Replace ASCII output files by FITS files in cssens (#1707)

Jürgen Knödlseder
1 parent b58507bd
Showing 11 changed files with 271 additions and 72 deletions
ChangeLog
NEWS
README.md
cscripts/csphasecrv.py
cscripts/cssens.par
cscripts/cssens.py
doc/source/admin/release_history/2.0.rst
doc/source/users/reference_manual/cssens.rst
examples/show_sensitivity.py
test/test_csphasecrv.py
test/test_cssens.py
-2021-05-07
+2021-05-11
         * Version 2.0.0 released
         ========================
@@ -17,6 +17,7 @@
         Replace std::tmpnam() by gammalib::tmpnam() function in ctselect (#3324)
         Remove implicit conversion warnings (#3324)
         Add support of rectangular On/Off regions (#3248)
+        Replace ASCII output files by FITS files (#1707)
 2021-04-28
 New Features and Important Changes in ctools 2.0.0
-7 May 2021
+11 May 2021
 Introduction
@@ -286,6 +286,8 @@ The test source position is now used as pointing direction, eventually offset
 in galactic latitude by a number of degrees specified using the "offset"
 parameter (#3555).
+The ASCII output file was replaced by a FITS file (#1707).
+
 csspec - Generates Spectral Energy Distribution
 -----------------------------------------------
@@ -358,7 +360,7 @@ regions (#3248).
 Examples
 --------
-None
+The show_sensitivity.py was adapted to deal also with FITS output files (#1707).
 caldb
 ctools information
 ==================
-* Version:             2.0.0.dev (7 May 2021)
+* Version:             2.0.0.dev (11 May 2021)
 * GammaLib dependency: 2.0.0.dev
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates spectra as a function of phase
 #
-# Copyright (C) 2017-2019 Rolf Buehler
+# Copyright (C) 2017-2021 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
@@ -88,12 +88,12 @@ class csphasecrv(ctools.csobservation):
         state = {'base'         : ctools.csobservation.__getstate__(self),
                  'srcname'      : self._srcname,
                  'phbins'       : self._phbins,
-                 'stacked'      : self._stacked,
                  'onoff'        : self._onoff,
+                 'stacked'      : self._stacked,
                  'fits'         : self._fits,
+                 'fitmodels'    : self._fitmodels,
                  'nthreads'     : self._nthreads,
-                 'excl_reg_map' : self._excl_reg_map,
-                 'fitmodels'    : self._fitmodels}
+                 'excl_reg_map' : self._excl_reg_map}
         # Return pickled dictionary
         return state
@@ -110,12 +110,12 @@ class csphasecrv(ctools.csobservation):
         ctools.csobservation.__setstate__(self, state['base'])
         self._srcname      = state['srcname']
         self._phbins       = state['phbins']
-        self._stacked      = state['stacked']
         self._onoff        = state['onoff']
+        self._stacked      = state['stacked']
         self._fits         = state['fits']
+        self._fitmodels    = state['fitmodels']
         self._nthreads     = state['nthreads']
         self._excl_reg_map = state['excl_reg_map']
-        self._fitmodels    = state['fitmodels']
         # Return
         return
@@ -291,7 +291,7 @@ class csphasecrv(ctools.csobservation):
                       'phmax': phmax,
                       'pars': pars,
                       'values': {}}
-            
+
             # Store fit results
             phname = str(phmin)+'-'+str(phmax)
@@ -310,14 +310,14 @@ class csphasecrv(ctools.csobservation):
             # 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
@@ -339,7 +339,7 @@ class csphasecrv(ctools.csobservation):
         # Save to fits
         self._fits.saveto(outfile, self._clobber())
-        
+
         # Return
         return
@@ -368,25 +368,33 @@ class csphasecrv(ctools.csobservation):
             # Log file name
             name = 'Phase [%s] file' % phname
             self._log_value(gammalib.NORMAL, name, outfile)
-        
+
         # Return
         return
-    def _phase_bin(self,phbin):
+    def _phase_bin(self, phbin):
+        """
+        Analyse one phase bin
+
+        Parameters
+        ----------
+        phbin : list of float
+            Phase bin limits
+        """
         # Write time bin into header
         self._log_header2(gammalib.TERSE, 'PHASE %f - %f' %
                           (phbin[0], phbin[1]))
         # Select events
-        select = ctools.ctselect(self.obs().copy())
+        select = ctools.ctselect(self.obs())
         select['emin'] = self['emin'].real()
         select['emax'] = self['emax'].real()
         select['tmin'] = 'UNDEFINED'
         select['tmax'] = 'UNDEFINED'
-        select['rad'] = 'UNDEFINED'
-        select['ra'] = 'UNDEFINED'
-        select['dec'] = 'UNDEFINED'
+        select['rad']  = 'UNDEFINED'
+        select['ra']   = 'UNDEFINED'
+        select['dec']  = 'UNDEFINED'
         select['expr'] = 'PHASE>' + str(phbin[0]) + ' && PHASE<' + str(phbin[1])
         select.run()
@@ -457,13 +465,13 @@ class csphasecrv(ctools.csobservation):
         # Switch screen logging on in debug mode
         if self._logDebug():
             self._log.cout(True)
-        
+
         # Get parameters
         self._get_parameters()
         # Write observation into logger
         self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
-        
+
         # Dictionary to save phase fitted models
         self._fitmodels = {}
@@ -507,7 +515,7 @@ class csphasecrv(ctools.csobservation):
         # Save results to XML files (one per phase bin)
         self._save_xml()
-        
+
         # Return
         return
@@ -569,7 +577,7 @@ class csphasecrv(ctools.csobservation):
         # Return
         return self._excl_reg_map
-        
+
 # ======================== #
 # Main routine entry point #
 # ======================== #
@@ -23,7 +23,7 @@ caldb,   s, a, prod2,,, &quot;Calibration database&quot;
 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"
+outfile, f, a, sensitivity.fits,,, "Output sensitivity FITS file"
 #
 # Script parameters
@@ -51,5 +51,3 @@ clobber,  b, h, yes,,, &quot;Overwrite existing output files with new output files?&quot;
 debug,    b, h, no,,, "Activate debugging mode?"
 mode,     s, h, ql,,, "Mode of automatic parameters"
 logfile,  f, h, cssens.log,,, "Log filename"
-
-
@@ -54,6 +54,7 @@ class cssens(ctools.csobservation):
         # Initialise class members
         self._ebounds     = gammalib.GEbounds()
         self._obs_ebounds = []
+        self._fits        = None
         self._srcname     = ''
         self._ra          = None
         self._dec         = None
@@ -78,6 +79,7 @@ class cssens(ctools.csobservation):
         state = {'base'         : ctools.csobservation.__getstate__(self),
                  'ebounds'      : self._ebounds,
                  'obs_ebounds'  : self._obs_ebounds,
+                 'fits'         : self._fits,
                  'srcname'      : self._srcname,
                  'ra'           : self._ra,
                  'dec'          : self._dec,
@@ -100,6 +102,7 @@ class cssens(ctools.csobservation):
         ctools.csobservation.__setstate__(self, state['base'])
         self._ebounds     = state['ebounds']
         self._obs_ebounds = state['obs_ebounds']
+        self._fits        = state['fits']
         self._srcname     = state['srcname']
         self._ra          = state['ra']
         self._dec         = state['dec']
@@ -145,7 +148,6 @@ class cssens(ctools.csobservation):
         self['sigma'].real()
         self['max_iter'].integer()
         self['type'].string()
-        self['outfile'].filename()
         self['edisp'].boolean()
         self['debug'].boolean()
@@ -154,6 +156,10 @@ class cssens(ctools.csobservation):
                                           gammalib.GEnergy(emin, 'TeV'),
                                           gammalib.GEnergy(emax, 'TeV'))
+        # Read ahead output parameters
+        if self._read_ahead():
+            self['outfile'].filename()
+
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
@@ -690,6 +696,88 @@ class cssens(ctools.csobservation):
         # Return results
         return result
+    def _create_fits(self, results):
+        """
+        Create FITS file from results
+
+        Parameters
+        ----------
+        results : list of dict
+            List of result dictionaries
+        """
+        # Create FITS table columns
+        nrows        = len(results)
+        e_mean       = gammalib.GFitsTableDoubleCol('E_MEAN', nrows)
+        e_min        = gammalib.GFitsTableDoubleCol('E_MIN', nrows)
+        e_max        = gammalib.GFitsTableDoubleCol('E_MAX', nrows)
+        flux_crab    = gammalib.GFitsTableDoubleCol('FLUX_CRAB', nrows)
+        flux_photon  = gammalib.GFitsTableDoubleCol('FLUX_PHOTON', nrows)
+        flux_energy  = gammalib.GFitsTableDoubleCol('FLUX_ENERGY', nrows)
+        sensitivity  = gammalib.GFitsTableDoubleCol('SENSITIVITY', nrows)
+        regcoeff     = gammalib.GFitsTableDoubleCol('REGRESSION_COEFF', nrows)
+        nevents      = gammalib.GFitsTableDoubleCol('NEVENTS', nrows)
+        npred        = gammalib.GFitsTableDoubleCol('NPRED', nrows)
+        e_mean.unit('TeV')
+        e_min.unit('TeV')
+        e_max.unit('TeV')
+        flux_crab.unit('')
+        flux_photon.unit('ph/cm2/s')
+        flux_energy.unit('erg/cm2/s')
+        sensitivity.unit('erg/cm2/s')
+        regcoeff.unit('')
+        nevents.unit('counts')
+        npred.unit('counts')
+
+        # Fill FITS table columns
+        for i, result in enumerate(results):
+            e_mean[i]      = math.pow(10.0, result['loge'])
+            e_min[i]       = result['emin']
+            e_max[i]       = result['emax']
+            flux_crab[i]   = result['crab_flux']
+            flux_photon[i] = result['photon_flux']
+            flux_energy[i] = result['energy_flux']
+            sensitivity[i] = result['sensitivity']
+            regcoeff[i]    = result['regcoeff']
+            nevents[i]     = result['nevents']
+            npred[i]       = result['npred']
+
+        # Initialise FITS Table with extension "SENSITIVITY"
+        table = gammalib.GFitsBinTable(nrows)
+        table.extname('SENSITIVITY')
+
+        # Add keywors for compatibility with gammalib.GMWLSpectrum
+        table.card('INSTRUME', 'CTA', 'Name of Instrument')
+        table.card('TELESCOP', 'CTA', 'Name of Telescope')
+
+        # Stamp header
+        self._stamp(table)
+
+        # Add script keywords
+        table.card('OBJECT',   self._srcname, 'Source for which sensitivity was estimated')
+        table.card('TYPE',     self['type'].string(), 'Sensitivity type')
+        table.card('SIGMA',    self['sigma'].real(), '[sigma] Sensitivity threshold')
+        table.card('MAX_ITER', self['max_iter'].integer(), 'Maximum number of iterations')
+        table.card('STAT',     self['statistic'].string(), 'Optimization statistic')
+
+        # Append filled columns to fits table
+        table.append(e_mean)
+        table.append(e_min)
+        table.append(e_max)
+        table.append(flux_crab)
+        table.append(flux_photon)
+        table.append(flux_energy)
+        table.append(sensitivity)
+        table.append(regcoeff)
+        table.append(nevents)
+        table.append(npred)
+
+        # Create the FITS file now
+        self._fits = gammalib.GFits()
+        self._fits.append(table)
+
+        # Return
+        return
+
     # Public methods
     def run(self):
@@ -708,10 +796,7 @@ class cssens(ctools.csobservation):
         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', 'regcoeff', 'nevents',
-                    'npred']
+        # Initialise results
         results  = []
         # Set test source model for this observation
@@ -745,20 +830,50 @@ class cssens(ctools.csobservation):
         else:
             for ieng in range(self._ebounds.size()):
-                #Run analysis in energy bin
+                # Run analysis in energy bin
                 result = self._e_bin(ieng)
                 # Append results
                 results.append(result)
-        # Write out trial result
-        for ieng, result in enumerate(results):
-            ioutils.write_csv_row(self['outfile'].filename().url(), ieng,
-                                  colnames, result)
+        # Create FITS file
+        self._create_fits(results)
         # Return
         return
+    def save(self):
+        """
+        Save sensitivity FITS file
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save sensitivity curve')
+
+        # Continue only if FITS file is valid
+        if self._fits != None:
+
+            # Get outmap parameter
+            outfile = self['outfile'].filename()
+
+            # Log file name
+            self._log_value(gammalib.NORMAL, 'Sensitivity file', outfile.url())
+
+            # Save sensitivity
+            self._fits.saveto(outfile, self['clobber'].boolean())
+
+        # Return
+        return
+
+    def sensitivity(self):
+        """
+        Return sensitivity FITS file
+
+        Returns:
+            FITS file containing sensitivity curve
+        """
+        # Return
+        return self._fits
+
 # ======================== #
 # Main routine entry point #
 .. _2.0:
-ctools 2.0 release (7 May 2021)
-===============================
+ctools 2.0 release (11 May 2021)
+================================
 Overview
 --------
@@ -44,6 +44,8 @@ Improvements
   Replace ``std::tmpnam()`` by ``gammalib::tmpnam()`` function in :ref:`ctselect`
 * [`3324 <https://cta-redmine.irap.omp.eu/issues/3324>`_] -
   Remove implicit conversion warnings
+* [`1707 <https://cta-redmine.irap.omp.eu/issues/1707>`_] -
+  Replace ASCII output files by FITS files
 New features
@@ -42,17 +42,11 @@ bin/threshold will be performed in parallel over as many processes as the number
 CPUs available on your machine. The maximum number of parallel processes can be set
 by the user through the ``nthreads`` hidden parameter.
-:ref:`cssens` will generate an ASCII file in comma-separated value (CSV) format
-containing the sensitivity as function of energy. The first row is a header row
-providing the column names. The following rows provide the mean logarithmic
-energy and the boundaries of the energy bin for which the sensitivity was
-computed. They also provide the flux within the energy bin in Crab units, in
-photons (ph/cm2/s) and in energy (erg/cm2/s). Finally, the sensitivity is given
-as the test source spectrum evaluated at the mean logarithmic energy multiplied
-by the energy squared (erg/cm2/s).
+:ref:`cssens` will generate a FITS file with a single extension that contains a
+binary table with the sensitivity as function of energy.
-The sensitivity can be displayed using the ``show_sensitivity.py`` script in the
-example folder. Matplotlib is required to execute the script.
+The sensitivity FITS file can be displayed using the ``show_sensitivity.py`` script
+in the example folder. Matplotlib is required to execute the script.
 General parameters
@@ -81,7 +75,7 @@ General parameters
     Average deadtime correction factor.
 ``outfile [file]``
-    ASCII file containing the sensitivity values.
+    Output sensitivity FITS file.
 ``(offset = 0.0) [real]``
     Offset angle of source in field of view (in degrees).	 
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Display differential sensitivity plot generated by cssens
 #
-# Copyright (C) 2015-2020 Juergen Knoedlseder
+# Copyright (C) 2015-2021 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
@@ -36,6 +36,29 @@ import cscripts
 # ===================== #
 def read_sensitivity(filename):
     """
+    Read sensitivity information
+
+    Parameters
+    ----------
+    filename : str
+        Name of FITS or CSV file
+    """
+    # Read sensitivity information dependent on file type
+    fname = gammalib.GFilename(filename)
+    if fname.is_fits():
+        sensitivity = read_sensitivity_fits(filename)
+    else:
+        sensitivity = read_sensitivity_csv(filename)
+
+    # Return sensitivity
+    return sensitivity
+
+
+# =================================== #
+# Read sensitivity data from CSV file #
+# =================================== #
+def read_sensitivity_csv(filename):
+    """
     Read sensitivity information from CSV file
     Parameters
@@ -94,6 +117,61 @@ def read_sensitivity(filename):
     return sensitivity
+# ==================================== #
+# Read sensitivity data from FITS file #
+# ==================================== #
+def read_sensitivity_fits(filename):
+    """
+    Read sensitivity information from FITS file
+
+    Parameters
+    ----------
+    filename : str
+        Name of CSV file
+
+    Returns
+    -------
+    sensitivity : dict
+        Dictionary with sensitivity information
+    """
+    # Open FITS file
+    fits = gammalib.GFits(filename)
+
+    # Get sensitivity table
+    table = fits.table('SENSITIVITY')
+
+    # Get relevant columns
+    c_e_mean      = table['E_MEAN']
+    c_e_min       = table['E_MIN']
+    c_e_max       = table['E_MAX']
+    c_sensitivity = table['SENSITIVITY']
+
+    # Get mode
+    mode = table.card('TYPE').string()
+
+    # Initialise vectors
+    energies      = []
+    sensitivities = []
+
+    # Fill vectors
+    nrows = table.nrows()
+    for row in range(nrows):
+        if mode == 'Differential':
+            energies.append(c_e_mean[row])
+        else:
+            energies.append(c_e_min[row])
+        sensitivities.append(c_sensitivity[row])
+
+
+    # Set dictionary
+    sensitivity = {'mode'        : mode,
+                   'energy'      : energies,
+                   'sensitivity' : sensitivities}
+
+    # Return
+    return sensitivity
+
+
 # ===================== #
 # Plot sensitivity data #
 # ===================== #
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csphasecrv script.
 #
-# Copyright (C) 2017-2020 Juergen Knoedlseder
+# Copyright (C) 2017-2021 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
@@ -33,7 +33,7 @@ class Test(test):
     This test class makes unit tests for the csphasecrv script by using it
     from the command line and from Python.
     """
-    
+
     # Constructor
     def __init__(self):
         """
@@ -280,7 +280,7 @@ class Test(test):
         pcrv['phbinalg'] = 'LIN'
         pcrv['phbins']   = 2
         pcrv['method']   = '3D'
-        pcrv['enumbins'] = 0
+        pcrv['enumbins'] = 0     # Unbinned analysis
         pcrv['emin']     = 1.0
         pcrv['emax']     = 100.0
         pcrv['outfile']  = 'csphasecrv_py6.fits'
@@ -378,6 +378,6 @@ class Test(test):
         # Close FITS file
         fits.close()
-        
+
         # Return
         return
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the cssens script
 #
-# Copyright (C) 2016-2018 Juergen Knoedlseder
+# Copyright (C) 2016-2021 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
@@ -77,7 +77,7 @@ class Test(test):
                      ' srcname="Crab"'+ \
                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                      ' duration=180.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \
-                     ' outfile="cssens_cmd1.dat"'+ \
+                     ' outfile="cssens_cmd1.fits"'+ \
                      ' logfile="cssens_cmd1.log" chatter=1'
         # Check if execution was successful
@@ -85,14 +85,14 @@ class Test(test):
              'Check successful execution from command line')
         # Check result file
-        self._check_result_file('cssens_cmd1.dat')
+        self._check_result_file('cssens_cmd1.dat', 1)
         # Setup cssens command
         cmd = cssens+' inmodel="model_that_does_not_exist.xml"'+ \
                      ' srcname="Crab"'+ \
                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                      ' duration=180.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \
-                     ' outfile="cssens_cmd2.dat"'+ \
+                     ' outfile="cssens_cmd2.fits"'+ \
                      ' logfile="cssens_cmd2.log" debug=yes chatter=2'
         # Check if execution failed
@@ -118,12 +118,12 @@ class Test(test):
         sens['srcname']  = 'Crab'
         sens['caldb']    = self._caldb
         sens['irf']      = self._irf
-        sens['outfile']  = 'cssens_py1.dat'
         sens['duration'] = 180.0
         sens['rad']      = 3.0
         sens['emin']     = 1.0
         sens['emax']     = 10.0
         sens['bins']     = 2
+        sens['outfile']  = 'cssens_py1.fits'
         sens['logfile']  = 'cssens_py1.log'
         sens['chatter']  = 4
@@ -132,8 +132,8 @@ class Test(test):
         sens.run()
         sens.save()
-        # Check pull distribution file
-        self._check_result_file('cssens_py1.dat',nrows=3)
+        # Check sensitivity file
+        self._check_result_file('cssens_py1.fits', 2)
         # Set-up cssens for integral sensitivity computation without
         # multiprocessing
@@ -143,13 +143,13 @@ class Test(test):
         sens['srcname']  = 'Crab'
         sens['caldb']    = self._caldb
         sens['irf']      = self._irf
-        sens['outfile']  = 'cssens_py2.dat'
         sens['duration'] = 180.0
         sens['rad']      = 3.0
         sens['emin']     = 1.0
         sens['emax']     = 10.0
         sens['bins']     = 1
         sens['type']     = 'Integral'
+        sens['outfile']  = 'cssens_py2.fits'
         sens['logfile']  = 'cssens_py2.log'
         sens['chatter']  = 4
         sens['nthreads'] = 1
@@ -157,8 +157,8 @@ class Test(test):
         # Execute cssens script
         sens.execute()
-        # Check pull distribution file
-        self._check_result_file('cssens_py2.dat')
+        # Check sensitivity file
+        self._check_result_file('cssens_py2.fits', 1)
         # Return
         return
@@ -178,12 +178,12 @@ class Test(test):
         sens['srcname']  = 'Crab'
         sens['caldb']    = self._caldb
         sens['irf']      = self._irf
-        sens['outfile']  = 'cssens_py1_pickle.dat'
         sens['duration'] = 180.0
         sens['rad']      = 3.0
         sens['emin']     = 1.0
         sens['emax']     = 10.0
         sens['bins']     = 2
+        sens['outfile']  = 'cssens_py1_pickle.fits'
         sens['logfile']  = 'cssens_py1_pickle.log'
         sens['chatter']  = 4
@@ -196,23 +196,24 @@ class Test(test):
         obj.save()
         # Check result file
-        self._check_result_file('cssens_py1_pickle.dat',nrows=3)
+        self._check_result_file('cssens_py1_pickle.fits', 2)
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename, nrows=2):
+    def _check_result_file(self, filename, nrows):
         """
         Check result file
         """
-        # Open result file as CSV file
-        results = gammalib.GCsv(filename, ',')
+        # Open result file as FITS file
+        fits  = gammalib.GFits(filename)
+        table = fits.table('SENSITIVITY')
         # Check dimensions
-        self.test_value(results.nrows(), nrows,
+        self.test_value(table.nrows(), nrows,
              'Check number of rows in sensitivity file')
-        self.test_value(results.ncols(), 10,
+        self.test_value(table.ncols(), 10,
              'Check number of columns in sensitivity file')
         # Return