Ion signal with SRIM

I am trying to collect the signal produced by the proton in the proportional counter. The working gas is a mixture of He3 and Ar. The energy of the proton is between 0.6MeV and 0.8MeV. I used SRIM to generate the files under the corresponding conditions.
But there was a problem when I tried to collect the signal with AvalancheMC::AvalancheElectron() to simulate the electron drift:

    while (track.GetCluster(xc, yc, zc, tc, nce, eDepc, extra)) {
        mc.AvalancheElectron(xc, yc, zc, tc, true);

But I can not get the ion signal. I tried using either ‘IonMobility_He+_He.txt’ or ‘IonMobility_Ar+_Ar.txt’ separately, but I don’t know if that’s true, because protons might ionize both He3 and Ar.

My Code:

#include <TApplication.h>
#include <TCanvas.h>
#include <TROOT.h>

#include "Garfield/AvalancheMC.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackSrim.hh"

#include <cstdlib>
#include <fstream>
#include <iostream>

using namespace Garfield;

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

    MediumMagboltz gas;
    const std::string path = std::getenv("GARFIELD_INSTALL");
    gas.LoadIonMobility(path + "/Data/IonMobility_He+_He.txt");

    ComponentAnalyticField cmp;

    const double rWire = 1e-3;
    const double rTube = 2.54 / 2;
    const double vWire = 1000.;
    const double vTube = 0.;
    cmp.AddWire(0., 0., 2 * rWire, vWire, "s", 30.);
    cmp.AddTube(rTube, vTube, 0, "t");

    Sensor sensor;
    sensor.AddElectrode(&cmp, "s");
    sensor.SetTimeWindow(-0.25, 0.5, 4000);

    TrackSrim track;

    AvalancheMC mc;

    TCanvas* cD = new TCanvas("cD", "", 600, 600);
    ViewCell cellView;
    ViewDrift driftView;

    int nc = 0, ne = 0, nce;
    double eDepSum = 0;
    double xc, yc, zc, tc, eDepc, extra;

    const double rTrack = 0.3;
    track.NewTrack(rTrack, rTrack, 0, 0, -1, 1, 0);
    while (track.GetCluster(xc, yc, zc, tc, nce, eDepc, extra)) {
        nc += 1;
        ne += nce;
        eDepSum += eDepc;

        mc.AvalancheElectron(xc, yc, zc, tc, true);
        // mc.SetHoleSignalScalingFactor(nce);
        // mc.AvalancheHole(xc, yc, zc, tc, true);

    printf("Cluster = %d, Total Count = %d, Energy Deposit = %.3f, W = %.3f\n", nc, ne, eDepSum, eDepSum / ne);

    driftView.Plot(true, false);

    ViewSignal signalView;
    signalView.PlotSignal("s", "tei");


    return 0;

helium3_625_ar_375_64atm_3000K.txt (103.4 KB)

proton_in_he3-argon(gas).txt (6.1 KB)

My Result:

Related Topic:

Welcome tp theROOT forum

I think @hschindl can help you.

AvalancheMC::AvalancheElectron does not simulate the ion tail. You would need to loop over the “electron endpoints” and simulate an ion drift line from the starting point of each electron trajectory. Alternatively, you could use DriftLineRKF which does simulate the ion tail (if you switch it on).

Thanks for your reply.

bool AvalancheElectron(const double x, const double y, const double z,
                         const double t, const bool hole = false);

I set hole to true, it should simulate the ions in the electron avalanche?

sorry, can you raise an example?

I add the code to the loop:

mc.AvalancheHole(xc, yc, zc, tc, true);

but it doesn’t work, and will raise an error:

AvalancheMC::GetVelocity: Error calculating hole velocity at (0.299768, 0.300192, -0.000018).
AvalancheMC::DriftLine: Abandoned the calculation.

and the result is similar to the previous one:

DriftLineRKF will give me a different result with AvalancheMC:

I use DriftLineRKF and change my code:

DriftLineRKF drift;
while (track.GetCluster(xc, yc, zc, tc, nce, eDepc, extra)) {

    drift.DriftIon(xc, yc, zc, tc);
    drift.DriftElectron(xc, yc, zc, tc);

I can see the tracks and signal of ions now, but it obviously does not take into account the contribution of ions in the electron-ion pair produced by the electron avalanche.

(because the electric field is weak, it takes a long time for heavier ions to be collected

another question is that the program only calculate Ar ions without considering possible He3 ions

No, this flag applies only to holes in semiconductors at the moment. But now that you mention it I could modify the code such that when this flag is set to true, the ion drift lines in the avalanche will also be simulated. I’ll put this on my to-do list…

To loop over the electron end points you can do something like this:

const unsigned int np = mc.GetNumberOfElectronEndpoints();
double xe1, ye1, ze1, te1, e1;
double xe2, ye2, ze2, te2, e2;
int status;
for (unsigned int j = 0; j < np; ++j) {
  mc.GetElectronEndpoint(j, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);
  mc.DriftIon(xe1, ye1, ze1, te1);

AvalancheHole only make sense for semiconductors; that’s why you get the error message.

DriftLineRKF does not take diffusion into account, that’s why the electron signal is “sharper” than with AvalancheMC. Also with AvalancheMC the avalanche size will fluctuate from avalanche to avalanche, while DriftLineRKF takes by default the average gain (computed by integrating the Townsend coefficient along the drift line).

The ion signal does look a bit strange. Can you provide a minimal working example that reproduces this plot?

Yes, that’s a good point; the treatment of ion transport is indeed a bit of an issue. The table of ion mobilities that you load at the beginning of the program should be understood as “effective” mobility of ions in the gas.

Sure, here is: (8.0 KB)

xe1, ye1, ze1 is the start point of the electron drift line, does it is the same point to the xc, yc or cluster.x?

and I got:

terminate called after throwing an instance of 'std::cad_alloc'


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) {
        mc.GetElectronEndpoint(i, xe1, ye1, ze1, te1, xe2, ye2, ze2, te2, status);
        mc.DriftIon(xe1, ye1, ze1, te1);

I think that the reason is the program does not take into account the contribution of ions in the electron-ion pair produced by the electron avalanche.
So I want to get the Gain of electron drift, and then driftIon in the endpoint of the electron drift:

while (track.GetCluster(xc, yc, zc, tc, nc, eDepc, extra)) {
    drift.DriftIon(xc, yc, zc, tc);
    drift.DriftElectron(xc, yc, zc, tc);

    double gain = drift.GetGain() * drift.GetLoss();
    drift.DriftIon(xc, yc, zc, tc);

the document shows that:

n_e= \exp(\int(\alpha-n)\mathrm{d}s)= \exp(\int\alpha\mathrm{d}s)\times\exp(-\int n\mathrm{d}s)=Gain\times Loss

but it always return 1 as the result of GetGain, the source sode:

double DriftLineRKF::GetGain(const double eps) const {
  if (m_status == StatusCalculationAbandoned) return 1.;
  return ComputeGain(m_x, m_particle, eps);

why m_status will equal to StatusCalculationAbandoned?

  • In DriftLineRKF::DriftLine, it shows that m_sensor may be not set.
  • In DriftLineRKF::Terminate, I don’t know exactly what it does

At the gas pressure and voltage settings you are using, the electric field is not high enough to start an avalanche. Also the gas table you are using only extends up to 10 kV/cm. At a pressure of 6.4 atm there is no multiplication yet at this field.
Can you check what are the parameters of the gas and the voltages that are used in the detector which you are simulating?
In any case, I think you’ll need to increase the range of electric fields in your gas table significantly. Something like 100 kV/cm is a typical value for atmospheric pressure; you will probably need to go a bit higher.


I got the reason, the voltage is high enough to start an avalanche. In the near-wire area, the electric filed will be up to 300kV/cm, so I regenerate the gas table extending up to 320kV/cm, and now I get the result:

According to Blanc’s Law, we know that:


where Kmix is the standard mobility of ion in gases mixture, xi is the mole fraction of the i-th component in the mixture and Ki is the standard mobility of ions in the i-th component.

Now I can get the mobility of He+ in He, Ar+ in He and Ar+ in Ar (how to get He+ in Ar is another question). So I need to know how many He+ and Ar+ are produced by the avalanche to determine the final ion mobility data, is there any method to get this?

I don’t think there is a way to get information from SRIM about the ion species that are produced by the projectile.
Also you would need to consider charge transfer reactions between the ions to determine which are the relevant types of ions in your case. For more details on this, see e. g. this paper:

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.