Document the differences between current LDMX Geant4 Bertini and upstream

[EinarElen]

The Bertini cascade is a crucial part of our simulations, we should make sure that the differences between what is in this repository and upstream Geant4 is clear to any users without requiring them to read through the Geant4 source code. The changes are detailed in Appendix A of https://arxiv.org/pdf/1808.05219.pdf but we know that some of these have been upstreamed into regular Geant4.

Some Bertini particle identifiers used below

• $p$: pro = 1
• $n$: neu = 2
• $\pi^+$: pip = 3
• $\pi^-$: pim = 5
• $\pi^0$: pi0 = 7
• $\gamma$: gam = 9
• $K^+$: kpl = 11
• $K^-$: kmi = 13
• $K^0$: k0 = 15
• $\bar{K}^0$: k0b = 17

[EinarElen]

Changes between LDMX 10.2.3 and Geant4 10.2.3 (NU = Not upstreamed), (LP = Not upstreamed but can be patched in ldmx-sw), (UP = Upstreamed)
From v10.2.3...LDMX.10.2.3_v0.5

• Use g4fsqrt instead of fsqrt in G4tgrEvaluator (irrelevant for Bertini)
• LP Change the maximum Bertini Cascade energy to 10 GeV and minimum QGS energy to 8 GeV
  This is something we can do directly in SimCore as part of https://github.com/LDMX-Software/SimCore/blob/trunk/include/SimCore/GammaPhysics.h
• NU Add a $\gamma + p / N \rightarrow \pi + p/n$ backscatter distribution to the G4TwoBodyAngularDist option which is applied when

// In 
// Final state multiplicity has been determined to be 2 in G4CascadeFinalStateAlgorithm::Configure 
// Otherwise, fs =0 
// In G4CascadeFinalStateAlgorithm::ChooseGenerators
    if((is == gam*pro || is == gam*neu) &&
       (kinds[0]==pro || kinds[0]==neu) && 
       (kinds[1]==pip || kinds[1]==pim || kinds[1]==pi0))
      kw = -1;   // NT: backscatter case

// In G4TwoBodyAngularDist::ChooseDist  
if (((is == gam*pro || is == gam*neu) &&
       (fs == pro*pim || pro*pi0 || neu*pi0 || neu*pip)) && kw== -1) {
    return gp_npiback;
  } 

Here, is is the initial state multiplied (bullet * target), fs is the final state multiplied, with kinds being the two resulting particles. So we are applying the new distribution when the initial state is $\gamma + p/n$ and the final state is $p + \pi^-$, $p + \pi^0$, $n + \pi^0$, or $n + \pi^+$ with the nucleon as the first final state particle (i.e. the nucleon is backscattering)

• Change $\gamma + N$ and $\gamma + P$ data bins in source/processes/hadronic/models/cascade/cascade/include/G4CascadeGamPChannel.hh and source/processes/hadronic/models/cascade/cascade/include/G4CascadeGamNChannel.hh to accout for the new two-body final states. The second part of the template here is the number of two-body final states

typedef G4CascadeData<30,6,6,4,5,6,7,7,7> data_t;
typedef G4CascadeData<30,8,6,4,5,6,7,7,7> data_t; 

... 
// G4CascadeData 
template <int NE,int N2,int N3,int N4,int N5,int N6,int N7,int N8=0,int N9=0>
  const G4int (&x2bfs)[N2][2];		// Initialized from file-scope inputs

The cross-section for this backscatter case is

  // p pi0 ( = n pi0 ) backscatter -- from integral of Anderson 1969 data.  Start at 2.4 GeV 
   {0.0,    0.0,    0.0,    0.0,    0.0,    0.0,   0.0,    0.0,    0.0,    0.0,
    0.0,    0.0,    0.0,    0.0,    0.0,    0.0,   0.0,    0.0,    0.0,    0.0,
   1.5e-4, 6.4e-5, 2.85e-5, 1.2e-5, 5.0e-6 , 2.1e-6, 9.6e-7,  3.6e-7,  1.5e-7, 6.4e-8}, 
  // n pi+ ( = p pi- ) backscatter -- from integral of Anderson 1969 data.  Start at 2.4 GeV
   {0.0,    0.0,    0.0,    0.0,    0.0,    0.0,   0.0,    0.0,    0.0,    0.0,
    0.0,    0.0,    0.0,    0.0,    0.0,    0.0,   0.0,    0.0,    0.0,    0.0,
   1.5e-4, 6.4e-5, 2.85e-5, 1.2e-5, 5.0e-6 , 2.1e-6, 9.6e-7,  3.6e-7,  1.5e-7, 6.4e-8}, 

• Add the two new final state possibilities to the list of possible outgoing particles

This change would be straight-forward to add to upstream Geant4, but might have been introduced elsewhere.

• UP Change the order of some $\gamma + N$ 2-body final states to have consistent kinematic distribution. This was a bug in Geant4 where the final nucleon was getting the kinematics of the other secondary in reactions like $\gamma + p \rightarrow p+ \pi^0$

• UP* Updated $\gamma + (NN)$ (quasi-deutron) cross-sections
  These were significantly overestimated in G4 10.2 for GeV-scale gammas and underetimated for lower. This change has been upstreamed in Geant4, but with a less conservative estimate

LDMX description:

• up to 0.56: match Kossov 2002 (10.1140/epja/i2002-10008-x) g1+g2
• above 0.56: based on exclusive deuteron dissociation data from CLAS Mirazita et al 2004 (10.1103/PhysRevC.70.014005), which covers ~0.5 to 3 GeV.
• 0.56 to 1.8: use highest central-value of differential cross-section from CLAS, in any angular bin, multiplied by 4pi. This overestimates TOTAL cross-section by a factor of 2-10 to avoid under-producing on the tails.
• 1.8 to 18.0: a rough fit to the 1-2 GeV CLAS data: (delta_t/s^8 * 10000 mb), extrapolated. Falls slower than the expected s^11/t scaling. Note like above, overestimates total xsec by a significant factor.

Geant4 11.2.0

• Quasideuteron absorption cross section is taken to be the deuteron photo-disintegration cross section ($\gamma + (pn) \rightarrow p + n$)
• Values taken from smooth curve drawn through data from 2.4 - 500 MeV,
• plus angle integration of JLab data from 0.5 - 3 GeV M. Mirazita et al., Phys. Rev. C 70, 014005 (2004)
• Points above 3 GeV are extrapolated.

• UP Add a potential thickness parameter to the cascade parameters and the G4NucleiModel, default to 1 fm. This is included in upstream Geant4, but is not runtime configurable like it is in LDMX's version.

The thickness is used to calculate $q_{perp}$ in addition to $q_{v}$ for deciding whether or not to reflect at a density boundary. $q_{perp} = \frac{2 \cdot p_{perp}^2 \cdot \xi}{r}$ with $\xi$ as the thickness and $p_{perp} = p^2 - p_r^2$. $q_v = dv^2 + 2 \cdot dv \cdot E + p_r^2$

Where $dv$ is the potential difference between the two regions that the particle is attempting move through (corresponding to the height of the wall seen by the particle at the boundary).

Boundary transition then carries out as follows

• If $q_v \leq 0, q_v + q_{perp} \leq 0$: Reflect particle

• If $q_v > 0$ Transition through the barrier

• Otherwise, perform transition by transverse KE and adjust $p_{perp}$ after the transition

• Change the dinucleon density in the G4NucleiModel based on https://arxiv.org/pdf/nucl-th/0301091.pdf, corresponding to a reduction in dinucleon density by ~1/7 and scale down the diproton/dineutron densities by a factor 1/2.

The effective number of quasi-deutrons, $P_{QD} \approx L \frac{Z(A-Z)}{A}$. Assuming that the density of quasi-dineutrons and diprotons are similar, we should expect that $P_{NN} = P_{PP} = \frac{1}{2}L\frac{Z(A-Z)}{A}$. We calculate the QD density as follows

• For each nuclear density zone $z$
  • Calculate number of QDs in this zone as $N_{QD,z} = V_{z} V_{Z} \rho_{p,z} \rho_{n,z}$
  • Sum these to get the total number of naive QDs $N_{QD}$
• Scale the density by $R_{NN} = \frac{N_{QD, LDA}}{N_{QD}} = \frac{L_{LDA} \frac{Z \cdot(Z-A) }{A} }{N_QD} $
  $P_{QD} = R_{NN} \rho_p \rho_n$
  $P_{nn} = \frac{R_{NN}}{2} \rho_n \rho_n$
  $P_{pp} = \frac{R_{NN}}{2} \rho_p \rho_p$

Below is the output from an example

 zone 0 : 1815.61 naive quasideuterons
 zone 1 : 599.455 naive quasideuterons
 zone 2 : 85.1433 naive quasideuterons
 zone 3 : 57.9201 naive quasideuterons
 zone 4 : 12.9922 naive quasideuterons
 zone 5 : 3.46688 naive quasideuterons
 total:  2574.59 naive quasideuterons
 vs. Levinger L = 9.11306 => 400.173 quasideuterons
 => dinucleon numbers being rescaled by 0.155432

• NU PN Algorithm changes
  See PN-Fix Methodology Changes Hadron Interactions #9

[EinarElen]

The following is a listing of all non-trivial changes to Bertini between 10.2.3 and latest

240 files changed, 19257 insertions(+), 4764 deletions(-)
What fun

Potentially relevant

• New collider for gamma on light targets (H, D, T, 3He) added in 2019 would show up in some of our PN events
  Irrelevant from LDMX perspective (~5% of EcalPN PN reactions are with hydrogen)
• Minor bugfix in G4NumIntTwoBodyAngDest, angular distribution bins were not zeroed out correctly (bug 2515)
• Major bugfix In G4TwoBodyAngDist, charge exchange was checking for the wrong type of reaction (bug 2516). This meant that the charge exchange interaction pi+ + n -> pi0 + p would not get the right angular distribution (if it would have gotten any)

if ((is == pim*pro && fs == pi0*neu) || (is == pip*neu && fs == pi0*pip) ||

should have been

if ((is == pim*pro && fs == pi0*neu) || (is == pip*neu && fs == pi0*pro) ||

• Major bugfix, LDMX's two-body final state swap bug has been corrected in several other final states than just gamma + - Pion-absorption bug in 10.2.3 Geant4 Bertini cascade #11 has been fixed
• Several channels have received updated many-body final states (especially strange pair-production channels)
  • p + p, p + n, n + n, pi0 + N, pi- + n, pi- + p, pi+ + n, pi+ + p have received updates to 6 to 9-body final states
  • K- + n, K- + p, K0b + n, K0b + p, K0 + n, and K0+ p Kaon induced reaction interactions (except K+?) have gotten 8/9 body final states added
• Minor bug in pi0 + N channel for 8-body final states, wasn't preserving quantum numbers
  Irrelevant changes
• The ABLA de-excitation module can be used to de-excite
• CascadeCoalescence removed clusterhas/triedclusters in order to reduce memory usage, instead handle each combination explicitly
• Bugfix in PreCompound interface related to internal conversion (+ a unit error), we don't use PreCompound