Merge branch 'integration' into 1920-fix-cslightcrv

# Conflicts: # NEWS

Merge branch 'integration' into 1920-fix-cslightcrv
# Conflicts: # NEWS
Jürgen Knödlseder
2 parents a2090093 fe07c6cf
Showing 13 changed files with 517 additions and 69 deletions
NEWS
cscripts/Makefile.am
cscripts/cslightcrv.py
cscripts/csroot2caldb.py
dev/macosx-package.sh
examples/Makefile.am
examples/README.md
examples/show_lightcurve.py
examples/show_spectrum.py
src/ctulimit/ctulimit.cpp
test/data/lightcurve.fits
test/test_cslightcrv.py
test/test_examples.py.in
 New Features and Important Changes in ctools 1.2.0
-12 December 2016
+26 January 2017
 Introduction
@@ -195,6 +195,9 @@ in particular the &quot;PhotonFlux&quot; and &quot;EnergyFlux&quot; 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
@@ -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 \
@@ -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'])
@@ -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/<tel>/<inst>/bcf/<obsid>
@@ -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.
+#!/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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+
+# ============================================= #
+# 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 \
+<license file="COPYING"/>\
+<title>'$CTOOLS'</title>\
+<organization>irap.omp.eu</organization>\
+<description>CTA Science Analysis Tools</description>\
+' "$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 <<EOM >$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
@@ -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
@@ -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.
+#! /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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+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()
@@ -35,7 +35,7 @@ except (ImportError, RuntimeError):
 # ============= #
 def plot_spectrum(filename, plotfile):
     """
-    Plot sensitivity data
+    Plot spectrum
     Parameters
     ----------
@@ -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&amp; min, const double&amp; 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) {
@@ -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)
@@ -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):
         """