Remove PhotonOutput1D, use regular output instead

[avvliet]

It might be useful at some point to allow for HDF5 output in the PhotonOutput1D module.

[TobiasWinchen]

Where is the difference to the regular Observer with the feature ObserverNucleusVeto with HDF5 Output? If there is none, we probably also should remove the PhotonOutputut1D as this is redundant to the text output.

[avvliet]

I agree it might be better to remove the PhotonOutput1D module and only use the regular output, also for photons, but at the moment that does not seem to work. If you run the file pasted below the photon output from 'obs1' will be empty, while the 'PhotonOutput1D' will give you a very long list of mainly electrons and positrons.

from crpropa import *

Run = 'PhotonOutputTestNucleusVeto'

neutrinos = False
photons = True

# module setup
EBL = IRB_Stecker16_upper
m = ModuleList()
m.add(SimplePropagation(1*kpc, 10*Mpc))
m.add(Redshift())
m.add(PhotoPionProduction(CMB, photons, neutrinos))
m.add(PhotoPionProduction(EBL, photons, neutrinos))
m.add(NuclearDecay(photons, photons, neutrinos))
m.add(ElectronPairProduction(CMB, photons))
m.add(ElectronPairProduction(EBL, photons))
m.add(MinimumEnergy(10.**8 * eV))
m.add(PhotonOutput1D('photon-input'+Run+'.txt'))

# observers
obs1 = Observer()
obs1.add(ObserverPoint())
obs1.add(ObserverNucleusVeto())  # we don't want nuclei here
output1 = TextOutput('photons-'+Run+'.txt', Output.Everything)
obs1.onDetection( output1 )
m.add(obs1)
obs2 = Observer()
obs2.add(ObserverPoint())
obs2.add(ObserverPhotonVeto()) # we don't want photons here
obs2.add(ObserverElectronVeto())
output2 = TextOutput('nucleons-'+Run+'.txt', Output.Everything)
obs2.onDetection( output2 )
m.add(obs2)

# source
source = Source()
source.add(SourceUniform1D(0, redshift2ComovingDistance(1)))
source.add(SourceRedshift1D())
source.add(SourcePowerLawSpectrum(1 * EeV, 100 * EeV, -1.0))
source.add(SourceParticleType(nucleusId(1, 1)))

# run simulation for 10 primaries and propagate secondary photons
m.setShowProgress(True)
m.run(source, 10, True)

[TobiasWinchen]

You cannot operate both, as the PhotonOutput deactivates the Photons. If I remove the PhotonOutput1D from your example I get some output.

[avvliet]

Ah yes, you're right. Sorry about that one. We should be a bit careful though. If you want to enter the photon output in, for example, DINT, you need the distance to the point of creation of the photon, and the energy of the photon at that point. I don't think it is possible to get that in the current standard output.

[TobiasWinchen]

Instead of using a ObserverPoint you can use ObserverDetectAll for the photons

• that will on the first step write out the particle. Then you get the creation point for further usage in e.g. DINT.

[avvliet]

Good point, yes that should work. But in that case you get the trajectory length (distance from source to interaction in 1D) not the distance from the creation to the observer as you would need for DINT. If you output everything, though, you can use 'X0' to compute this distance to the observer.

So yes, I completely agree, we should get rid of 'PhotonOutput1D'. The implementation of the regular output into DINT might not be so convenient as with PhotonOutput1D but it doesn't make sense to have these multiple ways of producing the output.

[TobiasWinchen]

Wouldn't X be the quantity you need?

[avvliet]

Ah yes, of course.

[TobiasWinchen]

I just added deprecation warnings. Could you update the examples on the wiki accordingly? We can then remove this shortly before releasing 3.2.

[saikatxdas]

When I use a regular observer with ObserverNucleusVeto() and pass the output to DINT, I get the error

DintPropagation: Energy too low -inf

But this doesn't happens when I use PhotonOutput1D(). I am using the script given in examples for secondary photons from proton propagation and have turned on the electrons in ElectronPairProduction. I have also included interactions with IRB as well. Here is the minimal script to reproduce the error:

from crpropa import *

# source setup
source = Source()
source.add(SourceParticleType(nucleusId(1, 1)))
source.add(SourcePowerLawSpectrum(10 * EeV, 100 * EeV, -2))
source.add(SourceUniform1D(3 * Mpc, 100.00001 * Mpc))

# setup module list for proton propagation
m = ModuleList()
m.add(SimplePropagation(0, 10 * Mpc))
m.add(MinimumEnergy(1 * EeV))

# observer
obs1 = Observer() # proton output
obs1.add( ObserverPoint() )
obs1.add( ObserverPhotonVeto() ) # we don't want photons here
obs1.onDetection( TextOutput('proton_output.txt', Output.Event1D) )
m.add(obs1)

obs2 = Observer() # photon output
obs2.add( ObserverDetectAll() ) # stores the photons at creation without propagating them
obs2.add( ObserverNucleusVeto() ) # we don't want hadrons here
out2 = TextOutput('photon_output.txt', Output.Event1D)

out2.disableAll()
out2.enable(Output.CurrentIdColumn)
out2.enable(Output.CurrentEnergyColumn)
out2.enable(Output.CurrentPositionColumn)
out2.enable(Output.CreatedIdColumn)
out2.enable(Output.CreatedEnergyColumn)
out2.enable(Output.SourceIdColumn)
out2.enable(Output.SourceEnergyColumn)
obs2.onDetection( out2 )
m.add(obs2)

# enable secondary electrons
m.add(ElectronPairProduction(CMB, True))
m.add(ElectronPairProduction(IRB, True))

# enable secondary photons
m.add(PhotoPionProduction(CMB, True))
m.add(PhotoPionProduction(IRB, True))

# run simulation
m.setShowProgress(True)
m.run(source, 100000, True)

crpropa.DintPropagation("photon_output.txt", "dint_out.dat", 4, 4, 1*crpropa.nG)

[rafaelab]

Hi Saikat, The output is in the wrong format. Check the required format here.

[saikatxdas]

Dear Rafael,

I have an older CRPropa version (October 2018) installed that doesn't have all the latest commits. However, I have made sure that the output format is the same for PhotonOutput1D and that with a regular observer. I see that the output of PhotonOutput1D is in exponential format, while that of the regular observer is in decimal format. Thus, a particle of energy 1e-7 EeV will be saved as 0.000000. How can I change this?

[rafaelab]

You can either download the latest version or change the source code of the functions that write the output (see https://github.com/CRPropa/CRPropa3/blob/master/src/module/TextOutput.cpp). Don't forget to recompile CRPropa after you implement your changes.

[saikatxdas]

Hi Rafael,

I have done a fresh installation of the latest version of CRPropa3. make test runs successfully without any failure. But when I try to propagate the photons using DINT, the program crashes with segmentation fault (core dumped) at the line crpropa.DintPropagation(...). When I try to debug using (gdb), I get the following:

(gdb) file python
Reading symbols from python...(no debugging symbols found)...done.
(gdb) run xphoton.py 
Starting program: /home/saikat/.virtualenvs/crpropa/bin/python xphoton.py
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
[New Thread 0x7ffff18bb700 (LWP 8524)]
[New Thread 0x7fffef0ba700 (LWP 8525)]
[New Thread 0x7fffec8b9700 (LWP 8526)]
crpropa::ModuleList: Number of Threads: 4
Run ModuleList
[New Thread 0x7fffe840f700 (LWP 8527)]
[New Thread 0x7fffe7c0e700 (LWP 8528)]
[New Thread 0x7fffe740d700 (LWP 8529)]
  Started Fri Apr 19 15:32:50 2019 : [ Finished ] 100%    Needed: 00:00:17  - Finished at Fri Apr 19 15:33:07 2019

Thread 1 "python" received signal SIGSEGV, Segmentation fault.
0x00007ffff618e9c3 in crpropa::DintPropagation(std::__cxx11::basic_string<char, 
std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char,
std::char_traits<char>, std::allocator<char> > const&, int, int, double, double) ()
   from /home/saikat/.virtualenvs/crpropa/lib/libcrpropa.so

I do not understand the problem. Can you please help? This is the script of the simulation that I am executing. Should I open a separate issue for this?

[rafaelab]

Hi, I understand that the PhotonOutput1D issue has been solved. The other issue reported doesn't belong to this thread and is not a CRPropa bug, so I will close the issue.