From 0ffbb73c6488a74dc3316a2daaeacc676f5e749c Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 10 Aug 2017 23:43:18 +0200 Subject: [PATCH] Reduce code duplication --- cscripts/cslightcrv.py | 86 ++++++++++---------------------------------------------------------------------------- cscripts/csphasecrv.py | 105 +++++++-------------------------------------------------------------------------------------------------- cscripts/obsutils.py | 79 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/ctool.i | 4 ++++ src/support/ctool.cpp | 39 +++++++++++++++++++++++++++++++++++++++ src/support/ctool.hpp | 3 +++ 6 files changed, 142 insertions(+), 174 deletions(-) diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index 13cb324..e5ff788 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -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 diff --git a/cscripts/csphasecrv.py b/cscripts/csphasecrv.py index 695d57c..7fc0b75 100755 --- a/cscripts/csphasecrv.py +++ b/cscripts/csphasecrv.py @@ -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(): diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index b66bfb5..4c17d3d 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -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 diff --git a/pyext/ctool.i b/pyext/ctool.i index 3f2e45f..3eadb0e 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -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 set_edisp(GObservations& obs, const bool& edisp) const; diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index ff17d96..3f6dc4b 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1156,6 +1156,45 @@ std::string ctool::set_outfile_name(const std::string& 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 diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index bf9cc5e..26a6aea 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -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, -- libgit2 0.21.4