Dear experts,
I’m currently working on protons crossing a plane parallel ionisation chamber geometry and studying the attachment. I have some problems when i try to understand the visualization of the drift lines.
- Why my proton (green line) stop in the chamber while the energy it loses is negligible ?
- Why does it look like there are 2 electrons, one collected by the electrode and the other attached, when the code give a total number of electrons of 1?
- In the second photo, i have a drift line of electron who stop but all the status are -7 (attachment) or -5 (lft the drift medium =? collected by the electrode) you can see the status printed in the terminal to the right.
One last question: is it possible to change the color of the positive or negative ion to distinguish them ?
Thank you for your help.
Pierre
// Make a medium
MediumMagboltz gas("o2", 21.1, "Ar", .9, "H2O", 2., "N2", 76.);
gas.SetTemperature(273.15 + 20);
gas.SetPressure(760.);
const std::string path = getenv("GARFIELD_INSTALL");
gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_N2+_N2.txt");
gas.LoadNegativeIonMobility(path + "/share/Garfield/Data/IonMobility_O2-_air.txt");
// Thickness of the gas gap [cm] and voltage [V]
constexpr double gap = .1;
constexpr double voltage = 500.;
// Drift canvas
const bool plotDrift = true;
TCanvas* cD = nullptr;
if (plotDrift) cD = new TCanvas("cD", "", 600, 600);
// Make a component
ComponentConstant cmp;
cmp.SetArea(0., -1., -1., gap, 1., 1.);
cmp.SetMedium(&gas);
double field = voltage / double(gap);
cmp.SetElectricField(field, 0., 0.);
// Make a sensor
Sensor sensor;
sensor.AddComponent(&cmp);
// DriftView class
ViewDrift driftView;
float ViewWidth = .02;
driftView.SetArea(0., -ViewWidth, -ViewWidth, gap, ViewWidth, ViewWidth);
if (plotDrift) driftView.SetCanvas(cD);
// AvalancheMicroscopic class (for electrons)
AvalancheMicroscopic aval;
aval.SetSensor(&sensor);
aval.EnablePlotting(&driftView);
aval.EnableExcitationMarkers(true);
// AvalancheMC class (for ions)
AvalancheMC drift;
drift.SetSensor(&sensor);
drift.EnablePlotting(&driftView);
// Track class
TrackHeed track;
track.SetSensor(&sensor);
track.SetParticle("proton");
//track.SetKineticEnergy(230.e6); // 230 MeV
track.SetGamma(3); // using MIP for more visibility in the drift view
constexpr bool verbose = true;
track.Initialise(&gas, verbose);
track.EnablePlotting(&driftView);
// total number of electrons
int nsum = 0;
// number of attached electrons
int asum = 0;
// energy loss
double esum = 0;
// Initial position and direction
double x0 = 0., y0 = 0., z0 = 0., t0 = 0.;
double dx0 = 1., dy0 = 0., dz0 = 0.;
track.NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);
for (const auto& cluster : track.GetClusters()) {
esum += cluster.energy;
for (const auto& electron : cluster.electrons) {
aval.AvalancheElectron(electron.x, electron.y, electron.z, electron.t, 0.1, 0, 0, 0);
// Loop over the electron avalanches.
for (int j = 0; j < aval.GetNumberOfElectronEndpoints(); ++j){
nsum ++;
double xe1, ye1, ze1, te1, e1;
double xe2, ye2, ze2, te2, e2;
int status;
aval.GetElectronEndpoint(j, xe1, ye1, ze1, te1, e1, xe2, ye2, ze2, te2, e2, status);
aval.DriftElectron(xe1, ye1, ze1, te1, .1, 0, 0, 0);
std::cout << "status: " << status << "\n";
drift.DriftIon(xe1, ye1, ye1, te1);
if (status == -7) {
asum ++;
drift.DriftNegativeIon(xe2, ye2, ye2, te2);
}
}
}
}
double attachRatio = (asum / (double)nsum);
// Print the results
std::cout << "\n" << "Gas: "<< "\n";
for (int i = 0; i < gas.GetNumberOfComponents(); ++i) {
std::string name;
double f;
gas.GetComponent(i, name, f);
std::cout << " " << name << " : " << f*100 << "\n";
}
std::cout << " " << gas.GetTemperature()-273.15 << " [°C]" << "\n";
std::cout << " " << gas.GetPressure() << " [Torr]" << "\n\n";
std::cout << "Gap: " << gap << " [cm]" << "\n";
std::cout << "Voltage: " << voltage << " [V]" << "\n";
std::cout << "Field: " << field << " [V/cm]" << "\n\n";
std::cout << "Energy of particle: " << track.GetEnergy() << " eV" << "\n";
std::cout << "Total energy loss: " << esum << " eV" << "\n\n";
std::cout << "Total electrons: " << nsum << "\n";
std::cout << "Attached electrons: " << asum << "\n";
std::cout << "Ratio of attached electrons: " << attachRatio << "\n";
std::cout << "Free electron fraction: " << 1 - attachRatio << "\n";
if (plotDrift) {
constexpr bool twod = true;
constexpr bool drawaxis = true;
driftView.Plot(twod, drawaxis);
}