Support decay in stepping loop #1991

amandalund · 2025-09-27T14:54:39Z

This adds initial support for decay in the stepping loop. Muon decay is currently the only decay channel supported, but adding additional channels should be straightforward now.

To easily support the different decay channel types (which may not have fixed daughter particle types, or may have an arbitrary number of daughters), this moves from using channel-specific data (e.g. MuDecayData) to the generic DecayChannelData and DecayTableData. Daughter particle IDs are passed to the interactor and are expected to be in a specific order (this should be validated in the DecayChannel constructor).

Changes include:

Add decay input structs
Import decay tables from Geant4 for particles that have the decay process enabled
Rename Process to InteractionProcess and add a new DecayProcess (used for all particle types that decay)
Add DecayChannel which extends from CoreStepActionInterface, responsible for validating decay table input and launching kernels for interactors
Add MuDecayChannel derived class
Refactor MuDecayInteractor to use generic decay data
Store decay tables and mapping from channel ID to action ID in physics data
Include decay as a step limiter, calculating cross section from decay constant
Sample decay process from among other discrete processes in select_discrete_interaction
If decay is selected, sample decay channel according to branching ratios

I haven't added any tests yet (only checked that running celer-g4 on a simple problem with muon decay imports the data and interacts by decay without any failures), but otherwise this should be ready for a first look.

See #1375

This reverts commit 757e9a1541c1efa4a118f62eca7cc22f8305781c.

- Rename Process to InteractionProcess and add DecayProcess - Add DecayChannel abstract base class - Add MuDecayChannel and executor - Add decay tables and channel - to - action mapping to physics data - Add DecayTableInserter helper class to build decay tables in physics data - Add decay to particle - processes in physics if present - Add decay_length method to ParticleTrackView - Add decay limiter in calc_physics_step_limit - Store decay process xs and sample from among other interaction processes - Sample decay channel using branching ratios if decay process is chosen - Cache sampled channel for accessing decay daughters in interactor - Refactor muon decay data and interactor to use generic decay data

github-actions · 2025-09-27T15:10:17Z

Test summary

4 576 files 7 305 suites 14m 23s ⏱️
1 993 tests 1 975 ✅ 18 💤 0 ❌
25 456 runs 25 361 ✅ 95 💤 0 ❌

Results for commit 08a2f1f.

♻️ This comment has been updated with latest results.

sethrj

I'm still working through the PR. It looks like currently we have one step action per decay channel per particle type? I think we'd be better off treating the decay physics input like a collection of models starting with those we support (as we do with the new optical surface physics), rather than the old ImportX style of making everything generic with an enum to tie it together. Each model would provide a set of particles it applies to.

So for muon decay, we might have optional<MuDecayChannel> in the DecayPhysics. We wouldn't need daughters since each decay channel knows what it can decay and how the particle decays. (If we really wanted to allow more granularity, we could have the MuDecayChannel input take a set of parent particles, which normally would be {mu-,mu+} or empty.) A more complex/general decay channel such as a three-body decay (subset of G4PhaseSpaceDecayChannel functionality) would have a different set of data underneath it: a vector of {particle -> daughters} etc. User-specified decay channel input would be a generic function that returns a vector of decay channels (or something).

As for branching ratios, it might make sense to give them as an unnormalized value or as a decay constant (1/time), which would get normalized among all the decay channels for a particular particle type. For the current capabilities, we only have one decay channel so the ratio doesn't even matter yet.

The DecayProcess would then construct the different decay channels, categorize them by particle ID, and determine branching ratios by relative decay constants/values. The decay channel itself would change (I think?) so that it might contain data from multiple particle types rather than a single one. For the typical case of particle/antiparticle (muon, tau) this is sensible because the only difference in the physics is the daughter particle IDs (I think?).

I'm going to keep reviewing, and I know that this suggestion (assuming I'm not missing something important) is a major change, so we don't have to implement it as part of this PR. But if we do want to go in that direction, we could at least delete aspects that aren't yet needed and will be deleted in the future (I'm thinking specifically the decay channel enum).

sethrj · 2025-09-27T16:01:11Z

src/celeritas/phys/ParticleTrackView.hh

+    CELER_EXPECT(this->mass() > zero_quantity());
+    CELER_EXPECT(this->decay_constant() != constants::stable_decay_constant);


Since momentum is always nonzero, I think we just get a (correct) value of infinity if either or both of these is zero.

sethrj · 2025-09-27T16:14:39Z

src/celeritas/ext/GeantImporter.cc

+        // Exclude the neutrino daughters for muon decay
+        //! \todo Remove once neutrinos are included in \c MuDecayInteractor
+        size_type num_daughters = channel.type == DecayChannelType::muon
+                                      ? 1
+                                      : g4_channel->GetNumberOfDaughters();


Should we warn about missing neutrinos? (And is there another process in which we're neglecting neutrinos as well, or am I thinking of the muon decay that Stefano already put in?) Also, are we sure the neutrino is always listed last? Maybe it'd be better test the PDG inside the loop and warn+skip if neutrino.

I think it's the right thing to omit them on the "input" side rather than pull them in and ignore them at the process level, but maybe it'd be a nice feature to be able to store a more complete physics as input have have the PhysicsParams class trim out deactivated/unsupported features, e.g. removing hadronic physics if hadron particles are disabled?

sethrj · 2025-09-28T18:19:10Z

src/celeritas/decay/DecayData.hh

+
+//---------------------------------------------------------------------------//
+/*!
+ * Decay channels for a particle type.


Add a note that these are stored by the PhysicsParamsData and used in the ProcessGroup (processes that apply to a particular particle type).

sethrj · 2025-09-28T18:37:14Z

src/celeritas/inp/Physics.hh

+ */
+struct DecayPhysics
+{
+    using DecayTable = std::vector<DecayChannel>;


I think VecDecayChannel since it's not really a table: there's no [row, column] access, and each vector of particle decay channels has additional data.

Does the order of this vector matter, or should it be a set<DecayChannel>?

I think I see where @amandalund is coming from with this name, since the PDG book has decay tables where the 1st column is the decay channel (so daughters here) and the 2nd column is the branching ratio. But Vec is definitely better in the code.

About the vector ordering: sorting channels by BR may be more efficient when selecting them. Specially in cases like pi0 -> 2gamma, where the BR is ~99%.

sethrj · 2025-09-29T10:28:18Z

src/celeritas/decay/DecayProcess.hh

+
+//---------------------------------------------------------------------------//
+/*!
+ * Process for decay.


A little more documentation here would be good: this single class manages all possible particle decay channels, each of which is analogous to a model for the action process.

sethrj · 2025-09-29T10:31:17Z

src/celeritas/inp/Physics.hh

+    DecayChannelType type{DecayChannelType::size_};
+    double branching_ratio{};
+    std::vector<PDGNumber> daughters;
+};


I'm not 100% sure I understand the input here. For muon decay, don't we know what the daughter types are? (Is this like specifying the secondary particle types for compton scattering?) Will other decay types need more data than just the branching ratios?

sethrj · 2025-09-29T10:34:34Z

src/celeritas/ext/GeantImporter.cc

            case G4ProcessType::fOptical:
                return (which & DataSelection::optical);
+            case G4ProcessType::fDecay:
+                return (which & (DataSelection::em | DataSelection::hadron));


Could you add a test to the GeantImporter to see what decay models it pulls out?

sethrj · 2025-09-29T10:48:48Z

src/celeritas/decay/DecayProcess.cc

+                {
+                    case DecayChannelType::muon:
+                        result.push_back(std::make_shared<MuDecayChannel>(
+                            *start_id++, input_));


Does this mean that for muons we'd create two separate step actions, one for each particle type?

whokion · 2025-09-29T15:27:16Z

@amandalund I am mostly off this week, but I can catch up on this MR later if a second review is needed.

amandalund · 2025-09-29T16:09:45Z

Thanks @sethrj for the comments! I'll look over the rest soon, but just to clarify it's one step action per channel, not per channel per particle type (I don't think we'd want that!)

@whokion no worries, enjoy the time off! If you have a chance when you're back I'd definitely appreciate your review to make sure I've got the physics right.

whokion

Overall, this looks reasonable to me, although the current implementation is somewhat over-engineered if we’re only supporting muon decay for now - of course, there is nothing wrong with being general from the start. Separately from this MR (perhaps by @stognini), we should add all daughters of the muon decay if we claim full support for it (in the context of the neutrino physics).

whokion · 2025-10-02T17:37:53Z

src/celeritas/phys/PhysicsStepUtils.hh

+            // Calculate the macroscopic cross section as the inverse of the
+            // particle's decay length
+            process_xs = 1 / particle.decay_length();


I assume this is a technical tweak treating decay as part of a macroscopic interaction. In reality, decay is independent of particle–matter interactions and is an intrinsic process governed solely by the particle’s proper time. A better approach would be to let the interaction process and the decay process compete independently, with the final proposed step length whichever occurs first. I understood this was part of the motivation for introducing InteractionProcess and DecayProcess separately (along with other variations). Mixing these two intrinsically different processes through a material-dependent macroscopic cross section is somewhat physically unnatural, since the proper time is independent of the material.

From a computational perspective I think it's better to combine the probabilities of a change in particle state into a single random number sample. But you do bring up a very good point about length versus time: is there any way this representation will interfere with assumptions about changing speed over the step length? Does Geant4's decay process factor in an average speed, or does it assume the speed is constant with the starting speed? Can we improve up on that? (Or must we avoid improvements to better agree with existing results?)

In Geant4, the particle velocity is assumed to remain constant during decay-in-flight. This approximation is generally valid, particularly for the muon weak decay, since the mean decay length of a high-energy muon is relatively long (on the order of $\sim 2.2 \mu {\rm sec} \cdot c$). As a result, a muon rarely decays inside dense materials, where particle-material interactions would typically limit the step length first. If the decay occurs quickly, the associated energy loss over that step is to be negligible—especially for muons, which lose only minimum ionizing energy—since it is much smaller than their (eloss) range in this case. Thus, muon decay is primarily relevant in sparse media, where the muon velocity can be treated as approximately constant until decay. The same assumption also holds for short-lived particles, of which lifetimes are so short that they decay essentially immediately in flight, before significant energy loss or material interactions occur.

In principle, your point remains valid in situations where particle decay competes with material interactions on comparable length scales. However, such cases are not typical for common unstable particles in HEP detector simulations and, when they do arise, of course, should be handled with additional care even though I have not seen any good example.

The application for muon decay being funded is muon catalyzed fusion which uses cold muons. So this may be an important effect we can model more accurately than geant4?

Not sure—probably @stognini can clarify. My assumption is that it would be undesirable for muons to produce Michel electrons during the catalyzed fusion process.

whokion · 2025-10-02T17:55:22Z

src/celeritas/phys/PhysicsTrackView.hh

+ * correspondence between the channel ID and the action ID. Multiple channel
+ * IDs can map to a single channel action.


Could you clarify how multiple channel IDs map to a single channel action? Does this imply that the action is only depending on the multiplicity of daughters?

stognini · 2025-10-13T19:13:10Z

src/celeritas/phys/detail/DecayTableInserter.cc

+    {
+        CELER_LOG(warning)
+            << "branching ratios for decay channels should sum to 1 "
+               "but instead sum tp "


Suggested change

"but instead sum tp "

"but instead sum to "

amandalund added 5 commits September 27, 2025 10:21

Import decay tables

2cf22fe

Refactor to store decay channels in particle input

8f9a328

Revert "Refactor to store decay channels in particle input"

ffa73cc

This reverts commit 757e9a1541c1efa4a118f62eca7cc22f8305781c.

Remove unsupported decay channel types

9a467b1

amandalund requested review from sethrj, stognini and whokion September 27, 2025 14:54

amandalund added enhancement New feature or request physics Particles, processes, and stepping algorithms labels Sep 27, 2025

amandalund added 2 commits September 27, 2025 19:13

Fix docs

c30626c

clang-tidy

08a2f1f

sethrj reviewed Sep 29, 2025

View reviewed changes

whokion reviewed Oct 2, 2025

View reviewed changes

stognini reviewed Oct 13, 2025

View reviewed changes

		CELER_EXPECT(this->mass() > zero_quantity());
		CELER_EXPECT(this->decay_constant() != constants::stable_decay_constant);

		* correspondence between the channel ID and the action ID. Multiple channel
		* IDs can map to a single channel action.

Support decay in stepping loop #1991

Are you sure you want to change the base?

Support decay in stepping loop #1991

Uh oh!

Conversation

amandalund commented Sep 27, 2025

Uh oh!

github-actions bot commented Sep 27, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Test summary

Uh oh!

sethrj left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

whokion commented Sep 29, 2025

Uh oh!

amandalund commented Sep 29, 2025

Uh oh!

whokion left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

github-actions bot commented Sep 27, 2025 •

edited

Loading