diff --git a/NEWS b/NEWS index 45e4473..3dc8aec 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.2.0 -12 December 2016 +26 January 2017 Introduction @@ -195,6 +195,9 @@ in particular the "PhotonFlux" and "EnergyFlux" names (#1832). The class now derives from ctlikelihood (#1828). +The tool now makes sure that the intensity parameter of the source of interest +is free in the model optimisation that preceeds the upper limit computation. + cterror - Likelihood profile error computation ---------------------------------------------- @@ -235,10 +238,14 @@ Add cssrcdetect script to detect sources in a sky map (#1617). The "ds9file" parameters were changed to "outds9file" to be consistent with the naming conventions. +Improve the upper limit computation in the cslightcrv script. Now all source +parameters, except the source intensity, are set to the initial values given +in the model definition file for upper limit computation. + Examples -------- -None +Add show_lightcurve script to display a light curve. caldb diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 685ec45..401a77f 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -166,6 +166,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/sensitivity_integral.dat \ $(top_srcdir)/test/data/skymap.fits \ $(top_srcdir)/test/data/spectrum.fits \ + $(top_srcdir)/test/data/lightcurve.fits \ $(top_srcdir)/test/data/ts.dat \ $(top_srcdir)/test/data/models_spatial.xml \ $(top_srcdir)/test/data/models_spectral.xml \ diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index 33a2bd6..a583f38 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -2,7 +2,7 @@ # ========================================================================== # Generates a lightcurve. # -# 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 @@ -69,7 +69,7 @@ class cslightcrv(ctools.cscript): """ # Set name self._name = 'cslightcrv' - self._version = '1.2.0' + self._version = '1.3.0' # Initialise some members self._srcname = '' @@ -335,18 +335,26 @@ class cslightcrv(ctools.cscript): table.append(ts) # Append upper limit column "UpperLimit" - ulimit = gammalib.GFitsTableDoubleCol('UpperLimit', nrows) - ulimit.unit('ph/cm2/s') + ul_diff = gammalib.GFitsTableDoubleCol('DiffUpperLimit', nrows) + ul_flux = gammalib.GFitsTableDoubleCol('FluxUpperLimit', nrows) + ul_eflux = gammalib.GFitsTableDoubleCol('EFluxUpperLimit', nrows) + ul_diff.unit('ph/cm2/s/MeV') + ul_flux.unit('ph/cm2/s') + ul_eflux.unit('erg/cm2/s') for i, result in enumerate(results): - ulimit[i] = result['ulimit'] - table.append(ulimit) + ul_diff[i] = result['ul_diff'] + ul_flux[i] = result['ul_flux'] + ul_eflux[i] = result['ul_eflux'] + table.append(ul_diff) + table.append(ul_flux) + table.append(ul_eflux) # Return table return table def _compute_ulimit(self, obs): """ - Computes upper flux limit + Computes upper limit Parameters ---------- @@ -355,35 +363,57 @@ class cslightcrv(ctools.cscript): Returns ------- - ulimit : float - Upper flux limit (-1.0 of not computed) + ul_diff, ul_flux, ul_eflux : tuple of float + Upper limits, -1.0 of not computed """ # Initialise upper flux limit - ulimit_value = -1.0 - + ul_diff = -1.0 + ul_flux = -1.0 + ul_eflux = -1.0 + # Perform computation only if requested if self['calc_ulim'].boolean(): # Write header in logger if self._logExplicit(): - self._log.header3('Computing upper flux limit') - - # Create upper limit object - ulimit = ctools.ctulimit(obs) - ulimit['srcname'] = self._srcname - ulimit['eref'] = 1.0 - - # Try to run upper limit and catch any exceptions + self._log.header3('Computing upper limit') + + # Create copy of observations + cpy_obs = obs.copy() + + # Fix parameters of source of interest in copy of observations. + # This assures that the original spectral parameters and position + # are used in the upper limit search. The ctulimit tool makes sure + # that the flux parameter is freed when needed. + source = cpy_obs.models()[self._srcname] + for par in source: + if par.is_free(): + par.fix() + + # Create instance of upper limit tool + ulimit = ctools.ctulimit(cpy_obs) + ulimit['srcname'] = self._srcname + ulimit['eref'] = 1.0 + ulimit['emin'] = self['emin'].real() + ulimit['emax'] = self['emax'].real() + ulimit['sigma_min'] = 0.0 + ulimit['sigma_max'] = 3.0 + + # Try to run upper limit tool and catch any exceptions try: ulimit.run() - ulimit_value = ulimit.flux_ulimit() + ul_diff = ulimit.diff_ulimit() + ul_flux = ulimit.flux_ulimit() + ul_eflux = ulimit.eflux_ulimit() except: if self._logTerse(): - self._log('Upper limit flux calculation failed.\n') - ulimit_value = -1.0 + self._log('Upper limit calculation failed.\n') + ul_diff = -1.0 + ul_flux = -1.0 + ul_eflux = -1.0 - # Return upper limit - return ulimit_value + # Return upper limit tuple + return ul_diff, ul_flux, ul_eflux def _bin_observation(self, obs): """ @@ -517,7 +547,9 @@ class cslightcrv(ctools.cscript): result = {'mjd': tmean.mjd(), 'e_mjd': twidth / gammalib.sec_in_day, 'ts': 0.0, - 'ulimit': 0.0, + 'ul_diff': 0.0, + 'ul_flux': 0.0, + 'ul_eflux': 0.0, 'pars': pars, 'values': {}} @@ -581,9 +613,12 @@ class cslightcrv(ctools.cscript): result['values']['e_'+par] = source[par].error() # Calculate upper limit (-1 if not computed) - ulimit_value = self._compute_ulimit(like.obs()) - if ulimit_value > 0.0: - result['ulimit'] = ulimit_value + #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs()) + ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs) + if ul_diff > 0.0: + result['ul_diff'] = ul_diff + result['ul_flux'] = ul_flux + result['ul_eflux'] = ul_eflux # Extract Test Statistic value if self['calc_ts'].boolean(): @@ -601,9 +636,13 @@ class cslightcrv(ctools.cscript): unit = source[par].unit() self._log_value(gammalib.NORMAL, par, str(value)+' +/- '+str(error)+' '+unit) - if result['ulimit'] > 0.0: + if ul_diff > 0.0: + self._log_value(gammalib.NORMAL, 'Upper flux limit', + str(result['ul_diff'])+' ph/cm2/s/MeV') + self._log_value(gammalib.NORMAL, 'Upper flux limit', + str(result['ul_flux'])+' ph/cm2/s') self._log_value(gammalib.NORMAL, 'Upper flux limit', - str(result['ulimit'])+' ph/cm2/s') + str(result['ul_eflux'])+' erg/cm2/s') if self['calc_ts'].boolean(): self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts']) diff --git a/cscripts/csroot2caldb.py b/cscripts/csroot2caldb.py index a2b15ba..867788f 100755 --- a/cscripts/csroot2caldb.py +++ b/cscripts/csroot2caldb.py @@ -146,7 +146,7 @@ class csroot2caldb(ctools.cscript): def _make_dirs(self, version, irf): """ Generate CALDB directory structure for one observation identifier - + The structure is given by data///bcf/ @@ -227,7 +227,7 @@ class csroot2caldb(ctools.cscript): def _open(self, irf, ds): """ Open existing or create new calibration - + The actual version will put all calibrations in the same file, although each part of the response function will have its own logical name. We can thus easily modify the script to put each calibration information @@ -322,7 +322,7 @@ class csroot2caldb(ctools.cscript): def _open_hdu(self, fits, extname, name, doc, irf): """ Open HDU - + Opens a FITS binary table with given "extname". If HDU does not exist in the FITS file it will be created and appended to the FITS file. @@ -534,10 +534,9 @@ class csroot2caldb(ctools.cscript): def _root2ea(self, tfile, irf, ds): """ Translate ROOT to CALDB effective area extension - + The following ROOT histograms are used: - EffectiveAreaEtrue_offaxis -> EFFAREA - - EffectiveArea_offaxis -> EFFAREA_RECO Parameters ---------- @@ -555,16 +554,13 @@ class csroot2caldb(ctools.cscript): # Get relevant ROOT histograms etrue = tfile.Get('EffectiveAreaEtrue_offaxis') - ereco = tfile.Get('EffectiveArea_offaxis') # If requested then normalize the 2D histogram on the on-axis 1D # histogram. This assures that the 2D histogram has the same on-axis # effective area dependence as the 1D histogram. if self['norm1d'].boolean(): etrue_1D = tfile.Get('EffectiveAreaEtrue') - ereco_1D = tfile.Get('EffectiveArea') self._renorm_onaxis(etrue, etrue_1D) - self._renorm_onaxis(ereco, ereco_1D) # If requested then rebin the effective area Etrue histogram. This # increases the number of bins by a factor 10. @@ -589,17 +585,6 @@ class csroot2caldb(ctools.cscript): for ieng in range(neng): value = etrue.GetBinContent(ieng+1,ioff+1) * eascale etrue.SetBinContent(ieng+1,ioff+1,value) - neng = ereco.GetXaxis().GetNbins() - noffset = ereco.GetYaxis().GetNbins() - for ioff in range(noffset): - for ieng in range(neng): - value = ereco.GetBinContent(ieng+1,ioff+1) * eascale - ereco.SetBinContent(ieng+1,ioff+1,value) - - # Set boundaries (use Ereco boundaries) - #bounds = self._make_2D(ereco, ds['HDU_EA'], None, 'm2') - #for b in bounds: - # irf['EA_BOUNDS'].append(b) # Write boundary keywords self._set_cif_keywords(ds['HDU_EA'], irf['EA_NAME'], @@ -608,16 +593,13 @@ class csroot2caldb(ctools.cscript): # Create "EFFAREA" data column self._make_2D(etrue, ds['HDU_EA'], 'EFFAREA', 'm2') - # Create "EFFAREA_RECO" data column - self._make_2D(ereco, ds['HDU_EA'], 'EFFAREA_RECO', 'm2') - # Return return def _root2psf(self, tfile, irf, ds): """ Translate ROOT to CALDB point spread function extension - + Parameters ---------- tfile : `~ROOT.TFile` @@ -646,7 +628,7 @@ class csroot2caldb(ctools.cscript): def _root2psf_gauss(self, tfile, irf, ds): """ Translate ROOT to CALDB point spread function extension - + The following ROOT histograms are used: - 1/(2*pi*SIGMA_1) -> SCALE - AngRes_offaxis -> SIGMA_1 (scaling: 1/0.8) @@ -729,7 +711,7 @@ class csroot2caldb(ctools.cscript): def _root2psf_king(self, tfile, irf, ds): """ Translate ROOT to CALDB point spread function extension - + The following ROOT histograms are used: - AngRes_offaxis - AngRes80_offaxis @@ -853,7 +835,7 @@ class csroot2caldb(ctools.cscript): def _root2edisp(self, tfile, irf, ds): """ Translate ROOT to CALDB energy dispersion extension - + The following ROOT histograms are used: - EestOverEtrue_offaxis -> MATRIX @@ -892,7 +874,7 @@ class csroot2caldb(ctools.cscript): def _root2bgd(self, tfile, irf, ds): """ Translate ROOT to CALDB background extension. - + The following ROOT histograms are used: - BGRatePerSqDeg_offaxis -> BGD @@ -945,7 +927,7 @@ class csroot2caldb(ctools.cscript): def _append_column_axis(self, hdu, name, unit, axis, log=False): """ Append column of ROOT axis values to HDU - + Parameters ---------- hdu : `~gammalib.GFitsHDU` @@ -991,7 +973,7 @@ class csroot2caldb(ctools.cscript): def _append_column_values(self, hdu, name, unit, values): """ Append column of values to HDU - + Parameters ---------- hdu : `~gammalib.GFitsHDU` @@ -1089,7 +1071,7 @@ class csroot2caldb(ctools.cscript): for ioff in range(noffset): value = hist2D_copy.GetBinContent(ieng+1,ioff+1) / onaxis hist2D_copy.SetBinContent(ieng+1,ioff+1,value) - + # Put back the energy dependence for ieng in range(neng): onaxis = hist1D.GetBinContent(ieng+1) @@ -1324,7 +1306,7 @@ class csroot2caldb(ctools.cscript): def _make_2D(self, array, hdu, name, unit, scale=1.0): """ Make 2D data column as function of energy and offset angle - + If the HDU has already the energy and offset angle columns, this method will simply add another data column. If name==None, the method will not append any data column. @@ -1381,7 +1363,7 @@ class csroot2caldb(ctools.cscript): def _make_3D(self, array, hdu, name, unit): """ Make 3D data column as function of DETX, DETY and energy - + If the HDU has already the energy and offset angle columns, this method will simply add another data column. If name==None, the method will not append any data column. @@ -1471,7 +1453,7 @@ class csroot2caldb(ctools.cscript): def _make_3D_migra(self, array, hdu, name, unit, scale=1.0): """ Make 3D data column as function of ETRUE, MIGRA and THETA - + If the HDU has already the energy and offset angle columns, this method will simply add another data column. If name==None, the method will not append any data column. diff --git a/dev/macosx-package.sh b/dev/macosx-package.sh new file mode 100755 index 0000000..0e146ea --- /dev/null +++ b/dev/macosx-package.sh @@ -0,0 +1,233 @@ +#!/bin/bash -f +# ========================================================================== +# ctools Mac OS X package creation +# +# Copyright (C) 2017 Sylvie Brau-Nogué +# +# 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 . +# +# ========================================================================== + +# ============================================= # +# Get ctools/GammaLib version from command line # +# ============================================= # +if [ $# -ne 1 ] ; then + echo "Please specify ctools/GammaLib version in the form x.y.z" + exit 1 +fi +VERSION=$1 + + +# =============================== # +# Set software component versions # +# =============================== # +CFITSIO=cfitsio3410 +NCURSES=ncurses-6.0 +READLINE=readline-7.0 +GAMMALIB=gammalib-$VERSION +CTOOLS=ctools-$VERSION + + +# ============== # +# Set parameters # +# ============== # +WRKDIR=$PWD/osxpackage +INSTALLDIR=$WRKDIR/install +SRCDIR=$WRKDIR/src +PKGDIR=$WRKDIR/pkg +PRODDIR=$WRKDIR/prod +PLISTFILE=$PKGDIR/$CTOOLS-components.plist +DISTFILE=$PRODDIR/$CTOOLS.dist +DMGFILE=$WRKDIR/$CTOOLS-macosx10.7.dmg + + +# ====================== # +# Clean package creation # +# ====================== # +rm -rf $INSTALLDIR +rm -rf $PKGDIR +rm -rf $PRODDIR +rm -rf $DMGFILE + + +# ============================= # +# Create package directory tree # +# ============================= # +mkdir -p $SRCDIR +mkdir -p $PKGDIR +mkdir -p $PRODDIR + + +# Make sure that wget is installed +# TODO: If not installed, execute "brew install wget" + + +# =============== # +# Install ncurses # +# =============== # +cd $SRCDIR +if [ ! -d "$NCURSES" ]; then + wget https://ftp.gnu.org/pub/gnu/ncurses/$NCURSES.tar.gz + tar xvfz $NCURSES.tar.gz +fi +cd $NCURSES +./configure --prefix=$INSTALLDIR --with-shared +make +make install + + +# ================ # +# Install readline # +# ================ # +cd $SRCDIR +if [ ! -d "$READLINE" ]; then + wget https://ftp.gnu.org/pub/gnu/readline/$READLINE.tar.gz + tar xvfz $READLINE.tar.gz +fi +cd $READLINE +./configure --prefix=$INSTALLDIR +make +make install + + +# =============== # +# Install cfitsio # +# =============== # +cd $SRCDIR +if [ ! -d "cfitsio" ]; then + wget http://heasarc.gsfc.nasa.gov/FTP/software/fitsio/c/$CFITSIO.tar.gz + tar xvfz $CFITSIO.tar.gz +fi +cd cfitsio +./configure --prefix=$INSTALLDIR +make shared +make install + + +# ================ # +# Install GammaLib # +# ================ # +cd $SRCDIR +if [ ! -d "$GAMMALIB" ]; then + wget http://cta.irap.omp.eu/ctools/releases/gammalib/$GAMMALIB.tar.gz + tar xvfz $GAMMALIB.tar.gz +fi +cd $GAMMALIB +./configure --prefix=$INSTALLDIR +make +make install + + +# ============== # +# Install ctools # +# ============== # +cd $SRCDIR +if [ ! -d "$CTOOLS" ]; then + wget http://cta.irap.omp.eu/ctools/releases/ctools/$CTOOLS.tar.gz + tar xvfz $CTOOLS.tar.gz +fi +cd $CTOOLS +./configure --prefix=$INSTALLDIR +make +make install + + +# ====================== # +# Build Mac OS X package # +# ====================== # +pkgbuild --analyze --root $INSTALLDIR $PLISTFILE +pkgbuild --identifier cta.irap.omp.eu.ctools.pkg \ + --version $VERSION \ + --root $INSTALLDIR \ + --install-location /usr/local/gamma \ + --component-plist $PLISTFILE \ + $PKGDIR/$CTOOLS.pkg + + +# ====================== # +# Build Mac OS X product # +# ====================== # +# Analyse existing package +productbuild --synthesize \ + --package $PKGDIR/$CTOOLS.pkg \ + $DISTFILE + +# Add information to dist file +sed -i '' '$ i \ +\ +'$CTOOLS'\ +irap.omp.eu\ +CTA Science Analysis Tools\ +' "$PRODDIR/$CTOOLS.dist" + +# Build product +productbuild --distribution $DISTFILE \ + --version $VERSION \ + --resources $SRCDIR/$CTOOLS \ + --package-path $PKGDIR \ + $PRODDIR/$CTOOLS.pkg + +# Add additional files to production folder +cp $SRCDIR/$CTOOLS/COPYING $PRODDIR/License.txt + +# Add ReadMe +/bin/cat <$PRODDIR/ReadMe.txt +This package will install + + ctools-$VERSION and gammalib-$VERSION + +for Mac OS X 10.7 or later. + +Installation requires approximately 260 MB of disk space. The package +will be installed in the directory /usr/local/gamma, hence you need +administrator privileges to do the install. + +Simply click on the installer package and follow the instructions. +The installation will take less than one minute. + +After installing the package, you have to initialise the ctools +environment. This will be done by invoking an initialisation script +that can be found in the /usr/local/gamma/bin directory. + +To invoke this script you need to add the following to your .bashrc +or $HOME/.profile script: + + export CTOOLS=/usr/local/gamma + source \$CTOOLS/bin/gammalib-init.sh + +If you use C shell or a variant then add the following to your +.cshrc or .tcshrc script: + + setenv CTOOLS /usr/local/gamma + source \$CTOOLS/bin/gammalib-init.csh + +More information on ctools can be found at +http://cta.irap.omp.eu/ctools/. + +If you encounter problems in installing or using this package, please +send your error report to jurgen.knodlseder@irap.omp.eu. In your report, +please specify your Mac OS X version and Python version and the problem +you're facing. +EOM + + + +# ==================== # +# Create Mac OS X disk # +# ==================== # +hdiutil create -volname $CTOOLS \ + -srcfolder $PRODDIR/$CTOOLS.pkg \ + -srcfolder $PRODDIR/License.txt \ + -srcfolder $PRODDIR/ReadMe.txt \ + $DMGFILE diff --git a/examples/Makefile.am b/examples/Makefile.am index ad47ac2..4e70286 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -39,6 +39,7 @@ scripts = README.md \ show_response.py \ show_sensitivity.py \ show_spectrum.py \ + show_lightcurve.py \ show_ts_distribution.py \ show_obs.py diff --git a/examples/README.md b/examples/README.md index 5aabb5b..6dba7ba 100644 --- a/examples/README.md +++ b/examples/README.md @@ -70,5 +70,8 @@ the ctools. `show_spectrum.py` > Displays spectrum generated by the csspec script. +`show_lightcurve.py` +> Displays lightcurve generated by the cslightcrv script. + `show_ts_distribution.py` > Displays Test Statistic histograms generated using the cstsdist script. diff --git a/examples/show_lightcurve.py b/examples/show_lightcurve.py new file mode 100755 index 0000000..ea0788d --- /dev/null +++ b/examples/show_lightcurve.py @@ -0,0 +1,144 @@ +#! /usr/bin/env python +# ========================================================================== +# Display lightcurve generated by cslightcrv +# +# 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 . +# +# ========================================================================== +import sys +import gammalib +import cscripts +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except (ImportError, RuntimeError): + print('This script needs the "matplotlib" module') + sys.exit() + + +# =============== # +# Plot lightcurve # +# =============== # +def plot_lightcurve(filename, plotfile): + """ + Plot lightcurve + + Parameters + ---------- + filename : str + Name of lightcurve FITS file + plotfile : str + Plot file name + """ + # Read spectrum file + fits = gammalib.GFits(filename) + table = fits.table(1) + c_mjd = table['MJD'] + c_emjd = table['e_MJD'] + c_flux = table['EnergyFlux'] + c_eflux = table['e_EnergyFlux'] + c_ts = table['TS'] + c_upper = table['EFluxUpperLimit'] + + # Initialise arrays to be filled + mjd = [] + e_mjd = [] + flux = [] + e_flux = [] + ul_mjd = [] + ul_e_mjd = [] + ul_flux = [] + ul_e_flux = [] + + # Loop over rows of the file + nrows = table.nrows() + for row in range(nrows): + + # Get Test Statistic, flux and flux error + ts = c_ts.real(row) + flx = c_flux.real(row) + e_flx = c_eflux.real(row) + + # If Test Statistic is larger than 9 and twice the flux error is + # smaller than the flux, then append flux point ... + if ts > 9.0 and 2.0*e_flx < flx: + mjd.append(c_mjd.real(row)) + e_mjd.append(c_emjd.real(row)) + flux.append(c_flux.real(row)) + e_flux.append(c_eflux.real(row)) + + # ... otherwise append upper limit + else: + ul_mjd.append(c_mjd.real(row)) + ul_e_mjd.append(c_emjd.real(row)) + ul_flux.append(c_upper.real(row)) + ul_e_flux.append(0.5*c_upper.real(row)) + + # Plot the spectrum + plt.figure() + plt.semilogy() + plt.grid() + plt.errorbar(mjd, flux, yerr=e_flux, xerr=[e_mjd, e_mjd], + fmt='ro') + plt.errorbar(ul_mjd, ul_flux, xerr=[ul_e_mjd, ul_e_mjd], + yerr=ul_e_flux, uplims=True, fmt='ro') + plt.xlabel('MJD (days)') + plt.ylabel(r'E $\times$ dN/dE (erg cm$^{-2}$ s$^{-1}$)') + + # Show figure + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return + + +# =============== # +# Show lightcurve # +# =============== # +def show_lightcurve(): + """ + Show lightcurve + """ + # Set usage string + usage = 'show_lightcurve.py [-p plotfile] [file]' + + # Set default options + options = [{'option': '-p', 'value': ''}] + + # Get arguments and options from command line arguments + args, options = cscripts.ioutils.get_args_options(options, usage) + + # Extract script parameters from options + plotfile = options[0]['value'] + + # Plot lightcurve + plot_lightcurve(args[0], plotfile) + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Show lightcurve + show_lightcurve() diff --git a/examples/show_spectrum.py b/examples/show_spectrum.py index f7ff462..18f18be 100755 --- a/examples/show_spectrum.py +++ b/examples/show_spectrum.py @@ -35,7 +35,7 @@ except (ImportError, RuntimeError): # ============= # def plot_spectrum(filename, plotfile): """ - Plot sensitivity data + Plot spectrum Parameters ---------- diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index 4eddd9c..7738a0f 100644 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -248,6 +248,9 @@ void ctulimit::run(void) // Write header into logger log_header1(TERSE, "Compute best-fit likelihood"); + // Make sure that requested model parameters is free + m_model_par->free(); + // Optimize and save best log-likelihood m_obs.optimize(m_opt); m_obs.errors(m_opt); @@ -610,6 +613,9 @@ void ctulimit::ulimit_bisection(const double& min, const double& max) if (std::abs(eval_mid) < m_tol) { break; } + if (std::abs(wrk_max-wrk_min) < m_tol) { + break; + } // Change boundaries for further iteration if (eval_mid > 0.0) { diff --git a/test/data/lightcurve.fits b/test/data/lightcurve.fits new file mode 100644 index 0000000..a1ef7b2 --- /dev/null +++ b/test/data/lightcurve.fits @@ -0,0 +1,15 @@ +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 = 80 / width of table in bytes NAXIS2 = 50 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 10 / number of fields in each row TTYPE1 = 'MJD ' / label for field 1 TFORM1 = '1D ' / data format of field: 8-byte DOUBLE TUNIT1 = 'days ' / physical unit of field TTYPE2 = 'e_MJD ' / label for field 2 TFORM2 = '1D ' / data format of field: 8-byte DOUBLE TUNIT2 = 'days ' / physical unit of field TTYPE3 = 'EnergyFlux' / label for field 3 TFORM3 = '1D ' / data format of field: 8-byte DOUBLE TUNIT3 = 'erg/cm2/s' / physical unit of field TTYPE4 = 'e_EnergyFlux' / label for field 4 TFORM4 = '1D ' / data format of field: 8-byte DOUBLE TUNIT4 = 'erg/cm2/s' / physical unit of field TTYPE5 = 'Index ' / label for field 5 TFORM5 = '1D ' / data format of field: 8-byte DOUBLE TTYPE6 = 'e_Index ' / label for field 6 TFORM6 = '1D ' / data format of field: 8-byte DOUBLE TTYPE7 = 'TS ' / label for field 7 TFORM7 = '1D ' / data format of field: 8-byte DOUBLE TTYPE8 = 'DiffUpperLimit' / label for field 8 TFORM8 = '1D ' / data format of field: 8-byte DOUBLE TUNIT8 = 'ph/cm2/s/MeV' / physical unit of field TTYPE9 = 'FluxUpperLimit' / label for field 9 TFORM9 = '1D ' / data format of field: 8-byte DOUBLE TUNIT9 = 'ph/cm2/s' / physical unit of field TTYPE10 = 'EFluxUpperLimit' / label for field 10 TFORM10 = '1D ' / data format of field: 8-byte DOUBLE TUNIT10 = 'erg/cm2/s' / physical unit of field EXTNAME = 'LIGHTCURVE' / name of this extension END @+?+Nc;mc_=DB(yǿFFmP<00|a~}= &=j@+?+Nc=-=ݪP0__YS?2@>v<\9F3=j;h=Д {@+ X?+Nk=էv{="ӣsYI?@ow@*s<}bT[P>ޠb=Ӝ_@+\(?+N=4Fvs~=ŃȜI`W??jf@yeȠNn!=(Lxh@+-?+Nc=6d=ėsVPN]?LL?z8@wȢ:<~5S>GNX(=BvԾ)@+/c?+Nk=bmOb=%^ɪ`?\0@zA6<~w>RQJr=( + @+?+Nc=#sD =%.ncY?֐@Ke` m=:dN@+i?+Nc=\'=~LGB?t @~, ʸ=l(S~@+ m:?+N[=ᰆCo=sS&bV?b@~z!0<w}> +^y=n1N]@+# +>?+Nc=rNb=iA?Rh@tU BV@<{q>K@Bu=F@+'@t?+Nk=*/=H6,/?]]X@|c<:> +vشP=2'RM@+*?+N=3=5sl e?ʍ0f@ҠN, +)a=g s@+.z?+Nc=n)vDZ*=ȍ_|aZ?o@C <,8 > .=զ4@+1~K?+Nk=i$=`2L8LA?DnX6{@{FY<k> ihi=lA]@+4N?+Nc=y w=A~l*? s`ѧ@uX0<=>yL}=hvW=ѰhF@+;?+N[=xAOn=M&g$39?<2@sq  ֆ=@+?%?+Nc="=05NibW?O)#6@}+p< fM> +j;"=rt@+B\)?+Nk=f~=+ș/i?bh@wA}I=S?@+E,`?+N=G =%B?ąD+,@p``<p,;> ]RP=.@+Ib?+Nc=HlQ=,^},T?.=@n@oE<}v>w>=r=hqXw@+L?+Nk=Np&=)խM|h1ZulA?@nu +tN45Q?̪@Nt6-{y!=ݿ%&@+aG?+N=P=73S?č#@V+k= p|A?_m@^5ڀH=kjP@+|/?+Nc=)=MdcQ?)K@w9"y+ 4= %Lط@+?+N=iǰ=FI > +3?Z@r6<_m >i=)t@+i7?+Nk='Ccs=λDt2?Ɲ*(@aT@4΋{<1G0b=N.=J5 @+?+Nc;^9 ^Ab\w*_@E%c-0C`r<3ô"Z= ~=BAN@+zH?+Nk;^9 ^AX@ENBOᦿC.<=Nbd=idb=*;U@+~?+Nc;^9 ^AS0O+@Ese<7ɾ=qG=c̪o@+N?+N;^9 ^=ܻ)UBe;ܼ=<1;h-= M1=ũ4u@+Q?+Nk;^9 ^= @8T-eBLEI.@<@UͰ8o=6=+<@+"""?+Nc;^9 ^==wBVbus.<<8@+=<(e=|]@+Y?+N[;^9 ^A%mf @F=yzuA<92x=7!=-C@+?+Nc;^9 ^By"I@FOjH<-]=Lu=9pX \ No newline at end of file diff --git a/test/test_cslightcrv.py b/test/test_cslightcrv.py index a5bcd71..2eaf4c9 100755 --- a/test/test_cslightcrv.py +++ b/test/test_cslightcrv.py @@ -2,7 +2,7 @@ # ========================================================================== # This scripts performs unit tests for the cslightcrv script. # -# 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 @@ -239,7 +239,8 @@ class Test(test): """ # Expected column names cols = ['MJD', 'e_MJD', 'Prefactor', 'e_Prefactor', - 'Index', 'e_Index', 'TS', 'UpperLimit'] + 'Index', 'e_Index', 'TS', 'DiffUpperLimit', + 'FluxUpperLimit', 'EFluxUpperLimit'] # Open FITS file fits = gammalib.GFits(filename) diff --git a/test/test_examples.py.in b/test/test_examples.py.in index 4c9cb27..6c6cab8 100755 --- a/test/test_examples.py.in +++ b/test/test_examples.py.in @@ -2,7 +2,7 @@ # ========================================================================== # Performs unit tests for the ctools example scripts # -# 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 @@ -128,6 +128,8 @@ class Test(gammalib.GPythonTestSuite): 'Test show_sensitivity.py') self.append(self.test_show_spectrum, 'Test show_spectrum.py') + self.append(self.test_show_lightcurve, + 'Test show_lightcurve.py') self.append(self.test_show_ts_distribution, 'Test ts_distribution.py') self.append(self.test_show_obs, @@ -401,6 +403,20 @@ class Test(gammalib.GPythonTestSuite): # Return return + # Test show_lightcurve + def test_show_lightcurve(self): + """ + Test show_lightcurve + """ + # Set arguments + args = '-p example_show_lightcurve.eps %s/lightcurve.fits' % self._datadir + + # Execute script + self._execute_python('show_lightcurve', args=args) + + # Return + return + # Test show_ts_distribution def test_show_ts_distribution(self): """