Extended source and astromodel : To normalized or not to normalized ?

[GallegoSav]

I see from the line here

get_integrated_extended_model is normalizing the spatial shape to sr^-1 but astromodel is actually already expecting the spatial model to have units of deg⁻2 , so to be already normalized. When declaring the model it will otherwise raise a warning like :

16:58:39 WARNING   2D template read from                                                        [functions_2D.py](file:///uni-mainz.de/homes/sgallego/software/COSItools/cosipy_israel/Pyenv_israelcosipy/lib/python3.11/site-packages/astromodels/functions/functions_2D.py):[749](file:///uni-mainz.de/homes/sgallego/software/COSItools/cosipy_israel/Pyenv_israelcosipy/lib/python3.11/site-packages/astromodels/functions/functions_2D.py#749)
                  /localscratch/sgallego/linkToXauron/MC/Response/BalloonResponse/1_809MeV/COSI                    
                  _DC2_26Al_R5000_z1000_M60_3deg_normalized.fits is normalized to                                  
                  0.0027415567780803775 (expected: 1)

For my case I am using a SpatialTemplate_2D for the spatial shape, using a fits file and I take care to normalized the image so that sum(pixel)*domega = 1 (what astromodel is expecting) .

#normalise the spatial distribution of the fits file to deg²⁻1 (asked by astromodels)
hdul = fits.open('/localscratch/sgallego/linkToXauron/MC/Response/BalloonResponse/1_809MeV/COSI_DC2_26Al_R5000_z1000_M60_3deg.fits', mode='update')
WCS = wcs.WCS(hdul[0].header)
area = wcs.utils.proj_plane_pixel_area(WCS)
dOmega = (area * u.deg * u.deg).to(u.sr).value

#normalize
hdul[0].data = hdul[0].data/(np.sum(hdul[0].data)* dOmega) 

This is removing the warning from astromodel when loading the model actually. And I think the spatial shape using predefined astromodel function like Gaussian_on_sphere give already the correct units of deg^-2.

All of this to say that I think we should remove the / u.sr in normalized_map = extendedmodel.spatial_shape(coords.l.deg, coords.b.deg) / u.sr.

What do you think @israelmcmc and @hiyoneda ?

[israelmcmc]

Yes, I think you are right @GallegoSav. That line shouldn't have the /us.r since, as described in the astromodels documention, the output from evaluation an extended source already has units of "per solid angle area", and they decided to use 1/deg2 (choosing between 1/sr and 1/deg2 is arbitrary).

Can you please verify that, once you remove the /u.sr, the final output flux_map has the expected units of 1/(cm2 s deg2)?

[GallegoSav]

@israelmcmc the unit method return nothing for every axes of the injected model

model_injected = injector.inject_model(model = model, make_spectrum_plot = False, make_psichi_plot = True ,data_save_path = None)
model_injected.project('PsiChi').unit

[israelmcmc]

Seems like a problem. The PsiChi axis shouldn't have units though, since it's a HealpixAxis. But the energy and Phi axes should.

It seems that the ExtendedSourceResponse checks the axes labels but not the units, which is likely the source of the issue (+some user error, I suppose):

cosipy/cosipy/response/ExtendedSourceResponse.py

Line 37 in 11ac4f5

if not list(self.axes.labels) == ['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi']:

[GallegoSav]

Actually model_injected.axes["Em"].unit and model_injected.axes["PsiChi"].unitgives the correct unit but notmodel_injected.project('Emi').unit`

[israelmcmc]

That's seems fine then. I think model_injected is already a histogram of the expected counts (unitless).

[hiyoneda]

@GallegoSav I think that it is OK. / u.sr is not for normalization, just for adding a unit. I believe that spatial shapes generated using the predefined astromodel function have a unit of sr^-1, but extendedmodel.spatial_shape returns a pure NumPy array without a unit. To make a generated histogram with a correct unit, I added it there.
I also paste my test code here.

from astromodels import *
import numpy as np
import healpy as hp
from histpy import Histogram, Axes, Axis, HealpixAxis

spectrum = Gaussian(F = 1e-3, mu = 1809, sigma = 1.8) # spectral model
emission_map = Disk_on_sphere(lon0 = 15, lat0 = -5, radius = 10) # spatial model

diffuse_model = ExtendedSource('disk_gaussian',spectral_shape=spectrum,spatial_shape=emission_map) # model

# axes for the function
energy_axis = Axis(np.array([1800.0, 1820.0]) * u.keV)
image_axis = HealpixAxis(nside = 32, coordsys='Galactic')

# test
npix = image_axis.npix
coords = image_axis.pix2skycoord(np.arange(npix))
np.sum(map_with_spectrum.spatial_shape(coords.l.deg, coords.b.deg)) #-> 982.2527326975237
np.sum(map_with_spectrum.spatial_shape(coords.l.deg, coords.b.deg)) * image_axis.pixarea().value #-> 1.0045 -> it means that map_with_spectrum.spatial_shape(coords.l.deg, coords.b.deg) has an unit of sr^-1, but ...
diffuse_model.spatial_shape(coords.l.deg, coords.b.deg) #-> yields a pure numpy array, so it is unitless 

[israelmcmc]

Interesting. Despite the astromodels documentation saying that the output is in 1/deg2, the test code does seem to indicate that it's in 1/sr 🤔

[GallegoSav]

@israelmcmc Should I close the issue then ?

[israelmcmc]

@GallegoSav Do you understand why you were getting the warning though? It seems that it's not because these units, but there seems to be still some normalization issue somewhere

[GallegoSav]

@israelmcmc Yes, when you provide your own spatial shape with the 2Dtemplate , Threeml is making this check :

 # test if the map is normalized as expected
            area = wcs.utils.proj_plane_pixel_area(self._wcs)
            dOmega = (area * u.deg * u.deg).to(u.sr).value
            total = self._map.sum() * dOmega

            if not np.isclose(total, 1, rtol=1e-2):
                log.warning(
                    "2D template read from {} is normalized to {} (expected: 1)".format(
                        self._fitsfile, total
                    )
                )

So that the units of the map you provide is well deg^-2 . I was having the warning because I was just normalizing the map to 1, not to self._map.sum() * dOmega

[hiyoneda]

OK, then it looks resolved, but please let me confirm what is done in the threeml and my understanding @GallegoSav @israelmcmc

Above the code, dOmega is converted to 'sr', and total is assumed to be normalized to 1. Then, should self._map.sum() (and self._map) have a unit of 'sr^-1'? Still, I am a bit confused about the unit of the output of the spatial models...

[GallegoSav]

@hiyoneda Right in the documentation it says deg^-2 but the check in the code clearly show they wants sr⁻1 so the unit of self._map is sr^-1

[israelmcmc]

I went through the rabbit hole, and so far I think it's an issue in the astromodels side. I opened this issue: threeML/astromodels#238

I'll leave this open while the astromodels developers answer and to see what we should to.