From 2bb02a07b75400fea7bcbf89d9b192a6a2c70ada Mon Sep 17 00:00:00 2001
From: Jürgen Knödlseder <jurgen.knodlseder@irap.omp.eu>
Date: Tue, 13 Jun 2017 11:14:26 +0200
Subject: [PATCH] Correct standard deviation computation in cssrcdetect and add constructor from skymap and method that returns model container

---
 cscripts/cssrcdetect.py | 61 ++++++++++++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 46 insertions(+), 15 deletions(-)

diff --git a/cscripts/cssrcdetect.py b/cscripts/cssrcdetect.py
index 1cc0ad5..8ebc8a3 100755
--- a/cscripts/cssrcdetect.py
+++ b/cscripts/cssrcdetect.py
@@ -45,6 +45,13 @@ class cssrcdetect(ctools.cscript):
         # Set protected members
         self._models = gammalib.GModels()
 
+        # Initialise sky map from constructor arguments
+        if len(argv) > 0 and isinstance(argv[0], gammalib.GSkyMap):
+            self._map = argv[0]
+            argv      = argv[1:]
+        else:
+            self._map = gammalib.GSkyMap()
+
         # Initialise application by calling the appropriate class constructor
         self._init_cscript(argv)
 
@@ -57,8 +64,9 @@ class cssrcdetect(ctools.cscript):
         """
         Get parameters from parfile
         """
-        # Query input parameters
-        self['inmap'].filename()
+        # Query input parameters if sky map is empty
+        if self._map.is_empty():
+            self['inmap'].filename()
 
         # Query further parameters
         self['srcmodel'].string()
@@ -228,7 +236,7 @@ class cssrcdetect(ctools.cscript):
         # Compute mean and standard deviation
         mean /= float(map.npix())
         std  /= float(map.npix())
-        std  -= mean
+        std  -= mean * mean
         std   = math.sqrt(std)
 
         # Return mean and standard deviation
@@ -334,6 +342,31 @@ class cssrcdetect(ctools.cscript):
         return model
 
 
+    def _load_skymap(self):
+        """
+        Load sky map
+
+        Returns
+        -------
+        map : `~gammalib.GSkyMap()`
+            Sky map
+        """
+        # Get skymap filename
+        inmap = self['inmap'].filename()
+
+        # Open sky map file
+        fits = gammalib.GFits(inmap)
+
+        # Extract primary extension as sky map
+        map = gammalib.GSkyMap(fits.image(0))
+
+        # Close sky map file
+        fits.close()
+
+        # Return
+        return map
+
+
     # Public methods
     def run(self):
         """
@@ -352,20 +385,12 @@ class cssrcdetect(ctools.cscript):
         # Write header
         self._log_header1(gammalib.NORMAL, 'Source detection')
 
-        # Get skymap filename
-        inmap = self['inmap'].filename()
-
-        # Open sky map file
-        fits = gammalib.GFits(inmap)
-
-        # Extract primary extension as counts map
-        counts = gammalib.GSkyMap(fits.image(0))
+        # If sky map is empty then load it from input parameter
+        if self._map.is_empty():
+            self._map = self._load_skymap()
 
         # Detect sources
-        self._detect_sources(counts)
-
-        # Close sky map file
-        fits.close()
+        self._detect_sources(self._map)
 
         # Write detected sources as models into logger
         self._log_models(gammalib.NORMAL, self._models, 'Detected source model')
@@ -433,6 +458,12 @@ class cssrcdetect(ctools.cscript):
         # Return
         return
 
+    def models(self):
+        """
+        Return model container
+        """
+        return self._models
+
     def execute(self):
         """
         Execute the script
--
libgit2 0.21.4