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
71 stars 69 forks source link

Tracking UHECRs corresponding to each secondary produced #333

Closed saikatxdas closed 2 years ago

saikatxdas commented 3 years ago

Hi,

I have a usage question. I am trying to track which EM secondaries are produced from which primary UHECR using ParticleCollector() but I end up getting SN0 column to be 448 for all the rows. I am using something like this.

# Observer for UHECRs
obs = Observer()
obs.add(ObserverPoint())
obs.add(ObserverElectronVeto())
output = TextOutput('cr.txt', Output.Event1D)
obs.onDetection(output)
sim.add(obs)

# Observer for EM particles
obs2 = Observer()
obs2.add(ObserverDetectAll())
obs2.add(ObserverNucleusVeto())
out2 = ParticleCollector()
obs2.onDetection(out2)
sim.add(obs2)

ObserverDetectAll() should detect the particles as soon as they are produced. But how come all of them have the same SN0 in out2?

lukasmerten commented 3 years ago

Could you please add a minimal example that is able to reproduce that behavior. With the provided information we cannot provide any help.

saikatxdas commented 3 years ago

Thank you for your reply. Yes, I have attached a working script emsim.txt

Here, the SN0 values generated seem erroneous. Because particles with the same SN0 must have the same E0. Here, they differ.

lukasmerten commented 3 years ago

I finally found some time to look into this.

Have you checked what happens if you use a TextOutput for your electrons, too? # Observer for EM particles out2 = TextOutput('em_part.txt', Output.Event1D) out2.enable(Output.SerialNumberColumn) obs2 = Observer() obs2.add(ObserverDetectAll()) obs2.add(ObserverNucleusVeto()) obs2.onDetection(out2)

In this case, I get around 8500 unique parent serial numbers (SN0 or SN1) for all observed electrons, which seems reasonable to me.

Do not forget to close both outputs after the simulation, when using jupyter notebooks or similar tools for your simulation.

lukasmerten commented 3 years ago

@saikat8018 could you try to explicitly set the constructor arguments of the ParticleCollector

ParticleCollector(int(1e6), False, False)

and keep the rest of your script the same? For me it is working in this case.

Which brings me to the question: Why does it make a difference if I set the default arguments explicitly? Does someone have an idea?

Here is also a smaller script that is able to reproduce this issue.

@TobiasWinchen and @adundovi is this related to #244?

lukasmerten commented 2 years ago

@saikatxdas can you confirm that my work-around by explicitely setting the arguments, does what you want?

saikatxdas commented 2 years ago

@lukasmerten Sorry for not following it up earlier. Yes, I have checked it now on CRPropa 3.1.7 and it works. Thank you for your help. This indeed works