How to 4D simulate photons?

[k4zuyoh1]

Hello,

I started working with CRPropa recently and I am having trouble trying to 4D simulate photons, as my output file comes empty.

I basicaly am using the 4D simulation code with Photon Propagation. I've also tried using PhotonOutput1D, with my output yet coming empty.

I'm running the code as follows:

from crpropa import *
import crpropa

# set up random turbulent field
turbSpectrum = SimpleTurbulenceSpectrum(1 * nG, 60 * kpc,800 * kpc,5./3.)
gridprops = GridProperties(Vector3d(0), 256, 30 * kpc)
Bfield = SimpleGridTurbulence(turbSpectrum, gridprops, 42)

# simulation setup
sim = ModuleList()
sim.add(PropagationCK(Bfield))
sim.add(Redshift())
sim.add(PhotoPionProduction(CMB()))
sim.add(PhotoPionProduction(IRB_Kneiske04()))
sim.add(PhotoDisintegration(CMB()))
sim.add(PhotoDisintegration(IRB_Kneiske04()))
sim.add(ElectronPairProduction(CMB()))
sim.add(ElectronPairProduction(IRB_Kneiske04()))
sim.add( EMPairProduction(CMB(), True) )
sim.add( EMPairProduction(IRB_Kneiske04(), True) )
sim.add( EMDoublePairProduction(CMB(), True) )
sim.add( EMDoublePairProduction(IRB_Kneiske04(), True) )
sim.add( EMInverseComptonScattering(IRB_Kneiske04(), True) )
sim.add( EMInverseComptonScattering(CMB(), True) )
sim.add( EMTripletPairProduction(CMB(), True) )
sim.add( EMTripletPairProduction(IRB_Kneiske04(), True) )
sim.add(NuclearDecay())
sim.add(MinimumEnergy(1 * EeV))
sim.add(MinimumRedshift(-0.1))

# periodic boundaries
extent = 256 * 30 * kpc  # size of the magnetic field grid
sim.add(PeriodicBox(Vector3d(-extent), Vector3d(2 * extent)))

# define the observer
obs1 = Observer()
obs1.add(ObserverSurface( Sphere(Vector3d(0.), 0.5 * Mpc)))
obs1.add( ObserverPhotonVeto() ) 
obs1.add( ObserverNeutrinoVeto() )
obs1.add(ObserverRedshiftWindow(-0.1, 0.1))
output = TextOutput('output.txt', Output.Event3D)
output.enable(output.RedshiftColumn)
obs1.onDetection(output)
sim.add(obs1)

# Definir observador de gamas 
obs2 = Observer()
obs2.add(ObserverSurface( Sphere(Vector3d(0.), 0.5 * Mpc))) 
obs2.add( ObserverDetectAll() ) 
obs2.add(ObserverInactiveVeto())
obs2.add( ObserverElectronVeto() )
obs2.add( ObserverNeutrinoVeto() )
obs2.add( ObserverNucleusVeto() ) 
obs2.add(ObserverRedshiftWindow(-0.1, 0.1))
out2 = TextOutput('output-gamma.txt', Output.Event1D)
out2.enable(Output.CreatedIdColumn) 
out2.enable(Output.CreatedEnergyColumn)
out2.enable(Output.CreatedPositionColumn)
out2.enable(Output.RedshiftColumn)
obs2.onDetection( out2 )
sim.add( obs2 )

# define the source(s)
source = Source()
source.add(SourcePosition(Vector3d(10, 0, 0) * Mpc))
source.add(SourceIsotropicEmission())
source.add(SourceParticleType(nucleusId(1, 1)))
source.add(SourcePowerLawSpectrum(1 * EeV, 200 * EeV, -2.3))
source.add(SourceRedshiftEvolution(1.5, 0.001, 3))

# run simulation
sim.setShowProgress(True)
sim.run(source, 10000)
#output.close()

Although my photon output is empty, the output for hadrons comes nicely. What am I doing wrong in this case?

[JulienDoerner]

hey @k4zuyoh1,

If I see it correct, your simulation does not contain any module that can produce gamma-rays. You should either add them in your source (if there is a gamma-ray emitter) or add the secondaries from the PPP module or the PD. But for both, the propagation length of 10 Mpc may be to short to produce a noticeable amount of photons.

[k4zuyoh1]

Hey @JulienDoerner ,

Thanks for you answer! I've given the PD and PPP modules the True/False booleans and now my gamma output is filled. But now I have another issue: if I try to propagate the photons of this output with DINT/EleCa, I receive the following error, even though I've enabled additional columns:

RuntimeError: DintElecaPropagation: Wrong header of input file. Use PhotonOutput1D or Event1D with additional columns enabled.

My output has D, z, ID, E and x columns, and even commenting my z and x columns, DINT/EleCa doesn't work.

My photons observer is as follows:

obs2 = Observer()
obs2.add(ObserverSurface( Sphere(Vector3d(0.), 0.5 * Mpc))) 
obs2.add( ObserverDetectAll() ) 
obs2.add( ObserverElectronVeto() )
obs2.add( ObserverNeutrinoVeto() )
obs2.add( ObserverNucleusVeto() ) 
obs2.add(ObserverRedshiftWindow(-0.1, 0.1))
out2 = TextOutput('output-gamma.txt', Output.Event1D)
out2.enable(Output.CreatedIdColumn) 
out2.enable(Output.CreatedEnergyColumn)
out2.enable(Output.CreatedPositionColumn)
out2.enable(Output.RedshiftColumn)
obs2.onDetection( out2 )
sim.add( obs2 )

and my propagation modules are:

sim = ModuleList()
sim.add(PropagationCK(Bfield))
sim.add(Redshift())
sim.add(PhotoPionProduction(CMB(), True, False, True, True))
sim.add(PhotoPionProduction(IRB_Kneiske04(), True, False, True, True))
sim.add(PhotoDisintegration(CMB(), havePhotons=True))
sim.add(PhotoDisintegration(IRB_Kneiske04(), havePhotons=True))
sim.add(ElectronPairProduction(CMB()))
sim.add(ElectronPairProduction(IRB_Kneiske04()))
sim.add( EMPairProduction(CMB(), True) )
sim.add( EMPairProduction(IRB_Kneiske04(), True) )
sim.add( EMDoublePairProduction(CMB(), True) )
sim.add( EMDoublePairProduction(IRB_Kneiske04(), True) )
sim.add( EMInverseComptonScattering(IRB_Kneiske04(), True) )
sim.add( EMInverseComptonScattering(CMB(), True) )
sim.add( EMTripletPairProduction(CMB(), True) )
sim.add( EMTripletPairProduction(IRB_Kneiske04(), True) )
sim.add(NuclearDecay(True, True, False))
sim.add(MinimumEnergy(1 * EeV))
sim.add(MinimumRedshift(-0.1))

[JulienDoerner]

DINT and EleCa are deprecated since version 3.2. (see changelog) The interactions of the photons with the background is already included in the EM* modules.

[lukasmerten]

I would guess that something in the header of the TextOutput file was changed. However, as @JulienDoerner said we are no longer providing any support for DINT and EleCa. We suggest to use the above mentioned modules together with the thinning option.

If there is problem with these (EM*)-modules please open a new issue.

[lukasmerten]

Just for reference: The example in doc/pages/example_notebooks/secondaries/photons.v4.ipynb works with DINT and EleCa if you disable the CandidateTagColumn. It was added lately and is in string and not float format.

out = Output()
out.disable(Output.CandidateTagColumn)