CRPropa / CRPropa3

CRPropa is a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles. https://crpropa.github.io/CRPropa3/
https://crpropa.desy.de
GNU General Public License v3.0
72 stars 69 forks source link

Associating Secondary Photons with Parent Protons in CRPropa Simulations #508

Closed Pier-mic closed 1 week ago

Pier-mic commented 2 weeks ago

Hello everyone,

I'm currently working on a simulation using CRPropa where I propagate primary protons and generate secondary photons. I would like to associate each produced photon with the specific proton that generated it. Essentially, I want to track the parent-progeny relationship between the primary protons and the secondary photons in my simulation data.

I've attempted to use the EventGroupID to achieve this association, but I'm not sure if I'm implementing it correctly. Below is a simplified version of my code:

from crpropa import *
import numpy as np

# Output file names
filename_primary = 'primary_protons.txt'
filename_photons = 'secondary_photons.txt'

# Source configuration
source = Source()
source.add(SourceParticleType(nucleusId(1, 1)))  # Protons
source.add(SourcePowerLawSpectrum(1e19 * eV, 1e20 * eV, -2))
source.add(SourceUniform1D(0, 100 * Mpc))

# Module list
modules = ModuleList()
modules.add(SimplePropagation(0, 10 * Mpc))
modules.add(MinimumEnergy(1e18 * eV))

# Observers for protons and photons
observer_protons = Observer()
observer_protons.add(Observer1D())
observer_protons.add(ObserverPhotonVeto())
observer_protons.add(ObserverElectronVeto())
output_protons = TextOutput(filename_primary, Output.Event1D)
output_protons.setEnergyScale(eV)
output_protons.enable(Output.WeightColumn)
output_protons.enable(Output.EventGroupIDColumn)  # Enable EventGroupID
observer_protons.onDetection(output_protons)

observer_photons = Observer()
observer_photons.add(Observer1D())
observer_photons.add(ObserverElectronVeto())
observer_photons.add(ObserverNucleusVeto())
output_photons = TextOutput(filename_photons, Output.Event1D)
output_photons.setEnergyScale(eV)
output_photons.enable(Output.WeightColumn)
output_photons.enable(Output.EventGroupIDColumn)  # Enable EventGroupID
observer_photons.onDetection(output_photons)

modules.add(observer_protons)
modules.add(observer_photons)

# Physical processes
cmb = CMB()
modules.add(PhotoPionProduction(cmb, True, False, False))

# Run the simulation
modules.run(source, 100000, True)

# Load data
data_protons = np.loadtxt(filename_primary, dtype=np.float64)
data_photons = np.loadtxt(filename_photons, dtype=np.float64)

# Attempt to associate photons with parent protons
proton_dict = {int(row[4]): row for row in data_protons}  # EventGroupID as key
for photon in data_photons:
    event_group_id = int(photon[4])
    parent_proton = proton_dict.get(event_group_id)
    if parent_proton is not None:
        # Process the photon and its parent proton
        pass

Despite this approach, I'm not confident that the EventGroupID is correctly linking the photons to their parent protons. Is there a recommended method or best practice in CRPropa for achieving this association? Any guidance or examples would be greatly appreciated.

Thank you for your help!

Best regards,

mohller commented 2 weeks ago

Hi, what is the intended goal? For the case you propose (1D) the result would be a distribution of the positions of production for the photons, and it might be simpler to obtain the distribution of photons as the exponential attenuation of protons with rate as provided by CRPropa. If the goal is to see the photons from the successive interaction of protons as they cascade in energy, this might achievable with this code, but you will need to correlate the positions of creation for photons and the positions in the protons output file. This is because in CRPropa at a given step size can produce many interactions and thus you can have multiple photons coming from one step for one proton.

One final note, in this example both photons and protons are only written to the output file when they reach the position Vector3D(0, 0, 0), and I would imagine this does not work for the stated goal.

Pier-mic commented 2 weeks ago

My goal is to produce an energy spectrum of secondary photons, generated by primary protons with an energy spectrum proportional to E^-3. To do this, I first need to produce an energy spectrum of primary protons proportional to E^-1, and then adjust the weights of the particles in the simulation to make it an E^-3 spectrum. For the protons, I managed to do this fairly easily, but the issue is that I need to modify the spectrum of secondary photons to ensure it is produced by primary protons with a spectrum proportional to E^-3. To achieve this, I need to know exactly which proton generated a given photon, and I don’t know how to do that.

Inviato da Outlook per iOShttps://aka.ms/o0ukef


Da: Leonel Morejon @.> Inviato: Thursday, November 7, 2024 12:36:07 PM A: CRPropa/CRPropa3 @.> Cc: Pier-mic @.>; Author @.> Oggetto: Re: [CRPropa/CRPropa3] Associating Secondary Photons with Parent Protons in CRPropa Simulations (Issue #508)

Hi, what is the intended goal? For the case you propose (1D) the result would be a distribution of the positions of production for the photons, and it might be simpler to obtain the distribution of photons as the exponential attenuation of protons with rate as provided by CRPropa. If the goal is to see the photons from the successive interaction of protons as they cascade in energy, this might achievable with this code, but you will need to correlate the positions of creation for photons and the positions in the protons output file. This is because in CRPropa at a given step size can produce many interactions and thus you can have multiple photons coming from one step for one proton.

One final note, in this example both photons and protons are only written to the output file when they reach the position Vector3D(0, 0, 0), and I would imagine this does not work for the stated goal.

— Reply to this email directly, view it on GitHubhttps://github.com/CRPropa/CRPropa3/issues/508#issuecomment-2462009429, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BMJAGC2XKWN3FPRSB6CVDBDZ7NGCPAVCNFSM6AAAAABRK52RDSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINRSGAYDSNBSHE. You are receiving this because you authored the thread.Message ID: @.***>

mohller commented 2 weeks ago

In this case I would recommend sampling directly from SOPHIA, the generator of photohadronic interactions that CRPropa employs.

Pier-mic commented 2 weeks ago

So, should I access the source code of SOPHIA and see how its libraries work?

Inviato da Outlook per iOShttps://aka.ms/o0ukef


Da: Leonel Morejon @.> Inviato: Thursday, November 7, 2024 4:12:20 PM A: CRPropa/CRPropa3 @.> Cc: Pier-mic @.>; Author @.> Oggetto: Re: [CRPropa/CRPropa3] Associating Secondary Photons with Parent Protons in CRPropa Simulations (Issue #508)

In this case I would recommend sampling directly from SOPHIA, the generator of photohadronic interactions that CRPropa employs.

— Reply to this email directly, view it on GitHubhttps://github.com/CRPropa/CRPropa3/issues/508#issuecomment-2462485676, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BMJAGC2UYWZ7ZZURM4E6E2TZ7N7NJAVCNFSM6AAAAABRK52RDSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINRSGQ4DKNRXGY. You are receiving this because you authored the thread.Message ID: @.***>

mohller commented 1 week ago

You can download it and use it as explained in its website https://www.uibk.ac.at/projects/he-cosmic-sources/tools/sophia/index.html.en