like what I typed in the topic, I want to get the totally number of the electrons generated after the proton passing through the gas ionization chamber ne, and in the meanwhile get the number that collected by the collecting plate ne’. Thus I can finally calculate the number of electrons which are recombined with ions or drifting somewhere else.
Now I can have the ne’ data, while to the number of ne I have no idea.
As discussed in the other thread, you can retrieve the coordinates of the end point of each electron trajectory using GetEndPoint and then count how many of the electrons ended on the collection electrode.
Hi,
Thanks a lot for your reply
I have calculate the number of electrons generated, the number is 4217,while the endpoint of the electrons is only more than 400.
Why is that?
Besides, how to simulate the incident proton beam that satisfies the Gaussian distribution.
heed.SetParticle("p");
heed.SetMomentum(120.e9);
heed.NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);
where x0, y0 are sampled from from a 2D Gaussian distribution,
see e.g. here: ROOT: tutorials/math/multivarGaus.C File Reference
Hi,Thanks a lot for your help.
I am trying to figure it out, while when I try to run the program, there are some warning to notice the .h file not existing or like this
Hi,
the compiler message is quite self-explanatory. Did you try “Gaussian2D” instead of “GaussianND”?
Maybe you need to add ROOT::MathMore to the list of libraries in your CMakeLists.txt?
add_executable(gem gem.C)
target_link_libraries(gem Garfield::Garfield ROOT::MathMore)
Yep, it works. Thanks a lot for your help.
But I don’t know how to set the Gaussian 2D parameter, what’s the sigmaX,sigmaY,rho,&x,&y mean? I tried to find some samples to read, but I can’t find more.
https://www.gnu.org/software/gsl/doc/html/randist.html#the-bivariate-gaussian-distribution
sigmaX, sigmaY seem to be the standard deviations in x and y, rho is the correlation coefficient; x and y are the random numbers you want to sample.
Hi,
Thanks a lot
But I don’t know how to get the random X and Y from the Gaussian 2D void.
const int MAX = 20;
unsigned int nesum = 0;
double genpars[2] = {0, 0};
for (int iEvent = 0; iEvent < MAX; ++iEvent) {
double xa = -13.7 + RndmUniform() * 6;
double ya = 18.4 + RndmUniform() * 1;
double za = -10.5 + RndmUniform() * 6;
genpars = rnd.Gaussian2D(2, 1, 0.9, xa, ya);
auto xc = genpars[0];
auto yc = genpars[1];
hX->Fill(xc);
hY->Fill(yc);
hZ->Fill(zc);
hXY->Fill(xc, yc);
hXZ->Fill(xc, zc);
hYZ->Fill(yc, zc);
sensor.NewSignal();
// Simulate an ion track.
tr.NewTrack(xc, yc, zc, 0., 0., -1., 0.);
// Loop over the clusters.
double tc, ec, ekin;
int ne = 0;
while (tr.GetCluster(xc, yc, zc, tc, ne, ec, ekin)) {
// Sum up the number of electron/ion pairs created.
nesum += ne;
// Simulate electron and ion drift lines starting
// from the cluster position.
// Scale the induced current by the number of electron/ion pairs
// in the cluster.
drift.SetElectronSignalScalingFactor(ne);
drift.DriftElectron(xc, yc, zc, tc);
double xf = 0., yf = 0., zf = 0., tf = 0.;
int status = 0;
drift.GetEndPoint(xf, yf, zf, tf, status);
std::cout << "Endpoint: (" << xf << ", " << yf << ", " << zf << ")\n";
hHitMap.Fill(xf, yf);
// drift.SetIonSignalScalingFactor(ne);
//drift.DriftIon(xc, yc, zc, tc);
}
double xa = 0., ya = 0.;
rnd.Gaussian2D(2, 1, 0.9, xa, ya);
should to the trick I think (but this is just my guess).
Yep,I have tried this way, while the result is like this:
I think maybe the range of the X and Y is :
Dimensions of the elementary block
-13.7044 < x < -7.70436 cm,
18.4927 < y < 19.4827 cm,
-10.5321 < z < -4.53211 cm,
-500 < V < 0.10642 V.
Periodicities
while the random xa and ya is out of this range, so the proton can’t inject into the effective area.
I have set the conditon
if (-13.704 < xc < -7.704) continue; if (18.492 < yc < 19.482) continue; while the result is the same.
GSLRandomEngine rnd;
rnd.Initialize();
TH1F* hX = new TH1F("hX", "hX;x;Counts", 100, -15, -5);
TH1F* hY = new TH1F("hY", "hY;y;Counts", 100, 18, 20);
TH1F* hZ = new TH1F("hZ", "hZ;z;Counts", 100, -11, -4);
TH2F* hXY = new TH2F("hXY", "hXY;x;y;Counts", 100, -13.7, -7.7, 100, 18.49, 19.48);
TH2F* hXZ = new TH2F("hXZ", "hXZ;x;z;Counts", 100, -13.7, 5-7.7, 100, -10.5, -4.5);
TH2F* hYZ = new TH2F("hYZ", "hYZ;y;z;Counts", 100, 18.49, 19.48, 100, -10.5, -4.5);
const int MAX = 13;
unsigned int nesum = 0;
double xc = -10.704;
double yc = 18.99;
for (int iEvent = 0; iEvent < MAX; ++iEvent) {
double zc = -10.5 + RndmUniform() * 6;
rnd.Gaussian2D(0.54, 0.0998, 0.9, xc, yc);
if ((-13.704 < xc < -7.704) && (18.492 < yc < 19.482))continue;
hX->Fill(xc);
hY->Fill(yc);
hZ->Fill(zc);
hXY->Fill(xc, yc);
hXZ->Fill(xc, zc);
hYZ->Fill(yc, zc);
sensor.NewSignal();
// Simulate an ion track.
tr.NewTrack(xc, yc, zc, 0., 0., -1., 0.);
// Loop over the clusters.
double tc, ec, ekin;
int ne = 0;
while (tr.GetCluster(xc, yc, zc, tc, ne, ec, ekin)) {
// Sum up the number of electron/ion pairs created.
nesum += ne;
// Simulate electron and ion drift lines starting
// from the cluster position.
// Scale the induced current by the number of electron/ion pairs
// in the cluster.
drift.SetElectronSignalScalingFactor(ne);
drift.DriftElectron(xc, yc, zc, tc);
double xf = 0., yf = 0., zf = 0., tf = 0.;
int status = 0;
drift.GetEndPoint(xf, yf, zf, tf, status);
std::cout << "Endpoint: (" << xf << ", " << yf << ", " << zf << ")\n";
hHitMap.Fill(xf, yf);
// drift.SetIonSignalScalingFactor(ne);
//drift.DriftIon(xc, yc, zc, tc);
}
}
std::cout << "nesum: (" << nesum << ")\n";
TCanvas c1;
hHitMap.Draw( );
//TCanvas c2;
//hEn.Draw("energy");
TCanvas* d = new TCanvas("d", "Multivariate gaussian random numbers");
d->Divide(3, 2);
d->cd(1);
hX->Draw();
d->cd(2);
hY->Draw();
d->cd(3);
hZ->Draw();
d->cd(4);
hXY->Draw("COL");
d->cd(5);
hXZ->Draw("COL");
d->cd(6);
hYZ->Draw("COL");
Sorry, what is your question?
my question is why my result always show the no medium at initial point? I want the result like the plot below, the proton injection is according to the Gaussian 2D distribution. the setting range of the X and Y obeys :
Dimensions of the elementary block
-13.7044 < x < -7.70436 cm,
18.4927 < y < 19.4827 cm,
-10.5321 < z < -4.53211 cm,
With this line:
if ((-13.704 < xc < -7.704) && (18.492 < yc < 19.482))continue;
you basically skip the xc and yc values you are interested in.
Please look up the use of continue statements in C++.
What you need is likely
if (!((-13.704 < xc < -7.704) && (18.492 < yc < 19.482)))continue;
you could have found yourself inserting more printouts in your code, which I strongly recommend.
For example:
std::cout<<"Launching new Track with coordinates (x,y,z) = ("<<xc<<","<<yc<<","<<zc<<")"<<std::endl;
tr.NewTrack(xc, yc, zc, 0., 0., -1., 0.);
These are not really garfield problems but more general C++ questions.
Yes, that’s really the problem is. when I use this code
if ((-13.704 < xc < -7.704) && (18.492 < yc < 19.482))continue;
it would jump out the loop cycle I just needed.
Thanks a lot for your kindly help, otherwise would you please tell me how to use the multithread computation