Bug #1415

errors with ctobssim simulation for diffuse source with Healpix file

Added by Yang Lili about 9 years ago. Updated over 8 years ago.

Status:ClosedStart date:02/03/2015
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.0.0
Duration:

Description

Use ctobssim for CTA simulation for diffuse sources with Healpix fits file, get the following error:

terminate called after throwing an instance of 'GException::invalid_argument’
what(): *** ERROR in GSkymap::pix2dir(GSkyPixel&): Invalid argument.
A 2-dimensional sky map pixel (298.425,-0.374402) is used to determine the sky direction for the 1-dimensional sky projection “HealPix”
Please specify a 1-dimensional sky map pixel.
Aborted (core dumped)

healpix.fits (200 KB) Yang Lili, 02/03/2015 01:57 PM

diffuse.xml Magnifier (1.06 KB) Yang Lili, 02/03/2015 01:57 PM

diffuse_src.xml Magnifier (725 Bytes) Knödlseder Jürgen, 02/04/2015 12:45 AM

map_car_gal.fits (1.98 MB) Knödlseder Jürgen, 02/04/2015 12:45 AM

wcs.fits (1.99 MB) Yang Lili, 02/04/2015 08:23 AM

wcs_map.png (226 KB) Knödlseder Jürgen, 02/05/2015 06:18 PM

Wcs_map

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen about 9 years ago

  • Assigned To set to Knödlseder Jürgen
  • Target version set to 1.0.0
Just for the record:
  • GammaLib needs the ORDERING keyword, the original file had a ORDER keyword. I should check whether both are “standard”
  • GammaLib needs the PIXTYPE keyword set to Healpix. Also should check if this is standard or if some other way should be okay to identify whether a file is Healpix.
  • GammaLib needs the COORDSYS keyword to know in which coordinate system the map is given. I see no real way around that problem.

Finally, I should put somewhere information about the expected HealPix format so that people can look up what is needed.

#2 Updated by Knödlseder Jürgen about 9 years ago

  • Status changed from New to In Progress
  • % Done changed from 0 to 10

Unfortunately there is some way to go as the actual code has been mainly geared towards WCS maps and not to the usage of HealPix maps.

So far I have identified the following:
  • The GModelSpatialDiffuseMap::mc needs to be adapted to work also on HealPix maps. I have some code working but I cannot yet make sense of the results. I first need to produce a test case.
  • The GSkymap::operator() sky direction operator is only implemented for WCS maps, HealPix does not work. This operator is needed for interpolating within a map. Implementation of this operator for HealPix needs nearest neighbour pixel determination, which in principle is easy in HealPix, but needs some reading to understand how to do it properly

#3 Updated by Knödlseder Jürgen about 9 years ago

Both issues have been fixed. The GSkymap::operator() sky direction operator now also works for HealPix maps, a working example is the conversion from HealPix into another projection:

$ python
>>> import gammalib
>>> healpix=gammalib.GSkymap("healpix.fits")
>>> map=gammalib.GSkymap("CAR","GAL",0.0,0.0,0.5,0.5,720,360)
>>> map += healpix
>>> map.save("map_car_gal.fits",True)

This code projects the HealPix map into a cartesian map and saves the result on disk.

The GModelSpatialDiffuseMap::mc has also been adapted. A working example is:

$ ctobssim
RA of pointing (degrees) (0-360) [128.84] 
Dec of pointing (degrees) (-90-90) [-45.18] 
Radius of FOV (degrees) (0-180) [20.0] 
Start time (MET in s) (0) [0.0] 
End time (MET in s) (0) [1800.0] 
Lower energy limit (TeV) (0) [0.1] 
Upper energy limit (TeV) (0) [100.0] 
Calibration database [dummy] 
Instrument response function [cta_dummy_irf] 
Model [diffuse_src.xml] 
Output event data file or observation definition file [events.fits]

using the attached XML file that contains only the source component. The intensity of the map has been changed so that a few 100 events get simulated. The simulation region is centred on Vela, the brightest source in the map. The simulated events obviously cluster around Vela.

Some more testing is needed to validate all changes.

#4 Updated by Yang Lili about 9 years ago

Thanks for your effort!! I think it would be good if you also put some information about what is needed for WCS format file. I am trying a wcs fits file with ctobssim, it gives the following errors,

  • ERROR encounterted in the execution of ctobssim. Run aborted ...
  • ERROR in GSkymap::alloc_wcs(GFitsImage*): Invalid World Coordinate System code (MOL). Projection code not known. Should be one of AIT/AZP/CAR/MER/TAN/STG.

Could you please indicate what is missing? Thank you in advance.

#5 Updated by Knödlseder Jürgen about 9 years ago

MOL is not yet supported. I will add support for this later.

#6 Updated by Knödlseder Jürgen about 9 years ago

  • % Done changed from 80 to 90

I added the MOL projection (class GWcsMOL). Will be merged into devel within the next hour.

#7 Updated by Knödlseder Jürgen about 9 years ago

Now merged into devel branch.

I created a follow-up action (#1416) to make sure that all missing WCS projections will be implemented at some point in time.

Is this issue now fixed from your side or do you still encounter problems somewhere?

#8 Updated by Yang Lili about 9 years ago

Using ctobssim with wcs fits file, find this error,

  • ERROR encounterted in the execution of ctobssim. Run aborted ...
  • ERROR in GWcsMOL::prj_x2s(int, int, int, int, double*, double*, double*, double*, int*): 2 (x,y) coordinates were invalid.

#9 Updated by Knödlseder Jürgen about 9 years ago

Yang Lili wrote:

Using ctobssim with wcs fits file, find this error,

  • ERROR encounterted in the execution of ctobssim. Run aborted ...
  • ERROR in GWcsMOL::prj_x2s(int, int, int, int, double*, double*, double*, double*, int*): 2 (x,y) coordinates were invalid.

I think something is wrong with that map, I even get an error when I open the file with fv. I see that you have 721 x 361 pixels of 0.5 deg, hence formally this map is bigger than the sky. I also see blank bands on the top and the bottom, which may indicate that these zones lie outside the range of valid coordinates. I don’t know how MOL is defined, but maybe you cannot map an entire sky with that projection.

Have you tried the original HealPix map with the next code? I guess it should work.

#10 Updated by Knödlseder Jürgen about 9 years ago

DS9 says “Off map” for the zones at the top and the bottom. Hence I think that ctobssim says basically the same thing, just in other words.

#11 Updated by Yang Lili about 9 years ago

Have you tried the original HealPix map with the next code? I guess it should work.

Yes, I have tried that file. It works fine.

#12 Updated by Knödlseder Jürgen over 8 years ago

  • Status changed from In Progress to Closed
  • % Done changed from 90 to 100

Close this now as problem was related to an invalid input map.

Also available in: Atom PDF