Incorrect CAR projection in gtobssim/gtbin

[anatolii-zenin]

I have encountered an issue when trying to simulate observation of an extended source with its morphology defined in a .fits file. The problem can be easily seen in this screenshot. The right part shows the input .fits file I provide to gtobssim, and the left part shows the output from gtbin. Unexpectedly, the shapes look different despite the fact their WCS values are the same, as can be seen in the lower part of the screenshot. The distorted image shown on the left changes if the source is moved to another position.

Is it a bug or do I do something wrong?

I use:
fermitools 1.2.1
fermitools-data 0.17

$ uname -a
Darwin zenin.local 17.5.0 Darwin Kernel Version 17.5.0: Mon Mar 5 22:24:32 PST 2018; root:xnu-4570.51.1~1/RELEASE_X86_64 x86_64

All the input and output files can be found here. The set of commands to reproduce the issue:

(fermi) $ gtobssim
WARNING: version mismatch between CFITSIO header (v3.43) and linked library (v3.41).
File of flux-style source definitions[ring_source_def.xml]
File containing list of source names[ring_source_list.txt]
Pointing history file[none]
Prefix for output files[test_ring]
Simulation time (seconds)[86400]
Simulation start time (seconds wrt MET 0)[INDEF]
Apply acceptance cone?[no]
Response functions[P8R3_SOURCE_V2]
Random number seed[293049]
added source "test_source"
Generating events for a simulation time of 86400 seconds....

(fermi) $ gtbin
WARNING: version mismatch between CFITSIO header (v3.43) and linked library (v3.41).
This is gtbin version HEAD
Type of output file (CCUBE|CMAP|LC|PHA1|PHA2|HEALPIX) [CMAP]
Event data file name[test_ring_events_0000.fits]
Output file name[ring_binned_cmap.fits]
Spacecraft data file name[NONE]
Size of the X axis in pixels[150]
Size of the Y axis in pixels[150]
Image scale (in degrees/pixel)[0.2]
Coordinate system (CEL - celestial, GAL -galactic) (CEL|GAL) [CEL]
First coordinate of image center in degrees (RA or galactic l)[270]
Second coordinate of image center in degrees (DEC or galactic b)[30]
Rotation angle of image axis, in degrees[0]
Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN:[CAR]
gtbin: WARNING: No spacecraft file: EXPOSURE keyword will be set equal to ontime.

[nmirabal]

This is likely an effect of gtbin when creating the counts map. Recall that you are simulating photon counts with a resolution in arcminutes. The input fits image will guide the simulation, but in the small count regime any photon excess that does not fall exactly within input ring will guide how the counts map is produced. If you look at the raw events in test_ring_events_0000.fits, or if you set set X,Y = 30 and Image scale = 1 in gtbin. The image will "look" rounder. But when the image is projected into a counts map with 0.1 image scale the result will not look exactly as you input file unless the photon count increases significantly. Again this is a simulation and the actual photon distribution will likely be slightly different and as a result the actual counts map could look rounder or more elliptical.

[akira-okumura]

The ring shape is obviously distorted. Please look at the East edge of the circle in the screen shot below. The circle touches the 30:00.0 grid in the top-right panel (input FITS file) but it does not in the other panels (top-left with low-energy photons, bottom-left with energies > 1 GeV and 10 times long obs. time). I believe it is not an "effect" but a bug.

this is a simulation and the actual photon distribution will likely be slightly different

I do not understand what you explain here. Simulators must simulate exactly what we mathematically expect.

[akira-okumura]

It happens with RA---AIT and DEC--AIT too. I first thought that the tool is handing the CAR projection improperly, but changing CAR to AIT did not help.

[anatolii-zenin]

This is what it looks like with RA = 180 deg and DEC = 70 deg:

However there is no distortion when I place the source it at (0.,0.) coordinates.

[nmirabal]

After discussing it with Seth Digel, we believe the issue is with the definition of the reference coordinates (specifically CRVAL2) in ring_source.fits. Wells et al.(1981) and Calabretta & Greisen (2002) are different -CAR conventions in the fits header. You can read more details in Section 2.8 of https://www.aanda.org/component/article?access=bibcode&bibcode=&bibcode=2002A%2526A...395.1077CFUL

The Fermitools use the Calabretta & Greisen (2002), as a result the reference pixel always needs to be on the equator CRVAL2=0. Then CRPIX2 needs to be compensated for the shift i.e. newCRPIX2 = originalCRPIX2 +(original CVAL2)*CDELT2 = 81.5.

[akira-okumura]

• The difference described in Calabretta & Greisen 2002 is difference between RA/DEC and RA---CAR/DEC--CAR. I think the latter is not introduced in Wells et al. (1981).
• Even if the Fermitools uses a particular WCS definition, the result obtained by gtbin should be consistent with the original diffuse map. Our results shows that gtobssim and gtbin use different WCS interpretations.
• The distortion appears when AIT is used as I wrote before.

[akira-okumura]

The Fermitools use the Calabretta & Greisen (2002), as a result the reference pixel always needs to be on the equator CRVAL2=0. Then CRPIX2 needs to be compensated for the shift i.e. newCRPIX2 = originalCRPIX2 +(original CVAL2)*CDELT2 = 81.5.

This is not the CAR definition in Calabretta & Greisen (2002). The user is allowed to use any value for CRVAL2. Please read the paper again.

[akira-okumura]

Where is the source file of the corresponding part of gtobssim, which reads the FITS image and coordinates headers?

[akira-okumura]

Please see the screen shot that I made with RA---AIT and DEC--AIT headers.

I believe gtobssim completely ignore the projection keywords such as AIT and CAR and it handles any FITS images as if only linear coordinates are used in FITS headers (Wells et al. 1981 style).

[anatolii-zenin]

The issue seems to behave exactly as Prof. Okumura described. Whenever I try to generate an event map with gtobssim, I get results like this:

The shape of this ring does not depend on its coordinates or the coordinates of the reference pixel, which means gtobssim ignores keywords related to projection and treats everything in linear coordinate system.

[Areustle]

This is a fairly tricky problem to solve. The author of the libraries gtobssim is depending upon here to write the resulting file is no longer with the mission. It would likely require substantial human effort to address this bug. We strongly suspect that there is a bug in how gtobbsim interfaces with the MapSources library, as the underlying cause.

An alternative approach, rather than throwing individual photons with gtobssim, would be to perform a binned analysis. Use gtmodel to create a model counts map then apply Poisson statistics to the resulting map. @eacharles can provide more information on that analysis.

[Areustle]

This would be substantially faster than using gtobssim, even without this projection bug.