Incorporating TRIM into Garfield Simulations

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

}

Hi,
in the program you sent there is a TrackSrim object created, but then it’s not used later on in the program. Instead the program uses TrackHeed (with default settings, i. e. a minimum-ionizing muon) to simulate tracks. Is that intentional?
To actually use TrackSrim, you would need to replace track by tr, and remove the GetElectronloop (with TrackSrim the coordinates of the electrons are identical to the “cluster” ones).

PS: I’ve seen your mail with the TRIM output file, but I’ve not yet had time to look into it in more detail.

Thank you for responding.
Yes, this code is the last version made that at least somewhat works. If that’s the case with Clusters and Elections in TrackSrim then I’ll be able to fix that at least. Thank you for that.

Understood about the TRIM file, I get that real life happens. I do want to be able to use it at some point but if this were easy anyone could do it.

Edited the Code to show what it looks like now. It looks better, I’ve picked up a logic error when it comes to actually showing the results but I’m certain it’s going through the SRIM properly otherwise now.

Still not complete, but feels like a big step forward has been done and thank you for that.

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