Commit a39a072768b16995c09eebed62804c8fd3bbe77b
1 parent
b3ae9f06
initial version
Showing
2 changed files
with
179 additions
and
0 deletions
inst/cta/test/test_irf_trafo.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ===========================================================================================# | ||
| 3 | +# Test coordinate transformations for IRF computation. | ||
| 4 | +# ===========================================================================================# | ||
| 5 | +from gammalib import * | ||
| 6 | +from math import * | ||
| 7 | + | ||
| 8 | + | ||
| 9 | +# ====================== # | ||
| 10 | +# Compute Omega0 # | ||
| 11 | +# ====================== # | ||
| 12 | +def omega0(): | ||
| 13 | + """ | ||
| 14 | + ... | ||
| 15 | + """ | ||
| 16 | + # Set pointing direction | ||
| 17 | + pointing = GSkyDir() | ||
| 18 | + pointing.radec_deg(10.0, 10.0) | ||
| 19 | + | ||
| 20 | + # Set model centre | ||
| 21 | + model = GSkyDir() | ||
| 22 | + model.radec_deg(12.0, 10.0) | ||
| 23 | + | ||
| 24 | + # Set measured photon direction | ||
| 25 | + p = GSkyDir() | ||
| 26 | + p.radec_deg(13.0, 9.0) | ||
| 27 | + | ||
| 28 | + # Determine angular distance between measured photon direction and model centre | ||
| 29 | + zeta = model.dist(p) | ||
| 30 | + | ||
| 31 | + # Determine angular distance between measured photon direction and pointing direction | ||
| 32 | + eta = pointing.dist(p) | ||
| 33 | + | ||
| 34 | + # Determine angular distance between model centre and pointing direction | ||
| 35 | + lamda = model.dist(pointing) | ||
| 36 | + | ||
| 37 | + # Compute Omega0 | ||
| 38 | + omega0 = 0.0 | ||
| 39 | + omega2 = 0.0 | ||
| 40 | + denom = sin(lamda) * sin(zeta) | ||
| 41 | + if denom != 0.0: | ||
| 42 | + arg = (cos(eta) - cos(lamda) * cos(zeta))/denom | ||
| 43 | + arg2 = cos(eta)/denom - 1./(tan(lamda) * tan(zeta)) | ||
| 44 | + omega0 = acos(arg) | ||
| 45 | + omega2 = acos(arg2) | ||
| 46 | + | ||
| 47 | + # Print results | ||
| 48 | + print "Omega0 ...: ", omega0*180.0/pi | ||
| 49 | + print "Omega0 ...: ", omega2*180.0/pi | ||
| 50 | + print "Zeta .....: ", zeta | ||
| 51 | + print "Eta ......: ", eta | ||
| 52 | + print "Lambda ...: ", lamda | ||
| 53 | + | ||
| 54 | + # Return | ||
| 55 | + return | ||
| 56 | + | ||
| 57 | + | ||
| 58 | + | ||
| 59 | +#==========================# | ||
| 60 | +# Main routine entry point # | ||
| 61 | +#==========================# | ||
| 62 | +if __name__ == '__main__': | ||
| 63 | + """ | ||
| 64 | + ... | ||
| 65 | + """ | ||
| 66 | + # ... | ||
| 67 | + omega0() |
inst/cta/test/test_npred_computation.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ========================================================================== | ||
| 3 | +# This script tests the Npred computation of the CTA response. | ||
| 4 | +# | ||
| 5 | +# Copyright (C) 2011 Jurgen Knodlseder | ||
| 6 | +# | ||
| 7 | +# This program is free software: you can redistribute it and/or modify | ||
| 8 | +# it under the terms of the GNU General Public License as published by | ||
| 9 | +# the Free Software Foundation, either version 3 of the License, or | ||
| 10 | +# (at your option) any later version. | ||
| 11 | +# | ||
| 12 | +# This program is distributed in the hope that it will be useful, | ||
| 13 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 14 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 15 | +# GNU General Public License for more details. | ||
| 16 | +# | ||
| 17 | +# You should have received a copy of the GNU General Public License | ||
| 18 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 19 | +# | ||
| 20 | +# ========================================================================== | ||
| 21 | +from gammalib import * | ||
| 22 | +import math | ||
| 23 | +import os | ||
| 24 | +import glob | ||
| 25 | +import sys | ||
| 26 | + | ||
| 27 | +# Optionally import matplotlib if available | ||
| 28 | +try: | ||
| 29 | + import matplotlib.pyplot as plt | ||
| 30 | + has_matplotlib = True | ||
| 31 | +except ImportError: | ||
| 32 | + print "WARNING: Matplotlib not available for plotting." | ||
| 33 | + has_matplotlib = False | ||
| 34 | + | ||
| 35 | + | ||
| 36 | +# ======================================= # | ||
| 37 | +# Compute NPsf as function of PSF position # | ||
| 38 | +# ======================================= # | ||
| 39 | +def npsf(caldb, irf, radius): | ||
| 40 | + """ | ||
| 41 | + """ | ||
| 42 | + # Setup CTA response | ||
| 43 | + rsp = GCTAResponse(irf, caldb) | ||
| 44 | + | ||
| 45 | + # Setup ROI | ||
| 46 | + instDir = GCTAInstDir() | ||
| 47 | + instDir.radec_deg(0.0, 0.0) | ||
| 48 | + roi = GCTARoi() | ||
| 49 | + roi.centre(instDir) | ||
| 50 | + roi.radius(radius) | ||
| 51 | + | ||
| 52 | + # Setup source energy | ||
| 53 | + srcEng = GEnergy() | ||
| 54 | + srcEng.TeV(0.1) | ||
| 55 | + | ||
| 56 | + # Setup dummy time and pointing | ||
| 57 | + srcTime = GTime() | ||
| 58 | + pnt = GCTAPointing() | ||
| 59 | + | ||
| 60 | + # Compute PSF as function of offset | ||
| 61 | + x = [] | ||
| 62 | + y = [] | ||
| 63 | + for i in range(150): | ||
| 64 | + | ||
| 65 | + # Compute offset angle | ||
| 66 | + offset = i*radius/100.0 | ||
| 67 | + | ||
| 68 | + # Set PSF position | ||
| 69 | + srcDir = GSkyDir() | ||
| 70 | + srcDir.radec_deg(0.0, offset) | ||
| 71 | + | ||
| 72 | + # Compute NPsf | ||
| 73 | + npsf = rsp.npsf(srcDir, srcEng.log10TeV(), srcTime, pnt, roi) | ||
| 74 | + | ||
| 75 | + # Store results | ||
| 76 | + x.append(offset) | ||
| 77 | + y.append(npsf) | ||
| 78 | + #print offset, npsf | ||
| 79 | + | ||
| 80 | + # Optionally plot results | ||
| 81 | + if has_matplotlib: | ||
| 82 | + | ||
| 83 | + # Create figure | ||
| 84 | + plt.figure(1) | ||
| 85 | + plt.title("NPsf") | ||
| 86 | + | ||
| 87 | + # Plot data | ||
| 88 | + plt.plot(x, y, 'r-') | ||
| 89 | + plt.plot([radius,radius], [0.0,1.0], 'b-') | ||
| 90 | + | ||
| 91 | + # Set axes | ||
| 92 | + plt.xlabel("Offset angle (deg)") | ||
| 93 | + plt.ylabel("Npsf") | ||
| 94 | + | ||
| 95 | + # Show plot | ||
| 96 | + plt.show() | ||
| 97 | + | ||
| 98 | + | ||
| 99 | + | ||
| 100 | + | ||
| 101 | +#==========================# | ||
| 102 | +# Main routine entry point # | ||
| 103 | +#==========================# | ||
| 104 | +if __name__ == '__main__': | ||
| 105 | + """ | ||
| 106 | + """ | ||
| 107 | + # Set response information | ||
| 108 | + caldb = "$GAMMALIB/share/caldb/cta" | ||
| 109 | + irf = "kb_E_50h_v3" | ||
| 110 | + | ||
| 111 | + # Compute NPsf | ||
| 112 | + npsf(caldb, irf, 1.0) |