Hello! I want to retrieve the signal of low energy ions and fission fragments (in the script I am attaching below,I want to track α-particles with kinetic energy of 5 MeV) passing through a Micromegas geometry.The code I am using is a combination of asacusa,DriftTube,SRIM and other examples.Retrieving the outputs Ι am attaching after the code,I had the following questions:
-
I see that,in most examples,functions DriftLineRKF,AvalancheMicroscopic,AvalancheMC are being used for signal calculation.In my script,I am using only DriftLineRKF as you can see. I want to know if this is the most efficient (or even a correct!) way to retrieve the signal or should I use a combination of AvalancheMicroscopic,AvalancheMC and DriftLineRKF?
-
What exactly is the definition of dx,dz,dz in TrackSrim::NewTrack?.Firstly I thought it was just a definition of direction using a beginning (x0,y0,z0) and an ending point (dx0,dy0,dz0) of the track but I noticed that maybe I misunderstood.Are they related to the angles? Can you make it clear to me?
-
Using DriftLineRKF::GetElectronSignal,I am always getting 0.Could you tell me what am I missing here? Also,I am wondering why the total signal on mesh2 is positive,even though the electrons are being collected there.Is it maybe due to the large number of positive ions in amplification region and so the total signal is positive?
-
Is the usage of DriftLineRKF::SetElectronSignalScalingFactor(nc) correct in the code (Ln 149) or should this be in/out of the loop?
Any answer to these or any suggestion/correction on the code would be Really helpful.
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"
using namespace std;
using namespace Garfield;
bool readTransferFunction(Sensor& sensor) {
std::ifstream infile;
infile.open("mdt_elx_delta.txt", std::ios::in);
if (!infile) {
std::cerr << "Could not read delta response function.\n";
return false;
}
std::vector<double> times;
std::vector<double> values;
while (!infile.eof()) {
double t = 0., f = 0.;
infile >> t >> f;
if (infile.eof() || infile.fail()) break;
times.push_back(1.e3 * t);
values.push_back(f);
}
infile.close();
sensor.SetTransferFunction(times, values);
return true;
}
int main(int argc, char * argv[]) {
TApplication app("app", &argc, argv);
// Set the gas mixture.
MediumMagboltz gas;
gas.LoadGasFile("donegas.gas");
const std::string path = std::getenv("GARFIELD_HOME");
gas.LoadIonMobility(path + "/Data/IonMobility_Ar+_Ar.txt");
// Build the Micromegas geometry.
// We use separate components for drift and amplification regions.
ComponentAnalyticField cmpD;
ComponentAnalyticField cmpA;
cmpD.SetMedium(&gas);
cmpA.SetMedium(&gas);
// NO magetic field [Tesla].
// 500 um amplification gap.
constexpr double yMesh = 0.005;
constexpr double vMesh = -300.;
cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
cmpA.AddPlaneY(0., 0., "pad");
// 6 mm drift gap.
constexpr double yDrift = yMesh + 0.6;
constexpr double vDrift = -800.;
cmpD.AddPlaneY(yDrift, vDrift, "drift");
cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
cmpA.AddReadout("mesh1");
cmpD.AddReadout("mesh2");
// Assemble a sensor.
Sensor sensor;
sensor.AddComponent(&cmpD);
sensor.AddComponent(&cmpA);
sensor.SetArea(-4.75, 0., -4.75, 4.75, yDrift, 4.75);
sensor.AddElectrode(&cmpA, "mesh1");
sensor.AddElectrode(&cmpD, "mesh2");
// Set the time window [ns] for the signal calculation.
const double tMin = 0.;
const double tMax = 220.;
const double tStep = 0.05;
const int nTimeBins = int((tMax - tMin) / tStep);
sensor.SetTimeWindow(0., tStep, nTimeBins);
// Set the delta reponse function.
if (!readTransferFunction(sensor)) return 0;
sensor.ClearSignal();
TrackSrim track;
track.SetSensor(&sensor);
// Read SRIM output from file.
const std::string file = "alpha in 93ar-7co2 0.01-5mev.txt";
if (!track.ReadFile(file)) {
std::cerr << "Reading SRIM file failed.\n";
return 0;
}
// Set the initial kinetic energy of the particle (in eV).
track.SetKineticEnergy(5e6);
// Set the W value and Fano factor of the gas.
track.SetWorkFunction(30.0);
track.SetFanoFactor(0.3);
// Set A and Z of the gas (not sure what's the correct mixing law). how to compute a and z in mixed gases
const double za = 0.93 * (18. / 40.) + 0.07 * (22. / 44.);
track.SetAtomicMassNumbers(22. / za, 22);
track.SetTargetClusterSize(500);
// RKF integration.
DriftLineRKF drift;
drift.SetSensor(&sensor);
//drift.SetGainFluctuationsPolya(0., 20000.);
drift.EnableSignalCalculation();
/*
TCanvas* cS = nullptr;
ViewSignal signalView;
constexpr bool plotSignal = true;
if (plotSignal) {
cS = new TCanvas("cS", "", 600, 600);
signalView.SetCanvas(cS);
signalView.SetSensor(&sensor);
signalView.SetLabelY("signal [fC]");
}
*/
std::ofstream outfile1;
outfile1.open("signal.txt");
const unsigned int nTracks = 200;
for (unsigned int j = 0; j < nTracks; ++j) {
sensor.ClearSignal();
track.NewTrack(0, yDrift, 0, 0, 1.023, 0.0, 0);
double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
int nc = 0;
while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
drift.SetElectronSignalScalingFactor(nc);
for (int k = 0; k < nc; ++k) {
drift.DriftElectron(xc, yc, zc, tc);
drift.DriftIon(xc, yc, zc, tc);
}
}
sensor.ConvoluteSignals();
const double qmesh1 = sensor.GetSignal("mesh1", nTimeBins - 1);
const double qmesh2 = sensor.GetSignal("mesh2", nTimeBins - 1);
const double qemesh1 = sensor.GetElectronSignal("mesh1", nTimeBins - 1);
const double qemesh2 = sensor.GetElectronSignal("mesh2", nTimeBins - 1);
outfile1 << qmesh2 << " " << qmesh1 << " " << qemesh1 << " " << qemesh2 << "\n";
}
outfile1.close();
app.Run(kTRUE);
}
And the outputs for the signal:
signal.txt (4.1 KB)
