Dear experts
I wrote a program to simulate the drift of muon in a gas chamber, but the NewTrack never finish… here are my codes:
#include "xt_simulation.hxx"
simXT::simXT(int HV, int nevents, int fileId, std::string primary, double momentum):
fGas(NULL), fGeo(NULL), fCompE_field(NULL), fSensor_e(NULL), fMyCanvas(NULL), fPrimaryParticle(primary), fPrimaryMomentum(momentum)
{
fFileId = fileId;
fHV = HV;
fNEvent = nevents;
fSize_box = 147.42; //cm
fMyCanvas = new TCanvas();
std::cout<<"Canvas created"<<std::endl;
fFile_out = new TFile(Form("root/output_%d.root",fileId), "RECREATE");
}
simXT::~simXT(){
}
void simXT::SetupGas(std::string gasFile, std::string ionMobilityFile){
fGas = new MediumMagboltz();
const double pressure = 750.;
const double temperature = 293.15;
fGas->LoadGasFile(gasFile.data());
fGas->LoadIonMobility(ionMobilityFile.data());
fGas->Initialise(true);
// fGas->PrintGas();
}
void simXT::SetupGeo(std::string geoFile){
fGeo= new GeometrySimple();
TFile* f_geo = new TFile(geoFile.c_str(), "read");
if(!f_geo) std::cerr<<"Can not find Geo file!"<<std::endl;
TTree* t_geo = (TTree*)f_geo->Get("t");
t_geo->SetBranchAddress("x", &wire_x);
t_geo->SetBranchAddress("y", &wire_y);
t_geo->SetBranchAddress("type", &wire_type);
SolidBox* box = new SolidBox(0., 0., 0., 10., 10., fSize_box);
fGeo->AddSolid(box,fGas);
fGeo->PrintSolids();
// fGeo->SetMedium(fGas); // Why need to set medium twice?
std::cout<<"\033[32m";
std::cout<<fGeo->IsInside(0, 0, 0.)<<std::endl;
std::cout<<"\033[0m";
fCompE_field = new ComponentAnalyticField();
fCompE_field->SetGeometry(fGeo);
for(int i=0; i<t_geo->GetEntries(); i++){
t_geo->GetEntry(i);
if(wire_type == 0) fCompE_field->AddWire(wire_x, wire_y, r_field, fHV*wire_type, "field");
if(wire_type == 1) fCompE_field->AddWire(wire_x, wire_y, r_sense, fHV*wire_type, "sense");
// std::cout<<"added wire at "<<wire_x<<" "<<wire_y<<std::endl;
}
fCompE_field->AddReadout("sense");
fCompE_field->SetMagneticField(0, 0, 0);
std::cout<<"Set geometry for EM field"<<std::endl;
//view geometry
ViewCell* cellView = new ViewCell();
cellView->SetCanvas(fMyCanvas);
cellView->EnableWireMarkers();
cellView->SetComponent(fCompE_field);
cellView->SetArea(-5, -5, -fSize_box/2, 5, 5, fSize_box/2);
cellView->Plot2d();
fFile_out->cd();
fMyCanvas->Write();
fFile_out->Close();
}
void simXT::SetupSensor(){
const double tMin = 0.;
const double tMax = 2.;
const double tStep = 0.01;
const int nSteps = (int)(tMax-tMin)/tStep;
fSensor_e = new Sensor();
fSensor_e->AddComponent(fCompE_field);
fSensor_e->SetArea(-5., -5., -fSize_box/2, 5., 5., fSize_box/2);
fSensor_e->AddElectrode(fCompE_field, "sense");
fSensor_e->SetTimeWindow(0., tStep, nSteps);
}
void simXT::MakeTrack(){
TrackHeed* track = new TrackHeed();
track->SetSensor(fSensor_e);
track->SetParticle(fPrimaryParticle);
track->SetMomentum(fPrimaryMomentum);
track->EnableDebugging();
DriftLineRKF* drift = new DriftLineRKF();
drift->SetSensor(fSensor_e);
track->NewTrack(0.5, 0.5, 0., 0., -1., -1., 0.);
//cluster coordinates
double xc = 0., yc = 0., zc = 0., tc = 0.;
//number of electrons produced in a collision
int nc = 0;
//energy loss in a collision
double ec = 0.;
//dummy variable
double extra = 0.;
std::cout<<"Getting clusters..."<<std::endl;
while(track->GetCluster(xc, yc, zc, tc, nc, ec, extra)){ // Why not use references but normal variables?
std::cout<<nc<<" of electrons generated"<<std::endl;
double del_x, del_y, del_z, del_t;
double ke;
double dummy;
track->GetElectron(0, del_x, del_y, del_z, del_t, ke, dummy, dummy, dummy);
drift->DriftElectron(del_x, del_y, del_z, del_t);
std::cout<<"electron at "<<del_x<<" "<<del_y<<" "<<del_z<<std::endl;
}
ViewDrift* draw_drift = new ViewDrift();
track->EnablePlotting(draw_drift);
drift->EnablePlotting(draw_drift);
draw_drift->Plot();
}
void simXT::DrawField(){
ViewField* draw_field = new ViewField();
draw_field->SetArea(-5, -5, -fSize_box/2, 5, 5, fSize_box);
draw_field->SetComponent(fCompE_field);
draw_field->SetSensor(fSensor_e);
draw_field->Plot("v", "COLZ");
}
int main(int argc, char** argv){
std::cout<<argc<<" arguments"<<std::endl;
if(argc<9){
std::cout<<"/xt_simulation <gasfile><ionMobfile><geoFile><HV><number_events><fileId><particle><momentum[eV]>"<<std::endl;
return 1;
}
std::string gasFile = argv[1];
std::string ionMobFile = argv[2];
std::string geoFile = argv[3];
int hv = atoi(argv[4]);
int nevents = atoi(argv[5]);
int fileId = atoi(argv[6]);
std::string primary_particle = argv[7];
double primary_momentum = atoi(argv[8]);
TApplication app("app", &argc, argv);
simXT *usercode = new simXT(hv, nevents, fileId, primary_particle, primary_momentum);
usercode->SetupGas(gasFile, ionMobFile);
usercode->SetupGeo(geoFile);
//usercode->DrawField();
usercode->SetupSensor();
usercode->MakeTrack();
app.Run();
return 0;
}
And I got no error message but the NewTrack just never stop as follows
MediumMagboltz::Mixer:
Lowest ionisation threshold in the mixture: 10.67 eV (iC4H10)
MediumMagboltz::Mixer:
Energy [eV] Collision Rate [ns-1]
2.50 1812.84
7.50 3044.33
12.50 3944.59
17.50 4276.87
22.50 4492.48
27.50 4591.58
32.50 4683.67
37.50 4787.23
GeometrySimple::PrintSolids:
1 solid
Index Type Medium
0 box He/iC4H10
1
ComponentAnalyticField::AddReadout:
Readout group sense comprises:
9 wires
Set geometry for EM field
Sensor::AddElectrode:
Added readout electrode "sense".
All signals are reset.
Sensor::SetTimeWindow: Resetting all signals.
TrackHeed::UpdateBoundingBox:
Bounding box dimensions:
x: 10 cm
y: 10 cm
z: 147.42 cm
TrackHeed::UpdateBoundingBox:
Center of bounding box:
x: 0 cm
y: 0 cm
z: 0 cm
TrackHeed::Initialise:
Database path: /Users/siyuan/Physics/garfield/install/share/Heed/database/
TrackHeed::Initialise:
Cluster density: nan cm-1
Stopping power (restricted): nan keV/cm
Stopping power (incl. tail): nan keV/cm
W value: 39.51 eV
Fano factor: 0.179
Min. ionization potential: 10.55 eV
TrackHeed::NewTrack:
Track starts at (0.5, 0.5, 0) at time 0
Direction: (-0.707107, -0.707107, 0)