Hi, I am trying to create a simulation of an isobutane-filled MWPC with two wire planes and a mylar cathode in between them. In this simulation, a 5-MeV alpha particle enters the volume of gas and leaves an ionization trail, which I consider primary electrons. I then want to take this trail of primary ionization electrons and propagate it through an MWPC, creating an avalanche of secondary electrons.
I am presently trying to make a movie for visualizing the process. I use HEED example to calculate the first ionization trail, save the positions and energies of resulting cluster primary electrons to a vector, and then feed this vector as primaries to code borrowed from single GEM movie example, which then calculates avalanches of secondary electrons.
The resulting simulation is very slow, which is perhaps unsurprising given the large number of particles I start with (a few hundred, which multiply to a few hundred thousand). More importantly, I am not sure I am approaching it correctly, and not fully understanding the HEED step of the process. I attach an image below of frame zero from my section of the code, which uses GEM example to produce an ionization movie. In my understanding, frame zero should represent the initial conditions fed to the GEM example. In the image, there are several sideways branches of ionization electrons visible, which suggests that there are already some secondary ionizations from electrons in the first step of the movie not resulting from the visible central alpha trace. I think I am likely misunderstanding how HEED works.
Am I interpreting this correctly? And if so, what would you recommend? I attach my full source code as well.
#include <iostream>
#include <cmath>
#include <vector>
#include <TFile.h>
#include <TCanvas.h>
#include <TH1I.h>
#include <TH2F.h>
#include <TString.h>
#include <TLatex.h>
#include <TApplication.h>
#include <TSystem.h>
//neBEM
#include "Garfield/SolidWire.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/MediumMylar.hh"
#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/ViewField.hh"
//Heed
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/ViewDrift.hh"
//Gem movie
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
using namespace Garfield;
int main(int argc, char* argv[]){
TApplication app("app", &argc, argv);
const float P = 5.250432; //7 mbar in torr
const float V_wire = 650.;
const int N_alpha = 1;
const double E_alpha = 5.; //MeV alpha particle
const double radius = 0.003; //Wire radius
const double halflength = 5.; //Wire half-length
const float pitch = 0.2; //Distance between wires in cm
//const float mylar_thickness = 0.00009;
const float mylar_thickness = 0.;
const float plot_z_height = 0.9;
const float wires1_height = 0.3;
const float wires2_height = 0.9;
const float cathode_height = 0.6;
const float cathode_eloss = 0.106451; //MeV, SRIM calculated energy loss of 5 MeV alpha in 0.9 um of Mylar
const bool gif = 1; //Boolean to animate avalanches
TFile* rootout = TFile::Open("MWPC.root", "recreate"); //ROOT output
//Output histograms
TH1I* h_de = new TH1I("h_de", "h_de", 700, 0, 7);
TH1I* h_ne = new TH1I("h_ne", "h_ne", 1e4, 0, 1e4);
TH1F* h_a_range = new TH1F("h_a_range", "h_a_range", 1000, 0, 1);
TH1D* h_ne_fin = new TH1D("h_ne_fin", "h_ne_fin", 1e6, 0, 1e8);
TH2F* h_E = new TH2F("h_E", "h_E", 1000, -5, 5, 1000, -5, 5);
//Set up gas
MediumMagboltz gas("ic4h10", 100.);
gas.SetPressure(P);
gas.SetTemperature(293.15); //Kelvin
gas.EnableThermalMotion(true);
//Geometry
GeometrySimple geo;
std::vector< std::vector<SolidWire> > Wires; //Vector to list wires
Wires.resize(2); //Column 0 is for x wires, column 1 is for y wires
float start_wire_pos = -1. * halflength; //Position in cm of first wire in cm in along +/-5 cm axis
float end_wire_pos = halflength; //Same as above but last wire
int N_wires = floor((end_wire_pos - start_wire_pos)) / pitch; //Number of wires fitting into that space
MediumConductor metal;
std::cout << std::endl << "Processing " << N_wires << " wires along x and y axis each." << std::endl << std::endl;
//Iteratively calculate x and y wire plane positions and add them to vector
for(int i = 0; i <= N_wires; i++){
float wire_pos = start_wire_pos + (float)i * pitch;
if(i == 1){std::cout << " Starting wire position " << wire_pos << std::endl;}
SolidWire wire1(0., wire_pos, wires1_height, radius, halflength, 1, 0, 0);
SolidWire wire2(wire_pos, 0., wires2_height, radius, halflength, 0, 1, 0);
wire1.SetBoundaryPotential(V_wire);
wire2.SetBoundaryPotential(V_wire);
Wires[0].push_back(wire1);
Wires[1].push_back(wire2);
}
//Add wires to geometry
for(int i = 0; i < (int)Wires[0].size(); i++){
geo.AddSolid(&Wires[0][i], &metal);
geo.AddSolid(&Wires[1][i], &metal);
}
MediumConductor mylar;
//MediumMylar mylar;
//Create a solid cathode plane made of Mylar between x and y wire planes
SolidBox cathode(0., 0., cathode_height, halflength, halflength, mylar_thickness / 2.);
//Cathode is grounded
cathode.SetBoundaryPotential(0.);
geo.SetMedium(&gas);
geo.AddSolid(&cathode, &mylar);
//Initialize neBEM calculation for electric fields
ComponentNeBem3d nebem;
nebem.SetGeometry(&geo);
nebem.SetTargetElementSize(1.);
nebem.UseLUInversion();
nebem.SetNumberOfThreads(18);
nebem.Initialise();
std::cout << std::endl;
//Calculate field and potential at point
double ex = 0., ey = 0., ez = 0., V = 0.;
int status = 0;
Medium* medium = nullptr;
//Make a 2D plot of electric field
for(int i = 0; i < h_E->GetNbinsX(); i++){
for(int j = 0; j < h_E->GetNbinsY(); j++){
nebem.ElectricField(-5. + 10. / (float)i, -5. + 10. / (float)j, plot_z_height, ex, ey, ez, V, medium, status);
h_E->SetBinContent(i + 1, j + 1, V);
h_E->SetBinError(i + 1, j + 1, 0.);
}
}
std::cout << std::endl << " Plotting electric fields " << std::endl << std::endl;
ViewField fieldView1(&nebem);
ViewField fieldView2(&nebem);
ViewField fieldView3(&nebem);
fieldView1.SetArea(-5., -5., 0., 5., 5., 1.);
fieldView2.SetArea(-5., -5., 0., 5., 5., 1.);
fieldView3.SetArea(-5., -5., 0., 5., 5., 1.);
fieldView1.SetVoltageRange(-2000., 2000.);
fieldView2.SetVoltageRange(-2000., 2000.);
fieldView3.SetVoltageRange(-2000., 2000.);
fieldView1.SetPlane(0., 0., 1., 0., 0., cathode_height);
fieldView1.PlotContour("e");
fieldView2.SetPlaneXZ();
fieldView2.PlotContour("e");
fieldView3.SetPlaneYZ();
fieldView3.PlotContour("e");
ViewGeometry geomView3d(&geo);
geomView3d.Plot();
//Gas file handling
gas.WriteGasFile("isobutane.gas");
//Track particle
Sensor sensor;
sensor.AddComponent(&nebem);
rootout->cd();
TCanvas canvas("c", "c", 600, 600);
canvas.SetLeftMargin(0.16);
double prev_z = 0.;
std::cout << std::endl << " Simulating avalanches for " << N_alpha << " alpha(s)" << std::endl << std::endl;
//Loop over desired number of alpha particles
for(int i = 0; i < N_alpha; i++){
TrackHeed heed; //HEED calculates ionizations along alpha particle track in gas
heed.SetSensor(&sensor);
heed.EnableElectricField();
heed.DisableDeltaElectronTransport();
heed.Initialise(&gas, 1);
//const float p_alpha = 1.e6 * sqrt(2. * 4. * 931.494 * E_alpha); //4 amu alpha momentum, momentum in eV/c
heed.SetParticle("ALPHA");
//heed.SetMomentum(p_alpha);
heed.SetKineticEnergy(1.e6 * E_alpha); //This is in units of eV
std::cout << std::endl << " Simulating " << E_alpha << "-MeV alpha with momentum " << heed.GetMomentum() / 1.e6 << " MeV/c (beta = " << heed.GetBeta() << "% of c, read-in energy " << heed.GetKineticEnergy() << " eV)" << std::endl << std::endl;
heed.NewTrack(0.1, 0.1, 0., 0., 0, 0, 1); //Particle goes in z direction
double esum = 0.; //Total energy loss
int nsum = 0; //Total number of electrons produced along the track
AvalancheMicroscopic aval(&sensor); //This is a list of initial particles for simulating avalanche cascades. For this example, we add all the electrons produced by alpha ionization here and propagate them later through the detector.
int n_clusters = 0;
rootout->cd();
TH1D* h_eloss = new TH1D(Form("h_eloss%d", i), Form("h_eloss%d", i), 1000, 0, 1);
for(const auto& cluster : heed.GetClusters()){
n_clusters++;
if(n_clusters == heed.GetClusters().size()){
h_a_range->Fill(cluster.z);
}
//Adjust alpha energy for energy lost to ionization
heed.SetKineticEnergy(heed.GetKineticEnergy() - cluster.energy);
//When alpha crosses the cathode, subtract SRIM-based energy loss in eV
if(prev_z < cathode_height && cluster.z >= cathode_height){
heed.SetKineticEnergy(heed.GetKineticEnergy() - 1.e6 * cathode_eloss);
}
prev_z = cluster.z;
h_eloss->SetBinContent(h_eloss->FindBin(cluster.z), heed.GetKineticEnergy() / 1.e3);
esum += cluster.energy;
nsum += cluster.electrons.size();
//Save all the ionization electrons for later cascade propagation
for(const auto& electron : cluster.electrons){
aval.AddElectron(electron.x, electron.y, electron.z, electron.t, electron.e);
} //End of loop over electrons in a cluster
} //End of loop over collision clusters
std::cout << std::endl << " Total energy loss " << esum << " eV" << std::endl;
std::cout << " Total secondary electrons " << nsum << std::endl << std::endl << std::endl;
h_de->Fill(esum / 1.e3);
h_ne->Fill(nsum);
std::cout << "Animating the cascade" << std::endl;
AvalancheMC drift(&sensor);
ViewDrift driftView;
aval.EnableExcitationMarkers(false);
aval.EnableIonisationMarkers(true);
aval.EnableAttachmentMarkers(false);
aval.EnablePlotting(&driftView);
drift.EnablePlotting(&driftView);
driftView.SetPlaneXZ();
driftView.SetArea(-5., 0., 5., 1.);
//driftView.SetArea(-1., 0., 1., 0.4);
driftView.SetCanvas(&canvas);
const int nFrames = 50; //Number of steps for gif of avalanches calculation
double tmin = 0.;
double dt = 0.3; //Time steps in avalanche in units of ns
drift.SetTimeSteps(dt);
rootout->cd();
TH1D* mult = new TH1D(Form("mult%d", i), Form("mult%d", i), nFrames, tmin, tmin + nFrames * dt);
for(int frame = 0; frame < nFrames; frame++){
std::cout << " Frame " << frame << std::endl;
std::cout << " Ne-: " << aval.GetElectrons().size() << std::endl;
//Make a histogram of total electrons inside the detector per unit time
mult->SetBinContent(mult->FindBin(tmin), aval.GetElectrons().size());
driftView.Clear();
if(!aval.GetElectrons().empty()){
aval.SetTimeWindow(tmin, tmin + dt);
drift.SetTimeWindow(tmin, tmin + dt);
aval.ResumeAvalanche();
drift.ResumeAvalanche();
canvas.cd();
canvas.Update();
gSystem->ProcessEvents();
tmin += dt;
driftView.Plot2d(true, true);
if(gif){
if(frame == nFrames - 1){
canvas.Print("MWPC_Movie.gif++");
}
else{
canvas.Print("MWPC_Movie.gif+3");
}
}
} //End of if avalanche is empty
else{frame = nFrames;}
} //End of loop over frames for a movie
h_ne_fin->Fill(aval.GetElectrons().size());
}
rootout->Write();
rootout->Close();
std::cout << std::endl << "Finished" << std::endl;
app.Run();
return 0;
}
