Adapt script to actual GammaLib interface

Jürgen Knödlseder
1 parent e3348b38
Showing 8 changed files with 85 additions and 63 deletions
inst/cta/test/benchmark_irf_computation.py
inst/cta/test/benchmark_ml_fitting.py
inst/cta/test/example_binned_ml_fit.py
inst/cta/test/example_make_model.py
inst/cta/test/example_sim_photons.py
inst/cta/test/example_sim_psf.py
inst/cta/test/test_model.py
inst/mwl/test/test_python.py
@@ -6,7 +6,7 @@
 # Monte-Carlo simulations. This script is used to validate the absolute
 # normalization of the CTA response computations.
 #
-# Copyright (C) 2014 Juergen Knoedlseder
+# Copyright (C) 2014-2015 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
@@ -117,7 +117,7 @@ def binned_irf(model, cntmap, irf, caldb, cntref):
 # =========================== #
 # CTA stacked IRF computation #
 # =========================== #
-def stacked_irf(model, cntmap, expcube, psfcube, cntref):
+def stacked_irf(model, cntmap, expcube, psfcube, bkgcube, cntref):
     """
     Perform stacked IRF computation.
     """
@@ -132,9 +132,10 @@ def stacked_irf(model, cntmap, expcube, psfcube, cntref):
     obs.load(cntmap)
     # Specify response for CTA observation
-    exposure = gammalib.GCTACubeExposure(expcube)
-    psf      = gammalib.GCTACubePsf(psfcube)
-    obs.response(exposure, psf)
+    exposure   = gammalib.GCTACubeExposure(expcube)
+    psf        = gammalib.GCTACubePsf(psfcube)
+    background = gammalib.GCTACubeBackground(bkgcube)
+    obs.response(exposure, psf, background)
     # Load model to describe the data from XML file
     models = gammalib.GModels(model)
@@ -179,6 +180,7 @@ if __name__ == &#39;__main__&#39;:
     caldb   = "../caldb"
     expcube = "data/expcube.fits"
     psfcube = "data/psfcube.fits"
+    bkgcube = "data/bkgcube.fits"
     # Set tests
     tests = ["point", "disk", "gauss", "shell", "ellipse", "diffuse"]
@@ -235,4 +237,4 @@ if __name__ == &#39;__main__&#39;:
         binned_irf(model, cntmap, irf, gammalib.GCaldb(caldb), cntref)
         # Perform stacked computation
-        stacked_irf(model, cntmap, expcube, psfcube, cntref)
+        stacked_irf(model, cntmap, expcube, psfcube, bkgcube, cntref)
@@ -3,7 +3,7 @@
 # This script performs a benchmark for maximum likelihood fitting of CTA
 # data.
 #
-# Copyright (C) 2012-2014 Juergen Knoedlseder
+# Copyright (C) 2012-2015 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
@@ -142,7 +142,7 @@ def binned_analysis(model, cntmap, irf, caldb):
 # ==================== #
 # CTA stacked analysis #
 # ==================== #
-def stacked_analysis(model, cntmap, expcube, psfcube):
+def stacked_analysis(model, cntmap, expcube, psfcube, bkgcube):
     """
     Perform stacked maximum likelihood fitting of CTA data using GammaLib
     classes.
@@ -166,9 +166,10 @@ def stacked_analysis(model, cntmap, expcube, psfcube):
     cta_obs.load(cntmap)
     # Specify response for CTA observation
-    exposure = gammalib.GCTACubeExposure(expcube)
-    psf      = gammalib.GCTACubePsf(psfcube)
-    cta_obs.response(exposure, psf)
+    exposure   = gammalib.GCTACubeExposure(expcube)
+    psf        = gammalib.GCTACubePsf(psfcube)
+    background = gammalib.GCTACubeBackground(bkgcube)
+    cta_obs.response(exposure, psf, background)
     # Append CTA observation to observation container
     obs.append(cta_obs)
@@ -217,6 +218,7 @@ if __name__ == &#39;__main__&#39;:
     caldb   = "../caldb"
     expcube = "data/expcube.fits"
     psfcube = "data/psfcube.fits"
+    bkgcube = "data/bkgcube.fits"
     # Set tests
     tests = ["point", "disk", "gauss", "shell", "ellipse", "diffuse"]
@@ -273,4 +275,4 @@ if __name__ == &#39;__main__&#39;:
         binned_analysis(model, cntmap, irf, gammalib.GCaldb(caldb))
         # Perform stacked analysis
-        stacked_analysis(model, cntmap, expcube, psfcube)
+        stacked_analysis(model, cntmap, expcube, psfcube, bkgcube)
@@ -115,7 +115,7 @@ def binned_analysis(model, cntmap, irf, caldb):
 # ==================== #
 # CTA stacked analysis #
 # ==================== #
-def stacked_analysis(model, cntmap, expcube, psfcube):
+def stacked_analysis(model, cntmap, expcube, psfcube, bkgcube):
     """
     Perform stacked maximum likelihood fitting of CTA data using GammaLib
     classes. The function performs two fits, one with the full model and
@@ -141,9 +141,10 @@ def stacked_analysis(model, cntmap, expcube, psfcube):
     cta_obs.load(cntmap)
     # Specify response for CTA observation
-    exposure = gammalib.GCTACubeExposure(expcube)
-    psf      = gammalib.GCTACubePsf(psfcube)
-    cta_obs.response(exposure, psf)
+    exposure   = gammalib.GCTACubeExposure(expcube)
+    psf        = gammalib.GCTACubePsf(psfcube)
+    background = gammalib.GCTACubeBackground(bkgcube)
+    cta_obs.response(exposure, psf, background)
     # Append CTA observation to observation container
     obs.append(cta_obs)
@@ -225,12 +226,13 @@ if __name__ == &#39;__main__&#39;:
     cntmap  = "data/crab_cntmap.fits"
     expcube = "data/expcube.fits"
     psfcube = "data/psfcube.fits"
+    bkgcube = "data/bkgcube.fits"
     # Perform binned analysis
     results_binned = binned_analysis(model, cntmap, irf, gammalib.GCaldb(caldb))
     # Perform stacked analysis
-    results_cube = stacked_analysis(model, cntmap, expcube, psfcube)
+    results_cube = stacked_analysis(model, cntmap, expcube, psfcube, bkgcube)
     # Print model results
     #print(result)
 #! /usr/bin/env python
-# ===========================================================================================#
+# ==========================================================================
 # This script makes a model map from an XML file.
-# ===========================================================================================#
-from gammalib import *
+#
+# Copyright (C) 2015 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 gammalib
 # ============================ #
@@ -13,16 +29,16 @@ def make_model_for_cntmap(cntmap, modname, xmlname, irf, caldb, clobber=True):
     Make a model counts map from an XML model.
     """
     # Allocate empty CTA observation
-    obs = GCTAObservation()
+    obs = gammalib.GCTAObservation()
     # Load counts map into CTA observation
-    obs.load_binned(cntmap)
+    obs.load(cntmap)
     # Set response
     obs.response(irf, caldb)
     # Load models from XML file
-    models = GModels(xmlname)
+    models = gammalib.GModels(xmlname)
     # Loop over all bins in counts map
     for event in obs.events():
@@ -44,18 +60,18 @@ if __name__ == &#39;__main__&#39;:
     Example illustrating how to create a model map for a counts map.
     """
     # Dump header
-    print
-    print "***************************************************"
-    print "* Create model map for counts map using XML model *"
-    print "***************************************************"
-    print "... please wait for a few seconds"
+    print("")
+    print("***************************************************")
+    print("* Create model map for counts map using XML model *")
+    print("***************************************************")
+    print("... please wait for a few seconds")
     # Set parameters
-    irf     = "kb_E_50h_v3"
-    caldb   = "../caldb"
+    irf     = "cta_dummy_irf"
+    caldb   = "./caldb"
     xmlname = "data/crab.xml"
     cntmap  = "data/crab_cntmap.fits"
     modname = "model.fits"
     # Make the model
-    make_model_for_cntmap(cntmap, modname, xmlname, irf, caldb)
+    make_model_for_cntmap(cntmap, modname, xmlname, irf, gammalib.GCaldb(caldb))
@@ -116,7 +116,7 @@ def show_photons(photons, xmlname, e_min, e_max, area, duration, ebins=30):
         # Plot data
         plt.loglog(energy, counts, 'ro')
-        plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
+        plt.errorbar(energy, counts, error, ecolor='r')
         # Plot model
         plt.plot(energy, model, 'b-')
@@ -89,7 +89,7 @@ def sim_psf(response, energy, r_max=0.8, rbins=1000, nmc=1000000):
         # Plot simulated data
         plt.plot(delta, counts, 'ro')
-        plt.errorbar(delta, counts, error, fmt=None, ecolor='r')
+        plt.errorbar(delta, counts, error, ecolor='r')
         # Plot PSF
         plt.plot(delta, psf, 'b-')
@@ -174,10 +174,10 @@ def observation(ra=83.6331, dec=22.0145, emin=0.1, emax=100.0,
     obs = gammalib.GCTAObservation()
     # Set response
-    #obs.response("cta_dummy_irf", GCaldb("../caldb"))
-    exposure = gammalib.GCTACubeExposure("data/expcube.fits")
-    psf      = gammalib.GCTACubePsf("data/psfcube.fits")
-    obs.response(exposure, psf)
+    exposure   = gammalib.GCTACubeExposure("data/expcube.fits")
+    psf        = gammalib.GCTACubePsf("data/psfcube.fits")
+    background = gammalib.GCTACubeBackground("data/bkgcube.fits")
+    obs.response(exposure, psf, background)
     # Set pointing
     dir = gammalib.GSkyDir()
@@ -201,9 +201,6 @@ def observation(ra=83.6331, dec=22.0145, emin=0.1, emax=100.0,
     obs.livetime(duration * deadc)
     obs.deadc(deadc)
-    # Optionally show observation
-    # print(obs)
-
     # Return observation
     return obs
 #! /usr/bin/env python
-from gammalib import *
-from math import *
+import gammalib
+import math
+#from gammalib import *
+#from math import *
 import os
@@ -17,16 +19,16 @@ def extract_data(filename):
     styles = ['ro', 'bo', 'go']
     # Open FITS file
-    fits = GFits(filename)
+    fits = gammalib.GFits(filename)
     # Loop over all table HDUs
     k = 0
     for extno in range(fits.size()):
-        if fits.at(extno).exttype() == GFitsHDU.HT_ASCII_TABLE or \
-           fits.at(extno).exttype() == GFitsHDU.HT_BIN_TABLE:
+        if fits.at(extno).exttype() == gammalib.GFitsHDU.HT_ASCII_TABLE or \
+           fits.at(extno).exttype() == gammalib.GFitsHDU.HT_BIN_TABLE:
             # Load spectrum
-            spectrum = GMWLSpectrum()
+            spectrum = gammalib.GMWLSpectrum()
             spectrum.load(filename, extno)
             # Add spectral points
@@ -42,7 +44,7 @@ def extract_data(filename):
                 y.append(yval * conv)
                 dy.append(yerr * conv)
             plt.loglog(x, y, styles[k])
-            plt.errorbar(x, y, dy, fmt=None, ecolor=ecolor[k])
+            plt.errorbar(x, y, dy, ecolor=ecolor[k])
             # Increment style counter
             k = k + 1
@@ -63,37 +65,38 @@ def fit_spectrum(filename, xmlname):
     Fit spectrum.
     """
     # Load observations
-    obs = GObservations()
-    comptel = GMWLObservation(filename, "COMPTEL")
-    lat = GMWLObservation(filename, "LAT")
-    hess = GMWLObservation(filename, "HESS")
+    obs     = gammalib.GObservations()
+    comptel = gammalib.GMWLObservation(filename, "COMPTEL")
+    lat     = gammalib.GMWLObservation(filename, "LAT")
+    hess    = gammalib.GMWLObservation(filename, "HESS")
     obs.append(comptel)
     obs.append(lat)
     obs.append(hess)
     # Setup model
-    models = GModels(xmlname)
+    models = gammalib.GModels(xmlname)
     obs.models(models)
     # Perform optimization
-    opt = GOptimizerLM()
+    opt = gammalib.GOptimizerLM()
     opt.max_iter(1000)
     obs.optimize(opt)
-    print obs
+    print(obs)
     # Get model values
     x = np.power(10., np.arange(-1., 8., 0.1))
     y = []
-    d = GSkyDir()
+    d = gammalib.GSkyDir()
     d.radec_deg(83.6331, 22.0145)
-    e = GEnergy()
-    t = GTime()
+    e = gammalib.GEnergy()
+    t = gammalib.GTime()
     for energy in x:
         e.MeV(energy)
         yval = 0.0
         for model in obs.models():
             if model.type() == "PointSource":
-                yval += model.value(GPhoton(d, e, t)) * energy * energy * 1.6021765e-6
+                yval += model.value(gammalib.GPhoton(d, e, t)) * \
+                        energy * energy * 1.6021765e-6
         y.append(yval)
     plt.plot(x, y)
@@ -106,10 +109,10 @@ if __name__ == &#39;__main__&#39;:
     Perform testing.
     """
     # Dump result
-    print
-    print "****************************"
-    print "* Python interface testing *"
-    print "****************************"
+    print("")
+    print("****************************")
+    print("* Python interface testing *")
+    print("****************************")
     # Only proceed if matplotlib is available
     try:
@@ -131,4 +134,4 @@ if __name__ == &#39;__main__&#39;:
         plt.show()
     except ImportError:
-        print "Matplotlib is not (correctly) installed on your system."
+        print("Matplotlib is not (correctly) installed on your system.")