DriftLineRKF lines ending inside a SolidWire


an issue with the electron drift simulation in Garfield++ has been plaguing me for the past few days.

I have been trying to simulate the passage of a muon track inside a drift chamber cell that I had set up:

  1. I defined the structure of the cell using the native geometry (SolidBox for plates and SolidWire for wires) and computed the electric field with ComponentNeBem3d.
  2. I then aimed at getting the signal induced on the central wire by a particle track, computing the drift lines with the DriftLineRKF class and plotting the induced signal and the simulated drift lines.

This is the code that I have been using:

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <sys/stat.h>
#include <math.h>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>

#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/SolidWire.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"

using namespace Garfield;

int main(int argc, char * argv[]){
    TApplication app("app", &argc, argv);

    //load the gas medium for the internal space
    MediumMagboltz gas;
    //define the conductor material for the strips and wire
    MediumConductor plane_medium;
    MediumConductor strip_medium;
    MediumConductor wire_medium;

    //Defining the finite plane geometry through a GeometrySimple object
    GeometrySimple geo;

    //geometric parameters
    //----whole cell----
    //half-length of the assembly along z [cm]
    const double hlength = 100.;
    //gap half-length
    const double half_gap = 0.1;    
    //half width of a cell
    const double cell_width = 1.;
    //half thickness of a cell
    const double cell_thick = 0.5*cell_width;

    //half dimensions
    const double stripWidth = 1*cell_width-half_gap;
    const double stripThick = 35e-4;

    //signal wire radius [cm]
    const double r_signal_wire=20.e-4/2;
    //field wire radius
    const double r_field_wire=100.e-4/2;

    //strip voltage
    const double vStrip = -1000.;
    //signal wire voltage
    const double v_signal_wire = +1500.;
    //field wire voltage
    const double v_field_wire = -3500.;

    //define the cathode strips
    SolidBox c_strip_up(0.,cell_thick-stripThick,0.,stripWidth,stripThick,hlength);
    SolidBox c_strip_down(0.,-cell_thick+stripThick,0.,stripWidth,stripThick,hlength);
    //define the signal wire
    SolidWire sens_wire(0.,0.,0.,r_signal_wire,hlength);
    //set this wire as the electrode
    //define the field wires: left side
    SolidWire field_wire_u_l(-cell_width,0.,0.,r_field_wire,hlength);
    //define the field wires: right side
    SolidWire field_wire_u_r(+cell_width,0.,0.,r_field_wire,hlength);
    //set the boundary potentials for the components

    //add wires strips and planes to the geometry
    geo.AddSolid(&c_strip_up, &strip_medium);
    geo.AddSolid(&c_strip_down, &strip_medium);

    //add the boundary volume filled with gas
    //Before moving on one needs to compute the field with neBEM
    ComponentNeBem3d nebem;
    //Use a Sensor object to interface to the transport classes
    Sensor sensor;
    //add the field and electrode components

    const double tstep = 0.5;
    const double tmin = 0.;
    const unsigned int nbins = 100.;
    //setup the binning for the signal calculation (start, bin width[ns],bin number)

    //Simulating the ionization with Heed
    //first a track object is defined
    std::cout<<"-- Particle track simulation --"<<std::endl;
    TrackHeed track;
    //compute the track
    //compute the electron drift lines
    DriftLineRKF drift;

    //set the initial position: middle-left side of the cell
    const double x0 = -cell_width/2, y0 = -cell_thick, z0 = 0., t0 = 0.;
    //set the track direction and define the track
    const double dx0 = 0., dy0 = 1., dz0 = 0.;

    //loop over the clusters
    double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
    int nc =0;
        //loop over the electrons in the cluster
        for(int k = 0; k<nc; ++k){
            //electron data vars. are defined
            double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
            double dx = 0., dy = 0., dz = 0.;
            double enx = 0., eny = 0., enz = 0., ent = 0.;
            int st = 0.;
            track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
            drift.DriftElectron(xe, ye, ze, te);
            drift.GetEndPoint(enx, eny, enz, ent, st);
            std::cout<<"x: "<<enx<<" y: "<<eny<<" z: "<<enz<<" t: "<<ent<<" stat: "<<st<<std::endl; 

    std::cout<<"Induced charge: "<<sensor.GetTotalInducedCharge("s")<<std::endl;

    //plot the induced current on the wire by the electrons with ViewSignal
    ViewSignal signalView;

    //Setting up the viewer of the simulated drift lines
    //create a new canvas
    TCanvas* cD = new TCanvas("cD"," ",1100,500);
    //the ViewDrift object plots the cluster coordinates and drift line points
    ViewDrift driftView;
    //set the object with which to plot and the canvas

    //set the object with which to plot and the canvas
    ViewGeometry geoView;
    constexpr bool twod = true;
    constexpr bool drawaxis = false;

    std::cout<<"End of program"<<std::endl;

    return 0;

Upon execution I don’t seem to be getting any error by the track and drift simulation, but the induced charge on the central electrode remains 0 and no drift lines are plotted when using ViewDrift::Plot().
I inspected the drift computations with DriftLineRKF::GetEndPoint() and found that the endpoints are inside the central wire volume and the status flag is set to -5 (so to the LeftDriftMedium status, if I’m not mistaken):

    Warning: Integrating the Townsend coefficients would lead to exponential overflow.
    Avalanche truncated.
x: -0.000597484 y: -0.000801882 z: -7.56021e-05 t: 124.356 stat: -5

I have tried many combinations of the Sensor and DriftLineRKF options, but I haven’t been able to solve the issue. I also noticed this forum thread dealing with AvalancheMC instead:
Electrons generated inside anode wire during avalanche.

What could be the issue with my code?
Any feedback that you might have would be greatly appreciated.

I am attaching the .gas file that I have been using as a .txt file.
ar_85_co2_15_1atm.txt.zip (7.3 KB)

1 Like

Welcome to the ROOT forum,

I guess @hschindl can help you.

It’s a few things. First of all, sensor.SetTimeWindow(0.,1000.,100) is way too long. Next, you need to EnablePlotting before the main loop. After you change to sensor.SetTimeWindow(0.,1.,300) and move the whole ViewDrift block up you get this:

So far so good. Why does the warning happen? Because there is overflow in DriftLineRKF. This avalanche is too big for poor DriftLineRKF. Why is induced charge zero? Because DriftLineRKF doesn’t support calculating induced charges. You can hack the system by just integrating the plot with signalView.GetHistogram()->Integral() but results will probably be incorrect with truncated avalanches.

1 Like

@hschindl should know.

Thank you so much for the help @mjag!
The drift line plot is working just fine now: I had managed to miss the correct placement of the ViewDrift block until now.
As for the signal plot, I am still getting a flat-line after setting your reasonable time window. Is any other modification needed?

In general, truncated avalanches could be avoided by computing the gas table for sufficiently high values of the fields, right?

As for the signal plot, I am still getting a flat-line after setting your reasonable time window. Is any other modification needed?

Hmm, not really. Works for me™. Maybe you changed something else in your code? Here’s your example with the changes I told you:
test.C (6.3 KB)
If that doesn’t work then I don’t know, maybe some Garfield++ version mismatch or something.

In general, truncated avalanches could be avoided by computing the gas table for sufficiently high values of the fields, right?

No. The limit is hardcoded to e^{30} \approx 10^{13} particles. Note that this happens very close to the wire. Wire is so thin that electric field becomes too large. Not sure if it’s Geiger region yet, but it sure is too much for Garfield++. I can’t suggest you modify the actual design of the detector so that it works with Garfield++ (lol) but when I increased r_signal_wire twice (to 20.e-4) then I got no warnings. Also, if it is indeed Geiger region, then you can have nasty effects like secondary UV photon generation, so it’s not easily quantifiable.

Well, this does seem quite puzzling: I tried executing your code, but got the usual empty result. Even after building Garfield again from the repo (I cloned the master branch) the output is the same. Is your version different?

What seems strange to me is that this type of signal plot works in the case of the Drift Tube example Examples/DriftTube · master · garfield / garfieldpp · GitLab.

I updated Garfield++ to recent master and I still can’t recreate your problem - I get non-zero signal. I don’t know how else I could help you. Maybe you could try running it on another machine?

I removed neBEMOut and rerun program. After repeating this a few times I did manage to get empty plots. Maybe you have some old field cached in neBEMOut? Try removing that directory and rerunning.

Unfortunately removing neBEMOut didn’t work either: I get empty plots even just after clearing the entire build directory.
For the sake of completeness, I made a final check with Sensor::EnableDebugging() and got messages like these:

Sensor::AddSignal:   Time: 159.02
  Step: 0.334589
  Charge: -0.995868
  Velocity: (0.00350158, -0.00629439, 7.85088e-10)
  Electrode s:
    Weighting potentials: 0, 0
    Induced charge: -0

I will try installing Garfield++ locally on my Mac and update you on the results.
Thanks for the suggestions in any case!

An update on the issue:

  • The installation that I had used so far was on Centos 7.9 building the project from the HEAD master repository.
  • I installed the HEAD master repository on Ubuntu-22 (with gcc 11.3.0). The code that you provided (test.C) produced the expected signal plot (that is non-empty).
  • Using instead Centos 7.9 and the installation from CVMFS with the commands listed on the website (shown below) I got empty signal plots.
source /cvmfs/sft.cern.ch/lcg/views/LCG_101/x86_64-centos7-gcc11-opt/setup.sh
source /cvmfs/sft.cern.ch/lcg/views/LCG_101/x86_64-centos7-gcc11-opt/share/Garfield/setupGarfield.sh

As our production is on Centos, this is quite concerning. What could be the reason behind the difference?

Hmm, no idea. Sounds like NeBem bug if you got signals with drift tube example. I never used CVMFS, so can’t really test it.

Sorry for being late to the discussion!
Can you try picking up a nightly build from CVMFS instead of LCG_101?

. /cvmfs/sft-nightlies.cern.ch/lcg/views/dev3/latest/x86_64-centos7-gcc11-opt/setup.sh
. /cvmfs/sft-nightlies.cern.ch/lcg/views/dev3/latest/x86_64-centos7-gcc11-opt/share/Garfield/setupGarfield.sh

Hello @hschindl.
Unfortunately picking up the nightly build does not seem to work on the remote CentOS that I was using: I still get the empty signal plot.

Just to let you know that I can reproduce the issue on Centos7 but I have not yet identified what causes it. Will do more debugging…

1 Like

sorry it took so long… The issue should be fixed now:

Let me know if it persists…

Thank you so much @hschindl!
The plots are indeed working on our CentOS machine.
I’ll mark your last answer as the solution.

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