Merge branch '1867-csresmap-significance-map' into integration

Jürgen Knödlseder
2 parents 2d5f1b30 4d08028c
Showing 4 changed files with 80 additions and 5 deletions
cscripts/csresmap.par
cscripts/csresmap.py
doc/source/users/reference_manual/csresmap.rst
test/test_csresmap.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|SIGNIFICANCE,, "Residual map computation algorithm"
 #
 # Standard parameters
@@ -239,6 +239,26 @@ class csresmap(ctools.cscript):
             self._resmap -= modelmap
             self._resmap /= modelmap.sqrt()
+        # Calculate significance from Li&Ma derivation
+        elif algorithm == 'SIGNIFICANCE':
+            signmap = (self._resmap - modelmap).sign()
+            logmap  = self._resmap/modelmap
+            # Masking pixels with zero counts. Put a small value >0 for log
+            # calculation and restore the zero 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')
@@ -94,9 +94,17 @@ General parameters
 ``binsz [real]``
     Cube bin size (in degrees/pixel).
-``(algorithm = SUBDIV) <SUB|SUBDIV|SUBDIVSQRT> [string]``
+``(algorithm = SUBDIV) <SUB|SUBDIV|SUBDIVSQRT|SIGNIF> [string]``
     Algorithm used to generate the residual map.
- 	 	 
+
+    SUB:         DATA - MODEL
+
+    SUBDIV:     (DATA - MODEL)/MODEL
+
+    SUBDIVSQRT: (DATA - MODEL)/sqrt(MODEL)
+
+    SIGNIF:     sign(DATA-MODEL) * sqrt( 2*( DATA*log(DATA/MODEL) + MODEL - DATA ) )
+     	 	 
 Standard parameters
 -------------------
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csresmap script.
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 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
@@ -33,7 +33,7 @@ class Test(test):
     This test class makes unit tests for the csresmap script by using it
     from the command line and from Python.
     """
-    
+
     # Constructor
     def __init__(self):
         """
@@ -165,6 +165,53 @@ 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 SIGNIFICANCE 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'] = 'SIGNIFICANCE'
+        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