Error reading CR 3D distribution

[ottaviofornieri]

Hi!

When trying to read the regular 3D CR distribution from DRAGON that comes with the HERMES release, I get this error: AttributeError: module 'pyhermes.cosmicrays' has no attribute 'Dragon3D'.

Thanks a lot in advance!

[adundovi]

Hi Ottavio,

thank for your report. It should be in the version I've pushed a week ago ( c98bc53dfe41c14bed12418c78542214d1d48c66 ). These are possible constructors for it:

1. pyhermes.cosmicrays.Dragon3D(arg0: pyhermes.PID)
2. pyhermes.cosmicrays.Dragon3D(arg0: List[pyhermes.PID])
3. pyhermes.cosmicrays.Dragon3D(arg0: str, arg1: List[pyhermes.PID])

Try installing the newest version in your old build directory:

git pull
make -j
make install

Tell me if it worked.

[ottaviofornieri]

Ciao Andrej! yes it worked thanks, I was sure I updated to the latest version already. However, with the command density = cr.getDensityPerEnergy(1, sun_pos) to read the particle distribution, I get:

TypeError: getDensityPerEnergy(): incompatible function arguments. The following argument types are supported:
    1. (self: pyhermes.cosmicrays.Dragon3D, arg0: pyhermes.QEnergy, arg1: pyhermes.Vector3QLength) -> pyhermes.QPDensityPerEnergy

Invoked with: <pyhermes.cosmicrays.Dragon3D object at 0x7f675d0ca2b0>, 1, <pyhermes.Vector3QLength object at 0x7f6721f010f0>

In principle, I can go ahead, but then I get negative values after the integration, so something is not read properly.

[carmeloevoli]

I think that 1 is not recognised since has no units. The first argument must be arg0: pyhermes.QEnergy Try with 1_GeV (or the equivalent in python).

[ottaviofornieri]

Hi Carmelo, thanks for the suggestion, but still cannot read it, although something has changed. I tried with the dimensional energy 1*GeV, and I get a different error.

IndexError                                Traceback (most recent call last)
<ipython-input-27-22bd50d15fbc> in <module>
     16 cr = cosmicrays.Dragon3D(fitsfile, [Proton])
     17 cr_list = [cr]
---> 18 density = cr.getDensityPerEnergy(1.0*GeV, sun_pos)
     19 #print(density)
     20 

IndexError: map::at

[adundovi]

Hi @ottaviofornieri, sorry for waiting. Since a Dragon3D file has well-defined energy bins, you cannot choose arbitrary energy as an argument for it. Or, to phrase it differently - we don't have a technique to easily interpolate CR density from a 4D grid (energy, pos_x, pos_y, pos_z). The code only interpolates CR density from the spatial part of the grid: (pos_x, pos_y, pos_z). That's why we have an energy iterator that "knows" which energy bins are present and possible as the first argument in getDensityPerEnergy(). For analytical models where you calculate the density for any given energy and position, this would not be the case, but for the grid data, it has to be like this (or try to interpolate a 4D grid, however, I don't see a use case with it). We should probably explain this in the documentation.

The following code will work:

from pyhermes import *

cr = cosmicrays.Dragon3D([Proton])

for e in cr.getEnergyAxis():
    print("Energy: ", e)
    density = cr.getDensityPerEnergy(e, Vector3QLength(8*units.kpc, 0*units.kpc, 0*units.kpc))
    print(density)

EDIT: ...and, as @carmeloevoli wrote, nothing works without appropriate units :-)

[ottaviofornieri]

Ciao @adundovi! thanks now this method works and I can ready the densities from my FITS file.

Thanks a lot!