Home | News | Documentation | Download

TrackHeed GetElectron: Index Out Of Range

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);
}

Hi,
there is probably a misunderstanding; I don’t think you need TrackHeed for what you want to do. To simulate the drift lines of the electrons produced by the 29-Si track, you can just do something like this:

Si29tr.NewTrack(xt, yt, zt, 0., sin(thetaSi * M_PI / 180), cos(thetaSi * M_PI / 180), 0.);
while (Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
  drift.SetElectronSignalScalingFactor(nc);
  drift.DriftElectron(xc, yc, zc, tc);
}

Also reading the code, I’m not sure if the while loop

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)) {
  // ..
}

is really what you intended. In the body of the loop, xc, yc, zc, tc, nc will have the values returned by H1tr.GetCluster(...), and the clusters from the Al and Si tracks will be “lost”.

Well I tried what you said, commenting out the unneeded code and making the adjustments as you suggested.

It does compile and seems to run, but I ran it for over an hour and a half and it wasn’t done. I suspect a memory leak or a logic error that loops everything to much somewhere. Here’s the resulting code.

#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);

      double xc, yc, zc, tc, ec, extra;
      int nc;

      // 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;
	  while (Al26tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
	    drift.SetElectronSignalScalingFactor(nc);
	    drift.DriftElectron(xc, yc, zc, tc);
	  }
	  
	  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.);
	  while (Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
	    drift.SetElectronSignalScalingFactor(nc);
	    drift.DriftElectron(xc, yc, zc, tc);
	  }
          ////////////////////////////////////////////////////////////////////////////////           
          //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 (H1tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
	      drift.SetElectronSignalScalingFactor(nc);
	      drift.DriftElectron(xc, yc, zc, tc);
	    }
	    /*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);
}

Just cleaned up the code a little bit. Like I said, it looks like it compiles and run okay but I’ve let it run for hours and it doesn’t finish.

Is there a leak somewhere that I’m missing? Or am I just not being patient enough?

#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);
    
    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;                        // 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);
    
    drift.EnablePlotting(&driftViewAl);             // Also need to specifically enable plotting drift
    drift.EnablePlotting(&driftViewSi);
    drift.EnablePlotting(&driftViewH);
    Al26tr.EnablePlotting(&driftViewAl);            // TRIM plotting 
    Si29tr.EnablePlotting(&driftViewSi);
    H1tr.EnablePlotting(&driftViewH);
    
    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);

      double xc, yc, zc, tc, ec, extra;
      int nc;

      // 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;
	  while (Al26tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
	    drift.SetElectronSignalScalingFactor(nc);
	    drift.DriftElectron(xc, yc, zc, tc);
	  }
	  
	  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.);
	  while (Si29tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
	    drift.SetElectronSignalScalingFactor(nc);
	    drift.DriftElectron(xc, yc, zc, tc);
	  }
          ////////////////////////////////////////////////////////////////////////////////           
          //Test of Track Output                                                               
          //Loop over the clusters.                                                            
	  
	  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 (H1tr.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
	      drift.SetElectronSignalScalingFactor(nc);
	      drift.DriftElectron(xc, yc, zc, tc);
	    }
	    
	  output << "\n \n";
	  }
	}
	std::cout << j << "%" << std::endl;
      }
      
      driftViewSi.SetColourTracks(kViolet);
      driftViewH.SetColourTracks(kBlue);
      ///////////////////////////////////////////////////////////////////////////////
      // 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);
      
    ///////////////////////////////////////////////////////////////////////////////////////////////
    // 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);
}

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.