Remove code duplication in csresmap and csresspec

Jürgen Knödlseder
1 parent 83ee3aec
Showing 3 changed files with 62 additions and 49 deletions
cscripts/csresmap.py
cscripts/csresspec.py
cscripts/obsutils.py
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a residual map.
 #
-# Copyright (C) 2014-2018 Michael Mayer
+# Copyright (C) 2014-2019 Michael Mayer
 #
 # 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
@@ -165,29 +165,8 @@ class csresmap(ctools.csobservation):
                 # Write header
                 self._log_header1(gammalib.TERSE, 'Generate binned map (ctbin)')
-                # Create countsmap
-                binning = ctools.ctbin(self.obs())
-                binning['xref']     = self['xref'].real()
-                binning['yref']     = self['yref'].real()
-                binning['proj']     = self['proj'].string()
-                binning['coordsys'] = self['coordsys'].string()
-                binning['ebinalg']  = self['ebinalg'].string()
-                binning['nxpix']    = self['nxpix'].integer()
-                binning['nypix']    = self['nypix'].integer()
-                binning['binsz']    = self['binsz'].real()
-                if self['ebinalg'].string() == 'FILE':
-                    binning['ebinfile'] = self['ebinfile'].filename().file()
-                else:
-                    binning['enumbins'] = self['enumbins'].integer()
-                    binning['emin']     = self['emin'].real()
-                    binning['emax']     = self['emax'].real()
-                binning['chatter']  = self['chatter'].integer()
-                binning['clobber']  = self['clobber'].boolean()
-                binning['debug']    = self['debug'].boolean()
-                binning.run()
-
-                # Retrieve counts cube
-                cta_counts_cube = binning.cube()
+                # Create counts cube
+                cta_counts_cube = obsutils.create_counts_cube(self, self.obs())
             # Assign GCTAEventCube to skymap
             countmap = cta_counts_cube.counts()
@@ -18,10 +18,10 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import sys
 import gammalib
 import ctools
 from cscripts import obsutils
-import sys
 # =============== #
@@ -672,35 +672,18 @@ class csresspec(ctools.csobservation):
         # Log processing header
         self._log_header1(gammalib.TERSE, 'Processing Observations')
-        #If observations are unbinned and we stack
+        # If observations are unbinned and we stack
         if self._stack and self.obs()[0].classname() == 'GCTAObservation':
+            # Log counts cube computation
             msg = 'Computing count cube from multiple unbinned observations'
             self._log_string(gammalib.NORMAL, msg)
-            # Build count cube
-            binning = ctools.ctbin(self.obs())
-            binning['xref']     = self['xref'].real()
-            binning['yref']     = self['yref'].real()
-            binning['proj']     = self['proj'].string()
-            binning['coordsys'] = self['coordsys'].string()
-            binning['ebinalg']  = self['ebinalg'].string()
-            binning['nxpix']    = self['nxpix'].integer()
-            binning['nypix']    = self['nypix'].integer()
-            binning['binsz']    = self['binsz'].real()
-            if self['ebinalg'].string() == 'FILE':
-                binning['ebinfile'] = self['ebinfile'].filename().file()
-            else:
-                binning['enumbins'] = self['enumbins'].integer()
-                binning['emin']     = self['emin'].real()
-                binning['emax']     = self['emax'].real()
-            binning['chatter'] = self['chatter'].integer()
-            binning['clobber'] = self['clobber'].boolean()
-            binning['debug']   = self['debug'].boolean()
-            binning.run()
-
-            #compute residuals using cube
-            table = self._residuals_3D(self.obs(),'',binning.cube())
+            # Create counts cube
+            cta_counts_cube = obsutils.create_counts_cube(self, self.obs())
+
+            # Compute residuals using cube
+            table = self._residuals_3D(self.obs(), '', cta_counts_cube)
             # Append results table to output file
             self._fits.append(table)
@@ -937,3 +937,54 @@ def residuals(cls, counts, model):
     # Return
     return residuals
+
+
+# ================= #
+# Create counts map #
+# ================= #
+def create_counts_cube(cls, obs):
+    """
+    Create counts cube from observations
+
+    Parameters
+    ----------
+    cls : `~ctools.cscript`
+        cscript class
+    obs : `~gammalib.GObservations`
+        Observation container
+
+    Returns
+    -------
+    cntcube : `~gammalib.GCTAEventCube`
+        Counts cube
+    """
+    # Initialise ctbin
+    ctbin = ctools.ctbin(obs)
+    
+    # Set parameters
+    ctbin['xref']     = cls['xref'].real()
+    ctbin['yref']     = cls['yref'].real()
+    ctbin['proj']     = cls['proj'].string()
+    ctbin['coordsys'] = cls['coordsys'].string()
+    ctbin['ebinalg']  = cls['ebinalg'].string()
+    ctbin['nxpix']    = cls['nxpix'].integer()
+    ctbin['nypix']    = cls['nypix'].integer()
+    ctbin['binsz']    = cls['binsz'].real()
+    if cls['ebinalg'].string() == 'FILE':
+        ctbin['ebinfile'] = cls['ebinfile'].filename().file()
+    else:
+        ctbin['enumbins'] = cls['enumbins'].integer()
+        ctbin['emin']     = cls['emin'].real()
+        ctbin['emax']     = cls['emax'].real()
+    ctbin['chatter']  = cls['chatter'].integer()
+    ctbin['clobber']  = cls['clobber'].boolean()
+    ctbin['debug']    = cls['debug'].boolean()
+
+    # Run ctbin
+    ctbin.run()
+
+    # Retrieve counts cube
+    cntcube = ctbin.cube().copy()
+
+    # Return counts cube
+    return cntcube