I have simulated protons passing through the gas ionization chamber , according to the result, I found that there is a big difference between the number of electrons generated by the proton and finally collected by the copper strips . when the proton number is 1000, the generated electrons would be 3759 while the collected electrons would be 320, nearly one tenth .when the proton number is 10000, the generated electrons would be 37235 while the collected electrons would be 3417.
this is the endpoint coordinates of the electrons collected by the copper strips
That’s the code:
const int MAX = 100000;
unsigned int nesum = 0;
unsigned int necollect = 0;
for (int iEvent = 0; iEvent < MAX; ++iEvent) {
double xc = RndmGaussian(-10.7044,0.4) ;
double yc = RndmGaussian(18.99,0.4) ;
double zc = RndmGaussian(-7.53,0.4) ;
//double xc = -13.704 + RndmUniform() * 6;
//double yc = 18.49 + RndmUniform() * 1;
//double zc = -10.5 + RndmUniform() * 6;
//rnd.Gaussian2D(sigmaxc, sigmayc, 0.9, xc, yc);
//if (!((-13.704 < xc < -7.704) && (18.492 < yc < 19.482)&&(-10.53 < zc < -4.53)))continue;
hX->Fill(xc);
hY->Fill(yc);
hZ->Fill(zc);
hXY->Fill(xc, yc);
hXZ->Fill(xc, zc);
hYZ->Fill(yc, zc);
sensor.NewSignal();
std::cout<<"Launching new Track with coordinates (x,y,z) = ("<<xc<<","<<yc<<","<<zc<<")"<<std::endl;
tr.NewTrack(xc, yc, zc, 0., 0., -1., 0.);
// Loop over the clusters.
double 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";
necollect += 1;
hHitMap.Fill(xf, zf);
// drift.SetIonSignalScalingFactor(ne);
//drift.DriftIon(xc, yc, zc, tc);
}
}
std::cout << "nesum: (" << nesum << ")\n";
std::cout << "necollect: (" << necollect << ")\n";
Hi,
the reason for the discrepancy is that necollect counts the number of clusters on a track, while nesum counts the number of electron/ion pairs along the track. With TrackSrim and TrackTrim you typically have more than one electron/ion pair in each cluster.
Hi, Thanks a lot.
So how should I get the number of the endpoint while not the cluster number?
necollect += ne;
You don’t apply any filter on the start or end point of the electron drift line, so necollect == nesum.
Hi, thanks a lot for your help and happy new year.
What means filter? what’s that work for ? How should I use it?
Sorry for being unclear. What I meant is that you simulate an electron drift line for every “cluster” along the track. You don’t check where the end point of the drift line is, so the number of “end points” is the same as the number of clusters.
So the problem is I don’t know how to check the end point of the drift line. So initially I just count the number of (xf, yf,zf)
// multiple tracks and calculate the average collected charge
#include <iostream>
#include <TApplication.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TH2D.h>
#include <TCanvas.h>
#include "Garfield/MediumSilicon.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackTrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include <ctime>
#include <cassert>
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"
#include "gsl/gsl_vector.h"
#include "gsl/gsl_version.h"
#include "Math/GSLRndmEngines.h"
#include "TFile.h"
using namespace Garfield;
using namespace ROOT::Math;
int main(int argc, char *argv[]) {
// Application
TApplication app("app", &argc, argv);
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "m");
fm.PrintRange();
fm.PrintMaterials();
fm.SetWeightingField("weight.lis", "readout");
double x0, y0, z0, x1, y1, z1;
fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);
// Define the medium.
MediumMagboltz gas;
gas.LoadGasFile("720Torr-NEW.gas");
fm.SetMedium(0, &gas);
ViewMedium view;
view.SetMedium(&gas);
Sensor sensor;
sensor.AddComponent(&fm);
sensor.AddElectrode(&fm, "readout");
// Set the time bins for the induced current.
const unsigned int nTimeBins = 800;
const double tmin = 0.;
//onst double tmax = 2000.;
const double tmax = 800.;
const double tstep = (tmax - tmin) / nTimeBins;
sensor.SetTimeWindow(tmin, tstep, nTimeBins);
// Read the TRIM output file.
TrackTrim tr;
const std::string filename = "EXYZ-720He.txt";
// Import the first 100 ions.
if (!tr.ReadFile(filename)) {
std::cerr << "Reading TRIM EXYZ file failed.\n";
return 1;
}
tr.SetWorkFunction(42);
tr.SetFanoFactor(0.17);
// Connect the track to a sensor.
tr.SetSensor(&sensor);
AvalancheMicroscopic aval;
aval.SetSensor(&sensor);
AvalancheMC drift;
drift.SetSensor(&sensor);
drift.SetDistanceSteps(2.e-4);
// 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, -10.5, -4.5);
// Make a histogram of the electron energy distribution.
//TH1D hEn("hEn","energy distribution", 100, 0., 10.);
// aval.EnableElectronEnergyHistogramming(&hEn);
// Simulate the avalanche
ViewDrift driftView;
constexpr bool plotDrift = true;
if (plotDrift) {
aval.EnablePlotting(&driftView);
drift.EnablePlotting(&driftView);
}
//GSLRandomEngine rnd;
//rnd.Initialize();
const int MAX = 1000;
unsigned int nesum = 0;
unsigned int necollect = 0;
const double tEnd = 500.0;
const int nsBins = 500;
for (int iEvent = 0; iEvent < MAX; ++iEvent) {
double xc = RndmGaussian(-10.7044,0.4) ;
double yc = RndmGaussian(18.99,0.4) ;
double zc = RndmGaussian(-7.53,0.4) ;
//double xc = -13.704 + RndmUniform() * 6;
//double yc = 18.49 + RndmUniform() * 1;
//double zc = -10.5 + RndmUniform() * 6;
//rnd.Gaussian2D(sigmaxc, sigmayc, 0.9, xc, yc);
if (xc < -13.704 || xc > -7.704 || yc < 18.492 || yc > 19.482) continue;
// multiple tracks and calculate the average collected charge
// Calculate the avalanche.
aval.AvalancheElectron(xc, yc, zc, 0., 0.1, 0., -1., 0.);
std::cout << "... avalanche complete with "
<< aval.GetNumberOfElectronEndpoints() << " electron tracks.\n";
for (const auto& electron : aval.GetElectrons()) {
const auto& p0 = electron.path[0];
drift.DriftIon(p0.x, p0.y, p0.z, p0.t);
}
}
// Extract the calculated signal.
double bscale = tEnd / nsBins; // time per bin
double sum = 0.; // to keep a running sum of the integrated signal
// Create ROOT histograms of the signal and a file in which to store them.
TFile* f = new TFile("avalanche_signals.root", "RECREATE");
TH1F* hS = new TH1F("hh", "hh", nsBins, 0, tEnd); // total signal
TH1F* hInt = new TH1F("hInt", "hInt", nsBins, 0, tEnd); // integrated signal
// Fill the histograms with the signals.
// Note that the signals will be in C/(ns*binWidth), and we will divide by e
// to give a signal in e/(ns*binWidth).
// The total signal is then the integral over all bins multiplied by the bin
// width in ns.
for (int i = 0; i < nsBins; i++) {
double wt = sensor.GetSignal("wtlel", i) / ElementaryCharge;
sum += wt;
hS->Fill(i * bscale, wt);
hInt->Fill(i * bscale, sum);
}
// Write the histograms to the TFile.
hS->Write();
hInt->Write();
f->Close();
// Plot the signal.
const bool plotSignal = false;
if (plotSignal) {
TCanvas* cSignal = new TCanvas("signal", "Signal");
ViewSignal* vSignal = new ViewSignal();
vSignal->SetSensor(&sensor);
vSignal->SetCanvas(cSignal);
vSignal->PlotSignal("wtlel");
}
app.Run(true);
}
When I use these codes trying to figure out the number of endpoints through aval.GetNumberOfElectronEndpoints().
The results showed like this:
I’m not sure I understand the question (it would be helpful if you could explain what you are trying to do), but from the output you sent it looks like there is no multiplication happening in your chamber (at least not at the specific starting point you are using). Which would mean that the total number of electrons is the same as the number of “primary” electrons deposited along the track.
Hi, what I want is to get the number of “end points”, while I don’t know how to check the end point of the drift line. So I tired many ways, I used the AvalancheMC or the DriftLineRKF method, while I really don’t understand what’s the difference
With DriftLineRKF you can use the method GetAvalancheSize to retrieve the average number of electrons (and ions) at the end of the drift line. From the results you posted using AvalancheMC it looks like there is no multiplication happening (and little attachment), so to first order you can probably assume that the number of electrons after the drift in the chamber is the same as the number of primary electrons.
Hi, Thanks a lot for your reply. I understand now . If the multiplication happened, the number of electrons would be more than the primary electrons,however I don’t know what’s the key element not to happen the avalanche. Besides I do not want the endpoints of the drift line, but the electrons dropped on the copper strips (the collecting plate). That means the final numbers of electrons collected by the copper strips (the collecting plate).
What I mean is that when the proton inject into the gas ionization chamber, the numbers of electrons generated must be different from the one collected in the end point. because of the electron drifting way to somewhere or the electrons and ions recombination. So if the the number of electrons after the drift in the chamber is the same as the number of primary electrons, this situation probably not match the real condition.
Hi,
I’m afraid electron-ion recombination is not taken into account in any of the methods (DriftLineRKF, AvalancheMicroscopic, AvalancheMC) implemented in Garfield++ for simulation electron or ion drift lines.
You can use the coordinates of the drift line “end point” to check whether the drift line ends up on the collection electrode.
Hi,
I tried this method trying to figure out the number of the electrons collected by electrode. While it still can’t work.
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);
if (19.49 < yf < 19.5 ) {
++plane1;}
else if (18.98 < yf < 19) {
++cathode;}
else if (18.49< yf < 18.5) {
++plane2;}
else {
++escape;}
std::cout << "Endpoint: (" << xf << ", " << yf << ", " << zf << ")\n";
hHitMap.Fill(xf, zf);
// drift.SetIonSignalScalingFactor(ne);
//drift.DriftIon(xc, yc, zc, tc);
}
Yep ,I know that the Garfield++ can’t solve the recombination problem, is there any other software can do it?
I have read the example
and used TrackHeed trying to get the number of electrons at different places. While it is still failed.
// multiple tracks and calculate the average collected charge
#include <iostream>
#include <TApplication.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TH2D.h>
#include <TCanvas.h>
#include "Garfield/MediumSilicon.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackTrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include <ctime>
#include <cassert>
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"
#include "gsl/gsl_vector.h"
#include "gsl/gsl_version.h"
#include "Math/GSLRndmEngines.h"
using namespace Garfield;
using namespace ROOT::Math;
int main(int argc, char *argv[]) {
// Application
TApplication app("app", &argc, argv);
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "m");
fm.PrintRange();
fm.PrintMaterials();
fm.SetWeightingField("weight.lis", "readout");
double x0, y0, z0, x1, y1, z1;
fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);
// Define the medium.
MediumMagboltz gas;
gas.LoadGasFile("720Torr-NEW.gas");
fm.SetMedium(0, &gas);
ViewMedium view;
view.SetMedium(&gas);
Sensor sensor;
sensor.AddComponent(&fm);
sensor.AddElectrode(&fm, "readout");
// Set the time bins for the induced current.
const unsigned int nTimeBins = 800;
const double tmin = 0.;
//onst double tmax = 2000.;
const double tmax = 800.;
const double tstep = (tmax - tmin) / nTimeBins;
sensor.SetTimeWindow(tmin, tstep, nTimeBins);
// Read the TRIM output file.
TrackTrim tr;
const std::string filename = "EXYZ-720He.txt";
// Import the first 100 ions.
if (!tr.ReadFile(filename)) {
std::cerr << "Reading TRIM EXYZ file failed.\n";
return 1;
}
tr.SetWorkFunction(42);
tr.SetFanoFactor(0.17);
// Connect the track to a sensor.
tr.SetSensor(&sensor);
// Set up the charged particle track.
TrackHeed track;
track.SetSensor(&sensor);
track.SetParticle("proton");
// Set the momentum [eV / c].
track.SetMomentum(1.e9);
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, -10.5, -4.5);
// Make a histogram of the electron energy distribution.
//TH1D hEn("hEn","energy distribution", 100, 0., 10.);
// aval.EnableElectronEnergyHistogramming(&hEn);
// Simulate the avalanche
//GSLRandomEngine rnd;
//rnd.Initialize();
TH1F* hX = new TH1F("hX", "hX;x;Counts", 100, -15, -5);
TH1F* hY = new TH1F("hY", "hY;y;Counts", 100, 18, 20);
TH1F* hZ = new TH1F("hZ", "hZ;z;Counts", 100, -11, -4);
TH2F* hXY = new TH2F("hXY", "hXY;x;y;Counts", 100, -13.7, -7.7, 100, 18.49, 19.48);
TH2F* hXZ = new TH2F("hXZ", "hXZ;x;z;Counts", 100, -13.7, 5-7.7, 100, -10.5, -4.5);
TH2F* hYZ = new TH2F("hYZ", "hYZ;y;z;Counts", 100, 18.49, 19.48, 100, -10.5, -4.5);
const int MAX = 1000;
unsigned int nesum = 0;
unsigned int necollect = 0;
int plane1 = 0, plane2 = 0, cathode = 0, escape = 0;
double tc;
for (int iEvent = 0; iEvent < MAX; ++iEvent) {
double xc = RndmGaussian(-10.7044,0.4) ;
double yc = RndmGaussian(18.99,0.4) ;
double zc = RndmGaussian(-7.53,0.4) ;
//double xc = -13.704 + RndmUniform() * 6;
//double yc = 18.49 + RndmUniform() * 1;
//double zc = -10.5 + RndmUniform() * 6;
//rnd.Gaussian2D(sigmaxc, sigmayc, 0.9, xc, yc);
if (xc < -13.704 || xc > -7.704 || yc < 18.492 || yc > 19.482) continue;
hX->Fill(xc);
hY->Fill(yc);
hZ->Fill(zc);
hXY->Fill(xc, yc);
hXZ->Fill(xc, zc);
hYZ->Fill(yc, zc);
sensor.NewSignal();
std::cout<<"Launching new Track with coordinates (x,y,z) = ("<<xc<<","<<yc<<","<<zc<<")"<<std::endl;
tr.NewTrack(xc, yc, zc, 0., 0., -1., 0.);
for (const auto& cluster : track.GetClusters()) {
// Retrieve the electrons of the cluster.
for (const auto& electron : cluster.electrons) {
//sensor.EnableComponent(0, true);
//sensor.EnableComponent(1, false);
const double xc = electron.x;
const double yc = electron.y;
const double zc = electron.z;
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);
//sensor.EnableComponent(1, true);
//sensor.EnableComponent(0, false);
if (yf > 19.49 && yf< 19.5 ) {
++plane1;}
else if (yf > 18.98 && yf < 19) {
++cathode;}
else if (yf > 18.49 && yf < 18.5) {
++plane2;}
else {
++escape;}
std::cout << "Endpoint: (" << xf << ", " << yf << ", " << zf << ")\n";
hHitMap.Fill(xf, zf);
//std::cout << "nesum: (" << nesum << ")\n";
}
}
// drift.SetIonSignalScalingFactor(ne);
//drift.DriftIon(xc, yc, zc, tc);
}
std::cout << "plane1: (" << plane1 << ")\n";
std::cout << "cathode: (" << cathode << ")\n";
std::cout << "plane2: (" << plane2 << ")\n";
std::cout << "escape: (" << escape << ")\n";
//std::cout << "necollect: (" << necollect << ")\n";
TCanvas c1= new TCanvas("c1", "electron collected numbers");
hHitMap.Draw( );
//TCanvas c2;
//hEn.Draw("energy");
TCanvas* d = new TCanvas("d", "Multivariate gaussian random numbers");
d->Divide(3, 2);
d->cd(1);
hX->Draw();
d->cd(2);
hY->Draw();
d->cd(3);
hZ->Draw();
d->cd(4);
hXY->Draw("COL");
d->cd(5);
hXZ->Draw("COL");
d->cd(6);
hYZ->Draw("COL");
app.Run(true);
}
Why did you switch to TrackHeed (instead of TrackTrim)?
How many clusters do you get along the track?
Hi,
Thanks a lot for your reply
There is no clusters get along the track using the TrackHeed.
I am just confused for not getting the number of electrons at different places using the if-else sentences through the Tracktrim method. So switching to the TrackHeed is also not wise.