Hi there,
I am trying to simulate a double gridded drift cell by using DriftRKF to drift electrons from a cathode towards an anode using forces caused by an electric field, which is generated by field shaping rings that the electrons will drift through. I am getting a segmentation violation when I use DriftRKF::DriftElectron(), but not when I use DriftIon(). Something that has been confusing me about DriftIon() though is that I get status codes for drift lines such as 67, 69, or 71, which I could not find defined in GarfieldConstants.hh, so was also wondering where I could find what these numbers mean.
Before I get the segmentation violation, I get the error: “DriftLineRKF::GetVelocity: Cannot retrieve drift velocity at (x, y, z)” which makes me wonder if the two issues are related, however I had that error previously when using DriftIon(), but without the segmentation violation.
The terminal says that the error is at line 165 of DriftLineRKF.cc, which uses the std::vector::back() function. Therefore, I think that the error could be that the vector which calls the back() function is empty, or something similar. Attached below is my code. I am currently using gdb to try to debug the issue - any pointers/advice would be greatly appreciated.
#include <iostream>
#include <fstream>
#include <TApplication.h>
#include <TCanvas.h>
#include <TH1F.h>
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/Random.hh"
using namespace Garfield;
int main(int argc, char * argv[]) {
TApplication app("app", &argc, argv);
//variables to change
//grid pitch (30 lpi or 70 lpi)
constexpr bool lpi30 = false;
//number of rings
constexpr int numRings = 20;
//drift length [cm]
constexpr double driftLength = 5.;
//inner diameter of the rings [cm]
constexpr double ringInnerDiam = 2.5;
//thickness of the rings [cm]
constexpr double ringThickness = 0.1;
//gap between the grid and the first ring [cm]
constexpr double egGap = 0.162;
// wire diameter [cm]
constexpr double gridDiam = 0.0035;
//electric field between the anode grid and cathode grid [V/cm]
constexpr double eGG = 400;
//electric field between the anode and the grid [V/cm]
constexpr double eAG = 800;
//electric field between the cathode and the grid [V/cm]
constexpr double eCG = 200;
//grid spacing [cm]
constexpr double gridGap = lpi30 ? 0.0846: 0.0363;
// number of grid wires
constexpr int numGridWire = (ringInnerDiam-ringThickness)/gridGap;
//spacing between the rings
constexpr double ringGap = driftLength/(numRings+1);
// voltage settings [V], lower grid set at ground
constexpr double vAnode = eAG*egGap;
constexpr double vLowerGrid = 0;
constexpr double vUpperGrid = -(driftLength*eGG);
constexpr double vCathode = vUpperGrid - eCG*egGap;
// yGrid placement [cm]
constexpr double yLowerGrid = egGap;
constexpr double yUpperGrid = driftLength + egGap;
//electrode placement [cm]
constexpr double yCathode = yUpperGrid + egGap;
constexpr double yAnode = 0.;
//setup the gas
MediumMagboltz gas;
//set the temperature [K] and pressure [Torr].
//100% xenon gas
gas.SetComposition("xe", 100.);
//temperature [K]
gas.SetTemperature(293.15);
//pressure [Torr]
gas.SetPressure(15*760.);
//initialize table of scattering rates (called internally when a collision rate is requested and the gas mixture or other parameters have changed)
gas.Initialise(false);
//number of electrons produced in initial avalanche
int initElectronTotal = 0;
//number of electrons that end up at the anode
int finalElectronTotal = 0;
//load the ion mobilities.
const std::string path = std::getenv("GARFIELD_INSTALL");
gas.LoadIonMobility(path + "/share/Garfield/Data/IonMobility_Xe+_P12_Xe.txt");
//setup field
ComponentAnalyticField cmp;
cmp.SetMedium(&gas);
//add electrodes
cmp.AddPlaneY(yAnode, vAnode, "anode");
cmp.AddPlaneY(yCathode, vCathode, "cathode");
Sensor sensor;
sensor.AddComponent(&cmp);
//add electrode that will measure induced current
sensor.AddElectrode(&cmp, "anode");
//set drift region
sensor.SetArea();
//setup binning for signal calculation [ns]
sensor.SetTimeWindow(0., 1., 100);
//add grids
for (int i=0; i<numGridWire; i++) {
cmp.AddWire(0+i*gridGap, yLowerGrid, gridDiam, vLowerGrid, "lowergrid");
cmp.AddWire(0+i*gridGap, yUpperGrid, gridDiam, vUpperGrid, "uppergrid");
}
//add rings
for (int i=0; i<numRings; i++) {
std::ostringstream oss;
oss << "ring" << i+1;
std::string label = oss.str();
double v =(i+1)*(vUpperGrid/(numRings+1));
cmp.AddWire(-ringThickness/2, yLowerGrid+(i+1)*ringGap, ringThickness, v, label);
cmp.AddWire(ringInnerDiam+ringThickness/2, yLowerGrid+(i+1)*ringGap, ringThickness, v, label);
}
//drift line rkf to drift electrons -> calculates path of electron or ion by numerical integration of drift velocity vector, well adapted to smooth fields (e.g. analytic potentials)
DriftLineRKF driftline;
driftline.SetSensor(&sensor);
//enable plotting drift lines
ViewDrift driftView;
driftline.EnablePlotting(&driftView);
//below: using DriftRKF to drift electrons, causes segmentation violation
const int nIons = 10;
int plane = 0;
for (int i = 0; i<nIons; i++) {
constexpr double r = 1.25;
const double angle = 0.05;
const double startX = RndmGaussian(1.25, 0.05);
const double startY = yCathode;
driftline.DriftElectron(startX, yCathode, 0, 0);
double x1, y1, z1, t1;
int status = 0;
driftline.GetEndPoint(x1, y1, z1, t1, status);
std::cout<<"("<<x1<<", "<<y1<<")"<<" status: "<<status<<std::endl;;
if (y1<=1.62&&x1>=0&&x1<ringInnerDiam) {
++plane;
}
}
float gridTransparency = static_cast<float>(plane)/nIons;
std::cout<<"Total initial electrons: "<<nIons<<std::endl;
std::cout<<"Total final electrons: "<<plane<<std::endl;
std::cout<<"Grid transparency: "<<gridTransparency<<std::endl;
ViewCell cellView0;
cellView0.SetComponent(&cmp);
cellView0.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
cellView0.EnableWireMarkers(false);
cellView0.Plot2d();
driftView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
driftView.SetCanvas(cellView0.GetCanvas());
driftView.Plot(true, false);
//plot induced current on anode
ViewSignal signalView;
signalView.SetSensor(&sensor);
TCanvas inducedCurrent("inducedCurrent", "Induced Current over Time", 800, 600);
signalView.SetCanvas(&inducedCurrent);
signalView.PlotSignal("cathode");
//plot electric field
TCanvas* canvas = new TCanvas("c1", "Field and Cell View", 800, 600);
ViewField fieldView;
fieldView.SetCanvas(canvas);
fieldView.SetComponent(&cmp);
fieldView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
fieldView.SetVoltageRange(vCathode, vAnode);
fieldView.PlotContour();
ViewCell cellView;
cellView.SetCanvas(canvas);
cellView.SetComponent(&cmp);
cellView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
cellView.EnableWireMarkers(false);
cellView.Plot2d();
//plot electric field lines
TCanvas* fieldLinesCanvas = new TCanvas("fieldLinesCanvas", "Electric Field Lines", 800, 600);
ViewField fieldLinesView;
fieldLinesView.SetCanvas(fieldLinesCanvas);
fieldLinesView.SetComponent(&cmp);
fieldLinesView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
constexpr int numFieldLines = 75;
constexpr double fieldGap = ringInnerDiam/numFieldLines;
std::vector<double> x0;
std::vector<double> y0;
std::vector<double> z0;
for (int i=0; i<numFieldLines; i++) {
x0.push_back(0+i*fieldGap);
y0.push_back(yCathode);
z0.push_back(0.0);
}
fieldLinesView.PlotFieldLines(x0, y0, z0, true, true, kOrange-3);
ViewCell cellView2;
cellView2.SetCanvas(fieldLinesCanvas);
cellView2.SetComponent(&cmp);
cellView2.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
cellView2.EnableWireMarkers(false);
cellView2.Plot2d();
app.Run(true);
}