Ion signal with SRIM

[YiHui_Liu]

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.SetElectronSignalScalingFactor(nce);
        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:

Summary

#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;
    gas.LoadGasFile("helium3_625_ar_375_64atm_3000K.gas");
    const std::string path = std::getenv("GARFIELD_INSTALL");
    gas.LoadIonMobility(path + "/Data/IonMobility_He+_He.txt");

    ComponentAnalyticField cmp;
    cmp.SetMedium(&gas);

    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");
    cmp.AddReadout("s");

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

    TrackSrim track;
    track.SetSensor(&sensor);
    track.ReadFile("proton_in_he3-argon(gas).txt");
    track.SetKineticEnergy(0.65e6);

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

    TCanvas* cD = new TCanvas("cD", "", 600, 600);
    ViewCell cellView;
    ViewDrift driftView;
    cellView.SetCanvas(cD);
    cellView.SetComponent(&cmp);
    driftView.SetCanvas(cD);
    mc.EnablePlotting(&driftView);
    track.EnablePlotting(&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.SetElectronSignalScalingFactor(nce);
        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);

    cD->Clear();
    cellView.Plot2d();
    driftView.Plot(true, false);

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

    app.Run(kTRUE);

    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:

• Signal calculation of low energy particles and fission fragments in a Micromegas geometry

[couet]

Welcome tp theROOT forum

I think @hschindl can help you.

[hschindl]

Hi,
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).

[YiHui_Liu]

Thanks for your reply.

@hschindl

AvalancheMC::AvalancheElectron does not simulate the ion tail.

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?

@hschindl

You would need to loop over the “electron endpoints”

sorry, can you raise an example?

@hschindl

simulate an ion drift line from the starting point of each electron trajectory

I add the code to the loop:

mc.SetHoleSignalScalingFactor(nce);
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:

@hschindl

Alternatively, you could use DriftLineRKF which does simulate the ion tail (if you switch it on).

DriftLineRKF will give me a different result with AvalancheMC:

[YiHui_Liu]

Hello,
I use DriftLineRKF and change my code:

DriftLineRKF drift;
drift.SetSensor(&sensor);
drift.SetMaximumStepSize(5.e-4);
drift.EnableIonTail();
drift.EnableSignalCalculation();
...
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);
}

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

[hschindl]

@YiHui_Liu

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

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.

@YiHui_Liu

DriftLineRKF will give me a different result with AvalancheMC

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).

[hschindl]

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

@YiHui_Liu

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

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.

[YiHui_Liu]

@hschindl

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

Sure, here is:
example.zip (8.0 KB)

[YiHui_Liu]

@hschindl

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);
}

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'
what():std::bad_alloc

with:

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);
    }
}

@hschindl

The ion signal does look a bit strange.

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.SetIonSignalScalingFactor(nc);
    drift.DriftIon(xc, yc, zc, tc);
    drift.SetElectronSignalScalingFactor(nc);
    drift.DriftElectron(xc, yc, zc, tc);

    double gain = drift.GetGain() * drift.GetLoss();
    drift.SetIonSignalScalingFactor(gain);
    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

[hschindl]

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.

[YiHui_Liu]

Thanks!

@hschindl

Also the gas table you are using only extends up to 10 kV/cm.

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:

@hschindl

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.

According to Blanc’s Law, we know that:

\frac{1}{K_{mix}}=\sum\limits_i\frac{x_i}{K_i}

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?

[hschindl]

Hi,
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:
https://doi.org/10.1088/1748-0221/10/07/P07004