Add csphasecrv unit test, use global ctools version and update parameter files

Jürgen Knödlseder
1 parent 1110a353
Showing 64 changed files with 644 additions and 889 deletions
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/cscaldb.py
cscripts/csebins.py
cscripts/csfindobs.py
cscripts/csiactcopy.py
cscripts/csiactdata.py
cscripts/csiactobs.py
cscripts/cslightcrv.par
cscripts/cslightcrv.py
cscripts/csmodelinfo.py
cscripts/csmodelmerge.py
cscripts/csmodelselect.py
cscripts/csobs2caldb.py
cscripts/csobsdef.py
cscripts/csobsinfo.py
cscripts/csobsselect.py
cscripts/csphasecrv.par
cscripts/csphasecrv.py
cscripts/cspull.par
@@ -32,6 +32,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
                 $(srcdir)/csobsdef.py \
                 $(srcdir)/csobsinfo.py \
                 $(srcdir)/csobsselect.py \
+                $(srcdir)/csphasecrv.py \
                 $(srcdir)/cspull.py \
                 $(srcdir)/csresmap.py \
                 $(srcdir)/csroot2caldb.py  \
@@ -63,6 +64,7 @@ parfiles = $(srcdir)/cscaldb.par \
            $(srcdir)/csobsdef.par \
            $(srcdir)/csobsinfo.par \
            $(srcdir)/csobsselect.par \
+           $(srcdir)/csphasecrv.par \
            $(srcdir)/cspull.par \
            $(srcdir)/csresmap.par \
            $(srcdir)/csroot2caldb.par \
@@ -91,6 +93,7 @@ linkscripts = cscaldb \
               csobsdef \
               csobsinfo \
               csobsselect \
+              csphasecrv \
               cspull \
               csresmap \
               csroot2caldb \
@@ -183,7 +186,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/iact_runlist_1.dat \
            $(top_srcdir)/test/data/iact_runlist_2.dat \
            $(top_srcdir)/test/data/model_temporal_phasecurve.xml \
-           $(top_srcdir)/test/data/model_temporal_phasecurve.fits
+           $(top_srcdir)/test/data/model_temporal_phasecurve.fits \
+           $(top_srcdir)/test/data/phased_events.fits.gz
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata
 testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
@@ -35,6 +35,7 @@ __all__ = [
     "csobsdef",
     "csobsinfo",
     "csobsselect",
+    "csphasecrv",
     "cspull",
     "csresmap",
     "csroot2caldb",
@@ -69,6 +70,7 @@ from cscripts.csobs2caldb   import csobs2caldb
 from cscripts.csobsdef      import csobsdef
 from cscripts.csobsinfo     import csobsinfo
 from cscripts.csobsselect   import csobsselect
+from cscripts.csphasecrv    import csphasecrv
 from cscripts.cspull        import cspull
 from cscripts.csresmap      import csresmap
 from cscripts.csroot2caldb  import csroot2caldb
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Show the content of the ctools calibration database
 #
-# Copyright (C) 2014-2016 Juergen Knoedlseder
+# Copyright (C) 2014-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
@@ -55,7 +55,7 @@ class cscaldb(ctools.cscript):
         """
         # Set name
         self._name    = 'cscaldb'
-        self._version = '1.1.0'
+        self._version = ctools.__version__
         # Initialise application by calling the appropriate class
         # constructor.
@@ -40,7 +40,7 @@ class csebins(ctools.cscript):
         """
         # Set name and version
         self._name    = 'csebins'
-        self._version = '1.0.0'
+        self._version = ctools.__version__
         # Set members
         self._ebounds = gammalib.GEbounds()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Find observations from an IACT data store
 #
-# Copyright (C) 2016 Michael Mayer
+# Copyright (C) 2016-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -40,7 +40,7 @@ class csfindobs(ctools.cscript):
         """
         # Set name
         self._name         = 'csfindobs'
-        self._version      = '1.2.0'
+        self._version      = ctools.__version__
         self._datapath     = os.getenv('VHEFITS','')
         self._prodname     = ''
         self._select_radec = True
@@ -44,7 +44,7 @@ class csiactcopy(ctools.cscript):
         """
         # Set name
         self._name          = 'csiactcopy'
-        self._version       = '1.2.0'
+        self._version       = ctools.__version__
         self._datapath      = os.getenv('VHEFITS','')
         self._remote_master = gammalib.GFilename()
         self._remote_base   = ''
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Inspect an IACT data storage
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -45,7 +45,7 @@ class csiactdata(ctools.cscript):
         """
         # Set name
         self._name        = 'csiactdata'
-        self._version     = '1.1.0'
+        self._version     = ctools.__version__
         self._datapath    = os.getenv('VHEFITS','')
         self._prodnames   = []
         self._master_indx = ''
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates an IACT observation definition XML file.
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -47,7 +47,7 @@ class csiactobs(ctools.cscript):
         """
         # Set name and version
         self._name    = 'csiactobs'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise some members
         self._obs              = gammalib.GObservations()
@@ -36,7 +36,7 @@ emin,      r, a, 0.1,,, &quot;Lower energy limit of events (TeV)&quot;
 emax,      r, a, 100.0,,, "Upper energy limit of events (TeV)"
 enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)"
 coordsys,  s, a, CEL,CEL|GAL,, "Coordinate System"
-proj,      s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+proj,      s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
@@ -69,7 +69,7 @@ class cslightcrv(ctools.cscript):
         """
         # Set name
         self._name    = 'cslightcrv'
-        self._version = '1.3.0'
+        self._version = ctools.__version__
         # Initialise some members
         self._srcname = ''
@@ -509,11 +509,7 @@ class cslightcrv(ctools.cscript):
         self._get_parameters()
         # Write observation into logger
-        if self._logTerse():
-            self._log('\n')
-            self._log.header1(gammalib.number('Input observation',len(self._obs)))
-            self._log(str(self._obs))
-            self._log('\n')
+        self._log_observations(gammalib.NORMAL, self._obs, 'Observation')
         # Get time boundaries
         tmin = self._tbins.tstart(0)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Shows the content of a model container
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -44,7 +44,7 @@ class csmodelinfo(ctools.cscript):
         """
         # Set name
         self._name    = 'csmodelinfo'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._models = gammalib.GModels()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Merge model definition XML files
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -43,7 +43,7 @@ class csmodelmerge(ctools.cscript):
         """
         # Set name
         self._name    = 'csmodelmerge'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._files      = None
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Selects model from a model definition XML file
 #
-# 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
@@ -38,7 +38,7 @@ class csmodelselect(ctools.cscript):
         """
         # Set name
         self._name    = 'csmodelselect'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise observation container from constructor arguments
         self._obs, argv = self._set_input_obs(argv)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Creates a calibration database entry for an IACT observation.
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -48,7 +48,7 @@ class csobs2caldb(ctools.cscript):
         """
         # Set name and version
         self._name    = 'csobs2caldb'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise members
         self._observation = gammalib.GCTAObservation()
@@ -86,7 +86,7 @@ class csobsdef(ctools.cscript):
         """
         # Set name and version
         self._name    = 'csobsdef'
-        self._version = '1.3.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._obs    = gammalib.GObservations()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Dump information about observation into log file
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -38,7 +38,7 @@ class csobsinfo(ctools.cscript):
         """
         # Set name
         self._name    = 'csobsinfo'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._obj_dir        = None
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Selects observations from an observation definition XML file
 #
-# 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
@@ -38,7 +38,7 @@ class csobsselect(ctools.cscript):
         """
         # Set name
         self._name    = 'csobsselect'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise observation container from constructor arguments
         self._obs, argv = self._set_input_obs(argv)
@@ -27,14 +27,14 @@ outfile,   f, a, phasecrv.fits,,, &quot;Output fits file name.&quot;
 #
 # Script parameters
 #==================
-tbinalg,   s, a, LIN,LIN|FILE,, "Algorithm for defining time bins"
+phbinalg,  s, a, LIN,LIN|FILE,, "Algorithm for defining phase bins"
 phbins,    i, a, 5,1,10000, "Number of phase bins"
 phbinfile, f, a, NONE,,, "File defining the phase binning"
 emin,      r, a, 0.1,,, "Lower energy limit of events (TeV)"
 emax,      r, a, 100.0,,, "Upper energy limit of events (TeV)"
 enumbins,  i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)"
 coordsys,  s, a, CEL,CEL|GAL,, "Coordinate System"
-proj,      s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+proj,      s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
@@ -44,6 +44,7 @@ binsz,     r, a, 0.02,,, &quot;Pixel size (deg/pixel)&quot;
 #
 # Standard parameters
 #====================
+publish, b, h, no,,, "Publish phase curve 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"
 #! /usr/bin/env python
 # ==========================================================================
-# Generates pulsar spectra as a function of phase
+# Generates spectra as a function of phase
 #
-# Copyright (C) 2014-2017 Rolf Buehler
+# Copyright (C) 2017 Rolf Buehler
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -29,22 +29,22 @@ from cscripts import obsutils
 # ================ #
 class csphasecrv(ctools.cscript):
     """
-    Generates a spectra in phase bins
+    Generates spectra in phase bins
     This script computes spectra by performing a maximum likelihood fit
-    using :doc:`ctlike` in a series of phase bins for pulsars.
-    The phase bins can be either specified in an ASCII file, as an interval
-    divided into equally sized phase bins. 
+    using :doc:`ctlike` in a series of phase bins for pulsars. The phase bins
+    can be either specified in an ASCII file or as an interval divided into
+    equally sized phase bins.
     Examples:
-            >>> phrv = csphasecrv()
-            >>> phrv.run()
+            >>> phcrv = csphasecrv()
+            >>> phcrv.run()
             >>> ... (querying for parameters) ...
-            >>> phrv = phrv.phasecrv()
+            >>> phcrv = phrv.phasecrv()
                 Generates phase fits and retrieves dictionary with best fit models.
-            >>> lcrv = csphasecrv(obs)
-            >>> lcrv.execute()
+            >>> phcrv = csphasecrv(obs)
+            >>> phcrv.execute()
             >>> ... (querying for parameters) ...
                 Generates phase fits from the observations 
                 and saves results in XML files.
@@ -57,13 +57,14 @@ class csphasecrv(ctools.cscript):
         """
         # Set name
         self._name    = 'csphasecrv'
-        self._version = '1.0.0'
+        self._version = ctools.__version__
-        # Initialise some members
+        # Initialise some members. Phases are stored in a nested list
+        # [[ph1min,ph1max], [ph2min,ph2max],..]
         self._srcname = ''
-        # Phases are stored in a nested list [[ph1min,ph1max], [ph2min,ph2max],..]
-        self._phbins  = [[0.0,1.0]]   
+        self._phbins  = [[0.0,1.0]]
         self._stacked = False
+        self._fits    = gammalib.GFits()
         # Initialise observation container from constructor arguments.
         self._obs, argv = self._set_input_obs(argv)
@@ -144,7 +145,7 @@ class csphasecrv(ctools.cscript):
         # Get algorithm to use for defining the time intervals. Possible
         # values are "FILE" or "LIN". This is enforced at
         # parameter file level, hence no checking is needed.
-        algorithm = self['tbinalg'].string()
+        algorithm = self['phbinalg'].string()
         # If the algorithm is "FILE" then handle a FITS or an ASCII file for
         # the time bin definition
@@ -173,8 +174,8 @@ class csphasecrv(ctools.cscript):
             # Compute time step and setup the GTIs
             ph_step = 1.0 / float(nbins)
             for i in range(nbins):
-                phmin = i *ph_step
-                phmax = (i+1)*ph_step
+                phmin = i     * ph_step
+                phmax = (i+1) * ph_step
                 phbins.append([phmin,phmax])
         # Return Good Time Intervals
@@ -252,6 +253,168 @@ class csphasecrv(ctools.cscript):
         # Return new oberservation container
         return new_obs
+    def _get_free_par_names(self):
+        """
+        Return list of free parameter names
+
+        Returns
+        -------
+        names : list of str
+            List of free parameter names.
+        """
+        # Initialise empty list of free parameter names
+        names = []
+
+        # Collect list of free parameter names
+        for par in self._obs.models()[self._srcname]:
+            if par.is_free():
+                names.append(par.name())
+
+        # Return names
+        return names
+
+    def _create_fits_table(self, results):
+        """
+        Creates FITS binary table containing light curve results
+
+        Parameters
+        ----------
+        results : list of dict
+            List of result dictionaries
+
+        Returns
+        -------
+        table : `~gammalib.GFitsBinTable`
+            FITS binary table containing light curve
+        """
+        # Determine number of rows in FITS table
+        nrows = len(results)
+
+        # Create FITS Table with extension "PHASECURVE"
+        table = gammalib.GFitsBinTable(nrows)
+        table.extname('PHASECURVE')
+
+        # Append phase columns "PHASE_MIN" and "PHASE_MAX"
+        phmin = gammalib.GFitsTableDoubleCol('PHASE_MIN', nrows)
+        phmax = gammalib.GFitsTableDoubleCol('PHASE_MAX', nrows)
+        for i, result in enumerate(results):
+            phmin[i] = result['phmin']
+            phmax[i] = result['phmax']
+        table.append(phmin)
+        table.append(phmax)
+
+        # Create parameter columns
+        for par in self._obs.models()[self._srcname]:
+            if par.is_free():
+                name   = par.name()
+                e_name = 'e_'+par.name()
+                col    = gammalib.GFitsTableDoubleCol(name, nrows)
+                e_col  = gammalib.GFitsTableDoubleCol(e_name, nrows)
+                col.unit(par.unit())
+                e_col.unit(par.unit())
+                for i, result in enumerate(results):
+                    col[i]   = result['values'][name]
+                    e_col[i] = result['values'][e_name]
+                table.append(col)
+                table.append(e_col)
+
+        # Return table
+        return table
+
+    def _create_fits(self):
+        """
+        Create FITS file object from fit results
+        """
+        # Initialise list of result dictionaries
+        results = []
+
+        # Get source parameters
+        pars = self._get_free_par_names()
+
+        # Loop over time bins
+        for i in range(len(self._phbins)):
+
+            # Get time boundaries
+            phmin = self._phbins[i][0]
+            phmax = self._phbins[i][1]
+
+            # Initialise result dictionary
+            result = {'phmin': phmin,
+                      'phmax': phmax,
+                      'pars': pars,
+                      'values': {}}
+            
+            # Store fit results
+            phname = str(phmin)+'-'+str(phmax)
+            source = self._fitmodels[phname][self._srcname]
+            for par in pars:
+                result['values'][par]      = source[par].value()
+                result['values']['e_'+par] = source[par].error()
+            
+            # Append result to list of dictionaries
+            results.append(result)
+            
+        # Create FITS table from results
+        table = self._create_fits_table(results)
+
+        # Create FITS file and append FITS table to FITS file
+        self._fits = gammalib.GFits()
+        self._fits.append(table)
+        
+        # Return
+        return
+
+    def _save_fits(self):
+        """
+        Saved phase fit results into a fits file
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save phase FITS file')
+
+        # Create FITS file
+        self._create_fits()
+
+        # Get file name
+        outfile = self['outfile'].filename()
+
+        # Log file name
+        self._log_value(gammalib.NORMAL, 'Phase curve file', outfile.url())
+
+        # Save to fits
+        self._fits.saveto(outfile, self._clobber())
+        
+        # Return
+        return
+
+    def _save_xml(self):
+        """
+        Save phase fit results into one XML model each bin
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save phase XML files')
+
+        # Get file name
+        outname = self['outfile'].filename()
+
+        # Loop over all phase bins
+        for phbin in self._phbins:
+
+            # Set file name for phase bin
+            phname        = str(phbin[0])+'-'+str(phbin[1])
+            replace_name  = '_phase_'+phname+'.xml'
+            outname_phase = (outname.file()).replace('.fits', replace_name)
+            outfile       = outname.path() + outname_phase
+
+            # Save XML file
+            self._fitmodels[phname].save(outfile)
+
+            # Log file name
+            name = 'Phase [%s] file' % phname
+            self._log_value(gammalib.NORMAL, name, outfile)
+        
+        # Return
+        return
+
     # Public methods
     def run(self):
@@ -264,24 +427,26 @@ class csphasecrv(ctools.cscript):
         # Get parameters
         self._get_parameters()
-        
+
         # Write observation into logger
-        if self._logTerse():
-            self._log('\n')
-            self._log.header1(gammalib.number('Input observation',len(self._obs)))
-            self._log(str(self._obs))
-            self._log('\n')
-                
-        #Store and clear observations to split them in phase
+        self._log_observations(gammalib.NORMAL, self._obs, 'Observation')
+        
+        # Store and clear observations to split them into phase bins
         orig_obs = self._obs.copy()
-        #self._obs.clear()
-        #Dictionary to save phase fitted models
+        # Dictionary to save phase fitted models
         self._fitmodels = {}
-        
-        #Loop over all phases
+
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Generate phase curve')
+
+        # Loop over all phases
         for phbin in self._phbins:
-            
+
+            # Write time bin into header
+            if self._logTerse():
+                self._log.header2('PHASE '+str(phbin[0])+'-'+str(phbin[1]))
+
             # Select events
             select = ctools.ctselect(orig_obs)
             select['emin'] = self['emin'].real()
@@ -291,13 +456,13 @@ class csphasecrv(ctools.cscript):
             select['rad']  = 'UNDEFINED'
             select['ra']   = 'UNDEFINED'
             select['dec']  = 'UNDEFINED'
-            select['expr'] = "PHASE>"+str(phbin[0])+" && PHASE<"+str(phbin[1])
+            select['expr'] = 'PHASE>'+str(phbin[0])+' && PHASE<'+str(phbin[1])
             select.run()
-            #Add phase to observation id
+            # Add phase to observation id
             for i in range(0,select.obs().size()):
                 oldid = select.obs()[i].id()
-                select.obs()[i].id(oldid+"_"+str(phbin[0])+"-"+str(phbin[1]))
+                select.obs()[i].id(oldid+'_'+str(phbin[0])+'-'+str(phbin[1]))
             obs = select.obs()
             # If a stacked analysis is requested then bin the events
@@ -318,15 +483,18 @@ class csphasecrv(ctools.cscript):
             # Renormalize models to phase selection
             # TODO move the scaling from the temporal to the spectral component
             for model in like.obs().models():
-                scaled_norm = model["Normalization"].value()/(phbin[1] - phbin[0])
-                model["Normalization"].value(scaled_norm)
+                scaled_norm = model['Normalization'].value()/(phbin[1] - phbin[0])
+                model['Normalization'].value(scaled_norm)
-            #Store fit model
-            self._fitmodels[str(phbin[0])+"-"+str(phbin[1])] = like.obs().models().copy()
+            # Store fit model
+            self._fitmodels[str(phbin[0])+'-'+str(phbin[1])] = like.obs().models().copy()
+
+        # Create FITS file
+        self._create_fits()
-            ## Store observations
-            #for ob in select.obs():
-                #self._obs.append(ob.copy())
+        # Optionally publish light curve
+        if self['publish'].boolean():
+            self.publish()
         # Return
         return
@@ -350,171 +518,57 @@ class csphasecrv(ctools.cscript):
         # Return
         return
-    def saveXML(self):
-        """
-        Save phase fits into one XML model each
-        """
-        # Write header
-        self._log_header1(gammalib.TERSE, 'Save phase XML')
-
-        # Get out filename
-        outfile = str()
-
-        # Save XMLs
-        for phbin in self._phbins:
-            phname = str(phbin[0])+"-"+str(phbin[1])
-            outname = self['outfile'].filename()
-            outname_ph = (outname.file()).replace(".fits","_ph_"+phname+".xml")
-            outfile = outname.path()+outname_ph
-            print outfile
-            self._fitmodels[phname].save(outfile)
-        
-        # Return
-        return
-    
     def save(self):
         """
         Saves results to fits and XML
         """
-        print "SAVING"
-        self.savefits()
-        self.saveXML()
+        # Save results to FITS
+        self._save_fits()
+
+        # Save results to XML files (one per phase bin)
+        self._save_xml()
         # Return
         return
-    def phaserv(self):
-        """
-        Return dictionary with best fit models
-        """
-        # Return
-        return self._fitmodels
-
-    def _get_free_par_names(self):
+    def publish(self, name=''):
         """
-        Return list of free parameter names
+        Publish phase curve
-        Returns
-        -------
-        names : list of str
-            List of free parameter names.
-        """
-        # Initialise empty list of free parameter names
-        names = []
-
-        # Collect list of free parameter names
-        for par in self._obs.models()[self._srcname]:
-            if par.is_free():
-                names.append(par.name())
-
-        # Return names
-        return names
-
-    def savefits(self):
-        """
-        Saved phase fit results to a fits file
+        Parameters
+        ----------
+        name : str, optional
+            Name of phase curve
         """
-        # Initialise list of result dictionaries
-        results = []
-
-        # Get source parameters
-        pars = self._get_free_par_names()
-
-        # Loop over time bins
-        for i in range(len(self._phbins)):
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Publish phase curve')
-            # Get time boundaries
-            phmin = self._phbins[i][0]
-            phmax = self._phbins[i][1]
+        # Continue only if light curve is valid
+        if self._fits.contains('PHASECURVE'):
-            # Write time bin into header
-            if self._logTerse():
-                self._log.header2('PHASE '+str(phmin)+'-'+str(phmax))
+            # Set default name is user name is empty
+            if not name:
+                user_name = self._name
+            else:
+                user_name = name
-            # Initialise result dictionary
-            result = {'phmin': phmin,
-                      'phmax': phmax,
-                      'pars': pars,
-                      'values': {}}
-                        # Extract parameter values
-            
-            #Store fit results
-            phname = str(phmin)+"-"+str(phmax)
-            source = self._fitmodels[phname][self._srcname]
-            for par in pars:
-                result['values'][par]      = source[par].value()
-                result['values']['e_'+par] = source[par].error()
-            
-            # Append result to list of dictionaries
-            results.append(result)
-            
-        # Create FITS table from results
-        table = self._create_fits_table(results)
+            # Log file name
+            self._log_value(gammalib.NORMAL, 'Phase curve name', user_name)
-        # Create FITS file and append FITS table to FITS file
-        self._fits = gammalib.GFits()
-        self._fits.append(table)
+            # Publish phase curve
+            self._fits.publish('PHASECURVE', user_name)
-        # Save to fits
-        outfile = self['outfile'].filename()
-        print outfile
-        self._fits.saveto(outfile, self._clobber())
-        
-        #Return
+        # Return
         return
-    def _create_fits_table(self, results):
+    def phasecurve(self):
         """
-        Creates FITS binary table containing light curve results
-
-        Parameters
-        ----------
-        results : list of dict
-            List of result dictionaries
-
-        Returns
-        -------
-        table : `~gammalib.GFitsBinTable`
-            FITS binary table containing light curve
+        Return dictionary with best fit models
         """
-        # Determine number of rows in FITS table
-        nrows = len(results)
-
-        # Create FITS Table with extension "LIGHTCURVE"
-        table = gammalib.GFitsBinTable(nrows)
-        table.extname('PHASESPEC')
-
-        # Append time columns "MJD" and "e_MJD"
-        phmin   = gammalib.GFitsTableDoubleCol('PHMIN', nrows)
-        phmax = gammalib.GFitsTableDoubleCol('PHMAX', nrows)
-        #mjd.unit('days')
-        #e_mjd.unit('days')
-        for i, result in enumerate(results):
-            phmin[i]   = result['phmin']
-            phmax[i]   = result['phmax']
-        table.append(phmin)
-        table.append(phmax)
-
-        # Create parameter columns
-        for par in self._obs.models()[self._srcname]:
-            if par.is_free():
-                name   = par.name()
-                e_name = 'e_'+par.name()
-                col    = gammalib.GFitsTableDoubleCol(name, nrows)
-                e_col  = gammalib.GFitsTableDoubleCol(e_name, nrows)
-                col.unit(par.unit())
-                e_col.unit(par.unit())
-                for i, result in enumerate(results):
-                    col[i]   = result['values'][name]
-                    e_col[i] = result['values'][e_name]
-                table.append(col)
-                table.append(e_col)
-
-        # Return table
-        return table
+        # Return
+        return self._fitmodels
-        
 # ======================== #
 # Main routine entry point #
 # ======================== #
@@ -31,7 +31,7 @@ ntrials,  i, a, 10,,, &quot;Number of trials&quot;
 ra,       r, a, 83.6331,0,360, "RA of pointing (deg)"
 dec,      r, a, 22.0145,-90,90, "Dec of pointing (deg)"
 coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
-proj,     s, a, TAN,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+proj,     s, a, TAN,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 emin,     r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,     r, a, 100.0,,, "Upper energy limit (TeV)"
 enumbins, i, a, 0,,, "Number of energy bins (0=unbinned)"
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This script generates the pull distribution for all model parameters.
 #
-# 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 cspull(ctools.cscript):
         """
         # Set name
         self._name    = 'cspull'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise some members
         self._edisp       = False
@@ -37,7 +37,7 @@ emax,      r, a, 100.0,,, &quot;Upper energy limit (TeV)&quot;
 enumbins,  i, h, 20,,, "Number of model cube energy bins"
 ebinfile,  f, a, NONE,,, "Name of the file containing the energy bin definition"
 coordsys,  s, a, CEL,CEL|GAL,, "Coordinate System"
-proj,      s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+proj,      s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a residual map.
 #
-# Copyright (C) 2014-2016 Michael Mayer
+# Copyright (C) 2014-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -38,7 +38,7 @@ class csresmap(ctools.cscript):
         """
         # Set name
         self._name    = 'csresmap'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._resmap       = None
@@ -50,7 +50,7 @@ class csroot2caldb(ctools.cscript):
         """
         # Set name
         self._name    = 'csroot2caldb'
-        self._version = '1.3.0'
+        self._version = ctools.__version__
         # Initialise application by calling the appropriate class constructor
         self._init_cscript(argv)
@@ -50,7 +50,7 @@ class cssens(ctools.cscript):
         """
         # Set name
         self._name    = 'cssens'
-        self._version = '1.3.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._ebounds     = gammalib.GEbounds()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a spectrum.
 #
-# Copyright (C) 2014-2016 Michael Mayer
+# Copyright (C) 2014-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -43,7 +43,7 @@ class csspec(ctools.cscript):
         """
         # Set name
         self._name    = 'csspec'
-        self._version = '1.1.0'
+        self._version = ctools.__version__
         # Initialise some members
         self._ebounds     = gammalib.GEbounds()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Detects sources in a sky map
 #
-# 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
@@ -40,7 +40,7 @@ class cssrcdetect(ctools.cscript):
         """
         # Set name
         self._name    = 'cssrcdetect'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Set protected members
         self._models = gammalib.GModels()
@@ -36,7 +36,7 @@ ntrials,  i, a, 10,,, &quot;Number of trials&quot;
 ra,       r, a, 83.63,0,360, "RA of pointing (deg)"
 dec,      r, a, 22.01,-90,90, "Dec of pointing (deg)"
 coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)"
-proj,     s, a, TAN,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+proj,     s, a, TAN,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 emin,     r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,     r, a, 100.0,,, "Upper energy limit (TeV)"
 enumbins, i, a, 0,,, "Number of energy bins (0=unbinned)"
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates the TS distribution for a particular model
 #
-# 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
@@ -41,7 +41,7 @@ class cstsdist(ctools.cscript):
         """
         # Set name
         self._name    = 'cstsdist'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Initialise some members
         self._srcname     = ''
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Merge Test Statistic maps
 #
-# Copyright (C) 2015-2016 Michael Mayer
+# Copyright (C) 2015-2017 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -40,7 +40,7 @@ class cstsmapmerge(ctools.cscript):
         """
         # Set name
         self._name    = 'cstsmapmerge'
-        self._version = '1.1.0'
+        self._version = ctools.__version__
         # Initialise class members
         self._files        = None
@@ -41,9 +41,9 @@ nxpix,     i, a, 200,,, &quot;Size of the X axis in pixels&quot;
 nypix,     i, a, 200,,, "Size of the Y axis in pixels"
 binsz,     r, a, 0.02,,, "Image scale (in degrees/pixel)"
 coordsys,  s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,      s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,      s, a, CAR, AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Standard parameters
@@ -47,7 +47,7 @@ class cstsmapsplit(ctools.cscript):
         # Set name and version
         self._name    = 'cstsmapsplit'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Set data members
         self._outmap       = gammalib.GFilename()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Compute a visibility cube
 #
-# 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
@@ -50,7 +50,7 @@ class csviscube(ctools.cscript):
         """
         # Set name
         self._name    = 'csviscube'
-        self._version = '1.0.0'
+        self._version = ctools.__version__
         # Initialise members
         self._cube = gammalib.GSkyMap()
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Executes analysis workflow
 #
-# 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
@@ -53,7 +53,7 @@ class csworkflow(ctools.cscript):
         """
         # Set name and version
         self._name    = 'csworkflow'
-        self._version = '1.2.0'
+        self._version = ctools.__version__
         # Set members
         self._workflow = gammalib.GXml()
@@ -75,7 +75,7 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
 ``xref [real]``
@@ -9,13 +9,16 @@ Computes phase dependent spectra for a given source.
 Synopsis
 --------
-This script computes spectra by performing a maximum likelihood fit
-using :doc:`ctlike` in a series of phase bins for pulsars.
-The phase bins can be either specified in an ASCII file, as an interval
-divided into equally sized time bins.
-The format of the ASCII file is one row per time bin, each specifying the
-start of stop value of the phase bin, separated by a whitespace. The phase 
-goes from 0.0 to 1.0 .
+This script computes spectra by performing a maximum likelihood fit using
+:doc:`ctlike` in a series of phase bins for pulsars. The phase bins can be
+either specified in an ASCII file or as an interval divided into equally sized
+phase bins. The format of the ASCII file is one row per phase bin, each
+specifying the start of stop value of the phase bin, separated by a whitespace.
+The phase goes from 0.0 to 1.0.
+
+On output the script writes the fitting results into a FITS file. The script
+also produces one XML file per phase bin that contains the fitting results for
+that bin.
 General parameters
@@ -25,10 +28,10 @@ General parameters
     Input event list or observation definition XML file.
 ``inmodel [file]``
-    Input model XML file.
+    Input model definition XML file.
 ``srcname [string]``
-    Name of the periodic source in the source model XML file
+    Name of the periodic source in the source model XML file.
 ``caldb [string]``
     Calibration database.
@@ -40,10 +43,11 @@ General parameters
     Applies energy dispersion to response computation.
 ``outfile [file]``
-    Name of the XML output file. The phase intervall will be automatically appended to the name.
+    Name of the XML output file. The phase interval will be automatically
+    appended to the name.
-``(tbinalg = LIN) <FILE|LIN> [string]``
-    Algorithm for defining time bins.
+``phbinalg <FILE|LIN> [string]``
+    Algorithm for defining phase bins.
 ``phbins [integer]``
     Number of phase bins.
@@ -63,7 +67,7 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
 ``xref [real]``
@@ -81,9 +85,13 @@ General parameters
 ``binsz [real]``
     Pixel size (in degrees/pixel).
+
 Standard parameters
 -------------------
+``(publish = no) [boolean]``
+    Specifies whether the phase curve should be published on VO Hub.
+
 ``(chatter = 2) [integer]``
     Verbosity of the executable:
      chatter = 0: no information will be logged
@@ -71,7 +71,7 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
 ``emin [real]``
@@ -76,7 +76,7 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
 ``xref [real]``
@@ -78,7 +78,7 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
 ``emin [real]``
@@ -80,16 +80,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``xref [real]``
     Right Ascension / Galactic longitude of map centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of map centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
-
-
 Standard parameters
 -------------------
@@ -65,15 +65,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+ 	 	 
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of cube centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
- 	 	 
 Standard parameters
 -------------------
@@ -105,15 +105,15 @@ General parameters
 ``coordsysL <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+ 	 	 
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of cube centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
- 	 	 
 Standard parameters
 -------------------
@@ -88,15 +88,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of cube centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
-
 ``(migramax = 2.0) [real]``
     Upper bound of ratio between reconstructed and true photon energy.
@@ -95,15 +95,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
 ``yref [real]``
-    Declination / Galactic latitude of cube centre (J2000, in degrees).
- 	 	 
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
- 	 	 
+    Declination / Galactic latitude of cube centre (J2000, in degrees). 	 	 
+
 Standard parameters
 -------------------
@@ -61,15 +61,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+ 	 	 
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of cube centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
- 	 	 
 Standard parameters
 -------------------
@@ -125,15 +125,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``xref [real]``
     Right Ascension / Galactic longitude of image centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of image centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
-
 Standard parameters
 -------------------
@@ -87,15 +87,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``xref [real]``
     Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of cube centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
-
 ``(amax = 0.3) [real]``
     Upper bound of angular separation between true and measued photon
     direction (in degrees).
@@ -65,7 +65,7 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
-``proj <AIT|AZP|CAR|MER|STG|TAN> [string]``
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
     Projection method.
 ``xref [real]``
@@ -78,15 +78,15 @@ General parameters
 ``coordsys <CEL|GAL> [string]``
     Coordinate system (CEL - celestial, GAL - galactic).
+``proj <AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN> [string]``
+    Projection method.
+
 ``xref [real]``
     Right Ascension / Galactic longitude of map centre (J2000, in degrees).
 ``yref [real]``
     Declination / Galactic latitude of map centre (J2000, in degrees).
-``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
-    Projection method.
-
 ``(binmin = -1) [integer]``
     First bin to compute.
@@ -43,13 +43,14 @@ cscripts
    cscaldb --- Lists available instrument response functions <cscaldb>
    csebins --- Generates energy boundaries for stacked analysis <csebins>
-   cslightcrv --- Computes lightcurve <cslightcrv>
+   cslightcrv --- Computes light curve <cslightcrv>
    csmodelinfo --- Shows model container content <csmodelinfo>
    csmodelmerge --- Merges several model containers into one file <csmodelmerge>
    csmodelselect --- Select models from model definition file <csmodelselect>
    csobsdef --- Generates observation definition file <csobsdef>
    csobsinfo --- Shows observation container content <csobsinfo>
    csobsselect --- Select observations from observation definition file <csobsselect>
+   csphasecrv --- Computes phase curve <csphasecrv>
    cspull --- Generates pull distribution <cspull>
    csresmap --- Generates residual map <csresmap>
    cssens --- Computes CTA sensitivity <cssens>
@@ -5,7 +5,6 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctbin.par                                                     #
-#  Version:   1.1.0                                                         #
 #                                                                           #
 #  Author:    Juergen Knoedlseder                                           #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -37,9 +36,9 @@ nxpix,    i, a, 200,,, &quot;Size of the X axis in pixels&quot;
 nypix,    i, a, 200,,, "Size of the Y axis in pixels"
 binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
 coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Standard parameters
@@ -41,10 +41,10 @@ usepnt,   b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 nxpix,    i, a, 200,,, "Size of the X axis in pixels"
 nypix,    i, a, 200,,, "Size of the Y axis in pixels"
 binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
-coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Standard parameters
@@ -5,7 +5,6 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctedispcube.par                                               #
-#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Maria Haupt                                                   #
 #             maria.haupt@desy.de                                           #
@@ -40,10 +39,10 @@ usepnt,   b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 nxpix,    i, a, 10,,, "Size of the X axis in pixels"
 nypix,    i, a, 10,,, "Size of the Y axis in pixels"
 binsz,    r, a, 1.0,,, "Image scale (in degrees/pixel)"
-coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Angular separation binning parameters
@@ -5,7 +5,6 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctexpcube.par                                                 #
-#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Juergen Knoedlseder                                           #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -40,10 +39,10 @@ usepnt,   b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 nxpix,    i, a, 200,,, "Size of the X axis in pixels"
 nypix,    i, a, 200,,, "Size of the Y axis in pixels"
 binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
-coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Standard parameters
@@ -5,7 +5,6 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctmapcube.par                                                 #
-#  Version:   1.0.0                                                         #
 #                                                                           #
 #  Author:    Juergen Knoedlseder                                           #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -37,9 +36,9 @@ nxpix,    i, a, 200,,, &quot;Size of the X axis in pixels&quot;
 nypix,    i, a, 200,,, "Size of the Y axis in pixels"
 binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
 coordsys, s, a, CEL,CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Standard parameters
@@ -54,10 +54,10 @@ usepnt,   b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 nxpix,    i, a, 200,,, "Size of the X axis in pixels"
 nypix,    i, a, 200,,, "Size of the Y axis in pixels"
 binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
-coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Standard parameters
@@ -5,7 +5,6 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctpsfcube.par                                                 #
-#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Chia-Chun Lu                                                  #
 #             c.lu@mpi-hd.mpg.de                                            #
@@ -40,10 +39,10 @@ usepnt,   b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 nxpix,    i, a, 10,,, "Size of the X axis in pixels"
 nypix,    i, a, 10,,, "Size of the Y axis in pixels"
 binsz,    r, a, 1.0,,, "Image scale (in degrees/pixel)"
-coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,     s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,     s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Angular separation binning parameters
@@ -5,7 +5,6 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctselect.par                                                  #
-#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Jurgen Knodlseder                                             #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -32,8 +31,8 @@ tmin,    r, a, 0.0,,, &quot;Start time (CTA MET in seconds)&quot;
 tmax,    r, a, 0.0,,, "End time (CTA MET in seconds)"
 emin,    r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
-expr,    s, h, "",,, "Additional selection expression"
-usethres,s, h, "NONE",NONE|DEFAULT|USER,, "Energy threshold type"
+expr,    s, h, ,,, "Additional selection expression"
+usethres,s, h, NONE,NONE|DEFAULT|USER,, "Energy threshold type"
 #
 # Standard parameters
@@ -35,7 +35,7 @@ nxpix,    i, a, 200,,, &quot;Size of the X axis in pixels&quot;
 nypix,    i, a, 200,,, "Size of the Y axis in pixels"
 binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
 coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
-proj,     s, a, CAR, AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+proj,     s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,     r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
@@ -40,10 +40,10 @@ usepnt,    b, h, no,,, &quot;Use pointing instead of xref/yref parameters?&quot;
 nxpix,     i, a, 200,,, "Size of the X axis in pixels"
 nypix,     i, a, 200,,, "Size of the Y axis in pixels"
 binsz,     r, a, 0.02,,, "Image scale (in degrees/pixel)"
-coordsys,  s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+coordsys,  s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+proj,      s, a, CAR, AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
 xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
 yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
-proj,      s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
 #
 # Parameters for splitting and speed purpose
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the csphasecrv script.
+#
+# Copyright (C) 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
+# 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 gammalib
+import cscripts
+from testing import test
+
+
+# ================================ #
+# Test class for csphasecrv script #
+# ================================ #
+class Test(test):
+    """
+    Test class for csphasecrv script
+
+    This test class makes unit tests for the csphasecrv script by using it
+    from the command line and from Python.
+    """
+    
+    # Constructor
+    def __init__(self):
+        """
+        Constructor
+        """
+        # Call base class constructor
+        test.__init__(self)
+
+        # Set test data
+        self._phased_events  = self._datadir + '/phased_events.fits'
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name('csphasecrv')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test csphasecrv on command line')
+        self.append(self._test_python, 'Test csphasecrv from Python')
+
+        # Return
+        return
+
+    # Test csphasecrv on command line
+    def _test_cmd(self):
+        """
+        Test csphasecrv on the command line.
+        """
+        # Set script name
+        csphasecrv = self._script('csphasecrv')
+
+        # Setup csphasecrv command
+        cmd = csphasecrv+' inobs="'+self._phased_events+'"'+ \
+                         ' inmodel="'+self._model+'" srcname="Crab"'+ \
+                         ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+                         ' phbinalg="LIN" phbins=5'+ \
+                         ' enumbins=0 emin=0.1 emax=100.0'+ \
+                         ' outfile="csphasecrv_cmd1.fits"'+ \
+                         ' logfile="csphasecrv_cmd1.log" chatter=1'
+
+        # Check if execution of wrong command fails
+        self.test_assert(self._execute('command_that_does_not_exist') != 0,
+             'Self test of test script')
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Check phase curve
+        self._check_phase_curve('csphasecrv_cmd1.fits', 5)
+
+        # Setup csphasecrv command
+        cmd = csphasecrv+' inobs="events_that_do_not_exist.fits"'+ \
+                         ' inmodel="'+self._model+'" srcname="Crab"'+ \
+                         ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+                         ' phbinalg="LIN" phbins=5'+ \
+                         ' enumbins=0 emin=0.1 emax=100.0'+ \
+                         ' outfile="csphasecrv_cmd2.fits"'+ \
+                         ' logfile="csphasecrv_cmd2.log" chatter=1'
+
+        # Check if execution failed
+        self.test_assert(self._execute(cmd) != 0,
+             'Check invalid input file when executed from command line')
+
+        # Return
+        return
+
+    # Test csphasecrv from Python
+    def _test_python(self):
+        """
+        Test csphasecrv from Python
+        """
+        # Set-up unbinned csphasecrv
+        pcrv = cscripts.csphasecrv()
+        pcrv['inobs']    = self._phased_events
+        pcrv['inmodel']  = self._model
+        pcrv['srcname']  = 'Crab'
+        pcrv['caldb']    = self._caldb
+        pcrv['irf']      = self._irf
+        pcrv['phbinalg'] = 'LIN'
+        pcrv['phbins']   = 5
+        pcrv['enumbins'] = 0
+        pcrv['emin']     = 0.1
+        pcrv['emax']     = 100.0
+        pcrv['outfile']  = 'csphasecrv_py1.fits'
+        pcrv['logfile']  = 'csphasecrv_py1.log'
+        pcrv['chatter']  = 2
+        pcrv['publish']  = True
+
+        # Run csphasecrv script and save phase curve
+        pcrv.logFileOpen()   # Make sure we get a log file
+        pcrv.run()
+        pcrv.save()
+
+        # Check phase curve
+        self._check_phase_curve('csphasecrv_py1.fits', 5)
+
+        # Now use FILE as time bin algorithm. For this we need first to
+        # create an ASCII file. We use now 2 phase bins. The ASCII file
+        # is saved into the file "csphasecrv_py2.dat".
+        csv     = gammalib.GCsv(2,2)
+        phmin   = 0.0
+        phdelta = 1.0 / float(csv.nrows())
+        for i in range(csv.nrows()):
+            csv[i,0] = '%.5f' % (phmin +  i    * phdelta)
+            csv[i,1] = '%.5f' % (phmin + (i+1) * phdelta)
+        csv.save('csphasecrv_py2.dat', ' ', True)
+
+        # Set-up unbinned csphasecrv
+        pcrv = cscripts.csphasecrv()
+        pcrv['inobs']     = self._phased_events
+        pcrv['inmodel']   = self._model
+        pcrv['srcname']   = 'Crab'
+        pcrv['caldb']     = self._caldb
+        pcrv['irf']       = self._irf
+        pcrv['phbinalg']  = 'FILE'
+        pcrv['phbinfile'] = 'csphasecrv_py2.dat'
+        pcrv['enumbins']  = 0
+        pcrv['emin']      = 0.1
+        pcrv['emax']      = 100.0
+        pcrv['outfile']   = 'csphasecrv_py2.fits'
+        pcrv['logfile']   = 'csphasecrv_py2.log'
+        pcrv['chatter']   = 3
+
+        # Execute csphasecrv script
+        pcrv.logFileOpen()   # Make sure we get a log file
+        pcrv.execute()
+
+        # Check phase curve
+        self._check_phase_curve('csphasecrv_py2.fits', csv.nrows())
+
+        # Now we setup an observation container on input. We attached the
+        # model to the observation container so that csphasecrv should
+        # no longer query for the parameter.
+        cta = gammalib.GCTAObservation(self._phased_events)
+        obs = gammalib.GObservations()
+        obs.append(cta)
+        obs.models(self._model)
+
+        # Set-up unbinned csphasecrv from observation container. Now use
+        # the GTI algorithm so that we test all timing algorithms.
+        pcrv = cscripts.csphasecrv(obs)
+        pcrv['srcname']  = 'Crab'
+        pcrv['caldb']    = self._caldb
+        pcrv['irf']      = self._irf
+        pcrv['phbinalg'] = 'LIN'
+        pcrv['phbins']   = 5
+        pcrv['enumbins'] = 0
+        pcrv['emin']     = 0.1
+        pcrv['emax']     = 100.0
+        pcrv['outfile']  = 'csphasecrv_py3.fits'
+        pcrv['logfile']  = 'csphasecrv_py3.log'
+        pcrv['chatter']  = 4
+
+        # Execute csphasecrv script
+        pcrv.logFileOpen()   # Make sure we get a log file
+        pcrv.execute()
+
+        # Check phase curve
+        self._check_phase_curve('csphasecrv_py3.fits', 5)
+
+        # Finally we set-up a binned csphasecrv
+        pcrv = cscripts.csphasecrv()
+        pcrv['inobs']    = self._phased_events
+        pcrv['inmodel']  = self._model
+        pcrv['srcname']  = 'Crab'
+        pcrv['caldb']    = self._caldb
+        pcrv['irf']      = self._irf
+        pcrv['phbinalg'] = 'LIN'
+        pcrv['phbins']   = 2
+        pcrv['emin']     = 0.1
+        pcrv['emax']     = 100.0
+        pcrv['enumbins'] = 10
+        pcrv['coordsys'] = 'CEL'
+        pcrv['proj']     = 'TAN'
+        pcrv['xref']     = 83.63
+        pcrv['yref']     = 22.01
+        pcrv['nxpix']    = 20
+        pcrv['nypix']    = 20
+        pcrv['binsz']    = 0.02
+        pcrv['outfile']  = 'csphasecrv_py4.fits'
+        pcrv['logfile']  = 'csphasecrv_py4.log'
+        pcrv['chatter']  = 4
+
+        # Execute csphasecrv script
+        pcrv.logFileOpen()   # Make sure we get a log file
+        pcrv.execute()
+
+        # Check phase curve
+        self._check_phase_curve('csphasecrv_py4.fits', 2)
+
+        # Return
+        return
+
+    # Check phase curve result file
+    def _check_phase_curve(self, filename, bins):
+        """
+        Check phase curve file
+        """
+        # Expected column names
+        cols = ['PHASE_MIN', 'PHASE_MAX', 'Prefactor', 'e_Prefactor',
+                'Index', 'e_Index']
+
+        # Open FITS file
+        fits = gammalib.GFits(filename)
+
+        # Check FITS file structure
+        self.test_value(fits.size(), 2,
+             'Check for 2 extensions in phase curve FITS file')
+        self.test_assert(fits.contains('PHASECURVE'),
+             'FITS file contains "PHASECURVE" extension')
+
+        # Get PHASECURVE table
+        table = fits['PHASECURVE']
+
+        # Check FITS table structure
+        self.test_value(table.ncols(), len(cols),
+             'Check for %d columns in phase curve FITS table' % len(cols))
+        self.test_value(table.nrows(), bins,
+             'Check for %d rows in phase curve FITS table' % bins)
+        for col in cols:
+            self.test_assert(table.contains(col),
+                 'FITS file contains "'+col+'" column')
+
+        # Close FITS file
+        fits.close()
+        
+        # Return
+        return
@@ -36,6 +36,7 @@ import test_csobs2caldb
 import test_csobsdef
 import test_csobsinfo
 import test_csobsselect
+import test_csphasecrv
 import test_cspull
 import test_csresmap
 import test_csroot2caldb
@@ -124,6 +125,7 @@ def test(installed=False, debug=False):
              test_csobsdef.Test(),
              test_csobsinfo.Test(),
              test_csobsselect.Test(),
+             test_csphasecrv.Test(),
              test_cspull.Test(),
              test_csresmap.Test(),
              test_csroot2caldb.Test(),