Pulse height tally #1881

cfichtlscherer · 2021-09-17T20:18:38Z

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

Update recognized thermal scattering materials

paulromano · 2021-09-21T02:39:51Z

@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.

gridley

OK, so, every nuclear engineer out there using OpenMC will agree this is an incredibly useful piece of functionality. Maybe it's why people are still hanging on to ol' MCNP, even. There are a few things that have to be addressed before merging this.

A lot of the tallying goes totally outside the normal tally system in OpenMC. The fact that this only allows filtering based on cell IDs is a real shortcoming. Ideally, you should be able to use any of the filter bin structures provided by OpenMC. They're actually not terribly hard to use, if you look at some of the other tally types.
This code will give unexpected results if the distribcell (considering Particle::cell_instance()) approach is used. But, that goes in the same area of complaint as not using the tally filtering system.

So, generally, I think a lot of this is in the right direction of what we'd want in OpenMC. Treating the energy deposition tally as a special case definitely makes sense, and your approach of storing the bins as a particle-private variable is definitely the way to go. The thing to change here, though, will be figuring out how to define these with the capability of using most of OpenMC's filters, or, at least defining this tally score such that only a very strict subset of the available filters is allowed.

I think you are not far from the ideal solution. @paulromano will have a better feel for the state of functionality here required before merging. Hopefully these comments will get you on the way towards that, but know it will almost surely take a few iterations to get this merged.

include/openmc/particle.h

gridley · 2021-09-21T03:24:27Z

include/openmc/simulation.h

@@ -47,6 +47,8 @@ extern const RegularMesh* ufs_mesh;
 extern vector<double> k_generation;
 extern vector<int64_t> work_index;

+extern std::vector<double> bins_pht;


You shouldn't make these extern declarations twice. It appears you should only have them in simulation.h. If you need them in settings.cpp, you should include simulation.h there.

I think that's part of the problem: my implementation doesn't yet work properly with OpenMC's tally functions. I think I need to work on how the tally output is determined in the end from the pht_storage array. Then we don't need to define external std::vector<double> bins_pht; nor external int cell_pht; here.
I think I was in too much of a hurry here overall and really wanted to create a reasonably running version.

gridley · 2021-09-21T03:33:50Z

src/particle.cpp

+void Particle::energy_delivert_in_cell()
+{
+  // Update pht_storage array
+  pht_storage()[coord(0).cell] += E_last() - E();


These new functions appear to be very short, and you could redesign the code so that you only need to call this once. You call it right before the return in the various branches of sample_photon_reaction branches. Rather than calling it at the end of each branch, can't you instead just call it after calling sample_photon_reaction?

Or, rather, I think since this would in principle be the same for each particle type, we'd probably want to call this in the event_collide function right around the score_collision_tally function. The problem here is that we would be assuming we only lose energy at collisions, and if we ever do charged particle transport correctly, that would be the wrong assumption. But maybe that bridge could be crossed at a later date with a cautionary TODO statement.

Also, the comment "update the pht_storage array" is not helpful at all because anyone with a basic knowledge of C++ can immediately tell that's what this line of code is doing. A more helpful comment would be something along the lines of "tally particle energy decrement to top-level universe cell". One thing I notice here is that this will probably not work with the distribcell approach (which would require taking cell_instance() into account in your pht_storage array).

Several comments here:

I like the idea of calling it only once after sample_photon_reactions a lot. Thanks for that idea.
I placed this function as you suggested in Particle::event_collide() in line 257 directly after the score_collision_tally. Since we don't want the electrons to account into this value, I added another if condition:

if (type() == ParticleType::photon){energy_delivered_in_cell();}

I also added a TODO. I think we should "cross this bridge later".

when you say these functions appear to be very short, do you think we still should work here with functions or simply call

pht_storage()[coord(0).cell] += E_last() - E();

?

I changed the comment to
// Adds the energy the particle loses during the collision.
but we can also use yours.

It currently will not work with the distribcell approach, which could be interesting for users. My thought was that in a first version this tally only works with the CellFilter and EnergyFilter. What do you think?

gridley · 2021-09-21T03:35:41Z

src/particle.cpp

+void Particle::killed_particle_energy_delivert()
+{
+  // Add the energy of the killed particles to the pht_storage array
+  pht_storage()[coord(0).cell] += E();


Similarly, we should probably not have this as a separate function and instead put it in Particle::event_death. Moreover, this will also not work with tallies based on cell_instance.

killed_particle_energy_delivered() is now called at the end of Particle::event_death. But we could just remove the entire function and simply call the single line

pht_storage()[coord(0).cell] += E();

Again, I need to think a bit on how to make this work with the cell_instance feature.

src/tallies/tally.cpp

gridley · 2021-09-21T03:48:54Z

src/tallies/tally.cpp

+              << std::endl;
+  }
+  if (exists == 1 && scores_.size() > 1) {
+    std::cout << "Error: The Pulse-Height Tally can currently not be used in "


OK, well, this somewhat resolves my previous comment about pulse heights. I don't think we should merge this if putting this limitation in, because this limitation is not due to any fundamental shortcoming of openmc but rather lack of time put into programming a design that can work with OpenMC's diverse, extensible tallying system.

Yes that's part of the "tally problem". I need to pass the pht_storage to the tally quite differently, I'll try to tackle that in the next few days, think I need a bit of time for that.

src/tallies/tally.cpp

src/tallies/tally_scoring.cpp

cfichtlscherer · 2021-09-21T08:51:08Z

@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.

…ge interaction.

… into pulse-height-tally

cfichtlscherer · 2021-09-21T15:34:25Z

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.

paulromano

Thanks again @cfichtlscherer for starting this PR. If you haven't already, please take a look at our code review criteria. For a PR with a new feature to be accepted, you'll need:

Changes in the Python API to be able to use the new feature
Accompanying tests (at the very least, a regression test for this capability) so we can ensure the feature isn't broken over time
Documentation of new user input (Python) options

Let us know if you have any questions about this and we'll do our best to point you in the right direction.

include/openmc/constants.h

paulromano · 2021-09-21T22:16:08Z

src/tallies/filter_energy.cpp

@@ -18,6 +19,7 @@ namespace openmc {
 void EnergyFilter::from_xml(pugi::xml_node node)
 {
  auto bins = get_node_array<double>(node, "bins");
+  simulation::bins_pht = bins; //! store them as a global variable for the pht


Agree with @gridley. I'd encourage you to think about the user interface for this feature a bit since that will influence how data is stored on the C++ side. For example, would the following be appropriate?

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

It may be the case that pulse-height tallies are different enough that they deserve their own class, in which case it might look something like:

ph_tally = openmc.PulseHeightTally() ph_tally.filters = [openmc.CellFilter(...), openmc.EnergyFilter(...)]

I'm not sure offhand whether this latter idea is necessary. We may have to think a bit more carefully about the overall design of different tally types.

src/tallies/tally.cpp

paulromano · 2021-09-21T22:21:03Z

src/tallies/tally.cpp

+  if (exists && !settings::photon_transport) {
+      fatal_error("Error: The Pulse-Height Tally works only with photon transport True.");
+  }


How difficult will it be to extend this capability to neutrons?

I have thought about how useful this feature would be. To simulate neutron detectors, I think the absorption tally is enough. Are there other applications? Would users be interested in this? I could imagine that fission makes the whole thing more complicated. But if you like, we can work on it. Maybe in an update when we have implemented the tally for photons?

There may be users interested in this, but thinking about it a bit, it does seem that it will be a little more complicated for a number of reasons so perhaps we should defer for now. Note that checking settings::photon_transport is actually not sufficient to determine that the simulation doesn't contain neutrons. settings::photon_transport is true when we run coupled neutron-photon simulations too. The most robust way would be to look through each external source distribution openmc::model::external_sources and see if it emits photons. However, I don't see a way to do that easily since we have different types of external sources (IndependentSource, FileSource, CustomSourceWrapper)... let me think more about how to handle that

cfichtlscherer · 2021-09-22T06:42:09Z

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 · 2021-09-22T20:21:31Z

@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 · 2021-09-22T20:50:01Z

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?

… into pulse-height-tally

shikhar413

@cfichtlscherer thank you for submitting this PR, I did some preliminary reading on PHT's earlier this year and I've finally gotten around to testing out your code. I hope you don't mind an additional set of eyes to look over your code.

It is really nice to see that the functionality for PHT's is actually not all that complicated, and it is exciting to see preliminary results for this code's potential use cases. After thinking about the design for this code, it might make sense to have a src/tallies/pulse_height_tally.cpp file that is included by src/tallies/tally.cpp where you store any PHT-related variables like pht_storage_, bins_pht_, cell_pht, etc. and also conduct any filter/score error-checking to make sure that the tally is defined properly by the user. That way, you can still define PulseHeightTally as a Tally instance but make it clear that it behaves differently from your typical tally.

Since there also aren't very many custom definitions in your implementation, I think it would also be helpful to store a boolean variable is_pht_simulation or something along those lines to wrap all of your newly defined functions in order to make it clear where this new functionality diverges from typical tally behavior. Feel free to reach out to me directly if you need any help with C++/OpenMC-specific questions

shikhar413 · 2021-10-12T05:38:51Z

src/physics.cpp

@@ -55,7 +55,6 @@ void collision(Particle& p)
    sample_positron_reaction(p);
    break;
  }
-


Undo this line

shikhar413 · 2021-10-12T05:43:59Z

src/tallies/filter_cell.cpp

@@ -5,6 +5,7 @@
 #include "openmc/capi.h"
 #include "openmc/cell.h"
 #include "openmc/error.h"
+#include "openmc/simulation.h"


shikhar413 · 2021-10-12T06:03:28Z

src/tallies/tally.cpp

+  bool exists = std::find(std::begin(scores_), std::end(scores_),
+                  TallyScore::SCORE_PULSE_HEIGHT) != std::end(scores_);
+  if (exists && scores_.size() > 1) {
+    fatal_error("Error: The Pulse-Height Tally can currently not be used in combination with other tallys.");


By this logic, this checks whether more than one score is defined for a single pulse height tally, so the error message should be "The Pulse Height Tally cannot currently be used with multiple scores"

shikhar413 · 2021-10-12T06:12:54Z

src/tallies/tally.cpp

+  if (exists && scores_.size() > 1) {
+    fatal_error("Error: The Pulse-Height Tally can currently not be used in combination with other tallys.");
+  }
+  if (exists && !settings::photon_transport) {


Add checks to make sure each filter is of type cell or energy, are there any other filters that might be relevant? In that case, I also think there should be a check to make sure there is at least one energy and cell filter defined in the tally.

One design choice would also be to enforce exactly one cell filter and one energy filter for this tally type, and users can define multiple pulse height tallies if there are multiple detectors in the system

shikhar413 · 2021-10-12T06:21:14Z

src/simulation.cpp

+        if (tally.scores_[0] == -17) {
+          score_pht_tally(p);
+        }
+      }


This loop could be optimized by storing the PHT id as pht_tally_id and having a boolean is_pht_simulation as global variables so that you aren't iterating through each tally in the system after each photon death

if (is_pht_simulation) score_pht_tally(pht_tally_id, p)

shikhar413 · 2021-10-12T06:29:45Z

include/openmc/settings.h

@@ -116,6 +116,8 @@ extern int trigger_batch_interval; //!< Batch interval for triggers
 extern "C" int verbosity;          //!< How verbose to make output
 extern double weight_cutoff;       //!< Weight cutoff for Russian roulette
 extern double weight_survive;      //!< Survival weight after Russian roulette
+extern std::vector<double> bins_pht; //! Used for the pulse-height tally
+extern int cell_pht;                 //! Used for the pulse-height tally


Where exactly is cell_pht initialized?

shikhar413 · 2021-10-12T08:32:58Z

src/particle.cpp

+    }
+
+    // Set birth cell attribute
+    if (cell_born() == C_NONE)


@cfichtlscherer, upon further inspection I think this logic can be simplified a bit since the code will only enter this block if a new source particle from the secondary bank is populated:

void Particle::remove_energy_of_secondary() { if (!exhaustive_find_cell(*this)) { mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); return; } // Set birth cell attribute cell_born() = coord(n_coord() - 1).cell; if (type() == ParticleType::photon && E() > settings::energy_cutoff[1]) pht_storage()[cell_born()] -= E(); }

shikhar413 · 2021-10-12T08:58:29Z

src/particle.cpp

@@ -351,6 +360,41 @@ void Particle::event_death()
    int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
    simulation::progeny_per_particle[offset] = n_progeny();
  }
+killed_particle_energy_delivered();


shikhar413 · 2021-10-12T09:35:20Z

src/particle.cpp

+  // TODO: The assumption behind the modelling of the pulse-height tally is that
+  // particles lose energy only in collisions, for a more comprehensive modeling
+  // of charged particles this has to be changed. 
+  if (type() == ParticleType::photon){energy_delivered_in_cell();} // score the pht only for photons


According to this logic, secondary photons that start with an energy less than the threshold will still be scored in this function but their initial energy does not get scored in remove_energy_of_secondary. Is this correct? Should the condition instead be:

if (type() == ParticleType::photon && E_last() > settings::energy_cutoff[1]) energy_delivered_in_cell();

cfichtlscherer · 2021-10-12T10:30:56Z

@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 · 2021-10-19T09:09:36Z

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 · 2022-09-09T08:35:40Z

@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 · 2023-04-12T13:00:29Z

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

cfichtlscherer · 2023-04-12T13:10:29Z

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

pshriwise and others added 2 commits September 17, 2021 13:20

Merge pull request openmc-dev#1877 from paulromano/decay-thermal-fix

974dd9e

Update recognized thermal scattering materials

Merge branch 'openmc-dev:develop' into pulse-height-tally

f826f19

gridley requested changes Sep 21, 2021

View reviewed changes

cfichtlscherer added 2 commits September 21, 2021 17:25

Most suggestions from @gridley included. Still open: tally, pht_stora…

60dbd36

…ge interaction.

Merge branch 'pulse-height-tally' of github.com:cfichtlscherer/openmc…

fd00996

… into pulse-height-tally

cfichtlscherer requested a review from nelsonag as a code owner September 21, 2021 15:26

cfichtlscherer requested a review from gridley September 21, 2021 16:56

paulromano requested changes Sep 21, 2021

View reviewed changes

minor bugs, 1001 -> -27

bfc0941

python input name pht -> pulse-height

7b1269b

cfichtlscherer added 6 commits September 23, 2021 08:14

python input name pht -> pulse-height

9d8f244

Pulse-Height from ReactionType=-27 -> TallyScore=-17

16c1bdd

Pulse-Height from ReactionType=-27 -> TallyScore=-17

03785e0

Merge branch 'pulse-height-tally' of github.com:cfichtlscherer/openmc…

5f2ef89

… into pulse-height-tally

regression test / small bugs

a3236d7

works with multiple batches

7328a4a

shikhar413 reviewed Oct 12, 2021

View reviewed changes

cfichtlscherer mentioned this pull request Dec 23, 2021

Tally filter for time #1932

Closed

cfichtlscherer added 3 commits February 23, 2022 08:49

search the pht error start

9e66d14

secondary if contidion kill energy

a78a374

update

de24a93

Merge branch 'openmc-dev:develop' into pulse-height-tally

137f829

cfichtlscherer mentioned this pull request Apr 3, 2023

Pulse height tally for photons #2452

Merged

cfichtlscherer closed this Apr 12, 2023

Pulse height tally #1881

Pulse height tally #1881

Conversation

cfichtlscherer commented Sep 17, 2021

paulromano commented Sep 21, 2021

gridley left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cfichtlscherer Sep 21, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cfichtlscherer commented Sep 21, 2021

cfichtlscherer commented Sep 21, 2021

paulromano left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cfichtlscherer commented Sep 22, 2021

paulromano commented Sep 22, 2021

cfichtlscherer commented Sep 22, 2021 • edited Loading

shikhar413 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cfichtlscherer commented Oct 12, 2021

cfichtlscherer commented Oct 19, 2021

Ebiwonjumi commented Sep 9, 2022

shimwell commented Apr 12, 2023

cfichtlscherer commented Apr 12, 2023

gridley left a comment •

edited

Loading

cfichtlscherer Sep 21, 2021 •

edited

Loading

cfichtlscherer commented Sep 22, 2021 •

edited

Loading