Reduce code duplication

Jürgen Knödlseder
1 parent ff3e3615
Showing 6 changed files with 142 additions and 174 deletions
cscripts/cslightcrv.py
cscripts/csphasecrv.py
cscripts/obsutils.py
pyext/ctool.i
src/support/ctool.cpp
src/support/ctool.hpp
@@ -106,30 +106,10 @@ class cslightcrv(ctools.cscript):
         # Get time boundaries
         self._tbins = self._create_tbounds()
  
-        # Set stacked analysis flag to True if the requested number of
-        # energy bins is positive. Otherwise an unbinned analysis will
-        # be done and the stacked analysis flag will be set to False.
-        if self['enumbins'].integer() > 0:
-            self._stacked = True
-        else:
-            self._stacked = False
-
-        # Make sure that remaining user parameters are queried now. We
-        # do not store the actual parameter values as we do not want
-        # too many instance attributes with enhances the maintenance
-        # costs.
-        self['emin'].real()
-        self['emax'].real()
-        if self._stacked:
-            self['coordsys'].string()
-            self['proj'].string()
-            self['xref'].real()
-            self['yref'].real()
-            self['nxpix'].integer()
-            self['nypix'].integer()
-            self['binsz'].real()
-
-        # Do the same for the hidden parameters, just in case
+        # Set stacked analysis flag and query relevant user parameters
+        self._stacked = self._is_stacked()
+
+        # Query the hidden parameters, just in case
         self['edisp'].boolean()
         self['calc_ulim'].boolean()
         self['calc_ts'].boolean()
@@ -428,57 +408,11 @@ class cslightcrv(ctools.cscript):
         obs : `~gammalib.GObservations`
             Observation container where the first observation is a binned observation
         """
-        # Header
-        if self._logExplicit():
-            self._log.header3('Binning events')
-
-        # Bin events
-        cntcube = ctools.ctbin(obs)
-        cntcube['usepnt']   = False
-        cntcube['ebinalg']  = 'LOG'
-        cntcube['xref']     = self['xref'].real()
-        cntcube['yref']     = self['yref'].real()
-        cntcube['binsz']    = self['binsz'].real()
-        cntcube['nxpix']    = self['nxpix'].integer()
-        cntcube['nypix']    = self['nypix'].integer()
-        cntcube['enumbins'] = self['enumbins'].integer()
-        cntcube['emin']     = self['emin'].real()
-        cntcube['emax']     = self['emax'].real()
-        cntcube['coordsys'] = self['coordsys'].string()
-        cntcube['proj']     = self['proj'].string()
-        cntcube.run()
-
-        # Header
-        if self._logExplicit():
-            self._log.header3('Creating stacked response')
-
-        # Get stacked response
-        response = obsutils.get_stacked_response(obs,
-                                                 self['xref'].real(),
-                                                 self['yref'].real(),
-                                                 binsz=self['binsz'].real(),
-                                                 nxpix=self['nxpix'].integer(),
-                                                 nypix=self['nypix'].integer(),
-                                                 emin=self['emin'].real(),
-                                                 emax=self['emax'].real(),
-                                                 enumbins=self['enumbins'].integer(),
-                                                 edisp=self['edisp'].boolean(),
-                                                 coordsys=self['coordsys'].string(),
-                                                 proj=self['proj'].string())
-
-        # Retrieve a new oberservation container
-        new_obs = cntcube.obs().copy()
-
-        # Set stacked response
-        if self['edisp'].boolean():
-            new_obs[0].response(response['expcube'], response['psfcube'],
-                                response['edispcube'], response['bkgcube'])
-        else:
-            new_obs[0].response(response['expcube'], response['psfcube'],
-                                response['bkgcube'])
-
-        # Get new models
-        models = response['models']
+        # Get new observation container with binned observation
+        new_obs = obsutils.get_stacked_obs(self, obs)
+        
+        # Extract models
+        models = new_obs.models()
  
         # Fix background models if required
         if self['fix_bkg'].boolean():
@@ -487,7 +421,7 @@ class cslightcrv(ctools.cscript):
                     for par in model:
                         par.fix()
  
-        # Set models for new oberservation container     
+        # Put back models
         new_obs.models(models)
  
         # Return new oberservation container
@@ -91,33 +91,13 @@ class csphasecrv(ctools.cscript):
         # Get source name
         self._srcname = self['srcname'].string()
  
-        # Get time boundaries
+        # Get phase boundaries
         self._phbins = self._create_tbounds()
  
-        # Set stacked analysis flag to True if the requested number of
-        # energy bins is positive. Otherwise an unbinned analysis will
-        # be done and the stacked analysis flag will be set to False.
-        if self['enumbins'].integer() > 0:
-            self._stacked = True
-        else:
-            self._stacked = False
-
-        # Make sure that remaining user parameters are queried now. We
-        # do not store the actual parameter values as we do not want
-        # too many instance attributes with enhances the maintenance
-        # costs.
-        self['emin'].real()
-        self['emax'].real()
-        if self._stacked:
-            self['coordsys'].string()
-            self['proj'].string()
-            self['xref'].real()
-            self['yref'].real()
-            self['nxpix'].integer()
-            self['nypix'].integer()
-            self['binsz'].real()
-
-        # Do the same for the hidden parameters, just in case
+        # Set stacked analysis flag and query relevant user parameters
+        self._stacked = self._is_stacked()
+
+        # Query the hidden parameters, just in case
         self['edisp'].boolean()
  
         # Read ahead output parameters
@@ -181,78 +161,6 @@ class csphasecrv(ctools.cscript):
         # Return Good Time Intervals
         return phbins
  
-    def _bin_observation(self, obs):
-        """
-        Bin an observation if a binned analysis was requested
-
-        Parameters
-        ----------
-        obs : `~gammalib.GObservations`
-            Observation container
-
-        Returns
-        -------
-        obs : `~gammalib.GObservations`
-            Observation container where the first observation is a binned observation
-        """
-        # Header
-        if self._logExplicit():
-            self._log.header3('Binning events')
-
-        # Bin events
-        cntcube = ctools.ctbin(obs)
-        cntcube['usepnt']   = False
-        cntcube['ebinalg']  = 'LOG'
-        cntcube['xref']     = self['xref'].real()
-        cntcube['yref']     = self['yref'].real()
-        cntcube['binsz']    = self['binsz'].real()
-        cntcube['nxpix']    = self['nxpix'].integer()
-        cntcube['nypix']    = self['nypix'].integer()
-        cntcube['enumbins'] = self['enumbins'].integer()
-        cntcube['emin']     = self['emin'].real()
-        cntcube['emax']     = self['emax'].real()
-        cntcube['coordsys'] = self['coordsys'].string()
-        cntcube['proj']     = self['proj'].string()
-        cntcube.run()
-
-        # Header
-        if self._logExplicit():
-            self._log.header3('Creating stacked response')
-
-        # Get stacked response
-        response = obsutils.get_stacked_response(obs,
-                                                 self['xref'].real(),
-                                                 self['yref'].real(),
-                                                 binsz=self['binsz'].real(),
-                                                 nxpix=self['nxpix'].integer(),
-                                                 nypix=self['nypix'].integer(),
-                                                 emin=self['emin'].real(),
-                                                 emax=self['emax'].real(),
-                                                 enumbins=self['enumbins'].integer(),
-                                                 edisp=self['edisp'].boolean(),
-                                                 coordsys=self['coordsys'].string(),
-                                                 proj=self['proj'].string())
-
-        # Retrieve a new oberservation container
-        new_obs = cntcube.obs().copy()
-
-        # Set stacked response
-        if self['edisp'].boolean():
-            new_obs[0].response(response['expcube'], response['psfcube'],
-                                response['edispcube'], response['bkgcube'])
-        else:
-            new_obs[0].response(response['expcube'], response['psfcube'],
-                                response['bkgcube'])
-
-        # Get new models
-        models = response['models']
-
-        # Set models for new oberservation container     
-        new_obs.models(models)
-
-        # Return new oberservation container
-        return new_obs
-
     def _get_free_par_names(self):
         """
         Return list of free parameter names
@@ -469,7 +377,8 @@ class csphasecrv(ctools.cscript):
             # and compute the stacked response functions and setup
             # an observation container with a single stacked observation.
             if self._stacked:
-                obs = self._bin_observation(select.obs())
+                #obs = self._bin_observation(select.obs())
+                obs = obsutils.get_stacked_obs(self, select.obs())
  
             # Header
             if self._logExplicit():
@@ -639,3 +639,82 @@ def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200,
  
     # Return response cubes
     return response
+
+
+# ================================= #
+# Get stacked observation container #
+# ================================= #
+def get_stacked_obs(cls, obs):
+    """
+    Bin an observation and return an observation container with a single
+    binned observation
+
+    Parameters
+    ----------
+    cls : `~ctools.cscript`
+        cscript class
+    obs : `~gammalib.GObservations`
+        Observation container
+
+    Returns
+    -------
+    obs : `~gammalib.GObservations`
+        Observation container where the first observation is a binned observation
+    """
+    # Write header
+    if cls._logExplicit():
+        cls._log.header3('Binning events')
+    
+    # Bin events
+    cntcube = ctools.ctbin(obs)
+    cntcube['usepnt']   = False
+    cntcube['ebinalg']  = 'LOG'
+    cntcube['xref']     = cls['xref'].real()
+    cntcube['yref']     = cls['yref'].real()
+    cntcube['binsz']    = cls['binsz'].real()
+    cntcube['nxpix']    = cls['nxpix'].integer()
+    cntcube['nypix']    = cls['nypix'].integer()
+    cntcube['enumbins'] = cls['enumbins'].integer()
+    cntcube['emin']     = cls['emin'].real()
+    cntcube['emax']     = cls['emax'].real()
+    cntcube['coordsys'] = cls['coordsys'].string()
+    cntcube['proj']     = cls['proj'].string()
+    cntcube.run()
+
+    # Write header
+    if cls._logExplicit():
+        cls._log.header3('Creating stacked response')
+
+    # Get stacked response
+    response = get_stacked_response(obs,
+                                    cls['xref'].real(),
+                                    cls['yref'].real(),
+                                    binsz=cls['binsz'].real(),
+                                    nxpix=cls['nxpix'].integer(),
+                                    nypix=cls['nypix'].integer(),
+                                    emin=cls['emin'].real(),
+                                    emax=cls['emax'].real(),
+                                    enumbins=cls['enumbins'].integer(),
+                                    edisp=cls['edisp'].boolean(),
+                                    coordsys=cls['coordsys'].string(),
+                                    proj=cls['proj'].string())
+
+    # Retrieve a new oberservation container
+    new_obs = cntcube.obs().copy()
+
+    # Set stacked response
+    if cls['edisp'].boolean():
+        new_obs[0].response(response['expcube'], response['psfcube'],
+                            response['edispcube'], response['bkgcube'])
+    else:
+        new_obs[0].response(response['expcube'], response['psfcube'],
+                            response['bkgcube'])
+
+    # Get new models
+    models = response['models']
+
+    # Set models for new oberservation container
+    new_obs.models(models)
+
+    # Return new oberservation container
+    return new_obs
@@ -99,6 +99,7 @@ public:
     %rename(_get_pointing)             get_pointing;
     %rename(_get_skydir)               get_skydir;
     %rename(_set_outfile_name)         set_outfile_name;
+    %rename(_is_stacked)               is_stacked;
     %rename(_set_response)             set_response;
     %rename(_set_edisp)                set_edisp;
     %rename(_restore_edisp)            restore_edisp;
@@ -143,6 +144,9 @@ public:
     GSkyDir           get_skydir(void);
     std::string       set_outfile_name(const std::string& filename);
  
+    // Protected methods that query user parameters
+    bool              is_stacked(void);
+
     // Protected support methods
     void              set_response(GObservations& obs);
     std::vector<bool> set_edisp(GObservations& obs, const bool& edisp) const;
@@ -1156,6 +1156,45 @@ std::string ctool::set_outfile_name(const std::string&amp; filename)
  
  
 /***********************************************************************//**
+ * @brief Query user parameters for stacked analysis
+ *
+ * @return True if a stacked analysis should be done.
+ *
+ * Queries the user parameters "emin", "emax" and "enumbins", and if
+ * "enumbins" is positive, signal that a stacked analysis is requested. In
+ * that case, also the "coordsys", "proj", "xref", "yref", "nxpix", "nypix"
+ * and "binsz" parameters are queried.
+ *
+ * Using this method assures that the parameters are always queried in the
+ * same order.
+ ***************************************************************************/
+bool ctool::is_stacked(void)
+{
+    // First query the minimum and maximum energies
+    (*this)["emin"].real();
+    (*this)["emax"].real();
+
+    // Now query the number of energy bins and set the stacked flag
+    bool stacked = ((*this)["enumbins"].integer() > 0) ? true : false;
+
+    // If a stacked analysis is requested then query the spatial definition
+    // of the cube
+    if (stacked) {
+        (*this)["coordsys"].string();
+        (*this)["proj"].string();
+        (*this)["xref"].real();
+        (*this)["yref"].real();
+        (*this)["nxpix"].integer();
+        (*this)["nypix"].integer();
+        (*this)["binsz"].real();
+    }
+
+    // Return stacked flag
+    return stacked;
+}
+
+
+/***********************************************************************//**
  * @brief Set response for all CTA observations in container
  *
  * @param[in,out] obs Observation container
@@ -102,6 +102,9 @@ protected:
     GSkyDir           get_skydir(void);
     std::string       set_outfile_name(const std::string& filename);
  
+    // Protected methods that query user parameters
+    bool              is_stacked(void);
+
     // Protected methods for logging
     void              log_observations(const GChatter&      chatter,
                                        const GObservations& obs,