Add coordsys, proj and binsz parameters to comobsres script

Jürgen Knödlseder
1 parent 10e2ae51
Showing 4 changed files with 39 additions and 7 deletions
doc/source/users/reference_manual/comobsres.rst
modules/comscripts/comobsres.par
modules/comscripts/comobsres.py
test/test_comobsres.py
@@ -36,10 +36,19 @@ General parameters
 ``(armmax = 3.0) [real]``
     Maximum ARM (deg).
  
+``coordsys <CEL|GAL> [string]``
+    Coordinate system (CEL - celestial, GAL - galactic).
+
+``proj <AIT|AZP|CAR|GLS|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``(margin = 30.0) [real]``
     Sky map margin (deg). This margin is added to the data cube size so that
     also source positions outside the range spanned by the data cube are considered.
  
+``(binsz = 1.0) [real]``
+    Sky map bin size (deg).
+
 ``(grouping = 1) [integer]``
     Number of Chi/Psi bins to group for residual computation.
  
@@ -27,7 +27,10 @@ outmap,    f, a, resmap.fits,,,&quot;Output residual map&quot;
 algorithm,  s, a, SIGNIFICANCE,SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE,, "Residual map computation algorithm"
 armmin,     r, h, -3.0,,, "Minimum ARM (deg)"
 armmax,     r, h,  3.0,,, "Maximum ARM (deg)"
+coordsys,   s, a, GAL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,       s, a, TAN,AIT|AZP|CAR|GLS|MOL|SFL|SIN|STG|TAN,, "Projection method"
 margin,     r, h, 30.0,,, "Sky map margin (deg)"
+binsz,      r, h, 1.0,0.1,10.0, "Bin size (deg)"
 grouping,   i, h, 1,1,80, "Number of Chi/Psi bins to group for residual computation"
 dri,        b, h, no,,, "Compute DRI residuals?"
 png,        b, h, no,,, "Generate residual histogram png files?"
@@ -145,7 +145,10 @@ class comobsres(ctools.csobservation):
         self['algorithm'].string()
         self['armmin'].real()
         self['armmax'].real()
+        self['coordsys'].string()
+        self['proj'].string()
         self['margin'].real()
+        self['binsz'].real()
  
         # Query PNG parameters
         if self['png'].boolean():
@@ -187,6 +190,11 @@ class comobsres(ctools.csobservation):
         glat_min = None
         glat_max = None
  
+        # Get parameters
+        proj     = self['proj'].string()
+        coordsys = self['coordsys'].string()
+        binsz    = self['binsz'].real()
+
         # Loop over all input observations
         for obs in inobs:
  
@@ -243,15 +251,22 @@ class comobsres(ctools.csobservation):
         # Compute width and number of pixels
         glon_width = glon_max - glon_min
         glat_width = glat_max - glat_min
-        nx         = int(glon_width+0.5)
-        ny         = int(glat_width+0.5)
+        nx         = int(glon_width/binsz+0.5)
+        ny         = int(glat_width/binsz+0.5)
         if nx > 360:
             nx = 360
         if ny > 180:
             ny = 180
  
+        # Translate to celestial coordinates if requested
+        if coordsys == 'CEL':
+            centre = gammalib.GSkyDir()
+            centre.lb_deg(glon, glat)
+            glon = centre.ra_deg()
+            glat = centre.dec_deg()
+
         # Allocate map
-        map = gammalib.GSkyMap('TAN', 'GAL', glon, glat, -1.0, 1.0, nx, ny)
+        map = gammalib.GSkyMap(proj, coordsys, glon, glat, -binsz, binsz, nx, ny)
  
         # Return map
         return map
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the comobsres script.
 #
-# Copyright (C) 2021 Juergen Knoedlseder
+# Copyright (C) 2021-2022 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -75,7 +75,8 @@ class Test(test):
  
         # Setup command
         cmd = comobsres+' inobs="'+self._obs+'" inmodel="'+self._model+'" '+ \
-                        'algorithm="SIGNIFICANCE" margin="0.0" '+ \
+                        'algorithm="SIGNIFICANCE" margin="0.0" binsz="5.0" '+ \
+                        'coordsys="GAL" proj="TAN" '+ \
                         'outfolder="comobsres_cmd1" outmap="comobsres_cmd1.fits" '+ \
                         'logfile="comobsres_cmd1.log" chatter=1'
  
@@ -88,7 +89,8 @@ class Test(test):
  
         # Setup command
         cmd = comobsres+' inobs="file_that_does_not_exist.xml" inmodel="'+self._model+'" '+ \
-                        'algorithm="SIGNIFICANCE" margin="0.0" '+ \
+                        'algorithm="SIGNIFICANCE" margin="0.0" binsz="5.0" '+ \
+                        'coordsys="GAL" proj="TAN" '+ \
                         'outfolder="comobsres_cmd2" outmap="comobsres_cmd2.fits" '+ \
                         'logfile="comobsres_cmd2.log" chatter=1'
  
@@ -112,7 +114,10 @@ class Test(test):
         add['inobs']     = self._obs
         add['inmodel']   = self._model
         add['algorithm'] = 'SIGNIFICANCE'
+        add['coordsys']  = 'GAL'
+        add['proj']      = 'TAN'
         add['margin']    = 0.0
+        add['binsz']     = 5.0
         add['outfolder'] = 'comobsres_py1'
         add['outmap']    = 'comobsres_py1.fits'
         add['logfile']   = 'comobsres_py1.log'
@@ -138,7 +143,7 @@ class Test(test):
         map = gammalib.GSkyMap(filename)
  
         # Check that there is one observation
-        self.test_value(map.npix(), 3025, 'Check for number of pixels in sky map')
+        self.test_value(map.npix(), 121, 'Check for number of pixels in sky map')
  
         # Return
         return