diff --git a/cscripts/csresmap.par b/cscripts/csresmap.par index 67564f2..6416350 100644 --- a/cscripts/csresmap.par +++ b/cscripts/csresmap.par @@ -43,7 +43,7 @@ yref, r, a, 22.01,-90,90, "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 diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index ff1e7cd..28db110 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -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') diff --git a/test/test_csresmap.py b/test/test_csresmap.py index 392e372..f828818 100755 --- a/test/test_csresmap.py +++ b/test/test_csresmap.py @@ -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 diff --git a/test/test_python_ctools.py b/test/test_python_ctools.py index 9f546da..782fafc 100755 --- a/test/test_python_ctools.py +++ b/test/test_python_ctools.py @@ -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')