Energy Conservation violation

akhilsadam · May 14, 2020, 5:12pm

I am firing 40MeV protons into a skull and brain phantom.

My energy deposition and secondary energies are lower than the primary energy of the proton for a specific interaction, where the proton seems to lose some energy.
Please view the following images.

Interaction : (38.107MeV) p + volume = (5.18MeV) p + n + 4 gammas + C11 + B11 + nu_e + volume + ~4MeV of energy deposit.

The netfraction plot shows the net change in energy between every step (the 40MeV peaks are due to new protons firing). Clearly the dip at around 879 is incorrect.
This plot is generated in my SteppingAction, where if the primary is a proton, I log the (change in primary energy + the energy deposit + the kinetic energy of all generated secondaries) every step.

The 5.18 MeV subpeak is due to my code misclassifying the scattered proton in the interaction as another primary.

The fractional plot is simply the (change in primary energy every step), also logged at the same place in the SteppingAction.

The red areas indicated normal decay trails, and the blue areas indicate new protons. The purple indicates the conservation violating interaction, while the black is the scattered proton and its decay (since my code considers it a new primary).

This is my log (also in Stepping Action):

The color scheme is the same as before.
Ep is primary energy, Es is secondary energy, and Edep is the GetTotalDeposit().
Clearly (in the purple area) the secondaries do not add up to ~ 38.107 MeV.
The black is the scattered proton.

I suspect my Physics List, so I have attached it for your reference.

Thank you.

akhilsadam · May 15, 2020, 3:09pm

Anyone?

Any sort of advice would be helpful.

akhilsadam · May 19, 2020, 2:55pm

Anyone?

Any sort of advice would be helpful.

donahuw2 · May 20, 2020, 5:39pm

I would like to state that I am not an expert here, and I am writing this as I find concerns. I think one thing that is confusing as a whole is your output format. Is that a list for each proton? Do the IDs correspond to a particular track or a step in the track? I assume all this data looks at the same track.

I am looking at this and there are other things that bother me about your output that you say is good but it is not. Looking at the very first, lines of your output we see that there was an energy Dep of 0.8, but the energy changes by 0.02.

I would ask are you pulling the data from the pre-step or post-step point? The data here is different and depending on what information you want, you need to select the correct one.

Are you sure that the TotalEnergyDeposit includes the non-ionizing component? I don’t know and I really don’t want to dig through the code base to find out. You should though (https://geant4.kek.jp/lxr/source/)

If you are worried about the physics list I would suggest running your simulation with just the G4EmStandardPhysics_opt4 or one of the preconstructed Hadron lists. I think that this would give you an idea if you are stacking anything strange on your physics lists.

One other thought, are you are missing the remaining energy of the primary proton. The primary proton may have an energy of 38 MeV, but are you sure it is fully stopping in the medium? is it possible you messed up the units on your geometry so it is a much shorter target than initially thought?

akhilsadam · May 20, 2020, 6:25pm

Every step in every track that has a proton primary. The SteppingAction is threadlocked, so each track is calculated sequentially.

This is since I am pulling data from pre-step point, so the energy change from the previous (unlisted) primary to the current primary deposits 0.8MeV. The Edep refers to the previous change.

Yes. I haven’t looked at the codebase, but this source mentions TotalEnergyDeposit = ionizing + non-ionizing.

github.com/ActarSimGroup/Actarsim

GEANT4 : GetEnergyDeposit

opened 05:07PM - 02 Dec 15 UTC

piotrkon

For getting the energy deposit in the gas we use now the function GetDeltaEnergy…Deposit(). It is indeed the energy difference between the two step point. The function called before was GetTotalEnergyDeposit() which gives the total energy deposit in the gas for the tracked particle. The problem of this previous function is that it doesn't take into account the contribution of the delta electrons produce to the energy deposit in the gas. Since we use this energy to calculate the number of electrons drifting to the pad plane, I'm wondering is it a good method to estimate the energy deposit in the gas? Actually there is another function called GetNonIonizingEnergyDeposit() which gives the energy deposit in the gas not by ionization. So what is recommended to get the energy deposit by ionization is GetTotalEnergyDeposit()-GetNonIonizingEnergyDeposit(). I looked at the contribution of the NonIonizingEnergyDeposit and also the difference between total and delta energy deposit.

I will try this and post if it does not fix the issue.

I am pretty sure that this is not the case, as I can see protons stopping using the visualizer, and I have a set of energy deposition plots.
As far as exiting the volume, it may be that these protons exit the world with this energy (since I do not have any volume based cuts), but that would be a set of crazy coincidences, since it would only happen via this interaction.

donahuw2 · May 20, 2020, 7:02pm

So if the data we are looking at is all of the steps for a single proton in your console output screen. Then: Step ID 874 resulted in a loss of 0.799 MeV
Step ID 875: the proton starts with an energy of 38.121 MeV -> but according to your explanation it should be 37.323 MeV (38.122-0.799)

I think that this is highlighting that you are not pulling the appropriate values from the steps based on what my understanding is.

akhilsadam · May 20, 2020, 7:23pm

No it is not all of the steps! I apologize if I am being unclear.
As in the graph, between every 40 MeV spike it is the same track, and all steps for that track.
My screenshot only shows part of the log, with the first few steps being a single track (the one in red)
Then the interaction (in purple), and then another proton (generated as a secondary in the previous interaction, but counted as a primary since it is a proton).

But for this particular set of steps it is the same proton.

Step ID 873: (not listed) had (38.122MeV + 0.799 MeV)
874: 38.122MeV
873: 38.121MeV.

My energy deposition lists the energy deposited in the previous step! (so the total energy at that step is primary energy + edep (listed below, but of previous step) + energy of all secondaries (listed below, but again of previous step).

The below is an old run (note the primary is still 40 MeV - the log used to output total energy (938+40) instead of kinetic energy)

akhilsadam · May 22, 2020, 3:23am

I used the Hadr04 list, and did not have any violations. So, technically problem solved.
However, I can no longer see any positrons being produced - at all! Is there a better physics list I can use?

civanch · May 25, 2020, 1:26pm

Hello,

Hadr04 is a specific example where only neutron physics HP is activated. If you have interest to study full physics - use Hadr01.

Check on energy conservation is not a trivial one - you need to take into account several processes, for example, target atom de-excitation. Also note, that HP neutron models by default do not conserve energy event by event, only in average. If any non-HP reference Physics List is used then energy conservation is fine. It is well tested for each modification of the code.

VI

akhilsadam · May 25, 2020, 1:29pm

Thank you, Mr. Ivantchenko.

I will look into the Hadr01 physics list and see if that solves the issue.