diff --git a/.gitignore b/.gitignore index 15c3b77..ad3dd70 100644 --- a/.gitignore +++ b/.gitignore @@ -82,6 +82,7 @@ test/*.reg test/*.trs test/*.txt test/*.eps +test/*.jpg test/test_cscripts.sh test/test_examples.py test/reports/*.xml diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 8b69c45..3f033fe 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -53,6 +53,7 @@ class cssens(ctools.cscript): # Initialise class members self._obs = gammalib.GObservations() self._ebounds = gammalib.GEbounds() + self._obs_ebounds = [] self._srcname = "" self._outfile = gammalib.GFilename() self._ra = None @@ -270,8 +271,41 @@ class cssens(ctools.cscript): # Loop over all observations in container for obs in self._obs: + # Get observation energy boundaries + obs_ebounds = obs.events().ebounds() + + # Get minimum and maximum energy + obs_emin = obs_ebounds.emin() + obs_emax = obs_ebounds.emax() + + # Case A: bin fully contained in observation ebounds + if obs_emin <= emin and obs_emax >= emax: + ebounds = gammalib.GEbounds(emin, emax) + + # Case B: bin fully outside of obs ebounds + elif emax < obs_emin or emin > obs_emax: + + # Set zero range (inspired by ctselect) + e0 = gammalib.GEnergy(0.0, "TeV") + ebounds = gammalib.GEbounds(e0, e0) + + # Case C: bin partly overlapping with observation ebounds + else: + + # Set energy range as obs ebounds were fully contained inside energy bin + set_emin = emin + set_emax = emax + + # Adjust energy bin to respect observation energy boundary + if emin < obs_emin: + set_emin = obs_emin + if emax > obs_emax: + set_emax = obs_emax + + #Set energy boundaries + ebounds = gammalib.GEbounds(set_emin, set_emax) + # Set energy boundaries - ebounds = gammalib.GEbounds(emin, emax) obs.events().ebounds(ebounds) # Return @@ -624,6 +658,10 @@ class cssens(ctools.cscript): for par in model: self._log(str(par)+"\n") + # Restore energy boundaries of observation container + for i, obs in enumerate(self._obs): + obs.events().ebounds(self._obs_ebounds[i]) + # Store result result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ 'crab_flux': crab_flux, 'photon_flux': photon_flux, \ @@ -645,7 +683,11 @@ class cssens(ctools.cscript): # Get parameters self._get_parameters() - + + # Loop over observations and store ebounds + for obs in self._obs: + self._obs_ebounds.append(obs.events().ebounds()) + # Initialise script colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', 'energy_flux', 'sensitivity'] diff --git a/examples/show_obs.py b/examples/show_obs.py new file mode 100755 index 0000000..3dd46a0 --- /dev/null +++ b/examples/show_obs.py @@ -0,0 +1,200 @@ +#! /usr/bin/env python +# ========================================================================== +# Display summary of observation definition XML file +# +# Copyright (C) 2016 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 +# 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 ctools +import cscripts +try: + import matplotlib.pyplot as plt + plt.figure() + plt.close() +except: + print('This script needs the "matplotlib" module') + sys.exit() + + +# ==================== # +# Run csobsinfo script # +# ==================== # +def run_csobsinfo(filename, ra=None, dec=None, debug=True): + """ + Run csobsinfo script + + Parameters + ---------- + filename : str + File name of observation definition XML file + ra : float, optional + Target Right Ascension (deg) + dec : float, optional + Target declination of pointing (deg) + debug : bool, optional + Switch on debugging in csobsinfo run + + Returns + ------- + info : `~cscripts.csobsinfo` + csobsinfo instance + """ + # Get observation definition XML filename + obsdef = gammalib.GFilename(filename) + + # Setup csobsinfo script + info = cscripts.csobsinfo() + info['inobs'] = obsdef.url() + info['debug'] = debug + + # Set offset + if ra != None and dec != None: + info['ra'] = ra + info['dec'] = dec + info['offset'] = True + else: + info['offset'] = False + + # Run csobsinfo + info.run() + + # Return + return info + + +# ================ # +# Plot information # +# ================ # +def plot_information(info, plotfile=''): + """ + Plot information + + Parameters + ---------- + info : `~cscripts.csobsinfo` + csobsinfo instance + plotfile : str, optional + Plot filename + """ + # Retrieve observation info + zeniths = info.zeniths() + azimuths = info.azimuths() + offsets = info.offsets() + ebounds = info.ebounds() + gti = info.gti() + + # Create figure with subplots + nrows = 2 + ncols = 2 + if info['offset'].boolean() and ebounds.size() > 0: + ncols += 1 + iplot = 1 + plt.figure() + plt.subplot(nrows, ncols, iplot) + + # Plot zenith angle distribution + zmin = min(zeniths) + zmax = max(zeniths) + plt.hist(zeniths, bins=30, range=(zmin, zmax), fc='blue') + plt.xlabel('Zenith Angle (deg)') + plt.ylabel('Frequency') + plt.title('Zenith angle distribution') + + # Plot azimuth angle distribution + iplot += 1 + plt.subplot(nrows, ncols, iplot) + amin = min(azimuths) + amax = max(azimuths) + plt.hist(azimuths, bins=30, range=(amin, amax), fc='blue') + plt.xlabel('Azimuth Angle (deg)') + plt.ylabel('Frequency') + plt.title('Azimuth distribution') + + # Plot offset if possible + if info['offset'].boolean(): + iplot += 1 + plt.subplot(nrows, ncols, iplot) + omin = min(offsets) + omax = max(offsets) + plt.hist(offsets, bins=30, range=(omin, omax), fc='blue') + plt.xlabel('Offset from (RA, DEC)=('+str(ra)+','+str(dec)+') (deg)') + plt.ylabel('Frequency') + plt.title('Offset distribution') + + # Plot energy thresholds if possible + if ebounds.size() > 0: + emin = [] + emax = [] + for i in range(ebounds.size()): + emin.append(ebounds.emin(i).log10TeV()) + emax.append(ebounds.emax(i).log10TeV()) + iplot += 1 + plt.subplot(nrows, ncols, iplot) + plt.hist(emin, bins=80, range=(-1.0, 2.0), fc='red', label='emin') + plt.hist(emax, bins=80, range=(-1.0, 2.0), fc='blue', label='emax') + plt.xlabel('Energy threshold (log10 (E/TeV))') + plt.ylabel('Frequency') + plt.legend(loc='upper left') + plt.title('Energy threshold') + + # Plot observation point in time wrt zenith angle + iplot += 1 + plt.subplot(nrows, ncols, iplot) + times = [] + for i in range(gti.size()): + tmean = gti.tstart(i) + 0.5*(gti.tstart(i)-gti.tstart(i)) + times.append(tmean.mjd()) + plt.plot(times, zeniths, 'o', lw=2.0, color='black') + plt.xlabel('Time (MJD)') + plt.ylabel('Zenith Angle (deg)') + plt.title('Observation time') + + # Show plots or save it into file + if len(plotfile) > 0: + plt.savefig(plotfile) + else: + plt.show() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Print usage message if wrong number of arguments is provided + if len(sys.argv) < 2: + sys.exit('Usage: inspect_obs.py inobs.xml [ra] [dec] [file]') + + # Extract optional parameters + plotfile = '' + ra = None + dec = None + if len(sys.argv) == 3: + plotfile = sys.argv[2] + elif len(sys.argv) == 4: + ra = float(sys.argv[2]) + dec = float(sys.argv[3]) + + # Run csobsinfo + info = run_csobsinfo(sys.argv[1], ra=ra, dec=dec) + + # Plot information + plot_information(info, plotfile=plotfile) diff --git a/test/test_examples.py.in b/test/test_examples.py.in index 54aa4bf..402ce30 100755 --- a/test/test_examples.py.in +++ b/test/test_examples.py.in @@ -108,6 +108,7 @@ class Test(gammalib.GPythonTestSuite): self.append(self.test_show_sensitivity, 'Test test_show_sensitivity.py') self.append(self.test_show_spectrum, 'Test test_show_spectrum.py') self.append(self.test_show_ts_distribution, 'Test show_ts_distribution.py') + self.append(self.test_show_obs, 'Test show_obs.py') # Return return @@ -128,8 +129,6 @@ class Test(gammalib.GPythonTestSuite): self._execute_python(script, args='gc', log='make_pointings_gc') self._execute_python(script, args='lmc', log='make_pointings_lmc') - #TODO: Do any testing - # Return return @@ -141,8 +140,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('generate_prod3_irfs') - #TODO: Do any testing - # Return return @@ -154,8 +151,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('make_spectrum', args='eps') - #TODO: Do any testing - # Return return @@ -167,8 +162,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('make_ts_distributions', args='-n 2 -e 0 -d 1800') - #TODO: Do any testing - # Return return @@ -180,8 +173,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('pipeline_binned_disk') - #TODO: Do any testing - # Return return @@ -193,8 +184,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('pipeline_binned_mem') - #TODO: Do any testing - # Return return @@ -206,8 +195,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('pipeline_stacked_disk') - #TODO: Do any testing - # Return return @@ -219,8 +206,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('pipeline_stacked_mem') - #TODO: Do any testing - # Return return @@ -232,8 +217,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('pipeline_unbinned_disk') - #TODO: Do any testing - # Return return @@ -245,8 +228,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('pipeline_unbinned_mem') - #TODO: Do any testing - # Return return @@ -260,8 +241,6 @@ class Test(gammalib.GPythonTestSuite): args='data/butterfly.txt ' 'example_show_butterfly.eps') - #TODO: Do any testing - # Return return @@ -274,8 +253,6 @@ class Test(gammalib.GPythonTestSuite): self._execute_python('show_model', args='data/crab.xml Crab example_show_model.eps') - #TODO: Do any testing - # Return return @@ -288,8 +265,6 @@ class Test(gammalib.GPythonTestSuite): self._execute_python('show_pha', args='data/pha.fits example_show_pha.eps') - #TODO: Do any testing - # Return return @@ -303,8 +278,6 @@ class Test(gammalib.GPythonTestSuite): args='data/pull.dat Pull_Crab_Prefactor ' 'example_show_pull_evolution.eps') - #TODO: Do any testing - # Return return @@ -318,8 +291,6 @@ class Test(gammalib.GPythonTestSuite): args='data/pull.dat Pull_Crab_Prefactor 50 ' 'example_show_pull_histogram.eps') - #TODO: Do any testing - # Return return @@ -331,8 +302,6 @@ class Test(gammalib.GPythonTestSuite): # Execute script self._execute_python('show_response', args='example_show_response.eps') - #TODO: Do any testing - # Return return @@ -346,8 +315,6 @@ class Test(gammalib.GPythonTestSuite): args='data/sensitivity.dat ' 'example_show_sensitivity.eps') - #TODO: Do any testing - # Return return @@ -361,8 +328,6 @@ class Test(gammalib.GPythonTestSuite): args='data/spectrum.fits ' 'example_show_spectrum.eps') - #TODO: Do any testing - # Return return @@ -376,7 +341,17 @@ class Test(gammalib.GPythonTestSuite): args='data/ts.dat -n 100 ' '-f example_show_ts_distribution.eps') - #TODO: Do any testing + # Return + return + + # Test show_obs + def test_show_obs(self): + """ + Test show_obs + """ + # Execute script + self._execute_python('show_obs', + args='obs_selected.xml example_show_obs.jpg') # Return return