From ddd8c148bb4b883264e7f2e41cb58d0683d285c7 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 24 May 2015 22:08:35 +0200 Subject: [PATCH] Adapt script to actual GammaLib interface --- inst/cta/test/benchmark_irf_computation.py | 14 ++++++++------ inst/cta/test/benchmark_ml_fitting.py | 14 ++++++++------ inst/cta/test/example_binned_ml_fit.py | 12 +++++++----- inst/cta/test/example_make_model.py | 44 ++++++++++++++++++++++++++++++-------------- inst/cta/test/example_sim_photons.py | 2 +- inst/cta/test/example_sim_psf.py | 2 +- inst/cta/test/test_model.py | 11 ++++------- inst/mwl/test/test_python.py | 49 ++++++++++++++++++++++++++----------------------- 8 files changed, 85 insertions(+), 63 deletions(-) diff --git a/inst/cta/test/benchmark_irf_computation.py b/inst/cta/test/benchmark_irf_computation.py index 2d6346f..0d5b989 100755 --- a/inst/cta/test/benchmark_irf_computation.py +++ b/inst/cta/test/benchmark_irf_computation.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__ == '__main__': 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__ == '__main__': 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) diff --git a/inst/cta/test/benchmark_ml_fitting.py b/inst/cta/test/benchmark_ml_fitting.py index 2172907..3060542 100755 --- a/inst/cta/test/benchmark_ml_fitting.py +++ b/inst/cta/test/benchmark_ml_fitting.py @@ -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__ == '__main__': 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__ == '__main__': binned_analysis(model, cntmap, irf, gammalib.GCaldb(caldb)) # Perform stacked analysis - stacked_analysis(model, cntmap, expcube, psfcube) + stacked_analysis(model, cntmap, expcube, psfcube, bkgcube) diff --git a/inst/cta/test/example_binned_ml_fit.py b/inst/cta/test/example_binned_ml_fit.py index 2fc15f9..e985016 100755 --- a/inst/cta/test/example_binned_ml_fit.py +++ b/inst/cta/test/example_binned_ml_fit.py @@ -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__ == '__main__': 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) diff --git a/inst/cta/test/example_make_model.py b/inst/cta/test/example_make_model.py index 1d06c2c..f3b655d 100755 --- a/inst/cta/test/example_make_model.py +++ b/inst/cta/test/example_make_model.py @@ -1,8 +1,24 @@ #! /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 . +# +# ========================================================================== +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__ == '__main__': 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)) diff --git a/inst/cta/test/example_sim_photons.py b/inst/cta/test/example_sim_photons.py index 79c69b4..1288413 100755 --- a/inst/cta/test/example_sim_photons.py +++ b/inst/cta/test/example_sim_photons.py @@ -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-') diff --git a/inst/cta/test/example_sim_psf.py b/inst/cta/test/example_sim_psf.py index 7d6c50c..ad033d0 100755 --- a/inst/cta/test/example_sim_psf.py +++ b/inst/cta/test/example_sim_psf.py @@ -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-') diff --git a/inst/cta/test/test_model.py b/inst/cta/test/test_model.py index e94fce9..2d06a54 100755 --- a/inst/cta/test/test_model.py +++ b/inst/cta/test/test_model.py @@ -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 diff --git a/inst/mwl/test/test_python.py b/inst/mwl/test/test_python.py index 49c347b..bb9a6e6 100755 --- a/inst/mwl/test/test_python.py +++ b/inst/mwl/test/test_python.py @@ -1,7 +1,9 @@ #! /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__ == '__main__': 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__ == '__main__': plt.show() except ImportError: - print "Matplotlib is not (correctly) installed on your system." + print("Matplotlib is not (correctly) installed on your system.") -- libgit2 0.21.4