“Hey Garfield++ users, I’m simulating an MWPC with neBEM2d but the charges (e⁻/ions) aren’t following the E-field lines as expected. Has anyone encountered similar issues?”
I’m trying to create an MWPC (anode wires and cathode strips) using neBEM and run DriftLineRKF, but I’m encountering some issues. Here’s my approach:
type or paste code here
- Use neBEM to create anode wires and cathode strips (neBEM.segment), and use componentAnalyticField to generate weighting fields with the strips as sensing electrodes.
- Use TrackHeed to create a proton that generates electron-ion pairs after passing through the MWPC, then use DriftLineRKF to drift the electrons.
- Plan to collect induced signals on the strips and use the centroid method to determine the proton’s incident position.
My questions are:
- From the results shown in viewcell, viewfield, and viewmedium, I’ve generated the correct electric field, weighting field on the strips, gas, and cell. But when I output the drift lines, they look strange - the electrons and positive ions aren’t moving along the electric field direction shown in viewfield. Why is this happening? What did I do wrong?
- When using DriftLineRKF, I need to disable avalanche, otherwise I get warnings. I’ve seen this issue mentioned in other posts. Can I solve this by increasing the electric field range in the gas file?
- Back to question 1 - I tried using componentAnalyticField to create the MWPC (Garfield++ has an example, but it uses a planar cathode), and in that case the electrons drift normally. So I think I made a mistake in setting up the neBEM class, but I’m not sure where. I’ll include my relevant code and output results.
#include <TCanvas.h>
#include <TApplication.h>
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/ComponentNeBem2d.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
using namespace Garfield;
int main(int argc, char * argv[]){
TApplication app("app", &argc, argv);
//Add gas
MediumMagboltz gas;
gas.LoadGasFile("ar_93_co2_7_3bar.gas");
gas.LoadIonMobility("IonMobility_Ar+_Ar.txt");
// parameters
const double wirediameter = 30.e-4;
const double wirePitch = 0.2;
const int numberOfwires = 5;
const double stripWidth = 0.15;
const double stripPitch = 0.025;
const double stripHeigh = 0.003;
const int numberOfstrips = 6;
const int ndiv = 20;
const double chamberWidth = (numberOfwires + 1) * wirePitch;
const double chamberHeigh = 0.6;
const double wireVoltage = 2600;
const double stripVoltage = 0.;
//Add field
ComponentNeBem2d nebem;
nebem.SetMedium(&gas);
nebem.SetNumberOfDivisions(ndiv);
// nebem.EnableDebugging();
//Add wires
for(int i = 0; i < numberOfwires; i++){
const double x = (numberOfwires/2) * wirePitch - i * wirePitch ;
nebem.AddWire(x, 0., wirediameter, wireVoltage);
}
//Add strips
for(int j= 0; j < numberOfstrips; j++){
const double x0 = chamberWidth/2 - j * (stripPitch + stripWidth);
const double x1 = chamberWidth/2 - j * (stripPitch + stripWidth) - stripWidth;
const double y00 = -chamberHeigh/2;
const double y01 = -chamberHeigh/2;
const double y10 = chamberHeigh/2;
const double y11 = chamberHeigh/2;
nebem.AddSegment(x0, y00, x1, y01, stripVoltage);
nebem.AddSegment(x0, y10, x1, y11, stripVoltage);
}
Sensor sensor;
sensor.AddComponent(&nebem);
sensor.SetTimeWindow(0.,50., 1000);
//Add weighing field
ComponentAnalyticField wfield;
for(int b = 0; b < numberOfwires; b++){
const double x = (numberOfwires/2) * wirePitch - b * wirePitch ;
wfield.AddWire(x, 0., wirediameter, wireVoltage);
}
wfield.AddPlaneY( chamberHeigh/2, stripVoltage);
wfield.AddPlaneY(-chamberHeigh/2, stripVoltage);
for(int k = 0; k < numberOfstrips; k++){
std::string top_strip_lable = "top_strip_" + std::to_string(k);
std::string bottom_strip_lable = "bottom_strip_" + std::to_string(k);
const double s0 = chamberWidth/2 - k * (stripPitch + stripWidth);
const double s1 = chamberWidth/2 - k * (stripPitch + stripWidth) - stripWidth;
wfield.AddStripOnPlaneY('z', chamberHeigh/2, s0, s1, top_strip_lable, 0.02);
wfield.AddStripOnPlaneY('z', -chamberHeigh/2, s0, s1, bottom_strip_lable, 0.02);
sensor.AddElectrode(&wfield, top_strip_lable);
sensor.AddElectrode(&wfield, bottom_strip_lable);
}
// View the field ...
ViewMedium mediumview(&gas);
ViewCell cellview(&nebem);
ViewField fieldview(&sensor);
ViewField wfieldview(&sensor);
wfieldview.SetComponent(&wfield);
// //Add particle
TrackHeed track;
track.SetSensor(&sensor);
const double KineticEnergyOfParticle = 170.e6;
track.SetParticle("p");
track.SetKineticEnergy(KineticEnergyOfParticle);
const int nParticles = 1;
// track.EnableElectricField();
// track.EnablePlotting(&driftview);
// constexpr bool DeltaElectron = true;
//Add drift
DriftLineRKF drift(&sensor);
//View the Drift
ViewDrift driftview;
drift.EnableAvalanche(false);
track.EnablePlotting(&driftview);
drift.EnablePlotting(&driftview);
const double x0 = 0.1;
const double y0 = chamberHeigh/2 - stripHeigh;
// const double y0 = 0.2;
const double z0 = 0.;
const double t0 = 0.;
for(int a = 0; a < nParticles; a++){
track.NewTrack(x0, y0, z0, t0, 0.,-1., 0.);
for(const auto& cluster : track.GetClusters()){
for(const auto& electron : cluster.electrons){
drift.DriftElectron(electron.x, electron.y, electron.z, electron.t);
}
// for(const auto& ion : cluster.ions){
// drift.DriftIon(ion.x, ion.y, ion.z, ion.t);
// }
}
}
TCanvas* c1 = new TCanvas("c1", "", 600, 600);
driftview.SetCanvas(c1);
driftview.SetArea(-chamberWidth/2,-chamberHeigh/2, chamberWidth/2, chamberHeigh/2);
driftview.Plot2d(true);
constexpr bool plotgas = true;
if(plotgas){
mediumview.PlotElectronTownsend();
}
constexpr bool plotField = true;
if(plotField){
fieldview.SetArea(-chamberWidth/2,-chamberHeigh/2, chamberWidth/2, chamberHeigh/2);
fieldview.PlotContour("v");
// fieldview.PlotProfile(0,-0.25,0.,0,0.25,0,"v",false);
}
constexpr bool plotcell = true;
if(plotcell){
cellview.SetCanvas(c1);
cellview.SetArea(-chamberWidth/2,-chamberHeigh/2, chamberWidth/2, chamberHeigh/2);
cellview.Plot2d();
}
constexpr bool plotwfield = true;
if(plotwfield){
// wfieldview.SetArea(-chamberWidth/2, -chamberHeigh/2, chamberWidth/2, -chamberHeigh/2 + 0.02);
// wfieldview.PlotContourWeightingField("bottom_strip_4", "v");
wfieldview.SetArea(-chamberWidth/2, chamberHeigh/2 - 0.02, chamberWidth/2, chamberHeigh/2);
wfieldview.PlotContourWeightingField("top_strip_4", "v");
}
// constexpr bool plotdrift = true;
// if(plotdrift){
// cellview.SetArea(-chamberWidth/2,-chamberHeigh/2, chamberWidth/2, chamberHeigh/2);
// driftview.Plot2d();
// }
if(plotwfield | plotField | plotcell ){
app.Run(true);
}
return 0;
}
