diff --git a/cscripts/csresmap.par b/cscripts/csresmap.par index ed0c42c..0841d94 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|SIGNIFICANCE,, "Residual map computation algorithm" # # Standard parameters diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index ae53697..a72ca50 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -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') diff --git a/doc/source/users/reference_manual/csresmap.rst b/doc/source/users/reference_manual/csresmap.rst index bdbed8e..2ed0619 100644 --- a/doc/source/users/reference_manual/csresmap.rst +++ b/doc/source/users/reference_manual/csresmap.rst @@ -94,9 +94,17 @@ General parameters ``binsz [real]`` Cube bin size (in degrees/pixel). -``(algorithm = SUBDIV) [string]`` +``(algorithm = SUBDIV) [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 ------------------- diff --git a/test/test_csresmap.py b/test/test_csresmap.py index 392e372..4f6e93b 100755 --- a/test/test_csresmap.py +++ b/test/test_csresmap.py @@ -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