Integrate code

Jürgen Knödlseder
1 parent 14632f27
Showing 9 changed files with 364 additions and 262 deletions
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/csresspec.py
cscripts/obsutils.py
examples/Makefile.am
examples/show_residuals.py
pyext/ctcubemask.i
src/ctcubemask/ctcubemask.hpp
test/test_csresspec.py
@@ -37,7 +37,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
                 $(srcdir)/csphasecrv.py \
                 $(srcdir)/cspull.py \
                 $(srcdir)/csresmap.py \
-		$(srcdir)/csresspec.py \
+                $(srcdir)/csresspec.py \
                 $(srcdir)/csroot2caldb.py  \
                 $(srcdir)/cssens.py \
                 $(srcdir)/csspec.py \
@@ -72,7 +72,7 @@ parfiles = $(srcdir)/cscaldb.par \
            $(srcdir)/csphasecrv.par \
            $(srcdir)/cspull.par \
            $(srcdir)/csresmap.par \
-	   $(srcdir)/csresspec.par \
+           $(srcdir)/csresspec.par \
            $(srcdir)/csroot2caldb.par \
            $(srcdir)/cssens.par \
            $(srcdir)/csspec.par \
@@ -104,7 +104,7 @@ linkscripts = cscaldb \
               csphasecrv \
               cspull \
               csresmap \
-	      csresspec \
+              csresspec \
               csroot2caldb \
               cssens \
               csspec \
 # ==========================================================================
 # cscripts Python module
 #
-# Copyright (C) 2015-2017 Juergen Knoedlseder
+# Copyright (C) 2015-2018 Juergen Knoedlseder
 #
 # 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
@@ -40,6 +40,7 @@ __all__ = [
     "csphasecrv",
     "cspull",
     "csresmap",
+    "csresspec",
     "csroot2caldb",
     "cssens",
     "csspec",
@@ -77,6 +78,7 @@ from cscripts.csphagen      import csphagen
 from cscripts.csphasecrv    import csphasecrv
 from cscripts.cspull        import cspull
 from cscripts.csresmap      import csresmap
+from cscripts.csresspec     import csresspec
 from cscripts.csroot2caldb  import csroot2caldb
 from cscripts.cssens        import cssens
 from cscripts.csspec        import csspec
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a residual spectrum.
 #
-# Copyright (C) 2017- Luigi Tibaldo
+# Copyright (C) 2017-2018 Luigi Tibaldo
 #
 # 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
@@ -29,7 +29,7 @@ import sys
 # =============== #
 class csresspec(ctools.csobservation):
     """
-    Generates a residual map
+    Generates a residual spectrum
     """
  
     # Constructor
@@ -43,9 +43,9 @@ class csresspec(ctools.csobservation):
  
         # Initialise class members
         self._use_maps = False
-        self._stack = False
-        self._mask = False
-        self._fits = gammalib.GFits()
+        self._stack    = False
+        self._mask     = False
+        self._fits     = gammalib.GFits()
  
         # Return
         return
@@ -65,8 +65,8 @@ class csresspec(ctools.csobservation):
         # Collect number of unbinned, binned and On/Off observations in
         # observation container
         n_unbinned = 0
-        n_binned = 0
-        n_onoff = 0
+        n_binned   = 0
+        n_onoff    = 0
         for obs in self.obs():
             if obs.classname() == 'GCTAObservation':
                 if obs.eventtype() == 'CountsCube':
@@ -75,21 +75,10 @@ class csresspec(ctools.csobservation):
                     n_unbinned += 1
             elif obs.classname() == 'GCTAOnOffObservation':
                 n_onoff += 1
-        n_cta = n_unbinned + n_binned + n_onoff
+        n_cta   = n_unbinned + n_binned + n_onoff
         n_other = self.obs().size() - n_cta
  
-        # Log census of input observations
-        self._log_value(gammalib.NORMAL, 'Unbinned CTA observations',
-                        n_unbinned)
-        self._log_value(gammalib.NORMAL, 'Binned CTA observations', n_binned)
-        self._log_value(gammalib.NORMAL, 'On/off CTA observations', n_onoff)
-        self._log_value(gammalib.NORMAL, 'Other observations', n_other)
-        if n_other > 0:
-            msg = 'WARNING: Only CTA observation can be handled, all non-CTA ' \
-                  + 'observations will be ignored.\n'
-            self._log_string(gammalib.TERSE, msg)
-
-        # Query wether to compute model for individual components
+        # Query whether to compute model for individual components
         components = self['components'].boolean()
  
         # If there is only one binned observation and no model for individual
@@ -109,18 +98,13 @@ class csresspec(ctools.csobservation):
                 self['emin'].real()
                 self['emax'].real()
                 self['enumbins'].integer()
-            msg = 'User defined energy binning will be used for %s unbinned observations.' % n_unbinned
-            self._log_string(gammalib.TERSE, msg)
             if n_cta > n_unbinned:
                 n_notunbin = n_cta - n_unbinned
-                msg = 'The intrinsic binning will be used for the remaining %s CTA observations.' % (
-                    n_notunbin)
-                self._log_string(gammalib.TERSE, msg)
  
         # If there is more than one observation, and observations are all
         # unbinned or all onoff query user to know if they wish stacked results
-        if self.obs().size() > 1 and (
-                        n_unbinned == self.obs().size() or n_onoff == self.obs().size()):
+        if self.obs().size() > 1 and \
+           (n_unbinned == self.obs().size() or n_onoff == self.obs().size()):
             self._stack = self['stack'].boolean()
             # If we are to stack event lists query parameters for cube creation
             if self._stack and n_unbinned == self.obs().size():
@@ -151,9 +135,7 @@ class csresspec(ctools.csobservation):
         # Unless all observations are On/Off, or we are using precomputed model
         # maps query whether to use energy dispersion
         if n_onoff == n_cta or self._use_maps:
-            msg = 'Energy dispersion is applied based on the input data/model ' \
-                  + 'and not according to the edisp parameter'
-            self._log_string(gammalib.TERSE, msg)
+            pass
         else:
             self['edisp'].boolean()
  
@@ -167,6 +149,35 @@ class csresspec(ctools.csobservation):
         # Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
  
+        # Write header for observation census
+        self._log_header1(gammalib.TERSE, 'Observation census')
+
+        # Log census of input observations
+        self._log_value(gammalib.NORMAL, 'Unbinned CTA observations', n_unbinned)
+        self._log_value(gammalib.NORMAL, 'Binned CTA observations', n_binned)
+        self._log_value(gammalib.NORMAL, 'On/off CTA observations', n_onoff)
+        self._log_value(gammalib.NORMAL, 'Other observations', n_other)
+        if n_other > 0:
+            msg = 'WARNING: Only CTA observation can be handled, all non-CTA ' \
+                  + 'observations will be ignored.'
+            self._log_string(gammalib.TERSE, msg)
+
+        # Log for unbinned observations
+        if n_unbinned != 0:
+            msg = ' User defined energy binning will be used for %d unbinned '\
+                  'observations.' % (n_unbinned)
+            self._log_string(gammalib.TERSE, msg)
+            if n_cta > n_unbinned:
+                msg = ' The intrinsic binning will be used for the remaining '\
+                      '%d CTA observations.' % (n_notunbin)
+                self._log_string(gammalib.TERSE, msg)
+
+        # Signal how energy dispersion is applied
+        if n_onoff == n_cta or self._use_maps:
+            msg = ' Energy dispersion is applied based on the input data/model ' \
+                  + 'and not according to the edisp parameter'
+            self._log_string(gammalib.TERSE, msg)
+
         # Return
         return
  
@@ -175,23 +186,26 @@ class csresspec(ctools.csobservation):
         Stack multiple observations and replace observation container
         with single stacked observation
  
-        :return: `~gammalib.GObservations' container with a single stacked observation
+        Returns
+        -------
+        obs : `~gammalib.GObservations'
+            Observation container with a single stacked observation
         """
         # Use first observation to determine the type and apply the
         # appropriate stacking method
  
         # If On/Off use the GCTAOnOffObservation constructor
         if self.obs()[0].classname() == 'GCTAOnOffObservation':
-            msg = 'Stacking %s On/Off observations.' % (self.obs().size())
+            msg = 'Stacking %d On/Off observations.' % (self.obs().size())
             self._log_string(gammalib.NORMAL, msg)
-            stacked = gammalib.GCTAOnOffObservation(self.obs())
+            stacked     = gammalib.GCTAOnOffObservation(self.obs())
             stacked_obs = gammalib.GObservations()
             stacked_obs.append(stacked)
             stacked_obs.models(self.obs().models())
  
         # If event list then bin observations
         elif self.obs()[0].classname() == 'GCTAObservation':
-            msg = 'Stacking %s event lists.' % (self.obs().size())
+            msg = 'Stacking %d event lists.' % (self.obs().size())
             self._log_string(gammalib.NORMAL, msg)
             stacked_obs = obsutils.get_stacked_obs(self, self.obs())
  
@@ -201,24 +215,31 @@ class csresspec(ctools.csobservation):
     def _bin_evlist(self, obs):
         """
         Turn single event list into counts cube
-        :param   obs: `~gammalib.GObservations' observation container with single event list
-        :return: `~gammalib.GObservations' binned observation container
-                 dict with event list ROI and energy range information
  
+        Parameters
+        ----------
+        obs : `~gammalib.GObservations()`
+            Observation container with single event list
+
+        Returns
+        -------
+        obs : `~gammalib.GObservations()`
+            Binned observation container dict with event list ROI and energy
+            range information
         """
         # Retrieve information about ROI in event list
         roi = obs[0].roi()
-        ra = roi.centre().dir().ra_deg()
+        ra  = roi.centre().dir().ra_deg()
         dec = roi.centre().dir().dec_deg()
         rad = roi.radius()
  
         # We will cover the whole ROI with 0.02 deg binning
-        npix = int(2 * rad / 0.02) + 1
+        npix = int(2.0 * rad / 0.02) + 1
+
+        # Log binning of events
+        self._log_string(gammalib.EXPLICIT, 'Binning events')
  
         # Bin events
-        msg = 'Binning events'
-        self._log_string(gammalib.EXPLICIT, msg)
-        #
         cntcube = ctools.ctbin(obs)
         cntcube['xref'] = ra
         cntcube['yref'] = dec
@@ -239,9 +260,10 @@ class csresspec(ctools.csobservation):
         # Retrieve the binned observation container
         binned_obs = cntcube.obs().copy()
  
+        # Log response computation
+        self._log_string(gammalib.EXPLICIT, 'Computing binned response')
+
         # Compute binned response
-        msg = 'Computing binned response'
-        self._log_string(gammalib.EXPLICIT, msg)
         response = obsutils.get_stacked_response(obs, ra, dec,
                                                  binsz=0.02,
                                                  nxpix=npix,
@@ -287,44 +309,68 @@ class csresspec(ctools.csobservation):
                      regfile='NONE'):
         """
         Mask an event cube and returns the masked cube
-        :param cube: `~gammalib.GCTAEventCube' Events cube
-        :param ra: double (str 'INDEF' for no selection on direction)
-        :param dec: double (str 'INDEF' for no selection on direction)
-        :param rad: double (str 'INDEF' for no selection on direction)
-        :param emin: double (str 'INDEF' for no selection on energy)
-        :param emax: double (str 'INDEF' for no selection on energy)
-        :param regfile: str with name of ds9 file with exclusion regions
-        :return: `~gammalib.GCTAEventCube'
-        """
  
+        Parameters
+        ----------
+        cube : `~gammalib.GCTAEventCube()`
+            Event cube
+        ra : float (str 'INDEF' for no selection on direction)
+            Right Ascension (deg)
+        dec : float (str 'INDEF' for no selection on direction)
+            Declination (deg)
+        rad : float (str 'INDEF' for no selection on direction)
+            Radius (deg)
+        emin : float (str 'INDEF' for no selection on energy)
+            Minimum energy (TeV)
+        emax : float (str 'INDEF' for no selection on energy)
+            Maximum energy (TeV)
+
+        Returns
+        -------
+        cube : `~gammalib.GCTAEventCube()`
+            Event cube
+        """
         # Turn cube into observation container to feed to ctcubemask
         obs = gammalib.GCTAObservation()
         obs.events(cube)
         obs_cont = gammalib.GObservations()
         obs_cont.append(obs)
  
-        # Use ctcubemask to mask
+        # Use ctcubemask to mask event cube pixels
         cubemask = ctools.ctcubemask(obs_cont)
-        cubemask['ra'] = ra
-        cubemask['dec'] = dec
-        cubemask['rad'] = rad
-        cubemask['emin'] = emin
-        cubemask['emax'] = emax
+        cubemask['ra']      = ra
+        cubemask['dec']     = dec
+        cubemask['rad']     = rad
+        cubemask['emin']    = emin
+        cubemask['emax']    = emax
         cubemask['regfile'] = regfile
         cubemask.run()
  
-        # Return
-        # Note: we return a copy to avoid memory leaks in SWIG
-        return cubemask.cube().copy()
+        # Extract copy of cube from observation container (copy is needed to
+        # avoid memory leaks in SWIG)
+        cube = cubemask.obs()[0].events().copy()
+
+        # Return cube
+        return cube
  
     def _cube_to_spectrum(self, cube, evlist_info):
         """
-        Derive from event cube a count spectrum. If data come from event list
-        use only the ROI and energy range of the original data. Apply user
-        defined mask if requested.
-        :param cube: `~gammalib.GCTAEventCube' input cube
-        :param evlist_info: dic with information on original event list
-        :return: `~gammalib.GNdarray' count spectrum
+        Derive from event cube a count spectrum
+        
+        If data come from event list use only the ROI and energy range of
+        the original data. Apply user defined mask if requested.
+
+        Parameters
+        ----------
+        cube : `~gammalib.GCTAEventCube()`
+            Event cube
+        evlist_info : dict
+            Dictionary with information on original event list
+
+        Returns
+        -------
+        array : `~gammalib.GNdarray()'
+            Counts spectrum
         """
         # If we started from event list mask the ROI only
         # for residual computation
@@ -350,7 +396,7 @@ class csresspec(ctools.csobservation):
         counts = counts.clip(0.)
  
         # Convert skymap into GNdarray count spectrum
-        counts = counts.total_counts()
+        counts = counts.counts()
  
         # Return
         return counts
@@ -358,30 +404,42 @@ class csresspec(ctools.csobservation):
     def _residuals_table(self, obs_id, ebounds, counts, model, residuals):
         """
         Create a Fits Table and store counts, model, and residuals
-        :param obs_id: str observation id
-        :param ebounds: `~gammalib.GEbounds' energy bounds
-        :param counts: `~gammalib.GNdarray'
-        :param model: `~gammalib.GNdarray'
-        :param residuals: `~gammalib.GNdarray'
-        :return: `~gammalib.GFitsBinTable'
+
+        Parameters
+        ----------
+        obs_id : str
+            Observation id
+        ebounds : `~gammalib.GEbounds()'
+            Energy boundaries
+        counts : `~gammalib.GNdarray()'
+            Counts spectrum
+        model : `~gammalib.GNdarray()'
+            Model spectrum
+        residuals : `~gammalib.GNdarray()'
+            Residual spectrum
+
+        Returns
+        -------
+        table : `~gammalib.GFitsBinTable()'
+            Residual spectrum as FITS binary table
         """
         # Create FITS table columns
-        nrows = ebounds.size()
-        energy_low = gammalib.GFitsTableDoubleCol('Emin', nrows)
+        nrows       = ebounds.size()
+        energy_low  = gammalib.GFitsTableDoubleCol('Emin', nrows)
         energy_high = gammalib.GFitsTableDoubleCol('Emax', nrows)
-        counts_col = gammalib.GFitsTableDoubleCol('Counts', nrows)
-        model_col = gammalib.GFitsTableDoubleCol('Model', nrows)
-        resid_col = gammalib.GFitsTableDoubleCol('Residuals', nrows)
+        counts_col  = gammalib.GFitsTableDoubleCol('Counts', nrows)
+        model_col   = gammalib.GFitsTableDoubleCol('Model', nrows)
+        resid_col   = gammalib.GFitsTableDoubleCol('Residuals', nrows)
         energy_low.unit('TeV')
         energy_high.unit('TeV')
  
         # Fill FITS table columns
         for i in range(nrows):
-            energy_low[i] = ebounds.emin(i).TeV()
+            energy_low[i]  = ebounds.emin(i).TeV()
             energy_high[i] = ebounds.emax(i).TeV()
-            counts_col[i] = counts[i]
-            model_col[i] = model[i]
-            resid_col[i] = residuals[i]
+            counts_col[i]  = counts[i]
+            model_col[i]   = model[i]
+            resid_col[i]   = residuals[i]
  
         # Initialise FITS Table with extension set to obs id
         table = gammalib.GFitsBinTable(nrows)
@@ -401,15 +459,26 @@ class csresspec(ctools.csobservation):
         table.append(model_col)
         table.append(resid_col)
  
+        # Return binary table
         return table
  
     def _append_column(self, table, name, data):
         """
         Append optional column to residual table
-        :param table: `~gammalib.GFitsBinTable'
-        :param name: str column name
-        :param data: `gammalib.GNdarray' data to be filled into new column
-        :return: `~gammalib.GFitsBinTable'
+
+        Parameters
+        ----------
+        table : `~gammalib.GFitsBinTable()'
+            FITS binary table
+        name : str
+            Column name
+        data : `~gammalib.GEbounds()'
+            Data to be filled into new column
+
+        Returns
+        -------
+        table : `~gammalib.GFitsBinTable()'
+            FITS binary table
         """
         # Check size compatibility
         if table.nrows() == data.size():
@@ -417,7 +486,7 @@ class csresspec(ctools.csobservation):
         # Otherwise throw error
         else:
             msg = 'csresspec._append_column: FITS table and data have ' \
-                  + 'incompatible size.'
+                  'incompatible size.'
             raise RuntimeError(msg)
  
         # Create column
@@ -433,6 +502,7 @@ class csresspec(ctools.csobservation):
         # Return modified table
         return table
  
+
     # Public methods
     def run(self):
         """
@@ -453,17 +523,21 @@ class csresspec(ctools.csobservation):
             self._log_header1(gammalib.TERSE, 'Stacking Observations')
             self.obs(self._stack_observations())
  
+        # Log processing header
         self._log_header1(gammalib.TERSE, 'Processing Observations')
+
         # Loop over observations and calculate residuals
         for s, observation in enumerate(self.obs()):
  
             # Retrieve and store obs id
             obs_id = observation.id()
+
             # If observation id is empty and there is more than one observation
             # replace with incremental number
             if obs_id == '' and self.obs().size() > 1:
                 obs_id = str(s)
  
+            # Log processing of observation
             if self.obs().size() > 1:
                 self._log_header2(gammalib.NORMAL,
                                   'Processing observation %s' % obs_id)
@@ -476,7 +550,7 @@ class csresspec(ctools.csobservation):
             # if 3D observation
             if obs[0].classname() == 'GCTAObservation':
  
-                ## Prepare Observations
+                # Prepare Observations
  
                 # If already binned set the evlist_info dictionary to have
                 # attribute was_list False
@@ -490,7 +564,7 @@ class csresspec(ctools.csobservation):
                     self._log_string(gammalib.NORMAL, msg)
                     obs, evlist_info = self._bin_evlist(obs)
  
-                ## Calculate Model and residuals
+                # Calculate Model and residuals
                 # If model cube is provided load it
                 if self._use_maps:
                     modcube = gammalib.GCTAEventCube(self['inmodel'].filename())
@@ -509,7 +583,7 @@ class csresspec(ctools.csobservation):
                 msg = 'Computing counts, model, and residual spectra'
                 self._log_string(gammalib.NORMAL, msg)
                 counts = self._cube_to_spectrum(cntcube, evlist_info)
-                model = self._cube_to_spectrum(modcube, evlist_info)
+                model  = self._cube_to_spectrum(modcube, evlist_info)
  
                 # Calculate residuals
                 residuals = obsutils.residuals(self, counts, model)
@@ -523,7 +597,7 @@ class csresspec(ctools.csobservation):
                 table = self._residuals_table(obs_id, ebounds, counts, model,
                                               residuals)
  
-                ## Calculate models of individual components if requested
+                # Calculate models of individual components if requested
                 if self['components'].boolean():
                     for component in obs.models():
                         self._log_value(gammalib.NORMAL,
@@ -551,7 +625,7 @@ class csresspec(ctools.csobservation):
  
                 onoff = obs[0]
  
-                ## Calculate Counts, Model and residuals
+                # Calculate Counts, Model and residuals
                 msg = 'Computing counts, model, and residual spectra'
                 self._log_string(gammalib.NORMAL, msg)
  
@@ -559,12 +633,11 @@ class csresspec(ctools.csobservation):
                 counts = onoff.on_spec().counts_spectrum()
  
                 # Model for On region
-                background = onoff.model_background(
-                    obs.models()).counts_spectrum()
-                alpha = onoff.on_spec().backscal_spectrum()
-                model = background.copy()
-                model *= alpha
-                model += onoff.model_gamma(obs.models()).counts_spectrum()
+                background = onoff.model_background(obs.models()).counts_spectrum()
+                alpha      = onoff.on_spec().backscal_spectrum()
+                model      = background.copy()
+                model     *= alpha
+                model     += onoff.model_gamma(obs.models()).counts_spectrum()
  
                 # On Residuals
                 residuals = obsutils.residuals(self, counts, model)
@@ -594,7 +667,7 @@ class csresspec(ctools.csobservation):
                 table = self._append_column(table, 'Residuals_Off',
                                             residuals_off)
  
-                ## Calculate models of individual components if requested
+                # Calculate models of individual components if requested
                 if self['components'].boolean():
                     for component in obs.models():
                         # If background pass
@@ -625,9 +698,12 @@ class csresspec(ctools.csobservation):
                     table = self._append_column(table, 'Background',
                                                 background)
  
-            ## Append results table to output file
+            # Append results table to output file
             self._fits.append(table)
  
+        # Return
+        return
+
     def save(self):
         """
         Save residuals
@@ -637,6 +713,7 @@ class csresspec(ctools.csobservation):
  
         # Continue only if FITS file is valid
         if self._fits != None:
+
             # Get outfile parameter
             outfile = self['outfile'].filename()
  
@@ -654,6 +731,7 @@ class csresspec(ctools.csobservation):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
+
     # Create instance of application
     app = csresspec(sys.argv)
  
 # ==========================================================================
 # Utility functions for observation handling
 #
-# Copyright (C) 2011-2017 Juergen Knoedlseder
+# Copyright (C) 2011-2018 Juergen Knoedlseder
 #
 # 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
@@ -829,22 +829,33 @@ def get_onoff_obs(cls, obs):
     # Return On/Off oberservation container
     return onoff_obs
  
+
 # ========================================= #
 # Calculate residuals from counts and model #
 # ========================================= #
-def residuals(cls,counts,model):
+def residuals(cls, counts, model):
     """
-    Calculate residuals given counts and models,
-    according to algorithm specified by user.
+    Calculate residuals given counts and models, according to algorithm
+    specified by user.
+
     Can handle GSkyMap or GNdarray objects
-    :param cls: `~ctools.cscript`
+
+
+    Parameters
+    ----------
+    cls : `~ctools.cscript`
         cscript class
-    :param counts: `~gammalib.GSkyMap/~gammalib.GNdarray' data counts
-    :param model: `~gammalib.GSkyMap/~gammalib.GNdarray' model counts
-    :return residuals: `~gammalib.GSkyMap/~gammalib.GNdarray' residuals
+    counts : `~gammalib.GSkyMap/~gammalib.GNdarray'
+        Data counts
+    model : `~gammalib.GSkyMap/~gammalib.GNdarray'
+        Model counts
+
+    Returns
+    -------
+    residuals : `~gammalib.GSkyMap/~gammalib.GNdarray'
+        Residuals
     """
-    # Find type of objects we are manipulating
-    # and set size to iterate later
+    # Find type of objects we are manipulating and set size to iterate later
  
     # If GNdarray
     if counts.classname() == 'GNdarray':
@@ -893,9 +904,8 @@ def residuals(cls,counts,model):
                 # significance^2 and save it ...
                 data_val = residuals[i]
                 if data_val > 0.0:
-                    log_val = math.log(data_val / model_val)
-                    residuals[i] = (
-                                      data_val * log_val) + model_val - data_val
+                    log_val      = math.log(data_val / model_val)
+                    residuals[i] = (data_val * log_val) + model_val - data_val
  
                 # ... otherwise compute the reduced value of the above
                 # expression. This is necessary to avoid computing log(0).
@@ -908,7 +918,7 @@ def residuals(cls,counts,model):
  
         # Compute significance map
         residuals *= 2.0
-        residuals = residuals.sqrt()
+        residuals  = residuals.sqrt()
         residuals *= sign
  
     # Raise exception if algorithm is unknown
@@ -916,4 +926,4 @@ def residuals(cls,counts,model):
         raise TypeError('Algorithm "' + algorithm + '" not known')
  
     # Return
-    return residuals
 \ No newline at end of file
+    return residuals
 #############################################################################
 # Makefile for ctools examples                                              #
 # ------------------------------------------------------------------------- #
-# Copyright (C) 2016-2017 Juergen Knoedlseder                               #
+# Copyright (C) 2016-2018 Juergen Knoedlseder                               #
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  This program is free software: you can redistribute it and/or modify     #
@@ -28,7 +28,7 @@ scripts = README.md \
           show_pha.py \
           show_pull_evolution.py \
           show_pull_histogram.py \
-	  show_residuals.py \
+          show_residuals.py \
           show_response.py \
           show_sensitivity.py \
           show_spectrum.py \
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Display residuals generated by csresspec
 #
-# Copyright (C) 2017- Luigi Tibaldo
+# Copyright (C) 2017-2018 Luigi Tibaldo
 #
 # 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
@@ -22,11 +22,9 @@ import sys
 import math
 import gammalib
 import cscripts
-
 try:
     import matplotlib.pyplot as plt
     from matplotlib import gridspec
-
     plt.figure()
     plt.close()
 except (ImportError, RuntimeError):
@@ -42,23 +40,41 @@ def fill_cmr(row, counts, model, resid, e_counts, e_resid, c_counts, c_model,
     """
     Helper function that fills in the counts, model, residuals lists and
     calculate relative errors for given table row
-    :param row: int row number
-    :param counts: list
-    :param model: list
-    :param resid: list
-    :param e_counts: list
-    :param e_resid: list
-    :param c_counts: `~gammalib.GFits.table.column'
-    :param c_model: `~gammalib.GFits.table.column'
-    :param c_resid: `~gammalib.GFits.table.column'
-    :param algorithm: string
-    :return: counts, model, resid, e_counts, e_resid
+
+    Parameters
+    ----------
+    row : int
+        Table row index
+    counts : list of float
+        Counts
+    model : list of float
+        Model
+    resid : list of float
+        Residual
+    e_counts : list of float
+        Counts error
+    e_resid : list of float
+        Residual error
+    c_counts : `~gammalib.GFitsTableCol()'
+        Counts FITS table column
+    c_model : `~gammalib.GFitsTableCol()'
+        Model FITS table column
+    c_resid : `~gammalib.GFitsTableCol()'
+        Residual FITS table column
+    algorithm : str
+        Algorithm
+
+    Returns
+    -------
+    counts, model, resid, e_counts, e_resid : tuple of float
+        Values
     """
+    # Extract values
     counts.append(c_counts.real(row))
     model.append(c_model.real(row))
     resid.append(c_resid.real(row))
  
-    # calculate count error
+    # Calculate count error
     err = math.sqrt(c_counts.real(row))
     if err == 1:
         # this prevents visualization problem in matplotlib
@@ -67,7 +83,7 @@ def fill_cmr(row, counts, model, resid, e_counts, e_resid, c_counts, c_model,
         pass
     e_counts.append(err)
  
-    # calculate residual error
+    # Calculate residual error
     if algorithm == 'SUB':
         e_resid.append(err)
     elif algorithm == 'SUBDIV':
@@ -77,6 +93,7 @@ def fill_cmr(row, counts, model, resid, e_counts, e_resid, c_counts, c_model,
     elif algorithm == 'SIGNIFICANCE':
         e_resid.append(1.)
  
+    # Return tuple
     return counts, model, resid, e_counts, e_resid
  
  
@@ -97,59 +114,63 @@ def plot_residuals(filename, plotfile, hdu):
         Number of observation to plot
     """
     # Read spectrum file    
-    fits = gammalib.GFits(filename)
-    table = fits.table(1 + hdu)
-    c_emin = table['Emin']
-    c_emax = table['Emax']
+    fits     = gammalib.GFits(filename)
+    table    = fits.table(1 + hdu)
+    c_emin   = table['Emin']
+    c_emax   = table['Emax']
     c_counts = table['Counts']
-    c_model = table['Model']
-    c_resid = table['Residuals']
+    c_model  = table['Model']
+    c_resid  = table['Residuals']
     try:
         c_counts_off = table['Counts_Off']
-        c_model_off = table['Model_Off']
-        c_resid_off = table['Residuals_Off']
-        is_onoff = True
+        c_model_off  = table['Model_Off']
+        c_resid_off  = table['Residuals_Off']
+        is_onoff     = True
     except:
         is_onoff = False
  
     # Initialise arrays to be filled
-    ebounds = []
-    ed_engs = []
-    eu_engs = []
-    em_engs = []
-    counts = []
-    e_counts = []
-    model = []
-    resid = []
-    e_resid = []
-    counts_off = []
+    ebounds      = []
+    ed_engs      = []
+    eu_engs      = []
+    em_engs      = []
+    counts       = []
+    e_counts     = []
+    model        = []
+    resid        = []
+    e_resid      = []
+    counts_off   = []
     e_counts_off = []
-    model_off = []
-    resid_off = []
-    e_resid_off = []
+    model_off    = []
+    resid_off    = []
+    e_resid_off  = []
  
     # Residual algorithm
     algorithm = table.header()['ALGORITHM'].string()
  
     # Loop over rows of the file
     nrows = table.nrows()
-    # add first energy boundary
+
+    # Add first energy boundary
     ebounds.append(c_emin.real(0))
     for row in range(nrows):
-        # boundaries
+
+        # Boundaries
         ebounds.append(c_emax.real(row))
-        # geometrical mean energy and errors
+
+        # Geometrical mean energy and errors
         emean = math.sqrt(c_emin.real(row) * c_emax.real(row))
         em_engs.append(emean)
         ed_engs.append(emean - c_emin.real(row))
         eu_engs.append(c_emax.real(row) - emean)
+
         # counts, model, residuals
         counts, model, resid, e_counts, e_resid = fill_cmr(row, counts, model,
                                                            resid, e_counts,
                                                            e_resid, c_counts,
                                                            c_model, c_resid,
                                                            algorithm)
-        # onoff
+        # On/off
         if is_onoff:
             counts_off, model_off, resid_off, e_counts_off, e_resid_off = fill_cmr(
                 row, counts_off,
@@ -168,7 +189,7 @@ def plot_residuals(filename, plotfile, hdu):
     # Plot
     axarr = []
     if is_onoff:
-        f = plt.figure(figsize=(8, 4))
+        f  = plt.figure(figsize=(8, 4))
         gs = gridspec.GridSpec(2, 2, height_ratios=[2, 1])
         gs.update(hspace=0, left=0.09, right=0.97)
         axarr.append(f.add_subplot(gs[0]))
@@ -176,7 +197,7 @@ def plot_residuals(filename, plotfile, hdu):
         axarr.append(f.add_subplot(gs[1], sharex=axarr[0]))
         axarr.append(f.add_subplot(gs[3], sharex=axarr[1]))
     else:
-        f = plt.figure()
+        f  = plt.figure()
         gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
         gs.update(hspace=0)
         axarr.append(f.add_subplot(gs[0]))
@@ -185,12 +206,14 @@ def plot_residuals(filename, plotfile, hdu):
     axarr[0].set_xscale('log')
     axarr[1].set_xscale('log')
     axarr[0].set_yscale('log')
+
     # Counts and model
     axarr[0].errorbar(em_engs, counts, yerr=e_counts, xerr=[ed_engs, eu_engs],
                       fmt='ko', capsize=0, linewidth=2, zorder=2, label='Data')
     axarr[0].step(ebounds, model, color='0.5', linewidth=2, zorder=1,
                   label='Model')
     axarr[0].set_ylabel('Counts')
+
     # Residuals
     axarr[1].errorbar(em_engs, resid, yerr=e_resid, xerr=[ed_engs, eu_engs],
                       fmt='ko', capsize=0, linewidth=2, zorder=2)
@@ -202,6 +225,7 @@ def plot_residuals(filename, plotfile, hdu):
         axarr[1].set_ylabel('Residuals (fraction)')
     elif algorithm == 'SUBDIVSQRT' or algorithm == 'SIGNIFICANCE':
         axarr[1].set_ylabel(r'Residuals ($\sigma$)')
+
     # On/Off
     if is_onoff:
         axarr[0].set_title('ON')
@@ -210,6 +234,7 @@ def plot_residuals(filename, plotfile, hdu):
         axarr[2].set_xscale('log')
         axarr[3].set_xscale('log')
         axarr[2].set_yscale('log')
+
         # Counts and model
         axarr[2].errorbar(em_engs, counts_off, yerr=e_counts_off,
                           xerr=[ed_engs, eu_engs],
@@ -217,6 +242,7 @@ def plot_residuals(filename, plotfile, hdu):
                           label='Data')
         axarr[2].step(ebounds, model_off, color='0.5', linewidth=2, zorder=1,
                       label='Model')
+
         # Residuals
         axarr[3].errorbar(em_engs, resid_off, yerr=e_resid_off,
                           xerr=[ed_engs, eu_engs],
@@ -284,5 +310,6 @@ def show_residuals():
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    # Show spectrum
+
+    # Show residuals
     show_residuals()
@@ -48,7 +48,6 @@ public:
     void run(void);
     void save(void);
     void publish(const std::string& name = "");
-    GCTAEventCube* cube(void);
 };
  
  
@@ -62,7 +62,6 @@ public:
     void run(void);
     void save(void);
     void publish(const std::string& name = "");
-    GCTAEventCube* cube(void) const;
  
 protected:
     // Protected methods
@@ -89,19 +88,4 @@ protected:
     bool                     m_select_energy; //!< Perform energy selection
 };
  
-/***********************************************************************//**
- * @brief Return events cube
- *
- * @return Reference to events cube correspoding to ith observation
- *
- * Returns a reference to the events cube.
- ***************************************************************************/
-inline
-GCTAEventCube* ctcubemask::cube(void) const
-{
-
-  return dynamic_cast<GCTAEventCube*>(const_cast<GEvents*>(m_obs[0]->events()));
-  
-}
-
 #endif /* CTCUBEMASK_HPP */
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csresspec script.
 #
-# Copyright (C) 2017- Luigi Tibaldo
+# Copyright (C) 2017-2018 Luigi Tibaldo
 #
 # 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
@@ -43,7 +43,7 @@ class Test(test):
         test.__init__(self)
  
         # Set members
-        self._inonoff = self._datadir + '/onoff_obs.xml'
+        self._inonoff     = self._datadir + '/onoff_obs.xml'
         self._onoff_model = self._datadir + '/onoff_model.xml'
  
         # Return
@@ -74,12 +74,12 @@ class Test(test):
  
         # Setup csresspec command
         cmd = csresspec + ' inobs="' + self._events + '"' + \
-              ' outfile="csresspec_cmd1.fits"' + \
-              ' inmodel="' + self._model + '"' + \
-              ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
-              ' emin=0.1 emax=100.0 enumbins=20' + \
-              ' mask=no algorithm="SUB"' + \
-              ' logfile="csresspec_cmd1.log" chatter=4'
+                          ' outfile="csresspec_cmd1.fits"' + \
+                          ' inmodel="' + self._model + '"' + \
+                          ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
+                          ' emin=0.1 emax=100.0 enumbins=20' + \
+                          ' mask=no algorithm="SUB"' + \
+                          ' logfile="csresspec_cmd1.log" chatter=4'
  
         # Check if execution of wrong command fails
         self.test_assert(self._execute('command_that_does_not_exist') != 0,
@@ -94,12 +94,12 @@ class Test(test):
  
         # Setup csresspec command
         cmd = csresspec + ' inobs="event_file_that_does_not_exist.fits"' + \
-              ' outfile="csresspec_cmd2.fits"' + \
-              ' inmodel="' + self._model + '"' + \
-              ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
-              ' emin=0.1 emax=100.0 enumbins=20' + \
-              ' mask=no algorithm="SUBDIVSQRT"' + \
-              ' logfile="csresspec_cmd2.log" debug=yes chatter=2'
+                          ' outfile="csresspec_cmd2.fits"' + \
+                          ' inmodel="' + self._model + '"' + \
+                          ' caldb="' + self._caldb + '" irf="' + self._irf + '"' + \
+                          ' emin=0.1 emax=100.0 enumbins=20' + \
+                          ' mask=no algorithm="SUBDIVSQRT"' + \
+                          ' logfile="csresspec_cmd2.log" debug=yes chatter=2'
  
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
@@ -118,18 +118,18 @@ class Test(test):
         """
         # Set-up csresspec for event list
         resspec = cscripts.csresspec()
-        resspec['inobs'] = self._events
-        resspec['outfile'] = 'csresspec_py1.fits'
-        resspec['inmodel'] = self._model
-        resspec['caldb'] = self._caldb
-        resspec['irf'] = self._irf
-        resspec['emin'] = 0.1
-        resspec['emax'] = 100.0
-        resspec['enumbins'] = 20
-        resspec['mask'] = False
+        resspec['inobs']     = self._events
+        resspec['outfile']   = 'csresspec_py1.fits'
+        resspec['inmodel']   = self._model
+        resspec['caldb']     = self._caldb
+        resspec['irf']       = self._irf
+        resspec['emin']      = 0.1
+        resspec['emax']      = 100.0
+        resspec['enumbins']  = 20
+        resspec['mask']      = False
         resspec['algorithm'] = 'SUB'
-        resspec['logfile'] = 'csresspec_py1.log'
-        resspec['chatter'] = 2
+        resspec['logfile']   = 'csresspec_py1.log'
+        resspec['chatter']   = 2
  
         # Run csresspec script
         resspec.logFileOpen()  # Make sure we get a log file
@@ -142,19 +142,19 @@ class Test(test):
         # Set-up csresspec for event list
         # with computation of individual model components
         resspec = cscripts.csresspec()
-        resspec['inobs'] = self._events
-        resspec['outfile'] = 'csresspec_py2.fits'
-        resspec['inmodel'] = self._model
-        resspec['caldb'] = self._caldb
-        resspec['irf'] = self._irf
-        resspec['emin'] = 0.1
-        resspec['emax'] = 100.0
-        resspec['enumbins'] = 20
-        resspec['mask'] = False
-        resspec['algorithm'] = 'SUBDIV'
+        resspec['inobs']      = self._events
+        resspec['outfile']    = 'csresspec_py2.fits'
+        resspec['inmodel']    = self._model
+        resspec['caldb']      = self._caldb
+        resspec['irf']        = self._irf
+        resspec['emin']       = 0.1
+        resspec['emax']       = 100.0
+        resspec['enumbins']   = 20
+        resspec['mask']       = False
+        resspec['algorithm']  = 'SUBDIV'
         resspec['components'] = True
-        resspec['logfile'] = 'csresspec_py2.log'
-        resspec['chatter'] = 2
+        resspec['logfile']    = 'csresspec_py2.log'
+        resspec['chatter']    = 2
  
         # Run csresspec script
         resspec.execute()
@@ -164,20 +164,20 @@ class Test(test):
  
         # Set-up csresspec for counts cube
         resspec = cscripts.csresspec()
-        resspec['inobs'] = self._cntcube
-        resspec['modcube'] = 'NONE'
-        resspec['expcube'] = 'NONE'
-        resspec['psfcube'] = 'NONE'
+        resspec['inobs']     = self._cntcube
+        resspec['modcube']   = 'NONE'
+        resspec['expcube']   = 'NONE'
+        resspec['psfcube']   = 'NONE'
         resspec['edispcube'] = 'NONE'
-        resspec['bkgcube'] = 'NONE'
-        resspec['caldb'] = self._caldb
-        resspec['irf'] = self._irf
-        resspec['inmodel'] = self._model
-        resspec['outfile'] = 'csresspec_py3.fits'
-        resspec['mask'] = False
+        resspec['bkgcube']   = 'NONE'
+        resspec['caldb']     = self._caldb
+        resspec['irf']       = self._irf
+        resspec['inmodel']   = self._model
+        resspec['outfile']   = 'csresspec_py3.fits'
+        resspec['mask']      = False
         resspec['algorithm'] = 'SUBDIVSQRT'
-        resspec['logfile'] = 'csresspec_py3.log'
-        resspec['chatter'] = 3
+        resspec['logfile']   = 'csresspec_py3.log'
+        resspec['chatter']   = 3
  
         # Run csresspec script
         resspec.execute()
@@ -188,24 +188,24 @@ class Test(test):
         # Set-up csresspec for counts cube
         # with source mask
         resspec = cscripts.csresspec()
-        resspec['inobs'] = self._cntcube
-        resspec['modcube'] = 'NONE'
-        resspec['expcube'] = 'NONE'
-        resspec['psfcube'] = 'NONE'
+        resspec['inobs']     = self._cntcube
+        resspec['modcube']   = 'NONE'
+        resspec['expcube']   = 'NONE'
+        resspec['psfcube']   = 'NONE'
         resspec['edispcube'] = 'NONE'
-        resspec['bkgcube'] = 'NONE'
-        resspec['caldb'] = self._caldb
-        resspec['irf'] = self._irf
-        resspec['inmodel'] = self._model
-        resspec['outfile'] = 'csresspec_py4.fits'
-        resspec['mask'] = True
-        resspec['ra'] = 83.63
-        resspec['dec'] = 22.01
-        resspec['rad'] = 0.2
-        resspec['regfile'] = 'NONE'
+        resspec['bkgcube']   = 'NONE'
+        resspec['caldb']     = self._caldb
+        resspec['irf']       = self._irf
+        resspec['inmodel']   = self._model
+        resspec['outfile']   = 'csresspec_py4.fits'
+        resspec['mask']      = True
+        resspec['ra']        = 83.63
+        resspec['dec']       = 22.01
+        resspec['rad']       = 0.2
+        resspec['regfile']   = 'NONE'
         resspec['algorithm'] = 'SUBDIVSQRT'
-        resspec['logfile'] = 'csresspec_py4.log'
-        resspec['chatter'] = 3
+        resspec['logfile']   = 'csresspec_py4.log'
+        resspec['chatter']   = 3
  
         # Run csresspec script
         resspec.execute()
@@ -215,13 +215,13 @@ class Test(test):
  
         # Set-up csresspec for On/Off observations
         resspec = cscripts.csresspec()
-        resspec['inobs'] = self._inonoff
-        resspec['inmodel'] = self._onoff_model
-        resspec['outfile'] = 'csresspec_py5.fits'
+        resspec['inobs']     = self._inonoff
+        resspec['inmodel']   = self._onoff_model
+        resspec['outfile']   = 'csresspec_py5.fits'
         resspec['algorithm'] = 'SIGNIFICANCE'
-        resspec['stack'] = False
-        resspec['logfile'] = 'csresspec_py5.log'
-        resspec['chatter'] = 3
+        resspec['stack']     = False
+        resspec['logfile']   = 'csresspec_py5.log'
+        resspec['chatter']   = 3
  
         # Run csresspec script
         resspec.execute()
@@ -232,14 +232,14 @@ class Test(test):
         # Set-up csresspec for On/Off observations
         # in stacked mode and with WSTAT background
         resspec = cscripts.csresspec()
-        resspec['inobs'] = self._inonoff
-        resspec['inmodel'] = self._onoff_model
-        resspec['outfile'] = 'csresspec_py6.fits'
+        resspec['inobs']     = self._inonoff
+        resspec['inmodel']   = self._onoff_model
+        resspec['outfile']   = 'csresspec_py6.fits'
         resspec['algorithm'] = 'SIGNIFICANCE'
         resspec['statistic'] = 'WSTAT'
-        resspec['stack'] = True
-        resspec['logfile'] = 'csresspec_py6.log'
-        resspec['chatter'] = 3
+        resspec['stack']     = True
+        resspec['logfile']   = 'csresspec_py6.log'
+        resspec['chatter']   = 3
  
         # Run csresspec script
         resspec.execute()
@@ -263,9 +263,11 @@ class Test(test):
  
         # Check first hdu with residuals
         table = fits[1]
+
         # Check number of rows = number of energy bins
         self.test_value(table.nrows(), enumbins,
                         'Check for %d rows in residuals' % enumbins)
+
         # Check number of columns, 2 for energy bounds + 3 for 3D or 6 for
         # On/Off + 1 for each model component if components = yes
         self.test_value(table.ncols(), ncols,