Merge branch 'devel' into pull-request-1982-add-phase-info-ctphase

# Conflicts: # src/Makefile.am

Merge branch 'devel' into pull-request-1982-add-phase-info-ctphase
# Conflicts: # src/Makefile.am
Jürgen Knödlseder
2 parents 70c82fc3 ee25e9a8
Showing 73 changed files with 2785 additions and 1366 deletions
.gitignore
AUTHORS
ChangeLog
NEWS
cscripts/csroot2caldb.py
cscripts/cssens.par
cscripts/cssens.py
doc/source/users/reference_manual/cssens.rst
doc/source/users/reference_manual/ctbutterfly.rst
doc/source/users/reference_manual/ctlike.rst
doc/source/users/reference_manual/ctprob.rst
doc/source/users/reference_manual/index.rst
pyext/Makefile.am
pyext/ctool.i
pyext/ctools/tools.i
pyext/ctphase.i
pyext/ctprob.i
src/Makefile.am
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
@@ -71,6 +71,8 @@ src/ctcubemask/ctcubemask
 src/ctbutterfly/ctbutterfly
 src/ctulimit/ctulimit
 src/cterror/cterror
+src/ctprob/ctprob
+src/ctphase/ctphase
 test/pfiles/
 test/caldb/
 test/prod3
@@ -2,6 +2,7 @@ Johan Bregeon (johan.bregeon@lupm.in2p3.fr)
 Rolf Buehler (rolf.buehler@desy.de)
 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
+Leonardo Di Venere (leonardo.divenere@ba.infn.it)
 Florent Forest (florent.forest@isae.fr)
 Lucie Gerard (lgerard@wgs18.zeuthen.desy.de)
 Maria Haupt (maria.haupt@desy.de)
-2017-04-05
+2017-04-11
  
         * Version 1.3.0 released
         ========================
  
+        Add ctprob tool (#1991)
+        ctlike stores covariance matrix in either CSV or FITS format (#1863)
+        Use correct ROOT histograms in csroot2caldb
+        Implement regression scheme in cssens and remove num_avg parameter (#1989)
+        Implement Gaussian error propagation in ctbutterfly (#1862)
         Copy cscripts into bin directory instead of creating a symbolic links (#1987)
         If code originates from Git, always compile the Python bindings (#1985)
         Raise exception instead of exiting with return code in Python tests (#1983)
 New Features and Important Changes in ctools 1.3.0.dev1
  
-5 April 2017
+11 April 2017
  
  
 Introduction
@@ -57,7 +57,10 @@ None
  
 ctlike - Maximum likelihood fitting
 -----------------------------------
-None
+The tool now stores the covariance matrix in either FITS or CSV format,
+depending on the file type specified in the file name. For file names
+terminating with ".fits" or ".fit" a FITS file will be written, otherwise
+a CSV file will be written (#1863).
  
  
 cttsmap - Generation of Test Statistic map
@@ -113,7 +116,11 @@ None
  
 ctbutterfly - Generate butterfly diagrams
 -----------------------------------------
-None
+The tool now uses Gaussian error propagation as the standard method for the
+computation of the butterfly diagram. The envelope method that existed before
+is still available for power law spectral models as an option. The "method"
+parameter has been added to allow the selection of the computation method
+(#1862).
  
  
 ctulimit - Compute upper limit
@@ -126,6 +133,15 @@ cterror - Likelihood profile error computation
 None
  
  
+ctprob - Event probability computation
+--------------------------------------
+The ctprob tool has been added that computates for each event the probability
+that the event originates from a given model component. The tool takes on input
+either an event list or an observation definition XML file. It appends to each
+event list a column for each model component that will contain the probability
+that the event originates from that component (#1991).
+
+
 cscripts
 --------
 The csebins script has been added that generates an energy boundaries
@@ -133,6 +149,15 @@ extension that can be used for a stacked analysis. The energy boundaries
 are set as function of the variation in the effective area and background
 rates of the Instrument Response Function (#1976).
  
+A regression scheme now replaces the averaging scheme in cssens as this
+allows more accurately the prediction of the flux value at the threshold
+Test Statistic. The num_avg parameter has been removed (#1989).
+
+The names of some ROOT histograms has been changed in the Monte Carlo data
+and the "EffectiveAreaEtrue_offaxis" needed to be changed to
+"EffectiveAreaEtrueNoTheta2cut_offaxis" and "EffectiveAreaEtrue" to
+"EffectiveAreaEtrueNoTheta2cut" in csroot2caldb.
+
  
 Examples
 --------
@@ -536,7 +536,8 @@ class csroot2caldb(ctools.cscript):
         Translate ROOT to CALDB effective area extension
  
         The following ROOT histograms are used:
-        - EffectiveAreaEtrue_offaxis -> EFFAREA
+        - EffectiveAreaEtrueNoTheta2cut_offaxis -> EFFAREA
+        - EffectiveAreaEtrueNoTheta2cut (for normalization)
  
         Parameters
         ----------
@@ -551,13 +552,13 @@ class csroot2caldb(ctools.cscript):
         self._log_header1(gammalib.TERSE, 'Generate effective area extension')
  
         # Get relevant ROOT histograms
-        etrue = tfile.Get('EffectiveAreaEtrue_offaxis')
+        etrue = tfile.Get('EffectiveAreaEtrueNoTheta2cut_offaxis')
  
         # If requested then normalize the 2D histogram on the on-axis 1D
         # histogram. This assures that the 2D histogram has the same on-axis
         # effective area dependence as the 1D histogram.
         if self['norm1d'].boolean():
-            etrue_1D = tfile.Get('EffectiveAreaEtrue')
+            etrue_1D = tfile.Get('EffectiveAreaEtrueNoTheta2cut')
             self._renorm_onaxis(etrue, etrue_1D)
  
         # If requested then rebin the effective area Etrue histogram. This
@@ -20,7 +20,7 @@ inobs,   f, h, NONE,,, &quot;Input event list, counts cube, or observation definition
 inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
 srcname, s, a, Crab,,, "Source name"
 caldb,   s, a, prod2,,, "Calibration database"
-irf,     s, a, South_0.5h,,, "Instrument response function"
+irf,     s, a, South_50h,,, "Instrument response function"
 edisp,   b, h, no,,, "Apply energy dispersion?"
 deadc,   r, h, 0.98,0,1, "Deadtime correction factor"
 outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file"
@@ -31,16 +31,15 @@ outfile, f, h, sensitivity.dat,,, &quot;Output sensitivity ASCII file&quot;
 offset,   r, h, 0.0,,, "Source offset angle (deg)"
 duration, r, a, 180000.0,,, "Effective exposure time (s)"
 rad,      r, a, 5.0,,, "Radius of ROI (deg)"
-emin,     r, a, 0.020,,, "Lower energy limit (TeV)"
-emax,     r, a, 200.0,,, "Upper energy limit (TeV)"
-bins,     i, a, 21,,, "Number of energy bins for differential sensitivity computation"
+emin,     r, a, 0.03163,,, "Lower energy limit (TeV)"
+emax,     r, a, 125.8925,,, "Upper energy limit (TeV)"
+bins,     i, a, 18,,, "Number of energy bins for differential sensitivity computation"
 enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)"
 npix,     i, h, 200,,, "Number of pixels for binned analysis"
 binsz,    r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)"
 type,     s, h, Differential,Differential|Integral,, "Sensitivity type"
 sigma,    r, h, 5.0,,, "Significance threshold (Gaussian sigma)"
 max_iter, i, h, 50,,, "Maximum number of iterations"
-num_avg,  i, h, 3,,, "Number of iterations for sliding average"
  
 #
 # Standard parameters
@@ -49,6 +48,6 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
 debug,   b, h, no,,, "Debugging mode activated"
 mode,    s, h, ql,,, "Mode of automatic parameters"
-logfile, f, h, cspull.log,,, "Log filename"
+logfile, f, h, cssens.log,,, "Log filename"
  
  
 #! /usr/bin/env python
 # ==========================================================================
-# Computes the array sensitivity using a test source
+# Computes the array sensitivity using the Test Statistic for a test source
 #
-# Copyright (C) 2011-2016 Juergen Knoedlseder
+# Copyright (C) 2011-2017 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,7 +40,7 @@ class cssens(ctools.cscript):
     and the flux level is determined that leads to a particular significance.
  
     The significance is determined using the Test Statistic, defined as twice
-    the likelihood difference between fitting with and  without the test source.
+    the likelihood difference between fitting with and without the test source.
     """
  
     # Constructor
@@ -50,31 +50,20 @@ class cssens(ctools.cscript):
         """
         # Set name
         self._name    = 'cssens'
-        self._version = '1.2.0'
+        self._version = '1.3.0'
  
         # Initialise class members
-        self._obs         = gammalib.GObservations()
         self._ebounds     = gammalib.GEbounds()
         self._obs_ebounds = []
         self._srcname     = ''
-        self._outfile     = gammalib.GFilename()
         self._ra          = None
         self._dec         = None
-        self._edisp       = False
-        self._emin        = 0.020
-        self._emax        = 200.0
-        self._bins        = 21
-        self._enumbins    = 0
-        self._npix        = 200
-        self._binsz       = 0.05
-        self._type        = 'Differential'
-        self._ts_thres    = 25.0
-        self._max_iter    = 50
-        self._num_avg     = 3
         self._log_clients = False
-        
-        # Initialise application by calling the appropriate base class
-        # constructor
+
+        # Initialise observation container from constructor arguments
+        self._obs, argv = self._set_input_obs(argv)
+
+        # Initialise application by calling the appropriate class constructor
         self._init_cscript(argv)
  
         # Return
@@ -87,8 +76,8 @@ class cssens(ctools.cscript):
         Get user parameters from parfile
         """
         # Set observation if not done before
-        if self._obs == None or self._obs.size() == 0:
-            self._obs = self._set_obs()
+        if self._obs.size() == 0:
+            self._obs = self._set_obs(self['emin'].real(), self['emax'].real())
  
         # Set models if we have none
         if self._obs.models().size() == 0:
@@ -98,31 +87,28 @@ class cssens(ctools.cscript):
         self._srcname = self['srcname'].string()
  
         # Read further parameters
-        self._outfile = self['outfile'].filename()
-        self._emin    = self['emin'].real()
-        self._emax    = self['emax'].real()
-        self._bins    = self['bins'].integer()
-
-        # Read parameters for binned if requested
-        self._enumbins = self['enumbins'].integer()
-        if not self._enumbins == 0:
-            self._npix  = self['npix'].integer()
-            self._binsz = self['binsz'].real()
-
-        # Read remaining parameters
-        self._edisp    = self['edisp'].boolean()
-        self._ts_thres = self['sigma'].real() * self['sigma'].real()
-        self._max_iter = self['max_iter'].integer()
-        self._num_avg  = self['num_avg'].integer()
-        self._type     = self['type'].string()
-
-        # Query remaining parameters
+        emin = self['emin'].real()
+        emax = self['emax'].real()
+        bins = self['bins'].integer()
+
+        # Query parameters for binned if requested
+        enumbins = self['enumbins'].integer()
+        if not enumbins == 0:
+            self['npix'].integer()
+            self['binsz'].real()
+
+        # Query input parameters
+        self['sigma'].real()
+        self['max_iter'].integer()
+        self['type'].string()
+        self['outfile'].filename()
+        self['edisp'].boolean()
         self['debug'].boolean()
  
         # Derive some parameters
-        self._ebounds = gammalib.GEbounds(self._bins,
-                                          gammalib.GEnergy(self._emin, 'TeV'),
-                                          gammalib.GEnergy(self._emax, 'TeV'))
+        self._ebounds = gammalib.GEbounds(bins,
+                                          gammalib.GEnergy(emin, 'TeV'),
+                                          gammalib.GEnergy(emax, 'TeV'))
  
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
@@ -130,20 +116,20 @@ class cssens(ctools.cscript):
         # Return
         return
  
-    def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
+    def _set_obs(self, emin, emax, lpnt=0.0, bpnt=0.0):
         """
         Set an observation container
  
         Parameters
         ----------
+        emin : float
+            Minimum energy (TeV)
+        emax : float
+            Maximum energy (TeV)
         lpnt : float, optional
             Galactic longitude of pointing (deg)
         bpnt : float, optional
             Galactic latitude of pointing (deg)
-        emin : float, optional
-            Minimum energy (TeV)
-        emax : float, optional
-            Maximum energy (TeV)
  
         Returns
         -------
@@ -205,23 +191,23 @@ class cssens(ctools.cscript):
  
             # Get energy boundaries of the observation
             obs_ebounds = self._obs_ebounds[i]
-    
+
             # Get minimum and maximum energy of the observation
             obs_emin = obs_ebounds.emin()
             obs_emax = obs_ebounds.emax()
-            
+
             # If [emin,emax] is fully contained in the observation energy range
             # the use [emin,emax] as energy boundaries
             if obs_emin <= emin and obs_emax >= emax:
                 ebounds = gammalib.GEbounds(emin, emax)
-            
+
             # ... otherwise, if [emin,emax] is completely outside the
             # observation energy range then set the energy boundaries to the
             # zero-width interval [0,0]
             elif emax < obs_emin or emin > obs_emax:
                 e0      = gammalib.GEnergy(0.0, 'TeV')
                 ebounds = gammalib.GEbounds(e0, e0)
-            
+
             # ... otherwise, if [emin,emax] overlaps partially with the
             # observation energy range then set the energy boundaries to the
             # overlapping part
@@ -229,10 +215,10 @@ class cssens(ctools.cscript):
                 # Set overlapping energy range
                 set_emin = max(emin, obs_emin)
                 set_emax = min(emax, obs_emax)
-                
+
                 # Set energy boundaries
                 ebounds = gammalib.GEbounds(set_emin, set_emax)
-            
+
             # Set the energy boundaries as the boundaries of the observation
             obs.events().ebounds(ebounds)
  
@@ -286,13 +272,27 @@ class cssens(ctools.cscript):
         # Set TeV->erg conversion factor
         tev2erg = 1.6021764
  
+        # Set parameters
+        ts_thres = self['sigma'].real() * self['sigma'].real()
+        max_iter = self['max_iter'].integer()
+        enumbins = self['enumbins'].integer()
+        if not enumbins == 0:
+            npix  = self['npix'].integer()
+            binsz = self['binsz'].real()
+        else:
+            npix  = 200
+            binsz = 0.05
+
+        # Set flux ratio precision required for convergence to 5%
+        ratio_precision = 0.05
+
         # Set energy boundaries
         self._set_obs_ebounds(emin, emax)
  
-        # Determine energy boundaries from first observation in the container
-        loge      = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
-        e_mean    = math.pow(10.0, loge)
-        erg_mean  = e_mean * tev2erg
+        # Determine mean energy for energy boundary
+        e_mean   = math.sqrt(emin.TeV()*emax.TeV())
+        loge     = math.log10(e_mean)
+        erg_mean = e_mean * tev2erg
  
         # Compute Crab unit. This is the factor with which the Prefactor needs
         # to be multiplied to get 1 Crab.
@@ -300,20 +300,27 @@ class cssens(ctools.cscript):
         src_flux  = test_model[self._srcname].spectral().flux(emin, emax)
         crab_unit = crab_flux/src_flux
  
-        # Write header and some introductory information into logger
+        # Initialise regression coefficient
+        regcoeff = 0.0
+
+        # Write header for energy bin
+        self._log_string(gammalib.TERSE, '')
         self._log_header2(gammalib.TERSE, 'Energies: '+str(emin)+' - '+str(emax))
+
+        # Write initial parameters
+        self._log_header3(gammalib.TERSE, 'Initial parameters')
         self._log_value(gammalib.TERSE, 'Crab flux', str(crab_flux)+' ph/cm2/s')
-        self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+
-                        ' ph/cm2/s')
+        self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+' ph/cm2/s')
         self._log_value(gammalib.TERSE, 'Crab unit factor', crab_unit)
  
         # Initialise loop
-        crab_flux_value   = []
-        photon_flux_value = []
-        energy_flux_value = []
-        sensitivity_value = []
-        iterations        = 0
-        test_crab_flux    = 0.1 # Initial test flux in Crab units (100 mCrab)
+        results        = []
+        iterations     = 0
+        test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab)
+
+        # Write header for iterations for terse chatter level
+        if self._logTerse():
+            self._log_header3(gammalib.TERSE, 'Iterations')
  
         # Loop until we break
         while True:
@@ -335,11 +342,11 @@ class cssens(ctools.cscript):
  
             # Simulate events for the models. "sim" holds an observation
             # container with observations containing the simulated events.
-            sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations,
-                               binsz=self._binsz, npix=self._npix,
+            sim = obsutils.sim(self._obs, nbins=enumbins, seed=iterations,
+                               binsz=binsz, npix=npix,
                                log=self._log_clients,
                                debug=self['debug'].boolean(),
-                               edisp=self._edisp)
+                               edisp=self['edisp'].boolean())
  
             # Determine number of events in simulation by summing the events
             # over all observations in the observation container
@@ -354,7 +361,7 @@ class cssens(ctools.cscript):
             # Fit test source to the simulated events in the observation
             # container
             fit = ctools.ctlike(sim)
-            fit['edisp']   = edisp=self._edisp
+            fit['edisp']   = self['edisp'].boolean()
             fit['debug']   = self['debug'].boolean()
             fit['chatter'] = self['chatter'].integer()
             fit.run()
@@ -383,29 +390,24 @@ class cssens(ctools.cscript):
             # Assess quality based on a comparison between Npred and Nevents
             quality = npred - nevents
  
-            # Write test source fit results into logger
-            self._log_header3(gammalib.EXPLICIT, 'Test source model fit')
-            self._log_value(gammalib.EXPLICIT, 'Test statistics', ts)
-            self._log_value(gammalib.EXPLICIT, 'log likelihood', logL)
-            self._log_value(gammalib.EXPLICIT, 'Number of predicted events', quality)
-            self._log_value(gammalib.EXPLICIT, 'Fit quality', npred)
-            for model in models:
-                self._log_value(gammalib.EXPLICIT, 'Model', model.name())
-                for par in model:
-                    self._log_string(gammalib.EXPLICIT, str(par))
-
-            # If TS was non-positive then increase the test flux and start
-            # over
+            # Write fit results into logger
+            name  = 'Iteration %d' % iterations
+            value = ('TS=%10.4f  Sim=%9.4f mCrab  Fit=%9.4f mCrab  '
+                     'Sens=%e erg/cm2/s' %
+                     (ts, test_crab_flux*1000.0, crab_flux*1000.0, sensitivity))
+            self._log_value(gammalib.TERSE, name, value)
+
+            # If TS was non-positive then increase the test flux and start over
             if ts <= 0.0:
  
                 # If the number of iterations was exceeded then stop
-                if (iterations >= self._max_iter):
+                if (iterations >= max_iter):
                     self._log_string(gammalib.TERSE,
-                         ' Test ended after %d iterations.' % self._max_iter)
+                         ' Test ended after %d iterations.' % max_iter)
                     break
  
-                # Increase test flux by a factor of 2
-                test_crab_flux = test_crab_flux * 2.0
+                # Increase test flux
+                test_crab_flux *= 3.0
  
                 # Signal start we start over
                 self._log_string(gammalib.EXPLICIT,
@@ -414,95 +416,64 @@ class cssens(ctools.cscript):
                 # ... and start over
                 continue
  
-            # Compute flux correction factor based on average TS. We
-            # assume here that the significance is proportional to the source
-            # flux, and that the significance scales as the square root of
-            # the Test Statistic
-            correct = math.sqrt(self._ts_thres/ts)
+            # Append result entry to result list
+            result = {'ts': ts, 'crab_flux': crab_flux,
+                                'photon_flux': photon_flux,
+                                'energy_flux': energy_flux}
+            results.append(result)
+
+            # Predict Crab flux at threshold TS using a linear regression of
+            # the log(TS) and log(crab_flux) values that have so far been
+            # computed. If not enough results are available than use a simple
+            # TS scaling relation.
+            if len(results) > 1:
+                pred_crab_flux, regcoeff = self._predict_flux(results, ts_thres)
+                correct                  = pred_crab_flux / crab_flux
+            else:
+                correct = math.sqrt(ts_thres/ts)
  
             # Compute extrapolated fluxes based on the flux correction factor
             crab_flux   = correct * crab_flux
             photon_flux = correct * photon_flux
             energy_flux = correct * energy_flux
             sensitivity = correct * sensitivity
-            crab_flux_value.append(crab_flux)
-            photon_flux_value.append(photon_flux)
-            energy_flux_value.append(energy_flux)
-            sensitivity_value.append(sensitivity)
-
-            # Write fit results
-            if self._logExplicit():
-                self._log_value(gammalib.EXPLICIT, 'Photon flux',
-                                str(photon_flux)+' ph/cm2/s')
-                self._log_value(gammalib.EXPLICIT, 'Energy flux',
-                                str(energy_flux)+' erg/cm2/s')
-                self._log_value(gammalib.EXPLICIT, 'Crab flux',
-                                str(crab_flux*1000.0)+' mCrab')
-                self._log_value(gammalib.EXPLICIT, 'Differential sensitivity',
-                                str(sensitivity)+' erg/cm2/s')
-                for model in models:
-                    self._log_value(gammalib.EXPLICIT, 'Model', model.name())
-                    for par in model:
-                        self._log_string(gammalib.EXPLICIT, str(par))
-            elif self._logTerse():
-                name  = 'Iteration %d' % iterations
-                value = ('TS=%.4f corr=%e  %e ph/cm2/s = %e erg/cm2/s = %.2f '
-                         'mCrab = %e erg/cm2/s' % (ts, correct, photon_flux,
-                         energy_flux, crab_flux*1000.0, sensitivity))
-                self._log_value(gammalib.TERSE, name, value)
-
-            # Compute sliding average of extrapolated fitted prefactor,
-            # photon and energy flux. This damps out fluctuations and
-            # improves convergence.
-            crab_flux   = 0.0
-            photon_flux = 0.0
-            energy_flux = 0.0
-            sensitivity = 0.0
-            num         = 0.0
-            for k in range(self._num_avg):
-                inx = len(crab_flux_value) - k - 1
-                if inx >= 0:
-                    crab_flux   += crab_flux_value[inx]
-                    photon_flux += photon_flux_value[inx]
-                    energy_flux += energy_flux_value[inx]
-                    sensitivity += sensitivity_value[inx]
-                    num      += 1.0
-            crab_flux   /= num
-            photon_flux /= num
-            energy_flux /= num
-            sensitivity /= num
-
-            # Compare average flux to last average
-            if iterations > self._num_avg:
+
+            # If we have at least 3 results then check if the flux determination
+            # at the TS threshold has converged
+            if len(results) > 3:
                 if test_crab_flux > 0:
+
+                    # Compute fractional change in the Crab flux between two
+                    # iterations
                     ratio = crab_flux/test_crab_flux
  
-                    # We have 2 convergence criteria:
-                    # 1. The average flux does not change
-                    # 2. The flux correction factor is small
-                    if ratio   >= 0.99 and ratio   <= 1.01 and \
-                       correct >= 0.9  and correct <= 1.1:
-                        if self._logTerse():
-                            self._log(' Converged ('+str(ratio)+')\n')
+                    # If fractional change is smaller than the required position
+                    # the iterations are stopped
+                    if ratio > 1.0-ratio_precision and \
+                       ratio < 1.0+ratio_precision:
+                        value = ('TS=%10.4f  Sim=%9.4f mCrab                  '
+                                 '     Sens=%e erg/cm2/s' %
+                                 (ts, crab_flux*1000.0, sensitivity))
+                        self._log_value(gammalib.TERSE, 'Converged result', value)
+                        self._log_value(gammalib.TERSE, 'Converged flux ratio', ratio)
+                        self._log_value(gammalib.TERSE, 'Regression coefficient',
+                                        regcoeff)
                         break
                 else:
-                    if self._logTerse():
-                        self._log(' Flux is zero.\n')
+                    self._log_value(gammalib.TERSE, 'Not converged', 'Flux is zero')
                     break
  
-            # Use average for next iteration
+            # Set test flux for next iteration
             test_crab_flux = crab_flux
  
             # Exit loop if number of trials exhausted
-            if (iterations >= self._max_iter):
-                if self._logTerse():
-                    self._log(' Test ended after '+str(self._max_iter)+
-                              ' iterations.\n')
+            if (iterations >= max_iter):
+                self._log_string(gammalib.TERSE,
+                                 ' Test ended after %d iterations.' % max_iter)
                 break
  
         # Write fit results into logger
         self._log_header3(gammalib.TERSE, 'Fit results')
-        self._log_value(gammalib.TERSE, 'Test statistics', ts)
         self._log_value(gammalib.TERSE, 'Photon flux',
                         str(photon_flux)+' ph/cm2/s')
         self._log_value(gammalib.TERSE, 'Energy flux',
@@ -523,16 +494,84 @@ class cssens(ctools.cscript):
         # Restore energy boundaries of observation container
         for i, obs in enumerate(self._obs):
             obs.events().ebounds(self._obs_ebounds[i])
-        
+
         # Store result
         result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
                   'crab_flux': crab_flux, 'photon_flux': photon_flux, \
                   'energy_flux': energy_flux, \
-                  'sensitivity': sensitivity}
+                  'sensitivity': sensitivity, 'regcoeff': regcoeff, \
+                  'nevents': nevents, 'npred': npred}
  
         # Return result
         return result
  
+    def _predict_flux(self, results, ts):
+        """
+        Predict Crab flux for a given TS value
+
+        The Crab flux at a given Test Statistic value is predicted by doing a
+        linear regression of the log(TS) and log(crab_flux) values in a results
+        list.
+
+        See https://en.wikipedia.org/wiki/Simple_linear_regression
+
+        Parameters
+        ----------
+        results : list of dict
+            List of results
+        ts : float
+            Test Statistic value
+
+        Returns
+        -------
+        crab_flux_prediction, regcoeff : tuple of float
+            Predicted Crab flux and regression coefficient
+        """
+        # Compute means and regression coefficient
+        mean_x  = 0.0
+        mean_y  = 0.0
+        mean_xy = 0.0
+        mean_xx = 0.0
+        mean_yy = 0.0
+        for result in results:
+            x        = math.log(float(result['ts']))
+            y        = math.log(float(result['crab_flux']))
+            mean_x  += x
+            mean_y  += y
+            mean_xy += x * y
+            mean_xx += x * x
+            mean_yy += y * y
+        norm     = 1.0 / float(len(results))
+        mean_x  *= norm
+        mean_y  *= norm
+        mean_xy *= norm
+        mean_xx *= norm
+        mean_yy *= norm
+        rxy      = (mean_xy - mean_x * mean_y) / \
+                   math.sqrt((mean_xx - mean_x*mean_x) *
+                             (mean_yy - mean_y*mean_y))
+        regcoeff = rxy*rxy
+
+        # Compute regression line slope
+        beta_nom   = 0.0
+        beta_denom = 0.0
+        for result in results:
+            x           = math.log(float(result['ts']))
+            y           = math.log(float(result['crab_flux']))
+            beta_nom   += (x - mean_x) * (y - mean_y)
+            beta_denom += (x - mean_x) * (x - mean_x)
+        beta = beta_nom / beta_denom
+
+        # Compute regression line offset
+        alpha = mean_y - beta * mean_x
+
+        # Predict Crab flux at TS threshold
+        log_ts_thres         = math.log(ts)
+        crab_flux_prediction = math.exp(alpha + beta * log_ts_thres)
+
+        # Return
+        return (crab_flux_prediction, regcoeff)
+
  
     # Public methods
     def run(self):
@@ -545,15 +584,16 @@ class cssens(ctools.cscript):
  
         # Get parameters
         self._get_parameters()
-        
+
         # Loop over observations and store a deep copy of the energy
         # boundaries for later use
         for obs in self._obs:
             self._obs_ebounds.append(obs.events().ebounds().copy())
-        
+
         # Initialise script
         colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux',
-                    'energy_flux', 'sensitivity']
+                    'energy_flux', 'sensitivity', 'regcoeff', 'nevents',
+                    'npred']
         results  = []
  
         # Set test source model for this observation
@@ -566,30 +606,35 @@ class cssens(ctools.cscript):
         # Write models into logger
         self._log_models(gammalib.NORMAL, models, 'Model')
  
+        # Get sensitivity type
+        sensitivity_type = self['type'].string()
+
         # Write header
         self._log_header1(gammalib.TERSE, 'Sensitivity determination')
-        self._log_value(gammalib.TERSE, 'Type', str(self._type))
+        self._log_value(gammalib.TERSE, 'Type', sensitivity_type)
  
         # Loop over energy bins
         for ieng in range(self._ebounds.size()):
  
             # Set energies
-            if self._type == 'Differential':
+            if sensitivity_type == 'Differential':
                 emin  = self._ebounds.emin(ieng)
                 emax  = self._ebounds.emax(ieng)
-            elif self._type == 'Integral':
+            elif sensitivity_type == 'Integral':
                 emin  = self._ebounds.emin(ieng)
                 emax  = self._ebounds.emax()
             else:
                 msg = ('Invalid sensitivity type "%s" encountered. Either '
-                       'specify "Differential" or "Integral".' % self._type)
+                       'specify "Differential" or "Integral".' %
+                       sensitivity_type)
                 raise RuntimeError(msg)
  
             # Determine sensitivity
             result = self._get_sensitivity(emin, emax, models)
  
             # Write out trial result
-            ioutils.write_csv_row(self._outfile.url(), ieng, colnames, result)
+            ioutils.write_csv_row(self['outfile'].filename().url(), ieng,
+                                  colnames, result)
  
             # Append results
             results.append(result)
@@ -100,9 +100,6 @@ General parameters
 ``(max_iter = 50) [integer]``
     Maximum number of iterations.
  
-``(num_avg = 3) [integer]``
-    Number of iterations for sliding average.
-
  
 Standard parameters
 -------------------
@@ -3,21 +3,25 @@
 ctbutterfly
 ===========
  
-Computes butterfly diagram for a power law model.
+Computes butterfly diagram for a given spectral model.
  
  
 Synopsis
 --------
  
-This tool calculates a butterfly diagram for a specific source with power law
-spectral model. The butterfly diagram is the envelope of all power law models
-that are within a given confidence limit compatible with the data. By default
-a confidence level of 68% is used, but the level can be adjusted using the 
-hidden ``confidence`` parameter. ctbutterfly computes this envelope by
-evaluating for each energy the minimum and maximum intensity of all power law
-models that fall within the error ellipse of the prefactor and index parameters.
-The error ellipse is derived from the covariance matrix of a maximum likelihood
-fit.
+This tool calculates a butterfly diagram for a specific source according to 
+its spectral model. The butterfly diagram is the envelope of all spectral
+models that are within a given confidence limit compatible with the data.
+The default method used for the calculation is Gaussian error propagation
+using the covariance matrix from a maximum likelihood fit. By default a
+confidence level of 68% is used, but the level can be adjusted using the
+hidden ``confidence`` parameter. For power law models, an alternative
+calculation method can be specified by setting the hidden parameter
+``method=ENVELOPE``. By using this method, the envelope is computed by
+evaluating for each energy the minimum and maximum intensity of all power
+law models that fall within the error ellipse of the prefactor and index
+parameters. The error ellipse is derived from the covariance matrix of a
+maximum likelihood fit.
  
 ctbutterfly assumes that the input model (parameter ``inmodel``) has been 
 adjusted using :doc:`ctlike` to the data, but if this is not the case you 
@@ -68,11 +72,14 @@ General parameters
  
 ``outfile [file]``
     Output butterfly ASCII file name.
-
+    
 ``(fit = no) [boolean]``
     Performs maximum likelihood fitting of input model ignoring any provided
     covariance matrix.
- 	 	 
+
+``(method = "GAUSSIAN") <GAUSSIAN|ENVELOPE> [string]``
+    Computation method.
+    
 ``(confidence = 0.68) [real]``
     Confidence level for error computation.
  
@@ -48,13 +48,13 @@ the spatial model parameters can be kept fixed during that computation
 ctlike generates an output model XML file that contains the values of the 
 best fitting model parameters. For all free parameters, an ``error`` attribute
 is added that provides the statistical uncertainty in the parameter estimate.
-In addition, the entire output of the covariance matrix to a separate FITS file
-can be queried (parameter ``outcovmat``). The parameter order given in the file
-corresponds to the arrangement of the covariance matrix entries. If computation
-of the Test Statistics has been requested for a model component, a ``ts``
-attribute providing the Test Statistics value is added. The output model can be
-used as an input model for other ctools.
+If computation of the Test Statistics has been requested for a model component,
+a ``ts`` attribute providing the Test Statistics value is added. The output model
+can be used as an input model for other ctools.
  
+In addition, the covariance matrix can be saved into a file by specifying a
+file name for the ``outcovmat`` parameter. If the file name ends with ``.fits``
+or ``.fit`` a FITS file is written, otherwise a CSV file is written.
  
  
 General parameters
@@ -89,7 +89,7 @@ General parameters
     the maximum likelihood fit.
  
 ``(outcovmat = NONE) [string]``
-    Output FITS file to store covariance matrix.
+    Output FITS or CSV file to store covariance matrix.
  
 ``(stat = POISSON) [string]``
     Fitting statistics (POISSON or GAUSSIAN; only affects binned analysis).
+.. _ctprob:
+
+ctprob
+========
+
+Compute event probability for a given model.
+
+
+Synopsis
+--------
+
+This tool takes an event list file (or an XML observation definition file) and
+a model as input and computes for each event the probability that it originates
+from a specific model component. The probability is computed by by evaluating
+the differential event probability for a given model and normalizing these
+probabilities so that the sum of the probabilities for all model components is
+unity for each event.
+
+The tool appends for each model component a single precision column to the event
+list that contains the event probabilities. The name of the column is build from
+the model name prefixed with ``PROB_``.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+    Input event list or observation definition XML file
+ 	 	 
+``outobs [file]``
+    Output event list or observation definition XML file.
+ 	 	 
+``(prefix = "prob_") [string]``
+    Prefix for output event lists in observation definition file.
+ 	 	 
+``inmodel [string]``
+    Input model definition XML file.
+
+``caldb [string]``
+    Calibration database.
+ 	 	 
+``irf [string]``
+    Instrument response function.
+ 	 	 
+``(edisp = no) [boolean]``
+    Applies energy dispersion to response computation.
+
+
+Standard parameters
+-------------------
+
+``(publish = no) [boolean]``
+    Specifies whether the event list(s) should be published on VO Hub.
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     chatter = 0: no information will be logged
+     
+     chatter = 1: only errors will be logged
+     
+     chatter = 2: errors and actions will be logged
+     
+     chatter = 3: report about the task execution
+     
+     chatter = 4: detailed report about the task execution
+ 	 	 
+``(clobber = yes) [boolean]``
+    Specifies whether existing output files should be overwritten.
+ 	 	 
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+ 	 	 
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = ctprob.log) [string]``
+    Name of log file.
+
+
+Related tools or scripts
+------------------------
+
+:doc:`ctphase`
+:doc:`csphasecrv`
+
@@ -27,6 +27,7 @@ ctools
    ctmodel --- Computes model counts cube <ctmodel>
    ctobssim --- Simulate observations <ctobssim>
    ctpsfcube --- Generates point spread function cube <ctpsfcube>
+   ctprob --- Computes event probability for a given model <ctprob>
    ctselect --- Selects event data <ctselect>
    ctskymap --- Generates sky map <ctskymap>
    cttsmap --- Generates Test Statistic map <cttsmap>
@@ -59,6 +59,7 @@ testscripts = ctools/tests/__init__.py \
               $(top_srcdir)/test/test_ctbutterfly.py \
               $(top_srcdir)/test/test_ctulimit.py \
               $(top_srcdir)/test/test_cterror.py \
+              $(top_srcdir)/test/test_ctprob.py \
               $(top_srcdir)/test/test_pipelines.py
  
 # Test data to be distributed in $(pyexecdir)/ctools/tests/data
 /***************************************************************************
  *                        ctool - ctool base class                         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2014-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -105,6 +105,8 @@ public:
     %rename(_insert_energy_boundaries) insert_energy_boundaries;
     %rename(_cube_layer_usage)         cube_layer_usage;
     %rename(_is_valid_filename)        is_valid_filename;
+    %rename(_get_gtiname)              get_gtiname;
+    %rename(_save_event_list)          save_event_list;
     %rename(_warn_too_few_energies)    warn_too_few_energies;
     %rename(_warn_xml_suffix)          warn_xml_suffix;
  
@@ -147,6 +149,13 @@ public:
     std::vector<bool> cube_layer_usage(const GEbounds& cube_ebounds,
                                        const GEbounds& list_ebounds) const;
     bool              is_valid_filename(const GFilename& filename) const;
+    std::string       get_gtiname(const std::string& filename,
+                                  const std::string& evtname) const;
+    void              save_event_list(const GCTAObservation* obs,
+                                      const std::string&     infile,
+                                      const std::string&     evtname,
+                                      const std::string&     gtiname,
+                                      const std::string&     outfile) const;
  
     // Protected warning strings
     std::string     warn_too_few_energies(const GEnergies& energies) const;
@@ -73,6 +73,7 @@ namespace std {
 %include "ctbutterfly.i"
 %include "ctulimit.i"
 %include "cterror.i"
+%include "ctprob.i"
  
 /* __ Test function ______________________________________________________ */
 %pythoncode %{
 /***************************************************************************
  *          ctphase - Append phase information to CTA events file          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Leonardo Di Venere                               *
+ *  copyright (C) 2017 by Joshua Cardenzana                                *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file ctphase.i
  * @brief Append phase information to CTA events file
- * @author Leonardo Di Venere
+ * @author Joshua Cardenzana
  */
 %{
 /* Put headers and other declarations here that are needed for compilation */
+/***************************************************************************
+ *           ctprob - Computes probability for a given model               *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2012-2016 by Leonardo Di Venere                          *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctprob.i
+ * @brief Computes probability for a given model
+ * @author Leonardo Di Venere
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctprob.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctprob
+ *
+ * @brief Data selection tool
+ ***************************************************************************/
+class ctprob : public ctobservation {
+public:
+    // Constructors and destructors
+    ctprob(void);
+    explicit ctprob(const GObservations& obs);
+    ctprob(int argc, char *argv[]);
+    ctprob(const ctprob& app);
+    virtual ~ctprob(void);
+
+    // Methods
+    void clear(void);
+    void run(void);
+    void save(void);
+    void publish(const std::string& name = "");
+};
+
+
+/***********************************************************************//**
+ * @brief Data selection tool Python extension
+ ***************************************************************************/
+%extend ctprob {
+    ctprob copy() {
+        return (*self);
+    }
+}
 #############################################################################
 # Makefile.am for ctools module                                             #
 # ------------------------------------------------------------------------- #
-# Copyright (C) 2016 Juergen Knoedlseder                                    #
+# Copyright (C) 2016-2017 Juergen Knoedlseder                               #
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  This program is free software: you can redistribute it and/or modify     #
@@ -41,7 +41,8 @@ SUBDIRS = support \
           ctbutterfly \
           ctulimit \
           cterror \
-          ctphase
+          ctphase \
+          ctprob
  
 # Scripts to be distributed
 scripts = ctools-setup \
@@ -70,7 +71,8 @@ libctools_la_LIBADD = support/libsupport.la \
                       ctbutterfly/libctbutterfly.la \
                       ctulimit/libctulimit.la \
                       cterror/libcterror.la \
-                      ctphase/libctphase.la
+                      ctphase/libctphase.la \
+                      ctprob/libctprob.la
  
 libctools_la_LDFLAGS = -version-info $(CTOOLS_LT_VERSION) @OPENMP_LIB_LDFLAGS@
  
@@ -55,7 +55,7 @@
  *
  * Constructs an empty event binning tool.
  ***************************************************************************/
-ctbin::ctbin(void) : ctobservation(CTBIN_NAME, CTBIN_VERSION)
+ctbin::ctbin(void) : ctobservation(CTBIN_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -73,7 +73,7 @@ ctbin::ctbin(void) : ctobservation(CTBIN_NAME, CTBIN_VERSION)
  * Constructs event binning tool from an observation container.
  ***************************************************************************/
 ctbin::ctbin(const GObservations& obs) :
-       ctobservation(CTBIN_NAME, CTBIN_VERSION, obs)
+       ctobservation(CTBIN_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -94,7 +94,7 @@ ctbin::ctbin(const GObservations&amp; obs) :
  * parameter setting.
  ***************************************************************************/
 ctbin::ctbin(int argc, char *argv[]) : 
-       ctobservation(CTBIN_NAME, CTBIN_VERSION, argc, argv)
+       ctobservation(CTBIN_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTBIN_NAME    "ctbin"
-#define CTBIN_VERSION "1.2.0"
+#define CTBIN_NAME "ctbin"
  
  
 /***********************************************************************//**
@@ -50,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
+ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -69,7 +69,7 @@ ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
  * observations container.
  ***************************************************************************/
 ctbkgcube::ctbkgcube(const GObservations& obs) :
-		ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, obs)
+		ctobservation(CTBKGCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -87,7 +87,7 @@ ctbkgcube::ctbkgcube(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctbkgcube::ctbkgcube(int argc, char *argv[]) :
-		ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, argc, argv)
+		ctobservation(CTBKGCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTBKGCUBE_NAME    "ctbkgcube"
-#define CTBKGCUBE_VERSION "1.1.0"
+#define CTBKGCUBE_NAME "ctbkgcube"
  
  
 /***********************************************************************//**
 /***************************************************************************
  *                 ctbutterfly - butterfly calculation tool                *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2016 by Michael Mayer                               *
+ *  copyright (C) 2014-2017 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -52,7 +52,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
+ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -71,7 +71,7 @@ ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERS
  * observations container.
  ***************************************************************************/
 ctbutterfly::ctbutterfly(const GObservations& obs) :
-             ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, obs)
+             ctlikelihood(CTBUTTERFLY_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ ctbutterfly::ctbutterfly(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctbutterfly::ctbutterfly(int argc, char *argv[]) :
-             ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
+             ctlikelihood(CTBUTTERFLY_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -243,11 +243,14 @@ void ctbutterfly::run(void)
         m_obs.optimize(m_opt);
         m_obs.errors(m_opt);
  
-        // Get covariance matrix by inverting the curvature matrix. Note that
-        // the curvature matrix is for the unscaled parameters, hence we have
-        // to scale them before we can proceed.
+        // Get covariance matrix.
         GObservations::likelihood likelihood = m_obs.function();
-        m_covariance = likelihood.curvature()->invert();
+        if (m_method == "GAUSSIAN") {
+            m_covariance = likelihood.curvature()->invert();
+        }
+        else {   
+            m_covariance = likelihood.covariance();
+        }
  
         // Write optimizer, optimised model and covariance matrix into logger
         log_string(NORMAL, m_opt.print(m_chatter));
@@ -271,157 +274,26 @@ void ctbutterfly::run(void)
         GObservations::likelihood likelihood = m_obs.function();
         likelihood.eval(pars);
  
-        // Get covariance matrix by inverting the curvature matrix. Note that
-        // the curvature matrix is for the unscaled parameters, hence we have
-        // to scale them before we can proceed.
-        m_covariance = likelihood.curvature()->invert();
+        // Get covariance matrix.
+        if (m_method == "GAUSSIAN") {
+            m_covariance = likelihood.curvature()->invert();
+        }
+        else {   
+            m_covariance = likelihood.covariance();
+        }
  
         // Write covariance matrix into logger
         log_string(EXPLICIT, m_covariance.print(m_chatter));
  
     } // endif: covariance computation needed
  
-    // Write header into logger
-    log_header1(TERSE, "Prepare butterfly computation");
-
-    // Find parameter indices. We do this by the memory location of the
-    // model parameters.
-    GModelSky* skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
-    GModelPar* par_prefactor = &(*skymodel)["Prefactor"];
-    GModelPar* par_index     = &(*skymodel)["Index"];
-    int        inx_prefactor = -1;
-    int        inx_index     = -1;
-    for (int i = 0; i < models.npars(); ++i) {
-        if (models.pars()[i] == par_prefactor) {
-            inx_prefactor = i;
-        }
-        if (models.pars()[i] == par_index) {
-            inx_index = i;
-        }
+    // Compute butterfly depending on selected method
+    if (m_method == "GAUSSIAN") {
+        gaussian_error_propagation(models);
+    }
+    else {
+        ellipsoid_boundary(models);
     }
-
-    // Store model values and scales
-    double  prefactor_mean  = (*skymodel)["Prefactor"].value();
-    double  index_mean      = (*skymodel)["Index"].value();
-    double  pivot_mean      = (*skymodel)["PivotEnergy"].value();
-    double  prefactor_scale = (*skymodel)["Prefactor"].scale();
-    double  index_scale     = (*skymodel)["Index"].scale();
-    GEnergy pivot(pivot_mean, "MeV");
-
-    // Write parameter indices into logger
-    log_value(NORMAL, "Prefactor", gammalib::str(prefactor_mean)+
-                      " ph/cm2/s/MeV ("+gammalib::str(inx_prefactor)+")");
-    log_value(NORMAL, "Index", gammalib::str(index_mean)+
-                      " ("+gammalib::str(inx_index)+")");
-    log_value(NORMAL, "Pivot energy", pivot.print());
-
-    // Extract covariance elements
-    double cov_pp = m_covariance(inx_prefactor,inx_prefactor);
-    double cov_pg = m_covariance(inx_prefactor,inx_index);
-    double cov_gg = m_covariance(inx_index,inx_index);
-
-    // Compute eigenvectors and eigenvalues
-    double  lambda1;
-    double  lambda2;
-    GVector vector1(2);
-    GVector vector2(2);
-    eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg,
-                 &lambda1, &lambda2, &vector1, &vector2);
-
-    // Write eigenvectors and eigenvalues into logger
-    log_value(NORMAL, "Eigenvalue 1", lambda1);
-    log_value(NORMAL, "Eigenvalue 2", lambda2);
-    log_value(NORMAL, "Eigenvector 1", vector1.print());
-    log_value(NORMAL, "Eigenvector 2", vector2.print());
-
-    // Confidence scaling
-    double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
-    double scale = (sigma*sigma);
-
-    // Write confidence level information into logger
-    log_value(NORMAL, "Confidence level", m_confidence);
-    log_value(NORMAL, "Corresponding scaling", scale);
-
-    // Compute minor and major axes of error ellipse
-    double major = scale*std::sqrt(lambda1);
-    double minor = scale*std::sqrt(lambda2);
-
-    // Write header into logger
-    log_header1(TERSE, "Generate butterfly");
-
-    // Initialise result arrays
-    m_energies.assign(m_ebounds.size(), 0.0);
-    m_intensities.assign(m_ebounds.size(), 0.0);
-    m_min_intensities.assign(m_ebounds.size(), 1e30);
-    m_max_intensities.assign(m_ebounds.size(), 0.0);
-
-    // Walk around the error ellipse
-    int    nt = 360;
-    double dt = gammalib::twopi/double(nt);
-    double t  = 0.0;
-    for (int i = 0; i < nt; ++i, t += dt) {
-
-        // Pre-compute sin and cosine
-        double sin_t = std::sin(t);
-        double cos_t = std::cos(t);
-
-        // Compute prefactor and index
-        double prefactor = (major * cos_t * vector1[0] +
-                            minor * sin_t * vector2[0]) * prefactor_scale +
-                           prefactor_mean;
-        double index     = (major * cos_t * vector1[1] +
-                            minor * sin_t * vector2[1]) * index_scale +
-                           index_mean;
-
-        // Write prefactor and index into logger
-        log_value(EXPLICIT, "Angle "+gammalib::str(t*gammalib::rad2deg),
-                  gammalib::str(prefactor)+" ph/cm2/s/MeV; "+
-                  gammalib::str(index));
-
-        // Setup power law model
-        GModelSpectralPlaw plaw(prefactor, index, pivot);
-
-        // Evaluate power law at each energy and extract minimum and
-        // maximum intensity
-        for (int k = 0; k < m_ebounds.size(); ++k) {
-
-            // Get the energy of current bin
-            GTime   time;
-            GEnergy energy = m_ebounds.elogmean(k);
-
-            // Evaluate power law
-            double intensity = plaw.eval(energy, time);
-
-            // Get extremes
-            if (intensity < m_min_intensities[k]) {
-                m_min_intensities[k] = intensity;
-            }
-            if (intensity > m_max_intensities[k]) {
-                m_max_intensities[k] = intensity;
-            }
-
-        } // endfor: looped over all energies
-
-    } // endfor: looped over all angles
-
-    // Setup power law model for mean prefactor and index
-    GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot);
-
-    // Compute energy and intensity vectors for mean power law
-    for (int k = 0; k < m_ebounds.size(); ++k) {
-
-        // Get the energy of current bin
-        GTime   time;
-        GEnergy energy = m_ebounds.elogmean(k);
-
-        // Evaluate power law
-        double intensity = plaw.eval(energy, time);
-
-        // Store results
-        m_energies[k]    = energy.MeV();
-        m_intensities[k] = intensity;
-
-    } // endfor: looped over energies
  
     // Write header into logger
     log_header1(TERSE, "Results");
@@ -431,11 +303,11 @@ void ctbutterfly::run(void)
  
         // Write results into logger
         log_value(NORMAL, "Intensity at "+gammalib::str(m_energies[k])+" MeV",
-                  gammalib::str(m_intensities[k])+" ("+
-                  gammalib::str(m_min_intensities[k])+", "+
-                  gammalib::str(m_max_intensities[k])+") ph/cm2/s/MeV");
+                gammalib::str(m_intensities[k])+" ("+
+                gammalib::str(m_min_intensities[k])+", "+
+                gammalib::str(m_max_intensities[k])+") ph/cm2/s/MeV");
  
-    } // endfor: looped over energies
+    } //endfor: loop over energy bin  
  
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -504,6 +376,7 @@ void ctbutterfly::init_members(void)
 {
     // Initialise members
     m_srcname.clear();
+    m_method.clear();
     m_confidence  = 0.68;
     m_max_iter    = 50;
     m_apply_edisp = false;
@@ -537,6 +410,7 @@ void ctbutterfly::copy_members(const ctbutterfly&amp; app)
 {
     // Copy attributes
     m_srcname     = app.m_srcname;
+    m_method      = app.m_method;
     m_confidence  = app.m_confidence;
     m_max_iter    = app.m_max_iter;
     m_apply_edisp = app.m_apply_edisp;
@@ -544,7 +418,7 @@ void ctbutterfly::copy_members(const ctbutterfly&amp; app)
     m_chatter     = app.m_chatter;
     m_ebounds     = app.m_ebounds;
     m_outfile     = app.m_outfile;
- 
+
     // Copy protected members
     m_covariance      = app.m_covariance;
     m_energies        = app.m_energies;
@@ -583,16 +457,13 @@ void ctbutterfly::get_parameters(void)
 {
     // Setup observations from "inobs" parameter
     setup_observations(m_obs);
-    
+
     // Setup models from "inmodel" parameter
     setup_models(m_obs, (*this)["srcname"].string());
  
     // Get name of test source
     m_srcname = (*this)["srcname"].string();
  
-    // Check model name and type
-    check_model();
-
     // Create energy boundaries from user parameters
     m_ebounds = create_ebounds();
  
@@ -600,18 +471,22 @@ void ctbutterfly::get_parameters(void)
     std::string matrixfilename = (*this)["matrix"].filename();
     if (matrixfilename != "NONE") {
         std::string msg = "Loading of matrix from file not implemented yet. "
-                          "Use filename = \"NONE\" to induce a recomputation of "
-                          "the matrix internally.";
+                "Use filename = \"NONE\" to induce a recomputation of "
+                "the matrix internally.";
         throw GException::feature_not_implemented(G_GET_PARAMETERS, msg);
         // m_covariance.load(matrixfilename);
     }
  
     // Get remaining parameters
     m_apply_edisp = (*this)["edisp"].boolean();
+    m_method      = (*this)["method"].string();
     m_confidence  = (*this)["confidence"].real();
     m_fit         = (*this)["fit"].boolean();
     m_chatter     = static_cast<GChatter>((*this)["chatter"].integer());
  
+    // Check model name and type
+    check_model();
+
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
     if (read_ahead()) {
@@ -627,6 +502,279 @@ void ctbutterfly::get_parameters(void)
  
  
 /***********************************************************************//**
+ * @brief Compute butterfly using the ellipsoid boundary method
+ ***************************************************************************/
+void ctbutterfly::ellipsoid_boundary(GModels& models)
+{
+    // Write header into logger
+    log_header1(TERSE, "Prepare butterfly computation");
+
+    // Find parameter indices. We do this by the memory location of the
+    // model parameters.
+    GModelSky* skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
+    GModelPar* par_prefactor = &(*skymodel)["Prefactor"];
+    GModelPar* par_index     = &(*skymodel)["Index"];
+    int        inx_prefactor = -1;
+    int        inx_index     = -1;
+    for (int i = 0; i < models.npars(); ++i) {
+        if (models.pars()[i] == par_prefactor) {
+            inx_prefactor = i;
+        }
+        if (models.pars()[i] == par_index) {
+            inx_index = i;
+        }
+    }
+
+    // Store model values
+    double  prefactor_mean  = (*skymodel)["Prefactor"].value();
+    double  index_mean      = (*skymodel)["Index"].value();
+    double  pivot_mean      = (*skymodel)["PivotEnergy"].value();
+    double  prefactor_scale = (*skymodel)["Prefactor"].scale();
+    double  index_scale     = (*skymodel)["Index"].scale();
+    GEnergy pivot(pivot_mean, "MeV");
+
+    // Write parameter indices into logger
+    log_value(NORMAL, "Prefactor", gammalib::str(prefactor_mean)+
+                                   " ph/cm2/s/MeV ("+
+                                   gammalib::str(inx_prefactor)+")");
+    log_value(NORMAL, "Index", gammalib::str(index_mean)+
+                               " ("+gammalib::str(inx_index)+")");
+    log_value(NORMAL, "Pivot energy", pivot.print());
+
+    // Extract covariance elements
+    double cov_pp = m_covariance(inx_prefactor,inx_prefactor);
+    double cov_pg = m_covariance(inx_prefactor,inx_index);
+    double cov_gg = m_covariance(inx_index,inx_index);
+
+    // Compute eigenvectors and eigenvalues
+    double  lambda1;
+    double  lambda2;
+    GVector vector1(2);
+    GVector vector2(2);
+    eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg,
+                 &lambda1, &lambda2, &vector1, &vector2);
+
+    // Write eigenvectors and eigenvalues into logger
+    log_value(NORMAL, "Eigenvalue 1", lambda1);
+    log_value(NORMAL, "Eigenvalue 2", lambda2);
+    log_value(NORMAL, "Eigenvector 1", vector1.print());
+    log_value(NORMAL, "Eigenvector 2", vector2.print());
+
+    // Confidence scaling, calculated from a Chi-Sqr. 
+    // distribution with 2 dof
+    double scale = std::sqrt( -2. * std::log(1.-m_confidence) );
+
+    // Write confidence level information into logger
+    log_value(NORMAL, "Confidence level", m_confidence);
+    log_value(NORMAL, "Corresponding scaling", scale);
+
+    // Compute minor and major axes of error ellipse
+    double major = scale*std::sqrt(lambda1);
+    double minor = scale*std::sqrt(lambda2);
+
+    // Write header into logger
+    log_header1(TERSE, "Generate butterfly");
+
+    // Initialise result arrays
+    m_energies.assign(m_ebounds.size(), 0.0);
+    m_intensities.assign(m_ebounds.size(), 0.0);
+    m_min_intensities.assign(m_ebounds.size(), 1e30);
+    m_max_intensities.assign(m_ebounds.size(), 0.0);
+
+    // Walk around the error ellipse
+    int    nt = 360;
+    double dt = gammalib::twopi/double(nt);
+    double t  = 0.0;
+    for (int i = 0; i < nt; ++i, t += dt) {
+
+        // Pre-compute sin and cosine
+        double sin_t = std::sin(t);
+        double cos_t = std::cos(t);
+
+        // Compute prefactor and index
+        double prefactor = (major * cos_t * vector1[0] +
+                           minor * sin_t * vector2[0]) * prefactor_scale +
+                           prefactor_mean;
+        double index     = (major * cos_t * vector1[1] +
+                           minor * sin_t * vector2[1]) * index_scale +
+                           index_mean;
+
+        // Write prefactor and index into logger
+        log_value(EXPLICIT, "Angle "+gammalib::str(t*gammalib::rad2deg),
+                gammalib::str(prefactor)+" ph/cm2/s/MeV; "+
+                gammalib::str(index));
+
+        // Setup power law model
+        GModelSpectralPlaw plaw(prefactor, index, pivot);
+
+        // Evaluate power law at each energy and extract minimum and
+        // maximum intensity
+        for (int k = 0; k < m_ebounds.size(); ++k) {
+
+            // Get the energy of current bin
+            GTime   time;
+            GEnergy energy = m_ebounds.elogmean(k);
+
+            // Evaluate power law
+            double intensity = plaw.eval(energy, time);
+
+            // Get extremes
+            if (intensity < m_min_intensities[k]) {
+                m_min_intensities[k] = intensity;
+            }
+            if (intensity > m_max_intensities[k]) {
+                m_max_intensities[k] = intensity;
+            }
+
+        } // endfor: looped over all energies
+
+    } // endfor: looped over all angles
+
+    // Setup power law model for mean prefactor and index
+    GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot);
+
+    // Compute energy and intensity vectors for mean power law
+    for (int k = 0; k < m_ebounds.size(); ++k) {
+
+        // Get the energy of current bin
+        GTime   time;
+        GEnergy energy = m_ebounds.elogmean(k);
+
+        // Evaluate power law
+        double intensity = plaw.eval(energy, time);
+
+        // Store results
+        m_energies[k]    = energy.MeV();
+        m_intensities[k] = intensity;
+
+    } // endfor: looped over energies
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Compute butterfly using Gaussian error propagation
+ *
+ * Computes the butterfly diagram using Gaussian error propagation. This
+ * works for any kind of model.
+ ***************************************************************************/
+void ctbutterfly::gaussian_error_propagation(GModels& models)
+{
+    // Write header into logger
+    log_header1(TERSE, "Prepare butterfly computation");
+
+    // Confidence scaling
+    double scale = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
+    
+    // Write confidence level information into logger
+    log_value(NORMAL, "Confidence level", m_confidence);
+    log_value(NORMAL, "Corresponding scaling", scale);
+
+    // Write header into logger
+    log_header1(TERSE, "Generate butterfly");
+ 
+    // Initialise result arrays
+    m_energies.clear();
+    m_intensities.clear();
+    m_min_intensities.clear();
+    m_max_intensities.clear();
+   
+    // Initialise dummy time to evaluate spectral model
+    GTime time = GTime();
+
+    // Initialise vector of gradients
+    GVector grad = GVector(m_obs.models().npars());
+
+    // Loop over energy bins
+    for (int i = 0; i < m_ebounds.size(); ++i) {
+
+        // Initialise number of gradients
+        int num_gradient = 0;
+
+        // Get the mean logarithmic energy of current bin
+        GEnergy energy = m_ebounds.elogmean(i);
+
+        // Initialise model flux value
+        double model_flux = 0.0;
+
+        // Loop over models
+        for (int j = 0; j < models.size(); ++j) {
+
+            // Check wether model is a skymodel
+            GModelSky* skymodel = dynamic_cast<GModelSky*>(models[j]);
+
+            // Yes ...
+            if (skymodel != NULL) {
+
+                // Set flux value of source of interest, evaluate gradients
+                if (skymodel->name() == m_srcname) {
+
+                    // Get spectral component of model
+                    GModelSpectral *spectral = skymodel->spectral();
+
+                    // Get model flux
+                    model_flux = spectral->eval(energy, time, true);
+
+                    // Loop over all model parameters
+                    for (int k = 0; k < skymodel->size(); ++k) {
+
+                        // Set gradient if we deal with spectral component
+                        bool parIsSpectral = false;
+                        for (int l = 0; l < spectral->size(); ++l) {
+                            if ((&((*spectral)[l])) == (&((*skymodel)[k]))) {
+                                grad[num_gradient] = (*spectral)[l].gradient();
+                                num_gradient++;
+                                parIsSpectral = true;
+                                break;
+                            }
+                        } // endfor: looped over spectral parameters
+
+                        // If parameter is spatial or temporal, count up
+                        if (!parIsSpectral) {
+                            num_gradient++;
+                        }
+
+                    } // endfor: looped over all model parameters
+
+                } // endif: model was source of interest
+                else {
+                    num_gradient += skymodel->size();
+                }
+
+            } // endif: model was sky model
+
+            // Skip other models (e.g. GModelData instances)
+            else {
+                num_gradient += models[j]->size();
+            }
+
+        } // endfor: Looped over models
+
+        // Multiply covariance to the gradient vector
+        GVector vector = m_covariance * grad;
+
+        // Get the error from the scalar product
+        double error = std::sqrt(grad * vector);
+
+        // Multiply with confidence scaling
+        error *= scale;
+
+        // Store flux, value and energy for saving
+        m_intensities.push_back(model_flux);
+        m_energies.push_back(energy.MeV());
+        m_min_intensities.push_back(model_flux-error);
+        m_max_intensities.push_back(model_flux+error);
+
+    } // endfor: Looped over energy bins
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Check if sky model is valid
  *
  * @exception GException::invalid_value
@@ -645,17 +793,22 @@ void ctbutterfly::check_model(void)
         GModelSky* model  = dynamic_cast<GModelSky*>(models[m_srcname]);
         if (model == NULL) {
             std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
-                              "Please specify the name of a sky model for "
-                              "butterfly computation.";
+                    "Please specify the name of a sky model for "
+                    "butterfly computation.";
             throw GException::invalid_value(G_CHECK_MODEL, msg);
         }
  
-        // Check that the spectral model is a power law
-        if (model->spectral()->type() != "PowerLaw") {
-            std::string msg = "\""+model->spectral()->type()+"\" cannot be "
-                              "used as spectral model for an butterfly "
-                              "computation. Please specify a power law model.";
-            throw GException::invalid_value(G_CHECK_MODEL, msg);
+        // Check that the spectral model is a power law for the "ENVELOPE"
+        // method
+        if (m_method == "ENVELOPE") {
+            if (model->spectral()->type() != "PowerLaw") {
+                std::string msg = "\""+model->spectral()->type()+"\" cannot be "
+                                  "used as spectral model for an butterfly "
+                                  "computation with method=ENVELOPE. Please "
+                                  "specify a power law model or switch to "
+                                  "method=GAUSSIAN.";
+                throw GException::invalid_value(G_CHECK_MODEL, msg);
+            }
         }
  
     } // endif: source model existed
@@ -663,8 +816,8 @@ void ctbutterfly::check_model(void)
     // Otherwise throw an exception
     else {
         std::string msg = "Source \""+m_srcname+"\" not found in model "
-                          "container. Please add a source with that name "
-                          "or check for possible typos.";
+                "container. Please add a source with that name "
+                "or check for possible typos.";
         throw GException::invalid_value(G_CHECK_MODEL, msg);
     }
  
@@ -689,13 +842,13 @@ void ctbutterfly::check_model(void)
  * http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html
  ***************************************************************************/
 void ctbutterfly::eigenvectors(const double& a,
-                               const double& b,
-                               const double& c,
-                               const double& d,
-                               double*       lambda1,
-                               double*       lambda2,
-                               GVector*      vector1,
-                               GVector*      vector2)
+        const double& b,
+        const double& c,
+        const double& d,
+        double*       lambda1,
+        double*       lambda2,
+        GVector*      vector1,
+        GVector*      vector2)
 {
     // Compute trace and determinant
     double trace = a + d;
@@ -730,7 +883,7 @@ void ctbutterfly::eigenvectors(const double&amp; a,
     // Normalize eigenvectors
     *vector1 /= norm(*vector1);
     *vector2 /= norm(*vector2);
-    
+
     // Return
     return;
 }
 /***************************************************************************
  *                ctbutterfly - butterfly calculation tool                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2016 by Michael Mayer                               *
+ *  copyright (C) 2014-2017 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -33,8 +33,7 @@
 #include "ctlikelihood.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTBUTTERFLY_NAME    "ctbutterfly"
-#define CTBUTTERFLY_VERSION "1.1.0"
+#define CTBUTTERFLY_NAME "ctbutterfly"
  
  
 /***********************************************************************//**
@@ -77,6 +76,10 @@ protected:
     void copy_members(const ctbutterfly& app);
     void free_members(void);
     void get_parameters(void);
+    void gaussian_error_propagation(GModels& models);
+    void ellipsoid_boundary(GModels& models);
+
+
     void check_model(void);
     void eigenvectors(const double& a,
                       const double& b,
@@ -89,6 +92,7 @@ protected:
  
     // User parameters
     std::string m_srcname;      //!< Name of source to compute butterfly
+    std::string m_method;       //!< Computation method
     double      m_confidence;   //!< Confidence level
     int         m_max_iter;     //!< Maximum number of iterations
     bool        m_apply_edisp;  //!< Apply energy dispersion?
@@ -27,6 +27,7 @@ caldb,      s, a, &quot;prod2&quot;,,, &quot;Calibration database&quot;
 irf,        s, a, "South_0.5h",,, "Instrument response function"
 edisp,      b, h, no,,, "Apply energy dispersion?"
 outfile,    f, a, "butterfly.txt",,, "Output ASCII file"
+method,     s, h, GAUSSIAN, GAUSSIAN|ENVELOPE,,"Butterfly computation method"
  
 #
 # Butterfly parameters
@@ -51,7 +51,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
+ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -70,7 +70,7 @@ ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION
  * an observation container.
  ***************************************************************************/
 ctcubemask::ctcubemask(const GObservations& obs) :
-            ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, obs)
+            ctobservation(CTCUBEMASK_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ ctcubemask::ctcubemask(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctcubemask::ctcubemask(int argc, char *argv[]) : 
-            ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv)
+            ctobservation(CTCUBEMASK_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -35,8 +35,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTCUBEMASK_NAME    "ctcubemask"
-#define CTCUBEMASK_VERSION "1.2.0"
+#define CTCUBEMASK_NAME "ctcubemask"
  
  
 /***********************************************************************//**
@@ -52,7 +52,7 @@
  *
  * Constructs an empty energy dispersion tool.
  ***************************************************************************/
-ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION)
+ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -71,7 +71,7 @@ ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VER
  * in an observation container @p obs.
  ***************************************************************************/
 ctedispcube::ctedispcube(const GObservations& obs) :
-             ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, obs)
+             ctobservation(CTEDISPCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -92,7 +92,7 @@ ctedispcube::ctedispcube(const GObservations&amp; obs) :
  * on the command line.
  ***************************************************************************/
 ctedispcube::ctedispcube(int argc, char *argv[]) :
-             ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, argc, argv)
+             ctobservation(CTEDISPCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTEDISPCUBE_NAME    "ctedispcube"
-#define CTEDISPCUBE_VERSION "1.1.0"
+#define CTEDISPCUBE_NAME "ctedispcube"
  
  
 /***********************************************************************//**
@@ -52,7 +52,7 @@
  *
  * Constructs an empty cterror tool.
  ***************************************************************************/
-cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION)
+cterror::cterror(void) : ctlikelihood(CTERROR_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -70,7 +70,7 @@ cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION)
  * Constructs cterror tool from an observations container.
  ***************************************************************************/
 cterror::cterror(const GObservations& obs) :
-         ctlikelihood(CTERROR_NAME, CTERROR_VERSION, obs)
+         ctlikelihood(CTERROR_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -90,7 +90,7 @@ cterror::cterror(const GObservations&amp; obs) :
  * parameters that are provided as command line arguments.
  ***************************************************************************/
 cterror::cterror(int argc, char *argv[]) :
-         ctlikelihood(CTERROR_NAME, CTERROR_VERSION, argc, argv)
+         ctlikelihood(CTERROR_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTERROR_NAME    "cterror"
-#define CTERROR_VERSION "1.2.0"
+#define CTERROR_NAME "cterror"
  
  
 /***********************************************************************//**
@@ -50,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION)
+ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -69,7 +69,7 @@ ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION)
  * observations container.
  ***************************************************************************/
 ctexpcube::ctexpcube(const GObservations& obs) :
-           ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, obs)
+           ctobservation(CTEXPCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -87,7 +87,7 @@ ctexpcube::ctexpcube(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctexpcube::ctexpcube(int argc, char *argv[]) :
-           ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, argc, argv)
+           ctobservation(CTEXPCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTEXPCUBE_NAME    "ctexpcube"
-#define CTEXPCUBE_VERSION "1.1.0"
+#define CTEXPCUBE_NAME "ctexpcube"
  
  
 /***********************************************************************//**
 /***************************************************************************
  *                ctlike - Maximum likelihood fitting tool                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -48,7 +48,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION)
+ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -66,7 +66,7 @@ ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION)
  * Constructs ctlike tool from an observations container.
  ***************************************************************************/
 ctlike::ctlike(const GObservations& obs) :
-        ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, obs)
+        ctlikelihood(CTLIKE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -86,7 +86,7 @@ ctlike::ctlike(const GObservations&amp; obs) :
  * that are provided as command line arguments.
  ***************************************************************************/
 ctlike::ctlike(int argc, char *argv[]) : 
-        ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, argc, argv)
+        ctlikelihood(CTLIKE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -321,8 +321,13 @@ void ctlike::run(void)
 /***********************************************************************//**
  * @brief Save results
  *
- * This method saves the fit results into a XML file and a FITS file. If
- * the filename parameters are "NONE", the files are not saved.
+ * This method saves the fit results and the covariance matrix. The fit
+ * results are written into a XML file while the covariance matrix is written
+ * into either a FITS or a CSV file, depending on the file type extension
+ * (an extension of `.fits` or `.fit` produce a FITS file, any other
+ * extension produces a CSV file).
+ *
+ * If the filenames are `NONE` no information is saved.
  ***************************************************************************/
 void ctlike::save(void)
 {
 /***************************************************************************
  *                ctlike - Maximum likelihood fitting tool                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTLIKE_NAME    "ctlike"
-#define CTLIKE_VERSION "1.2.0"
+#define CTLIKE_NAME "ctlike"
  
  
 /***********************************************************************//**
@@ -68,17 +67,17 @@ protected:
     double reoptimize_lm(void);
  
     // User parameters
-    bool          m_refit;           //!< Refitting
-    GFilename     m_outmodel;        //!< Source model output XML file name
-    GFilename     m_outcovmat;       //!< Covariance matrix output FITS file name
-    bool          m_apply_edisp;     //!< Apply energy dispersion?
-    bool          m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation?
-    GChatter      m_chatter;         //!< Chattiness
+    bool      m_refit;           //!< Refitting
+    GFilename m_outmodel;        //!< Source model output XML file name
+    GFilename m_outcovmat;       //!< Covariance matrix output file name
+    bool      m_apply_edisp;     //!< Apply energy dispersion?
+    bool      m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation?
+    GChatter  m_chatter;         //!< Chattiness
  
     // Members
-    int           m_max_iter;   //!< Maximum number of iterations
-    int           m_max_stall;  //!< Maximum number of stalls
-    double        m_logL;       //!< Maximum log likelihood
+    int       m_max_iter;        //!< Maximum number of iterations
+    int       m_max_stall;       //!< Maximum number of stalls
+    double    m_logL;            //!< Maximum log likelihood
 };
  
 #endif /* CTLIKE_HPP */
@@ -50,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION)
+ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -67,7 +67,7 @@ ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION)
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctmapcube::ctmapcube(int argc, char *argv[]) :
-           ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION, argc, argv)
+           ctool(CTMAPCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctool.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTMAPCUBE_NAME    "ctmapcube"
-#define CTMAPCUBE_VERSION "1.3.0"
+#define CTMAPCUBE_NAME "ctmapcube"
  
  
 /***********************************************************************//**
@@ -53,7 +53,7 @@ const GEnergy g_energy_margin(1.0e-12, &quot;TeV&quot;);
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION)
+ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -72,7 +72,7 @@ ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION)
  * observations container.
  ***************************************************************************/
 ctmodel::ctmodel(const GObservations& obs) :
-         ctobservation(CTMODEL_NAME, CTMODEL_VERSION, obs)
+         ctobservation(CTMODEL_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -90,7 +90,7 @@ ctmodel::ctmodel(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctmodel::ctmodel(int argc, char *argv[]) :
-         ctobservation(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
+         ctobservation(CTMODEL_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTMODEL_NAME    "ctmodel"
-#define CTMODEL_VERSION "1.1.0"
+#define CTMODEL_NAME "ctmodel"
  
  
 /***********************************************************************//**
@@ -67,7 +67,7 @@ const double g_roi_margin = 0.5;      //!&lt; Simulation radius margin (degrees)
  *
  * Constructs an empty ctobssim tool.
  ***************************************************************************/
-ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION)
+ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -85,7 +85,7 @@ ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION)
  * Constructs ctobssim tool from an observation container.
  ***************************************************************************/
 ctobssim::ctobssim(const GObservations& obs) :
-          ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION, obs)
+          ctobservation(CTOBSSIM_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -105,7 +105,7 @@ ctobssim::ctobssim(const GObservations&amp; obs) :
  * setting.
  ***************************************************************************/
 ctobssim::ctobssim(int argc, char *argv[]) :
-          ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
+          ctobservation(CTOBSSIM_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -1471,7 +1471,8 @@ void ctobssim::simulate_background(GCTAObservation* obs,
  
                 // Debug code: signal that we step into the response MC method
                 #if defined(G_BACKGROUND_DEBUG)
-                std::cout << "ctobssim::simulate_background: model->mc in" << std::endl;
+                std::cout << "ctobssim::simulate_background: model->mc in"
+                          << std::endl;
                 #endif
  
                 // Get simulated CTA event list. Note that this method
@@ -1481,7 +1482,8 @@ void ctobssim::simulate_background(GCTAObservation* obs,
  
                 // Debug code: signal that we came back
                 #if defined(G_BACKGROUND_DEBUG)
-                std::cout << "ctobssim::simulate_background: model->mc out" << std::endl;
+                std::cout << "ctobssim::simulate_background: model->mc out"
+                          << std::endl;
                 #endif
  
                 // Continue only if we got a CTA event list
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTOBSSIM_NAME    "ctobssim"
-#define CTOBSSIM_VERSION "1.2.0"
+#define CTOBSSIM_NAME "ctobssim"
  
  
 /***********************************************************************//**
 #############################################################################
 # Makefile.am for ctphase tool                                              #
 # ------------------------------------------------------------------------- #
-# Copyright (C) 2017 Leonardo Di Venere                                     #
+# Copyright (C) 2017 Joshua Cardenzana                                      #
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  This program is free software: you can redistribute it and/or modify     #
 /***************************************************************************
  *          ctphase - Append phase information to CTA events file          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Leonardo Di Venere                               *
+ *  copyright (C) 2017 by Joshua Cardenzana                                *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file ctphase.cpp
  * @brief Append phase information to CTA events file
- * @author Leonardo Di Venere
+ * @author Joshua Cardenzana
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -33,15 +33,12 @@
 #include "GTools.hpp"
  
 /* __ Method name definitions ____________________________________________ */
-#define G_RUN                                               "ctphase::run()"
-#define G_PHASE_EVENTS          "ctphase::phase_events(GCTAObservation*, "\
-                                  "std::string&, std::string&, std::string&)"
-#define G_SET_EBOUNDS    "ctphase::set_ebounds(GCTAObservation*, GEbounds&)"
+#define G_GET_PARAMETERS                          "ctphase::get_parameters()"
+#define G_PHASE_EVENTS              "ctphase::phase_events(GCTAObservation*)"
  
 /* __ Debug definitions __________________________________________________ */
  
 /* __ Coding definitions _________________________________________________ */
-//#define G_USE_MKSTEMP               //!< Use mkstemp for temporary filename
  
  
 /*==========================================================================
@@ -53,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, CTPHASE_VERSION)
+ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -72,7 +69,7 @@ ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, CTPHASE_VERSION)
  * provided in an observation container.
  ***************************************************************************/
 ctphase::ctphase(const GObservations& obs) :
-          ctobservation(CTPHASE_NAME, CTPHASE_VERSION, obs)
+         ctobservation(CTPHASE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -91,7 +88,7 @@ ctphase::ctphase(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctphase::ctphase(int argc, char *argv[]) : 
-          ctobservation(CTPHASE_NAME, CTPHASE_VERSION, argc, argv)
+         ctobservation(CTPHASE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -203,10 +200,11 @@ void ctphase::clear(void)
  
  
 /***********************************************************************//**
- * @brief Append phase information to event data
+ * @brief Runs the ctprob tool
  *
  * This method reads in the application parameters and loops over all
- * observations that were found to append phase information to the data.
+ * unbinned observations to compute the phase information for each event
+ * based on an input model.
  ***************************************************************************/
 void ctphase::run(void)
 {
@@ -222,118 +220,22 @@ void ctphase::run(void)
     log_observations(NORMAL, m_obs, "Input observation");
  
     // Write header into logger
-    log_header1(TERSE, "Appending phase informaiton");
+    log_header1(TERSE, "Compute event phase");
  
     // Initialise counters
     int n_observations = 0;
  
-    // Loop over all observation in the container
-    for (int i = 0; i < m_obs.size(); ++i) {
-
-        // Write header for the current observation
-        log_header3(TERSE, get_obs_header(m_obs[i]));
-
-        // Initialise event input and output filenames and the event
-        // and GTI extension names
-        m_infiles.push_back("");
-        m_evtname.push_back(gammalib::extname_cta_events);
-        m_gtiname.push_back(gammalib::extname_gti);
-
-        // Get CTA observation
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Skip observation if it's not CTA
-        if (obs == NULL) {
-            std::string msg = " Skipping "+m_obs[i]->instrument()+
-                              " observation";
-            log_string(NORMAL, msg);
-            continue;
-        }
-
-        // Skip observation if we have a binned observation
-        if (obs->eventtype() == "CountsCube") {
-            std::string msg = " Skipping binned "+m_obs[i]->instrument()+
-                              " observation";
-            log_string(NORMAL, msg);
-            continue;
-        }
+    // Loop over all unbinned CTA observations in the container
+    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+         obs = next_unbinned_observation()) {
  
         // Increment counter
         n_observations++;
  
-        // Save event file name (for possible saving)
-        m_infiles[i] = obs->eventfile();
-
-        // Extract event and GTI extension names from input FITS file
-        GFilename fname(m_infiles[i]);
-        if (fname.has_extname()) {
-            m_evtname[i] = fname.extname();
-        }
-        m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]);
-
-        // Write input file information into logger
-        log_value(NORMAL, "Input filename", m_infiles[i]);
-        log_value(NORMAL, "Event extension name", m_evtname[i]);
-        log_value(NORMAL, "GTI extension name", m_gtiname[i]);
-
-        // Fall through in case that the event file is empty
-        if (obs->events()->size() == 0) {
-            log_string(NORMAL, " Warning: No events in event file \""+
-                       m_infiles[i]+"\". Event phasing skipped.");
-            continue;
-        }
+        // Compute event phase
+        phase_events(obs);
  
-        // If we have an input file then check it
-        if (!m_infiles[i].empty()) {
-            std::string message = check_infile(m_infiles[i], m_evtname[i]);
-            if (!message.empty()) {
-                throw GException::invalid_value(G_RUN, message);
-            }
-        }
-
-        // Get temporary file name
-        #if G_USE_MKSTEMP
-        char tpl[]  = "ctphaseXXXXXX";
-        int  fileid = mkstemp(tpl);
-        std::string filename(tpl);
-        #else
-        std::string filename = std::tmpnam(NULL);
-        #endif
-
-        // Save observation in temporary file. We add here the events and
-        // GTI extension name so that the GCTAObservation::save method can
-        // use this information for writing the proper extension names into
-        // the temporary file
-        obs->save(filename+"["+m_evtname[i]+";"+ m_gtiname[i]+"]", true);
-
-        // Log saved FITS file.
-        if (logExplicit()) {
-            GFits tmpfile(filename);
-            log.header3("FITS file content of temporary file");
-            log << tmpfile << std::endl;
-            tmpfile.close();
-        }
-
-        // If we have a temporary file then check it
-        if (!filename.empty()) {
-            std::string message = check_infile(filename, m_evtname[i]);
-            if (!message.empty()) {
-                throw GException::invalid_value(G_RUN, message);
-            }
-        }
-
-        // Load observation from temporary file and compute phase
-        phase_events(obs, filename, m_evtname[i], m_gtiname[i]);
-
-        // Close temporary file
-        #if G_USE_MKSTEMP
-        close(fileid);
-        #endif
-
-        // Remove temporary file
-        std::remove(filename.c_str());
-
-    } // endfor: looped over all observations
+    } // endfor: looped over observations
  
     // If more than a single observation has been handled then make sure that
     // an XML file will be used for storage
@@ -380,12 +282,12 @@ void ctphase::save(void)
  
     // Case A: Save event file(s) and XML metadata information
     if (m_use_xml) {
-        save_xml();
+        save_events_xml();
     }
  
     // Case B: Save event file as FITS file
     else {
-        save_fits();
+        save_events_fits();
     }
  
     // Return
@@ -459,15 +361,10 @@ void ctphase::publish(const std::string&amp; name)
 void ctphase::init_members(void)
 {
     // Initialise parameters
-    m_outobs.clear();
-    m_prefix.clear();
-    m_phase.clear();
     m_chatter = static_cast<GChatter>(2);
  
     // Initialise protected members
-    m_infiles.clear();
-    m_evtname.clear();
-    m_gtiname.clear();
+    m_phase.clear();
  
     // Return
     return;
@@ -482,15 +379,10 @@ void ctphase::init_members(void)
 void ctphase::copy_members(const ctphase& app)
 {
     // Copy parameters
-    m_outobs   = app.m_outobs;
-    m_prefix   = app.m_prefix;
-    m_chatter  = app.m_chatter;
+    m_chatter = app.m_chatter;
  
     // Copy protected members
-    m_infiles       = app.m_infiles;
-    m_evtname       = app.m_evtname;
-    m_gtiname       = app.m_gtiname;
-    m_phase         = app.m_phase;
+    m_phase = app.m_phase;
  
     // Return
     return;
@@ -502,8 +394,6 @@ void ctphase::copy_members(const ctphase&amp; app)
  ***************************************************************************/
 void ctphase::free_members(void)
 {
-    m_phase.clear();
-    
     // Return
     return;
 }
@@ -526,381 +416,81 @@ void ctphase::get_parameters(void)
     // information and do not accept counts cubes.
     setup_observations(m_obs, false, true, false);
  
-    // Check whether a model XML file and source name were supplied
-    if (!((*this)["inmodel"].filename().file().empty()) &&
-        !((*this)["source"].string().empty())) {
-        
-        // Get the information from an xml file
-        read_phase_info_from_xml();
-    }
-    
-    // Check for sanity of frequency and phase information
-    else if ((*this)["p0"].is_valid() &&
-             (*this)["mjd"].is_valid() &&
-             (*this)["f0"].is_valid()&&
-             (*this)["f1"].is_valid()&&
-             (*this)["f2"].is_valid()) {
-        
-        // The following creates an XML object responsible for initializing
-        // 'm_phase'. This was done as the default values for the parameters in
-        // 'GModelTemporalPhaseCurve' are such that the ranges are predefined.
-        // Initializing 'm_phase' from an XML object circumvents this restriction.
-        
-        GXmlElement norm;
-        norm.name("parameter");
-        norm.attribute("name", "Normalization");
-        norm.attribute("value", "1");
-        norm.attribute("min", "0");
-        norm.attribute("max", "1000");
-        
-        GXmlElement mjd;
-        mjd.name("parameter");
-        mjd.attribute("name", "MJD");
-        mjd.attribute("value", (*this)["mjd"].value());
-        mjd.attribute("min", (*this)["mjd"].min());
-        mjd.attribute("max", (*this)["mjd"].max());
-        
-        GXmlElement p0;
-        p0.name("parameter");
-        p0.attribute("name", "Phase");
-        p0.attribute("value", (*this)["p0"].value());
-        p0.attribute("min", (*this)["p0"].min());
-        p0.attribute("max", (*this)["p0"].max());
-        
-        GXmlElement f0;
-        f0.name("parameter");
-        f0.attribute("name", "F0");
-        f0.attribute("value", (*this)["f0"].value());
-        f0.attribute("min", (*this)["f0"].min());
-        f0.attribute("max", (*this)["f0"].max());
-        
-        GXmlElement f1;
-        f1.name("parameter");
-        f1.attribute("name", "F1");
-        f1.attribute("value", (*this)["f1"].value());
-        f1.attribute("min", (*this)["f1"].min());
-        f1.attribute("max", (*this)["f1"].max());
-        
-        GXmlElement f2;
-        f2.name("parameter");
-        f2.attribute("name", "F2");
-        f2.attribute("value", (*this)["f2"].value());
-        f2.attribute("min", (*this)["f2"].min());
-        f2.attribute("max", (*this)["f2"].max());
-        
-        GXmlElement info;
-        info.append(norm);
-        info.append(mjd);
-        info.append(p0);
-        info.append(f0);
-        info.append(f1);
-        info.append(f2);
-        
-        info.attribute("file","");
-        
-        m_phase = GModelTemporalPhaseCurve();
-        m_phase.read(info);
-    }
-    else {
-        GException::invalid_value("ctphase::get_parameters()",
-                                  "Invalid parameter value for p0, mjd, f0, f1, or f2.");
-    }
-
-    // Get other User parameters
-    m_chatter  = static_cast<GChatter>((*this)["chatter"].integer());
-
-    // Optionally read ahead parameters so that they get correctly
-    // dumped into the log file
-    if (read_ahead()) {
-        m_outobs = (*this)["outobs"].filename();
-        m_prefix = (*this)["prefix"].string();
-    }
-
-    // Write parameters into logger
-    log_parameters(TERSE);
-
-    // Return
-    return;
-}
-
-/*************************************************************************
- * @brief Read phase information from XML file
- *
- * @exception GException::par_error
- *            Source does not have a temporal component with phase information.
- ***************************************************************************/
-void ctphase::read_phase_info_from_xml(void)
-{
-    // Load the models from the xml file
-    GModels in_models((*this)["inmodel"].filename());
-    
-    // Get the temporal component of this source
-    std::string src_name = (*this)["source"].string();
-    GModelSky* model = dynamic_cast<GModelSky*>(in_models[src_name]);
-    GModelTemporalPhaseCurve* temporal_model = dynamic_cast<GModelTemporalPhaseCurve*>(model->temporal());
-    
-    // Make sure the model has a temporal component with phase info
-    if (temporal_model != NULL) {
-        
-        // Copy the phase model
-        m_phase = GModelTemporalPhaseCurve(*temporal_model);
-        
-    } else {
-        // Throw an error as there is no temporal component for this source
-        throw GException::par_error("ctphase::get_parameters(void)",
-                                    "source",
-                                    "\"" + src_name + "\" source model " +
-                                    "does not have a temporal component with "+
-                                    "phase information!");
-    }
-}
-
-
-/***********************************************************************//**
- * @brief Define phase for events
- *
- * @param[in,out] obs CTA observation.
- * @param[in] filename File name.
- * @param[in] evtname Event extension name.
- * @param[in] gtiname GTI extension name.
- *
- * @exception GException::invalid_value
- *            No events extension found in FITS file.
- *
- * Append phase information to events from a FITS file by adding a "PHASE"
- * column to the FITS file.
- ***************************************************************************/
-void ctphase::phase_events(GCTAObservation*   obs,
-                             const std::string& filename,
-                             const std::string& evtname,
-                             const std::string& gtiname)
-{
-    // Get CTA event list pointer
-    GCTAEventList* list =
-        static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
-    
-    GCTAEventAtom* event = (*list)[0];
-    
-    // Print a warning into the log file if there are some precision concerns
-    double dt = std::abs(event->time() - m_phase.mjd());
-    if (((m_phase.f1() * dt) > 1.0e8) ||
-        (m_phase.f2() * dt * dt) > 1.0e8) {
-        log_string(NORMAL,"*** WARNING ***");
-        log_string(NORMAL,"   Values supplied for reference MJD and/or f1 and/or f2");
-        log_string(NORMAL,"   are large and may result in numerical precision issues.");
-        log_string(NORMAL,"   Consider using an ephemeris derived from a date closer");
-        log_string(NORMAL,"   to the data being analyzed.");
-        log_value(NORMAL, "     Observation ID", obs->obs_id());
-        log_value(NORMAL, "     Reference MJD", m_phase.mjd().mjd());
-        log_value(NORMAL, "     First event MJD", event->time().mjd());
-        log_string(NORMAL,"*** WARNING ***");
-    }
-    
-    for (int evnt=0; evnt<list->size(); evnt++) {
-        // Get the next event
-        event = (*list)[evnt];
-        
-        // Get the event time
-        GTime time = event->time();
-        
-        // Set the event phase
-        event->phase(m_phase.phase(time));
-    }
-    
-    // Make sure the file knows to save the phase information
-    if (!list->has_phase()) {
-        list->has_phase(true);
-    }
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Check input filename
- *
- * @param[in] filename File name.
- * @param[in] evtname Event extension name.
- *
- * This method checks if the input FITS file is correct.
- ***************************************************************************/
-std::string ctphase::check_infile(const std::string& filename,
-                                  const std::string& evtname) const
-{
-    // Initialise message string
-    std::string message = "";
-
-    // Open FITS file
-    GFits fits(filename);
-
-    // Check for existence of events extensions
-    if (!fits.contains(evtname)) {
-        message = "No \""+evtname+"\" extension found in input file \""+
-                  filename + "\".";
-    }
-
-    // ... otherwise check column names
-    else {
+    // If a model definition XML file is provided the
+    GFilename inmodel = (*this)["inmodel"].filename();
+    if (is_valid_filename(inmodel)) {
  
-        // Get pointer to FITS table
-        GFitsTable* table = fits.table(evtname);
+        // Load the models from the XML file
+        GModels models(inmodel);
  
-        // Initialise list of missing columns
-        std::vector<std::string> missing;
+        // Get the source component
+        std::string source = (*this)["source"].string();
  
-        // Check for existence of TIME column
-        if (!table->contains("TIME")) {
-            missing.push_back("TIME");
+        // Throw an exception if source is not found in the model container
+        if (!models.contains(source)) {
+            std::string msg = "Source model \""+source+"\" not found in model "
+                              "definition XML file \""+inmodel.url()+"\".";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
         }
  
-        // Check for existence of ENERGY column
-        if (!table->contains("ENERGY")) {
-            missing.push_back("ENERGY");
-        }
+        // Get pointer to source model
+        GModelSky* sky = dynamic_cast<GModelSky*>(models[source]);
  
-        // Check for existence of RA column
-        if (!table->contains("RA")) {
-            missing.push_back("RA");
+        // Throw an exception if model is not a sky model
+        if (sky == NULL) {
+            std::string msg = "Source model \""+source+"\" is not a sky model.";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
         }
  
-        // Check for existence of DEC column
-        if (!table->contains("DEC")) {
-            missing.push_back("DEC");
-        }
+        // Get pointer to temporal phase curve
+        GModelTemporalPhaseCurve* phasecurve =
+              dynamic_cast<GModelTemporalPhaseCurve*>(sky->temporal());
  
-        // Set error message for missing columns
-        if (!missing.empty()) {
-            message = "The following columns are missing in the "
-                      "\""+evtname+"\" extension of input file \""+
-                      filename + "\": ";
-            for (int i = 0; i < missing.size(); ++i) {
-                message += "\"" + missing[i] + "\"";
-                if (i < missing.size()-1) {
-                    message += ", ";
-                }
-            }
+        // Throw an exception if model has not a temporal phase curve
+        // component
+        if (phasecurve == NULL) {
+            std::string msg = "Source model \""+source+"\" does not have a "
+                              "temporal component with phase information.";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
         }
  
-    } // endelse: checked column names
-
-    // Return
-    return message;
-}
-
-
-/***********************************************************************//**
- * @brief Set output file name.
- *
- * @param[in] filename Input file name.
- *
- * Converts an input file name into an output filename by prepending the
- * prefix stored in the member m_prefix to the input file name. Any path as
- * well as extension will be stripped from the input file name. Also a
- * trailing ".gz" will be stripped as one cannot write into gzipped files.
- ***************************************************************************/
-std::string ctphase::set_outfile_name(const std::string& filename) const
-{
-    // Create filename
-    GFilename fname(filename);
-
-    // Split input filename without any extensions into path elements
-    std::vector<std::string> elements = gammalib::split(fname.url(), "/");
-
-    // The last path element is the filename
-    std::string outname = m_prefix + elements[elements.size()-1];
-
-    // Strip any ".gz"
-    outname = gammalib::strip_chars(outname, ".gz");
-    
-    // Return output filename
-    return outname;
-}
-
-
-/***********************************************************************//**
- * @brief Get Good Time Intervals extension name
- *
- * @param[in] filename Input file name.
- * @param[in] evtname Events extension name.
- *
- * Extracts the Good Time Intervals extension name from the event file. We
- * do this by loading the events and accessing the Good Time Intervals
- * extension name using the GCTAEventList::gtiname() method. If the file name
- * is empty, the method returns `GTI`.
- ***************************************************************************/
-std::string ctphase::get_gtiname(const std::string& filename,
-                                  const std::string& evtname) const
-{
-    // Initialise GTI name
-    std::string gtiname = gammalib::extname_gti;
-
-    // Continue only if the filename is not empty
-    if (!filename.empty()) {
+        // Set phase curve
+        m_phase = GModelTemporalPhaseCurve(*phasecurve);
  
-        // Load events
-        GCTAEventList events(filename+"["+evtname+"]");
+    } // endif: model definition XML file provided
  
-        // Get GTI name
-        gtiname = events.gtiname();
-
-    } // endif: filename was not empty
-
-    // Return GTI name
-    return (gtiname);
-}
-
-
-/***********************************************************************//**
- * @brief Save event list in FITS format.
- *
- * Save the event list as a FITS file. The file name of the FITS file is
- * specified by the "outobs" parameter.
- ***************************************************************************/
-void ctphase::save_fits(void)
-{
-    // Save only if there are observations
-    if (m_obs.size() > 0) {
-
-        // Get output filename
-        m_outobs = (*this)["outobs"].filename();
-
-        // Get CTA observation from observation container
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
-
-        // Save only if it's a CTA observation
-        if (obs != NULL) {
-    
-            // Save only if file name is non-empty
-            if (m_infiles[0].length() > 0) {
+    // Otherwise read phase curve information from the parameter files
+    else {
  
-                // Create file name object
-                GFilename fname(m_outobs);
+        // Read phase curve user parameters
+        GFilename phasecurve = (*this)["phasecurve"].filename();
+        double    mjd_value  = (*this)["mjd"].real();
+        double    phase      = (*this)["phase"].real();
+        double    f0         = (*this)["f0"].real();
+        double    f1         = (*this)["f1"].real();
+        double    f2         = (*this)["f2"].real();
  
-                // Extract filename and event extension name
-                std::string outfile = fname.url();
+        // Set reference time
+        GTime mjd;
+        mjd.mjd(mjd_value);
  
-                // Append event extension name. We handle here the possibility
-                // to write the events into a different extension.
-                if (fname.has_extname()) {
-                    outfile += "["+fname.extname()+"]";
-                }
-                else {
-                    outfile += "["+m_evtname[0]+"]";
-                }
+        // Set phase curve
+        m_phase = GModelTemporalPhaseCurve(phasecurve, mjd, phase, f0, f1, f2);
  
-                // Log filename
-                log_value(NORMAL, "Event list file", outfile);
+    } // endelse: read phase curve information from User parameters
  
-                // Save event list
-                save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0],
-                                outfile);
-
-            } // endif: filename was non empty
+    // Get other User parameters
+    m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
  
-        } // endif: observation was CTA observation
+    // Optionally read ahead parameters so that they get correctly
+    // dumped into the log file
+    if (read_ahead()) {
+        (*this)["outobs"].filename();
+        (*this)["prefix"].string();
+    }
  
-    } // endif: there were observations
+    // Write parameters into logger
+    log_parameters(TERSE);
  
     // Return
     return;
@@ -908,169 +498,68 @@ void ctphase::save_fits(void)
  
  
 /***********************************************************************//**
- * @brief Save event list(s) in XML format.
+ * @brief Compute event phase for an observation
  *
- * Save the event list(s) into FITS files and write the file path information
- * into a XML file. The filename of the XML file is specified by the outfile
- * parameter, the filename(s) of the event lists are built by prepending a
- * prefix to the input event list filenames. Any path present in the input
- * filename will be stripped, i.e. the event list(s) will be written in the
- * local working directory (unless a path is specified in the prefix).
+ * @param[in,out] obs CTA observation.
+ *
+ * @exception GException::invalid_value
+ *            No events extension found in FITS file.
+ *
+ * Append phase information to events from a FITS file by adding a "PHASE"
+ * column to the FITS file.
  ***************************************************************************/
-void ctphase::save_xml(void)
+void ctphase::phase_events(GCTAObservation* obs)
 {
-    // Get output filename and prefix
-    m_outobs = (*this)["outobs"].filename();
-    m_prefix = (*this)["prefix"].string();
-
-    // Issue warning if output filename has no .xml suffix
-    log_string(TERSE, warn_xml_suffix(m_outobs));
-
-    // Loop over all observation in the container
-    for (int i = 0; i < m_obs.size(); ++i) {
-
-        // Get CTA observation
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Skip observations that are no CTA observations
-        if (obs == NULL) {
-            continue;
-        }
+    // Write header into logger
+    log_header3(NORMAL, "Compute event phases");
+
+    // Make sure that the observation holds a CTA event list. If this
+    // is not the case then throw an exception.
+    GCTAEventList* events = dynamic_cast<GCTAEventList*>(const_cast<GEvents*>
+                            (obs->events()));
+    if (events == NULL) {
+        std::string msg = "CTA Observation does not contain an event list. An "
+                          "event list is needed to compute the event phases.";
+        throw GException::invalid_value(G_PHASE_EVENTS, msg);
+    }
  
-        // Skip observations that have empty names
-        if (m_infiles[i].length() == 0) {
-            continue;
+    // Continue only if there are events
+    if (events->size() > 0) {
+
+        // Print a warning into the log file if there are some precision concerns
+        GCTAEventAtom* event = (*events)[0];
+        double         dt    = std::abs(event->time() - m_phase.mjd());
+        if (((m_phase.f1() * dt) > 1.0e8) ||
+             (m_phase.f2() * dt * dt) > 1.0e8) {
+            log_string(NORMAL,"*** WARNING ***");
+            log_string(NORMAL,"   Values supplied for reference MJD and/or f1 and/or f2");
+            log_string(NORMAL,"   are large and may result in numerical precision issues.");
+            log_string(NORMAL,"   Consider using an ephemeris derived from a date closer");
+            log_string(NORMAL,"   to the data being analyzed.");
+            log_value(NORMAL, "     Observation ID", obs->obs_id());
+            log_value(NORMAL, "     Reference MJD", m_phase.mjd().mjd());
+            log_value(NORMAL, "     First event MJD", event->time().mjd());
+            log_string(NORMAL,"*** WARNING ***");
         }
  
-        // Set event output file name
-        std::string outfile = set_outfile_name(m_infiles[i]);
-
-        // Append event extension name
-        outfile += "["+m_evtname[i]+"]";
-
-        // Log filename
-        log_value(NORMAL, "Event list file", outfile);
-
-        // Store output file name in observation
-        obs->eventfile(outfile);
-
-        // Save event list
-        save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i],
-                        outfile);
-
-    } // endfor: looped over observations
-
-    // Write observation definition XML file name into logger
-    log_value(NORMAL, "Obs. definition file", m_outobs);
-
-    // Save observations in XML file
-    m_obs.save(m_outobs);
-
-    // Return
-    return;
-}
-
+        // Compute the phase for all events
+        for (int i = 0; i < events->size(); ++i) {
  
-/***********************************************************************//**
- * @brief Save event list into FITS file
- *
- * @param[in] obs Pointer to CTA observation.
- * @param[in] infile Input file name.
- * @param[in] evtname Event extension name.
- * @param[in] gtiname GTI extension name.
- * @param[in] outfile Output file name.
- *
- * Saves an event list and the corresponding Good Time Intervals into a FITS
- * file and copy all others extensions from the input file to the output
- * file.
- *
- * If an extension name is specified in the @p outfile argument, the events
- * and eventually also the Good Time Intervals will be extracted from the
- * argument and used for writing the events. The format is
- *
- *      <filename>[<event extension name;GTI extension name>]
- *
- * where <filename> needs to be replaced by the name of the FITS file,
- * and <event extension name;GTI extension name> by the name of the events
- * and Good Time Intervals extensions. For example 
- *
- *      myfits.fits[EVENTS1;GTI1]
- *
- * will write the updated events into the "EVENTS1" extension and the
- * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
- * file. If the Good Time Intervals extension name is skipped, e.g.
- *
- *      myfits.fits[EVENTS1]
- *
- * the original extension name for the Good Time Intervals will be kept.
- * Analogously, only the Good Time Intervals extension name can be changed
- * by specifying
- *
- *      myfits.fits[;GTI1]
- *
- * In none of the cases will the original events and Good Time Intervals be
- * copied over to the output file.
- ***************************************************************************/
-void ctphase::save_event_list(const GCTAObservation* obs,
-                               const std::string&     infile,
-                               const std::string&     evtname,
-                               const std::string&     gtiname,
-                               const std::string&     outfile) const
-{
-    // Save only if we have an event list
-    if (obs->eventtype() == "EventList") {
-
-        // Set output FITS file event extension names
-        GFilename   outname(outfile);
-        std::string outevt = evtname;
-        std::string outgti = gtiname;
-        if (outname.has_extname()) {
-            std::vector<std::string> extnames =
-                       gammalib::split(outname.extname(), ";");
-            if (extnames.size() > 0) {
-                std::string extname = gammalib::strip_whitespace(extnames[0]);
-                if (!extname.empty()) {
-                    outevt = extname;
-                }
-            }
-            if (extnames.size() > 1) {
-                std::string extname = gammalib::strip_whitespace(extnames[1]);
-                if (!extname.empty()) {
-                    outgti = extname;
-                }
-            }
-        }
+            // Get pointer to the event
+            GCTAEventAtom* event = (*events)[i];
  
-        // Create output FITS file
-        GFits outfits;
-
-        // Write observation into FITS file
-        obs->write(outfits, outevt, outgti);
-
-        // Copy all extensions other than evtname and gtiname extensions
-        // from the input to the output event list. The evtname and
-        // gtiname extensions are written by the save method, all others
-        // that may eventually be present have to be copied over
-        // explicitly.
-        GFits infits(infile);
-        for (int extno = 1; extno < infits.size(); ++extno) {
-            GFitsHDU* hdu = infits.at(extno);
-            if (hdu->extname() != evtname &&
-                hdu->extname() != gtiname &&
-                hdu->extname() != outevt  &&
-                hdu->extname() != outgti) {
-                outfits.append(*hdu);
-            }
-        }
-
-        // Close input file
-        infits.close();
+            // Get the event time
+            GTime time = event->time();
+        
+            // Set the event phase
+            event->phase(m_phase.phase(time));
  
-        // Save file to disk and close it (we need both operations)
-        outfits.saveto(outname.url(), clobber());
-        outfits.close();
+        } // endfor: looped over events
  
-    } // endif: observation was unbinned
+    } // endif: there were events
+    
+    // Signal that phase information is available
+    events->has_phase(true);
  
     // Return
     return;
 /***************************************************************************
  *          ctphase - Append phase information to CTA events file          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Leonardo Di Venere                               *
+ *  copyright (C) 2017 by Joshua Cardenzana                                *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file ctphase.hpp
  * @brief Append phase information to CTA events file
- * @author Leonardo Di Venere
+ * @author Joshua Cardenzana
  */
  
 #ifndef CTPHASE_HPP
@@ -35,8 +35,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTPHASE_NAME    "ctphase"
-#define CTPHASE_VERSION "0.0.1"
+#define CTPHASE_NAME "ctphase"
  
  
 /***********************************************************************//**
@@ -65,38 +64,17 @@ public:
  
 protected:
     // Protected methods
-    void        init_members(void);
-    void        copy_members(const ctphase& app);
-    void        free_members(void);
-    void        get_parameters(void);
-    void        read_phase_info_from_xml(void);
-    void        phase_events(GCTAObservation* obs,
-                             const std::string& filename,
-                             const std::string& evtname,
-                             const std::string& gtiname);
-    std::string check_infile(const std::string& filename,
-                             const std::string& evtname) const;
-    std::string set_outfile_name(const std::string& filename) const;
-    std::string get_gtiname(const std::string& filename,
-                            const std::string& evtname) const;
-    void        save_fits(void);
-    void        save_xml(void);
-    void        save_event_list(const GCTAObservation* obs,
-                                const std::string&     infile,
-                                const std::string&     evtname,
-                                const std::string&     gtiname,
-                                const std::string&     outfile) const;
+    void init_members(void);
+    void copy_members(const ctphase& app);
+    void free_members(void);
+    void get_parameters(void);
+    void phase_events(GCTAObservation* obs);
  
     // User parameters
-    std::string m_outobs;               //!< Output event list or XML file
-    std::string m_prefix;               //!< Prefix for multiple event lists
-    GModelTemporalPhaseCurve m_phase;   //!< Phase information object
-    GChatter    m_chatter;              //!< Chattiness
+    GChatter m_chatter;                 //!< Chattiness
  
     // Protected members
-    std::vector<std::string> m_infiles;       //!< Input event filenames
-    std::vector<std::string> m_evtname;       //!< Event extension names
-    std::vector<std::string> m_gtiname;       //!< GTI extension names
+    GModelTemporalPhaseCurve m_phase;   //!< Phase curve model
 };
  
  
@@ -5,11 +5,10 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctphase.par                                                   #
-#  Version:   0.0.1                                                         #
 #                                                                           #
-#  Author:    Leonardo Di Venere                                            #
-#             leonardo.divenere@ba.infn.it                                  #
-#             INFN                                                          #
+#  Author:    Joshua Cardenzana                                             #
+#             jcardenzana@irap.omp.eu                                       #
+#             IRAP                                                          #
 #                                                                           #
 #  Purpose:   Parameter file for the ctphase application                    #
 #                                                                           #
@@ -17,20 +16,21 @@
 #
 # General parameters
 #===================
-inobs,    f, a,events.fits,,, "Input event list or observation definition XML file"
-outobs,   f, a,phased_events.fits,,, "Output event list or observation definition XML file"
-prefix,   s, h, "phased_",,, "Prefix for event lists in observation definition XML file"
-inmodel,  f, h, ,,, "Model file defining the phase information"
-source,   s, h, ,,, "Source in model file to use for phase information"
+inobs,   f, a, "events.fits",,, "Input event list or observation definition XML file"
+outobs,  f, a, "phased_events.fits",,, "Output event list or observation definition XML file"
+prefix,  s, h, "phased_",,, "Prefix for event lists in observation definition XML file"
+inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Model file defining the phase information"
+source,  s, a, "Crab",,, "Source in model file to use for phase information"
  
 #
 # Phase information parameters
 #=====================
-mjd,    r, a, 51544.5,,, "Reference MJD for phase information"
-p0,     r, a, 0.0,0.0,1.0, "Phase at reference MJD"
-f0,     r, a, 1,,, "Frequency at reference MJD [Hz]"
-f1,     r, a, 0.0,,, "First freq. derivative at reference MJD [s^-1]"
-f2,     r, a, 0.0,,, "Second freq. derivative at reference MJD [s^-2]"
+phasecurve, f, a, phase.fits,,, "Phase curve FITS file"
+mjd,        r, a, 51544.5,,, "Reference MJD for phase information"
+phase,      r, a, 0.0,0.0,1.0, "Phase at reference MJD"
+f0,         r, a, 1.0,,, "Frequency at reference MJD (Hz)"
+f1,         r, a, 0.1,,, "First freq. derivative at reference MJD (s^-1)"
+f2,         r, a, 0.01,,, "Second freq. derivative at reference MJD (s^-2)"
  
 #
 # Standard parameters
 /***************************************************************************
  *          ctphase - Append phase information to CTA events file          *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2017 by Leonardo Di Venere                               *
+ *  copyright (C) 2017 by Joshua Cardenzana                                *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file ctphase/main.cpp
  * @brief Append phase information to CTA events file
- * @author Leonardo Di Venere
+ * @author Joshua Cardenzana
  */
  
 /* __ Includes ___________________________________________________________ */
+#############################################################################
+# Makefile.am for ctprob tool                                               #
+# ------------------------------------------------------------------------- #
+# Copyright (C) 2017 Leonardo Di Venere                                     #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  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     #
+#  the Free Software Foundation, either version 3 of the License, or        #
+#  (at your option) any later version.                                      #
+#                                                                           #
+#  This program is distributed in the hope that it will be useful,          #
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of           #
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            #
+#  GNU General Public License for more details.                             #
+#                                                                           #
+#  You should have received a copy of the GNU General Public License        #
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.    #
+#                                                                           #
+#############################################################################
+
+# Set include directories
+INCLUDES = -I$(top_srcdir)/src/support
+
+# Program name
+bin_PROGRAMS = ctprob
+
+# Sources
+sources = ctprob.cpp ctprob.hpp
+
+# Libtool library
+noinst_LTLIBRARIES = libctprob.la
+libctprob_la_SOURCES = $(sources)
+libctprob_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@
+
+# Program definition
+ctprob_SOURCES = main.cpp
+ctprob_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
+ctprob_CPPFLAGS = @CPPFLAGS@
+ctprob_LDADD = libctprob.la $(top_builddir)/src/support/libsupport.la
+
+# Parameter files
+pfilesdir = $(prefix)/syspfiles
+pfiles_DATA = $(srcdir)/ctprob.par
+
+# Extra files to be distributed in the tar file
+EXTRA_DIST = $(pfiles_DATA)
+/***************************************************************************
+ *          ctprob - Computes event probability for a given model          *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2017 by Leonardo Di Venere                               *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctprob.cpp
+ * @brief Computes event probability for a given model
+ * @author Leonardo Di Venere
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdlib>
+#include "ctprob.hpp"
+#include "GTools.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_GET_OBS                                         "ctprob::get_obs()"
+#define G_EVALUATE_PROB      "ctprob::evaluate_probability(GCTAObservation*)"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+ctprob::ctprob(void) : ctobservation(CTPROB_NAME, VERSION)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Observations constructor
+ *
+ * param[in] obs Observation container.
+ *
+ * Creates an instance of the class that is initialised using the information
+ * provided in an observation container.
+ ***************************************************************************/
+ctprob::ctprob(const GObservations& obs) :
+        ctobservation(CTPROB_NAME, VERSION, obs)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+
+
+/***********************************************************************//**
+ * @brief Command line constructor
+ *
+ * @param[in] argc Number of arguments in command line.
+ * @param[in] argv Array of command line arguments.
+ ***************************************************************************/
+ctprob::ctprob(int argc, char *argv[]) : 
+        ctobservation(CTPROB_NAME, VERSION, argc, argv)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+ctprob::ctprob(const ctprob& app) : ctobservation(app)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(app);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+ctprob::~ctprob(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app Application.
+ * @return Application.
+ ***************************************************************************/
+ctprob& ctprob::operator=(const ctprob& app)
+{
+    // Execute only if object is not identical
+    if (this != &app) {
+
+        // Copy base class members
+        this->ctobservation::operator=(app);
+
+        // Free members
+        free_members();
+
+        // Initialise members
+        init_members();
+
+        // Copy members
+        copy_members(app);
+
+    } // endif: object was not identical
+
+    // Return this object
+    return *this;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                            Public methods                               =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Clear ctprob tool
+ *
+ * Clears ctprob tool.
+ ***************************************************************************/
+void ctprob::clear(void)
+{
+    // Free members
+    free_members();
+    this->ctobservation::free_members();
+    this->ctool::free_members();
+
+    // Clear base class (needed to conserve tool name and version)
+    this->GApplication::clear();
+
+    // Initialise members
+    this->ctool::init_members();
+    this->ctobservation::init_members();
+    init_members();
+
+    // Write header into logger
+    log_header();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Runs the ctprob tool
+ *
+ * This method reads in the application parameters and loops over all
+ * unbinned observations to compute the probability for each event that it
+ * comes from one of the model components.
+ ***************************************************************************/
+void ctprob::run(void)
+{
+    // Switch screen logging on in debug mode
+    if (logDebug()) {
+        log.cout(true);
+    }
+
+    // Get parameters
+    get_parameters();
+
+    // Set energy dispersion flags of all CTA observations and save old
+    // values in save_edisp vector
+    std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
+
+    // Write input observation container into logger
+    log_observations(NORMAL, m_obs, "Input observation");
+
+    // Write header into logger
+    log_header1(TERSE, "Compute event probabilities");
+
+    // Initialise counters
+    int n_observations = 0;
+
+    // Loop over all unbinned CTA observations in the container
+    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+         obs = next_unbinned_observation()) {
+
+        // Increment counter
+        n_observations++;
+
+        // Evaluate event probabilities
+        evaluate_probability(obs);
+
+    } // endfor: looped over observations
+
+    // If more than a single observation has been handled then make sure that
+    // an XML file will be used for storage
+    if (n_observations > 1) {
+        m_use_xml = true;
+    }
+
+    // Write observation(s) into logger
+    log_observations(NORMAL, m_obs, "Output observation");
+
+    // Restore energy dispersion flags of all CTA observations
+    restore_edisp(m_obs, save_edisp);
+
+    // Optionally publish event list(s)
+    if ((*this)["publish"].boolean()) {
+        publish();
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save the selected event list(s)
+ *
+ * This method saves the selected event list(s) into FITS file(s). There are
+ * two modes, depending on the m_use_xml flag.
+ *
+ * If m_use_xml is true, all selected event list(s) will be saved into FITS
+ * files, where the output filenames are constructued from the input
+ * filenames by prepending the m_prefix string to name. Any path information
+ * will be stripped form the input name, hence event files will be written
+ * into the local working directory (unless some path information is present
+ * in the prefix). In addition, an XML file will be created that gathers
+ * the filename information for the selected event list(s). If an XML file
+ * was present on input, all metadata information will be copied from this
+ * input file.
+ *
+ * If m_use_xml is false, the selected event list will be saved into a FITS
+ * file.
+ ***************************************************************************/
+void ctprob::save(void)
+{
+    // Write header into logger
+    log_header1(TERSE, gammalib::number("Save event list", m_obs.size()));
+
+    // Case A: Save event file(s) and XML metadata information
+    if (m_use_xml) {
+        save_xml();
+    }
+
+    // Case B: Save event file as FITS file
+    else {
+        save_fits();
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Publish event lists
+ *
+ * @param[in] name Event list name.
+ ***************************************************************************/
+void ctprob::publish(const std::string& name)
+{
+    // Write header into logger
+    log_header1(TERSE, gammalib::number("Publish event list", m_obs.size()));
+
+    // Loop over all observation in the container
+    for (int i = 0; i < m_obs.size(); ++i) {
+
+        // Get CTA observation
+        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
+
+        // Handle only CTA observations
+        if (obs != NULL) {
+
+            // Continue only if there is an event list
+            if (obs->events()->size() != 0) {
+
+                // Set default name if user name is empty
+                std::string user_name(name);
+                if (user_name.empty()) {
+                    user_name = CTPROB_NAME;
+                }
+
+                // If there are several event lists then add an index
+                if (m_use_xml) {
+                    user_name += gammalib::str(i);
+                }
+
+                // Write event list name into logger
+                log_value(NORMAL, "Event list name", user_name);
+
+                // Write events into in-memory FITS file
+                GFits fits;
+                obs->write(fits);
+
+                // Publish
+                fits.publish(gammalib::extname_cta_events, user_name);
+
+            } // endif: there were events
+
+        } // endif: observation was a CTA observation
+
+    } // endfor: looped over observations
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void ctprob::init_members(void)
+{
+    // Initialise parameters
+    m_outobs.clear();
+    m_prefix.clear();
+    m_chatter = static_cast<GChatter>(2);
+
+    // Initialise protected members
+    m_infiles.clear();
+    m_evtname.clear();
+    m_gtiname.clear();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctprob::copy_members(const ctprob& app)
+{
+    // Copy parameters
+    m_outobs   = app.m_outobs;
+    m_prefix   = app.m_prefix;
+    m_chatter  = app.m_chatter;
+
+    // Copy protected members
+    m_infiles       = app.m_infiles;
+    m_evtname       = app.m_evtname;
+    m_gtiname       = app.m_gtiname;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctprob::free_members(void)
+{
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * Get all user parameters from parameter file or (if required) by querying
+ * the user. Times are assumed to be in the native CTA MJD format.
+ *
+ * This method also loads observations if no observations are yet allocated.
+ * Observations are either loaded from a single CTA even list, or from a
+ * XML file using the metadata information that is stored in that file.
+ ***************************************************************************/
+void ctprob::get_parameters(void)
+{
+    // Setup observations from "inobs" parameter. Do not request response
+    // information and do not accept counts cubes.
+    if (m_obs.size() == 0) {
+        get_obs();
+    }
+    // ... otherwise add response information and energy boundaries in case
+    // that they are missing
+    else {
+        setup_observations(m_obs);
+    }
+
+    // Read model definition file if required
+    if (m_obs.models().size() == 0) {
+
+        // Get model filename
+        std::string inmodel = (*this)["inmodel"].filename();
+
+        // Load models from file
+        m_obs.models(inmodel);
+
+    } // endif: there were no models
+
+    // Get energy dispersion flag parameters
+    m_apply_edisp = (*this)["edisp"].boolean();
+
+
+    // Get remaining parameters
+    m_publish = (*this)["publish"].boolean();
+    m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
+
+    // Optionally read ahead parameters so that they get correctly
+    // dumped into the log file
+    if (read_ahead()) {
+        m_outobs = (*this)["outobs"].filename();
+        m_prefix = (*this)["prefix"].string();
+    }
+
+    // Write parameters into logger
+    log_parameters(TERSE);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get observation container
+ *
+ * @exception GException::invalid_value
+ *            Invalid FITS file.
+ *
+ * Get an observation container according to the user parameters. The method
+ * supports loading of a individual FITS file or an observation definition
+ * file in XML format.
+ ***************************************************************************/
+void ctprob::get_obs(void)
+{
+    // Get the filename from the input parameters
+    std::string filename = (*this)["inobs"].filename();
+
+    // If no observation definition file throw an exception 
+    if (!is_valid_filename(filename)) {
+      throw GException::invalid_value(G_GET_OBS, "Input event list is not valid.");
+    } // endif: filename was not valid
+
+    // ... otherwise we have a file name
+    else {
+
+        // If file is a FITS file then create an empty CTA observation
+        // and load file into observation
+        if (GFilename(filename).is_fits()) {
+
+            // Allocate empty CTA observation
+            GCTAObservation cta;
+
+            // Load data
+            cta.load(filename);
+
+            // Set response
+            set_obs_response(&cta);
+
+            // Append observation to container
+            m_obs.append(cta);
+
+            // Signal that no XML file should be used for storage
+            m_use_xml = false;
+
+        }
+
+        // ... otherwise load file into observation container
+        else {
+
+            // Load observations from XML file
+            m_obs.load(filename);
+
+            // For all observations that have no response, set the response
+            // from the task parameters
+            set_response(m_obs);
+
+            // Set observation boundary parameters (emin, emax, rad)
+            set_obs_bounds(m_obs);
+
+            // Signal that XML file should be used for storage
+            m_use_xml = true;
+
+        } // endelse: file was an XML file
+
+    }
+
+    // Return
+    return;
+
+}
+
+
+
+/***********************************************************************//**
+ * @brief Evaluate probability for events
+ *
+ * @param[in,out] obs CTA observation.
+ *
+ * @exception GException::invalid_value
+ *            Observation does not contain an event list.
+ *
+ * Computes for all events the probability that an event arises from a
+ * specific model component. The computation is done by by evaluating the
+ * differential event probability for a given model and normalizing these
+ * probabilities so that the sum of the probabilities for all model
+ * components is unity for each event.
+ *
+ * The method appends for each model component a single precision column to
+ * the event list that contains the event probabilities. The name of the
+ * column is build from the model name prefixed with `PROB_`.
+ ***************************************************************************/
+void ctprob::evaluate_probability(GCTAObservation* obs)
+{
+    // Write header into logger
+    log_header3(NORMAL, "Evaluate events probabilities");
+
+    // Make sure that the observation holds a CTA event list. If this
+    // is not the case then throw an exception.
+    GCTAEventList* events = dynamic_cast<GCTAEventList*>(const_cast<GEvents*>
+                            (obs->events()));
+    if (events == NULL) {
+        std::string msg = "CTA Observation does not contain an event list. An "
+                          "event list is needed to compute the event "
+                          "probabilities.";
+        throw GException::invalid_value(G_EVALUATE_PROB, msg);
+    }
+
+    // Continue only if there are model components
+    int nmodels = m_obs.models().size();
+    if (nmodels > 0) {
+
+        // Allocate single precision FITS table columns for all model components
+        std::vector<GFitsTableFloatCol*> columns;
+        for (int k = 0; k < nmodels; ++k) {
+            const std::string   colname = "PROB_" + m_obs.models()[k] ->name();
+            GFitsTableFloatCol* col     =
+                                new GFitsTableFloatCol(colname, events->size());
+            columns.push_back(col);
+        }
+
+        // Loop over all events in the event list
+        for (int i = 0; i < events->size(); ++i) {
+
+            // Initialise total probability
+            double total = 0.0;
+
+            // Get point to event
+            const GEvent* event = (*events)[i];
+
+            // Initialise probability vector
+            std::vector<double> probabilities(nmodels, 0.0);
+
+            // Loop over models
+            for (int k = 0; k < nmodels; ++k) {
+                double value      =  m_obs.models()[k]->eval(*event, *obs);
+                probabilities[k]  = value;
+                total            += value;
+            }
+
+            // Normalize probabilities to unity
+            if (total > 0.0) {
+                for (int k = 0; k < nmodels; ++k) {
+                    probabilities[k] /= total;
+                }
+            }
+
+            // Store probabilities in FITS columns
+            for (int k = 0; k < nmodels; ++k) {
+                (*(columns[k]))(i) = probabilities[k];
+            }
+	
+        } // endfor: looped over all events
+
+        // Append columns to event list
+        for (int k = 0; k < nmodels; ++k) {
+            events->append_column(*columns[k]);
+        }
+
+    } // endif: there were model components
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set output file name.
+ *
+ * @param[in] filename Input file name.
+ *
+ * Converts an input file name into an output filename by prepending the
+ * prefix stored in the member m_prefix to the input file name. Any path as
+ * well as extension will be stripped from the input file name. Also a
+ * trailing ".gz" will be stripped as one cannot write into gzipped files.
+ ***************************************************************************/
+std::string ctprob::set_outfile_name(const std::string& filename) const
+{
+    // Create filename
+    GFilename fname(filename);
+
+    // Split input filename without any extensions into path elements
+    std::vector<std::string> elements = gammalib::split(fname.url(), "/");
+
+    // The last path element is the filename
+    std::string outname = m_prefix + elements[elements.size()-1];
+
+    // Strip any ".gz"
+    outname = gammalib::strip_chars(outname, ".gz");
+    
+    // Return output filename
+    return outname;
+}
+
+
+/***********************************************************************//**
+ * @brief Save event list in FITS format.
+ *
+ * Save the event list as a FITS file. The file name of the FITS file is
+ * specified by the `outobs` parameter.
+ ***************************************************************************/
+void ctprob::save_fits(void)
+{
+    // Save only if there are observations
+    if (m_obs.size() > 0) {
+
+        // Get output filename
+        m_outobs = (*this)["outobs"].filename();
+
+        // Loop over all unbinned CTA observations in the container
+        for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+             obs = next_unbinned_observation()) {
+
+            // Get input file name and default extension names
+            std::string infile  = obs->eventfile();
+            std::string evtname = gammalib::extname_cta_events;
+            std::string gtiname = gammalib::extname_gti;
+
+            // Create file name object
+            GFilename fname(m_outobs);
+
+            // Extract filename and event extension name
+            std::string outfile = fname.url();
+
+            // Append event extension name. We handle here the possibility
+            // to write the events into a different extension.
+            if (fname.has_extname()) {
+                outfile += "["+fname.extname()+"]";
+            }
+            else {
+                outfile += "["+gammalib::extname_cta_events+"]";
+            }
+
+            // Log filename
+            log_value(NORMAL, "Event list file", outfile);
+            
+            // Save event list
+            save_event_list(obs, infile, evtname, gtiname, outfile);
+
+            // Exit the loop (if this method is called we should only
+            // have a single unbinned observation in the observation
+            // container
+            break;
+
+        } // endfor: looped over all unbinned observations
+
+    } // endif: there were observations
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save event list(s) in XML format.
+ *
+ * Save the event list(s) into FITS files and write the file path information
+ * into a XML file. The filename of the XML file is specified by the `outobs`
+ * parameter, the filename(s) of the event lists are built by prepending a
+ * prefix to the input event list filenames. Any path present in the input
+ * filename will be stripped, i.e. the event list(s) will be written in the
+ * local working directory (unless a path is specified in the prefix).
+ ***************************************************************************/
+void ctprob::save_xml(void)
+{
+    // Get output filename and prefix
+    m_outobs = (*this)["outobs"].filename();
+    m_prefix = (*this)["prefix"].string();
+
+    // Issue warning if output filename has no .xml suffix
+    log_string(TERSE, warn_xml_suffix(m_outobs));
+
+    // Loop over all unbinned CTA observations in the container
+    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+         obs = next_unbinned_observation()) {
+
+        // Get input file name and default extension names
+        std::string infile  = obs->eventfile();
+        std::string evtname = gammalib::extname_cta_events;
+        std::string gtiname = gammalib::extname_gti;
+
+        // Extract event and GTI extension names from input FITS file
+        GFilename fname(infile);
+        if (fname.has_extname()) {
+            evtname = fname.extname();
+        }
+        gtiname = get_gtiname(fname.url(), evtname);
+
+        // Set event output file name
+        std::string outfile = set_outfile_name(infile);
+
+        // Append event extension name
+        outfile += "["+evtname+"]";
+
+        // Log filename
+        log_value(NORMAL, "Event list file", outfile);
+
+        // Store output file name in observation
+        obs->eventfile(outfile);
+
+        // Save event list
+        save_event_list(obs, infile, evtname, gtiname, outfile);
+
+    } // endfor: looped over observations
+
+    // Write observation definition XML file name into logger
+    log_value(NORMAL, "Obs. definition file", m_outobs);
+
+    // Save observations in XML file
+    m_obs.save(m_outobs);
+
+    // Return
+    return;
+}
+/***************************************************************************
+ *          ctprob - Computes event probability for a given model          *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2017 by Leonardo Di Venere                               *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctprob.hpp
+ * @brief Computes event probability for a given model
+ * @author Leonardo Di Venere
+ */
+
+#ifndef CTPROB_HPP
+#define CTPROB_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include <vector>
+#include <string>
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+#include "ctobservation.hpp"
+
+/* __Definitions _________________________________________________________ */
+#define CTPROB_NAME "ctprob"
+
+
+/***********************************************************************//**
+ * @class ctprob
+ *
+ * @brief Data selection tool
+ ***************************************************************************/
+class ctprob : public ctobservation {
+
+public:
+    // Constructors and destructors
+    ctprob(void);
+    explicit ctprob(const GObservations& obs);
+    ctprob(int argc, char *argv[]);
+    ctprob(const ctprob& app);
+    virtual ~ctprob(void);
+
+    // Operators
+    ctprob& operator=(const ctprob& app);
+
+    // Methods
+    void clear(void);
+    void run(void);
+    void save(void);
+    void publish(const std::string& name = "");
+
+protected:
+    // Protected methods
+    void        init_members(void);
+    void        copy_members(const ctprob& app);
+    void        free_members(void);
+    void        get_parameters(void);
+    void        get_obs(void);
+    void        evaluate_probability(GCTAObservation* obs);
+    std::string set_outfile_name(const std::string& filename) const;
+    void        save_fits(void);
+    void        save_xml(void);
+
+    // User parameters
+    std::string m_outobs;      //!< Output event list or XML file
+    std::string m_prefix;      //!< Prefix for multiple event lists
+    bool        m_apply_edisp; //!< Apply energy dispersion?
+    bool        m_publish;     //!< Publish event list?
+    GChatter    m_chatter;     //!< Chattiness
+
+    // Protected members
+    std::vector<std::string> m_infiles;  //!< Input event filenames
+    std::vector<std::string> m_evtname;  //!< Event extension names
+    std::vector<std::string> m_gtiname;  //!< GTI extension names
+};
+
+#endif /* CTPROB_HPP */
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      ctprob.par                                                    #
+#                                                                           #
+#  Author:    Leonardo Di Venere                                            #
+#             leonardo.divenere@ba.infn.it                                  #
+#             INFN                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the ctprob application                     #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+inobs,   f, a, "events.fits",,, "Input event list or observation definition XML file"
+outobs,  f, a, "prob_events.fits",,, "Output event list or observation definition XML file"
+prefix,  s, h, "prob_",,, "Prefix for event lists in observation definition XML file"
+inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
+caldb,   s, a, "prod2",,, "Calibration database"
+irf,     s, a, "South_0.5h",,, "Instrument response function"
+edisp,   b, h, no,,, "Apply energy dispersion?"
+
+#
+# Standard parameters
+#====================
+publish, b, h, no,,, "Publish event list on VO Hub?"
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug,   b, h, no,,, "Debugging mode activated"
+mode,    s, h, "ql",,, "Mode of automatic parameters"
+logfile, f, h, "ctprob.log",,, "Log filename"
+/***************************************************************************
+ *          ctprob - Computes event probability for a given model          *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2017 by Leonardo Di Venere                               *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file main.cpp
+ * @brief Computes event probability for a given model
+ * @author Leonardo Di Venere
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "support.hpp"
+#include "ctprob.hpp"
+
+
+/***********************************************************************//**
+ * @brief Main entry point of application
+ *
+ * @param[in] argc Number of command line arguments.
+ * @param[in] argv Command line arguments.
+ *
+ * This is the main entry point of the ctprob application. It allocates a
+ * ctprob object and executes the application. Any exceptions that occur
+ * will be catched and corresponding error messages written in the
+ * application logger and into the standard output.
+ ***************************************************************************/
+int main (int argc, char *argv[])
+{
+    // Initialise return code
+    int rc = 1;
+
+    // Initialise pointer on application
+    ctprob* application = NULL;
+
+    // Execute application
+    try {
+
+        // Create instance of application
+        application = new ctprob(argc, argv);
+
+        // Execute ctool
+        rc = execute_ctool(application);
+
+    }
+
+    catch (std::exception &e) {
+
+        // Report exception
+        report_ctool_failure("ctprob", e.what());
+
+    }
+
+    // Delete application
+    delete application;
+
+    // Return
+    return rc;
+}
@@ -49,7 +49,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
+ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -68,7 +68,7 @@ ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
  * observations container.
  ***************************************************************************/
 ctpsfcube::ctpsfcube(const GObservations& obs) :
-           ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, obs)
+           ctobservation(CTPSFCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -86,7 +86,7 @@ ctpsfcube::ctpsfcube(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctpsfcube::ctpsfcube(int argc, char *argv[]) :
-           ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, argc, argv)
+           ctobservation(CTPSFCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTPSFCUBE_NAME    "ctpsfcube"
-#define CTPSFCUBE_VERSION "1.1.0"
+#define CTPSFCUBE_NAME "ctpsfcube"
  
  
 /***********************************************************************//**
@@ -53,7 +53,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, CTSELECT_VERSION)
+ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -72,7 +72,7 @@ ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, CTSELECT_VERSION)
  * provided in an observation container.
  ***************************************************************************/
 ctselect::ctselect(const GObservations& obs) :
-          ctobservation(CTSELECT_NAME, CTSELECT_VERSION, obs)
+          ctobservation(CTSELECT_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -91,7 +91,7 @@ ctselect::ctselect(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctselect::ctselect(int argc, char *argv[]) : 
-          ctobservation(CTSELECT_NAME, CTSELECT_VERSION, argc, argv)
+          ctobservation(CTSELECT_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -1175,39 +1175,6 @@ std::string ctselect::set_outfile_name(const std::string&amp; filename) const
  
  
 /***********************************************************************//**
- * @brief Get Good Time Intervals extension name
- *
- * @param[in] filename Input file name.
- * @param[in] evtname Events extension name.
- *
- * Extracts the Good Time Intervals extension name from the event file. We
- * do this by loading the events and accessing the Good Time Intervals
- * extension name using the GCTAEventList::gtiname() method. If the file name
- * is empty, the method returns `GTI`.
- ***************************************************************************/
-std::string ctselect::get_gtiname(const std::string& filename,
-                                  const std::string& evtname) const
-{
-    // Initialise GTI name
-    std::string gtiname = gammalib::extname_gti;
-
-    // Continue only if the filename is not empty
-    if (!filename.empty()) {
-
-        // Load events
-        GCTAEventList events(filename+"["+evtname+"]");
-
-        // Get GTI name
-        gtiname = events.gtiname();
-
-    }
-
-    // Return GTI name
-    return (gtiname);
-}
-
-
-/***********************************************************************//**
  * @brief Save event list in FITS format.
  *
  * Save the event list as a FITS file. The file name of the FITS file is
@@ -1325,109 +1292,3 @@ void ctselect::save_xml(void)
     // Return
     return;
 }
-
-
-/***********************************************************************//**
- * @brief Save event list into FITS file
- *
- * @param[in] obs Pointer to CTA observation.
- * @param[in] infile Input file name.
- * @param[in] evtname Event extension name.
- * @param[in] gtiname GTI extension name.
- * @param[in] outfile Output file name.
- *
- * Saves an event list and the corresponding Good Time Intervals into a FITS
- * file and copy all others extensions from the input file to the output
- * file.
- *
- * If an extension name is specified in the @p outfile argument, the events
- * and eventually also the Good Time Intervals will be extracted from the
- * argument and used for writing the events. The format is
- *
- *      <filename>[<event extension name;GTI extension name>]
- *
- * where <filename> needs to be replaced by the name of the FITS file,
- * and <event extension name;GTI extension name> by the name of the events
- * and Good Time Intervals extensions. For example 
- *
- *      myfits.fits[EVENTS1;GTI1]
- *
- * will write the selected events into the "EVENTS1" extension and the
- * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
- * file. If the Good Time Intervals extension name is skipped, e.g.
- *
- *      myfits.fits[EVENTS1]
- *
- * the original extension name for the Good Time Intervals will be kept.
- * Analogously, only the Good Time Intervals extension name can be changed
- * by specifying
- *
- *      myfits.fits[;GTI1]
- *
- * In none of the cases will the original events and Good Time Intervals be
- * copied over to the output file.
- ***************************************************************************/
-void ctselect::save_event_list(const GCTAObservation* obs,
-                               const std::string&     infile,
-                               const std::string&     evtname,
-                               const std::string&     gtiname,
-                               const std::string&     outfile) const
-{
-    // Save only if we have an event list
-    if (obs->eventtype() == "EventList") {
-
-        // Set output FITS file event extension names
-        GFilename   outname(outfile);
-        std::string outevt = evtname;
-        std::string outgti = gtiname;
-        if (outname.has_extname()) {
-            std::vector<std::string> extnames =
-                       gammalib::split(outname.extname(), ";");
-            if (extnames.size() > 0) {
-                std::string extname = gammalib::strip_whitespace(extnames[0]);
-                if (!extname.empty()) {
-                    outevt = extname;
-                }
-            }
-            if (extnames.size() > 1) {
-                std::string extname = gammalib::strip_whitespace(extnames[1]);
-                if (!extname.empty()) {
-                    outgti = extname;
-                }
-            }
-        }
-
-        // Create output FITS file
-        GFits outfits;
-
-        // Write observation into FITS file
-        obs->write(outfits, outevt, outgti);
-
-        // Copy all extensions other than evtname and gtiname extensions
-        // from the input to the output event list. The evtname and
-        // gtiname extensions are written by the save method, all others
-        // that may eventually be present have to be copied over
-        // explicitly.
-        GFits infits(infile);
-        for (int extno = 1; extno < infits.size(); ++extno) {
-            GFitsHDU* hdu = infits.at(extno);
-            if (hdu->extname() != evtname &&
-                hdu->extname() != gtiname &&
-                hdu->extname() != outevt  &&
-                hdu->extname() != outgti) {
-                outfits.append(*hdu);
-            }
-        }
-
-        // Close input file
-        infits.close();
-
-        // Save file to disk and close it (we need both operations)
-        outfits.saveto(outname.url(), clobber());
-        outfits.close();
-
-    } // endif: observation was unbinned
-
-    // Return
-    return;
-}
 /***************************************************************************
  *                      ctselect - Data selection tool                     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -35,8 +35,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTSELECT_NAME    "ctselect"
-#define CTSELECT_VERSION "1.2.0"
+#define CTSELECT_NAME "ctselect"
  
  
 /***********************************************************************//**
@@ -78,15 +77,8 @@ protected:
     std::string check_infile(const std::string& filename,
                              const std::string& evtname) const;
     std::string set_outfile_name(const std::string& filename) const;
-    std::string get_gtiname(const std::string& filename,
-                            const std::string& evtname) const;
     void        save_fits(void);
     void        save_xml(void);
-    void        save_event_list(const GCTAObservation* obs,
-                                const std::string&     infile,
-                                const std::string&     evtname,
-                                const std::string&     gtiname,
-                                const std::string&     outfile) const;
  
     // User parameters
     std::string m_outobs;     //!< Output event list or XML file
@@ -52,7 +52,7 @@
  *
  * Constructs an empty sky mapping tool.
  ***************************************************************************/
-ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION)
+ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -70,7 +70,7 @@ ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION)
  * Constructs sky mapping tool from an observation container.
  ***************************************************************************/
 ctskymap::ctskymap(const GObservations& obs) :
-          ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION, obs)
+          ctobservation(CTSKYMAP_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -91,7 +91,7 @@ ctskymap::ctskymap(const GObservations&amp; obs) :
  * parameter setting.
  ***************************************************************************/
 ctskymap::ctskymap(int argc, char *argv[]) : 
-          ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION, argc, argv)
+          ctobservation(CTSKYMAP_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -32,8 +32,7 @@
 #include "ctobservation.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTSKYMAP_NAME    "ctskymap"
-#define CTSKYMAP_VERSION "1.3.0"
+#define CTSKYMAP_NAME "ctskymap"
  
  
 /***********************************************************************//**
@@ -51,7 +51,7 @@
  *
  * Constructs an empty cttsmap tool.
  ***************************************************************************/
-cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION)
+cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -69,7 +69,7 @@ cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION)
  * Constructs cttsmap tool from an observations container.
  ***************************************************************************/
 cttsmap::cttsmap(const GObservations& obs) :
-         ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, obs)
+         ctlikelihood(CTTSMAP_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ cttsmap::cttsmap(const GObservations&amp; obs) :
  * that are provided as command line arguments.
  ***************************************************************************/
 cttsmap::cttsmap(int argc, char *argv[]) :
-         ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
+         ctlikelihood(CTTSMAP_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTTSMAP_NAME    "cttsmap"
-#define CTTSMAP_VERSION "1.2.0"
+#define CTTSMAP_NAME "cttsmap"
  
  
 /***********************************************************************//**
@@ -53,7 +53,7 @@
  *
  * Constructs an empty ctulimit tool.
  ***************************************************************************/
-ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION)
+ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -71,7 +71,7 @@ ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION)
  * Constructs ctulimit tool from an observations container.
  ***************************************************************************/
 ctulimit::ctulimit(const GObservations& obs) :
-          ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, obs)
+          ctlikelihood(CTULIMIT_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -92,7 +92,7 @@ ctulimit::ctulimit(const GObservations&amp; obs) :
  * parameters that are provided as command line arguments.
  ***************************************************************************/
 ctulimit::ctulimit(int argc, char *argv[]) :
-          ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
+          ctlikelihood(CTULIMIT_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
  
 /* __Definitions _________________________________________________________ */
-#define CTULIMIT_NAME    "ctulimit"
-#define CTULIMIT_VERSION "1.2.0"
+#define CTULIMIT_NAME "ctulimit"
  
  
 /***********************************************************************//**
 /***************************************************************************
  *             ctobservation - Base class for observation tools            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2016 by Juergen Knoedlseder                              *
+ *  copyright (C) 2016-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -326,3 +326,123 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const
     // Return next CTA observation
     return obs;
 }
+
+
+/***********************************************************************//**
+ * @brief Save event list in FITS format.
+ *
+ * Save the event list as a FITS file. The file name of the FITS file is
+ * specified by the `outobs` parameter.
+ ***************************************************************************/
+void ctobservation::save_events_fits(void)
+{
+    // Save only if there are observations
+    if (m_obs.size() > 0) {
+
+        // Get output filename
+        std::string outobs = (*this)["outobs"].filename();
+
+        // Loop over all unbinned CTA observations in the container
+        for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+             obs = next_unbinned_observation()) {
+
+            // Get input file name and default extension names
+            std::string infile  = obs->eventfile();
+            std::string evtname = gammalib::extname_cta_events;
+            std::string gtiname = gammalib::extname_gti;
+
+            // Create file name object
+            GFilename fname(outobs);
+
+            // Extract filename and event extension name
+            std::string outfile = fname.url();
+
+            // Append event extension name. We handle here the possibility
+            // to write the events into a different extension.
+            if (fname.has_extname()) {
+                outfile += "["+fname.extname()+"]";
+            }
+            else {
+                outfile += "["+gammalib::extname_cta_events+"]";
+            }
+
+            // Log filename
+            log_value(NORMAL, "Event list file", outfile);
+            
+            // Save event list
+            save_event_list(obs, infile, evtname, gtiname, outfile);
+
+            // Exit the loop (if this method is called we should only
+            // have a single unbinned observation in the observation
+            // container
+            break;
+
+        } // endfor: looped over all unbinned observations
+
+    } // endif: there were observations
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save event list(s) in XML format.
+ *
+ * Save the event list(s) into FITS files and write the file path information
+ * into a XML file. The filename of the XML file is specified by the `outobs`
+ * parameter, the filename(s) of the event lists are built by prepending a
+ * prefix to the input event list filenames. Any path present in the input
+ * filename will be stripped, i.e. the event list(s) will be written in the
+ * local working directory (unless a path is specified in the prefix).
+ ***************************************************************************/
+void ctobservation::save_events_xml(void)
+{
+    // Get output filename and prefix
+    std::string outobs = (*this)["outobs"].filename();
+
+    // Issue warning if output filename has no .xml suffix
+    log_string(TERSE, warn_xml_suffix(outobs));
+
+    // Loop over all unbinned CTA observations in the container
+    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+         obs = next_unbinned_observation()) {
+
+        // Get input file name and default extension names
+        std::string infile  = obs->eventfile();
+        std::string evtname = gammalib::extname_cta_events;
+        std::string gtiname = gammalib::extname_gti;
+
+        // Extract event and GTI extension names from input FITS file
+        GFilename fname(infile);
+        if (fname.has_extname()) {
+            evtname = fname.extname();
+        }
+        gtiname = get_gtiname(fname.url(), evtname);
+
+        // Set event output file name
+        std::string outfile = set_outfile_name(infile);
+
+        // Append event extension name
+        outfile += "["+evtname+"]";
+
+        // Log filename
+        log_value(NORMAL, "Event list file", outfile);
+
+        // Store output file name in observation
+        obs->eventfile(outfile);
+
+        // Save event list
+        save_event_list(obs, infile, evtname, gtiname, outfile);
+
+    } // endfor: looped over observations
+
+    // Write observation definition XML file name into logger
+    log_value(NORMAL, "Obs. definition file", outobs);
+
+    // Save observations in XML file
+    m_obs.save(outobs);
+
+    // Return
+    return;
+}
 /***************************************************************************
  *             ctobservation - Base class for observation tools            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2016 by Juergen Knoedlseder                              *
+ *  copyright (C) 2016-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -78,6 +78,8 @@ protected:
     void init_members(void);
     void copy_members(const ctobservation& app);
     void free_members(void);
+    void save_events_fits(void);
+    void save_events_xml(void);
  
     // Protected members
     GObservations m_obs;            //!< Observation container
@@ -1028,6 +1028,41 @@ GCTARoi ctool::get_roi(void)
  
  
 /***********************************************************************//**
+ * @brief Set output file name.
+ *
+ * @param[in] filename Input file name.
+ *
+ * Required user parameters:
+ *      prefix
+ *
+ * Converts an input file name into an output filename by prepending the
+ * prefix specified by the `prefix` parameter to the input file name. Any
+ * path as well as extension will be stripped from the input file name. Also
+ * a trailing `.gz` will be stripped as one cannot write into gzipped files.
+ ***************************************************************************/
+std::string ctool::set_outfile_name(const std::string& filename)
+{
+    // Get prefix User parameter
+    std::string prefix = (*this)["prefix"].string();
+
+    // Create filename
+    GFilename fname(filename);
+
+    // Split input filename without any extensions into path elements
+    std::vector<std::string> elements = gammalib::split(fname.url(), "/");
+
+    // The last path element is the filename
+    std::string outname = prefix + elements[elements.size()-1];
+
+    // Strip any ".gz"
+    outname = gammalib::strip_chars(outname, ".gz");
+    
+    // Return output filename
+    return outname;
+}
+
+
+/***********************************************************************//**
  * @brief Set response for all CTA observations in container
  *
  * @param[in,out] obs Observation container
@@ -1730,6 +1765,145 @@ std::vector&lt;bool&gt; ctool::cube_layer_usage(const GEbounds&amp; cube_ebounds,
  
  
 /***********************************************************************//**
+ * @brief Save event list into FITS file
+ *
+ * @param[in] obs Pointer to CTA observation.
+ * @param[in] infile Input file name.
+ * @param[in] evtname Event extension name.
+ * @param[in] gtiname GTI extension name.
+ * @param[in] outfile Output file name.
+ *
+ * Saves an event list and the corresponding Good Time Intervals into a FITS
+ * file and copy all others extensions from the input file to the output
+ * file.
+ *
+ * If an extension name is specified in the @p outfile argument, the events
+ * and eventually also the Good Time Intervals will be extracted from the
+ * argument and used for writing the events. The format is
+ *
+ *      <filename>[<event extension name;GTI extension name>]
+ *
+ * where <filename> needs to be replaced by the name of the FITS file,
+ * and <event extension name;GTI extension name> by the name of the events
+ * and Good Time Intervals extensions. For example 
+ *
+ *      myfits.fits[EVENTS1;GTI1]
+ *
+ * will write the selected events into the `EVENTS1` extension and the
+ * Good Time Intervals into the `GTI1` extension of the `myfits.fits` FITS
+ * file. If the Good Time Intervals extension name is skipped, e.g.
+ *
+ *      myfits.fits[EVENTS1]
+ *
+ * the original extension name for the Good Time Intervals will be kept.
+ * Analogously, only the Good Time Intervals extension name can be changed
+ * by specifying
+ *
+ *      myfits.fits[;GTI1]
+ *
+ * In none of the cases will the original events and Good Time Intervals be
+ * copied over to the output file.
+ ***************************************************************************/
+void ctool::save_event_list(const GCTAObservation* obs,
+                            const std::string&     infile,
+                            const std::string&     evtname,
+                            const std::string&     gtiname,
+                            const std::string&     outfile) const
+{
+    // Save only if we have an event list
+    if (obs->eventtype() == "EventList") {
+
+        // Set output FITS file event extension names
+        GFilename   outname(outfile);
+        std::string outevt = evtname;
+        std::string outgti = gtiname;
+        if (outname.has_extname()) {
+            std::vector<std::string> extnames =
+                       gammalib::split(outname.extname(), ";");
+            if (extnames.size() > 0) {
+                std::string extname = gammalib::strip_whitespace(extnames[0]);
+                if (!extname.empty()) {
+                    outevt = extname;
+                }
+            }
+            if (extnames.size() > 1) {
+                std::string extname = gammalib::strip_whitespace(extnames[1]);
+                if (!extname.empty()) {
+                    outgti = extname;
+                }
+            }
+        }
+
+        // Create output FITS file
+        GFits outfits;
+
+        // Write observation into FITS file
+        obs->write(outfits, outevt, outgti);
+
+        // Copy all extensions other than evtname and gtiname extensions
+        // from the input to the output event list. The evtname and
+        // gtiname extensions are written by the save method, all others
+        // that may eventually be present have to be copied over
+        // explicitly.
+        GFits infits(infile);
+        for (int extno = 1; extno < infits.size(); ++extno) {
+            GFitsHDU* hdu = infits.at(extno);
+            if (hdu->extname() != evtname &&
+                hdu->extname() != gtiname &&
+                hdu->extname() != outevt  &&
+                hdu->extname() != outgti) {
+                outfits.append(*hdu);
+            }
+        }
+
+        // Close input file
+        infits.close();
+
+        // Save file to disk and close it (we need both operations)
+        outfits.saveto(outname.url(), clobber());
+        outfits.close();
+
+    } // endif: observation was unbinned
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get Good Time Intervals extension name
+ *
+ * @param[in] filename Input file name.
+ * @param[in] evtname Events extension name.
+ *
+ * Extracts the Good Time Intervals extension name from the event file. We
+ * do this by loading the events and accessing the Good Time Intervals
+ * extension name using the GCTAEventList::gtiname() method. If the file name
+ * is empty, the method returns `GTI`.
+ ***************************************************************************/
+std::string ctool::get_gtiname(const std::string& filename,
+                               const std::string& evtname) const
+{
+    // Initialise GTI name
+    std::string gtiname = gammalib::extname_gti;
+
+    // Continue only if the filename is not empty
+    if (!filename.empty()) {
+
+        // Load events
+        GCTAEventList events(filename+"["+evtname+"]");
+
+        // Get GTI name
+        gtiname = events.gtiname();
+
+    }
+
+    // Return GTI name
+    return (gtiname);
+}
+
+
+/***********************************************************************//**
  * @brief Set warning string if there are too few energies
  *
  * @param[in] energies Energies.
 /***************************************************************************
  *                        ctool - ctool base class                         *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2014-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -104,6 +104,7 @@ protected:
  
     // Protected methods that extract user parameters
     GCTARoi           get_roi(void);
+    std::string       set_outfile_name(const std::string& filename);
  
     // Protected support methods
     void              set_response(GObservations& obs);
@@ -121,6 +122,13 @@ protected:
     std::vector<bool> cube_layer_usage(const GEbounds& cube_ebounds,
                                        const GEbounds& list_ebounds) const;
     bool              is_valid_filename(const GFilename& filename) const;
+    std::string       get_gtiname(const std::string& filename,
+                                  const std::string& evtname) const;
+    void              save_event_list(const GCTAObservation* obs,
+                                      const std::string&     infile,
+                                      const std::string&     evtname,
+                                      const std::string&     gtiname,
+                                      const std::string&     outfile) const;
  
     // Protected warning strings
     std::string     warn_too_few_energies(const GEnergies& energies) const;
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
+      <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
+    </spatialModel>
+    <temporal type="PhaseCurve" file="model_temporal_phasecurve.fits">
+      <parameter name="Normalization" scale="1" value="1.0"     min="0.0" max="1000.0"   free="0"/>
+      <parameter name="MJD"           scale="1" value="51544.5" min="0.0" max="100000.0" free="0"/>
+      <parameter name="Phase"         scale="1" value="0.0"     min="0.0" max="1.0"      free="0"/>
+      <parameter name="F0"            scale="1" value="1.0"     min="0.0" max="1000.0"   free="0"/>
+      <parameter name="F1"            scale="1" value="0.1"     min="0.0" max="1000.0"   free="0"/>
+      <parameter name="F2"            scale="1" value="0.01"    min="0.0" max="1000.0"   free="0"/>
+    </temporal>
+  </source>
+</source_library>
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the cssens script
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 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
@@ -168,7 +168,7 @@ class Test(test):
         # Check dimensions
         self.test_value(results.nrows(), 2,
              'Check number of rows in sensitivity file')
-        self.test_value(results.ncols(), 7,
+        self.test_value(results.ncols(), 10,
              'Check number of columns in sensitivity file')
  
         # Return
@@ -168,6 +168,9 @@ class Test(test):
         cpy_butterfly['outfile'] = 'ctbutterfly_py3.dat'
         cpy_butterfly['logfile'] = 'ctbutterfly_py3.log'
         cpy_butterfly['chatter'] = 4
+
+        # Execute ctbutterfly tool
+        cpy_butterfly.logFileOpen()   # Make sure we get a log file
         cpy_butterfly.execute()
  
         # Check result file
 #! /usr/bin/env python
 # ==========================================================================
-# This scripts performs unit tests for the ctselect tool.
+# This scripts performs unit tests for the ctphase tool.
 #
-# Copyright (C) 2014-2016 Juergen Knoedlseder
+# Copyright (C) 2017 by Joshua Cardenzana
 #
 # 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,8 @@ class Test(test):
         test.__init__(self)
  
         # Set test data
+        self._model_file     = self._datadir + '/model_temporal_phasecurve.xml'
+        self._phase_file     = self._datadir + '/model_temporal_phasecurve.fits'
         self._invalid_events = self._datadir + '/invalid_event_list.fits'
  
         # Return
@@ -70,10 +72,9 @@ class Test(test):
  
         # Setup ctphase command
         cmd = ctphase+' inobs="'+self._events+'"'+ \
-                       ' outobs="ctphase_cmd1.fits"'+ \
-                       ' p0=0.1 mjd=51544.5'+ \
-                       ' f0=0.1 f1=0. f2=0.'+ \
-                       ' logfile="ctphase_cmd1.log" chatter=1'
+                      ' outobs="ctphase_cmd1.fits"'+ \
+                      ' inmodel="'+self._model_file+'" source="Crab"'+ \
+                      ' logfile="ctphase_cmd1.log" chatter=1'
  
         # Check if execution of wrong command fails
         self.test_assert(self._execute('command_that_does_not_exist') != 0,
@@ -88,10 +89,9 @@ class Test(test):
  
         # Setup ctphase command
         cmd = ctphase+' inobs="event_file_that_does_not_exist.fits"'+ \
-                       ' outobs="ctphase_cmd2.fits"'+ \
-                       ' p0=0.1 mjd=51544.5'+ \
-                       ' f0=0.1 f1=0. f2=0.'+ \
-                       ' logfile="ctphase_cmd2.log" chatter=1'
+                      ' outobs="ctphase_cmd2.fits"'+ \
+                      ' inmodel="'+self._model_file+'" source="Crab"'+ \
+                      ' logfile="ctphase_cmd2.log" chatter=1'
  
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
@@ -124,6 +124,8 @@ class Test(test):
         self._check_obs(phase.obs(), nobs=0)
  
         # Check that saving does nothing
+        phase['inmodel'] = self._model_file
+        phase['source']  = 'Crab'
         phase['outobs']  = 'ctphase_py0.fits'
         phase['logfile'] = 'ctphase_py0.log'
         phase.logFileOpen()
@@ -136,11 +138,8 @@ class Test(test):
  
         # Now set ctphase parameters
         phase['inobs']   = self._events
-        phase['p0']      = 0.
-        phase['mjd']     = 51544.5
-        phase['f0']      = 1.
-        phase['f1']      = 0.0
-        phase['f2']      = 0.
+        phase['inmodel'] = self._model_file
+        phase['source']  = 'Crab'
         phase['outobs']  = 'ctphase_py1.fits'
         phase['logfile'] = 'ctphase_py1.log'
         phase['chatter'] = 2
@@ -151,34 +150,50 @@ class Test(test):
         phase.save()
  
         # Check result file
-        self._check_result_file('ctphase_py1.fits', nevents=6141)
+        self._check_result_file('ctphase_py1.fits')
+
+        # Now run ctphase tool without a model file
+        phase['inmodel']    = 'NONE'
+        phase['source']     = 'NONE'
+        phase['phasecurve'] = self._phase_file
+        phase['mjd']        = 51544.5
+        phase['phase']      = 0.0
+        phase['f0']         = 1.0
+        phase['f1']         = 0.1
+        phase['f2']         = 0.01
+        phase['outobs']     = 'ctphase_py2.fits'
+        phase['logfile']    = 'ctphase_py2.log'
+        phase['chatter']    = 3
+
+        # Execute ctphase tool
+        phase.logFileOpen()   # Make sure we get a log file
+        phase.execute()
+
+        # Check result file
+        self._check_result_file('ctphase_py2.fits')
  
         # Copy ctphase tool
         cpy_phase = phase.copy()
  
-        # Execute copy of ctphase tool again, now with a negative value
-        # of initial phase to test for invalid parameter value
-        self.test_try('Test invalid parameter value')
-        try:
-            cpy_phase['mjd']     = 51544.5
-            cpy_phase['f0']      = 1.
-            cpy_phase['f1']      = 0.2
-            cpy_phase['f2']      = 0.
-            cpy_phase['outobs']  = 'ctphase_py2.fits'
-            cpy_phase['logfile'] = 'ctphase_py2.log'
-            cpy_phase['chatter'] = 3
-            cpy_phase['publish'] = True
-            cpy_phase.logFileOpen()  # Needed to get a new log file
-            cpy_phase['p0']      = -0.5
-            cpy_phase.execute()
-            self.test_try_failure('Exception not thrown')
-        except ValueError:
-            self.test_try_success()
+        # Clear phase tool
+        phase.clear()
+
+        # Execute copy
+        cpy_phase['outobs']  = 'ctphase_py3.fits'
+        cpy_phase['logfile'] = 'ctphase_py3.log'
+        cpy_phase['chatter'] = 4
+        cpy_phase['publish'] = True
+        cpy_phase.logFileOpen()  # Needed to get a new log file
+        cpy_phase.execute()
  
+        # Check result file
+        self._check_result_file('ctphase_py3.fits')
+
+        # Return
         return
  
     # Check result file
-    def _check_result_file(self, filename, nevents=591, dec=22.01, rad=3.0):
+    def _check_result_file(self, filename, nevents=6141):
         """
         Check result file
  
@@ -188,14 +203,10 @@ class Test(test):
             Event list file name
         nevents : int, optional
             Expected number of events
-        dec : float, optional
-            Expected Declination (deg)
-        rad : float, optional
-            Expected radius (deg)
         """
         # Open result file
-        events = gammalib.GCTAEventList(filename)
-        phases = [ events[i].phase() for i in range(events.size())]
+        events      = gammalib.GCTAEventList(filename)
+        phases      = [events[i].phase() for i in range(events.size())]
         test_result = int(sum(phases)>0)
         self.test_value(test_result, 1, 'Check number of events')
  
@@ -213,8 +224,6 @@ class Test(test):
             Models
         nobs : int, optional
             Expected number of observations
-        nevents : int, optional
-            Expected number of events
         """
         # Check number of observations
         self.test_value(obs.size(), nobs, 'Check number of observations')
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the ctprob tool.
+#
+# Copyright (C) 2017 by Leonardo Di Venere
+#
+# 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
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import os
+import gammalib
+import ctools
+from testing import test
+
+
+# ============================ #
+# Test class for ctprob tool #
+# ============================ #
+class Test(test):
+    """
+    Test class for ctprob tool
+    """
+    # Constructor
+    def __init__(self):
+        """
+        Constructor
+        """
+        # Call base class constructor
+        test.__init__(self)
+
+        # Set test data
+        self._invalid_events = self._datadir + '/invalid_event_list.fits'
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions
+        """
+        # Set test name
+        self.name('ctprob')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test ctprob on command line')
+        self.append(self._test_python, 'Test ctprob from Python')
+
+        # Return
+        return
+
+    # Test ctprob on command line
+    def _test_cmd(self):
+        """
+        Test ctprob on the command line
+        """
+        # Set tool name
+        ctprob = self._tool('ctprob')
+
+        # Setup ctprob command
+        cmd = ctprob+' inobs="'+self._events+'"'+ \
+                     ' outobs="ctprob_cmd1.fits"'+ \
+                     ' inmodel='+self._model+ \
+                     ' caldb='+self._caldb+' irf='+self._irf+ \
+                     ' edisp=no'+ \
+                     ' logfile="ctprob_cmd1.log" chatter=1'
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Check result file
+        self._check_result_file('ctprob_cmd1.fits')
+
+        # Setup ctprob command
+        cmd = ctprob+' inobs="event_file_that_does_not_exist.fits"'+ \
+                     ' outobs="ctprob_cmd2.fits"'+ \
+                     ' inmodel='+self._model+ \
+                     ' caldb='+self._caldb+' irf='+self._irf+ \
+                     ' edisp=no'+ \
+                     ' logfile="ctprob_cmd2.log" chatter=1'
+
+        # Check if execution failed
+        self.test_assert(self._execute(cmd) != 0,
+             'Check invalid input file when executed from command line')
+
+        # Setup ctprob --help option
+        cmd = ctprob+' --help'
+
+        # Check if execution was successful in case that the CTOOLS
+        # environment variable was set or failed otherwise
+        if 'CTOOLS' in os.environ:
+            self.test_value(self._execute(cmd), 0,
+                 'Check successful execution with --help option')
+        else:
+            self.test_assert(self._execute(cmd) != 0,
+                 'Check execution failure with --help option')
+
+        # Return
+        return
+
+    # Test ctprob from Python
+    def _test_python(self):
+        """
+        Test ctprob from Python
+        """
+        # Allocate empty ctprob tool
+        prob = ctools.ctprob()
+
+        # Check that empty ctprob tool holds an empty observation
+        self._check_obs(prob.obs(), nobs=0)
+
+        # Check that saving does nothing
+        prob['outobs']  = 'ctprob_py0.fits'
+        prob['logfile'] = 'ctprob_py0.log'
+        prob.logFileOpen()
+        prob.save()
+        self.test_assert(not os.path.isfile('ctprob_py0.fits'),
+             'Check that no event list has been created')
+
+        # Check that clearing does not lead to an exception or segfault
+        prob.clear()
+
+        # Now set ctprob parameters
+        prob['inobs']   = self._events
+        prob['outobs']  = 'ctprob_py1.fits'
+        prob['inmodel'] = self._model
+        prob['caldb']   = self._caldb
+        prob['irf']     = self._irf
+        prob['edisp']   = False
+        prob['logfile'] = 'ctprob_py1.log'
+        prob['chatter'] = 2
+
+        # Run ctprob tool
+        prob.logFileOpen()   # Make sure we get a log file
+        prob.run()
+        prob.save()
+
+        # Check result file
+        self._check_result_file('ctprob_py1.fits')
+
+        # Setup ctprob tool for an invalid event file
+        prob = ctools.ctprob()
+        prob['inobs']   = self._invalid_events
+        prob['outobs']  = 'ctprob_py2.fits'
+        prob['inmodel'] = self._model
+        prob['caldb']   = self._caldb
+        prob['irf']     = self._irf
+        prob['edisp']   = False
+        prob['logfile'] = 'ctprob_py2.log'
+        prob['chatter'] = 3
+
+        # Execute tool
+        prob.logFileOpen()  # Needed to get a new log file
+        self.test_try('Test invalid event file')
+        try:
+            prob.execute()
+            self.test_try_failure('Exception not thrown')
+        except RuntimeError:
+            self.test_try_success()
+
+        # Setup ctprob with a list of observations and test prefix
+        prob = ctools.ctprob()
+        prob['inobs']   = self._datadir + '/obs_unbinned.xml'
+        prob['prefix']  = 'ctprob_py3_'
+        prob['outobs']  = 'ctprob_py3_obs.xml'
+        prob['inmodel'] = self._model
+        prob['caldb']   = self._caldb
+        prob['irf']     = self._irf
+        prob['edisp']   = True
+        prob['logfile'] = 'ctprob_py3.log'
+        prob['chatter'] = 4
+
+        # Run ctprob tool
+        prob.logFileOpen()   # Make sure we get a log file
+        prob.execute()
+
+        # Check result file
+        self._check_result_file('ctprob_py3_crab_events.fits')
+
+        # Return
+        return
+
+    # Check result file
+    def _check_result_file(self, filename, nevents=6141):
+        """
+        Check result file
+
+        Parameters
+        ----------
+        filename : str
+            Event list file name
+        nevents : int, optional
+            Expected number of events
+        """
+        # Open result file
+        events = gammalib.GCTAEventList(filename)
+
+        # Check event list
+        self._check_events(events, nevents=nevents)
+
+        # Check that probability columns exist
+        colname1 = 'PROB_Background'
+        colname2 = 'PROB_Crab'
+        self._check_column(filename, colname1)
+        self._check_column(filename, colname2)        
+
+        # Check that the probability values are correctly normalized
+        self._check_normalization(filename, [colname1,colname2])
+
+        # Return
+        return
+
+    # Check observation and event list
+    def _check_obs(self, obs, nobs=1, nevents=6141):
+        """
+        Check observation and event list
+
+        Parameters
+        ----------
+        obs : `~gammalib.GObservations`
+            Models
+        nobs : int, optional
+            Expected number of observations
+        nevents : int, optional
+            Expected number of events
+        """
+        # Check number of observations
+        self.test_value(obs.size(), nobs, 'Check number of observations')
+
+        # If there is an observation then check the event list of the first
+        # one
+        if obs.size() > 0:
+            self._check_events(obs[0].events(), nevents=nevents)
+
+        # Return
+        return
+
+    # Check events
+    def _check_events(self, events, nevents=6141):
+        """
+        Check event list
+
+        Parameters
+        ----------
+        events : `~gammalib.GCTAEventList`
+            Event list
+        nevents : int, optional
+            Expected number of events
+        """
+        # Check event list
+        self.test_value(events.size(), nevents, 'Check number of events')
+
+        # Return
+        return
+
+
+    def _check_column(self, filename, colname):
+        """
+        Check that column colname exists.
+
+        Parameters
+        ----------
+        filename : str
+            Event list file name
+        colname : str
+            Name of the column to check
+        """
+        # Check that column exists
+        gfits = gammalib.GFits(filename)
+        gtable = gfits[1]
+        check = int(gtable.contains(colname))
+        self.test_value(check, 1, 'Check column '+colname)
+
+        # Return
+        return
+
+
+    def _check_normalization(self, filename, colnames):
+        """
+        Check that the probability values in columns identified by colnames are
+        correctly normalized.
+
+        Parameters
+        ----------
+        filename : str
+            Event list file name
+        colnames : list
+            Names of the columns to check
+        """
+        # Check that the probability values are correctly normalized
+
+        # Open fits table
+        gfits  = gammalib.GFits(filename)
+        gtable = gfits[1]
+        ncols  = gtable.ncols()
+        nevt   = gtable.nrows()
+        srcs   = []
+
+        # Loop over colnames
+        for colname in colnames:
+
+            # Get the correct column corresponding to colname from the table
+            for icol in range(ncols):
+                gcol = gtable[icol] 
+                if gcol.name()==colname:
+                    break
+
+            # Read and save column content 
+            src1 = []
+            for ievt in range(nevt):
+                src1.append(gcol[ievt])
+            srcs.extend(src1)
+
+        # Compute sum
+        tot = sum(srcs)
+        self.test_value(int(tot+0.5), nevt, 'Check that probability columns '
+                                            'are normalized')
+
+        # Return
+        return
+
@@ -40,6 +40,7 @@ import test_ctcubemask
 import test_ctbutterfly
 import test_ctulimit
 import test_cterror
+import test_ctprob
 import test_pipelines
  
  
@@ -118,6 +119,7 @@ def test(installed=False):
              test_ctbutterfly.Test(),
              test_ctulimit.Test(),
              test_cterror.Test(),
+             test_ctprob.Test(),
              test_pipelines.Test()]
  
     # Allocate test suite container