diff --git a/ChangeLog b/ChangeLog index 25233f6..336dc42 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -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 diff --git a/NEWS b/NEWS index ea0c796..2c162c9 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ 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 diff --git a/README.md b/README.md index 8465e42..63f1a13 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ 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/) diff --git a/cscripts/csphasecrv.py b/cscripts/csphasecrv.py index 587cd31..85febfb 100755 --- a/cscripts/csphasecrv.py +++ b/cscripts/csphasecrv.py @@ -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 # # ======================== # diff --git a/cscripts/cssens.par b/cscripts/cssens.par index 2eb2c64..fbf139f 100755 --- a/cscripts/cssens.par +++ b/cscripts/cssens.par @@ -23,7 +23,7 @@ caldb, s, a, prod2,,, "Calibration database" 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,,, "Overwrite existing output files with new output files?" debug, b, h, no,,, "Activate debugging mode?" mode, s, h, ql,,, "Mode of automatic parameters" logfile, f, h, cssens.log,,, "Log filename" - - diff --git a/cscripts/cssens.py b/cscripts/cssens.py index fb17e0d..f057967 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -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 # diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 25ffe35..485550d 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,7 +1,7 @@ .. _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 `_] - Remove implicit conversion warnings +* [`1707 `_] - + Replace ASCII output files by FITS files New features diff --git a/doc/source/users/reference_manual/cssens.rst b/doc/source/users/reference_manual/cssens.rst index 533592e..9443837 100755 --- a/doc/source/users/reference_manual/cssens.rst +++ b/doc/source/users/reference_manual/cssens.rst @@ -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). diff --git a/examples/show_sensitivity.py b/examples/show_sensitivity.py index ff4b94c..dbff23f 100755 --- a/examples/show_sensitivity.py +++ b/examples/show_sensitivity.py @@ -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 # # ===================== # diff --git a/test/test_csphasecrv.py b/test/test_csphasecrv.py index 67f4e06..f22a222 100755 --- a/test/test_csphasecrv.py +++ b/test/test_csphasecrv.py @@ -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 diff --git a/test/test_cssens.py b/test/test_cssens.py index 61fa1d6..4897742 100755 --- a/test/test_cssens.py +++ b/test/test_cssens.py @@ -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