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";
}
}