diff --git a/inst/cta/test/test_irf_trafo.py b/inst/cta/test/test_irf_trafo.py new file mode 100755 index 0000000..15c84e7 --- /dev/null +++ b/inst/cta/test/test_irf_trafo.py @@ -0,0 +1,67 @@ +#! /usr/bin/env python +# ===========================================================================================# +# Test coordinate transformations for IRF computation. +# ===========================================================================================# +from gammalib import * +from math import * + + +# ====================== # +# Compute Omega0 # +# ====================== # +def omega0(): + """ + ... + """ + # Set pointing direction + pointing = GSkyDir() + pointing.radec_deg(10.0, 10.0) + + # Set model centre + model = GSkyDir() + model.radec_deg(12.0, 10.0) + + # Set measured photon direction + p = GSkyDir() + p.radec_deg(13.0, 9.0) + + # Determine angular distance between measured photon direction and model centre + zeta = model.dist(p) + + # Determine angular distance between measured photon direction and pointing direction + eta = pointing.dist(p) + + # Determine angular distance between model centre and pointing direction + lamda = model.dist(pointing) + + # Compute Omega0 + omega0 = 0.0 + omega2 = 0.0 + denom = sin(lamda) * sin(zeta) + if denom != 0.0: + arg = (cos(eta) - cos(lamda) * cos(zeta))/denom + arg2 = cos(eta)/denom - 1./(tan(lamda) * tan(zeta)) + omega0 = acos(arg) + omega2 = acos(arg2) + + # Print results + print "Omega0 ...: ", omega0*180.0/pi + print "Omega0 ...: ", omega2*180.0/pi + print "Zeta .....: ", zeta + print "Eta ......: ", eta + print "Lambda ...: ", lamda + + # Return + return + + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + """ + ... + """ + # ... + omega0() diff --git a/inst/cta/test/test_npred_computation.py b/inst/cta/test/test_npred_computation.py new file mode 100755 index 0000000..890ba71 --- /dev/null +++ b/inst/cta/test/test_npred_computation.py @@ -0,0 +1,112 @@ +#! /usr/bin/env python +# ========================================================================== +# This script tests the Npred computation of the CTA response. +# +# Copyright (C) 2011 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 . +# +# ========================================================================== +from gammalib import * +import math +import os +import glob +import sys + +# Optionally import matplotlib if available +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except ImportError: + print "WARNING: Matplotlib not available for plotting." + has_matplotlib = False + + +# ======================================= # +# Compute NPsf as function of PSF position # +# ======================================= # +def npsf(caldb, irf, radius): + """ + """ + # Setup CTA response + rsp = GCTAResponse(irf, caldb) + + # Setup ROI + instDir = GCTAInstDir() + instDir.radec_deg(0.0, 0.0) + roi = GCTARoi() + roi.centre(instDir) + roi.radius(radius) + + # Setup source energy + srcEng = GEnergy() + srcEng.TeV(0.1) + + # Setup dummy time and pointing + srcTime = GTime() + pnt = GCTAPointing() + + # Compute PSF as function of offset + x = [] + y = [] + for i in range(150): + + # Compute offset angle + offset = i*radius/100.0 + + # Set PSF position + srcDir = GSkyDir() + srcDir.radec_deg(0.0, offset) + + # Compute NPsf + npsf = rsp.npsf(srcDir, srcEng.log10TeV(), srcTime, pnt, roi) + + # Store results + x.append(offset) + y.append(npsf) + #print offset, npsf + + # Optionally plot results + if has_matplotlib: + + # Create figure + plt.figure(1) + plt.title("NPsf") + + # Plot data + plt.plot(x, y, 'r-') + plt.plot([radius,radius], [0.0,1.0], 'b-') + + # Set axes + plt.xlabel("Offset angle (deg)") + plt.ylabel("Npsf") + + # Show plot + plt.show() + + + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + """ + """ + # Set response information + caldb = "$GAMMALIB/share/caldb/cta" + irf = "kb_E_50h_v3" + + # Compute NPsf + npsf(caldb, irf, 1.0)