How to find the electrons' deposition position coordinates on the collected copper strip

[XiCheng]

After the further coding in the below

DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.SetMaximumStepSize(0.01);
  drift.DriftElectron(-11.7, 19.4, -7.53, 0.);
  drift.GetDriftTime();
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
 // Create a 2D histogram with 100 bins between -1. and 1. (same range and binning for x and y).
  TH2F hHitMap("hHitMap", "Hit map", 1000, -14., -6., 1000, -11., -3.);
   
  constexpr unsigned int nTracks = 50;
  unsigned int nesum = 0;
  for (unsigned int i = 0; i < nTracks; ++i) {
    sensor.NewSignal();
    // Simulate an ion track.
    tr.NewTrack(-11.7, 19.2, -7.53, 0., 0., -1., 0.);
    // Loop over the clusters.
    double xc, yc, zc, 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);
      //drift.SetIonSignalScalingFactor(ne);
      //drift.DriftIon(xc, yc, zc, tc);
      aval.AvalancheElectron(xc, yc, zc, tc, ec, 0, 0, 0);    
      for (const auto& electron : aval.GetElectrons()) {
        const auto& endPoint = electron.path.back();
        hHitMap.Fill(endPoint.x, endPoint.y);
      }
    }
  }
  
  TCanvas c1;
  hHitMap.Draw("coordinate");

the result is like this

what’s the meaning of the sentence “rate at …is not included in the current table”
and how to set the min-max scope to the right range

[hschindl]

• During initialisation, MediumMagboltz loads a table of the electron-molecule scattering rates as function of electron energy. By default, the upper limit of this table is 40 eV. If (during the simulation of an avalanche) an electron exceeds this energy, then the upper limit is increased and the table is recomputed, and a message “Increasing energy range to …” is printed. The message therefore doesn’t necessarily mean that there is a problem, but it does indicate that the electric field in your detector is quite high (not sure if this is what you expect).
• I’m not familiar with the “Error in min-max scope” warning but I guess it’s related to the line hHitMap.Draw("coordinate"). I don’t think TH2F has a draw option “coordinate” (see ROOT: THistPainter Class Reference).
• In the loop you simulate electron drift lines both using DriftLineRKF (line drift.DriftElectron(xc, yc, zc, tc);) and using AvalancheMicroscopic (line avalanche.AvalancheElectron(xc, yc, zc, tc, ec, 0, 0, 0);). I don’t think this makes sense. Also you should not use ec (which is the total energy deposited in the cluster) as the initial energy of the electron. If you want to simulate the electron drift lines using DriftLineRKF (which is what I would recommend in this case), you can do something like this:

for (unsigned int i = 0; i < nTracks; ++i) {
  sensor.NewSignal();
  // Simulate an ion track.
  tr.NewTrack(-11.7, 19.2, -7.53, 0., 0., -1., 0.);
  // Loop over the clusters.
  double xc, yc, zc, 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);
    hHitMap.Fill(xf, yf, ne);
    // drift.SetIonSignalScalingFactor(ne);
    // drift.DriftIon(xc, yc, zc, tc);
  }
}

If you really want to simulate the electron avalanche “microscopically” (i. e. using AvalancheMicroscopic you would need to do something like this:

for (unsigned int i = 0; i < nTracks; ++i) {
    sensor.NewSignal();
    // Simulate an ion track.
    tr.NewTrack(-11.7, 19.2, -7.53, 0., 0., -1., 0.);
    // Loop over the clusters.
    double xc, yc, zc, 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;
      for (int j = 0; j < ne; ++j) {
        aval.AvalancheElectron(xc, yc, zc, tc, 0.1, 0, 0, 0);
        for (const auto& electron : aval.GetElectrons()) {
          const auto& endPoint = electron.path.back();
          hHitMap.Fill(endPoint.x, endPoint.y);
        }
      }
    }
  }

[XiCheng]

Hi,That’s really kind of you to explain extraordinary detail to me.
Actually I don’t understand what’s the difference between the drift lines and the avalanche “microscopically”. I just want to get these results:

1. there are many copper strips to collect the electrons, I need to know the coordinate information, the electron finally dropped to which strip.
   2.I want to know the distribution of the electrons collected on every copper strips, some might have many ,others may less. like the plot below: the X coordinate is about the position, Y coordinate is about the number of electron collected. the middle position of the strips obviously collect the most electrons.
   
   

[XiCheng]

When I use the drift line to calculate. the result is always blank, that seems no electrons collected in the copper strip range. I have changed the coordinate direction from X and Y to X and Z or even X,Y,Z three dimension view. But it’s all like this:

now my code is like this:

  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
 // Create a 2D histogram with 100 bins between -1. and 1. (same range and binning for x and y).
  TH2F hHitMap("hHitMap", "Hit map", 1000, -14., -6., 1000, -11., -3.);
   // Make a histogram of the electron energy distribution.
  //TH1D hEn("hEn","energy distribution", 100, 0., 10.);
 // aval.EnableElectronEnergyHistogramming(&hEn);
    // Simulate the avalanche

  constexpr unsigned int nTracks =200;
  unsigned int nesum = 0;
  for (unsigned int i = 0; i < nTracks; ++i) {
    sensor.NewSignal();
  // Simulate an ion track.
    tr.NewTrack(-11.7, 19.4, -7.53, 0., 0., -1., 0.);
  // Loop over the clusters.
    double xc, yc, zc, 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);
      hHitMap.Fill(xf, zf);
    // drift.SetIonSignalScalingFactor(ne);
    // drift.DriftIon(xc, yc, zc, tc);
    }
  }
  TCanvas c1;
  hHitMap.Draw( );

[XiCheng]

When I simulate the electron avalanche “microscopically” , then the result would be stopped and showed the killed word.

[hschindl]

Hi,
so based what you describe in your earlier post, I propose you use DriftLineRKF to simulate the electron drift lines (not AvalancheMicroscopic).

Now to find out where the electrons go, I suggest you print out a few of the end points. Just add a line like

std::cout << "Endpoint: (" << xf << ", " << yf << ", " << zf << ")\n";

after

drift.GetEndPoint(xf, yf, zf, tf, status);

This should give you an idea where the electrons end up. You could also use ViewDrift to plot the drift lines by the way.

[XiCheng]

Hi,Thanks again for your reply.
I found out that the program cannot retrieve drift velocity at initial position, however the initial position seems noting wrong. Is that the reason for the EXYZ txt?

Sensor::SetTimeWindow: Resetting all signals.
TrackTrim::ReadFile: Read energy vs position for 12408 ions.
    Initial kinetic energy set to 100000 keV. Mass number: 1.01.
Endpoint: (-10.9999, 18.9827, -8.00042)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-11.000000, 18.840090, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-11.000000, 18.840090, -8.000000).
Endpoint: (-10.9999, 18.9827, -8.00045)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-11.000000, 18.722540, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-11.000000, 18.722540, -8.000000).
Endpoint: (-10.9999, 18.9827, -8.00042)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-11.000000, 18.678880, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-11.000000, 18.678880, -8.000000).
Endpoint: (-10.9999, 18.9827, -8.00042)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-10.999998, 18.579940, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-10.999998, 18.579940, -8.000000).
Endpoint: (-10.9999, 18.9827, -8.00046)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-11.000000, 18.626680, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-11.000000, 18.626680, -8.000000).
Endpoint: (-10.9999, 18.9827, -8.00043)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-10.999992, 18.522990, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-10.999992, 18.522990, -8.000000).
Endpoint: (-11, 18.9827, -8.00032)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-10.999995, 18.610670, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-10.999995, 18.610670, -8.000000).
Endpoint: (-10.9999, 18.9827, -8.00045)
DriftLineRKF::GetVelocity:
    Cannot retrieve drift velocity at (-10.999995, 18.519970, -8.000000).
DriftLineRKF::DriftLine:
    Cannot retrieve drift velocity at initial position (-10.999995, 18.519970, -8.000000).

[hschindl]

Did you load a gas file?

[XiCheng]

Yep,I have loaded the gas file

 // Define the medium.
  MediumMagboltz gas;
  gas.LoadGasFile("740Torr-NEW.gas");
  fm.SetMedium(0, &gas);
  ViewMedium view;
  view.SetMedium(&gas);

[XiCheng]

Now is two results, one is killed ,another is “cannot retrieve drift velocity at initial position”.
I have tried other gas file at different pressure, the same result “cannot retrieve drift velocity at initial position” comes. when I changed the TH2F to TH3F to see the three dimension plot, however it shows killed.

[XiCheng]

I have changed the coordinate information in the TH2F hHitMap(“hHitMap”, “Hit map”, 1000, -13.7, -7.7, 1000, 18.4, 19.4); and tr.NewTrack(-11.7, 19.4, -8, 0., 0., -1., 0.); there are some endpoint, and only one point in the plot.

 DriftLineRKF drift;
  drift.SetSensor(&sensor);
  drift.SetMaximumStepSize(0.01);
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);
 // Create a 2D histogram with 100 bins between -1. and 1. (same range and binning for x and y).
  TH2F hHitMap("hHitMap", "Hit map", 1000, -13.7, -7.7, 1000, 18.4, 19.4);
   // Make a histogram of the electron energy distribution.
  //TH1D hEn("hEn","energy distribution", 100, 0., 10.);
 // aval.EnableElectronEnergyHistogramming(&hEn);
    // Simulate the avalanche

  constexpr unsigned int nTracks =20;
  unsigned int nesum = 0;
  for (unsigned int i = 0; i < nTracks; ++i) {
    sensor.NewSignal();
  // Simulate an ion track.
    tr.NewTrack(-11.7, 19.4, -8, 0., 0., -1., 0.);
  // Loop over the clusters.
    double xc, yc, zc, 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);
    }
  }

[XiCheng]

@XiCheng

// Simulate an ion track.
    tr.NewTrack(-11.7, 19.4, -8, 0., 0., -1., 0.);

I think the basic reason is from this, the proton injected coordinate is not a position fixed but a better be a random coordinate in the range. So that would be a distribution of the endpoint but not only one just loop over and repeated in the( -11.7, 19.4, -8,) rejection point. Is that right? How can I realize it?

[hschindl]

To sample a random number between xmax and xmin you can do something like this:

double xt = xmin + RndmUniform() * (xmax - xmin);

But I think it would still be good to understand why you get the message Cannot get drift velocity at initial position. What does your geometry look like? Is there any reason why this point would not be inside the mesh of your field map?

[XiCheng]

@hschindl

What does your geometry look like?

the geometry is like this:

the copper strips has two directions, what I simulate in the program is only one direction like I plot in the below.

[XiCheng]

@hschindl

But I think it would still be good to understand why you get the message Cannot get drift velocity at initial position.

I use the random incident coordinates, when the track is 200, the distribution on the copper strips is like this:

when the track is 2000, the electrons collected in the strips are like this:

It is a chaotic distribution not a gaussian distribution I expected. So maybe it is wrong result

[XiCheng]

@hschindl

Is there any reason why this point would not be inside the mesh of your field map?

yep, that’s really weird for me, all the coordinate is set in the range of the dimensions. why the result shows no medium in some point or cannot get drift velocity at initial position. I don’t figure it out yet.