Ekin in EnergyDepositionSD gives kinetic energy after, not before first interaction

[jvavrek]

Jack Silano (LLNL) and I have found that in the EnergyDepositionSD class, the kinetic energy appears to be tallied incorrectly. Based on the documentation, Ekin should be the kinetic energy of the incident particle, but based on experimentation, the kinetic energy is that after the primary's first interaction.

For example, firing 2500 keV photons at the ZeroDegree detector, a verbose output of EnergyDepositionSD::ProcessHits() gives the following:

Time: 4.93493, PARTICLE: gamma, Process: compt, Track KE: 1783.58, Pre-step KE: 2500, Post-step KE: 1783.58, EDEP: 1.2228
Time: 4.9849, PARTICLE: gamma, Process: compt, Track KE: 314.978, Pre-step KE: 1783.58, Post-step KE: 314.978, EDEP: 0.0324
Time: 5.07972, PARTICLE: gamma, Process: Transportation, Track KE: 314.978, Pre-step KE: 314.978, Post-step KE: 314.978, EDEP: 0
Time: 4.98643, PARTICLE: e-, Process: eIoni, Track KE: 541.52, Pre-step KE: 1468.57, Post-step KE: 541.52, EDEP: 927.052
Time: 4.98683, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 541.52, Post-step KE: 0, EDEP: 541.52
Time: 4.93551, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 715.195, Post-step KE: 0, EDEP: 715.195
 Final Edep = 2185.02
 Final Ekin = 1783.58

Here we can see that the final value of Ekin = 1783.58 comes from the outgoing kinetic energy after the first Compton interaction at 2500 keV. Note also that this (quite often) causes Edep > Ekin.

A simple fix might be to change
hit->SetKineticEnergy(track->GetKineticEnergy());
to
hit->SetKineticEnergy(aStep->GetPreStepPoint()->GetKineticEnergy());

This gives a better picture of what the true incident kinetic energy is, but doesn't perfectly solve the Edep > Ekin problem, as it doesn't account for pileup within the same event—multiple secondary photons depositing energy in the detector, e.g.:

Time: 4.97272, PARTICLE: gamma, Process: compt, Track KE: 347.45, Pre-step KE: 759.784, Post-step KE: 347.45, deltaE: 412.334, EDEP: 0, CreatorProcess: primary
Time: 4.97661, PARTICLE: gamma, Process: compt, Track KE: 240.124, Pre-step KE: 347.45, Post-step KE: 240.124, deltaE: 107.326, EDEP: 1.2228, CreatorProcess: primary
Time: 5.00702, PARTICLE: gamma, Process: compt, Track KE: 216.2, Pre-step KE: 240.124, Post-step KE: 216.2, deltaE: 23.9242, EDEP: 11.067, CreatorProcess: primary
Time: 5.02753, PARTICLE: gamma, Process: phot, Track KE: 0, Pre-step KE: 216.2, Post-step KE: 0, deltaE: 216.2, EDEP: 11.067, CreatorProcess: primary
Time: 5.02764, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 205.133, Post-step KE: 0, deltaE: 205.133, EDEP: 205.133, CreatorProcess: phot
Time: 5.00702, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 12.8572, Post-step KE: 0, deltaE: 12.8572, EDEP: 12.8572, CreatorProcess: compt
Time: 4.97666, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 106.103, Post-step KE: 0, deltaE: 106.103, EDEP: 106.103, CreatorProcess: compt
Time: 4.973, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 412.334, Post-step KE: 0, deltaE: 412.334, EDEP: 412.334, CreatorProcess: compt
Time: 4.9251, PARTICLE: gamma, Process: phot, Track KE: 0, Pre-step KE: 72.2999, Post-step KE: 0, deltaE: 72.2999, EDEP: 11.067, CreatorProcess: eBrem
Time: 4.92513, PARTICLE: e-, Process: eIoni, Track KE: 0, Pre-step KE: 61.2329, Post-step KE: 0, deltaE: 61.2329, EDEP: 61.2329, CreatorProcess: phot
 Final Edep = 832.084
 Final Ekin = 759.784
 Edep 832.084 > 759.784 by 72.2999 !

Here two photons of energies 759.784 keV and 72.2999 keV both deposit their full energies in the detector, but only 759.784 keV is set as the value of Ekin.

Either way, I think the change recommended above is still a more intuitive measure of the 'Ekin' despite the pileup issue.

[u-eff-gee]

Dear jvavrek,

Thank you for pointing out this issue and for the thorough documentation.

Indeed, the readout of the kinetic energy was in contradiction to what is written in the README. I applied your proposed fix to all detector classes. Now, the value of Ekin should make more sense.

So far, I did not notice this issue, because I usually do not read out Ekin when using an EnergyDepositionSD. This is because the EnergyDepositionSD stores information only at the level of G4Events (analysisManager->FillNtupleDColumn() is only called at the time of EnergyDepositionSD::EndOfEvent. Compare this to ParticleSD, where the ROOT file is filled at the level of G4Steps). Of course, several different particles may contribute to the total energy deposition, which is the quantity of interest for the EnergyDepositionSD. However, it only stores information about the first step (at least now it does) of the first particle in the event, since the actual filling of the ROOT file is done at the end of the event. I figured that the Ekin, position, momentum or particle information that an EnergyDepositionSD gives you would not be very useful anyway.
This comes from my background in Nuclear Resonance Fluorescence with single-crystal HPGes, where the only information accessible at the end is the total energy deposition by all particles in an event.
The pileup effect of the EnergyDepositionSD which you mention is caused by the same limitations of the sensitive detector. It can also be reproduced by shooting multiple particles in the same event. If the AngularCorrelationGenerator is used with the following macro commands

##################################
# Define the first step
##################################

# /angcorr/particle initiates the creation of a new cascade step. 
# The macro will fail if this is not the first command in a step definition.
/angcorr/particle e+ 
/angcorr/energy 0.511 MeV
/angcorr/direction 1. 0. 0.
/angcorr/polarization 1. 0. 0.

##################################
# Define the second step
##################################

/angcorr/particle e-
/angcorr/energy 0.511 MeV

/angcorr/relativeangle 180. deg
/angcorr/polarization 1. 0. 0.

, an electron and a positron can be emitted in opposite directions inside a detector. They will be part of the same event. Similar to your pileup effect, the EnergyDepositionSD records

*    Row   *      edep *      ekin *  particle *    volume *
************************************************************
*        0 * 1.5329989 *     0.511 *        11 *         0 *

, because the electron is emitted first. Note that it records an energy deposition which is even higher than the sum of the kinetical energies of the two emitted particles, because the positron was annihilated.

Actually, no sensitive detector exists at the moment that can store all the information that tracking/verbose gives you at the G4Step level. If you want to simulate gamma-ray tracking inside a detector in the GRETINA/AGATA style, then this would not be possible. In principle, ParticleSD can read out information about every hit point, but at the moment it does that only for the first step of a particle inside the sensitive detector. Particles are discriminated by their track ID and event ID. To have a 'true' tracking detector, one would need to remove the

if (trackID != getCurrentTrackID() || eventID != getCurrentEventID())

condition in ParticleSD::ProcessHits.

[u-eff-gee]

I will try to make the limitations of the detector classes more clear in the documentation and implement a real tracking detector. For the NRF experiments usually performed at HIγS, this may be overpowered, but it is nice to have it for debugging.

[jvavrek]

@uga-uga thanks for the update. I agree that a full tracking detector would be overpowered for NRF. I've found that the initial kinetic energy can be a useful diagnostic for efficiencies and fluxes, though, so it's nice to have directly alongside the energy deposition.