Add muon decay interactor

[stognini]

This PR adds the muon decay interactor. Muons can only decay via one channel, either

$$\mu^- \longrightarrow e^- \overline{\nu}_e \nu_\mu$$

or

$$\mu^+ \longrightarrow e^+ \nu_e \overline{\nu}_\mu .$$

An important caveat is that the Geant4 physics manual defines the maximum possible sampled energy for the electron to be $$E_\text{max} = m_\mu / 2$$, while the source code defines $$E_\text{max} = m_\mu / 2 - m_e$$. This difference leads to a small missing energy, as the total energy of secondaries on a muon decay at rest is not equal to the muon rest energy in the case of the source code definition.

[stognini]

Sampling loop

The MuDecayInteractor uses two for loops with a hardcoded limit of 1000 iterations, instead of usual do {} while's, to sample the secondary energy. A test with 10k iterations leads to only up to 50 samplings, as the plot shows, being far below the 1000 threshold. I kept the sampling equivalent to Geant4, but we can move to a while loop if that is preferred.

$$E_\text{max}$$ definition difference

Using $$E_\text{max} = m_\mu / 2$$ (the physics manual one) I do recover the muon mass on a decay at rest (as shown in the MuDecayInteractor.test). Nevertheless, since our sampling limit changes, our sampled energies are slightly shifted to higher values for all 3 secondaries. Plots below show the shift for the electron. Top plot uses the definition used in the source code, while the bottom uses the one defined in the physics manual.

[github-actions]

Test summary

 3 385 files   5 208 suites   3m 53s ⏱️
 1 580 tests  1 552 ✅ 28 💤 0 ❌
17 413 runs  17 347 ✅ 66 💤 0 ❌

Results for commit 4eaff7c.

♻️ This comment has been updated with latest results.

[amandalund]

Looking great @stognini! Still have to go over most of this, but since we were talking about this point on the call, I guess my preference would be to just create the electron secondary until we have a use case for the neutrinos. The extra bit of computation won't be a big deal, but the two extra secondaries per decay (that we can't track anyway) could increase the secondary/initializer memory requirements a bit. And since you've already implemented it and we'll have it in the history, we can easily just revert the change when we need it. (I'll let @sethrj and @whokion weigh in, though.)

[whokion]

I do agree with @amandalund for now (unless we want to add neutrino-nucleus processes in a near-term future or want to send these neutrinos back to CPU so that Geant4 can handle them in the case that users add neutrino processes in their physics list, of which data transfer needs to be supported in a similar way when we add lepto-/photo-nuclear processes). Of course, please share any other users cases that you are considering.

Regarding to the maximum energy of the daughter electron from the muon decay, I think that the formulas in the P.M. is an approximation (i.e, assuming that secondaries are massless in the V-A calculation except the primary muon). The implementation in Geant4 may be more practically/experically correct. You may still recover the energy conservation in the rest frame of the muon following correct kinematics between (e + n_e) + n_muon and muon as necessary - I am surprised that Geant4 does not calculate them correctly. Anyway, I will also review and give more detailed comments later. Thanks for a good addition!

[whokion]

@stognini Please address the first round of review comments. Thanks.

[whokion]

@stognini Thanks for the update. Please see the second set of comments (which are mostly minors).

[whokion]

This may be safe without a loop check as the sampling is well-defined. Nevertheless, have you thought about combining two do-while loops into one (i.e., 2D sampling with a composition and rejection technique as their energies are correlated)?

[whokion]

As we are not creating neutrinos, we may need a comment for not calculating cos_theta between the decay electron and its neutrinos (i.e., electron neutrino direction w.r.t electron). As discussed, I guess that the electron direction can be simply sampled by
charged_lep_dir = from_spherical(UniformRealDist(-1, 1)(rng), UniformRealDist(0, 2 * constants::pi));
without EulerRotation.

[amandalund]

Looking good @stognini! A few small suggestions:

[amandalund]

The RejectionSampler could be used here (if the sampling is not refactored according to Soon's suggestion).

[stognini]

So this is a weird detail that happens in the sampler that doesn't allow me to use any of the standard helper classes (that have reasonable assumptions, such as asserting that the acceptance value is non-negative). Here is the problem:

The electron_nu_energy_frac range is [0, sample_max_), which is [0, 1.000023). Now, the f_ value set in the RejectionSampler is f_ = electron_nu_energy_frac * (1 - electron_nu_energy_frac). So that means that for the upper limit (sample_max_ = 1.000023), f_ < 0, which leads to an assertion error.

[sethrj]

I agree we should try to refactor this. Do you know the underlying PDFs to be sampled? It would be really good to explicitly document those as well.

[amandalund]

I also don't think sampling the energy fraction in [0, sample_max) makes sense... won't that be equivalent to sampling it in [0, 1), just slightly less efficient?

[amandalund]

Some thoughts on how this might be refactored: it looks like the inner rejection loop samples the $\nu_e$ energy fraction from the PDF $f(x) = 1/6 x (1 - x)$, $x \in [0, 1)$. This has the CDF $F(x) = 3x^2 - 2x^3$, so to use the inverse transform method we'd have to invert that cubic function.

This PDF is also the beta distribution $\text{B}(2, 2)$. To sample from a beta distribution $\text{B}(\alpha, \beta)$, we can draw two independent samples from a gamma distribution, $X \sim \Gamma(\alpha, \theta)$ and $Y \sim \Gamma(\beta, \theta)$: then $\frac{X}{X + Y} \sim \text{B}(\alpha, \beta)$.

Our GammaDistribution implements the algorithm described in Marsaglia and Tsang which uses a rejection method with a very high (>95%) acceptance rate, but also has nested loops and draw samples from a normal distribution. For this simple case where we'd just need samples from $\Gamma(2, 1)$, we could use this technique which would just require generating two exponentially distributed RVs: $-\sum_{k=1}^n \ln U_k \sim \Gamma(n, 1)$.

I think this would at least be more efficient than the current rejection method.

[sethrj]

@stognini can you add a decay section to doc/implementation and include the new decay class? I guess we'd want it to have its own section at the same level as EM/optical as opposed to a subsection of EM or living inside "core"?

[sethrj]

I think this should be a note ... no one is going to be doing EM physics and neutrino physics inside the same calculation since the cross sections are so absurdly different.

[stognini]

They won't be doing neutrino physics per se, but vertex quantities won't be conserved, and reconstruction/filtering for analyses might depend directly on verification of conservation of E and p to reconstruct a mother particle from an observed signal. So I thought it was reasonable to have a warning.

[sethrj]

Does Geant4 actually do this? Would the user have to have custom physics/stepper code to pull out the neutrinos and save their momentum and kill them?

[stognini]

Geant4 just returns all three particles. So if the user doesn't want the neutrinos to be transported and/or use their initial state, then yes, the user would likely have that coded in their G4VUserTrackingAction implementation.

[sethrj]

Maybe this function should just be electron momentum? Not passing in the mass (just used shared_.electron_mass) and instead passing in Energy{energy_frac * max_energy_}.

[stognini]

This was originally as is because it was the same function for all three secondaries. I kept it isolated for future refactoring if we need to add the neutrinos, but it could technically just be in the main interactor body, since it is just used once.

[sethrj]

If we're going to shoehorn the decay process into the regular "interactor", this is exactly equivalent to absorption, right? So let's just delete this and the from_decay.

And in fact, you need to do this, because the InteractionApplier tries to change the direction if
(result.action != Interaction::Action::absorbed), and it only kills the particle if it's absorbed.

So we should probably just delete this and use "from_absorption" . If you (or anyone else) think it's too confusing to mark it as absorbed, then just use

Suggested change

	decay, //!< N-body decay
	decay = absorbed, //!< N-body decay

and replace your function below with

CELER_FORCEINLINE_FUNCTION Interaction from_decay() { return from_absorption(); }

[stognini]

But we still end up with result.action == absorbed in this case, right? Because having it tagged as decay in our actions is important. For example: in a perfect vacuum, a particle will never undergo absorption, but it will always decay if not stable.

[sethrj]

Currently we don't store the interaction result in any way or provide a way for the user to interrogate it. It's only used to update the post-step data in the InteractionApplier. I think that's OK, since we should be able to use the post-step action ID for users to interrogate that sort of thing.

[sethrj]

Maybe we should change the name from absorbed to killed? But then that might be ambiguous with a particle that had an error and had to be killed.

[whokion]

Sorry that I used a wrong account for the earlier comment by gxbert (now deleted). Recovering/rewording it here (but shorten): It is confusing that decay = absorbed as the decay process is a material independent (in a sense that the particle never interacts with a material and never be physically absorbed). The action decay = absorbed will be ambiguous if the muon capture process is added later which can be considered as absorbed and physically totally different process even though multiple processes can be categorized to the same action. Nevertheless, a question (to @stognini ) is whether there is any specific treatment for decay products from others (for InteractionApplier) which may justify the addition.

[stognini]

My initial enum addition was a starting point to expand the InteractionApplier when the decay process class comes in, as well as its inclusion in the stepping loop on future PRs. So we would have a particle killed if it is absorbed or underwent decay. For now, our actions are internal, but I assume that at some point we will be able to provide access to actions to the end user, similar to querying processes in G4, for analysis purposes. So one could, for example, store only data from a given decay, if their focus is, say, understanding signals from a given channel.

[sethrj]

We can always query processes by looking at the post-step action ID, since each of those correlates to a model (decay or interaction) that the track undergoes. I don't think we'll want to expand that further. I think instead we may want to make the "interaction result" less about the physics and more about the programming implementation: whether to update the direction and/or kill the particle, e.g.

[stognini]

Fair point. Should we rename absorbed to a catch-all word, such as stopped? This would cover all cases: scatters, stops (absorption, decay), remains unchanged or fails during sampling.