search

fermi-lat / Likelihood

BSD 3-Clause "New" or "Revised" License
2 stars 1 forks source link

gtmodel mangles headers for CAR projection #123

Open kerrm opened 2 months ago

kerrm commented 2 months ago

I made a model map using gtmodel with the CAR projection:

gtmodel srcmaps=ext_4FGL_J1409.1-6121e_ccube_CAR.fits srcmdl=ext_4FGL_J1409.1-6121e.xml irfs=P8R3_SOURCE_V3 evtype=32 expcube=lt_summed.fits bexpmap=ext_4FGL_J1409.1-6121e_ecube_CAR.fits outtype=CCUBE edisp_bins=-2

The header of the resulting FITS file, pasted below, is correct "at the top" but is later clobbered by what seems to be the header from the counts cube. (Repeated cards start at 2nd instance of "SIMPLE".) The 2nd entry for BITPIX is 32, and when you use astropy to load this file, it **very incorrectly interprets the map data as integers. ds9, on the other hand, uses the first instance (BITPIX=-32) and thus correctly loads the data as floats.

Now, I don't know what the "correct" behavior is with regards to duplicate cards as far as the FITS standard goes, but it's pretty clear that the 2nd values are incorrect and shouldn't be there, so that's the bug.

SIMPLE  =                    T / File conforms to NOST standard                 
BITPIX  =                  -32 / Bits per pixel                                 
NAXIS   =                    3 / No data is associated with this header         
NAXIS1  =                  200 / Length of data axis 1                          
NAXIS2  =                  200 / Length of data axis 2                          
NAXIS3  =                   16 / Length of data axis 3                          
EXTEND  =                    T / Extensions may be present                      
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
CTYPE1  = 'RA---CAR'           / RA---%%%, %%% is the projection, e.g., AIT     
CRPIX1  =                100.5 / Reference pixel                                
CRVAL1  =              212.294 / RA at the reference pixel                      
CDELT1  =                 -0.1 / X-axis incr per pixel at ref pixel (deg)       
CUNIT1  = 'deg     '           / Physical unit of X-axis                        
CTYPE2  = 'DEC--CAR'           / DEC---%%%, %%% is the projection, e.g., AIT    
CRPIX2  =                100.5 / Reference pixel                                
CRVAL2  =              -61.353 / DEC at the reference pixel                     
CDELT2  =                  0.1 / Y-axis incr per pixel at ref pixel (deg)       
CUNIT2  = 'deg     '           / Physical unit of Y-axis                        
CTYPE3  = 'photon energy'      / log_MeV                                        
CRPIX3  =                   1. / Reference pixel                                
CRVAL3  =                1000. / energy at the reference pixel                  
CDELT3  =     333.521432163324 / z-axis logrithmic incr per pixel               
CUNIT3  = 'log_MeV '           / Physical unit of Y-axis                        
CROTA2  =                   0. / Image rotation (deg)                           
DATE    = '2024-04-14T18:06:14' /                                               
TELESCOP= 'GLAST   '           / Name of telescope generating data              
INSTRUME= 'LAT     '           / Name of instrument generating data             
DATE-OBS= '        '           / Start Date and Time of the observation (UTC)   
DATE-END= '        '           / End Date and Time of the observation (UTC)     
EQUINOX =                2000. / Equinox of RA & DEC specifications             
CREATOR = 'gtmodel '           / Software and version creating file             
HISTORY                   $Id: LatCountsMapTemplate,v 1.2 2004/09/24 03:54:20 jc
HISTORY hiang E                                                                 
CHECKSUM= '0000000000000000'   / HDU checksum updated 2024-04-14T22:06:14       
DATASUM = '1912753061'         / data unit checksum updated 2024-04-14T22:06:14 
SIMPLE  =                    T / File conforms to NOST standard                 
BITPIX  =                   32 / Bits per pixel                                 
EXTEND  =                    T / Extensions may be present                      
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
CTYPE1  = 'RA---CAR'           / RA---%%%, %%% represents the projection method 
CRPIX1  =                100.5 / Reference pixel                                
CRVAL1  =              212.294 / RA at the reference pixel                      
CDELT1  =                 -0.1 / X-axis incr per pixel of physical coord at posi
CUNIT1  = 'deg     '           / Physical unit of X-axis                        
CTYPE2  = 'DEC--CAR'           / DEC---%%%, %%% represents the projection method
CRPIX2  =                100.5 / Reference pixel                                
CRVAL2  =              -61.353 / DEC at the reference pixel                     
CDELT2  =                  0.1 / Y-axis incr per pixel of physical coord at posi
CUNIT2  = 'deg     '           / Physical unit of Y-axis                        
CROTA2  =                   0. / Image rotation (deg)                           
CTYPE3  = 'Energy  '           / Energy                                         
CRPIX3  =                    1 / Reference pixel                                
CRVAL3  =                1000. / Energy at the reference pixel                  
CDELT3  =     333.521432163324 / Z-axis incr per pixel of physical coord at posi
CUNIT3  = 'MeV     '           / Physical unit of Z-axis                        
DATE    = '2024-04-14T16:56:15' / file creation date (YYYY-MM-DDThh:mm:ss U     
FILENAME= 'ext_4FGL_J1409.1-6121e_ccube_CAR.fits' /                             
TELESCOP= 'GLAST   '           / name of telescope generating data              
INSTRUME= 'LAT     '           / name of instrument generating data             
DATE-OBS= '2008-08-04T15:43:36.0000' / start date and time of the observation (U
DATE-END= '2022-08-02T21:53:14.9999' / end date and time of the observation (UTC
TIMEUNIT= 's       '           / units for the time related keywords            
TIMEZERO=                   0. / clock correction                               
TIMESYS = 'TT      '           / type of time system that is used               
TIMEREF = 'LOCAL   '           / reference frame used for times                 
CLOCKAPP=                    F / whether a clock drift correction has been appli
GPS_OUT =                    F / whether GPS time was unavailable at any time du
EQUINOX =                2000. / Equinox of RA & DEC specifications             
OBSERVER= 'Peter Michelson'    / GLAST/LAT PI                                   
CREATOR = 'gtbin   '           / Software and version creating file             
HISTORY                   LatCountMapTemplate,v 1.3 2005/04/05 21:06:39 peachey 
HISTORY Exp                                                                     
CHECKSUM= '7bXYAZXV8aXVAYXV'   / HDU checksum updated 2024-04-14T22:06:14       
DATASUM = '351929  '           / data unit checksum updated 2024-04-14T20:56:15
HISTORY The following history was copied from input files by gtbin              
HISTORY ------------------------------------------------------------------------
HISTORY BEGIN history copied from /tank/kerrm/fermi_data/fgl_diffuse/data_4fgl/p
HISTORY sf3/ft1_psf3_zmax90.fits[EVENTS]                                        
HISTORY ------------------------------------------------------------------------
HISTORY Input merit file: /lustre/ki/pfs/fermi_scratch/glastmp/P305-FITS.mergeRu
HISTORY n/0.0/IN/r0239557414_v301_merit.root                                    
HISTORY Filter string: (((FswGamState==0)&&((GltGemSummary&0x60)==0)&&(TkrNumTra
HISTORY cks > 0))&&((-log10(1.0-WP8CTAllProb)) >= (0.010000)*(EvtJointLogEnergy 
HISTORY < 1.250000) + ((EvtJointLogEnergy <= 1.750000)*((0.010000)*(1.0)+(0.0000
HISTORY 00)*(pow((EvtJointLogEnergy-1.250000)/0.500000,1))+(0.018669)*(pow((EvtJ
HISTORY ointLogEnergy-1.250000)/0.500000,2)))+((EvtJointLogEnergy > 1.750000)*(E
HISTORY vtJointLogEnergy <= 2.250000))*((0.028669)*(1.0)+(0.037338)*(pow((EvtJoi
HISTORY ntLogEnergy-1.750000)/0.500000,1))+(-0.017111)*(pow((EvtJointLogEnergy-1
HISTORY .750000)/0.500000,2)))+((EvtJointLogEnergy > 2.250000)*(EvtJointLogEnerg
HISTORY y <= 2.750000))*((0.048897)*(1.0)+(0.003117)*(pow((EvtJointLogEnergy-2.2
HISTORY 50000)/0.500000,1))+(0.001967)*(pow((EvtJointLogEnergy-2.250000)/0.50000
HISTORY 0,2)))+((EvtJointLogEnergy > 2.750000)*(EvtJointLogEnergy <= 3.250000))*
HISTORY ((0.053980)*(1.0)+(0.007050)*(pow((EvtJointLogEnergy-2.750000)/0.500000,
HISTORY 1))+(-0.003525)*(pow((EvtJointLogEnergy-2.750000)/0.500000,2)))+((EvtJoi
HISTORY ntLogEnergy > 3.250000)*(EvtJointLogEnergy <= 3.750000))*((0.057505)*(1.
HISTORY 0)+(0.000000)*(pow((EvtJointLogEnergy-3.250000)/0.500000,1))+(0.121963)*
HISTORY (pow((EvtJointLogEnergy-3.250000)/0.500000,2)))+((EvtJointLogEnergy > 3.
HISTORY 750000)*(EvtJointLogEnergy <= 4.250000))*((0.179468)*(1.0)+(0.243925)*(p
HISTORY ow((EvtJointLogEnergy-3.750000)/0.500000,1))+(0.493075)*(pow((EvtJointLo
HISTORY gEnergy-3.750000)/0.500000,2)))+((EvtJointLogEnergy > 4.250000)*(EvtJoin
HISTORY tLogEnergy <= 4.750000))*((0.916468)*(1.0)+(1.230076)*(pow((EvtJointLogE
HISTORY nergy-4.250000)/0.500000,1))+(-0.501532)*(pow((EvtJointLogEnergy-4.25000
HISTORY 0)/0.500000,2)))+(EvtJointLogEnergy > 4.750000)*((1.645012)*(1.0)+(0.227
HISTORY 011)*(pow((EvtJointLogEnergy-4.750000)/0.500000,1))+(0.029483)*(pow((Evt
HISTORY JointLogEnergy-4.750000)/0.500000,2))))*(EvtJointLogEnergy >= 1.250000 &
HISTORY & EvtJointLogEnergy <= 5.750000) + (2.216967)*(EvtJointLogEnergy > 5.750
HISTORY 000))) || (((-log10(1.0-WP8CTCalTkrProb)) >=  -0.035931+ (   0.27029)*po
HISTORY w(EvtJointLogEnergy,1)+(  -0.23349)*pow(EvtJointLogEnergy,2)+ (   0.0786
HISTORY 7)*pow(EvtJointLogEnergy,3)+ (-0.005295)*pow(EvtJointLogEnergy,4))&&((Fs
HISTORY wGamState==0)&&((GltGemSummary&0x60)==0)&&(TkrNumTracks > 0)&&EvtCalCsIR
HISTORY Ln>4 && ((log10(max(CalTrackAngle,1E-4))) <= (0.529795)*(EvtJointLogEner
HISTORY gy < 3.000000)+ ((1.0)*((0.529795)*(1.0)+(-1.379791)*(pow((EvtJointLogEn
HISTORY ergy-3.000000)/0.916667,1))+(0.583401)*(pow((EvtJointLogEnergy-3.000000)
HISTORY /0.916667,2))+(-0.075555)*(pow((EvtJointLogEnergy-3.000000)/0.916667,3))
HISTORY ))*(EvtJointLogEnergy>= 3.000000 && EvtJointLogEnergy <= 5.750000) + (-0
HISTORY .398962)*(EvtJointLogEnergy > 5.750000)))&&(WP8CTPSFTail > 0.05 && WP8CT
HISTORY BestEnergyProb > 0.1)) && (EvtElapsedTime >= 239557417)  && (EvtElapsedT
HISTORY ime <= 239558070)                                                       
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /lustre/ki/pfs/fermi_scratch/glastmp/P305-FITS.mergeRun/0.0/foo.fit[EVEN
HISTORY TS][gtifilter("/lustre/ki/pfs/fermi_scratch/glastmp/P305-FITS.mergeRun/0
HISTORY .0/OUT/gll_xp_p305_r0239557414_v305.fit_tempgti")]                      
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /lustre/ki/pfs/fermi_scratch/glastmp/P305-FITS.mergeRun/0.0/OUT/gll_xp_p
HISTORY 305_r0239557414_v305.fit[EVENTS][((EVENT_CLASS&o200) != o0) && gtifilter
HISTORY ()]                                                                     
HISTORY Filter string: ((EVENT_CLASS&o200) != o0) && gtifilter()                
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /lustre/ki/pfs/fermi_scratch/glastmp/P305-FITS.mergeRun/0.0/OUT/foo.fits
HISTORY [EVENTS][gtifilter("/lustre/ki/pfs/fermi_scratch/glastmp/P305-FITS.merge
HISTORY Run/0.0/OUT/gll_ph_p305_r0239557414_v305.fit_tempgti")]                 
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /glast/Data/Flight/Reprocess/P305/ft1/gll_ph_p305_r0239557414_v305.fit[E
HISTORY VENTS][((EVENT_CLASS&o200) != o0) && 10 < ENERGY && ENERGY <= 2000000 &&
HISTORY 239330000 < TIME && TIME <= 239593000 && 0 < ZENITH_ANGLE && ZENITH_ANG 
HISTORY LE <= 105 && gtifilter()]                                               
HISTORY Filter string: ((EVENT_CLASS&o200) != o0) && 10 < ENERGY && ENERGY <= 20
HISTORY 00000 && 239330000 < TIME && TIME <= 239593000 && 0 < ZENITH_ANGLE && ZE
HISTORY NITH_ANGLE <= 105 && gtifilter()                                        
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /tmp/jballet/toto.fits[EVENTS][((EVENT_CLASS&o200) != o0) && 10 < ENERGY
HISTORY && ENERGY <= 2000000 && 0 < TIME && TIME <= 1000000000 && 0 < ZENITH_AN 
HISTORY GLE && ZENITH_ANGLE <= 105 && gtifilter()]                              
HISTORY Filter string: ((EVENT_CLASS&o200) != o0) && 10 < ENERGY && ENERGY <= 20
HISTORY 00000 && 0 < TIME && TIME <= 1000000000 && 0 < ZENITH_ANGLE && ZENITH_AN
HISTORY GLE <= 105 && gtifilter()                                               
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /tmp/jballet/P305_Source_001_01.fits_select[EVENTS][gtifilter("/tmp/jbal
HISTORY let/P305_Source_001_01.fits_tempgti")]                                  
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY /tmp/jballet/P305_Source_001_01.fits[EVENTS][((EVENT_CLASS&o200) != o0) 
HISTORY && 10 < ENERGY && ENERGY <= 2000000 && 0 < TIME && TIME <= 1000000000 &&
HISTORY 0 < ZENITH_ANGLE && ZENITH_ANGLE <= 105 && gtifilter()]                 
HISTORY Filter string: ((EVENT_CLASS&o200) != o0) && 10 < ENERGY && ENERGY <= 20
HISTORY 00000 && 0 < TIME && TIME <= 1000000000 && 0 < ZENITH_ANGLE && ZENITH_AN
HISTORY GLE <= 105 && gtifilter()                                               
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY ft1/P305_Source_001_zmax105.fits[EVENTS][((EVENT_CLASS&o200) != o0) && 5
HISTORY 0 < ENERGY && ENERGY <= 1000000 && 0 < ZENITH_ANGLE && ZENITH_ANGLE <= 9
HISTORY 0 && gtifilter()]                                                       
HISTORY Filter string: ((EVENT_CLASS&o200) != o0) && 50 < ENERGY && ENERGY <= 10
HISTORY 00000 && 0 < ZENITH_ANGLE && ZENITH_ANGLE <= 90 && gtifilter()          
HISTORY CFITSIO used the following filtering expression to create this table:   
HISTORY ../ft1_zmax90.fits[EVENTS][((EVENT_CLASS&o200) != o0) && ((EVENT_TYPE&o4
HISTORY 0) != o0) && 1000 < ENERGY && ENERGY <= 10000 && 0 < TIME && TIME <= 681
HISTORY 170000 && 0 < ZENITH_ANGLE && ZENITH_ANGLE <= 90 && gtifilter()]        
HISTORY Filter string: ((EVENT_CLASS&o200) != o0) && ((EVENT_TYPE&o40) != o0) &&
HISTORY 1000 < ENERGY && ENERGY <= 10000 && 0 < TIME && TIME <= 681170000 && 0  
HISTORY < ZENITH_ANGLE && ZENITH_ANGLE <= 90 && gtifilter()                     
HISTORY ------------------------------------------------------------------------
HISTORY END copied history                                                      
HISTORY ------------------------------------------------------------------------
DSTYP1  = 'TIME    '                                                            
DSUNI1  = 's       '                                                            
DSVAL1  = 'TABLE   '                                                            
DSREF1  = ':GTI    '                                                            
DSTYP2  = 'ENERGY  '                                                            
DSUNI2  = 'MeV     '                                                            
DSVAL2  = '1000:10000'                                                          
DSTYP3  = 'ZENITH_ANGLE'                                                        
DSUNI3  = 'deg     '                                                            
DSVAL3  = '0:90    '                                                            
DSTYP4  = 'IRF_VERSION'                                                         
DSUNI4  = 'DIMENSIONLESS'                                                       
DSVAL4  = 'P8R3_SOURCE_V3 (PSF)'                                                
DSTYP5  = 'BIT_MASK(EVENT_CLASS,128,P8R3)'                                      
DSUNI5  = 'DIMENSIONLESS'                                                       
DSVAL5  = '1:1     '                                                            
DSTYP6  = 'BIT_MASK(EVENT_TYPE,32,P8R3)'                                        
DSUNI6  = 'DIMENSIONLESS'                                                       
DSVAL6  = '1:1     '                                                            
NDSKEYS =                    6
kerrm commented 2 months ago

PS -- in case it's useful, when you use HEALPix projection, the header is correct.

Areustle commented 2 months ago

Thanks @kerrm , I'll look into this. Have you seen the behavior with any other tools?

kerrm commented 2 months ago

I haven't seen it in other tools -- this was very ad hoc, though, and I only noticed it because BITPIX BIT me.

Philippe kindof reproduced the issue with the same version of fermitools (2.2.0). For him, the header is still mangled, but the 2nd value of BITPIX is -32, as it should be! Ultimately the real problem is the repeated values, which for these keywords breaks FITS compliance in any case.