Commit e69cd7444bfc8da7b010647b7f56de10a8493cc2
1 parent
89694521
Add script to show CTA responses using standard interfaces.
Showing
1 changed file
with
292 additions
and
0 deletions
examples/show_response.py
0 → 100755
| 1 | +#! /usr/bin/env python | ||
| 2 | +# ========================================================================== | ||
| 3 | +# This script displays one or several CTA response functions. | ||
| 4 | +# | ||
| 5 | +# Required 3rd party modules: | ||
| 6 | +# - matplotlib | ||
| 7 | +# - numpy | ||
| 8 | +# | ||
| 9 | +# Copyright (C) 2014 Juergen Knoedlseder | ||
| 10 | +# | ||
| 11 | +# This program is free software: you can redistribute it and/or modify | ||
| 12 | +# it under the terms of the GNU General Public License as published by | ||
| 13 | +# the Free Software Foundation, either version 3 of the License, or | ||
| 14 | +# (at your option) any later version. | ||
| 15 | +# | ||
| 16 | +# This program is distributed in the hope that it will be useful, | ||
| 17 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 18 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 19 | +# GNU General Public License for more details. | ||
| 20 | +# | ||
| 21 | +# You should have received a copy of the GNU General Public License | ||
| 22 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. | ||
| 23 | +# | ||
| 24 | +# ========================================================================== | ||
| 25 | +import sys | ||
| 26 | +import math | ||
| 27 | +import gammalib | ||
| 28 | +import matplotlib.pyplot as plt | ||
| 29 | + | ||
| 30 | + | ||
| 31 | +# ================================== # | ||
| 32 | +# Compute Eq. (17) of Li & Ma (1983) # | ||
| 33 | +# ================================== # | ||
| 34 | +def sigma_lima(Non, Noff, alpha=0.2): | ||
| 35 | + """ | ||
| 36 | + Compute Eq. (17) of Li & Ma (1983). | ||
| 37 | + | ||
| 38 | + Parameters: | ||
| 39 | + Non - Number of on counts | ||
| 40 | + Noff - Number of off counts | ||
| 41 | + Keywords: | ||
| 42 | + alpha - Ratio of on-to-off exposure | ||
| 43 | + """ | ||
| 44 | + # Compute sensitivity | ||
| 45 | + alpha1 = alpha + 1.0 | ||
| 46 | + sum = Non + Noff | ||
| 47 | + arg1 = Non / sum | ||
| 48 | + arg2 = Noff / sum | ||
| 49 | + term1 = Non * math.log((alpha1/alpha)*arg1) | ||
| 50 | + term2 = Noff * math.log(alpha1*arg2) | ||
| 51 | + sigma = math.sqrt(2.0 * (term1 + term2)) | ||
| 52 | + | ||
| 53 | + # Return sensitivity | ||
| 54 | + return sigma | ||
| 55 | + | ||
| 56 | + | ||
| 57 | +# ======================================== # | ||
| 58 | +# Solve Eq. (17) of Li & Ma (1983) for Non # | ||
| 59 | +# ======================================== # | ||
| 60 | +def Non_lima(sigma, Noff, alpha=0.2): | ||
| 61 | + """ | ||
| 62 | + Solve Eq. (17) of Li & Ma (1983) for Non | ||
| 63 | + | ||
| 64 | + Parameters: | ||
| 65 | + sigma - Required significance | ||
| 66 | + Noff - Number of off counts | ||
| 67 | + Keywords: | ||
| 68 | + alpha - Ratio of on-to-off exposure | ||
| 69 | + """ | ||
| 70 | + # Set initial guess for Non | ||
| 71 | + S = 25.0 | ||
| 72 | + Non = S + alpha*Noff | ||
| 73 | + | ||
| 74 | + # Compute initial sigma estimate | ||
| 75 | + s = sigma_lima(Non, Noff, alpha=alpha) | ||
| 76 | + delta = s - sigma | ||
| 77 | + | ||
| 78 | + # Iterate | ||
| 79 | + while abs(delta) > 1.0e-6: | ||
| 80 | + S *= sigma/s | ||
| 81 | + Non = S + alpha*Noff | ||
| 82 | + s = sigma_lima(Non, Noff, alpha=alpha) | ||
| 83 | + delta = s - sigma | ||
| 84 | + | ||
| 85 | + # Return | ||
| 86 | + return Non | ||
| 87 | + | ||
| 88 | + | ||
| 89 | +# ======================= # | ||
| 90 | +# Show one effective area # | ||
| 91 | +# ======================= # | ||
| 92 | +def show_one_effective_area(rsp, name, color="r"): | ||
| 93 | + """ | ||
| 94 | + Show one effective area. | ||
| 95 | + """ | ||
| 96 | + # Set to figure 1 | ||
| 97 | + plt.figure(1) | ||
| 98 | + | ||
| 99 | + # Generate logE(TeV) vector | ||
| 100 | + nbins = 20 | ||
| 101 | + logE = [-1.7+0.2*i for i in range(nbins)] | ||
| 102 | + E = [math.pow(10.0, e) for e in logE] | ||
| 103 | + | ||
| 104 | + # Generate effective area vector | ||
| 105 | + aeff = [0.0 for i in range(nbins)] | ||
| 106 | + for i in range(nbins): | ||
| 107 | + aeff[i] = rsp.aeff()(logE[i]) | ||
| 108 | + | ||
| 109 | + # Plot data | ||
| 110 | + plt.loglog(E, aeff, color+'-', label=name) | ||
| 111 | + plt.loglog(E, aeff, color+'o') | ||
| 112 | + | ||
| 113 | + # Set axes | ||
| 114 | + plt.xlabel("Energy [TeV]") | ||
| 115 | + plt.ylabel("Effective area [cm2]") | ||
| 116 | + | ||
| 117 | + # Set legend | ||
| 118 | + plt.legend(loc="lower right") | ||
| 119 | + | ||
| 120 | + # Return | ||
| 121 | + return | ||
| 122 | + | ||
| 123 | + | ||
| 124 | +# ======================== # | ||
| 125 | +# Show one background rate # | ||
| 126 | +# ======================== # | ||
| 127 | +def show_one_background_rate(rsp, name, color="r"): | ||
| 128 | + """ | ||
| 129 | + Show one background rate. | ||
| 130 | + """ | ||
| 131 | + # Set to figure 2 | ||
| 132 | + plt.figure(2) | ||
| 133 | + | ||
| 134 | + # Generate logE(TeV) vector | ||
| 135 | + nbins = 20 | ||
| 136 | + logE = [-1.7+0.2*i for i in range(nbins)] | ||
| 137 | + E = [math.pow(10.0, e) for e in logE] | ||
| 138 | + | ||
| 139 | + # Generate on-axis background rate vector | ||
| 140 | + bgrate = [0.0 for i in range(nbins)] | ||
| 141 | + for i in range(nbins): | ||
| 142 | + bgrate[i] = rsp.background()(logE[i], 0.0, 0.0) | ||
| 143 | + if bgrate[i] < 1.0e-10: | ||
| 144 | + bgrate[i] = 0.0 | ||
| 145 | + | ||
| 146 | + # Plot data | ||
| 147 | + plt.loglog(E, bgrate, color+'-', label=name) | ||
| 148 | + plt.loglog(E, bgrate, color+'o') | ||
| 149 | + | ||
| 150 | + # Set axes | ||
| 151 | + plt.xlabel("Energy [TeV]") | ||
| 152 | + plt.ylabel("Background rate [events/s/MeV/sr]") | ||
| 153 | + | ||
| 154 | + # Set legend | ||
| 155 | + plt.legend(loc="lower left") | ||
| 156 | + | ||
| 157 | + # Return | ||
| 158 | + return | ||
| 159 | + | ||
| 160 | + | ||
| 161 | +# ==================== # | ||
| 162 | +# Show one sensitivity # | ||
| 163 | +# ==================== # | ||
| 164 | +def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sigma=5.0): | ||
| 165 | + """ | ||
| 166 | + Show one sensitivity. | ||
| 167 | + """ | ||
| 168 | + # Set constants | ||
| 169 | + r68_to_sigma = 0.6624305 | ||
| 170 | + | ||
| 171 | + # Set to figure 3 | ||
| 172 | + plt.figure(3) | ||
| 173 | + | ||
| 174 | + # Generate logE(TeV) vector | ||
| 175 | + nbins = 20 | ||
| 176 | + logE = [-1.7+0.2*i for i in range(nbins)] | ||
| 177 | + E = [math.pow(10.0, e) for e in logE] | ||
| 178 | + | ||
| 179 | + # Generate sensitivity vector | ||
| 180 | + flux = [0.0 for i in range(nbins)] | ||
| 181 | + for i in range(nbins): | ||
| 182 | + ewidth = (math.pow(10.0, logE[i]+0.1) - \ | ||
| 183 | + math.pow(10.0, logE[i]-0.1))*1.0e6 | ||
| 184 | + aeff = rsp.aeff()(logE[i]) | ||
| 185 | + r68 = rsp.psf().delta_max(logE[i]) / 5.0 / r68_to_sigma | ||
| 186 | + solidangle = gammalib.twopi * (1.0 - math.cos(r68)) | ||
| 187 | + bgd_counts = rsp.background()(logE[i], 0.0, 0.0) * duration * ewidth * solidangle | ||
| 188 | + Noff = bgd_counts/alpha | ||
| 189 | + if Noff > 0: | ||
| 190 | + Non = Non_lima(sigma, Noff) | ||
| 191 | + src_counts = Non - bgd_counts | ||
| 192 | + flux[i] = src_counts / (aeff * duration) * E[i]*1.0e6 * gammalib.MeV2erg | ||
| 193 | + | ||
| 194 | + # Plot data | ||
| 195 | + plt.loglog(E, flux, color+'-', label=name) | ||
| 196 | + plt.loglog(E, flux, color+'o') | ||
| 197 | + | ||
| 198 | + # Set axes | ||
| 199 | + plt.xlabel("Energy [TeV]") | ||
| 200 | + plt.ylabel("Sensitivity [erg/cm2/s]") | ||
| 201 | + | ||
| 202 | + # Set legend | ||
| 203 | + plt.legend(loc="upper right") | ||
| 204 | + | ||
| 205 | + # Return | ||
| 206 | + return | ||
| 207 | + | ||
| 208 | + | ||
| 209 | +# ==================== # | ||
| 210 | +# Show one sensitivity # | ||
| 211 | +# ==================== # | ||
| 212 | +def show_one_response(rspname, dbname, name, color="r"): | ||
| 213 | + """ | ||
| 214 | + Show one response. | ||
| 215 | + """ | ||
| 216 | + # Set-up calibration database | ||
| 217 | + caldb = gammalib.GCaldb() | ||
| 218 | + if gammalib.dir_exists(dbname): | ||
| 219 | + caldb.rootdir(dbname) | ||
| 220 | + else: | ||
| 221 | + caldb.open("cta", dbname) | ||
| 222 | + | ||
| 223 | + # Load response function | ||
| 224 | + rsp = gammalib.GCTAResponse(rspname, caldb) | ||
| 225 | + | ||
| 226 | + # Show effective area | ||
| 227 | + show_one_effective_area(rsp, name, color=color) | ||
| 228 | + | ||
| 229 | + # Show background rate | ||
| 230 | + show_one_background_rate(rsp, name, color=color) | ||
| 231 | + | ||
| 232 | + # Show sensitivity | ||
| 233 | + show_one_sensitivity(rsp, name, color=color) | ||
| 234 | + | ||
| 235 | + # Return | ||
| 236 | + return | ||
| 237 | + | ||
| 238 | + | ||
| 239 | +# ======================== # | ||
| 240 | +# Main routine entry point # | ||
| 241 | +# ======================== # | ||
| 242 | +if __name__ == '__main__': | ||
| 243 | + """ | ||
| 244 | + Display spectrum in PHA format. | ||
| 245 | + """ | ||
| 246 | + # Print usage information | ||
| 247 | + #usage = "Usage: show_pha filename" \ | ||
| 248 | + # " [-n bins] [-c column] [-t title] [-p plot]" | ||
| 249 | + #if len(sys.argv) < 2: | ||
| 250 | + # print(usage) | ||
| 251 | + # sys.exit() | ||
| 252 | + | ||
| 253 | + # Extract parameters | ||
| 254 | + #filename = sys.argv[1] | ||
| 255 | + | ||
| 256 | + # Create figures | ||
| 257 | + plt.figure(1) | ||
| 258 | + plt.title("Effective area") | ||
| 259 | + plt.figure(2) | ||
| 260 | + plt.title("Background rate") | ||
| 261 | + | ||
| 262 | + # Set response dictionary | ||
| 263 | + rsps = [ | ||
| 264 | + {'dbname': "/usr/local/gamma/share/caldb/data/cta/dummy", \ | ||
| 265 | + 'rspname': "cta_dummy_irf", \ | ||
| 266 | + 'name': "Prod1 (Array E, MPIK)", \ | ||
| 267 | + 'color': "b"}, \ | ||
| 268 | + {'dbname': "e", \ | ||
| 269 | + 'rspname': "IFAE20120510_50h", \ | ||
| 270 | + 'name': "Prod1 (Array E, IFAE)", \ | ||
| 271 | + 'color': "g"}, \ | ||
| 272 | + {'dbname': "aar", \ | ||
| 273 | + 'rspname': "DESY20140105_50h", \ | ||
| 274 | + 'name': "Prod2 (Aar, DESY)", \ | ||
| 275 | + 'color': "r"} \ | ||
| 276 | + ] | ||
| 277 | + | ||
| 278 | + # Loop over all responses | ||
| 279 | + for rsp in rsps: | ||
| 280 | + | ||
| 281 | + # Get attributes | ||
| 282 | + dbname = rsp['dbname'] | ||
| 283 | + rspname = rsp['rspname'] | ||
| 284 | + color = rsp['color'] | ||
| 285 | + name = rsp['name'] | ||
| 286 | + | ||
| 287 | + # Show response | ||
| 288 | + show_one_response(rspname, dbname, name, color=color) | ||
| 289 | + | ||
| 290 | + # Show plot | ||
| 291 | + plt.show() | ||
| 292 | + | ||
| 0 | \ No newline at end of file | 293 | \ No newline at end of file |