add SIGNIF option in csresmap. added tests for DIFFSQRT and SIGNIF algorithms

ldivenere
1 parent 4b7d6042
Showing 4 changed files with 71 additions and 2 deletions
cscripts/csresmap.par
cscripts/csresmap.py
test/test_csresmap.py
test/test_python_ctools.py
@@ -43,7 +43,7 @@ yref,      r, a, 22.01,-90,90, &quot;Second coordinate of image center in degrees (DE
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
 nypix,     i, a, 200,,, "Size of the Y axis in pixels"
 binsz,     r, a, 0.02,,, "Pixel size (deg/pixel)"
-algorithm, s, a, SUBDIV,SUB|SUBDIV|SUBDIVSQRT,, "Residual map computation algorithm"
+algorithm, s, a, SUBDIV,SUB|SUBDIV|SUBDIVSQRT|SIGNIF,, "Residual map computation algorithm"
  
 #
 # Standard parameters
@@ -241,6 +241,27 @@ class csresmap(ctools.cscript):
             self._resmap -= modelmap
             self._resmap /= modelmap.sqrt()
  
+        # Calculate significance from Li&Ma derivation
+        elif algorithm == 'SIGNIF':
+            signmap = (self._resmap - modelmap).sign()
+            logmap = self._resmap/modelmap
+            # Masking pixels with 0 counts. Put a small value >0 for log calculation
+            # and restore the 0. value afterwards
+            mask = []
+            for i in range(logmap.npix()):
+                # Testing pixel i and map 0 (there is one map only after stacking)
+                if logmap[i,0] == 0:  
+                    mask.append(i)
+                    logmap[i,0] = 1e-10
+            logmap = logmap.log()
+            for i in mask:
+                logmap[i,0] = 0.
+
+            signif_squared_map = (self._resmap*logmap) + modelmap  - self._resmap
+            signif_squared_map *= 2.
+            unsigned_resmap = signif_squared_map.sqrt()
+            self._resmap = unsigned_resmap * signmap 
+
         # Raise exception if algorithm is unknown
         else:
             raise TypeError('Algorithm "'+algorithm+'" not known')
@@ -165,6 +165,54 @@ class Test(test):
         # Check pull distribution file
         self._check_result_file('csresmap_py2.fits', nx=200, ny=200)
  
+
+        # Run csresmap with SUBDIVSQRT algorithm
+        resmap = cscripts.csresmap()
+        resmap['inobs']     = self._cntcube
+        resmap['modcube']   = 'NONE'
+        resmap['expcube']   = 'NONE'
+        resmap['psfcube']   = 'NONE'
+        resmap['edispcube'] = 'NONE'
+        resmap['bkgcube']   = 'NONE'
+        resmap['caldb']     = self._caldb
+        resmap['irf']       = self._irf
+        resmap['inmodel']   = self._model
+        resmap['outmap']    = 'csresmap_py3.fits'
+        resmap['algorithm'] = 'SUBDIVSQRT'
+        resmap['logfile']   = 'csresmap_py3.log'
+        resmap['chatter']   = 3
+        resmap['publish']   = True
+
+        # Run csresmap script
+        resmap.execute()
+
+        # Check pull distribution file
+        self._check_result_file('csresmap_py3.fits', nx=200, ny=200)
+
+        # Run csresmap with SIGNIF algorithm
+        resmap = cscripts.csresmap()
+        resmap['inobs']     = self._cntcube
+        resmap['modcube']   = 'NONE'
+        resmap['expcube']   = 'NONE'
+        resmap['psfcube']   = 'NONE'
+        resmap['edispcube'] = 'NONE'
+        resmap['bkgcube']   = 'NONE'
+        resmap['caldb']     = self._caldb
+        resmap['irf']       = self._irf
+        resmap['inmodel']   = self._model
+        resmap['outmap']    = 'csresmap_py4.fits'
+        resmap['algorithm'] = 'SIGNIF'
+        resmap['logfile']   = 'csresmap_py4.log'
+        resmap['chatter']   = 3
+        resmap['publish']   = True
+
+        # Run csresmap script
+        resmap.execute()
+
+        # Check pull distribution file
+        self._check_result_file('csresmap_py4.fits', nx=200, ny=200)
+
+
         # Return
         return
  
@@ -117,7 +117,7 @@ def test(installed=False):
              test_ctulimit.Test(),
              test_cterror.Test(),
              test_pipelines.Test()]
-
+    
     # Allocate test suite container
     suites = gammalib.GTestSuites('ctools unit testing')