Hello, in my stepping action I get the deposited energy integrated during the run…should be possible to get the deposited energy by a single particle?
I mean…I run a primary electron and it produces some secondary particles…I want to know how much deposited energy EACH ONE of the secondary particles deposites…how should I do?
1 Like
Several ways. You can create an N-tuple with three branches: eventID, trackID and Edeposit. In your stepping action or SD, accumulate the Edeposit for each hit into the correct N-tuple row. If you don’t want to reach in and modify the N-tuple content, you could do it with something like a std::map<G4int, G4double>, keying on the track ID. Keep the eventID as well, and when that changes, clear the map and start over.
Your best bet is to probably implement this stuff in your EventAction, with a public member function you can call from your stepping action or SD. In the EventAction, you can have begin-event and end-event functions that take care of saving the map to an N-tuple, then clearing it.
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???
would it make sense to filter for <1MeV photons/gammas hitting your detector instead of energy deposit per particle?
EDIT there was a misunderstanding…the right cut is Edep>1MeV…
So @mkelsey or @weller , please, do you know how to get the deposited energy by each secondary partlicle in the scintillator to make this cut Edep>1MeV?
As I wrote here Deposited energy by a single particle - #3 by faca87 I filled the N-tuple in the stepping action with trackId, Edep, and eventID in the scintillato volume,…but now…how can see the total deposited energy by each secondary particle???
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);
}
Thank you
maybe you can get inspiration from example B1? I think they do exactly what you need there 
The energy deposited is collected step by step for a selected volume
in B1::SteppingAction and accumulated event by event in B1::EventAction.At end of event, the value acummulated in B1::EventAction is added in B1::RunAction
and summed over the whole run (see B1::EventAction::EndOfevent()).
Hello @weller I use the B1 code …for example the deposited energy in my scintillator is (in the Stepping action)
if (volume == fScoringVolumeScint) {
G4double edepStepScint = step->GetTotalEnergyDeposit();
fEventAction->AddEdepScint(edepStepScint);
}
but this code measures the deposited energy during the event …it means that if the primary particle A produces two secondary particles B,C in the target and the secondary particles B,C enter in the scintillator the deposited energy by B and C is summed and GEANT4 says that it is the deposited energy by A
on the contrary I need to know the deposited energy by B and the deposited energy by C.…because I need to know the deposited energy by each particle entering in the scintillator!
1 Like
hm. can you make use of the trackID to keep track of the individual secondaries, and add up the deposited energy e.g. in a vector?
@weller do you have an example to do it? I wrote the code
> if(volume == fScoringVolumeScint) {
> 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);
> }
to get the track and the deposited energy…but I don’t know how to sum up the single secondary partyicle entering in the scintillator
If you want the total energy deposited by everything within your scintillator, for a single particle entering from the outside, that should be extremely easy.
If you want the total of everything, then you don’t need a bunch of conditionals. If you only want hits in a specific volume, that’s exactly what an SD is good at.
If you only want one particle entering your detector, then just throw one particle per event.
This very simple, very basic setup is demonstrated in most of the examples, starting with B1.
Hello @mkelsey
I simulated 6.25*10^8 primary electrons hitting a target in a vacuum chamber. In front of the vacuum chamber there is a scintillator, then some of the secondary particles produced in the target and outgoing the vacuum chamber can hit the scintillator. By this kinetic energy spectrum you see that
24431 e- (called h1 in the graph), 4155 e+ (called h2) and 121174 gamma (called h3) enter in the scintillator…then I’ve 24431+4155+121174=149760 particles entering in the scintillator.
I’m able to make the spectrum of the deposited energy in the scintillator , i.e. this one (based on the B1 example)
but as you can read in the statbox the entries of this spectrum are 6.25*10^8…because this spectrum refers to the primary electrons… i.e. if a primary particle A produces some daughters particles, the sum of the deposited energy by the daughters is shown as the deposited energy by A!
On the contrary, I need the 3 deposited energy spectra:
a. the first one refering to the 24431 electrons entering in the scintillator (i.e. the entries of the graph must be the 24431 electrons entering in the scintillator)
b. the second one refering to the 4155 positrons entering in the scintillator (i.e. the entries of the graph must be the 4155 positrons entering in the scintillator)
c. the last one refering to 121174 photons entering in the scintillator (i.e. the entries of the graph must be the 121174 photons entering in the scintillator)
well, the suggestions are on the table, I fear it is up to you to implement them…
additional idea that is probably easier to do in geant4 but more difficult in the post processing:
store the parent ID and the trackID, and the particle type of each energy deposit event with deposited energy > 0 in a ntuple (plus of course the deposited energy), and then do post processing by summing all ntuples that have the same parentID, trackID and particle type, and then accumulate the sum for the respective histogram.
can you help me with the code please? I really don’t know how to implement it…
I’m stuck here
if(volume == fScoringVolumeScint) {
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);
}
because I don’t know how tu sum over the secondaries entering in the scintillator instead over the primary particle. Can you show me an input of the code please?
like so? of course when defining the ntuple, account for the extra integer column 
(trackid, pdg and eventID could also be defined as integer columns)
the other approaches I would have to learn myself, probably looking through all the examples and check if something alike was already implemented, then copy pasting code blocks that seem useful, trying things, failing, trying again, failing again, until it works…
@weller I filled that Ntupla as @mkelsey said me to do in his first message, but now I don’t know how to continue…because now I’ve just the deposited energy for each step, not for each secondary particle entering in the scintillator!
I see that you suggested to add the parentID culomn, but my problem is that I don’t know how to relate the deposited energy to each secondary particle…
Unfortunately I searched both in the forum and in the Google, and I didn’t find discussions about the depositerò energy by eaxh secondary particle…I guess that isn’t a common study the depositerò energy by each secondary particle, then there aren’t examples to do it…
the parentID will give you the ID of the parent particle in that event, e.g., if it is a secondary directly produced from the primary particle, the parentID will be 0. All secondary particles with the same parent ID, same trackID and the same event ID are in fact the same particle. So in your n-tuple, you can now look for all rows with unique set of parentID, trackID, eventID and particle type, and then summarize all the energy deposit entries with that unique combination.
just realized that maybe the parentID is redundant, because the trackID already is uniquely identifying the secondary. so just read the above statement, and strip away all the parentIDs 
From the history of your posts, you probably use root to create the histograms!? So now it is more a question of root. Sorry I cannot be of more help here, I have never used root.
Hello @weller thank you
but when the particles enter in the scintillator, they produce other particles…then by the code
if(volume == fScoringVolumeScint) {
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);
// add this:
analysisManager->FillNtupleIColumn(9,4, step->GetTrack()->GetParentID());
analysisManager->AddNtupleRow(9);
}
I get the parent ID, the trackID and the eventID of particles produced in the scintillator, not the parentID , trackID and eventID of the particles entering in the scintillator…am I wrong?
Example the primary electron A interacts with the target producing 3 secondary particles B,C,D. Particles B,C outgo the vacuum chamber and enter in the scintillator producing other particles E,F…
then this code, will give me the parentID,the trackID and the eventID of E,F …in addition, moving in the scintillatore E,F can produce other particles (G,H,I) …then event if by ROOT I can sum the energies for a particular parentID, trackID and eventID I will get the energies of E,F,G,H,I not the energies of B,C…
am I worng?
Yes, I produce Root TTree by Geant4…I will think about ROOT later…now my problem is the previous one (i.e. interacting with the scintillator particles entering in it produce other particles, then I get the trackID,parentID and eventID of particles produced in the scintillator instead of particles entering in it)
how often do you have more than 1 particle of the same type entering the scintillator in 1 single event (with 1 event being 1 primary particle!)? and can you distinguish between “energy deposited by 2 electrons entering the scintillator simultaneously” and “energy deposited by 1 single electron”?
I don’t thinck it’s easy to see…because I’ve the plots of the steps and they are a lot (in this case I just simulate 110^6 primary particles to try…but I’ve to simulate 6.2510^8)
TrackID
EventID
I guess that to do it’s needed to relate the deposited energy to the particle who deposited it…but I don’t know how to do it…
i think this plot is good news for you, the way I understand it. of your 178 entries in that simulation, 174 are in the trackID == 2 bin, and only 4 are higher order - that is >97% of all events have 1 or 0 particles ending up in the scintillator. I would also not be surprised, if the very few events that give you trackID 5 or the one event with trackID 4 does not simultaneously have the trackID 2 particle in the scintillator. and even more: you want to distinguish between particle types anyways, so the changes that both daugthers are of the same kind is even more unlikely.
trackID 2 means it is a direct secondary, and not higher order (not a grand-daughter particle of the parent). also it means that there are dominantly single secondaries produced, and not multiple at once, because then trackID 3 and 4 would also be more present.
so in conclusion: I would use the inspiration from example B1 to implement the accumulation of deposited energy per event, do that individually for electrons, positrons and photons, and at the end of each event, fill one of the three histograms/ntuples with the total energy deposit per particle type per event.
fyi: it does not matter how many particles you want to simulate. the number of primary particles per event should be 1 in this case. just create 6.25e8 events…
1 Like