rafaelab
commented
3 years ago Hi,
There is no reason to expect the spectral shape to remain the same if interactions are included. There will be energy losses, production of secondary particles, etc. Your figure shows exactly that (the GZK cut-off, in this case).
The memory problem you are encountering is because many secondaries are produced by interaction, quickly exceeding the memory of most machines. You can, for instance, deactivate photons and electrons to solve this (at least in part). Or you can run many simulations with fewer events and combine the results afterwards.
Hi, I am trying to plot cosmic rays at the source and at the observer for a 1D simulation using the following code and I am also attaching the result I am getting. I have chosen a power law -2 at the source. Without any interaction module added, the source and observer spectra are exactly overlapping but when I add the interaction modules the source spectrum is completely modified. Shouldn't the source spectrum remain the same? Or am I doing any mistake?
Also, what is the maximum number of particles CRPropa can simulate? In my case, I am able to simulate successfully about 1e8 particles but beyond this, the program automatically gets killed. I am using a system with 8GB RAM and 1.8 GHz 8 core processor.
from crpropa import *
photons=True neutrinos=True electrons =False
sim = ModuleList() sim.add( SimplePropagation() ) sim.add( Redshift() ) sim.add(PhotoPionProduction(CMB(), photons, neutrinos, electrons)) sim.add( ElectronPairProduction(CMB(), True) ) sim.add( MinimumEnergy( 0.1 * EeV) )
observer and output
obs = Observer() obs.add( ObserverPoint() ) output = TextOutput('testProton.txt', Output.Event1D) output.disableAll() output.enable(Output.SourceIdColumn)
output.enable(Output.SourceEnergyColumn)
output.enable(Output.CurrentIdColumn)
output.enable(Output.CurrentEnergyColumn) obs.onDetection( output ) sim.add( obs )
nucleusId(1, 1)
source = Source() source.add(SourceParticleType(nucleusId(1, 1))) source.add(SourcePowerLawSpectrum(1 EeV, 10000 EeV, -2.0)) source.add(SourcePosition(Vector3d(50,0,0) *Mpc)) source.add( SourceRedshift1D() )
run simulation
sim.setShowProgress(True) sim.run(source, 100000, True)
from pylab import *
import pylab as t
output.close() figure(figsize=(6,6)) a = loadtxt("testProton.txt") E = logspace(17,22,20) idx0 = a[:,2] == 1000010010 sprotons = a[idx0,3]1e18 idx = a[:,0] == 1000010010 protons = a[idx,1]1e18
bincenter = (E[1:] -E[:-1])/2 + E[:-1] binwidth = E[1:] -E[:-1]
data,bins = histogram(sprotons,E) J = data/binwidth plot(bincenter, J * bincenter**2,label="source Protons, 1e5")
data1,bins = histogram(protons,E) J1 = data1/binwidth plot(bincenter, J1 * bincenter**2,label="Protons propagated")
grid()
loglog() xlim(3e18, 3e21) ylim(1e20, 1e26) legend(loc="lower right") xlabel("Energy [eV]") ylabel ('$E^2 dN/dE$ [a.u.]') show()