Pulse height tally

[cfichtlscherer]

Hi everyone,

this is my first git pull request. Further, I am not really experienced in C++. I have done my best, but if I have made mistakes, I apologize.

For simulating gamma spectroscopy applications, a so-called pulse-height tally is needed. Here I would like to present the implementation of the pulse-height tally for photons in OpenMC. This tally calculates how much energy is deposited in a cell by an entire particle history (means a particle and its progenies).

The basic idea is that the particle gets a new attribute pht_storage_ which is an array. In this array, the deposited energy of the particle for every cell is stored. After every Compton scattering, pair production, and photoelectric event, or the death of the particle since it falls under the cutoff energy, this array is updated. It also needs to take care to remove the energy of created secondary particles in these reactions (to prevent double counting). If the particle and all its progenies are dead, we take the value of the index of the analyzed cell for the pulse-height tally value.

There are currently several limitations. It only runs in a single thread mode, can not be combined with other tallies, and can only tally a single volume per simulation.

Details of an extensive validation and information on the tally can be found in the attached paper.

We validated the results with an analytical model in which only simplified versions of the reactions are modeled. For reproducing these results, I attach the two files (physics.cpp and run_shuttleworth.py - since GitHub does not allow uploading these file endings, I changed them to txt). The physics.cpp file contains the simplified reaction versions. To run this version, it needs to be swapped with the physics.cpp file in OpenMC and after compiling, the script run_shuttleworth.py will produce the analytical results.

Further, I attach a small example that can model the spectrum of a Cs-137 point source in a simple NaI scintilation detector. When you compare these simulation results with experimental data, you realize that the experimental data is broadened. I have an additional python script for doing this, but we could also discuss implementing it directly into the OpenMC code.

Looking forward to your feedback.

Best Christopher

Implementation and Validation of the Pulse-Height Tally in OpenMC.pdf example_cs137.txt physics.txt run_shuttleworth.txt

[paulromano]

@cfichtlscherer Thanks so much for your contribution! Sorry I haven't had a chance to review this yet. I'll try to take a closer look in the next few days.

[cfichtlscherer]

@gridley, thanks a lot for all your comments! I try to go through them and create a new commit.

I would be very happy to make the pulse-height tally functionality a part of a future OpenMC version. I believe it could be interesting for many potential users.

[cfichtlscherer]

Here already the smaller changes and suggestions incorporated. Thanks again to @gridley. I think it's much better already. I need some time for the tally filtering system. But I am happy for all suggestions and comments.

[cfichtlscherer]

Dear Paul, thanks a lot for your reply and your comments. That I try to answer all below. I will go through the code review criteria again. Sorry if I made some mistakes here / did not follow the correct procedure.

• Regarding the Python API. Currently, you can use the tally via Python (as in the attached example). Would this be enough? Sorry, I have the feeling I miss something here.

tallies = openmc.Tallies()
tally = openmc.Tally(name='pulse-height tally')
tally.filters = [openmc.CellFilter(crystal), openmc.EnergyFilter(np.linspace(0, 2e6, 101)]
tally.scores = ['pht']
tallies.append(tally)
tallies.export_to_xml()

• I will create a regression test by simulating several pulse-height values in some settings.
• With documentation, do you mean something for the user guide and maybe a jupyter notebook example (in this example, we could also show how to broaden the spectrum afterward)?

[paulromano]

@cfichtlscherer Ah yes, now I see that it is usable from Python as is. Let's change the score name from 'pht' to 'pulse-height' just to be a little more consistent with other score names.

[cfichtlscherer]

I changed it. I think we should do it like the first version you suggested (that's how I currently have it implemented):

ph_tally = openmc.Tally()
ph_tally.filters = [openmc.CellFilter(...), openmc.EnergyFilter(...)]
ph_tally.scores = ['pulse-height']

I think its more consistent with the other tallys. If that's okay with you, I will mark this comment above as resolved?

[cfichtlscherer]

@shikhar413, Thanks a lot! I am actually very happy with an additional set of eyes to look over this code. I hope we can finish this soon. I am a bit busy at the moment, but I will try to answer you later in more detail.

[cfichtlscherer]

I am sorry, the last weeks were pretty busy, I am on vacation right now. From the 1st of November, I will have more time for working on this again. @shikhar413 I am sorry, but I will also go through all your comments beginning of November. Thanks already again! And I will get in contact with you. :-)

[Ebiwonjumi]

@cfichtlscherer Sorry for bringing up an issue related to the pht implementation. What I'm really interested in is the broadening of the tally spectrum as it applies to OpenMC flux tallies. You mentioned in your first post that you have python script for this. I want to ask is it that the python script broadens the final simulation result or the individual particle contributions to the tally? Most importantly, how can I use your implementation for gaussian energy broadening of a flux tally? not for pht. Please let me know if you are still working on this.

Best, Bamidele

[shimwell]

Just wondering if this PR has been replaced by #2452 and if we should close it. @cfichtlscherer

[cfichtlscherer]

Yes that's right. Thank you for remembering. I will close it.