MapMaker max_offset cutout size

[cdeil]

I tried using MapMaker to make an all-sky AIT map for CTA DC1 today.

One issue I noticed is that the max_offset isn't used in the way I had hoped / expected.

As can be seen in the screenshot below, what happens is that for most parts of the image where WCS distortions occur, the cutout is chosen too small, the cutout images are truncated.

One thing we could do is change the algorithm that computes cutouts to something that doesn't truncate. E.g. we could compute a sky distance image, threshold it, and then compute the bounding box for that mask. Or we could create a polygon that traces the max offset cone, and then compute the cutout from the polygon.

The sky image distance image based method would even give "fully correct" results in cases where the FOV wraps to the opposite side of the image. That does happen in all-sky analysis, and others have encountered problems with such cases (see here or here). Of course, in that case, the bounding box would be huge, spanning almost the entire image. But a bit slow and correct is better, no?
The computation of the sky image distance images and bounding boxes could be optimised if it turns out to be a real slowdown.

So bottom line: I suggest to change the cutout slice computation to be distance-imaged based for now to avoid this truncation effect.

We could change to completely skycoord-array based eventually, then no cutout box is needed at all, since here in MapMaker no PSF is involved. Then it would work for HEALPix directly.

@adonath - Could it be that we have similar truncation issues in the model image computations, where cutouts are used? We probably need to add tests for exotic cases some day?

@adonath @registerrier @AtreyeeS - Thoughts?

[registerrier]

Hi @cdeil, I would say that this is a use case that is not supported by MapMaker. One motivation for the cutout is to have the ability to perform operations on individual observations such as ring background, FoV normalization and possibly to store the results on observation per observation basis.

That said, computing a set of valid pixels based on some region selection, and fill the maps only for those pixels can be a good approach. I think initially it worked like this rather than by multiplying by a mask.

Regarding model evaluation, MapEvaluator does not really support strongly distorted projections, because of PSF convolution and solid angle calculation.

[cdeil]

this is a use case that is not supported by MapMaker

Surely making an all-sky map with Gammapy should be a supported use case?
If we don't offer correct cutouts if max_offset is given, then the user has to basically script the MapMaker again, i.e. 100 lines.

One motivation for the cutout is to have the ability to perform operations on individual observations such as ring background, FoV normalization and possibly to store the results on observation per observation basis.

Well we already have this: https://docs.gammapy.org/0.11/api/gammapy.cube.MapMakerRing.html
I'm not happy about this, maintaining multiple map makers with a lot of copy & pasted code and docstrings is a pain. Like I've said before, my suggestion here would be to change to one maker per quantity, and then having multiple specialised makers for background estimation is fine, that's where the 90% of the complexity of code and map making in Gammapy will be.

That said, computing a set of valid pixels based on some region selection, and fill the maps only for those pixels can be a good approach.

It could work like that for most cases, and HEALPix would work out of the box. Except if background estimation needs a convolution operation, like the ring.

I think initially it worked like this rather than by multiplying by a mask.

I don't think so, already the first version or MapMaker, and the old BasicImageEstimator used cutouts

• https://docs.gammapy.org/0.8/api/gammapy.cube.MapMaker.html
• https://docs.gammapy.org/0.7/api/gammapy.image.BasicImageEstimator.html

Regarding model evaluation, MapEvaluator does not really support strongly distorted projections, because of PSF convolution and solid angle calculation.

Yes, at least for PSF-convolved a 2D image is needed. Eventually I think being able to make an all-sky model image is also a supported use case for Gammapy. So that's an argument against changing to 1D pixel arrays, and to keep using 2D cutouts.

---

For now, what I've done is to change to offset_max="10 deg" (see here) - that fixes the cutout truncation issue and for my use case I don't really need an offset_max. That is, I hope it fixes the issue - it's very slow, it takes 1-2 hours to make an all-sky map for the 1DC data.

One concrete proposal I still have: add a "correct" cutout method to [gammapy.maps.WcsGeom.cutout] (https://docs.gammapy.org/0.11/api/gammapy.maps.WcsGeom.html#gammapy.maps.WcsGeom.cutout) which implements what I suggest above: make distance image and find cutout that contains all pixels within the given radius. This could then be the default in MapMaker by changing one line. It would give correct results for all-sky, and for most maps we have so far would be almost the same - except to add a pixel row or two at the edge where currently we truncate too early.

[cdeil]

This is what happened with max_offset="10 deg".

Events were only simulated out to 5 deg, and background IRF extrapolation gives incorrect background at large offsets (black ring = negative excess = background too high).

So I'll go back to max_offset="5 deg" for now, accepting the clipping as the lesser evil.

Another issue is that excess is overall positive, i.e. background underestimated in some of the AGN fields. That looks like a bug in our solid angle computation. I think https://docs.gammapy.org/0.11/_modules/gammapy/maps/wcs.html#WcsGeom.solid_angle is incorrect if pixels aren't approx squares on the sphere, but parallelograms like in the case of AIT?
Another possible bug (or a second one) is that the northern background IRFs aren't correct somehow. This image is with 10 bins / decade for the background computation, and above 300 GeV, so energy integration should be ok.

Scripts and files are here: https://github.com/gammasky/cta-dc1-gps-analysis/blob/master/make_new.py

I wanted to try ctools and compare, starting with a counts image which might already differ (see here. But at least ctools 1.5 doesn't seem to work at all with AIT if there are non-sphere pixels in the image (see here. Still, comparisons could be made, e.g. focusing on some AGN FOV where AIT distortion is present, just a test with 1 or 2 runs should reveal all issues / differences.

[cdeil]

@adonath - Thanks for fixing the solid angle computation in #2276 .

I re-ran the CTA DC-1 all-sky map with the script here and the result is here and here.

It looks much better now, but there are still some visible issues for the fields of view at the edge of the all-sky image (also see screenshot pasted below). There is FOV cutout truncation (the main issue reported and discussed here), but also an issue of bright excess lines at the edge of the all-sky image, probably due to partially inside / outside pixels being filled with counts, but not with background, because there they get nan because one of the corners is outside. Maybe we should also zero out counts there, but then we'd waste time when an image in the common case where image have no such "outside" regions? There's also other smaller background estimation errors, at the edge of AGN exposures, or above the GC in the 1 TeV map, there the background is overestimated for unknown reasons.

I'll leave this issue open for now - would like to improve FOV cutout at least, and maybe consider doing something about the edge of the all-sky image to get a good result there also.

[adonath]

In the new data reduction scheme the cutout parameter is separated from the offset-max. See here for an example: https://docs.gammapy.org/dev/cube/index.html#combining-data-reduction-steps. With this separation it should be possible to resolve the issue described here. Needs verification but not before v0.15. Re-labelling to 1.1.

[registerrier]

Is there any left to do here? @adonath ?

[Astro-Kirsty]

Any update here @adonath?

[adonath]

Thanks @Astro-Kirsty, I think we can close the issue for now. Maybe we we will compute an all-sky map for the next CTA data challenge. If we still find problems, we can open new issues.