Distance travelled D by the photons inconsistent.

[mcroisonnier]

Describe the bug

I simulate photons only as they propagate through the intergalactic medium, with EBL, CMB and Intergalactic magnetic field.

1. When I simulate an observer at the surface of a sphere with a radius D_source, some of the photons (not a lot, but still) in the resulting data frame have a distance travelled D_travelled lower than D_source. I don't understand how this is possible
2. I also observe another strange behaviour; Some of the photons which did not interact ( same SN0 and SN, and same E0 and E) are travelling a longer distance than D_source. I also don't get which scenario allows this.

Is the bug related to the physics in CRPropa?
not sure

To Reproduce
randomSeed = 42
Brms=1e-8nG
lMin = 1Mpc
lMax=25Mpc
sIndex=5./3.
turbSpectrum = SimpleTurbulenceSpectrum(Brms, lMin, lMax, sIndex)
gridprops = GridProperties(Vector3d(0), 100, 0.5Mpc) #origin, N of grid points, spacing # 100 x 0.5 = 50 Mpc
BField = SimpleGridTurbulence(turbSpectrum, gridprops, randomSeed)

dumpGrid(BField.getGrid(), 'myfield.dat')
dumpGridToTxt(BField.getGrid(), 'myfield.txt')

vgrid=Grid3f(gridprops)
loadGrid(vgrid, 'myfield.dat')
loadGridFromTxt(vgrid, 'myfield.txt')

dsrc = redshift2ComovingDistance(0.03) * u.m
dsrc_mpc = dsrc.to(u.Mpc).value

electrons = True
photons = True
thinning = 0.90

cmb = CMB()
ebl = IRB_Gilmore12()

propagation = PropagationCK(BField)
propagation.setMinimumStep(1e-6 *parsec)

sim = ModuleList()
sim.add(propagation)
sim.add(Redshift())
sim.add(ElectronPairProduction(CMB()))
sim.add(PhotoDisintegration(CMB()))
sim.add(EMPairProduction(cmb, electrons, thinning))
sim.add(EMPairProduction(ebl, electrons, thinning))
sim.add(EMInverseComptonScattering(cmb, photons, thinning))
sim.add(EMInverseComptonScattering(ebl, photons, thinning))
sim.add(MinimumEnergy(10 * GeV))
sim.add(MaximumTrajectoryLength(1.02*dsrc_mpc * Mpc))

extent = np.ceil(dsrc_mpc / 50 * Mpc)
sim.add(PeriodicBox(Vector3d(0), Vector3d(extent50Mpc)))

E = 100 #Energy TeV

obs1 = Observer()
obs1.add(ObserverSurface(Sphere(Vector3d(0.), dsrc_mpc* Mpc )))
out1 = TextOutput(f"Mrk421_B1e-17_E{E:.0e}TeV.txt", Output.Everything)
out1.setEnergyScale(eV)
out1.setLengthScale(kpc)
out1.enable(Output.CreatedPositionColumn)
out1.enable(Output.CreatedDirectionColumn)
out1.enable(Output.WeightColumn)
out1.enable(Output.RedshiftColumn)
out1.disable(Output.CandidateTagColumn)
obs1.onDetection(out1)
sim.add(obs1)

source = Source()
source.add(SourcePosition(Vector3d(0, 0, 0)))
source.add(SourceEmissionCone(Vector3d(1,0,0)Mpc,5np.pi /180) )
source.add(SourceParticleType(22))
source.add(SourceEnergy(E*TeV))

sim.run(source, 100000, True)
out1.close()

Expected behavior

1. In this case, I expect D_travelled to be equal to or superior than D_source.
2. In this case, I expect D_travelled to be equal to D_source

System (please complete the following information):

• Operating system: Ubuntu 22.04.4 LTS
• Laptop or computing cluster? laptop
• Python version: Python 3.10.17
• Swig version: SWIG Version 4.3.1
• Cmake version: cmake version 4.0.2

[lukasmerten]

From the current set-up it's hard to tell which of the several modules might cause a potential error. Please try to narrow it further down by, e.g., removing one module after the other until the potential issue disappears. For example see what happens if you remove the redshift evolution. Can you pin it down to one interaction etc. pp.

Otherwise:
2) Particles need to overshoot the boundary to be detected. So I would expect that even without any interactions, the trajectory length at detection is larger than the observer radius at least by a small bit.

1. Might be connected to the periodic boundary and how the secondary position is sample from the parent particle position.

PS:
You don't need ElectronPairPorduction and Photodisintegration in this purely leptonic scenario.