Support leptonuclear and "extra" EM physics

[sethrj]

Celeritas doesn't currently support the models and processes implemented as "extra" physics. This shows up in #1756 as fewer steps with slightly less energy deposition, and it is non-negligible in ATLAS.

Suggested work plan for rpogressive capability:

• [ ] Log a warning if any processes that apply to Celeritas-managed particle types are ignored
• [ ] Load energy, material-dependent cross sections from processes that aren't implemented and combine into an "absorption" process to locally deposit energy (preserving total cross section of the particle type)
• [ ] Instead of just absorbing, reconstruct a G4Track from the current state (see https://github.com/celeritas-project/celeritas/issues/1894), get the interaction length for the G4 "extra" processes that apply, sample proportionally based on those, and call "DoIt" on the sampled process (?) and complete the rest of the G4 stepping loop on CPU, including hits and stacking G4 secondaries
• ~~Implement the process in Celeritas but offload the secondary to Geant4~~

Particle Type	Process	Model/Implementation	Notes
Gamma (γ)	Photonuclear	G4HadronInelasticProcess	Uses G4GammaNuclearXS
		G4LowEGammaNuclearModel	Low energy (< 200 MeV)
		G4CascadeInterface (Bertini)	Moderate energies
		G4TheoFSGenerator + G4QGSModel	High energies with string fragmentation
Electron (e⁻)	Electronuclear	G4ElectronNuclearProcess	Uses G4ElectroVDNuclearModel
Positron (e⁺)	Positronuclear	G4PositronNuclearProcess	Uses G4ElectroVDNuclearModel
Muon (μ⁺/μ⁻)	Muon-nuclear	G4MuonNuclearProcess	Uses G4MuonVDNuclearModel

[whokion]

The implementation of the associated lepto-nuclear models (including sampling of final states) should be outside the current scope, as it requires a substantial development effort. For example, adding the gamma–nuclear interaction of a Bertini-like cascade model also requires implementing intra-nuclear interactions for all other hadrons and the following de-excitation models — the same applies to QGS. As an alternative, the final state sampling could be performed on the CPU, while these additional physics processes compete for steps on the GPU with m.f.p. In such cases, the primary gamma or e⁺/e⁻ could be sent back to the CPU so that the corresponding DoIt is executed by Geant4. I will start to implement cross-section calculations for these processes on the GPU as hard-wired models, since the cross-section evaluations appear to be performed on the fly. This will serve as a necessary first step while we continue to explore better approaches or a systematic, long-term extension of physics support.

[sethrj]

Thanks @whokion , I've crossed out the last bullet point to reflect the feasibility of doing the sampling on GPU. It sounds like your alternate suggestion is the previous bullet point?

Regarding the cross sections: are they smooth/predictable enough that they could be precomputed and tabularized?

[whokion]

Yes, the last bullet point is likely the most practical approach. The cross sections for photo- and lepto-nuclear processes are parameterized per element (or per isotope for some light elements) as a function of incident energy, typically represented by tabulated (array) data. I’ll review G4HepEM to see whether it use any useful shortcuts (though I suspect it may not) and then outline a more detailed implementation plan before beginning the actual coding. My plan is to implement the parameterized cross sections as simply data structure as possible, without compromising accuracy, especially in the resonance regions. We can later consider alternative approaches, such as exploring AI/ML-based methods, once the initial implementation is in place if relevant.

[sethrj]

If it's simply array data, that's great. After you establish the requirements from the geant4 side let's meet to sketch out an implementation plan. I definitely want to take the opporunity to do this process differently:

1. Use inp models rather than ImportData to represent the processes we support; this is an easy way to also specify the cross sections as nonuniform grids (which is what we want).
2. Include generic nonuniform grid calculations as part of the physics calculator utility, so we can replace the "hardwire" models with tabulated lookups.
3. Use the same core classes to be used in #1894

[whokion]

1. Use `inp` models rather than `ImportData` to represent the processes we support; this is an easy way to also specify the cross sections as nonuniform grids (which is what we want).

These array data are not the cross sections themselves but are used to compute analytical or parameterized element cross sections at a given energy. In any case, nonuniform grids may still be useful.

[sethrj]

Aha, I see. That makes sense. Thanks!

[sethrj]

@whokion Looking at how G4HepEm does this for grammas: https://github.com/mnovak42/g4hepem/blob/0b1c7ecea142e5fa0867d39faaaf7f360f7fce24/G4HepEm/G4HepEmInit/src/G4HepEmGammaTableBuilder.cc#L123-L124

It seems G4HepEm is only approximating the cross sections by tabulating them, not actually calculating them on the fly. It simply adds them into the total. Likewise it tabulates on a grid for electrons:

https://github.com/mnovak42/g4hepem/blob/0b1c7ecea142e5fa0867d39faaaf7f360f7fce24/G4HepEm/G4HepEmInit/src/G4HepEmElectronTableBuilder.cc#L434-L498

Can we do this before trying on-the-fly?

Do we know in what regions the leptonuclear physics matters most, and is it possible for you to plot the cross sections as generated by Geant4 vs the ones used by G4HepEm? (I suppose the latter being just a subset of discrete points?) Maybe we should have, as a first step, a utility for loading cross sections from these processes?

Important: since the nuclear cross sections don't seem to be subject to secondary production cuts (there's no continuous processes?), the data is correlated with G4Material not G4MaterialCutsCouple, aka our GeoMatId not PhysMatId.

I think a starting point will be:

diff --git a/src/celeritas/inp/Physics.hh b/src/celeritas/inp/Physics.hh
index 89659979e..a27a59401 100644
--- a/src/celeritas/inp/Physics.hh
+++ b/src/celeritas/inp/Physics.hh
@@ -13,6 +13,7 @@
 #include "corecel/Types.hh"
 #include "corecel/io/Label.hh"
 #include "celeritas/Types.hh"
+#include "celeritas/inp/Grid.hh"
 #include "celeritas/phys/AtomicNumber.hh"
 
 #include "PhysicsProcess.hh"
@@ -95,6 +96,35 @@ struct OpticalPhysics
     }
 };
 
+/*!
+ * Store total cross sections for sending back to Geant4.
+ */
+struct NuclearOnloadProcess
+{
+    //! Total cross sections in log E for each GeoMatId
+    using MaterialGrid = std::vector<inp::UniformGrid>;
+
+    MaterialGrid xs;
+
+    std::string g4process;  //!< TODO: can make variant with integer ID
+
+    //! Whether gamma-nuclear physics is enabled
+    explicit operator bool() const { return !xs.empty(); }
+};
+
+struct EmHadronPhysics
+{
+    NuclearOnloadProcess gamma_nuclear;
+    NuclearOnloadProcess electron_nuclear;
+    NuclearOnloadProcess positron_nuclear;
+
+    //! Whether em-hadronic physics is enabled
+    explicit operator bool() const
+    {
+        return gamma_nuclear || electron_nuclear || positron_nuclear;
+    }
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Set up physics options.

[whokion]

@sethrj Exactly! As a first step, my plan was to visualize the photo- and lepto-nuclear cross sections from Geant4 and evaluate their smoothness to determine whether tabulating these cross sections is justified. These cross sections typically have four different regions (see [attached plot)

and are material-dependent. In principle, the G4HepEM approach may be still valid if the grid spacing for the tabulated cross sections is on the same order as the characteristic dimension of the cross-section parameters.

[sethrj]

Ah nice! That's a great reference.

[sethrj]

• Hit reconstruction will contain ActionId which can encode a process if needed (gamma vs electron)
• Call ApplyYourself which will create secondaries (with associated G4HadProjectile and G4Nucleus target), modify G4Track (possibly kill, possibly decrease energy: updated final state?)
• Use stacking capability to push the updated track if it's not absorbed
• Use stacking capability to push secondaries added to the track