Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents d6e2e852 ca35e885
Showing 7 changed files with 118 additions and 49 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
@@ -27,7 +27,10 @@ srcmodel,  s, a, POINT,POINT,, &quot;Source model type&quot;
 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)"
+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?"
  
@@ -51,6 +51,9 @@ class cssrcdetect(ctools.cscript):
         else:
             self._map = gammalib.GSkyMap()
  
+        # Initialise other members
+        self._map_dirs = []
+
         # Return
         return
  
@@ -73,6 +76,12 @@ class cssrcdetect(ctools.cscript):
         self['fit_pos'].boolean()
         self['fit_shape'].boolean()
  
+        # Query the smoothing parameters
+        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()
@@ -100,20 +109,13 @@ class cssrcdetect(ctools.cscript):
             self._log_header3(gammalib.NORMAL, 'Iteration '+str(i+1))
  
             # Get map moments
-            mean, std = self._map_moments(counts)
+            mean, std = self._map_moments(counts, self['avgrad'].real())
  
-            # Log map moments
-            self._log_value(gammalib.NORMAL, 'Map mean', mean)
-            self._log_value(gammalib.NORMAL, 'Map standard deviation', std)
- 
             # Compute threshold
-            threshold = mean + self['threshold'].real() * std
-
-            # Log map threshold
-            self._log_value(gammalib.NORMAL, 'Map threshold', threshold)
+            sigmap = (counts - mean)/std
  
             # Get maximum value and corresponding sky direction
-            value, pos = self._find_maximum(counts, threshold=threshold)
+            value, pos = self._find_maximum(sigmap)
  
             # If maximum found then log maximum and add model
             if pos is not None:
@@ -124,12 +126,12 @@ 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)
  
                 # Remove maximum from map
-                counts = self._remove_maximum(counts, pos, mean,
+                counts = self._remove_maximum(counts, pos, mean(pos),
                                               radius=self['exclrad'].real())
  
             # ... otherwise log that no maximum was found and break iterations
@@ -145,7 +147,7 @@ class cssrcdetect(ctools.cscript):
         # Return
         return
  
-    def _find_maximum(self, skymap, threshold=0.0):
+    def _find_maximum(self, skymap):
         """
         Find maximum in a sky map
  
@@ -153,8 +155,6 @@ class cssrcdetect(ctools.cscript):
         ----------
         skymap : `~gammalib.GSkyMap()`
             Sky map
-        threshold : float, optional
-            Threshold for maximum value
  
         Returns
         -------
@@ -162,14 +162,14 @@ class cssrcdetect(ctools.cscript):
             Maximum sky map value and corresponding sky direction
         """
         # Initialise maximum pixel value and sky direction
-        value = threshold
+        value = self['threshold'].real()
         pos   = None
  
         # Loop over all pixels and find maximum
         for i in range(skymap.npix()):
             if skymap[i] > value:
                 value = skymap[i]
-                pos   = skymap.inx2dir(i)
+                pos   = self._map_dirs[i]
  
         # Return sky direction of maximum
         return value, pos
@@ -196,17 +196,20 @@ class cssrcdetect(ctools.cscript):
         """
         # Copy skymap
         skymap_copy = skymap.copy()
-        
+
+        # Cache the cosine of the radius
+        cos_radius = math.cos(math.radians(radius))
+
         # Loop over all pixels and find maximum
         for i in range(skymap_copy.npix()):
-            skymap_dir = skymap_copy.inx2dir(i)
-            if skymap_dir.dist_deg(pos) < radius:
+            skymap_dir = self._map_dirs[i]
+            if skymap_dir.cos_dist(pos) > cos_radius:
                 skymap_copy[i] = value
  
         # Return copy of map
         return skymap_copy
  
-    def _map_moments(self, skymap):
+    def _map_moments(self, skymap, radius):
         """
         Determine moments of sky map pixels
  
@@ -214,26 +217,26 @@ class cssrcdetect(ctools.cscript):
         ----------
         skymap : `~gammalib.GSkyMap()`
             Sky map
+        radius : float
+            radius (deg) for pixel consideration
  
         Returns
         -------
-        mean, std : tuple of float
-            Mean and standard deviation of pixel values
+        mean, std : tuple of GSkyMap
+            Mean and standard deviation of pixel values within a given radius
         """
-        # Initialise mean and standard deviation
-        mean = 0.0
-        std  = 0.0
+        # Copy the input skymap
+        mean = gammalib.GSkyMap(skymap)
+        std  = gammalib.GSkyMap(skymap)
+        std *= std
  
-        # Loop over all pixels
-        for i in range(skymap.npix()):
-            mean += skymap[i]
-            std  += skymap[i]*skymap[i]
+        # Convolve by disk to get bin-by-bin mean
+        mean.smooth('DISK', radius)
+        std.smooth('DISK', radius)
  
-        # Compute mean and standard deviation
-        mean /= float(skymap.npix())
-        std  /= float(skymap.npix())
-        std  -= mean * mean
-        std   = math.sqrt(std)
+        # Compute the standard deviation for each pixel
+        std = std - (mean*mean)
+        std = std.sqrt()
  
         # Return mean and standard deviation
         return mean, std
@@ -338,6 +341,30 @@ class cssrcdetect(ctools.cscript):
         return model
  
  
+    def _smooth_skymap(self, skymap):
+        """
+        Smooth the input sky map if a valid kernel was supplied
+
+        Parameters
+        ----------
+        skymap : `~gammalib.GSkyMap()`
+            Sky map
+        """
+        # Log header
+        self._log_header3(gammalib.NORMAL, 'Smoothing Skymap')
+
+        # Make sure the smoothing kernel is not 'NONE'
+        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.')
+
+        return
+
+
     def _load_skymap(self):
         """
         Load sky map
@@ -360,7 +387,18 @@ class cssrcdetect(ctools.cscript):
         fits.close()
  
         # Return
-        return skymap
+        return skymap.extract(0)
+
+
+    def _cache_map_dirs(self, skymap):
+        """
+        Cache map pixel positions to save some computation time
+        """
+        # Setup the skymap directions
+        if not skymap.is_empty():
+            self._map_dirs = [skymap.inx2dir(i) for i in range(skymap.npix())]
+        else:
+            self._map_dirs = []
  
  
     # Public methods
@@ -378,13 +416,18 @@ class cssrcdetect(ctools.cscript):
         # Get parameters
         self._get_parameters()
  
-        # Write header
-        self._log_header1(gammalib.NORMAL, 'Source detection')
-
         # If sky map is empty then load it from input parameter
         if self._map.is_empty():
             self._map = self._load_skymap()
  
+        self._cache_map_dirs(self._map)
+
+        # Smooth the map
+        self._smooth_skymap(self._map)
+
+        # Write header
+        self._log_header1(gammalib.NORMAL, 'Source detection')
+
         # Detect sources
         self._detect_sources(self._map)
  
@@ -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>`_] -
@@ -51,6 +51,17 @@ General parameters
 ``(maxsrcs = 20) [integer]``
     Maximum number of sources that should be detected.
  
+``(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.
+
+``corr_rad [real]``
+    Correlation kernel radius (degrees).
+
+``(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
     detection (degrees).
@@ -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,7 +75,8 @@ class Test(test):
         cmd = cssrcdetect+' inmap="'+self._map+'"'+ \
                           ' outmodel="cssrcdetect_cmd1.xml"'+ \
                           ' outds9file="cssrcdetect_cmd1.reg"'+ \
-                          ' 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
@@ -89,7 +90,8 @@ class Test(test):
         cmd = cssrcdetect+' inmap="skymap_that_does_not_exist"'+ \
                           ' outmodel="cssrcdetect_cmd1.xml"'+ \
                           ' outds9file="cssrcdetect_cmd1.reg"'+ \
-                          ' 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
@@ -129,9 +131,10 @@ class Test(test):
         srcdetect['outds9file'] = 'cssrcdetect_py1.reg'
         srcdetect['srcmodel']   = 'POINT'
         srcdetect['bkgmodel']   = 'NONE'
-        srcdetect['threshold']  = 10.0
+        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
@@ -146,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
@@ -168,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