Update example scripts.

Jürgen Knödlseder
1 parent aa19272d
Showing 6 changed files with 36 additions and 1026 deletions
examples/Makefile.am
examples/README
examples/benchmark_cta_analysis.py
examples/make_binned_analysis.py
examples/make_cube_analysis.py
examples/make_unbinned_analysis.py
@@ -2,14 +2,21 @@
  
 # Scripts to be distributed
 scripts = README \
-          pipeline_stacked_mem.py \
-          make_binned_analysis.py \
+          benchmark_cta_analysis.py \
           make_pull_at_sensitivity_limit.py \
+          make_spectrum.py \
           make_ts_distributions.py \
-          make_unbinned_analysis.py \
+          pipeline_binned_disk.py \
+          pipeline_binned_mem.py \
+          pipeline_stacked_disk.py \
+          pipeline_stacked_mem.py \
+          pipeline_unbinned_disk.py \
+          pipeline_unbinned_mem.py \
+          show_butterfly.py \
+          show_pha.py \
+          show_response.py \
           show_pull_histogram.py \
-          show_ts_distribution.py \
-          show_butterfly.py
+          show_ts_distribution.py
  
 # Set directories
 scriptsdir   = $(datadir)/examples/python
@@ -2,6 +2,28 @@ This folder contains Python example scripts that illustrate the
 usage of the ctools.
  
 benchmark_cta_analysis.py
+  This script peforms a benchmark of the various CTA analysis types
+  (unbinned, binned, stacked), and if matplotlib is installed, creates
+  a plot of the benchmark results.
+
+make_pull_at_sensitivity_limit.py
+  This script creates pull distributions as function of energy
+  for a test source that is at the 5 sigma sensitivity limit for
+  50h of effective observing time. It takes as input file the
+  output file generated using cspull.
+  If the Python processing module is installed the script can make
+  use of multiple cores and runs in parallel.
+
+make_spectrum.py
+  This script shows how to make a spectrum using obsutils.
+
+make_ts_distributions.py
+  This script creates TS distributions as function of energy for
+  an off observation of 50h effective observing time.
+  If the Python processing module is installed the script can make
+  use of multiple cores and runs in parallel.
+
+benchmark_cta_analysis.py
   Performs performance benchmarks of the various CTA analysis types.
  
 pipeline_binned_disk.py
@@ -28,39 +50,6 @@ pipeline_unbinned_mem.py
   Example pipeline for simulating events and performing an unbinned
   analysis in memory. No files will be written to disk.
  
-
-
-
-make_binned_analysis.py
-  Example Python pipelines for binned analysis.
-  The first pipeline writes intermediate results to the disk while
-  the second pipeline keeps the results in memory.
-  If matplotlib is installed on the system the fitted spectrum
-  and the observed offset angle distribution is displayed in
-  counts space.
-
-make_pull_at_sensitivity_limit.py
-  This script creates pull distributions as function of energy
-  for a test source that is at the 5 sigma sensitivity limit for
-  50h of effective observing time. It takes as input file the
-  output file generated using cspull.
-  If the Python processing module is installed the script can make
-  use of multiple cores and runs in parallel.
-
-make_ts_distributions.py
-  This script creates TS distributions as function of energy for
-  an off observation of 50h effective observing time.
-  If the Python processing module is installed the script can make
-  use of multiple cores and runs in parallel.
-
-make_unbinned_analysis.py
-  Example Python pipelines for unbinned analysis.
-  The first pipeline writes intermediate results to the disk while
-  the second pipeline keeps the results in memory.
-  If matplotlib is installed on the system the fitted spectrum
-  and the observed offset angle distribution is displayed in
-  counts space.
-
 show_butterfly.py
   Displays a butterfly generated using the ctbutterfly tool.
  
 #! /usr/bin/env python
 # ==========================================================================
 # This script peforms a benchmark of the various CTA analysis types
-# (unbinned, binned, cube-style), and if matplotlib is installed, creates
+# (unbinned, binned, stacked), and if matplotlib is installed, creates
 # a plot of the benchmark results.
 #
-# Copyright (C) 2014 Jurgen Knodlseder
+# Copyright (C) 2014-2015 Jurgen Knodlseder
 #
 # 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
-#! /usr/bin/env python
-# ==========================================================================
-# This script illustrates how to build a binned analysis pipeline using the
-# ctools. If matplotlib is installed, a counts spectrum and an offset
-# histogram will be displayed.
-#
-# Copyright (C) 2011-2014 Jurgen Knodlseder
-#
-# 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 gammalib
-import ctools
-import os
-import glob
-import sys
-import math
-try:
-    import matplotlib.pyplot as plt
-    has_matplotlib = True
-except:
-    has_matplotlib = False
-
-
-# ============================= #
-# Analysis pipeline - version 1 #
-# ============================= #
-def pipeline_v1():
-    """
-    Binned analysis pipeline - save intermediate results.
-    
-    This function implements an analysis pipeline that successively calls
-    ctobssim, ctbin and ctlike by saving the intermediate results as FITS
-    files on disk. This replicates the way how the analysis would be done
-    in a ftools-like approach.
-    """
-    # Set script parameters
-    model_name  = "${CTOOLS}/share/models/crab.xml"
-    events_name = "events.fits"
-    cntmap_name = "cntmap.fits"
-    result_name = "results.xml"
-    caldb       = "${CALDB}/data/cta/dummy/bcf"
-    irf         = "cta_dummy_irf"
-    ra          =   83.63
-    dec         =   22.01
-    rad_sim     =   10.0
-    tstart      =    0.0
-    tstop       = 1800.0
-    emin        =    0.1
-    emax        =  100.0
-    enumbins    =   20
-    nxpix       =  200
-    nypix       =  200
-    binsz       =    0.02
-    coordsys    = "CEL"
-    proj        = "CAR"
-
-    # Initialise timing
-    wall_seconds = 0.0
-    cpu_seconds  = 0.0
-
-    # Simulate events
-    sim = ctools.ctobssim()
-    sim.logFileOpen()  # We need this to explicitely open the log file in Python mode
-    sim["inmodel"].filename(model_name)
-    sim["outevents"].filename(events_name)
-    sim["caldb"].string(caldb)
-    sim["irf"].string(irf)
-    sim["ra"].real(ra)
-    sim["dec"].real(dec)
-    sim["rad"].real(rad_sim)
-    sim["tmin"].real(tstart)
-    sim["tmax"].real(tstop)
-    sim["emin"].real(emin)
-    sim["emax"].real(emax)
-    sim.execute()
-    sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n")
-
-    # Update timing
-    wall_seconds += sim.telapse()
-    cpu_seconds  += sim.celapse()
-    
-    # Bin events into counts map
-    bin = ctools.ctbin()
-    bin.logFileOpen()  # We need this to explicitely open the log file in Python mode
-    bin["inobs"].filename(events_name)
-    bin["outcube"].filename(cntmap_name)
-    bin["ebinalg"].string("LOG")
-    bin["emin"].real(emin)
-    bin["emax"].real(emax)
-    bin["enumbins"].integer(enumbins)
-    bin["nxpix"].integer(nxpix)
-    bin["nypix"].integer(nypix)
-    bin["binsz"].real(binsz)
-    bin["coordsys"].string(coordsys)
-    bin["xref"].real(ra)
-    bin["yref"].real(dec)
-    bin["proj"].string(proj)
-    bin.execute()
-    sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n")
-
-    # Update timing
-    wall_seconds += bin.telapse()
-    cpu_seconds  += bin.celapse()
-    
-    # Perform maximum likelihood fitting
-    like = ctools.ctlike()
-    like.logFileOpen()  # We need this to explicitely open the log file in Python mode
-    like["inobs"].filename(cntmap_name)
-    like["inmodel"].filename(model_name)
-    like["outmodel"].filename(result_name)
-    like["caldb"].string(caldb)
-    like["irf"].string(irf)
-    like.execute()
-    sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n")
-
-    # Update timing
-    wall_seconds += like.telapse()
-    cpu_seconds  += like.celapse()
-    
-    # Show total times
-    sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n")
-    sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n")
-
-    # Return
-    return
-
-
-# ============================= #
-# Analysis pipeline - version 2 #
-# ============================= #
-def pipeline_v2():
-    """
-    Binned analysis pipeline - keep intermediate results in memory.
-    
-    This function implements an analysis pipeline that successively calls
-    ctobssim, ctbin and ctlike without saving the intermediate results as
-    FITS files on disk. All data are hold in memory.
-    
-    At the end, results are plotted (if matplotlib is installed)
-    """
-    # Set script parameters
-    model_name  = "${CTOOLS}/share/models/crab.xml"
-    caldb       = "${CALDB}/data/cta/dummy/bcf"
-    irf         = "cta_dummy_irf"
-    ra          =   83.63
-    dec         =   22.01
-    rad_sim     =   10.0
-    tstart      =    0.0
-    tstop       = 1800.0
-    emin        =    0.1
-    emax        =  100.0
-    enumbins    =   20
-    nxpix       =  200
-    nypix       =  200
-    binsz       =    0.02
-    coordsys    = "CEL"
-    proj        = "CAR"
-
-    # Initialise timing
-    wall_seconds = 0.0
-    cpu_seconds  = 0.0
-
-    # Simulate events
-    sim = ctools.ctobssim()
-    sim["inmodel"].filename(model_name)
-    sim["caldb"].string(caldb)
-    sim["irf"].string(irf)
-    sim["ra"].real(ra)
-    sim["dec"].real(dec)
-    sim["rad"].real(rad_sim)
-    sim["tmin"].real(tstart)
-    sim["tmax"].real(tstop)
-    sim["emin"].real(emin)
-    sim["emax"].real(emax)
-    sim.run()
-    sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n")
-
-    # Update timing
-    wall_seconds += sim.telapse()
-    cpu_seconds  += sim.celapse()
-
-    # Bin events into counts map
-    bin = ctools.ctbin(sim.obs())
-    bin["ebinalg"].string("LOG")
-    bin["emin"].real(emin)
-    bin["emax"].real(emax)
-    bin["enumbins"].integer(enumbins)
-    bin["nxpix"].integer(nxpix)
-    bin["nypix"].integer(nypix)
-    bin["binsz"].real(binsz)
-    bin["coordsys"].string(coordsys)
-    bin["xref"].real(ra)
-    bin["yref"].real(dec)
-    bin["proj"].string(proj)
-    bin.run()
-    sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n")
-
-    # Update timing
-    wall_seconds += bin.telapse()
-    cpu_seconds  += bin.celapse()
-
-    # Perform maximum likelihood fitting
-    like = ctools.ctlike(bin.obs())
-    like.run()
-    sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n")
-
-    # Update timing
-    wall_seconds += like.telapse()
-    cpu_seconds  += like.celapse()
-        
-    # Show total times
-    sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n")
-    sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n")
-
-    # Show model fitting results
-    #sys.stdout.write(like.obs().models()+"\n")
-
-    # Plot counts
-    if has_matplotlib:
-        plot_counts(bin.obs())
-    else:
-        sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n")
-    
-    # Return
-    return
-
-
-# =========== #
-# Plot counts #
-# =========== #
-def plot_counts(observations):
-    """
-    Plot counts.
-    """
-    # Set legend fontsize
-    params = {'legend.fontsize': 10}
-    plt.rcParams.update(params)
-
-    # Set plot styles
-    styles = ['b-', 'g-', 'y-', 'n-']
-
-    # Dump header
-    sys.stdout.write("\n")
-    sys.stdout.write("Make plots (using matplotlib):\n")
-    sys.stdout.write("==============================\n")
-    
-    # Create figure 1
-    plt.figure(1,figsize=(12,6))
-    plt.subplots_adjust(hspace=.7)
-
-    # Create subplot 1
-    plt.subplot(121)
-    plt.title("Spectrum (summed over all pixels)")
-    
-    # Loop over observations
-    for obs in observations:
-    
-        # Get event cube
-        cube = obs.events()
-        
-        # Create energy axis
-        energy  = []
-        ebounds = cube.ebounds()
-        for i in range(ebounds.size()):
-            energy.append(ebounds.elogmean(i).TeV())
-
-        # Create spectrum
-        sys.stdout.write("Extract data:\n")
-        counts = [0.0 for i in range(ebounds.size())]
-        for bin in cube:
-            index         = ebounds.index(bin.energy())
-            counts[index] = counts[index] + bin.counts()
-
-        # Create error bars
-        error = [math.sqrt(c) for c in counts]
-
-        # Plot spectrum
-        plt.loglog(energy, counts, 'ro', label='data')
-        plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
-
-        # Extract models
-        sys.stdout.write("Extract models:\n")
-        sum_model = [0.0 for i in range(ebounds.size())]
-        for k, m in enumerate(observations.models()):
-            sys.stdout.write("- "+m.name()+"\n")
-            model = [0.0 for i in range(ebounds.size())]
-            for bin in cube:
-                index        = ebounds.index(bin.energy())
-                prob         = m.eval(bin, obs)
-                model[index] = model[index] + prob * bin.size()
-            for i in range(ebounds.size()):
-                sum_model[i] = sum_model[i] + model[i]
-            plt.loglog(energy, model, styles[k], label=m.name())
-        plt.loglog(energy, sum_model, 'r-', label='total')
-
-    # Put labels
-    plt.xlabel("Energy (TeV)")
-    plt.ylabel("Counts")
-    plt.legend(loc="lower left")
-
-    # Create subplot 2
-    plt.subplot(122)
-    plt.title("Offset (summed over all energies)")
-
-    # Set Crab direction
-    crab = gammalib.GSkyDir()
-    crab.radec_deg(83.63, 22.01)
-    
-    # Loop over observations
-    for obs in observations:
-    
-        # Get event cube
-        cube = obs.events()
-        
-        # Create offset histogram
-        sys.stdout.write("Extract data:\n")
-        nx       = 30
-        doffset2 = 0.01
-        offset2  = [(i+0.5)*doffset2 for i in range(nx)]
-        counts   = [0.0  for i in range(nx)]
-        for bin in cube:
-            off   = gammalib.GCTAInstDir(bin.dir()).dir().dist_deg(crab)
-            off2  = off*off
-            index = int(off2/doffset2)
-            if index < nx:
-                counts[index] = counts[index] + bin.counts()
-
-        # Create error bars
-        error = [math.sqrt(c) for c in counts]
-
-        # Plot distribution
-        plt.semilogy(offset2, counts, 'ro', label='data')
-
-        # Extract models
-        sys.stdout.write("Extract models:\n")
-        sum_model = [0.0 for i in range(nx)]
-        for k, m in enumerate(observations.models()):
-            sys.stdout.write("- "+m.name()+"\n")
-            model = [0.0 for i in range(nx)]
-            for bin in cube:
-                off   = gammalib.GCTAInstDir(bin.dir()).dir().dist_deg(crab)
-                off2  = off*off
-                index = int(off2/doffset2)
-                if index < nx:
-                    prob         = m.eval(bin, obs)
-                    model[index] = model[index] + prob * bin.size()
-            for i in range(nx):
-                sum_model[i] = sum_model[i] + model[i]
-            plt.plot(offset2, model, styles[k], label=m.name())
-        plt.plot(offset2, sum_model, 'r-', label='total')
-        plt.ylim(ymin=0.1)
-
-    # Put labels
-    plt.xlabel("Offset (deg^2)")
-    plt.ylabel("Counts")
-    plt.legend(loc="upper right")
-
-    # Show counts spectra
-    plt.show()
-
-    # Return
-    return
-
-
-#==========================#
-# Main routine entry point #
-#==========================#
-if __name__ == '__main__':
-    """
-    CTA binned analysis script.
-    
-    This script implements the analysis steps to perform a binned analysis.
-    The steps include:
-    - Event simulation
-    - Event binning
-    - Binned maximum likelihood fitting
-    
-    Two variations of the script are implemented. One that calls the classes
-    like executables, storing the intermediate results on disk, and a second
-    that uses the results of the precedent class directly without storing
-    any data on disk.
-    
-    The second variation will also show a counts spectrum of the data after
-    fitting (only if matplotlib is available on your system). Note that the
-    counts spectrum is summed over the entire data space, i.e. no event selection
-    around the Crab has been performed.
-    """
-    # Dump header
-    sys.stdout.write("*************************************\n")
-    sys.stdout.write("*    CTA binned analysis scripts    *\n")
-    sys.stdout.write("*************************************\n")
-    sys.stdout.write("... this script will take 1-2 minutes\n")
-
-    # Remove any existing result files
-    list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
-    for files in list:
-        for file in files:
-            os.remove(file)
-    
-    # Save intermediate results on disk
-    sys.stdout.write("\n")
-    sys.stdout.write("Executable analysis pipeline:\n")
-    sys.stdout.write("=============================\n")
-    pipeline_v1()
-
-    # Do analysis in memory
-    sys.stdout.write("\n")
-    sys.stdout.write("In memory analysis pipeline:\n")
-    sys.stdout.write("============================\n")
-    pipeline_v2()
-    sys.stdout.write("\n")
-#! /usr/bin/env python
-# ==========================================================================
-# This script illustrates how to perform a stacked cube analysis.
-#
-# Copyright (C) 2014 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 ctools
-import gammalib
-from ctools import obsutils
-
-
-# ================= #
-# Analysis pipeline #
-# ================= #
-def pipeline(caldb = "dummy", \
-             irf   = "cta_dummy_irf", \
-             model = "$CTOOLS/share/models/crab.xml"):
-    """
-    Analysis pipeline for a stacked cube analysis. The pipeline writes
-    intermediate results to FITS files on disk. Also log files of all
-    individual ctools will be written.
-    
-    Keywords:
-     caldb - Calibration database repository
-     irf   - Name of instrument response
-     model - Source model
-    """
-    # Set pointing pattern (4 pointings)
-    pointings = obsutils.set_obs_patterns("four", offset=2.0)
-
-    # Set CTA observations (default: 1800 sec per pointing)
-    obs = obsutils.set_obs_list(pointings, caldb=caldb, irf=irf)
-
-    # Load models
-    models = gammalib.GModels(model)
-    obs.models(models)
-
-    # Run ctobssim (simulates 4 event lists)
-    sim = ctools.ctobssim(obs)
-    sim.logFileOpen()
-    sim.run()
-
-    # Run ctbin (builds stacked event cube)
-    bin = ctools.ctbin(sim.obs())
-    bin.logFileOpen()
-    bin["outcube"].filename("cntmap.fits")
-    bin["ebinalg"].string("LOG")
-    bin["emin"].real(0.1)
-    bin["emax"].real(100.0)
-    bin["enumbins"].integer(20)
-    bin["usepnt"].boolean(False)
-    bin["nxpix"].integer(300)
-    bin["nypix"].integer(300)
-    bin["binsz"].real(0.02)
-    bin["coordsys"].string("CEL")
-    bin["proj"].string("CAR")
-    bin["xref"].real(83.6331)
-    bin["yref"].real(22.0145)
-    bin.execute()
-
-    # Run ctexpcube (builds exposure cube)
-    expcube = ctools.ctexpcube(sim.obs())
-    expcube.logFileOpen()
-    expcube["incube"].filename("cntmap.fits")
-    expcube["outcube"].filename("expcube.fits")
-    expcube["caldb"].string(caldb)
-    expcube["irf"].string(irf)
-    expcube.execute()
-
-    # Run ctpsfcube (builds PSF cube)
-    psfcube = ctools.ctpsfcube(sim.obs())
-    psfcube.logFileOpen()
-    psfcube["incube"].filename("NONE")
-    psfcube["outcube"].filename("psfcube.fits")
-    psfcube["caldb"].string(caldb)
-    psfcube["irf"].string(irf)
-    psfcube["ebinalg"].string("LOG")
-    psfcube["emin"].real(0.1)
-    psfcube["emax"].real(100.0)
-    psfcube["enumbins"].integer(20)
-    psfcube["nxpix"].integer(5)
-    psfcube["nypix"].integer(5)
-    psfcube["binsz"].real(2.0)
-    psfcube["coordsys"].string("CEL")
-    psfcube["proj"].string("CAR")
-    psfcube["xref"].real(83.6331)
-    psfcube["yref"].real(22.0145)
-    psfcube["amax"].real(0.3)
-    psfcube["anumbins"].integer(60)
-    psfcube.execute()
-
-    # Run ctbkgcube (builds background cube)
-    bkgcube = ctools.ctbkgcube(sim.obs())
-    bkgcube.logFileOpen()
-    bkgcube["incube"].filename("cntmap.fits")
-    bkgcube["inmodel"].filename(model)
-    bkgcube["outcube"].filename("bkgcube.fits")
-    bkgcube["outmodel"].filename("bkgcube.xml")
-    bkgcube.execute()
-
-    # Create observation container for stacked cube analysis
-    obs.clear()
-    cta = gammalib.GCTAObservation("cntmap.fits", "expcube.fits", "psfcube.fits")
-    obs.append(cta)
-    obs.models(gammalib.GModels("bkgcube.xml"))
-
-    # Run ctlike
-    like = ctools.ctlike(obs)
-    like.logFileOpen()
-    like["debug"].boolean(True)
-    like["outmodel"].filename("cube_results.xml")
-    like.execute()
-
-    # Return
-    return
-
-
-#==========================#
-# Main routine entry point #
-#==========================#
-if __name__ == '__main__':
-    """
-    """
-    # Dump header
-    print("*************************************")
-    print("*     CTA stacked cube analysis     *")
-    print("*************************************")
-
-    # Run pipeline
-    pipeline()
-    
 \ No newline at end of file
-#! /usr/bin/env python
-# ==========================================================================
-# This script illustrates how to build an unbinned analysis pipeline with
-# the ctools. If matplotlib is installed, the spectrum will be displayed
-# on the screen.
-#
-# Required 3rd party modules:
-# None
-#
-# Copyright (C) 2011-2015 Jurgen Knodlseder
-#
-# 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 gammalib
-import ctools
-import math
-import os
-import glob
-import sys
-try:
-    import matplotlib.pyplot as plt
-    has_matplotlib = True
-except:
-    has_matplotlib = False
-
-
-# ============================= #
-# Analysis pipeline - version 1 #
-# ============================= #
-def pipeline_v1():
-    """
-    Unbinned analysis pipeline - save intermediate results.
-    
-    This function implements an analysis pipeline that successively calls
-    ctobssim, ctselect and ctlike by saving the intermediate results as FITS
-    files on disk. This replicates the way how the analysis would be done
-    in a ftools-like approach.
-    """
-    # Set script parameters
-    model_name           = "${CTOOLS}/share/models/crab.xml"
-    events_name          = "events.fits"
-    selected_events_name = "selected_events.fits"
-    result_name          = "results.xml"
-    caldb                = "${CALDB}/data/cta/dummy/bcf"
-    irf                  = "cta_dummy_irf"
-    ra                   =   83.63
-    dec                  =   22.01
-    rad_sim              =   10.0
-    tstart               =    0.0
-    tstop                = 1800.0
-    emin                 =    0.1
-    emax                 =  100.0
-    rad_select           =    3.0
-
-    # Initialise timing
-    wall_seconds = 0.0
-    cpu_seconds  = 0.0
-
-    # Simulate events
-    sim = ctools.ctobssim()
-    sim.logFileOpen()
-    sim["inmodel"].filename(model_name)
-    sim["outevents"].filename(events_name)
-    sim["caldb"].string(caldb)
-    sim["irf"].string(irf)
-    sim["ra"].real(ra)
-    sim["dec"].real(dec)
-    sim["rad"].real(rad_sim)
-    sim["tmin"].real(tstart)
-    sim["tmax"].real(tstop)
-    sim["emin"].real(emin)
-    sim["emax"].real(emax)
-    sim.execute()
-    print("Simulated events ("+str(sim.celapse())+" CPU seconds)")
-    
-    # Update timing
-    wall_seconds += sim.telapse()
-    cpu_seconds  += sim.celapse()
-    
-    # Select events
-    select = ctools.ctselect()
-    select.logFileOpen()
-    select["inobs"].filename(events_name)
-    select["outobs"].filename(selected_events_name)
-    select["ra"].real(ra)
-    select["dec"].real(dec)
-    select["rad"].real(rad_select)
-    select["tmin"].real(tstart)
-    select["tmax"].real(tstop)
-    select["emin"].real(emin)
-    select["emax"].real(emax)
-    select.execute()
-    print("Select event ("+str(select.celapse())+" CPU seconds)")
-
-    # Update timing
-    wall_seconds += select.telapse()
-    cpu_seconds  += select.celapse()
-    
-    # Perform maximum likelihood fitting
-    like = ctools.ctlike()
-    like.logFileOpen()
-    like["inobs"].filename(selected_events_name)
-    like["inmodel"].filename(model_name)
-    like["outmodel"].filename(result_name)
-    like["caldb"].string(caldb)
-    like["irf"].string(irf)
-    like.execute()
-    print("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)")
-
-    # Update timing
-    wall_seconds += like.telapse()
-    cpu_seconds  += like.celapse()
-    
-    # Show total times
-    print("Total wall time elapsed: "+str(wall_seconds)+" seconds")
-    print("Total CPU time used ...: "+str(cpu_seconds)+" seconds")
-
-    # Return
-    return
-
-
-# ============================= #
-# Analysis pipeline - version 2 #
-# ============================= #
-def pipeline_v2(show_data):
-    """
-    Unbinned analysis pipeline - keep intermediate results in memory.
-    
-    This function implements an analysis pipeline that successively calls
-    ctobssim, ctselect and ctlike without saving the intermediate results as
-    FITS files on disk. All data is only hold in memory.
-    """
-    # Set script parameters
-    model_name  = "${CTOOLS}/share/models/crab.xml"
-    caldb       = "${CALDB}/data/cta/dummy/bcf"
-    irf         = "cta_dummy_irf"
-    ra          =   83.63
-    dec         =   22.01
-    rad_sim     =   10.0
-    tstart      =    0.0
-    tstop       = 1800.0
-    emin        =    0.1
-    emax        =  100.0
-    rad_select  =    3.0
-
-    # Initialise timing
-    wall_seconds = 0.0
-    cpu_seconds  = 0.0
-
-    # Simulate events
-    sim = ctools.ctobssim()
-    sim["inmodel"].filename(model_name)
-    sim["caldb"].string(caldb)
-    sim["irf"].string(irf)
-    sim["ra"].real(ra)
-    sim["dec"].real(dec)
-    sim["rad"].real(rad_sim)
-    sim["tmin"].real(tstart)
-    sim["tmax"].real(tstop)
-    sim["emin"].real(emin)
-    sim["emax"].real(emax)
-    sim.run()
-    print("Simulated events ("+str(sim.celapse())+" CPU seconds)")
-
-    # Update timing
-    wall_seconds += sim.telapse()
-    cpu_seconds  += sim.celapse()
-
-    # Select events
-    select = ctools.ctselect(sim.obs())
-    select["ra"].real(ra)
-    select["dec"].real(dec)
-    select["rad"].real(rad_select)
-    select["tmin"].real(tstart)
-    select["tmax"].real(tstop)
-    select["emin"].real(emin)
-    select["emax"].real(emax)
-    select.run()
-    print("Select event ("+str(select.celapse())+" CPU seconds)")
-
-    # Update timing
-    wall_seconds += select.telapse()
-    cpu_seconds  += select.celapse()
-
-    # Perform maximum likelihood fitting
-    like = ctools.ctlike(select.obs())
-    like.run()
-    print("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)")
-
-    # Update timing
-    wall_seconds += like.telapse()
-    cpu_seconds  += like.celapse()
-        
-    # Show total times
-    print("Total wall time elapsed: "+str(wall_seconds)+" seconds")
-    print("Total CPU time used ...: "+str(cpu_seconds)+" seconds")
-
-    # Optionally plot counts
-    if show_data:
-        if has_matplotlib:
-            plot_counts(like.obs())
-        else:
-            sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n")
-    
-    # Return
-    return
-
-
-# =========== #
-# Plot counts #
-# =========== #
-def plot_counts(observations):
-    """
-    Plot counts.
-    """
-    # Set legend fontsize
-    params = {'legend.fontsize': 10}
-    plt.rcParams.update(params)
-
-    # Set plot styles
-    styles = ['b-', 'g-', 'y-', 'n-']
-
-    # Dump header
-    print("")
-    print("Make plots (using matplotlib):")
-    print("==============================")
-    
-    # Create figure 1
-    plt.figure(1,figsize=(12,6))
-    plt.subplots_adjust(hspace=.7)
-
-    # Create subplot 1
-    plt.subplot(121)
-    plt.title("Spectrum (summed over all pixels)")
-    
-    # Loop over observations
-    for obs in observations:
-    
-        # Get event list
-        list = obs.events()
-        
-        # Create energy axis
-        ebounds = gammalib.GEbounds(10, gammalib.GEnergy(0.1, "TeV"), gammalib.GEnergy(100.0, "TeV"))
-        energy  = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
-
-        # Create spectrum
-        print("Extract data:")
-        counts = [0.0 for i in range(ebounds.size())]
-        print(list.size())
-        for atom in list:
-            index         = ebounds.index(atom.energy())
-            counts[index] = counts[index] + 1.0
-
-        # Create error bars
-        error = [math.sqrt(c) for c in counts]
-
-        # Plot spectrum
-        plt.loglog(energy, counts, 'ro', label='data')
-        #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
-
-        # Extract models
-        #print("Extract models:")
-        #sum_model = [0.0 for i in range(ebounds.size())]
-        #for k, m in enumerate(observations.models()):
-        #   print("- "+m.name())
-        #   model = [0.0 for i in range(ebounds.size())]
-        #   for atom in list:
-        #       index        = ebounds.index(atom.energy())
-        #       prob         = m.eval(atom, obs)
-        #       model[index] = model[index] + prob * atom.size()
-        #   for i in range(ebounds.size()):
-        #       sum_model[i] = sum_model[i] + model[i]
-        #   plt.loglog(energy, model, styles[k], label=m.name())
-        #plt.loglog(energy, sum_model, 'r-', label='total')
-
-    # Put labels
-    plt.xlabel("Energy (TeV)")
-    plt.ylabel("Counts")
-    plt.legend(loc="lower left")
-
-    # Create subplot 2
-    plt.subplot(122)
-    plt.title("Offset (summed over all energies)")
-
-    # Set Crab direction
-    crab = gammalib.GSkyDir()
-    crab.radec_deg(83.63, 22.01)
-    
-    # Loop over observations
-    for obs in observations:
-    
-        # Get event cube
-        list = obs.events()
-        
-        # Create offset histogram
-        print "Extract data:"
-        nx       = 30
-        doffset2 = 0.01
-        offset2  = [(i+0.5)*doffset2 for i in range(nx)]
-        counts   = [0.0  for i in range(nx)]
-        for atom in list:
-            off   = gammalib.GCTAInstDir(atom.dir()).dir().dist_deg(crab)
-            off2  = off*off
-            index = int(off2/doffset2)
-            if index < nx:
-                counts[index] = counts[index] + 1.0
-
-        # Create error bars
-        error = [math.sqrt(c) for c in counts]
-
-        # Plot distribution
-        plt.semilogy(offset2, counts, 'ro', label='data')
-
-        # Extract models
-        print "Extract models:"
-        sum_model = [0.0 for i in range(nx)]
-        for k, m in enumerate(observations.models()):
-            print "- "+m.name()
-            model = [0.0 for i in range(nx)]
-            for atom in list:
-                off   = gammalib.GCTAInstDir(atom.dir()).dir().dist_deg(crab)
-                off2  = off*off
-                index = int(off2/doffset2)
-                if index < nx:
-                    prob         = m.eval(atom, obs)
-                    model[index] = model[index] + prob * atom.size()
-            for i in range(nx):
-                sum_model[i] = sum_model[i] + model[i]
-            #plt.plot(offset2, model, styles[k], label=m.name())
-        #plt.plot(offset2, sum_model, 'r-', label='total')
-        #plt.ylim(ymin=0.1)
-
-    # Put labels
-    plt.xlabel("Offset (deg^2)")
-    plt.ylabel("Counts")
-    plt.legend(loc="upper right")
-
-    # Show counts spectra
-    plt.show()
-
-    # Return
-    return
-
-
-#==========================#
-# Main routine entry point #
-#==========================#
-if __name__ == '__main__':
-    """
-    CTA unbinned analysis script.
-    
-    This script implements the analysis steps to perform an unbinned analysis.
-    The steps include:
-    - Event simulation
-    - Event selection
-    - Binned maximum likelihood fitting
-    
-    Two variations of the script are implemented. One that calls the classes
-    like executables, storing the intermediate results on disk, and a second
-    that uses the results of the precedent class directly without storing
-    any data on disk.
-    """
-    # Initialise flags
-    need_help = False
-    show_data = True
-    
-    # Test for command line arguments
-    print sys.argv[0]
-    if (len(sys.argv) > 1):
-        if sys.argv[1] == "-h":
-            need_help = True
-        elif sys.argv[1] == "-noshow":
-            show_data = False
-        else:
-            need_help = True            
-
-    # Print help if needed and exit
-    if need_help:
-        print "Usage: example_binned.py [OPTIONS]"
-        print "     -h       Display this usage message"
-        print "     -noshow  Do not show data using matplotlib"
-        sys.exit()
-
-    # Dump header
-    print "*********************************"
-    print "* CTA unbinned analysis scripts *"
-    print "*********************************"
-    
-    # Remove any existing result files
-    list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
-    for files in list:
-        for file in files:
-            os.remove(file)
-    
-    # Save intermediate results on disk
-    print ""
-    print "Executable analysis pipeline:"
-    print "============================="
-    pipeline_v1()
-    
-    # Do analysis in memory
-    print ""
-    print "In memory analysis pipeline:"
-    print "============================"
-    pipeline_v2(show_data)
-    print ""