diff --git a/NEWS b/NEWS index 5a73fd4..d403b55 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.2.0 -1 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 ---------------------------------------------- @@ -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 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 ac44d82..6a6b160 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.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']) 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/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): """