DriftLineRKF looping indefinitely in MPGD implementation

[gzunica]

Dear Experts,

I am trying to simulate drift lines with DriftLineRKF in a proof of concept MPGD-like implementation.

The geometry consists of two plates (cathode on top, anode on bottom) with a thin dielectric medium on top of the bottom plate. One readout wire and two “field” anodic wires are resting on top of the glass. The overall gap thickness is 5mm.

I have tested multiple combinations of “SetMaximumStepSize” and "SetIntegrationAccuracy” to no avail. More specifically, with some voltage combinations the computation converges, with some others it gets stuck an a loop trying to optimise the step size.

I attach a minimal example.

Is there something that i am missing?

Thank you in advance!

using namespace Garfield;

int main(int argc, char* argv[]) {


  MediumMagboltz gas;//("Ar", 87.5, "CO2", 12.5);
  gas.LoadGasFile("/afs/cern.ch/work/g/gzunica/garfield/Examples/GasFile/ar_93_co2_7.gas");
  gas.LoadIonMobility("/afs/cern.ch/work/g/gzunica/garfield/Data/IonMobility_Ar+_Ar.txt");
  gas.SetPressure(760.);
  gas.SetTemperature(293.15);


  MediumConductor Cu;
  MediumPlastic glass;
  glass.SetDielectricConstant(10.0);
  MediumPlastic kapton;
  kapton.SetDielectricConstant(3.0);

  // Geometry.
  GeometrySimple geo;

  const double cat_v=-2000.;
  const double sense_v=+10.;
  const double field_v=-1700.;
  const double backplane_v=-1500.;
  

  SolidBox box1(0, 0, 0.5, 3.2,2.5, 0.0);
  box1.SetBoundaryPotential(cat_v);
  SolidBox box2(0, 0, -0.105, 3.2,2.5, 0.1);
  box2.SetBoundaryDielectric();
  SolidBox box3(0, 0, -0.205, 3.2,2.5, 0.0);
  box3.SetBoundaryPotential(backplane_v);
 
   geo.AddSolid(&box1, &Cu);
   geo.AddSolid(&box2, &glass);
   geo.AddSolid(&box3, &Cu);
   

  const double radiusS = 0.0015; //0.0015
  const double radiusF = 0.005;
  const double halflength = 2.5;
  double passo = 0.4;


  SolidWire wireS(xposS, 0.0, -0.0035, radiusS, halflength, 0, 1, 0); 
  wireS.SetBoundaryPotential(sense_v);
  wireS.SetLabel("readout");
  geo.AddSolid(&wireS, &Cu);


  SolidWire wireF1(xposF, 0.0, 0.0, radiusF, halflength, 0, 1, 0);
  wireF1.SetBoundaryPotential(field_v);
  SolidWire wireF2(xposF+passo, 0.0, 0.0, radiusF, halflength, 0, 1, 0); 
  wireF2.SetBoundaryPotential(field_v);
  geo.AddSolid(&wireF1, &Cu);
  geo.AddSolid(&wireF2, &Cu);


  geo.PrintSolids();
  geo.SetMedium(&gas);

  ComponentNeBem3d nebem;
  nebem.SetGeometry(&geo);
  nebem.SetTargetElementSize(0.2);
  nebem.SetMinMaxNumberOfElements(1, 100); 

  nebem.Initialise();

  Medium* medium = nullptr;
  double ex = 0., ey = 0., ez = 0., v = 0.;
  int status = 0;
  nebem.ElectricField(0, 0, 0, ex, ey, ez, v, medium, status);
  std::printf("E = (%15.8f, %15.8f %15.8f), V = %15.8f, status = %d\n", ex, ey,
              ez, v, status);


 Sensor sensor(&nebem);
 sensor.AddElectrode(&nebem, "readout");


 Shaper shaper(3,0.5,1,"unipolar");
 sensor.SetTransferFunction(shaper);

 sensor.SetTimeWindow(0., 0.5, 200);


  DriftLineRKF driftRKF(&sensor);
  driftRKF.SetIntegrationAccuracy(1e-6);
  driftRKF.SetMaximumStepSize(0.0005);
  driftRKF.EnableIonTail();
  driftRKF.EnableAvalanche();
  //driftRKF.EnableDebugging();


  ViewDrift driftView;
  driftView.SetPlaneXZ();
  driftView.SetArea(-5.,-5.,-0.5,5.,5.,0.6);
  heed.EnablePlotting(&driftView);
  driftRKF.EnablePlotting(&driftView);




  ViewSignal signalView(&sensor);


   std::vector<double> xe;
   std::vector<double> ye;
   std::vector<double> ze;
   xe.push_back(0.);

   ye.push_back(0.);

   ze.push_back(0.1);

   for(int i=0; i<xe.size();i++){

       driftRKF.DriftElectron(xe.at(i), ye.at(i), ze.at(i), 0);
       double ne=0;
       double ni=0;

         driftRKF.GetAvalancheSize(ne,ni);
         std::cout<< "ne " << ne << ", ni "<< ni <<"\n";
     }
}

[bellenot]

Welcome to the ROOT Forum!
I’m sure @hschindl can help

[hschindl]

Hi,
what are the values of xposS and xposF?

[gzunica]

Hello,

those are the positions on the X axis of the wires. The “readout” wire should be at x=0 while the field wires should be at x=±0.2 cm.
I just realised i have deleted their definition while posting the script, apologies!

[hschindl]

Hi,
sorry for the delay! The issue is that the ion drift lines jump across the box electrodes. I should think of a better way to prevent/fix this, but as a work-around, you can/should set explicitly the active region of the sensor, something like this:

sensor.SetArea(-3., -3., -0.205, 3., 3., 0.5);

Alternatively, you could also add a SolidBox for the gas volume and remove the line geo.SetMedium(&gas);.

Not directly related to your question, but:

• I assume the glass layer in your detector is weakly conductive? If it is, that means for the purposes of calculating the electrostatic field in your detector, it is a conductor, while for the purposes of calculating the prompt component of the weighting potential, it is an insulator.
• The value of the max. step size limit (5 μm) seems a bit small, I would try increasing it a bit to speed up the calculation.

[gzunica]

Indeed using sensor.SetArea() did the solve the issue!

As for the glass layer, i’m not sure i understand.
I’m using the standard “MediumPlastic” definition, explicitly setting the dielectric constant.
I had assumed it’s non driftable and non ionisable (as for the MediumPlastic.h), although it would be useful if this “glass” could be tweaked to have some weak conductivity (for ions at least). Or better still if it could mimick more of a “resistive” medium.

[hschindl]

Hi,
regarding the glass layer, what I meant is that in neBEM a solid can only be either a conductor (at fixed potential) or a dielectric. If you describe the glass layer as a dielectric (i. e. an insulator), the electric field will be “wrong” because even though the resistivity is high, from an electrostatic point of view it is nevertheless a conductor. For calculating the prompt weighting field, treating the glass layer as a dielectric is correct (and the time constant is typically so high that you can neglect the contribution from the delayed weighting field).

[gzunica]

Hello,

indeed my aim is to have a feeling of how the electric field behaves when “crossing” the glass. From what i understood neBEM should be handling the effect of different mediums on the electric field, is this correct?

Just to be clear, the glass was not intended to “block” the electric field, but rather to stop particles drifting to the backplane electrode.