Update example scripts

Jürgen Knödlseder
1 parent bbe311d4
Showing 17 changed files with 126 additions and 123 deletions
examples/benchmark_cta_analysis.py
examples/make_pointings.py
examples/make_pull_at_sensitivity_limit.py
examples/make_spectrum.py
examples/make_survey.py
examples/make_ts_distributions.py
examples/pipeline_binned_disk.py
examples/show_butterfly.py
examples/show_model.py
examples/show_pha.py
examples/show_pull_evolution.py
examples/show_pull_histogram.py
examples/show_response.py
examples/show_sensitivity.py
examples/show_spectrum.py
examples/show_ts_distribution.py
examples/simulate_events.py
@@ -299,8 +299,7 @@ def stacked_pipeline(model_name, duration):
 # Main routine entry point #
 #==========================#
 if __name__ == '__main__':
-    """
-    """
+
     # Initialise parameters
     title      = "CTA analysis benchmark"
     model_name = "${CTOOLS}/share/models/crab.xml"
@@ -336,9 +336,7 @@ def write_obsdef(filename, obsdef):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    CTA pointing pattern generation.
-    """
+
     # Initialise flags
     need_help = False
  
@@ -135,11 +135,7 @@ def create_model(flux, index=-2.48, fitidx=False):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Create pull distributions for the energy bands and fluxes given in the
-    file specified in the argument. The file needs to be created using the
-    cssens script.
-    """
+
     # Get input arguments
     usage = """
     make_pull_at_sensitivity_limit filename [max_threads]
@@ -167,9 +167,7 @@ def plot_spectrum(spectrum):
 # Main routine entry point #
 #==========================#
 if __name__ == '__main__':
-    """
-    Generate a Crab spectrum.
-    """
+
     # Generate spectrum
     spectrum = make_spectrum()
  
@@ -303,9 +303,7 @@ def survey_gplane(lrange=10, lstep=2):
 # Main routine entry point #
 #==========================#
 if __name__ == '__main__':
-    """
-    CTA survey simulation and analysis script.
-    """
+
     # Initialise flags
     need_help = False
  
@@ -94,9 +94,7 @@ def create_ts(loge, emin, emax, ntrials=100, duration=180000.0,
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Create TS distribution in a number of energy bands.
-    """
+
     # Get input arguments
     usage = "make_ts_distributions [-n ntrials] [-e enumbins] [-m max_threads]"
     if len(sys.argv) < 1:
 #! /usr/bin/env python
 # ==========================================================================
-# This script illustrates how to perform a binned CTA analysis based on
-# simulated CTA data. You may use and adapt this script to implement your
-# own pipeline.
+# Perform binned CTA analysis based of simulated CTA data.
 #
-# Usage:
-#   ./pipeline_binned_disk.py
-#
-# ==========================================================================
-#
-# Copyright (C) 2015 Juergen Knoedlseder
+# Copyright (C) 2015-2016 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
@@ -33,25 +26,41 @@ from cscripts import obsutils
 # ================== #
 # Setup observations #
 # ================== #
-def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5,
+def setup_observations(pattern='four', ra=83.63, dec=22.01, offset=1.5,
                        emin=0.1, emax=100.0, rad=5.0, duration=1800.0,
-                       deadc=0.95,
-                       caldb="prod2", irf="South_50h"):
+                       deadc=0.95, caldb='prod2', irf='South_0.5h'):
     """
     Returns an observation container.
  
-    Keywords:
-     pattern   - Pointing pattern, either "single" or "four"
-     ra        - RA of pattern centre [deg] (default: 83.6331)
-     dec       - DEC of pattern centre [deg] (default: 22.0145)
-     offset    - Offset between observations of pattern [deg] (default: 1.5)
-     emin      - Minimum energy [TeV] (default: 0.1)
-     emax      - Maximum energy [TeV] (default: 100.0)
-     rad       - ROI radius used for analysis [deg] (default: 5.0)
-     duration  - Duration of one CTA observation [seconds] (default: 1800.0)
-     deadc     - Deadtime correction factor (default: 0.95)
-     caldb     - Calibration database path (default: "dummy")
-     irf       - Instrument response function (default: cta_dummy_irf)
+    Parameters
+    ----------
+    pattern : str, optional
+        Pointing pattern, either 'single' or 'four'
+    ra : float, optional
+        Right Ascension of pattern centre [deg]
+    dec : float, optional
+        Delication of pattern centre [deg]
+    offset : float, optional
+        Offset between observations of pattern [deg]
+    emin : float, optional
+        Minimum energy [TeV]
+    emax : float, optional
+        Maximum energy [TeV]
+    rad : float, optional
+        ROI radius used for analysis [deg]
+    duration : float, optional
+        Duration of one CTA observation [seconds]
+    deadc : float, optional
+        Deadtime correction factor
+    caldb : str, optional
+        Calibration database path
+    irf : str, optional
+        Instrument response function
+
+    Returns
+    -------
+    obs : `~gammalib.GObservations`
+        Observation container
     """
     # Set list of observations
     obs_def_list = obsutils.set_obs_patterns(pattern,
@@ -75,12 +84,21 @@ def setup_observations(pattern=&quot;four&quot;, ra=83.63, dec=22.01, offset=1.5,
 # =========== #
 # Setup model #
 # =========== #
-def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
+def setup_model(obs, model='${CTOOLS}/share/models/crab.xml'):
     """
     Setup model for analysis.
  
-    Keywords:
-     model - Model Xml file
+    Parameters
+    ----------
+    obs : `~gammalib.GObservations`
+        Observation container
+    model : str, optional
+        Model definition XML file
+
+    Returns
+    -------
+    obs : `~gammalib.GObservations`
+        Observation container
     """
     # Append model from file to observation container
     obs.models(gammalib.GModels(model))
@@ -94,56 +112,71 @@ def setup_model(obs, model=&quot;${CTOOLS}/share/models/crab.xml&quot;):
 # ================================ #
 def run_pipeline(obs, emin=0.1, emax=100.0,
                  enumbins=20, nxpix=200, nypix=200, binsz=0.02,
-                 coordsys="CEL", proj="CAR",
-                 model="${CTOOLS}/share/models/crab.xml",
-                 caldb="prod2", irf="South_50h",
+                 coordsys='CEL', proj='CAR',
+                 model='${CTOOLS}/share/models/crab.xml',
+                 caldb='prod2', irf='South_0.5h',
                  debug=False):
     """
     Simulation and binned analysis pipeline.
  
-    Keywords:
-     emin     - Minimum energy of cube [TeV] (default: 0.1)
-     emax     - Maximum energy of cube [TeV] (default: 100.0)
-     enumbins - Number of energy bins in cube (default: 20)
-     nxpix    - Number of RA pixels in cube (default: 200)
-     nypix    - Number of DEC pixels in cube (default: 200)
-     binsz    - Spatial cube bin size [deg] (default: 0.02)
-     coordsys - Cube coordinate system (CEL or GAL)
-     proj     - Cube World Coordinate System (WCS) projection
-     model    - Model Xml file
-     caldb    - Calibration database path (default: "dummy")
-     irf      - Instrument response function (default: cta_dummy_irf)
-     debug    - Enable debugging (default: False)
+    Parameters
+    ----------
+    obs : `~gammalib.GObservations`
+        Observation container
+    emin : float, optional
+        Minimum energy [TeV]
+    emax : float, optional
+        Maximum energy [TeV]
+    enumbins : int, optional
+        Number of energy bins
+    nxpix : int, optional
+        Number of pixels in X axis
+    nypix : int, optional
+        Number of pixels in Y axis
+    binsz : float, optional
+        Pixel size [deg/pixel]
+    coordsys : str, optional
+        Coordinate system
+    proj : str, optional
+        Coordinate projection
+    model : str, optional
+        Model definition XML file
+    caldb : str, optional
+        Calibration database path
+    irf : str, optional
+        Instrument response function
+    debug : bool, optional
+        Debug function
     """
     # Simulate events
     sim = ctools.ctobssim(obs)
-    sim["debug"] = debug
-    sim["outevents"] = "obs.xml"
+    sim['debug']     = debug
+    sim['outevents'] = 'obs.xml'
     sim.execute()
  
     # Bin events by looping over all observations in the container
-    sim_obs = gammalib.GObservations("obs.xml")
+    sim_obs = gammalib.GObservations('obs.xml')
     obs     = gammalib.GObservations()
     for run in sim_obs:
  
         # Get event filename and set counts cube filename
-        eventfile = run.eventfile()
-        cubefile  = "cube_"+eventfile
+        eventfile = run.eventfile().url()
+        cubefile  = 'cube_'+eventfile
  
         # Bin events for that observation
         bin = ctools.ctbin()
-        bin["inobs"]    = eventfile
-        bin["outcube"]  = cubefile
-        bin["ebinalg"]  = "LOG"
-        bin["emin"]     = emin
-        bin["emax"]     = emax
-        bin["enumbins"] = enumbins
-        bin["nxpix"]    = nxpix
-        bin["nypix"]    = nypix
-        bin["binsz"]    = binsz
-        bin["coordsys"] = coordsys
-        bin["usepnt"]   = True
-        bin["proj"]     = proj
+        bin['inobs']    = eventfile
+        bin['outcube']  = cubefile
+        bin['ebinalg']  = 'LOG'
+        bin['emin']     = emin
+        bin['emax']     = emax
+        bin['enumbins'] = enumbins
+        bin['nxpix']    = nxpix
+        bin['nypix']    = nypix
+        bin['binsz']    = binsz
+        bin['coordsys'] = coordsys
+        bin['usepnt']   = True
+        bin['proj']     = proj
         bin.execute()
  
         # Set observation ID
@@ -156,19 +189,19 @@ def run_pipeline(obs, emin=0.1, emax=100.0,
     # Save XML file
     xml = gammalib.GXml()
     obs.write(xml)
-    xml.save("obs_cube.xml")
+    xml.save('obs_cube.xml')
  
     # Perform maximum likelihood fitting
     like = ctools.ctlike()
-    like["inobs"]    = "obs_cube.xml"
-    like["inmodel"]  = model
-    like["outmodel"] = "fit_results.xml"
-    like["expcube"]  = "NONE"
-    like["psfcube"]  = "NONE"
-    like["bkgcube"]  = "NONE"
-    like["caldb"]    = caldb
-    like["irf"]      = irf
-    like["debug"]    = True # Switch this always on for results in console
+    like['inobs']    = 'obs_cube.xml'
+    like['inmodel']  = model
+    like['outmodel'] = 'fit_results.xml'
+    like['expcube']  = 'NONE'
+    like['psfcube']  = 'NONE'
+    like['bkgcube']  = 'NONE'
+    like['caldb']    = caldb
+    like['irf']      = irf
+    like['debug']    = True # Switch this always on for results in console
     like.execute()
  
     # Return
@@ -181,9 +214,9 @@ def run_pipeline(obs, emin=0.1, emax=100.0,
 if __name__ == '__main__':
  
     # Dump header
-    print("********************************************")
-    print("*       CTA binned analysis pipeline       *")
-    print("********************************************")
+    print('********************************************')
+    print('*       CTA binned analysis pipeline       *')
+    print('********************************************')
  
     # Setup observations
     obs = setup_observations()
@@ -30,8 +30,7 @@ except:
 # Script entry #
 # ============ #    
 if __name__ == "__main__":
-    """
-    """
+
     # Check for given butterfly file
     if not len(sys.argv) == 2:
         sys.exit("Usage: show_butterfly.py butterfly.txt")
@@ -82,8 +82,7 @@ def plot_spectrum(model, emin=0.01, emax=100.0, enumbins=100):
 # Script entry  #
 # ============= #    
 if __name__ == "__main__":
-    """
-    """
+
     # Check for model filename
     if len(sys.argv) < 2:
         sys.exit("Usage: show_model.py model.xml [index]")
@@ -79,9 +79,7 @@ def show_pha(pha):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Display spectrum in PHA format.
-    """
+
     # Print usage information
     usage = "Usage: show_pha filename" \
             " [-n bins] [-c column] [-t title] [-p plot]"
@@ -81,9 +81,7 @@ def read_pull(filename, parname):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Display pull evolution generated by cspull using matplotlib.
-    """
+
     # Print usage information
     usage = "Usage: show_pull_evolution filename parname"
     if len(sys.argv) < 3:
@@ -81,9 +81,7 @@ def read_pull(filename, parname):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Display pull histogram generated by cspull using matplotlib.
-    """
+
     # Print usage information
     usage = "Usage: show_pull_histogram filename parname [bins]"
     if len(sys.argv) < 3:
@@ -252,9 +252,7 @@ def show_one_response(rspname, dbname, name, rootdir=None, color=&quot;r&quot;):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Display response information.
-    """
+
     # Create figures
     plt.figure(1)
     plt.title("Effective area")
@@ -119,8 +119,7 @@ def show_sensitivity(sensitivity, filename):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    """
+
     # Print usage information
     usage = "Usage: show_sensitivity [filename]"
     if len(sys.argv) < 1:
@@ -31,8 +31,6 @@ except:
 # Script entry #
 # ============ #    
 if __name__ == "__main__":
-    """
-    """
  
     # Check for given spectrum file
     if not len(sys.argv) == 2:
@@ -219,9 +219,7 @@ def dist_pdf(values, nbins, title):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Display TS distribution generated using cstsdist using matplotlib.
-    """
+
     # Print usage information
     usage = "Usage: show_ts_distribution filename" \
             " [-n bins] [-c column] [-t title] [-p plot]"
@@ -239,17 +239,15 @@ def show_events(events, xmlname, duration, emin, emax, ebins=30):
 # Main routine entry point #
 #==========================#
 if __name__ == '__main__':
-    """
-    Simulate events.
-    """
+
     # Dump header
-    print("")
-    print("*******************")
-    print("* Simulate events *")
-    print("*******************")
+    print('')
+    print('*******************')
+    print('* Simulate events *')
+    print('*******************')
  
     # Set XML names
-    xmlname = "model_bgd.xml"
+    xmlname = 'model_bgd.xml'
  
     # Set simulation parameters
     duration = 360000.0 # Observation duration (s)