Hello @mkelsey
I tried to follow the first way…i.e.
In the steppingAction I added:
const G4Event* evt = G4RunManager::GetRunManager()->GetCurrentEvent();
if(NextVol && ThisVol->GetName()=="PlasticScintillator" && NextVol->GetName()=="PlasticScintillator") {
int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
analysisManager->FillNtupleDColumn(9,0, pdg);
int trackID = step->GetTrack()->GetTrackID();
analysisManager->FillNtupleDColumn(9,1, trackID);
int evtID = evt->GetEventID();
analysisManager->FillNtupleDColumn(9,2, evtID);
double partEdep = step->GetTotalEnergyDeposit();
analysisManager->FillNtupleDColumn(9,3, partEdep);
analysisManager->AddNtupleRow(9);
}
I used the filter
if(NextVol && ThisVol->GetName()=="PlasticScintillator" && NextVol->GetName()=="PlasticScintillator") {
B1SteppingAction.cc (43.3 KB)
because I need only particles depositing energy in my scintillator. Is the code right or did I make mistakes??
Now I’ve the plots
This is the trackID
This is the eventID
and this is the deposited energy
I explain you my goal…
In november 2021 we made some measurement of radiation hitting a scintillator connected to an ArduSiPM. The ArduSiPM measured a counting rate of 4.410^5 particles/second … I simulated the experimental setup by Geant4, but I got a counting rate of 1.510^7 particles/second… obiouvsly the measurement by ArduSiPM is limited by many factors… One of these factors is that ArduSiPM can detect only particles depositing energy <1MeV in the scintillator…therefore, from the total particles entering in the scintillator that I save by the code
if(NextVol && ThisVol->GetName()=="World" && NextVol->GetName()=="PlasticScintillator") {
int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
analysisManager->FillNtupleDColumn(6,0, pdg);
double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
analysisManager->FillNtupleDColumn(6,1, kinEnergy);
double MomDirx = step->GetTrack()->GetMomentumDirection().x();
analysisManager->FillNtupleDColumn(6,2, MomDirx);
double MomDiry = step->GetTrack()->GetMomentumDirection().y();
analysisManager->FillNtupleDColumn(6,3, MomDiry);
double MomDirz = step->GetTrack()->GetMomentumDirection().z();
analysisManager->FillNtupleDColumn(6,4, MomDirz);
double Vertx = step->GetTrack()->GetVertexPosition().x();
analysisManager->FillNtupleDColumn(6,5, Vertx);
double Verty = step->GetTrack()->GetVertexPosition().y();
analysisManager->FillNtupleDColumn(6,6, Verty);
double Vertz = step->GetTrack()->GetVertexPosition().z();
analysisManager->FillNtupleDColumn(6,7, Vertz);
double Mom = step->GetTrack()->GetDynamicParticle()->GetTotalMomentum();
analysisManager->FillNtupleDColumn(6,8, Mom);
double Dirxsp = step->GetTrack()->GetPosition().x();
analysisManager->FillNtupleDColumn(6,9, Dirxsp);
double Dirysp = step->GetTrack()->GetPosition().y();
analysisManager->FillNtupleDColumn(6,10, Dirysp);
double Dirzsp = step->GetTrack()->GetPosition().z();
analysisManager->FillNtupleDColumn(6,11, Dirzsp);
double Momx = step->GetTrack()->GetDynamicParticle()->GetMomentum().x();
analysisManager->FillNtupleDColumn(6,12, Momx);
double Momy = step->GetTrack()->GetDynamicParticle()->GetMomentum().y();
analysisManager->FillNtupleDColumn(6,13, Momy);
double Momz = step->GetTrack()->GetDynamicParticle()->GetMomentum().z();
analysisManager->FillNtupleDColumn(6,14, Momz);
double xprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().x())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
analysisManager->FillNtupleDColumn(6,15, xprime);
double yprime=atan((step->GetTrack()->GetDynamicParticle()->GetMomentum().y())/(step->GetTrack()->GetDynamicParticle()->GetMomentum().z()));
analysisManager->FillNtupleDColumn(6,16, yprime);
analysisManager->AddNtupleRow(6);
}
this is the energy spectrum of particles entering in the scintillator (here I’ve only 38 particles because I runned only 3*10^6 primary electrons)
Now from the 38 particles entering in the scintillator, I’ve to select only particles depositing an energy >1MeV in the scintillator…I created the TTree with trackID, eventID and deposited energy (hoping that I did it in the right way)…now how can understand how many particles they are???