Join GitHub today
GitHub is home to over 28 million developers working together to host and review code, manage projects, and build software together.Sign up
Change IRF extrapolation behaviour #1195
Last week @bkhelifi reported that CTA DC-1 AGN spectrum analysis failed miserably.
The core of the issue was identified by @registerrier to lie in negative effective area and energy dispersion from the evaluate method in our IRF classes:
As demonstrated there, what we currently do is to place the interpolation nodes at the bin centers, i.e. at 0.5, 1.5, 2.5, ... deg in FOV offset, and thus since the AGN were placed at offset = 0.0 deg, we extrapolate the IRF values from the nodes at 0.5 and 1.5 deg, and this can lead to negative EDISP and predicted counts and analysis.
This can be fixed in several ways, e.g. by padding with an extra row of interpolation nodes at offset = 0 deg with the same values as at 0.5 deg on IRF read, or by implementing a custom IRF interpolator that has axis-specific behaviour and knows about energy and offset axes and how to extend.
However, the simplest solution is to clip the IRF values to >= 0, i.e. to replace negative values with zero. Doing that is also exactly what Gammalib / ctools does (linear extrapolation followed by clipping to >=0), so at least for CTA DC-1 where simulation was done by ctools, it means that we should be able to analyse the data and get pretty good results. I think there will still be issues from IRF and sampling noise when doing analysis, but that remains to be seen / studies in the coming weeks, it's not clear at the moment.
So I plan to make a pull request shortly where I introduce the clip to >= 0 shortly, applying it to all IRFs (AEFF, EDISP, BKG, table PSF)
I added the clipping to >= 0 in NDData.evaluate, attached as a commit here: a19e769
Locally I see three test fails, because results in high-level analyses that evaluate IRFs changed very slightly: https://gist.github.com/cdeil/50b29850f074af1c898e1c663f7dac88
I plan to add some unit tests here for NDData and IRF evaluation that establish the behaviour concerning extrapolation and clipping to >= 0, but the fix seems OK to me. If we have cases later where we need to have negative outputs from NDData.evaluate, we could add a configure option to allow it in those cases later, no?