Segmentation Violation When Using DriftRKF::DriftElectron()

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, 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]
  //pressure [Torr]
  //initialize table of scattering rates (called internally when a collision rate is requested and the gas mixture or other parameters have changed)

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

  //add electrodes
  cmp.AddPlaneY(yAnode, vAnode, "anode");
  cmp.AddPlaneY(yCathode, vCathode, "cathode");
  Sensor sensor;
  //add electrode that will measure induced current
  sensor.AddElectrode(&cmp, "anode");
  //set drift region
  //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;

  //enable plotting drift lines
  ViewDrift 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) {

  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.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
  driftView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
  driftView.Plot(true, false);
  //plot induced current on anode
  ViewSignal signalView;
  TCanvas inducedCurrent("inducedCurrent", "Induced Current over Time", 800, 600);
  //plot electric field
  TCanvas* canvas = new TCanvas("c1", "Field and Cell View", 800, 600);
  ViewField fieldView;
  fieldView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);
  fieldView.SetVoltageRange(vCathode, vAnode);

  ViewCell cellView;
  cellView.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);

  //plot electric field lines
  TCanvas* fieldLinesCanvas = new TCanvas("fieldLinesCanvas", "Electric Field Lines", 800, 600);

  ViewField fieldLinesView;
  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++) {
  fieldLinesView.PlotFieldLines(x0, y0, z0, true, true, kOrange-3);

  ViewCell cellView2;
  cellView2.SetArea(-0.5, -0.5, ringThickness + ringInnerDiam + 0.5, yCathode + 0.5);


DriftLineRKF needs a table of the macroscopic drift velocity as function of the electric field. For ions, this is what you import when calling LoadIonMobility. For electron, you’ll need to load a gas table using LoadGasFile (pre-computed using Magboltz).

Thank you! That solved the issue.

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