Add show_lightcurve example script

Jürgen Knödlseder
1 parent be2350c0
Showing 10 changed files with 262 additions and 36 deletions
NEWS
cscripts/Makefile.am
cscripts/cslightcrv.py
examples/Makefile.am
examples/README.md
examples/show_lightcurve.py
examples/show_spectrum.py
test/data/lightcurve.fits
test/test_cslightcrv.py
test/test_examples.py.in
 New Features and Important Changes in ctools 1.2.0
-1 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
 ----------------------------------------------
@@ -232,11 +235,14 @@ definition XML file (#1881).
 Add cssrcdetect script to detect sources in a sky map (#1617).
+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.spectral()[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.spectral()[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'])
@@ -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
     ----------
@@ -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):
         """