Get energy deposit by processes

Hello, I am trying to determine the total energy deposition within a volume and to discern the specific contributions from the distinct interaction processes. But I do not understand my output and maybe I am missing some conceptual understanding.

I am using the UserSteppingAction to accumulate the energy depositions.
For the total energy deposit I do:

G4double deltaEdep = step->GetTotalEnergyDeposit();

and keep adding it to an inital value of 0. Pretty standard I think.

To get the energy deposition due to the specific processes I followed extended example OpNovice2:

const std::vector<const G4Track*>* secondaries = step->GetSecondaryInCurrentStep();

for(auto sec : *secondaries)
    G4String creator_process = sec->GetCreatorProcess()->GetProcessName();
    G4double edep = sec->GetKineticEnergy();

    detector->AddEnergyDepositToProcess(creator_process, edep);


Where AddEnergyDepositToProcess is just a method which adds the energy to a map where the key is the creator process.

When I run my simulation for one event I get the following output for the energies:

   detector: total energy: 340.759 MeV
Process: Cerenkov, Energy Deposit:   104.3 keV
Process: annihil, Energy Deposit: 15.2449 MeV
Process: compt, Energy Deposit: 82.7012 MeV
Process: conv, Energy Deposit: 1.47501 GeV
Process: eBrem, Energy Deposit: 1.88445 GeV
Process: eIoni, Energy Deposit: 42.4892 MeV
Process: phot, Energy Deposit: 23.3796 MeV

which does not make much sense to me. Maybe my code is wrong but before trying to fix it I want to make sure that this is the right approach or am I failing to understand something?

I think this approach makes sense. What part of your output does not make sense? The energy values? The processes?

So I thought the sum of the processes should be the total energy deposited? That is clearly not the case.

What energy are you expecting?

It is possible that your processes are overlapping. So electrons may deposit energy as eBrem, but then the generated X-rays will also deposit energy via eIoni

I am expecting the sum of all processes to be the total energy, ie. 340.759 MeV. How can the energy lost to eBrem possibly be greater than 1GeV? Sure, secondary electrons can deposit energy as eBrem and subsequent X-rays could ionize but in the end the sum of all processes should be the total energy.

Expecting as in, what is the initial energy of your particles? If you are shooting a single 100 MeV particle than clearly something is wrong. If you are shooting 10’s of 1 GeV particles then the energies aren’t off.

If your logic is adding the energy that is lost in the electron when creating the x-ray, then also counting the energy lost from ionization from the x-ray, the energy will be counted twice

It is only one event, i.e. a single electron with an energy of 800Mev.

Let me explain how I understand it:
In each step secondaries are created by a specific process, e.g. eBrem. So, X MeV are transferred to the x-rays in this example. I imagine that these x-rays also generate steps and transfer Y MeV to eIon. So X is added to my eBrem and Y to my eIon… Aaaaaah okay, as I am writing this I understand my fallacy and what you mean by your last comment. That was quite stupid.
So somehow I have to filter only the processes my primary particle undergoes, and that should sum up to the total energy lost, is that right?
It seems that something like this is done in extended/electromagnetic/TestEm18?

Yes, that’s basically correct. However…

You’re not going to be able to separate out the different processes that cause energy loss. The only process which is stored for each step is the discrete (“PostStepDoIt”) process which limited the step length.

Generally (at least for charged particles) the energy loss (G4Step::GetTotalEnergyDeposit()) is determined mainly by the continuous processes that are applied along the step (like dE/dx), combined with any energy that “would have” gone into secondaries that were suppressed by production cuts.

So your map above is very useful for telling you about the tracked secondaries, and how they originated. If you accumulate the initial kinetic energy of those secondaries, that will give you information about how the initial primary energy got distributed to secondaries.

Thanks for the answer, but I am confused about what you are trying to tell me regarding what is possible and what is not possible.

When I look at example TestEm18’s SteppingAction:

if (Edep > 0.) fEventaction->SumEnergyDeposited(trackID, Edep);


void EventAction::SumEnergyDeposited(G4int trackID, G4double edep)
  if (trackID == 1) fEdepPrimary  += edep;
  else fEdepSecondary += edep;

This groups the energy deposit into energy loss due to
a = interactions of the primary particle which is the energy continuously deposited along the primary track.
b = interaction of the secondary particles which is the energy continuously deposited along the secondary tracks.

In the example the total energy loss E_loss by the primary particle is calculated as
E_loss = a. + (energy transferred to secondaries). This tells me that a = (energyloss due to simple excitation) + (energy that would have gone into secondaries). And the latter is because of the default production cuts? If so, where can I read about these default cuts, as I have a hard time understanding why they are implemented. My current understanding ist that you need those to avoid computational effort (I do not know much about Monte-Carlo simulations). Please comment or correct me on what I said until here.

I won’t show the code but in the example it is furthermore demonstrated how to compute the energy transferred to secondaries due to a specific process by the primary particle. I guess THIS is what a I was asking for in my initial post!

Then in the end I get confused again, since in the example the total energy deposit is calculated as
E_tot = a + b. But I’d say that E_tot should be equal to the total energy loss of the primary particle, since a + b could be greater than E_loss which does not make sense regarding energy conservation.

So as you might see I am quite confused and maybe got things mixed up. I would appreciate some enlightenment maybe with example TestEm18 as a basis for explanation and some concepts explained in simple words as referring to the simulation objects and how they are interwined often does not help me as a beginner. Thank you very much!