Signal with Srim(using AvalancheMC)

My code is as follows. I hope to collect electron and ion signals from sensors.

int main(int argc, char * argv[]) {
TApplication app("app", &argc, argv);
MediumMagboltz gas;

gas.LoadGasFile("ar_85_co2_15.gas");
const std::string path = std::getenv("GARFIELD_HOME");
gas.LoadIonMobility(path + "/Data/IonMobility_Ar+_Ar.txt");


ComponentAnalyticField cmp;
cmp.SetMedium(&gas);

//
const double rWire = 30.e-4;//cm
const double rTube = 3.5;//cm
//
const double vWire = 1100.;
const double vTube = 0.;
//
cmp.AddWire(0 , 0, 2* rWire, vWire, "s");
cmp.AddTube(rTube, vTube, 0 , "t");
cmp.AddReadout("s");

//
Sensor sensor;
sensor.AddComponent(&cmp);
sensor.AddElectrode(&cmp, "s");
const double tstep = 0.1;     
const double tmin = 0;     
const double tmax=20000;             
const unsigned int nbins = (tmax-tmin)/tstep;   
sensor.SetTimeWindow(tmin, tstep, nbins);

sensor.ClearSignal();
 //----alpha in -----
// Read the SRIM output file. 
TrackSrim track(&sensor);
const std::string file = "Li_in_Ar_CO2_gas.txt";
if (!track.ReadFile(file)) {
std::cerr << "Reading Li_in_Ar_CO2_gas.txt file failed.\n";
return 1;
}

// Set the initial kinetic energy of the particle (in eV).
 track.SetKineticEnergy(0.84e6);
track.SetWorkFunction(27.7);      
track.SetFanoFactor(0.172);      
const double za = 0.85 * (18. / 39.948) +0.15 * (22. / 44.0095);  // Set A and Z of the gas     (not sure what's the correct mixing law).
track.SetAtomicMassNumbers(22. / za, 22);
 // Specify how many electrons we want to be grouped to a cluster.

track.EnableTransverseStraggling(false);
track.EnableLongitudinalStraggling(false);


AvalancheMC mc;
mc.SetSensor(&sensor);
mc.EnableSignalCalculation();


const double rTrack =3.0;     
const double projectilex= rTrack;
const double projectiley= sqrt(rTube * rTube - rTrack * rTrack);
const double projectilez= 0.;  

const double projectiledx=0.;    
const double projectiledy=-1.; 
const double projectiledz=0.; 

const double projectileT=0.;  
//const unsigned int nTracks = 1;  

track.NewTrack(projectilex, projectiley, projectilez, projectileT, projectiledx, projectiledy,       projectiledz);  

for (const auto& cluster : track.GetClusters()) {
//mc.SetElectronSignalScalingFactor(cluster.n);
mc.AvalancheElectron(cluster.x, cluster.y, cluster.z, cluster.t);
int status;
const unsigned int np = mc.GetNumberOfElectronEndpoints();
double xe1, ye1, ze1, te1, xe2, ye2, ze2, te2;
for (unsigned int i = 0; i < np; ++i) {
  mc.GetElectronEndpoint(i, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);
  mc.DriftIon(xe1, ye1, ze1, te1);
}
}

ViewSignal signalView;
signalView.SetSensor(&sensor);
signalView.PlotSignal("s");

app.Run();
return 0;
}

But I encountered the following problem:

Hi,
calling DriftIon resets the containers of electron and ion “endpoints” stored in the AvalancheMC object. The easiest solution is probably to create two AvalancheMC objects, one for simulating the electron avalanche and one for transporting the ions.

Hi!
Thank you for your reply.I change the code to this way, and now the program can run:

AvalancheMC mc;
mc.SetSensor(&sensor);
mc.EnableSignalCalculation();
track.NewTrack(projectilex, projectiley, projectilez, projectileT, projectiledx, projectiledy,       
projectiledz);  
for (const auto& cluster : track.GetClusters()) {
mc.AvalancheElectron(cluster.x, cluster.y, cluster.z, cluster.t);
int status;
const unsigned int np = mc.GetNumberOfElectronEndpoints();
double xe1, ye1, ze1, te1, xe2, ye2, ze2, te2;
for (unsigned int i = 0; i < np; ++i) {
  AvalancheMC mcion;
  mcion.SetSensor(&sensor);
  mcion.EnableSignalCalculation();
  mc.GetElectronEndpoint(i, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);
  mcion.DriftIon(xe1, ye1, ze1, te1);
}

But I still have a few questions:
1.The calculation of program results is very slow and may take 3-4 hours.So I want to add parallel computing code as follows:

  omp_set_num_threads(62);
  #pragma omp parallel for
  for (const auto& cluster : track.GetClusters()) {
  mc.AvalancheElectron(cluster.x, cluster.y, cluster.z, cluster.t);
  int status;
  const unsigned int np = mc.GetNumberOfElectronEndpoints();
  double xe1, ye1, ze1, te1, xe2, ye2, ze2, te2;
  for (unsigned int i = 0; i < np; ++i) {
    AvalancheMC mcion;
    mcion.SetSensor(&sensor);
    mcion.EnableSignalCalculation();
    mc.GetElectronEndpoint(i, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);
    mcion.DriftIon(xe1, ye1, ze1, te1);
  }

But I’m not sure if this can cause electronic drift in parallel to the clusters in the track.

2.Do I need to configure SetTargetClusterSize(),or let it remain default?What impact will this cause on the results?

Hi,

  1. You will have to create a separate AvalancheMC object for each cluster (i. e. inside the for loop). Note that in TrackSrim, clusters typically have more than one electron, so you’ll probably want to scale the induced signal (both for the electrons and the ions) by the number of electrons in the cluster (SetElectronSignalScalingFactor). Also you don’t need to call EnableSignalCalculation explicitly; signals are calculated by default now if you have electrodes defined.
  2. In TrackSrim the grouping of electrons to a cluster is a bit arbitrary. Please consult the documentation for more details.

Thanks a lot! I still have some questions to ask.I want to simulate the current output signal under the same parameters multiple times and observe the waveform characteristics of these output signals. At a working voltage of 1300V, due to the significant electron multiplication, both AvalancheMC and AvalancheMicroscopic methods require 5 hours of simulation time for only one iteration. I have tried using the DriftLineRKF method, but the output signal is very sharp, losing a lot of information about the leading and trailing edges of the current signal.
I want to know what I should do if I want to retain more current waveform information but simulate the results as quickly as possible.

Dear @Migu
Since you are simulating the response to alpha particles, I would believe you have a lot of primary electrons that are amplified in your detector. In what kind of detector are you simulating your alpha particles? A MPGD, a wire detector, an RPC, a TPC? Depending on the geometry (especially on the possibility of electrons being lost on insulating materials inside your detector) you could think of simulating the initial avalanche with Microscopic Tracking, and once your primary electron was multiplied with a factor 50-100 you pass your electrons to DriftLineRKF. Alternatively - depending on your detector geometry - you could think of simulating only the first half of your avalanche, and extrapolate to the full avalanche. This was done in the past for micromegas [1] and worked out quite well.
greets
Piet

[1] 10.1016/j.nima.2014.11.014

Dear @Piet
I am simulating the current signal from DriftTube.

Depending on the geometry (especially on the possibility of electrons being lost on insulating materials inside your detector) you could think of simulating the initial avalanche with Microscopic Tracking, and once your primary electron was multiplied with a factor 50-100 you pass your electrons to DriftLineRKF .

I will try the simulation in the way you mentioned above!
But when I only used the DriftLineRKF method before(1300V), I encountered the following problem:
Integrating the Townsend coefficients would lead to exponential overflow.
The code is as follows:

DriftLineRKF drift;
drift.SetSensor(&sensor);
drift.SetMaximumStepSize(5.e-4);
drift.EnableIonTail();
drift.EnableSignalCalculation();
track.NewTrack(projectilex, projectiley, projectilez, projectileT, projectiledx, projectiledy,       
projectiledz);  
while (track.GetCluster(xc, yc, zc, tc, nce, eDepc, extra)) {
    drift.SetIonSignalScalingFactor(nce);
    drift.DriftIon(xc, yc, zc, tc);
    drift.SetElectronSignalScalingFactor(nce);
    drift.DriftElectron(xc, yc, zc, tc);
}

Best regards
Migu