Improve interpolation method and edge case handling in SpecFromDat:

[krishnatejavedula]

• Use 'fill_value' in interp1d for better extrapolation.
• Remove normalization of dataFlux to use raw data.

[israelmcmc]

Thanks for working on this @krishnatejavedula.

About the normalization:
I think this is a good change. We shouldn't renormalize the user input based on the sum of the y-values. Since it happens behind the scenes, it's error-prone. In addition, it makes the result depend on the bin size, which is not very intuitive. Let's just leave the built-in value K to be the normalization in front of the values provided by the user.

However, this change will break some tutorials. I think only the extended source model fit tutorial. Can you make the appropriate change such that it returns the same result? It's probably only a matter of changing the values in the input file to account with a constant in front.

About the fill_value:
I think it should be kept at 0. Let’s leave it to the user how far the spectrum needs to be extended. If they need to extrapolate it with a constant (which is usually not the case), they can add extra points at the end.

[krishnatejavedula]

The above changes will remove the line dataFlux = dataFlux / sum(dataFlux) #normalized since the dat file has one point per keV from SpecFromDat. This change will affect the results in some notebooks (example: diffuse_511_spectral_fit.ipynb). To rectify this the data files used in the notebooks can be modified by using the following code:

# Load data excluding comments, skip footer and header
data = np.genfromtxt('OPsSpectrum.dat', comments="#", usecols=(0, 1, 2), skip_footer=1, skip_header=5)

# Extract columns
dataEn = data[:, 1]  # Energy column
dataFlux = data[:, 2]  # Flux column

# Calculate the normalization factor
normalization_factor = np.sum(dataFlux)

# Normalize dataFlux
dataFlux_normalized = dataFlux / normalization_factor

dataFlux / np.sum(dataFlux)
# Combine data into new format
updated_data = np.column_stack((data[:, 0], dataEn, dataFlux_normalized))

# Header and comments
header = "IP LIN "
comments = "# # Format: DP <energy in keV> <shape of O-Ps [XX/keV]> "

# Save the data to a new file
np.savetxt('adjusted_OPsSpectrum.dat', updated_data, fmt='%f', comments='# ' + comments, header=header)```

This will normalize the data and save it to a new file which can be used in the notebook to achieve similar results.

[ckarwin]

I'm a bit confused about the changes being made here. I am trying to understand why dataFlux was initially normalized, and what the motivation is for changing it. Ultimately, the evaluate function interpolates the spectrum, and scales it by the normalization constant K, which I believe is in line with the 3ML convention described here, and probably why it was originally formatted this way. Can you please clarify the motivation for changing this, and if things will still be compatible with the 3ML conventions?

[israelmcmc]

@ckarwin Since I suggested this change, let me explain it.

I think there are 3 sensible ways to define the normalization K:

1. The simplest, K is just a unitless value that multiplies all flux points from the file. It's up to the user to calculate the final physical flux cofrectly. That's what 3ML does in TemplateModel. As a side point, I think we should replace SpecFromData with TemplateModel once we figure out how to make it work with cosipy.
2. The second, which I think is in the functions in the link you shared, is for K to be the differential flux at a given reference energy. In this case, the user needs to provide the reference energy in addition to the data file, which now only defines the shape of the spectrum. You would need to interpolate at the reference energy and divide by that value.
3. Lastly, K could be the integrated flux within an energy range. The user needs to provide a definition of this range or assume that it corresponds to min/max energy values in the file.

The current version that this PR is replacing is close to option 3, but only when the points are sampled every 1 keV. The normalization was the sum of the y values, without taking into account the bin size. Since it is bin and unit-dependent, I suggested changing it to option 1, which is in line with astromodel's TemplateModel.

[ckarwin]

@israelmcmc Thanks for the clarification. It makes sense. I see now that the way this was implemented in example 3 of the extended source spectral fit tutorial for DC2 was not completely correct. I agree that the first approach you described is probably the most straight forward way. It would be helpful to specify the requirements for the file type and corresponding normalization factor K in the doc string of SpecFromDat.

I guess then that modifying the input spectral files for the DC2 tutorial will be a quick and easy patch to reproduce the originally published results. But I think the files still need to be changed directly, and there is no need to supply the code for users to do it themselves. We can revisit the details of the extended source fit later on (and in future examples).

I also note that some of the automatic tests are failing now, but it seems related to the issue that is fixed with PR #192.

[israelmcmc]

Thanks, @ckarwin

It would be helpful to specify the requirements for the file type and corresponding normalization factor K in the doc string of SpecFromDat.

Yes, I agree. @krishnatejavedula can you add this to the docstring, please? The file energy units should be keV, and flux units 1/cm2/s/kev. The normalization K is unitless.

We can revisit the details of the extended source fit later on (and in future examples).

Indeed. I opened an issue about it: #194. I think it's a "good first issue" and labeled it as such.

I also note that some of the automatic tests are failing now, but it seems related to the issue that is fixed with PR #192.

@krishnatejavedula I just merged #192. Can you sync this PR with the develop branch in order to fix the unit tests, please?

[israelmcmc]

The latest of @krishnatejavedula commit (which I asked him to) made me realize that I screwed this up. I hadn't realize that the .dat file was using the pre-defined MEGAlib format. From the cosima manual:

That is, the normalization constant K has units of total particle flux (1/cm2/s), integrated over the range of the spectrum. The units of the values of the spectrum should be proportional to ~1/keV (that is, accounting for the bins size), but the overall normalization depends on K.

Andreas had reasons to do it this way for simulations, and although I'm not sure it's the best definition for the analysis, I think we should stick to the previously-defined file format definition, otherwise, it will be really error-prone.

Summarizing, the spectrum points, in 1/cm2/s/keV should be:

$$ F_i = K \frac{y_i \Delta E_i}{\sum_i y_i \Delta E_i} $$

Since all energy bins in the current file are equal, the current code works. However, it should be fixed to work in the general case. Also, we should make sure the docstring matches the cosima documentation.

@ckarwin Can you please double-check that I got it right this time, since you've used this spectrum file format before?
@krishnatejavedula Does this make sense? Can you please make the appropriate changes? You will also need to adjust your input to the SourceInjector so it matches the analytic function @Yong2Sheng was using. I'm sorry for the extra work due to my mistake.

[ckarwin]

@israelmcmc Ok yes, if we are going to stick with the same convention as MEGAlib, then that should help to make things less error prone. Everything you said makes sense for the most part, but I don't think the differential flux ($F_i$) should have $\Delta E_i$ in the numerator. Please see my reasoning below, which also outlines what I think the function needs to do:

The spectral file has the form:

E_1 y_1
E_2 y_2
.    .
.    .
.    .

where

$[y_i] \propto 1/\mathrm{keV}$.

The normalization of the spectrum is given by $K$, which is the integral of the differential flux over the energy range defined in the file, with

$[K] = \mathrm{\frac{ph}{cm^2s}}$.

So the function should take the input spectral values from the file and normalize them by the integral:
$y_i \rightarrow y_i' = \frac{y_i}{y_\mathrm{int}}$,

where

$y_\mathrm{int}=\sum_i y_i \Delta E_i$.

Then the differential flux is given by

$F_i = Ky_i' = k\bigg( \frac{y_i}{\sum_i y_i \Delta E_i} \bigg) $
$\implies [F] = (\mathrm{\frac{ph}{cm^2s}) (\frac{1}{keV}})$

And we have

$\int_{l}^{h} F(E)\mathrm{d}E = K \sum_i y_i' \Delta E_i$
$= K\frac{1}{y_\mathrm{int}}\sum_i y_i \Delta E_i$
$= K (\frac{y_\mathrm{int}}{y_\mathrm{int}})$
$= K$.

The function would return the differential flux as a function of energy. Do you know what threeml expects?

If this is the convention we are using, then I think the units of $K$ in the current tutorial are not correct, since they are $\mathrm{\frac{1}{cm^2 \ s \ keV}}$.

Please let me know if this all makes sense.

[krishnatejavedula]

@ckarwin @israelmcmc I think this makes sense. I was trying to figure out the units for the differential flux and K too. With $\Delta E$ in the numerator, that would give us a differential flux of $\text{ph} \cdot \text{cm}^{-2} \cdot \text{s}^{-1}$, which is not correct. As per my understanding, it is the amount of energy passing through a unit area per unit time per unit energy interval ($\text{ph} \cdot \text{cm}^{-2} \cdot \text{s}^{-1} \cdot \text{keV}^{-1}$).

Using $F = K y_i' = K \left( \frac{y_i}{\sum_i y_i \Delta E_i} \right)$, we get a K with units $\text{ph} \cdot \text{cm}^{-2} \cdot \text{s}^{-1}$. This would have to be corrected in the notebook too (currently it is /u.cm/u.cm/u.s/u.keV).

This would also mean in the SpecFromDat file, the normalization of flux should be done using dataFlux = dataFlux / sum(dataFlux * ewidths).

[israelmcmc]

Thank you, both. No delta E in the numerator sounds indeed correct. @krishnatejavedula can you please implement this and, as a test, make sure the source injector returns what's expected? I think 3ML expects the output of evaluate() to be differential flux, but I'm on my phone right now.

[krishnatejavedula]

@israelmcmc I believe the output of the evaluate() to be differential flux too. I have already tested the changes with the source injector. I think it's the $\Delta E$ in the numerator is what was causing the issues in its output I mentioned on slack. Rectifying that fixed it and the results are as expected.

[israelmcmc]

Great! Thanks

[israelmcmc]

@krishnatejavedula I opened a PR into your branch for this PR adding a test and syncing with the develop branch. Please take a look. When you incorporate those changes we can merge this PR, assuming the tests will pass.

[israelmcmc]

I'll close and reopen this PR to see if the kick fixes codecov.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 37.05%. Comparing base (1bdbdd9) to head (13c9fd0).
Report is 3 commits behind head on develop.

Files	Coverage Δ
cosipy/threeml/custom_functions.py	27.55% <100.00%> (+6.92%)	⬆️

[israelmcmc]

It worked. I am merging this now.