Implementing methods for building background IRFs from observed DL2 event lists

[chaimain]

Hello everyone,
While discussing with @moralejo, @mstrzys and @rlopezcoto regarding creating a background model for LST-1 data, we thought it would be good to have a function to separate the source events from a given DL2 events table and then use pyirf's background_2d function on the remaining background events to generate an observed data based background model for creating skymaps and further analyses.

Will it be ok to have such a function in pyirf itself?

[maxnoe]

I don't really see why the current function isn't already enough for that. What is missing?

Can't you just give the table of observed events to the current function?

[chaimain]

I do not mean to change the current function, but to have a separate function, which selects the background region from an ON/Wobble observation, based on a method of selection, the user would like to use, and then the user can combine it with the current background function to continue with their analysis.

[maxnoe]

You mean extracting counts in the off regions? That's more a job for gammapy, e.g. the reflected regions.

For point-like analysis, I don't think any background IRF is needed. For Full-Enclosure, the story is different, and then you have to provide background rates for the full FoV, which requires a measurement without a source in the field of view or code that somehow fills the region where you expect the source with something else.

Is this what you are talking about?

[moralejo]

We'd like to also have the possibility of making a background model from the very observation that you intend to analyze, but somehow masking certain areas where known sources exist. Either that, or from dedicated off observations, we think building models from real data is useful because for long observations it would be hard, perhaps impossible (due to strong background suppression), to do it based on MC alone.

[maxnoe]

@moralejo This is exactly what https://docs.gammapy.org/dev/makers/reflected.html is for.

What would be gained by duplicating this kind of functionality in pyirf, when the analysis of the DL3 data is done in gammapy anyways?

[mstrzys]

Yes, it is for the full-enclosure IRFs. No, you do not necessarily need the Off measurements in case of wobble observations. The idea is to add functionalities like the ones we use in MAGIC. I think this is what you refer to in your second to last sentence, Max.

Sure, one is free to add this code at different places (pyirf, DL2toDL3, gammapy). During our discussion we considered pyirf to be most plausible place for them. For sure it would avoid code duplication, as the functions would be available to any other analysis program (not just gammapy) and could be used by other telescope types as well (not just LST). In any way, I cannot see any harm in having such a wrapper or do you see possible damage arising from it?

[maxnoe]

I see this kind of functionality in the scope of the science tools, not here. And I think large parts of that are already implemented in those.

Maybe @adonath can comment from the gammapy point of view.

See for example gammapy background makers:

• https://docs.gammapy.org/dev/makers/ring.html
• https://docs.gammapy.org/dev/makers/reflected.html

I suspect ctools has equivalent functionality.

[mstrzys]

Well, as the open gamma ray format states Background can either be considered an IRF or not. The open gamma ray data format considers the former and thus you have a functionality for it in pyirf, don't you?

Later on I can image that one would optimise/cross-check e.g. a MC obtained/mathematical description of the acceptance function with real data before storing it to file. Both would require having the data based map at this earlier stage. Furthermore, there are background methods which need to be applied before the cuts and thus can never be applied at the science tool stages (or a rework of pyirf without cuts is needed). Again, there seems to be rather potential than harm.

[chaimain]

The gammapy background maker functions do help with the dataset reduction and providing the appropriate exposure to proceed with the spectral analysis, but we would like to create a general background model, not just from the OFF observations (as given in this gammapy tutorial), but also from the ON/Wobble observations, and then use the similar methods as in the gammapy background maker functions, or something different, to select the OFF counts. And for this, I think it is better to have a function to separate such OFF counts before generating the DL3 file and then generating a background model based on a large dataset of such observations.

Maybe @adonath can provide clarification on this from gammapy's perspective.

[adonath]

I read through the comments so far and had a hard time to understand what is planned to achieve and what is currently missing functionality for this. In general it seemed to me there was a misunderstanding / confusion of methods. So maybe I start with a short overview of what I understood so far:

• @chaimain From your initial post it seemed like what you wanted to achieve is, what I know as "template background": the background template is created from DL2 event files and the events that don't pass G/H selection cuts. This methods is described e.g. in https://arxiv.org/pdf/astro-ph/0310025.pdf However in addition you would like to exclude certain regions from the background estimate. This method can only be implemented on the pyirf level, as it requires access to DL2 event lists.
• The other currently used method is what I would call "FoV background": the background IRF template / acceptance is created from multiple off observations in zenith bands, but again from DL3 events passing G/H selection cuts. This is what is demonstrated in Gammapy notebook: https://docs.gammapy.org/0.18.2/tutorials/background_model.html However from the ST perspective, we prefer if this model was just delivered as an "IRF" and created by pyirf. This is how it is e.g. handled in HESS. This method is e.g. described in https://arxiv.org/pdf/1910.08088.pdf.
• The remaining methods are ring and reflected regions background and on/off observations, where the background is estimated from a single off observation or off regions from DL3 event lists or binned data. Those methods are implemented in Gammapy. See e.g. https://arxiv.org/abs/astro-ph/0610959

@moralejo @chaimain @mstrzys From what I understood you would like to create the template from DL2 event files and events that don't pass G/H selection cuts, but in addition have the possibility to exclude certain regions (option 1 above). Is this correct?

[maxnoe]

From your initial post it seemed like what you wanted to achieve is, what I know as "template background": the background template is created from DL2 event files and the events that don't pass G/H selection cuts

I am confused, of what use is a background model of events that do not pass the cuts?

I think we are talking about 2).

[adonath]

@maxnoe This is much simplified and I don't know the details of the method, but with certain corrections I think one can basically create a spatial template just as in method 2, but from the events classified as hadrons. The normalisation / total rate is then adjusted to predict the correct number of background events outside exclusion regions for the events classified as gammas.

[kosack]

I am confused, of what use is a background model of events that do not pass the cuts?

That's another method entirely - using data from an off-region in parameter space in the same sky region as the on-region. I don't think that is what is being asked for here: just ways to make a background model. The background model is really part of the science tools, I think, it doesn't really fit in the scope of pyirf (though of course both might be used together).

https://www.aanda.org/articles/aa/abs/2003/40/aah4546/aah4546.html

[mstrzys]

Generally, one can see the background as part of the IRFs or not. I think the general confusion comes from the fact that we deal with irreducible background from hadrons (rather part of IRFs) and diffuse gamma-ray emission (certainly part of the science tools). This is a somewhat metaphysical discussion and up to personal taste. Since the gamma ray data format does provide definitions of the background and pyirf has functionalities to include it, to us it sounded plausible to place these functionalities here.

Our aim is a data based full-enclosure background model possibly avoiding assumption such as radial symmetry. As Axel pointed out there are several methods:

1. slight extension of the reflected region background like e.g. the Wobble/Blind maps from S. Klepser (see fig 1. https://doi.org/10.1051/0004-6361/201833139)
2. template background (1. from Axel's list)
3. FoV background sampling from Off data (2. from Axel's list)
4. an analytical model, which can be fitted to the former 3

I imagine to start with 1 and 4, and add 3. 2 certainly requires additional investigations. These methods have certain strength and weaknesses but at them end aim for the same thing and provide the same output. I am open to your suggestions where the full-enclosure background model shall be added, but there should be one place for it, shouldn't it?

[maxnoe]

@mstrzys Yes, and I don't have strong opinion on where, the argumentation that it should already be part of the DL3 files would point to pyirf.

So if you want to make a pull request implementing any of the mentioned methods, that would be very welcome.

[kosack]

I fully agree this would be nice thing to have, but this all requires spatial regions, binning in WCS coordinates, etc. I think it really is a Science Tool, not an IRF tool. From a CTA point of view, the background model software is a product of SUSS (science tools, quick look, etc) , not DPPS (calibration reconstruction and IRFs). The background models themselves are expected to be in IKC (likely), as they require quite a bit of work and sufficient real data.

Generally making a background model is a large and complex project that requires DL3 events (not DL2) and all the functionality of the science tools. If you want to prototype this, I would suggest making it a separate project that uses pyirf if needed, but is not part of it. This would be a very useful tool, but doesn't belong inside PyIRF. We could create a repo in the CTA GitLab for this project.

Also, you might want to wait for a week or so to see which of the Science Tools are chosen by CTA!

[mstrzys]

@kosack , @maxnoe , @adonath thanks for your comments. Let me just summarize them and please feel free to correct me if I got anything wrong.

It seems we cannot find an agreement regarding the existing packages. Based on CTA guidelines Karl suggests the background models to be part of the science tools or at least that they belong to the same data level. On the other hand on behalf of gammapy Axel mentions that the full-enclosure background is preferred to be provided with the IRFs. Hence, they seem to fall somewhere in between.

In that case, the best long-term solution seem to be a separate package as you suggest, Karl. In that case we make use of some functionalities from pyirf such as cuts, binning, I/O etc. Maybe for the time being we can prototype here and as soon as it working we will get back to you, Karl, regarding an own repo. What do you think about it?

@maxnoe Does your thumb-up mean that you support a separate package or were you referring to another part of Karl's comment, Max?

Side note:
I am afraid that on the long-run it might be problematic to base the background tools on DL3 data only. There are methods which need a wider data set e.g. the template background method and who knows what people come up with in the future. It can work in DL3 only if events are flagged rather than cut, but this is somewhat clashing with the idea of reduced event lists. Thus the software should be flexible to switch between both. If the background comes as a separate package then would be anyway would be and this this part of the discussion can be postponed for the moment.

[adonath]

@kosack I'm not completely familiar with the CTA perspective, but to me it's not 100% obvious that the background modelling should only be part of SUSS. I fully agree that making a background model is a large and complex project, but this in turn means it cannot be the responsibility of science users (at least not a large part) to create such a model.

On the other hand it is clear that part of the background modelling, such as reflected, ring and on/off analysis, always remains on the science tools side. But even the aforementioned methods typically require some knowledge of the hadronic acceptance of the instrument (here I'm using acceptance / background model in a synonymous way). Putting the responsibility to always create an acceptance / background model from off observations in the hand of science users is a recipe for disaster I think. Instead there should be one "standard" model that is delivered with the IRFs (including systematics / priors) and that can be used as a "baseline" background model for any kind of analysis. Whether this will be derived from OFF observations or DL2 templates I don't know. The residual hadronic background I typically consider as a property of the instrument / reconstruction method. I also sometimes like to call it "instrumental background", as compared to "astrophysical background".

In general users are free to create their own instrumental background models as well (but should have a good reason to do so) e.g. from OFF observations and for these cases it would be nice to provide users with a set of tools. Here I agree again that the place where those tools live is not obvious. One can certainly rely on some functionality of the ST e.g. for the access of DL3 events, but typically the models are binned in camera coordinates and not WCS coordinates. As a compromise I agree these tools could live in an independent repository for now. However (and this is my main point) I still think a "baseline" instrumental background model, created e.g. as described here, should be delivered with the IRFs and could be seen in the responsibility of pyirf.

@mstrzys Parametric background models are currently not supported in the ST (at least not in Gammapy) and I doubt this will ever be the case. For H.E.S.S. (see here) the shape of the instrumental background is just too complex to be described by a simple parametric model: it changes from a center filled to ring like shape with energy, it is not radially symmetric etc. However one typically relies on the template and then applies first order parametric corrections in norm, gradient and spectrum of the model on the ST level.

I'll add @lmohrmann to the discussion as he also has a lot of experience with background models for H.E.S.S. and delivering them in the current GADF format.

[maxnoe]

@mstrzys

Maybe for the time being we can prototype here and as soon as it working we will get back to you, Karl, regarding an own repo. What do you think about it?

That doesn't make sense IMHO. Just start im another repo, include pyirf and gammapy as a dependency. I don't have a strong opinion on whether it should go in the cta github orga or the aswg group in the cta gitlab, both are fine.

If you want, I can create the repo here, if you want to use the cta gitlab, you should be able to create a repo in the aswg group yourself or ask orel

[mstrzys]

@maxnoe fine, then I consider your aforementioned suggestion of a pull request outdated and we go independent straightaway. Since, most of the software resides here, I think for the time being the background model tools should do as well. So, could you please create us a repo here - no opinion about the name, but lets go with pybkgmodel for consistency with pyirf.

@adonath thanks for this detailed comment; yep, I could not agree more. Particularly the point that the background model should not be a task for the analyzer; he should have the possibility to generate his own, but a default option should be there - at least if CTA wants to achieve the convenience of an fermi-like analysis. That the science tools will contain some kinds of background tools is also clear as a simple "aperture photometry" approaches need DL3 files only. This - as you rightly point out - should however not conflict with the tools we have in mind as they could serve what you call baseline background model.
Regarding the analytical background description, I can imagine that is rather complicated for an array such as HESS (for MAGIC it is the same). However, for single telescopes like the LST the story might be different - let's see.

[lmohrmann]

Thanks for adding me, @adonath. I completely agree with you that CTA will need to provide with the DL3 data a "default" background model, however that is being generated. I also think that some users might want to generate their own models, and there should be some way of achieving this. Since at least for some methods (template background, as was mentioned) access to DL2 data is needed, I would also tend to argue that developing such methods fully within the framework of the Science Tools does not seem a very good idea. Whether this should be included in pyIRF or put into a dedicated repo, I don't have a strong opinion about.

[maxnoe]

@mstrzys sorry for the delay, here is the repo: https://github.com/cta-observatory/pybkgmodel