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

Jürgen Knödlseder
1 parent d04d0a7a
Showing 14 changed files with 254 additions and 61 deletions
NEWS
cscripts/Makefile.am
cscripts/cspull.py
cscripts/cstsdist.par
cscripts/cstsdist.py
doc/source/users/reference_manual/cspull.rst
doc/source/users/reference_manual/cstsdist.rst
doc/source/users/reference_manual/ctbutterfly.rst
examples/show_pull_evolution.py
examples/show_pull_histogram.py
examples/show_ts_distribution.py
test/data/ts.fits
test/test_cstsdist.py
test/test_examples.py.in
@@ -307,7 +307,7 @@ None
  
 cstsdist - Computes Test Statistic distribution
 -----------------------------------------------
-None
+The ASCII output file was replaced by a FITS file (#1707).
  
  
 cstsmapmerge - Merge Test Statistic maps
@@ -361,9 +361,9 @@ regions (#3248).
  
 Examples
 --------
-The show_butterfly, show_pull_evolution.py, show_pull_histogram.py and
-show_sensitivity.py scripts were adapted to deal also with FITS output files
-(#1707).
+The show_butterfly, show_pull_evolution.py, show_pull_histogram.py,
+show_sensitivity.py and show_ts_distribution.py scripts were adapted to deal
+also with FITS output files (#1707).
  
  
 caldb
@@ -214,6 +214,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/residual_onoff.fits \
            $(top_srcdir)/test/data/lightcurve.fits \
            $(top_srcdir)/test/data/ts.dat \
+           $(top_srcdir)/test/data/ts.fits \
            $(top_srcdir)/test/data/models_spatial.xml \
            $(top_srcdir)/test/data/models_spectral.xml \
            $(top_srcdir)/test/data/model_filefunction.txt \
@@ -476,7 +476,6 @@ class cspull(ctools.csobservation):
             # If multiprocessing was used then recover results and put them
             # into the log file
             if self._nthreads > 1:
-                #result = poolresults[i][0]
                 results.append(poolresults[i][0])
                 self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
  
@@ -27,7 +27,7 @@ caldb,     s, a, prod2,,, &quot;Calibration database&quot;
 irf,       s, a, South_0.5h,,, "Instrument response function"
 edisp,     b, h, no,,, "Apply energy dispersion?"
 deadc,     r, h, 0.98,0,1, "Deadtime correction factor"
-outfile,   f, a, ts.dat,,, "Output Test Statistic distribution file"
+outfile,   f, a, ts.fits,,, "Output Test Statistic distribution FITS file"
  
 #
 # Script parameters
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates the TS distribution for a particular model
 #
-# Copyright (C) 2011-2018 Juergen Knoedlseder
+# Copyright (C) 2011-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
@@ -45,6 +45,7 @@ class cstsdist(ctools.csobservation):
  
         # Initialise some members
         self._srcname     = ''
+        self._fits        = None
         self._log_clients = False
         self._model       = None
         self._nthreads    = 0
@@ -65,6 +66,7 @@ class cstsdist(ctools.csobservation):
         # Set pickled dictionary
         state = {'base'        : ctools.csobservation.__getstate__(self),
                  'srcname'     : self._srcname,
+                 'fits'        : self._fits,
                  'log_clients' : self._log_clients,
                  'model'       : self._model,
                  'nthreads'    : self._nthreads}
@@ -83,6 +85,7 @@ class cstsdist(ctools.csobservation):
         """
         ctools.csobservation.__setstate__(self, state['base'])
         self._srcname     = state['srcname']
+        self._fits        = state['fits']
         self._log_clients = state['log_clients']
         self._model       = state['model']
         self._nthreads    = state['nthreads']
@@ -113,7 +116,10 @@ class cstsdist(ctools.csobservation):
         self['edisp'].boolean()
         self['ntrials'].integer()
         self['debug'].boolean()
-        self['outfile'].filename()
+
+        # Read ahead output parameters
+        if self._read_ahead():
+            self['outfile'].filename()
  
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
@@ -261,14 +267,14 @@ class cstsdist(ctools.csobservation):
                 if par.is_free():
  
                     # Set parameter name
-                    name = model_name+"_"+par.name()
+                    name = model_name + '_' + par.name()
  
                     # Append value
                     colnames.append(name)
                     values[name] = par.value()
  
                     # Append error
-                    name = 'e_'+name
+                    name = name+'_error'
                     colnames.append(name)
                     values[name] = par.error()
  
@@ -278,6 +284,75 @@ class cstsdist(ctools.csobservation):
         # Return
         return result
  
+    def _create_fits(self, results):
+        """
+        Create FITS file from results
+
+        Parameters
+        ----------
+        results : list of dict
+            List of result dictionaries
+        """
+        # Gather headers for parameter columns
+        headers = []
+        for colname in results[0]['colnames']:
+            if colname != 'TS' and colname != 'Nevents' and \
+               colname != 'Npred':
+                headers.append(colname)
+
+        # Create FITS table columns
+        nrows   = len(results)
+        ts      = gammalib.GFitsTableDoubleCol('TS', nrows)
+        nevents = gammalib.GFitsTableDoubleCol('NEVENTS', nrows)
+        npred   = gammalib.GFitsTableDoubleCol('NPRED', nrows)
+        ts.unit('')
+        nevents.unit('counts')
+        npred.unit('counts')
+        columns = []
+        for header in headers:
+            name   = gammalib.toupper(header)
+            column = gammalib.GFitsTableDoubleCol(name, nrows)
+            column.unit('')
+            columns.append(column)
+
+        # Fill FITS table columns
+        for i, result in enumerate(results):
+            ts[i]      = result['values']['TS']
+            nevents[i] = result['values']['Nevents']
+            npred[i]   = result['values']['Npred']
+            for k, column in enumerate(columns):
+                column[i] = result['values'][headers[k]]
+
+        # Initialise FITS Table with extension "TS_DISTRIBUTION"
+        table = gammalib.GFitsBinTable(nrows)
+        table.extname('TS_DISTRIBUTION')
+
+        # 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('NTRIALS', self['ntrials'].integer(),  'Number of trials')
+        table.card('STAT',    self['statistic'].string(), 'Optimization statistic')
+        table.card('EDISP',   self['edisp'].boolean(),    'Use energy dispersion?')
+
+        # Append filled columns to fits table
+        table.append(ts)
+        table.append(nevents)
+        table.append(npred)
+        for column in columns:
+            table.append(column)
+
+        # Create the FITS file now
+        self._fits = gammalib.GFits()
+        self._fits.append(table)
+
+        # Return
+        return
+
  
     # Public methods
     def run(self):
@@ -306,6 +381,9 @@ class cstsdist(ctools.csobservation):
         # Get number of trials
         ntrials = self['ntrials'].integer()
  
+        # Initialise results
+        results  = []
+
         # If more than a single thread is requested then use multiprocessing
         if self._nthreads > 1:
             args        = [(self, '_trial', i) for i in range(ntrials)]
@@ -317,16 +395,20 @@ class cstsdist(ctools.csobservation):
             # If multiprocessing was used then recover results and put them
             # into the log file
             if self._nthreads > 1:
-                result = poolresults[i][0]
+                results.append(poolresults[i][0])
                 self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
  
             # ... otherwise make a trial now
             else:
+
+                # Run trial
                 result = self._trial(i)
  
-            # Write out trial result
-            ioutils.write_csv_row(self['outfile'].filename().url(), i,
-                                  result['colnames'], result['values'])
+                # Append results
+                results.append(result)
+
+        # Create FITS file
+        self._create_fits(results)
  
         # Return
         return
@@ -346,6 +428,38 @@ class cstsdist(ctools.csobservation):
         # Return
         return
  
+    def save(self):
+        """
+        Save TS distribution FITS file
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save TS distribution')
+
+        # 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, 'TS distribution file', outfile.url())
+
+            # Save TS distribution
+            self._fits.saveto(outfile, self['clobber'].boolean())
+
+        # Return
+        return
+
+    def ts_distribution(self):
+        """
+        Return TS distribution FITS file
+
+        Returns:
+            FITS file containing TS distribution
+        """
+        # Return
+        return self._fits
+
  
 # ======================== #
 # Main routine entry point #
@@ -74,7 +74,7 @@ General parameters
     Average deadtime correction factor.
  
 ``outfile [file]``
-    ASCII file containing the individual pull values.
+    Output pull distribution FITS ASCII file.
  
 ``ntrials [integer]``
     Number of samples for generating the pull distribution.
@@ -20,11 +20,10 @@ be performed in parallel over as many processes as the number of CPUs available 
 machine. The maximum number of parallel processes can be set by the user through the
 ``nthreads`` hidden parameter.
  
-:ref:`cstsdist` will create an ASCII file in comma-separated value (CSV) format,
-containing one row per TS computation. The first row is a header row providing
-the column names. The following rows give the TS value, the log-likelihood 
-values of the fit with or without the source, the number of observed and 
-fitted events, as well as the values and errors for all fitted parameters.
+:ref:`cstsdist` will create on output a FITS file containing a table with one row
+per TS computation. The columns in the table comprise the TS value, the number
+of observed and fitted events, and the values and errors for all fitted
+parameters.
  
 From the output file, TS distribution plots can be generated using for
 example the ``show_ts_distribution.py`` script in the examples folder. The
@@ -74,7 +73,7 @@ General parameters
     Average deadtime correction factor.
  
 ``outfile [file]``
-    Output ASCII file containing the TS distribution values.
+    Output TS distribution FITS file.
  
 ``ntrials [integer]``
     Number of Monte Carlo samples.
@@ -68,7 +68,7 @@ General parameters
     Applies energy dispersion to response computation.
  
 ``outfile [file]``
-    Output butterfly ASCII file name.
+    Output butterfly FITS file.
  
 ``(fit = no) [boolean]``
     Performs maximum likelihood fitting of input model ignoring any provided
@@ -51,8 +51,8 @@ def read_pull_distribution(filename, parname):
  
     Returns
     -------
-    pulls : dict
-        Dictionary with pull distribution
+    pulls : list of floats
+        List with pull distribution
     """
     # Read pull distribution dependent on file type
     fname = gammalib.GFilename(filename)
@@ -81,8 +81,8 @@ def read_pull_distribution_fits(filename, parname):
  
     Returns
     -------
-    pulls : dict
-        Dictionary with pull distribution
+    pulls : list of floats
+        List with pull distribution
     """
     # Open FITS file
     fits = gammalib.GFits(filename)
@@ -51,8 +51,8 @@ def read_pull_distribution(filename, parname):
  
     Returns
     -------
-    pulls : dict
-        Dictionary with pull distribution
+    pulls : list of floats
+        List with pull distribution
     """
     # Read pull distribution dependent on file type
     fname = gammalib.GFilename(filename)
@@ -81,8 +81,8 @@ def read_pull_distribution_fits(filename, parname):
  
     Returns
     -------
-    pulls : dict
-        Dictionary with pull distribution
+    pulls : list of floats
+        List with pull distribution
     """
     # Open FITS file
     fits = gammalib.GFits(filename)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Display TS distribution generated by cstsdist
 #
-# Copyright (C) 2011-2020 Jurgen Knodlseder
+# Copyright (C) 2011-2021 Jurgen Knodlseder
 #
 # 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,13 +33,84 @@ try:
 except ImportError:
     print('This script needs the "numpy" module')
     sys.exit()
+import gammalib
 import cscripts
  
  
 # ==================== #
 # Read TS distribution #
 # ==================== #
-def read_ts(filename, tsname='TS'):
+def read_ts_distribution(filename):
+    """
+    Read TS distribution
+
+    Parameters
+    ----------
+    filename : str
+        Name of FITS or CSV file
+    parname : str
+        Parameter name
+
+    Returns
+    -------
+    ts : list of floats
+        List with TS distribution
+    """
+    # Read TS distribution dependent on file type
+    fname = gammalib.GFilename(filename)
+    if fname.is_fits():
+        ts = read_ts_distribution_fits(filename)
+    else:
+        ts = read_ts_distribution_csv(filename)
+
+    # Return ts
+    return ts
+
+
+# =================================== #
+# Read TS distribution from FITS file #
+# =================================== #
+def read_ts_distribution_fits(filename):
+    """
+    Read TS distribution from FITS file
+
+    Parameters
+    ----------
+    filename : str
+        Name of FITS file
+    parname : str
+        Parameter name
+
+    Returns
+    -------
+    ts : list of floats
+        List with TS distribution
+    """
+    # Open FITS file
+    fits = gammalib.GFits(filename)
+
+    # Get TS distribution table
+    table = fits.table('TS_DISTRIBUTION')
+
+    # Get relevant column
+    column = table['TS']
+
+    # Initialise vector
+    ts = []
+
+    # Fill vectors
+    nrows = table.nrows()
+    for row in range(nrows):
+        ts.append(column[row])
+
+    # Return TS distribution
+    return ts
+
+
+# ==================== #
+# Read TS distribution #
+# ==================== #
+def read_ts_distribution_csv(filename, tsname='TS'):
     """
     Read TS distribution values from CSV file
  
@@ -52,8 +123,8 @@ def read_ts(filename, tsname=&#39;TS&#39;):
  
     Returns
     -------
-    a : `numpy.array`
-        Array of TS values
+    values : list of floats
+        List with TS distribution
     """
     # Initialise list
     values = []
@@ -82,11 +153,8 @@ def read_ts(filename, tsname=&#39;TS&#39;):
         # Flag that first row has been passed
         first = False
  
-    # Create numpy array
-    a = np.array(values)
-
-    # Return array
-    return a
+    # Return TS values
+    return values
  
  
 # =========== #
@@ -309,9 +377,11 @@ def show_ts_distribution():
     disttype = options[4]['value']
     plotfile = options[5]['value']
  
-    # Read values from CSV file
-    values = read_ts(args[0], tsname=tsname)
-    print('%d values read.' % len(values))
+    # Read TS distribution from file
+    ts = read_ts_distribution(args[0])
+
+    # Create Numpy array
+    values = np.array(ts)
  
     # Show histogram
     if disttype == 'pdf':
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the cstsdist 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
@@ -83,7 +83,7 @@ class Test(test):
                        ' inmodel="'+self._model+'" srcname="Crab"'+ \
                        ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                        ' ntrials=2 '+ \
-                       ' outfile="cstsdist_cmd1.dat"'+ \
+                       ' outfile="cstsdist_cmd1.fits"'+ \
                        ' logfile="cstsdist_cmd1.log" chatter=2'
  
         # Check if execution was successful
@@ -91,14 +91,14 @@ class Test(test):
              'Check successful execution from command line')
  
         # Check result file
-        self._check_result_file('cstsdist_cmd1.dat')
+        self._check_result_file('cstsdist_cmd1.fits')
  
         # Setup cstsdist command
         cmd = cstsdist+' inobs="event_file_that_does_not_exist.xml"'+ \
                        ' inmodel="'+self._model+'" srcname="Crab"'+ \
                        ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                        ' ntrials=2 '+ \
-                       ' outfile="cstsdist_cmd2.dat"'+ \
+                       ' outfile="cstsdist_cmd2.fits"'+ \
                        ' logfile="cstsdist_cmd2.log" chatter=2'
  
         # Check if execution failed
@@ -124,7 +124,7 @@ class Test(test):
         tsdist['caldb']   = self._caldb
         tsdist['irf']     = self._irf
         tsdist['ntrials'] = 2
-        tsdist['outfile'] = 'cstsdist_py1.dat'
+        tsdist['outfile'] = 'cstsdist_py1.fits'
         tsdist['logfile'] = 'cstsdist_py1.log'
         tsdist['chatter'] = 3
  
@@ -134,7 +134,7 @@ class Test(test):
         tsdist.save()
  
         # Check TS distribution file
-        self._check_result_file('cstsdist_py1.dat')
+        self._check_result_file('cstsdist_py1.fits')
  
         # Return
         return
@@ -154,7 +154,7 @@ class Test(test):
         tsdist['caldb']   = self._caldb
         tsdist['irf']     = self._irf
         tsdist['ntrials'] = 2
-        tsdist['outfile'] = 'cstsdist_py2.dat'
+        tsdist['outfile'] = 'cstsdist_py2.fits'
         tsdist['logfile'] = 'cstsdist_py2.log'
         tsdist['chatter'] = 4
  
@@ -163,7 +163,7 @@ class Test(test):
         tsdist.execute()
  
         # Check TS distribution file
-        self._check_result_file('cstsdist_py2.dat')
+        self._check_result_file('cstsdist_py2.fits')
  
         # Return
         return
@@ -178,7 +178,7 @@ class Test(test):
         tsdist['inmodel'] = self._stacked_model
         tsdist['srcname'] = 'Crab'
         tsdist['ntrials'] = 2
-        tsdist['outfile'] = 'cstsdist_py3.dat'
+        tsdist['outfile'] = 'cstsdist_py3.fits'
         tsdist['logfile'] = 'cstsdist_py3.log'
         tsdist['chatter'] = 4
  
@@ -187,7 +187,7 @@ class Test(test):
         tsdist.execute()
  
         # Check TS distribution file
-        self._check_result_file('cstsdist_py3.dat')
+        self._check_result_file('cstsdist_py3.fits')
  
         # Return
         return
@@ -209,7 +209,7 @@ class Test(test):
         tsdist['irf']      = self._irf
         tsdist['ntrials']  = 2
         tsdist['nthreads'] = 1 # This tests at the same time no multiprocessing
-        tsdist['outfile']  = 'cstsdist_py1_pickle.dat'
+        tsdist['outfile']  = 'cstsdist_py1_pickle.fits'
         tsdist['logfile']  = 'cstsdist_py1_pickle.log'
         tsdist['chatter']  = 3
  
@@ -222,13 +222,13 @@ class Test(test):
         obj.save()
  
         # Check TS distribution file
-        self._check_result_file('cstsdist_py1_pickle.dat')
+        self._check_result_file('cstsdist_py1_pickle.fits')
  
         # Return
         return
  
     # Check result file
-    def _check_result_file(self, filename, nrows=3, ncols=9):
+    def _check_result_file(self, filename, nrows=2, ncols=9):
         """
         Check result file
  
@@ -241,12 +241,13 @@ class Test(test):
         ncols : int, optional
             Required number of columns
         """
-        # Open result file as CSV file
-        results = gammalib.GCsv(filename, ',')
+        # Open pull file as FITS file
+        fits  = gammalib.GFits(filename)
+        table = fits.table('TS_DISTRIBUTION')
  
         # Check dimensions
-        self.test_value(results.nrows(), nrows, 'Check rows in TS file')
-        self.test_value(results.ncols(), ncols, 'Check columns in TS file')
+        self.test_value(table.nrows(), nrows, 'Check rows in TS file')
+        self.test_value(table.ncols(), ncols, 'Check columns in TS file')
  
         # Return
         return
@@ -460,14 +460,21 @@ class Test(gammalib.GPythonTestSuite):
         Test show_ts_distribution
         """
         # Set arguments
-        args = ('-n 100 -p example_show_ts_distribution1.eps %s/ts.dat' %
+        args = ('-n 100 -p example_show_ts_distribution_csv1.eps %s/ts.dat' %
                 self._datadir)
  
         # Execute script
         self._execute_python('show_ts_distribution', args=args)
  
         # Set arguments
-        args = ('-n 100 -t cdf -p example_show_ts_distribution2.eps %s/ts.dat' %
+        args = ('-n 100 -t cdf -p example_show_ts_distribution_csv2.eps %s/ts.dat' %
+                self._datadir)
+
+        # Execute script
+        self._execute_python('show_ts_distribution', args=args)
+
+        # Set arguments
+        args = ('-n 100 -p example_show_ts_distribution_fits.eps %s/ts.fits' %
                 self._datadir)
  
         # Execute script