From 9dd6c7493b90309626e6ee3e4fd7853143fb6321 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 19 May 2021 17:26:33 +0200 Subject: [PATCH] Replace ASCII output files by FITS files in cstsdist (#1707) --- NEWS | 8 ++++---- cscripts/Makefile.am | 1 + cscripts/cspull.py | 1 - cscripts/cstsdist.par | 2 +- cscripts/cstsdist.py | 130 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------- doc/source/users/reference_manual/cspull.rst | 2 +- doc/source/users/reference_manual/cstsdist.rst | 11 +++++------ doc/source/users/reference_manual/ctbutterfly.rst | 2 +- examples/show_pull_evolution.py | 8 ++++---- examples/show_pull_histogram.py | 8 ++++---- examples/show_ts_distribution.py | 94 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------ test/data/ts.fits | 2 ++ test/test_cstsdist.py | 35 ++++++++++++++++++----------------- test/test_examples.py.in | 11 +++++++++-- 14 files changed, 254 insertions(+), 61 deletions(-) create mode 100644 test/data/ts.fits diff --git a/NEWS b/NEWS index a92a8c2..8f6d1c0 100755 --- a/NEWS +++ b/NEWS @@ -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 diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 7cf54b0..b13f5bf 100755 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -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 \ diff --git a/cscripts/cspull.py b/cscripts/cspull.py index cc86e97..32f9e86 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -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) diff --git a/cscripts/cstsdist.par b/cscripts/cstsdist.par index d38ee3c..92c8801 100755 --- a/cscripts/cstsdist.par +++ b/cscripts/cstsdist.par @@ -27,7 +27,7 @@ caldb, s, a, prod2,,, "Calibration database" 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 diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index 61fade5..58a69be 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -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 # diff --git a/doc/source/users/reference_manual/cspull.rst b/doc/source/users/reference_manual/cspull.rst index 7ee15d1..e73f1ee 100755 --- a/doc/source/users/reference_manual/cspull.rst +++ b/doc/source/users/reference_manual/cspull.rst @@ -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. diff --git a/doc/source/users/reference_manual/cstsdist.rst b/doc/source/users/reference_manual/cstsdist.rst index 05238c9..535f050 100755 --- a/doc/source/users/reference_manual/cstsdist.rst +++ b/doc/source/users/reference_manual/cstsdist.rst @@ -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. diff --git a/doc/source/users/reference_manual/ctbutterfly.rst b/doc/source/users/reference_manual/ctbutterfly.rst index 3c46a7d..5be96d8 100755 --- a/doc/source/users/reference_manual/ctbutterfly.rst +++ b/doc/source/users/reference_manual/ctbutterfly.rst @@ -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 diff --git a/examples/show_pull_evolution.py b/examples/show_pull_evolution.py index d012743..e15df05 100755 --- a/examples/show_pull_evolution.py +++ b/examples/show_pull_evolution.py @@ -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) diff --git a/examples/show_pull_histogram.py b/examples/show_pull_histogram.py index 2328672..1e9844e 100755 --- a/examples/show_pull_histogram.py +++ b/examples/show_pull_histogram.py @@ -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) diff --git a/examples/show_ts_distribution.py b/examples/show_ts_distribution.py index 5bf4998..dbfaa31 100755 --- a/examples/show_ts_distribution.py +++ b/examples/show_ts_distribution.py @@ -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='TS'): 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='TS'): # 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': diff --git a/test/data/ts.fits b/test/data/ts.fits new file mode 100644 index 0000000..338c034 --- /dev/null +++ b/test/data/ts.fits @@ -0,0 +1,2 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = 8 / number of bits per data pixel NAXIS = 0 / number of data axes EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 72 / width of table in bytes NAXIS2 = 10 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 9 / number of table fields TTYPE1 = 'TS ' / label for field 1 TFORM1 = '1D ' / data format of field 1 TTYPE2 = 'NEVENTS ' / label for field 2 TFORM2 = '1D ' / data format of field 2 TUNIT2 = 'counts ' / physical unit of field2 TTYPE3 = 'NPRED ' / label for field 3 TFORM3 = '1D ' / data format of field 3 TUNIT3 = 'counts ' / physical unit of field3 TTYPE4 = 'CRAB_PREFACTOR' / label for field 4 TFORM4 = '1D ' / data format of field 4 TTYPE5 = 'CRAB_PREFACTOR_ERROR' / label for field 5 TFORM5 = '1D ' / data format of field 5 TTYPE6 = 'CTABACKGROUNDMODEL_PREFACTOR' / label for field 6 TFORM6 = '1D ' / data format of field 6 TTYPE7 = 'CTABACKGROUNDMODEL_PREFACTOR_ERROR' / label for field 7 TFORM7 = '1D ' / data format of field 7 TTYPE8 = 'CTABACKGROUNDMODEL_INDEX' / label for field 8 TFORM8 = '1D ' / data format of field 8 TTYPE9 = 'CTABACKGROUNDMODEL_INDEX_ERROR' / label for field 9 TFORM9 = '1D ' / data format of field 9 EXTNAME = 'TS_DISTRIBUTION' / name of this extension INSTRUME= 'CTA ' / Name of Instrument TELESCOP= 'CTA ' / Name of Telescope CREATOR = 'cstsdist v2.0.0.dev' / Program which created the file NTRIALS = 10 / Number of trials STAT = 'DEFAULT ' / Optimization statistic EDISP = F / Use energy dispersion? TUNIT1 = ' ' / physical unit of field1 TUNIT4 = ' ' / physical unit of field4 TUNIT5 = ' ' / physical unit of field5 TUNIT6 = ' ' / physical unit of field6 TUNIT7 = ' ' / physical unit of field7 TUNIT8 = ' ' / physical unit of field8 TUNIT9 = ' ' / physical unit of field9 END @CNsܭ@@O<$5Zy?”?0 UHfS.?kEN@=@@<ώv