Revised units for parameter K (ph/cm2/s) in extended source fit tutorial to match the units in the SpecFromDat file

[krishnatejavedula]

Addresses issue #194

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 80.46%. Comparing base (02c7bc6) to head (9b386be).
Report is 39 commits behind head on develop.

Files with missing lines	Coverage Δ
cosipy/threeml/custom_functions.py	40.59% <ø> (ø)

... and 9 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[israelmcmc]

Hi @krishnatejavedula, thanks for working on this. First, I verified that you are changing only this line

-    "K          = [0.00046, 0.0011, 0.0027, 4.5e-3]/u.cm/u.cm/u.s/u.keV\n",
+    "K          = [0.00046, 0.0011, 0.0027, 4.5e-3]/u.cm/u.cm/u.s\n",

the rest are just images in the notebook that are actually the same, as expected.

I think there's something still off. Looking at SpecFromDat the K should indeed be integral flux --i.e. 1/cm2/s, since we perform this normalization operation:

dataFlux = dataFlux  / np.sum(dataFlux * ewidths)

However, the evaluate() should still have units of flux density (per unit of energy) but this is not what happens. For example:

from cosipy.threeml.custom_functions import  SpecFromDat
import astropy.units as u

spec = SpecFromDat(dat="</path/to/cosipy>/docs/tutorials/spectral_fits/extended_source_fit/OPsSpectrum.dat")

K = 1/u.cm/u.cm/u.s

spec.evaluate([1]*u.MeV, K)

returns [7.35916986e-06] 1 / (cm2 s).

The issue is that when you divide by ewidths the units should also be included. As per the .dat documentation, ewidths has units of keV, the missing unit. In summary, I agree with your change in the notebook, but this should be accompanied by the appropriate change in SpecFromDat. Overall, although the result will be identical as it was before, the units and meaning of each variable will make sense.

I'm actually not sure how this will affect other code using SpecFromData. I think it's mostly the SourceInjector, so you should check. I tried to check the SourceInjector tutorial but it doesn't have units!:

spectrum = SpecFromDat(K=1/18, dat="crab_spec.dat")

I think that should be revisited as well.

[krishnatejavedula]

Hello @israelmcmc I have updated the SpecFromDat class to ensure correct unit consistency in the evaluate() method. I've checked the notebooks that use this function, and the change does not affect their functionality.

[israelmcmc]

Closing and reopening to re-run the tests.

[israelmcmc]

@krishnatejavedula I'm getting this error

UnitConversionError: '1 / (s cm2)' (particle flux) and '1 / (keV s cm2)' are not convertible

when executing the extended source notebook, here:

Cell In[34], line 6
      4 c = SkyCoord(l=0*u.deg, b=0*u.deg, frame='galactic')
      5 c_icrs = c.transform_to('icrs')
----> 6 ModelCentralPoint = PointSource('centralPoint', ra = c_icrs.ra.deg, dec = c_icrs.dec.deg, spectral_shape=SpecCentralPoint)
      7 ModelNarrowBulge = ExtendedSource('narrowBulge',spectral_shape=SpecNarrowBulge,spatial_shape=MapNarrowBulge)
      8 ModelBroadBulge = ExtendedSource('broadBulge',spectral_shape=SpecBroadBulge,spatial_shape=MapBroadBulge)

Can you please try to run it again after syncing with develop? Maybe it's something that changed since you opened this PR (sorry for just getting to it!)

[krishnatejavedula]

Hi @israelmcmc, it's okay. I am checking it now.

[krishnatejavedula]

Hello @israelmcmc.

I have reverted the changes I made in this PR and updated the docstring for SpecFromDat so that the units of K match those used in the notebooks.

We may need to revisit the units of K in the future, but for now, this resolves the issues we were encountering.