Geant4 Failing to Match with MCNP6 and NIST for Simple Test of Energy Deposited on a Cesium Iodide Cube

jpjpjp4747 · 2020-05-19T18:42:23.691Z

Thanks for your observation! I am using Geant 10.2. Perhaps there is some sort of physical process that must be manually enabled that is causing this flattening in the curve.

Aleksei · 2020-05-20T15:04:40.717Z

Hi again! Maybe there should be included scintillation process to CsI material to get good comparision? Because by default scintillation process is not taken into account.

jpjpjp4747 · 2020-05-20T17:01:27.596Z

Hi, I think that is certainly a possibility. I will investigate this over the next few days, and post my findings here. Thanks for the suggestion!

jpjpjp4747 · 2020-05-21T15:33:10.562Z

Hi Aleksei,

I wanted to point out that I have tried conducting this same simulation but with lead, and am receiving similar results, where the energy deposited per photon at higher energies is still not linear. This leads me to believe that the issue may not be scintillation, and I am also aware that scintillation in Geant is not the most straight-forward process to implement. Lastly, I wanted to point out that the plot which was obtained in MCNP uses an F6 tally, which is a method that measures energy deposited based on collisional heating, which I do not think includes scintillation.

Aleksei · 2020-05-21T17:11:39.426Z

Hi! I agree with you. Yesterday I icluded scintillation process in that simulation and doposited energy from scintillation was very small and I think it is not worth to take into account this praction of energy. Could you please show the plot from MCNP of deposited energy in case of using 1 cudic centimeter lead to compare that with geant4?

jpjpjp4747 · 2020-05-25T16:32:19.033Z

Hi Aleksei, apologies for the delayed response. Unfortunately, I have no experience using MCNP. The MCNP plot which I am comparing with Geant was created a few years ago by another lab member. It may be worth noting that using the online NIST database for various attenuation coefficients, one can make a plot of the energy deposited per photon as a function of photon energy using the NIST coefficients for Cesium Iodide. I have attached one that I have made myself below. Note, E_0 is the photon energy, u is the mass-energy absorption coefficient from NIST at that energy, and x = 1 since I am interested in a cube of 1 cm thickness. NIST only has these coefficients up to 20 MeV on their website, so I have only gone up to 20 MeV, though it still seems that the “correct” distribution is in fact a linear one, at least for CsI. Furthermore, I recently read a paper comparing Geant and MCNP which stated that MCNP internally uses NIST coefficients for some of its energy calculations, which would explain why it matches nicely with the linear distribution. I am confident that Geant should be able to match this since it is used so widely, I just cannot seem to figure out what is missing.
nistplot

jpjpjp4747 · 2020-05-28T02:12:27.892Z

For anyone interested, i found this paper which matches Geant4 data to NIST coefficients successfully even at higher energies, but using Geant 6.2.

https://ieeexplore.ieee.org/abstract/document/1495783?casa_token=Lf-LmxdP-ecAAAAA:6MUzIwOKHXi-4yd6wtpVaX6OFDDVFtBVqnHy-_TSRhyn6oiisBZv7_sqJPibhqGuNXOK-Bph

jrmadsen · 2020-06-07T01:52:21.881Z

Try using the physics list constructor for fiber grained control of which em physics list is used instead of the pre-defined physics list. Here is an example of using the physics list constructor: https://gitlab.cern.ch/geant4/geant4/-/blob/master/examples/extended/parallel/ThreadsafeScorers/src/TSPhysicsList.cc#L105

Also, make sure you are adding the G4StepLimiter to a volume and enabling it in the physics list – in order for it to actually be applied, you have to do both.

jpjpjp4747 · 2020-06-07T22:00:13.055Z

Thanks so much for the suggestion. I will try implementing this and let you know if it solves my issue.

Aleksei · 2020-06-12T11:19:36.364Z

Hi! I start to think that there is a some bug in the program. Because this process is so common that it should work in any physics list. Maybe it worth to open a bug report? Especially because you found the article where simulation data and experemental have a good agreement. Or I hope someone from developers will tell what wrong.

jpjpjp4747 · 2020-06-12T16:02:31.732Z

Hi jrmadsen, I tried implementing your physics list in my program, using the step limiter as you have and using only the EMStandardphysics constructor, and yet I still am not getting as much energy deposited on the CsI cube per photon as I should be compared with MCNP6, I am still stuck at about 1.6 MeV/ph which is about what I was at with my pre-defined physics. I very much appreciate your help, I am thinking of reaching out to a developer.

jpjpjp4747 · 2020-06-12T17:41:44.476Z

Hi Aleksei, as per your suggestion, I have submitted a bug report. I will update this post once they respond.

jrmadsen · 2020-06-15T17:32:12.500Z

Vladimir and Daren are the maintainers for the EM physics. They probably can provide you with some insight. @civanch @dsawkey

jpjpjp4747 · 2020-06-15T18:22:36.391Z

Thank you. Hopefully they will respond here.

dsawkey · 2020-06-16T14:10:33.179Z

I don’t think the first plot shows what you think it does. The detector volume is only 1 cc and many secondaries will leave the volume. Values of the mass energy absorption coefficient assume charge secondaries deposit energy in the material.

Try some numbers. At 50 MeV NIST shows mu en/rho = 0.04 cm2/g and mu/rho = 0.07 cm2/g. With rho = 4.5 g/cm2 this gives mu_en = 0.18 /cm. With dx = 1 cm and E = 50 MeV,

dE = - mu_en E dx = 9 MeV

In your first (MCNP) plot, dE at 50 MeV is 9 MeV. Thus, the MCNP is not taking account any secondaries leaving the detector. Geant4 is telling you there are secondaries leaving.

jpjpjp4747 · 2020-06-16T18:06:51.263Z

Thank you so much for your explanation, it is tremendously appreciated. In response to your insight, I have added a very high production cut to gammas to try and increase the energy deposited in my Cesium Iodide cube as much as possible and was able to obtain the following plot. valplotwithcuts Fortunately it is modestly linear, though this is still off from the original plot in MCNP by about a factor of 1/2 at most data points. I will continue investigating setting cuts and see if I can further get the MeV deposited/photon to increase.

wwols · 2020-06-18T15:27:21.329Z

Hi jpjpjp4747,

I agree with dsaweky, the first plot from MCNP6 is really suspicious.

dsawkey:

I don’t think the first plot shows what you think it does. The detector volume is only 1 cc and many secondaries will leave the volume. Values of the mass energy absorption coefficient assume charge secondaries deposit energy in the material.

It is a bit optimistic to expect that a 1 cm3 crystal will be able to contain an electromagnetic shower at 150 MeV. Ok, it’s just impossible. The Geant4 plot is much more logical than the linear dependence obtained with MCNP. At low energies, <500 keV the photon interactions are dominated by the photoelectric effect and most of events will result in a full energy deposition in a 1 cm3 CsI. When you move up with energy, Compton scattering starts to dominate. At 1 MeV you will have mostly partial energy deposition due to Compton scattering. Additionally, CsI is not very high density material, it won’t have an impressive photofraction in general (compared to PWO, BGO, LSO etc.). At around 10 MeV pair production kicks in, and your energy deposition per photon increases (and this is exactly what your G4 plot shows). Above 20 MeV things start to level out because a small 1 cm3 crystal cannot contain an electromagnetic shower which develops.

You can make your plot more linear by making CsI crystal larger. An infinitely large crystal will eventually stop any energy photon and you will get your linear dependence. I would advise to return to the standard em physics list and increase the crystal size to 1 meter.

Cheers,
Weronika

jpjpjp4747 · 2020-06-18T18:55:28.754Z

Hi Weronika,

Thank you for your thorough explanation. I have been looking into it, and it seems like the MCNP “tallies” which use collisional heating to measure energy deposition, such as the F6 and *F8 tallies, assume that “All energy transferred to electrons is assumed to be deposited locally” (MCNP Manual). This is in agreement with dsawkey’s statement that MCNP, as well as the mass energy absorption coefficients from NIST, make this same assumption (that no secondary charged particles leave the volume, they deposit all of their energy locally). If this is unphysical as you describe, do you have any ideas as to why people using MCNP would use this assumption for their dose calculations?

Best regards,

Joseph

wwols · 2020-06-18T20:14:44.819Z

Hi Joseph,

It’s difficult to guess why they did it, if I don’t know anything about the application. Maybe they were not interested in real energy deposition in the crystal, maybe it was used for some physics tests or calibration. Maybe detector was supposed to work in coincidence and something else was important (triggering, coincidences etc.).

But you have to remember, that you will have energy leakage not only from secondary electrons, but also from secondary photons (x-ray fluorescence, 511 keV from positron annihilation), and finally, you will still have primary photon leakage. In Compton scattering your primary photon will escape the crystal (with reduced energy of course), and only some part of energy will be deposited as a secondary electron. So even if you will force all electrons to deposit energy locally in the scintillator, you still have neutral particles escaping (photons). If you will force everything to deposit energy locally, then what’s the point of the simulation?

Cheers,
Weronika

civanch · 2020-06-22T07:48:03.540Z

Hello all,

Thank you Weronika for this clear explanation. I only would suggest to check $G4INSTALL/examples/extended/electromagnetic/TestEm0, which allows check cross sections and attenuation coefficients without real simulation and compare results with NIST tables.

VI