Howdy.
I’m back. So with the help of the last thread I made a lot of progress and think I’m just one bug away from getting it how I want.
What this program does is model the results of an 26-Al ion beam hitting a helium target and tracking the resulting 29-Si + p reaction results. All of those parts are in, I just also want to see the resulting Electron Drift while I’m at it.
The code that follows appears to do all that as well as it can, but I get a “TrackHeed GetElectron: Index Out of Range” error each time I try to run the loop when I run it. I know what an Index Out of Range error is, I don’t understand why I’m seeing it.
Any help would be appreciated.
#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/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/TrackElectron.hh"
#include "Garfield/TrackTrim.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 TRIM
TrackHeed etrack;
TrackTrim Al26tr, Si29tr, H1tr;
Al26tr.SetSensor(&sensor);
Si29tr.SetSensor(&sensor);
H1tr.SetSensor(&sensor);
etrack.SetSensor(&sensor);
//const std::string file = "Aluminum-26 in No. 340 Isobutane (gas).txt";
//tr.ReadFile(file);
const unsigned int nIons = 10000;
const unsigned int nSkip = 0;
if ((!Al26tr.ReadFile("Al26EXYZ.txt", nIons, nSkip)) || (!Si29tr.ReadFile("Si29EXYZ.txt", nIons, nSkip)) || (!H1tr.ReadFile("H1EXYZ.txt", nIons, nSkip))) {
std::cerr << "Reading TRIM EXYZ files failed.\n";
return 1;
}
std::vector<double> AlSelect = { 15.696e6 , 13.362e6 , 34.725e6 , 32.583e6 };
/*
Kinetic Energy -
Ecm = 2.5 MeV , 0.72 MeV/u -> T = 18.73e6
Ecm = 5 MeV , 1.441 MeV/u -> T = 37.46e6
0.72 MeV/u, 1.12 um -> el = 0.162 -> AlEkin = 15.696e6
1.441 MeV/u, 1.12 um -> el = 0.073 -> AlEkin = 34.725e6
0.72 MeV/u, 2 um -> el = 0.287 -> AlEkin = 13.362e6
1.441 MeV/u, 2 um -> el = 0.130 -> AlEkin = 32.583e6
This is a crude approach, but good enough for now.
*/
int a, b, c;
std::cout << "Select Beam Energy - 1.) 0.72 MeV/u or 2.) 1.441 MeV/u: ";
std::cin >> a ;
std::cout << "Select Mylar Thickness - 1.) 1.12 um or 2.) 2 um: ";
std::cin >> b ;
double AlEkin, SiEkin, HEkin;
if ( (a == 1) && (b == 1) ) AlEkin = AlSelect[0];
if ( (a == 1) && (b == 2) ) AlEkin = AlSelect[1];
if ( (a == 2) && (b == 1) ) AlEkin = AlSelect[2];
if ( (a == 2) && (b == 2) ) AlEkin = AlSelect[3];
Al26tr.SetWorkFunction(30); // Work Function, just set to 30 eV
Si29tr.SetWorkFunction(30);
H1tr.SetWorkFunction(30);
///////////////////////////////////////////////////////////////////////////////////////////////
// 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
drift.EnableSignalCalculation();
///////////////////////////////////////////////////////////////////////////////////////////////
// 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
// 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 driftViewAl, driftViewSi, driftViewH, driftViewe; // To visualise electron drift lines
driftViewAl.SetCanvas(cD); // Set canvas to cD
driftViewAl.SetArea(xmin, ymin, xmax, ymax); // Restrict area to be plotted
driftViewSi.SetCanvas(cD); // Set canvas to cD
driftViewSi.SetArea(xmin, ymin, xmax, ymax); // Restrict area to be plotted
driftViewH.SetCanvas(cD);
driftViewH.SetArea(xmin, ymin, xmax, ymax);
driftViewe.SetCanvas(cD);
driftViewe.SetArea(xmin, ymin, xmax, ymax);
drift.EnablePlotting(&driftViewAl); // Also need to specifically enable plotting drift
drift.EnablePlotting(&driftViewSi);
drift.EnablePlotting(&driftViewH);
drift.EnablePlotting(&driftViewe);
Al26tr.EnablePlotting(&driftViewAl); // TRIM plotting
Si29tr.EnablePlotting(&driftViewSi);
H1tr.EnablePlotting(&driftViewH);
etrack.EnablePlotting(&driftViewe);
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";
///////////////////////////////////////////////////////////////////////////////////////////////
std::cout << "Select Si-Excitation Energy - 1.) 0 MeV, 2.) 4 MeV, or 3.) 8 MeV: ";
std::cin >> c ;
std::vector<float> SiAngToE, HAngToE;
//Silicon
if ( (AlEkin == AlSelect[0]) && (c == 1) ) SiAngToE.insert( SiAngToE.end(), { 0.089, 0.181, 16.485, 7.02 });
else if ( (AlEkin == AlSelect[0]) && (c == 2) ) SiAngToE.insert( SiAngToE.end(), { 0.116, 0.140, 15.15, 4.59 });
else if ( (AlEkin == AlSelect[0]) && (c == 3) ){
std::cout << "System Kinematics Not Possible" ;
return 0;
}
else if ( (AlEkin == AlSelect[1]) && (c == 1) ) SiAngToE.insert( SiAngToE.end(), { 0.086, 0.171, 13.98, 7.02 });
else if ( (AlEkin == AlSelect[1]) && (c == 2) ) SiAngToE.insert( SiAngToE.end(), { 0.111, 0.131, 12.67, 4.59 });
else if ( (AlEkin == AlSelect[1]) && (c == 3) ){
std::cout << "System Kinematics Not Possible" ;
return 0;
}
else if ( (AlEkin == AlSelect[2]) && (c == 1) ) SiAngToE.insert( SiAngToE.end(), { 0.214, 0.342, 34.94, 5.69 });
else if ( (AlEkin == AlSelect[2]) && (c == 2) ) SiAngToE.insert( SiAngToE.end(), { 0.251, 0.275, 33.35, 4.26 });
else if ( (AlEkin == AlSelect[2]) && (c == 3) ) SiAngToE.insert( SiAngToE.end(), { 0.314, 0.072, 31.08, 1.86 });
else if ( (AlEkin == AlSelect[3]) && (c == 1) ) SiAngToE.insert( SiAngToE.end(), { 0.216, 0.346, 32.51, 5.69 });
else if ( (AlEkin == AlSelect[3]) && (c == 2) ) SiAngToE.insert( SiAngToE.end(), { 0.254, 0.278, 30.90, 4.26 });
else if ( (AlEkin == AlSelect[3]) && (c == 3) ) SiAngToE.insert( SiAngToE.end(), { 0.317, 0.072, 28.60, 1.86 });
else SiAngToE.insert( SiAngToE.end(), { 0.00, 0.00, 0.00, 0.00});
//Proton
if ( (AlEkin == AlSelect[0]) && (c == 1) ) HAngToE.insert( HAngToE.end(), { 0.000609, 0.00801, 8.416, 18.06 });
else if ( (AlEkin == AlSelect[0]) && (c == 2) ) HAngToE.insert( HAngToE.end(), { 0.000498, 0.00232, 3.918, 18.06 });
else if ( (AlEkin == AlSelect[0]) && (c == 3) ){
std::cout << "System Kinematics Not Possible" ;
return 0;
}
else if ( (AlEkin == AlSelect[1]) && (c == 1) ) HAngToE.insert( HAngToE.end(), { 0.000504, 0.00815, 6.142, 18.06 });
else if ( (AlEkin == AlSelect[1]) && (c == 2) ) HAngToE.insert( HAngToE.end(), { 0.000363, 0.00193, 2.421, 18.06 });
else if ( (AlEkin == AlSelect[1]) && (c == 3) ){
std::cout << "System Kinematics Not Possible" ;
return 0;
}
else if ( (AlEkin == AlSelect[2]) && (c == 1) ) HAngToE.insert( HAngToE.end(), { 0.00141, 0.00606, 13.567, 18.06 });
else if ( (AlEkin == AlSelect[2]) && (c == 2) ) HAngToE.insert( HAngToE.end(), { 0.00110, 0.00694, 8.516, 18.06 });
else if ( (AlEkin == AlSelect[2]) && (c == 3) ) HAngToE.insert( HAngToE.end(), { 0.000850, 0.00193, 3.314, 18.06 });
else if ( (AlEkin == AlSelect[3]) && (c == 1) ) HAngToE.insert( HAngToE.end(), { 0.00133, 0.00649, 11.06, 18.06 });
else if ( (AlEkin == AlSelect[3]) && (c == 2) ) HAngToE.insert( HAngToE.end(), { 0.000927, 0.00725, 6.230, 18.06 });
else if ( (AlEkin == AlSelect[3]) && (c == 3) ) HAngToE.insert( HAngToE.end(), { 0.000606, 0.00159, 1.973, 18.06 });
else HAngToE.insert( HAngToE.end(), { 0.00, 0.00, 0.00, 0.00});
double thetaSi, thetaH;
///////////////////////////////////////////////////////////////////
//Silicon-29 Ion Angular Spread and Chamber Energy
//SiEkin = c + b*theta - a*theta*theta, theta_limit = d
//(a, b, c, d) calculated using LISE++ and Catkin
//Chamber maxes at theta = 18.06 degrees
//////////////////////////////////////////////////////////////////
Al26tr.SetKineticEnergy(AlEkin);
// Main loop
int n1 = 100; //100 Electrons per row
int n2 = 100; //100 rows
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
// Simulate a track.
const double xt = ((xmin + xmax) / 2) + (100*j + k) * (xmax - xmin) / 2500000; //y-track; 2500000 goes to 1 mm width
const double yt = ymin; //y-track
Al26tr.NewTrack(xt, yt, 0., 0., 0., 1., 0.); //y-track;
thetaSi = -1 * SiAngToE[3] + (100*j + k) * SiAngToE[3] / 5000;
SiEkin = (SiAngToE[2] + SiAngToE[1]*abs(thetaSi) - SiAngToE[0]*thetaSi*thetaSi) * 1e6;
Si29tr.SetKineticEnergy(SiEkin);
Si29tr.NewTrack(xt + 23*tan( thetaSi * M_PI / 180), yt, 0., 0., sin( thetaSi * M_PI / 180),cos( thetaSi * M_PI / 180), 0.);
////////////////////////////////////////////////////////////////////////////////
//Test of Track Output
//Loop over the clusters.
double xc, yc, zc, tc, ec, extra;
int nc;
double ee, dx, dy, dz;
if ( (j % 10) == 1) {
thetaH = -1 * HAngToE[3] + (100*j) * HAngToE[3] / 5000; //got rid of k and added an if loop to clean up the view
HEkin = (HAngToE[2] - HAngToE[1]*abs(thetaH) - HAngToE[0]*thetaH*thetaH) * 1e6;
H1tr.SetKineticEnergy(HEkin);
H1tr.NewTrack(xt + 23*tan( thetaH * M_PI / 180), yt, 0., 0., sin( thetaH * M_PI / 180), cos( thetaH * M_PI / 180), 0.);
while (Al26tr.GetCluster(xc, yc, zc, tc, nc, ec, extra) && Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, extra) && H1tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
etrack.NewTrack(xc, yc, zc, tc, 0, 0, 0);
etrack.GetElectron(j, xc, yc, zc, tc, ee, dx, dy, dz);
drift.DriftElectron(xc, yc, zc, tc);
}
}
if ( ( j % 10) != 1){
while (Al26tr.GetCluster(xc, yc, zc, tc, nc, ec, extra) && Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
etrack.NewTrack(xc, yc, zc, tc, 0, 0, 0);
etrack.GetElectron(j, xc, yc, zc, tc, ee, dx, dy, dz);
drift.DriftElectron(xc, yc, zc, tc);
}
}
output << "\n \n";
}
std::cout << j << "%" << std::endl;
}
driftViewSi.SetColourTracks(kViolet);
driftViewH.SetColourTracks(kBlue);
driftViewe.SetColourTracks(kOrange);
///////////////////////////////////////////////////////////////////////////////
// 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;
driftViewAl.Plot(twod, drawaxis); // Plot to plot drift lines (set up earlier) - NOTE: when drawaxis is true, cell gets overwritten... FIX(?)
driftViewSi.Plot(twod, drawaxis);
driftViewH.Plot(twod, drawaxis);
driftViewe.Plot(twod, drawaxis); //Test
///////////////////////////////////////////////////////////////////////////////////////////////
// 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);
}