Change names of cssrcdetect robustness parameters (#2498)

Jürgen Knödlseder
1 parent 77884dc8
Showing 7 changed files with 55 additions and 52 deletions
ChangeLog
NEWS
cscripts/cssrcdetect.par
cscripts/cssrcdetect.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/cssrcdetect.rst
test/test_cssrcdetect.py
-2018-09-07
+2018-09-10
         * Version 1.6.0 released
         ========================
@@ -9,6 +9,7 @@
         Use full RMF true energy range for On/Off analysis in csspec (#2656)
         Defer import of ROOT module in csroot2caldb to run method (#2652)
         Transparently update parameter files (#2513)
+        Make cssrcdetect more robust (#2498)
         Change formula for number of true energy bins in csphagen (#2490)
         Remove ctools::is_valid_filename() method
         Handle zero model for SUBDIV method in show_response.py (#2456)
 New Features and Important Changes in ctools 1.6.0
-7 September 2018
+10 September 2018
 Introduction
@@ -335,7 +335,12 @@ certain amount (#2233).
 cssrcdetect - Detects sources in sky map
 ----------------------------------------
-None
+The script was made more robust by optionally convolving the input sky map
+with a smoothing kernel with a specified correlation radius (options are a
+"DISK" or a "GAUSSIAN" kernel, default is now "DISK"). In addition, the mean
+and standard deviation of the pixel values are now computed locally over a
+given radius, controlled through the "avgrad" parameter. This accomodates for
+counts variations in the field of view (#2498).
 cstsdist - Computes Test Statistic distribution
@@ -23,16 +23,16 @@ outds9file, f, a, ds9.reg,,, &quot;Output DS9 region file&quot;
 #
 # Script parameters
 #==================
-srcmodel,     s, a, POINT,POINT,, "Source model type"
-bkgmodel,     s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type"
-threshold,    r, a, 5.0,,, "Detection threshold (Gaussian sigma)"
-maxsrcs,      i, h, 20,1,1000, "Maximum number of sources"
-exclrad,      r, h, 0.2,,, "Exclusion radius (deg)"
-fit_pos,      b, h, yes,,, "Fit source positions?"
-fit_shape,    b, h, yes,,, "Fit source shapes?"
-momentradius, r, a, 200,,, "Radius (deg) for computing map mean and std-deviation values"
-smoothkernel, s, a, NONE,NONE|DISK|GAUSSIAN,,"Smoothing kernel to use for smoothing the input sky map"
-smoothparam,  r, a, 0.05,,,"Smoothing kernel parameters (sigma for GAUS, radius for DISK)"
+srcmodel,  s, a, POINT,POINT,, "Source model type"
+bkgmodel,  s, a, NONE,NONE|IRF|AEFF|CUBE,, "Background model type"
+threshold, r, a, 5.0,,, "Detection threshold (Gaussian sigma)"
+maxsrcs,   i, h, 20,1,1000, "Maximum number of sources"
+avgrad,    r, h, 1.0,,, "Averaging radius (deg) for significance computation"
+corr_rad,  r, a, 0.1,,, "Correlation kernel radius (deg)"
+corr_kern, s, h, DISK,NONE|DISK|GAUSSIAN,, "Correlation kernel to apply to input sky map"
+exclrad,   r, h, 0.3,,, "Exclusion radius (deg)"
+fit_pos,   b, h, yes,,, "Fit source positions?"
+fit_shape, b, h, yes,,, "Fit source shapes?"
 #
 # Standard parameters
@@ -77,15 +77,15 @@ class cssrcdetect(ctools.cscript):
         self['fit_shape'].boolean()
         # Query the smoothing parameters
-        self['smoothkernel'].string()
-        if self['smoothkernel'].string().upper() != 'NONE':
-            self['smoothparam'].real()
+        self['avgrad'].real()
+        self['corr_kern'].string()
+        if self['corr_kern'].string().upper() != 'NONE':
+            self['corr_rad'].real()
         # Query ahead output model filename
         if self._read_ahead():
             self['outmodel'].filename()
             self['outds9file'].filename()
-            self['momentradius'].real()
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
@@ -109,11 +109,11 @@ class cssrcdetect(ctools.cscript):
             self._log_header3(gammalib.NORMAL, 'Iteration '+str(i+1))
             # Get map moments
-            mean, std = self._map_moments(counts, self['momentradius'].real())
- 
+            mean, std = self._map_moments(counts, self['avgrad'].real())
+
             # Compute threshold
             sigmap = (counts - mean)/std
-            
+
             # Get maximum value and corresponding sky direction
             value, pos = self._find_maximum(sigmap)
@@ -126,7 +126,7 @@ class cssrcdetect(ctools.cscript):
                 # Log maximum
                 self._log_value(gammalib.NORMAL, 'Map maximum', str(value))
                 self._log_value(gammalib.NORMAL, name+' position', str(pos))
-            
+
                 # Add model
                 self._add_model(pos, name)
@@ -196,7 +196,7 @@ class cssrcdetect(ctools.cscript):
         """
         # Copy skymap
         skymap_copy = skymap.copy()
-        
+
         # Cache the cosine of the radius
         cos_radius = math.cos(math.radians(radius))
@@ -354,16 +354,13 @@ class cssrcdetect(ctools.cscript):
         self._log_header3(gammalib.NORMAL, 'Smoothing Skymap')
         # Make sure the smoothing kernel is not 'NONE'
-        if self['smoothkernel'].string().upper() != 'NONE':
-            self._log_value(gammalib.NORMAL, 'Kernel',
-                            self['smoothkernel'].string())
-            self._log_value(gammalib.NORMAL, 'Parameter',
-                            self['smoothparam'].real())
-            skymap.smooth(self['smoothkernel'].string(),
-                          self['smoothparam'].real())
+        if self['corr_kern'].string().upper() != 'NONE':
+            self._log_value(gammalib.NORMAL, 'Kernel', self['corr_kern'].string())
+            self._log_value(gammalib.NORMAL, 'Parameter', self['corr_rad'].real())
+            skymap.smooth(self['corr_kern'].string(), self['corr_rad'].real())
         else:
             self._log_string(gammalib.NORMAL, 
-                        'Smoothing kernel is NONE, smoothing will be ignored.')
+                'Smoothing kernel is "NONE", smoothing will be ignored.')
         return
@@ -43,6 +43,8 @@ Improvements
   Use full ``RMF`` true energy range for On/Off analysis in :ref:`csspec`
 * [`2513 <https://cta-redmine.irap.omp.eu/issues/2513>`_] -
   Transparently update parameter files
+* [`2498 <https://cta-redmine.irap.omp.eu/issues/2498>`_] -
+  Make :ref:`cssrcdetect` more robust
 * [`2490 <https://cta-redmine.irap.omp.eu/issues/2490>`_] -
   Change formula for number of true energy bins in :ref:`csphagen`
 * [`2421 <https://cta-redmine.irap.omp.eu/issues/2421>`_] -
@@ -48,18 +48,19 @@ General parameters
     Detection threshold (Gaussian sigma). Only sources above this detection
     threshold will be added to the model definiton XML file.
-``momentradius [real]``
-    Radius (deg) for computing map mean and std-deviation values.
+``(maxsrcs = 20) [integer]``
+    Maximum number of sources that should be detected.
-``smoothkernel <NONE|DISK|GAUSSIAN> [string]``
-    Smoothing kernel to use for smoothing the input sky map.
+``(avgrad = 1.0) [real]``
+    Averaging radius for significance computation (degrees). The script computes
+    the mean pixel value and standard deviation value over this radius and
+    uses the result to estimate the excess significance.
-``smoothparam [real]``
-    Smoothing kernel parameters (sigma for GAUS, radius for DISK, both in 
-    degrees).
+``corr_rad [real]``
+    Correlation kernel radius (degrees).
-``(maxsrcs = 20) [integer]``
-    Maximum number of sources that should be detected.
+``(corr_kern = DISK) <NONE|DISK|GAUSSIAN> [string]``
+    Correlation kernel to use for smoothing the input sky map.
 ``(exclrad = 0.2) [real]``
     Radius around a detected source that is excluded from further source
@@ -33,7 +33,7 @@ class Test(test):
     This test class makes unit tests for the cssrcdetect script by using it
     from the command line and from Python.
     """
-    
+
     # Constructor
     def __init__(self):
         """
@@ -75,9 +75,8 @@ class Test(test):
         cmd = cssrcdetect+' inmap="'+self._map+'"'+ \
                           ' outmodel="cssrcdetect_cmd1.xml"'+ \
                           ' outds9file="cssrcdetect_cmd1.reg"'+ \
-                          ' momentradius=200.0'+ \
-                          ' smoothkernel=NONE smoothparam=0.1'+ \
-                          ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \
+                          ' corr_rad=0.1'+ \
+                          ' srcmodel="POINT" bkgmodel="NONE" threshold=8.0'+ \
                           ' logfile="cssrcdetect_cmd1.log" chatter=1'
         # Check if execution was successful
@@ -91,9 +90,8 @@ class Test(test):
         cmd = cssrcdetect+' inmap="skymap_that_does_not_exist"'+ \
                           ' outmodel="cssrcdetect_cmd1.xml"'+ \
                           ' outds9file="cssrcdetect_cmd1.reg"'+ \
-                          ' momentradius=200.0'+ \
-                          ' smoothkernel=NONE smoothparam=0.1'+ \
-                          ' srcmodel="POINT" bkgmodel="NONE" threshold=10.0'+ \
+                          ' corr_rad=0.1'+ \
+                          ' srcmodel="POINT" bkgmodel="NONE" threshold=8.0'+ \
                           ' logfile="cssrcdetect_cmd2.log" debug=yes chatter=2'
         # Check if execution failed
@@ -133,12 +131,10 @@ class Test(test):
         srcdetect['outds9file'] = 'cssrcdetect_py1.reg'
         srcdetect['srcmodel']   = 'POINT'
         srcdetect['bkgmodel']   = 'NONE'
-        srcdetect['threshold']  = 10.0
-        srcdetect['momentradius'] = 200.0
-        srcdetect['smoothkernel'] = "NONE"
-        srcdetect['smoothparam']  = 0.1
+        srcdetect['threshold']  = 8.0
+        srcdetect['corr_rad']   = 0.1
         srcdetect['logfile']    = 'cssrcdetect_py1.log'
-        srcdetect['chatter']    = 2
+        srcdetect['chatter']    = 4
         # Run cssrcdetect script and save models
         srcdetect.logFileOpen()   # Make sure we get a log file
@@ -153,6 +149,7 @@ class Test(test):
         srcdetect['outmodel']   = 'cssrcdetect_py2.xml'
         srcdetect['outds9file'] = 'cssrcdetect_py2.reg'
         srcdetect['bkgmodel']   = 'IRF'
+        srcdetect['threshold']  = 8.0
         srcdetect['logfile']    = 'cssrcdetect_py2.log'
         srcdetect['chatter']    = 3
         srcdetect.logFileOpen()  # Needed to get a new log file
@@ -175,6 +172,6 @@ class Test(test):
         # Check number of models
         self.test_value(models.size(), number, 'Check for '+str(number)+
                         ' models in XML file')
-        
+
         # Return
         return