Rare anomalous production of superheavy nuclei

[DavidWalz]

Exception in crpropa::ModuleList::run: 
nuclearMass: nucleus not found 86, 44

Exception in crpropa::ModuleList::run: 
nuclearMass: nucleus not found 114, 55

simulation setup

ModuleList
  SimplePropagation: Step size = 10 - 10000 kpc
  FutureRedshift: h0 = 0.673, omegaL = 0.685, omegaM = 0.315
  PhotoPionProduction: CMB
  PhotoPionProduction: IRB (Dominguez 2011)
  PhotoDisintegration: CMB
  PhotoDisintegration: IRB (Dominguez 2011)
  NuclearDecay
  ElectronPairProduction: CMB
  ElectronPairProduction: IRB (Dominguez 2011)
  ...
crpropa::ModuleList: Number of Threads: 1
Run ModuleList

This is reported by @davidwit to occur very rarely.

[DavidWalz]

with Gilmore EBL

Exception in crpropa::ModuleList::run: 
nuclearMass: nucleus not found 108, 53

[adundovi]

In one scenario I inject nuclei heavier than 56, thus, for this case, I needed to adjust the following line in ParticleMass.cpp:

if ((A < 1) or (A > 56) or (Z < 0) or (Z > 26) or (Z > A))
    return A*amu;

is there any reason why it shouldn't be a default behaviour? Instead of the error-throwing...

[DavidWalz]

We can extend nuclearMass and remove the error throw like you propose.
Please note however

1. decay and photodisintegration data is only provided for Z < 27, N < 31.
   PhototoDisintegration quietly skips, NuclearDecay needs a if ((Z > 26) or (N > 30)) return; in line

   CRPropa3/src/module/NuclearDecay.cpp

   Line 82 in 2f81224
   or it will treat cases Z < 26, N > 30 falsely.

2. I'd still like to solve the underlying problem

[adundovi]

@Froehliche-Kernschmelze this could be related to SOPHIA and #225.

[marcus-wirtz-snkeos]

Maybe it helps somehow: for me this error occured lately with helium propagation of an E^-1 spectrum above 10 EeV and an propagation length of 160 Mpc. All interactions were turned on:

sim = crp.ModuleList()
sim.add(crp.SimplePropagation(1*crp.kpc, 10*crp.Mpc))
sim.add(crp.Redshift())
sim.add(crp.PhotoPionProduction(crp.CMB))
sim.add(crp.PhotoPionProduction(crp.IRB))
sim.add(crp.PhotoDisintegration(crp.CMB))
sim.add(crp.PhotoDisintegration(crp.IRB))
sim.add(crp.NuclearDecay())
sim.add(crp.ElectronPairProduction(crp.CMB))
sim.add(crp.ElectronPairProduction(crp.IRB))
sim.add(crp.MinimumEnergy(10*crp.EeV))

The log:
[WARNING] nuclearMass: nuclear mass not found in the mass table for A = 66, Z = 32. Approximated value used A * amu - Z * m_e instead.
crpropa::ModuleList: Signal 15 (SIGINT/SIGTERM) received

For me it's a bit rare that this heavy nuclei can be produced by a helium candidate:
crp.SourceParticleType(crp.nucleusId(4, 2))