Hello,
So this program uses Garfield++ to simulate a ionization chamber filled with Isobutane gas at 58 torr and releases Al-26 Ions into it from the bottom. I’ve been trying to improve the simulations by incorporating SRIM and TRIM into it as well.
It’s not working well. I do have a TRIM relevent gas file called COLLISION.txt that I’m trying to incorporate, just trying to figure out how to do it right. This is the last version of the program that at least somewhat works.
Any thoughts you all might have would be appreciated.
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <math.h>
#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/TrackSrim.hh"
using namespace Garfield;
int main(int argc, char* argv[]) {
TApplication app("app", &argc, argv);
srand(time(NULL));
///////////////////////////////////////////////////////////////////////////////////////////////
// Make a gas medium.
MediumMagboltz gas;
gas.LoadGasFile("isobutane.gas");
///////////////////////////////////////////////////////////////////////////////////////////////
// Make a component with analytic electric field.
ComponentAnalyticField cmp;
cmp.SetMedium(&gas); // Attach gas to cmp
// Set parameters
const int nElec = 19; // Number of electrodes - assumes CACACACA.... //(19)
const int nWires = 50; // Wires per electrode //(50)
const double dWire = 0.0007*2.54; // Diameter of each wire [cm]
const double elecspac = 0.497 * 2.54; // Spacing between electrodes [cm]
const double wirespac = 0.2; // Spacing between wires [cm]
//Only have one of each Cathode and Anode not commented out at a time
//const double vCathode = -250.; // Potential of cathodes [V]
//const double vAnode = 250.; // Potential of anodes [V]
const double vCathode = 0.; // Potential of cathodes [V]
const double vAnode = 500.; // Potential of anodes [V] //!! 500ish?
const double chamberRad = (nWires * wirespac / 2.) +2.5; // Radius of chamber [cm] //!!
double vWire;
std::string label;
// Create electrodes
bool cathode = true; // First electrode is a cathode
for (int j = 0; j < nElec; j++) {
// Set y positions
double ywire = elecspac * j;
// Set potential and label
if (cathode) {
vWire = vCathode;
label = "C"; // Elements with same label are grouped together
}
else if (!cathode) {
vWire = vAnode;
label = "A";
}
// Create wires
for (int i = -nWires / 2; i < nWires / 2; i++) {
double xwire = wirespac * i + wirespac/2.;
cmp.AddWire(xwire, ywire, dWire, vWire, label);
}
// Set 'cathode' for next electrode
cathode = !cathode; // Alternate CACACA...
}
// Define area (for sensor and outer casing)
const double xmin = -chamberRad;
const double xmax = chamberRad;
const double ymin = -elecspac; // !! Should actually be size of dead region
const double ymax = nElec * elecspac; // !!?
const double zmin = -chamberRad;
const double zmax = chamberRad;
const unsigned int nx = 1000; //Should give point density, Default is 200
const unsigned int ny = 1000; //Should give point density, Default is 200
// Create outer casing - model w 4 planes
cmp.AddPlaneX(xmin, 0., "P1");
cmp.AddPlaneY(ymin, 0., "P2");
cmp.AddPlaneX(xmax, 0., "P3");
cmp.AddPlaneY(ymax, 0., "P4");
// Request calculation of the weighting field.
cmp.AddReadout("A"); // Only anodes read out
///////////////////////////////////////////////////////////////////////////////////////////////
// Make a sensor.
Sensor sensor;
sensor.AddComponent(&cmp); // Signal from cmp
// Set area // Area the sensor considers
sensor.SetArea(xmin, ymin, zmin, xmax, ymax, zmax);
sensor.ClearSignal();
///////////////////////////////////////////////////////////////////////////////////////////////
// Set up SRIM
TrackSrim tr;
tr.SetSensor(&sensor);
const std::string file = "Aluminum-26 in No. 340 Isobutane (gas).txt";
tr.ReadFile(file);
double ekin = 4.e6;
tr.SetKineticEnergy(ekin);
///////////////////////////////////////////////////////////////////////////////////////////////
// Set up drift lines
std::cout << "RKF integration\n";
DriftLineRKF drift;
drift.SetMaximumStepSize(0.1 * wirespac); // Setting Max Step Size; Experimental to get rid of "Hooks"
drift.SetSensor(&sensor); // Attach to the sensor
///////////////////////////////////////////////////////////////////////////////////////////////
// Set canvas
std::cout << "Set canvas for drift\n";
TCanvas* cD = nullptr; // Canvas cD will be used to plot the electron drift lines
cD = new TCanvas("cD", "", 600, 600); // Sets stuff up but won't be plotted yet
//SetNumberOfSamples2d(nx, ny); // Set number of points, default 200
// Defined outside loop to allow it to be used after loop
ViewCell cellView; // To visualise outline of cell i.e. tube, wire etc
cellView.SetComponent(&cmp); // Viewing cmp
ViewDrift driftView; // To visualise electron drift lines
driftView.SetCanvas(cD); // Set canvas to cD
driftView.SetArea(xmin, ymin, xmax, ymax); // Restrict area to be plotted
drift.EnablePlotting(&driftView); // Also need to specifically enable plotting drift
tr.EnablePlotting(&driftView); // Test, SRIM plotting
cellView.SetCanvas(cD); // cellView will be plotted on same canvas
cellView.SetArea(xmin, ymin, xmax, ymax);
///////////////////////////////////////////////////////////////////////////////////////////////
// Open output file
std::ofstream output;
output.open( "out.txt" );
output << "# xi yi wire time \n";
///////////////////////////////////////////////////////////////////////////////////////////////
// Main loop
int n1 = 2; // 10 Electrons per row
//int n1 = 100; // 100 Electrons per row
int n2 = 2;//100; // Number of rows - "10"
//int n2 = 100; // Number of rows - "100"
for (int j = 0; j < n2; ++j) { // Each row //~20mins per row of 200 electrons
for (int k = 0; k < n1+1; ++k) { // Each electron in the row
// Set initial positions
//double xei = double(rand()) / double(RAND_MAX) * (xmax - xmin) + xmin;
//double yei = double(rand()) / double(RAND_MAX) * (ymax - ymin) + ymin;
//double zei = double(rand()) / double(RAND_MAX) * (zmax - zmin) + zmin;
//double xei = k * 2. * chamberRad / n1 + xmin;
//double yei = j * (ymax - ymin) / n2 + ymin; //(j * 0.1 + 0.05) * elecspac;
//double yei = j * (ymax - ymin) / n2 + (3 * ymin / 2); //(j * 0.1 + 0.05) * elecspac for adjusted middle positions
//double zei = 0.;
//double tei = 0.;
//std::cout << xei << " " << yei << " " << zei << " " << tei << "\n";
// Simulate a track.
//const double xt = xmin;
const double xt = ((xmin + xmax) / 2) + (10*j + k) * (xmax - xmin) / 100; //y-track
//const double yt = ymin + (10*j + k) * (ymax - ymin) / (n1 * n2);
const double yt = ymin; //y-track
//track.NewTrack(xt, yt, 0., 0., 1., 0., 0.);
tr.NewTrack(xt, yt, 0., 0., 0., 1., 0.); //y-track;
////////////////////////////////////////////////////////////////////////////////
//Test of Track Output
// Loop over the clusters.
double xc, yc, zc, tc, ec, extra;
int nc;
while (tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
for (int j = 0; j < nc; ++j) {
tr.GetCluster(xc, yc, zc, tc, nc, ec, ekin);
}
}
}
output << "\n \n";
}
///////////////////////////////////////////////////////////////////////////////
// Close output file
output.close();
// Plot drift
std::cout << "Plot drift\n";
cD->Clear();
cellView.Plot2d(); // Plot2d to plot cell (set up earlier)
constexpr bool twod = true;
constexpr bool drawaxis = false;
driftView.Plot(twod, drawaxis); // Plot to plot drift lines (set up earlier) - NOTE: when drawaxis is true, cell gets overwritten... FIX(?)
///////////////////////////////////////////////////////////////////////////////////////////////
// Set up fieldView to plot the electric potential
std::cout << "Set up fieldView\n";
TCanvas* cF = nullptr; // Canvas cF will be used to plot potential/field
ViewField fieldView; // ViewField object manages this
cF = new TCanvas("cF", "", 700, 600); // Doesn't rely on the main simulation loop so we can set up then plot immediately
fieldView.SetCanvas(cF);
fieldView.SetSensor(&sensor); // ViewField gets the field/potentials from sensor (sensor is an interface between a lot of things)
fieldView.SetPlane(0., 0., 1., 0., 0., 0.); // Set plane to plot - xy plane
cellView.SetCanvas(cF); // Also plot cell on same plot (this is why cellView was set up BEFORE main loop)
cellView.SetPlane(0., 0., 1., 0., 0., 0.); // Set plane to plot - xy plane
cellView.SetArea(xmin, ymin, xmax, ymax); // Restrict area (so that it lines up with area set by sensor)
// Plot field
std::cout << "Plot field\n";
//fieldView.Plot("e","CONT2"); // ARR looks promising but too many cells
fieldView.PlotContour("v"); // Plot "e" for field magnitude, "v" for potential (see table 4.2, pg55)
cellView.Plot2d(); // Also plot cell
///////////////////////////////////////////////////////////////////////////////////////////////
app.Run(kTRUE);
}