Add show_phases.py example script and unit test

Jürgen Knödlseder
1 parent 96d2ba22
Showing 5 changed files with 158 additions and 11 deletions
examples/Makefile.am
examples/README.md
examples/show_phases.py
test/data/obs_unbinned.xml
test/test_examples.py.in
@@ -41,6 +41,7 @@ scripts = README.md \
           show_sensitivity.py \
           show_spectrum.py \
           show_lightcurve.py \
+          show_phases.py \
           show_ts_distribution.py \
           show_obs.py \
           show_pointings.py
@@ -4,7 +4,7 @@ ctools examples
 This folder contains Python example scripts that illustrate the usage of
 the ctools.
  
-`generate_prod3_irfs/py`
+`generate_prod3_irfs.py`
 >  Generates Prod3 Instrument Response Functions from FITS files.
  
 `make_pointings.py`
@@ -48,12 +48,27 @@ the ctools.
 `show_butterfly.py`
 >  Displays a butterfly generated using the ctbutterfly tool.
  
+`show_irf.py`
+>  Displays an Instrument Response Function.
+
+`show_lightcurve.py`
+>  Displays lightcurve generated by the cslightcrv script.
+
 `show_model.py`
 >  Displays a spectral model.
  
+`show_obs.py`
+>  Displays an observation definition XML file.
+
 `show_pha.py`
 >  Displays a spectrum in PHA format.
  
+`show_phases.py`
+>  Displays event phases in an event list of observation definition XML file.
+
+`show_pointings.py`
+>  Displays pointings in an observation definition XML file.
+
 `show_pull_evolution.py`
 >  Displays the evolution of the mean and rms pull as function of number of
 >  trials.
@@ -70,8 +85,5 @@ 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 phase information from event list
+#
+# 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 phases #
+# =========== #
+def plot_phases(filename, plotfile, nphases=40):
+    """
+    Plot phases
+
+    Parameters
+    ----------
+    filename : str
+        Name of lightcurve FITS file
+    plotfile : str
+        Plot file name
+    """
+    # Read observation container. If an exception occurs then try loading the
+    # file as an event list and build an observation container.
+    try:
+        obs = gammalib.GObservations(filename)
+    except:
+        obs = gammalib.GObservations()
+        obs.append(gammalib.GCTAObservation(filename))
+
+    # Initialise phases
+    phases = []
+
+    # Loop over all observations
+    for run in obs:
+
+        # Get events
+        events = run.events()
+
+        # Loop over all events
+        for event in events:
+
+            # Get phase
+            phase = event.phase()
+
+            # Collect phase
+            phases.append(phase)
+            phases.append(phase+1.0)
+
+    # Plot phase
+    plt.figure()
+    plt.hist(phases,bins=nphases,range=(0.0,2.0),facecolor='r')
+    plt.xlabel('Phase')
+    plt.ylabel('Events')
+
+    # Show figure
+    if len(plotfile) > 0:
+        plt.savefig(plotfile)
+    else:
+        plt.show()
+
+    # Return
+    return
+
+
+# =========== #
+# Show phases #
+# =========== #
+def show_phases():
+    """
+    Show phases
+    """
+    # Set usage string
+    usage = 'show_phases.py [-p plotfile] [file]'
+
+    # Set default options
+    options = [{'option': '-n', 'value': '20'},
+               {'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
+    nphases  = int(options[0]['value'])*2
+    plotfile = options[1]['value']
+
+    # Plot phases
+    plot_phases(args[0], plotfile, nphases=nphases)
+
+    # Return
+    return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+    # Show phases
+    show_phases()
 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
 <observation_list title="observation list">
   <observation name="" id="00001" instrument="CTA">
-    <parameter name="EventList"           file="crab_events.fits" />
+    <parameter name="EventList"           file="$(TEST_DATA)/crab_events.fits" />
     <parameter name="EffectiveArea"       file="$(TEST_DATA)/irf_file.fits" />
     <parameter name="PointSpreadFunction" file="$(TEST_DATA)/irf_file.fits" />
     <parameter name="EnergyDispersion"    file="$(TEST_DATA)/irf_file.fits" />
@@ -120,6 +120,8 @@ class Test(gammalib.GPythonTestSuite):
                     'Test show_model.py')
         self.append(self.test_show_pha,
                     'Test show_pha.py')
+        self.append(self.test_show_phases,
+                    'Test show_phases.py')
         self.append(self.test_show_pull_evolution,
                     'Test show_pull_evolution.py')
         self.append(self.test_show_pull_histogram,
@@ -334,6 +336,20 @@ class Test(gammalib.GPythonTestSuite):
         # Return
         return
  
+    # Test show_phases
+    def test_show_phases(self):
+        """
+        Test show_phases
+        """
+        # Set arguments
+        args = '-p example_show_phases.eps %s/phased_events.fits' % self._datadir
+
+        # Execute script
+        self._execute_python('show_phases', args=args)
+
+        # Return
+        return
+
     # Test show_pull_evolution
     def test_show_pull_evolution(self):
         """
@@ -526,12 +542,6 @@ if __name__ == &#39;__main__&#39;:
     # Save test results
     suites.save('reports/examples.xml')
  
-    # Set return code
-    if success:
-        rc = 0
-    else:
-        rc = 1
-
     # Raise an exception in case of failure
     if not success:
         raise RuntimeError('At least one error occured during the test.')