Veronika Vodeb / ctools

Commit e16a786351dad4a157db455ed25590ae945dcf69

Authored by Jürgen Knödlseder
1 parent 72358246

Add IRF display script

Inline Side-by-side
Showing 1 changed file with 548 additions and 0 deletions
examples/show_irf.py 0 → 100755
  View file @e16a786
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Display Instrument Response Function
  4 +#
  5 +# Copyright (C) 2017 Juergen Knoedlseder
  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 +import sys
  22 +import math
  23 +import gammalib
  24 +import cscripts
  25 +try:
  26 + import matplotlib.pyplot as plt
  27 + plt.figure()
  28 + plt.close()
  29 +except (ImportError, RuntimeError):
  30 + print('This script needs the "matplotlib" module')
  31 + sys.exit()
  32 +
  33 +
  34 +# ========= #
  35 +# Plot Aeff #
  36 +# ========= #
  37 +def plot_aeff(sub, aeff, emin=None, emax=None, tmin=None, tmax=None,
  38 + nengs=100, nthetas=100):
  39 + """
  40 + Plot effective area
  41 +
  42 + Parameters
  43 + ----------
  44 + sub : figure
  45 + Subplot
  46 + aeff : `~gammalib.GCTAAeff2d`
  47 + Instrument Response Function
  48 + emin : float, optional
  49 + Minimum energy (TeV)
  50 + emax : float, optional
  51 + Maximum energy (TeV)
  52 + tmin : float, optional
  53 + Minimum offset angle (deg)
  54 + tmax : float, optional
  55 + Maximum offset angle (deg)
  56 + nengs : int, optional
  57 + Number of energies
  58 + nthetas : int, optional
  59 + Number of offset angles
  60 + """
  61 + # Determine energy range
  62 + ieng = aeff.table().axis('ENERG')
  63 + neng = aeff.table().axis_bins(ieng)
  64 + if emin == None:
  65 + emin = aeff.table().axis_lo(ieng, 0)
  66 + if emax == None:
  67 + emax = aeff.table().axis_hi(ieng, neng-1)
  68 +
  69 + # Determine offset angle range
  70 + itheta = aeff.table().axis('THETA')
  71 + ntheta = aeff.table().axis_bins(itheta)
  72 + if tmin == None:
  73 + tmin = aeff.table().axis_lo(itheta, 0)
  74 + if tmax == None:
  75 + tmax = aeff.table().axis_hi(itheta, ntheta-1)
  76 +
  77 + # Use log energies
  78 + emin = math.log10(emin)
  79 + emax = math.log10(emax)
  80 +
  81 + # Set axes
  82 + denergy = (emax - emin)/(nengs-1)
  83 + dtheta = (tmax - tmin)/(nthetas-1)
  84 + logenergies = [emin+i*denergy for i in range(nengs)]
  85 + thetas = [tmax-i*dtheta for i in range(nthetas)]
  86 +
  87 + # Initialise image
  88 + image = []
  89 +
  90 + # Loop over offset angles
  91 + for theta in thetas:
  92 +
  93 + # Initialise row
  94 + row = []
  95 +
  96 + # Loop over energies
  97 + for logenergy in logenergies:
  98 +
  99 + # Get effective area value
  100 + value = aeff(logenergy, theta*gammalib.deg2rad)
  101 +
  102 + # Append value
  103 + row.append(value)
  104 +
  105 + # Append row
  106 + image.append(row)
  107 +
  108 + # Plot image
  109 + sub.imshow(image, extent=[emin,emax,tmin,tmax], aspect=0.5)
  110 +
  111 + # Show boundary contours
  112 + contours = sub.contour(logenergies, thetas, image, [0.0], colors=('white'))
  113 + sub.clabel(contours, inline=1, fontsize=8)
  114 +
  115 + # Plot title and axis
  116 + sub.set_title('Effective area')
  117 + sub.set_xlabel('log10(E/TeV)')
  118 + sub.set_ylabel('Offset angle (deg)')
  119 +
  120 + # Return
  121 + return
  122 +
  123 +
  124 +# ======== #
  125 +# Plot PSF #
  126 +# ======== #
  127 +def plot_psf(sub, psf, emin=None, emax=None, tmin=None, tmax=None,
  128 + nengs=100, nthetas=100):
  129 + """
  130 + Plot Point Spread Function
  131 +
  132 + Parameters
  133 + ----------
  134 + sub : figure
  135 + Subplot
  136 + psf : `~gammalib.GCTAPsf2D`
  137 + Instrument Response Function
  138 + emin : float, optional
  139 + Minimum energy (TeV)
  140 + emax : float, optional
  141 + Maximum energy (TeV)
  142 + tmin : float, optional
  143 + Minimum offset angle (deg)
  144 + tmax : float, optional
  145 + Maximum offset angle (deg)
  146 + nengs : int, optional
  147 + Number of energies
  148 + nthetas : int, optional
  149 + Number of offset angles
  150 + """
  151 + # Determine energy range
  152 + ieng = psf.table().axis('ENERG')
  153 + neng = psf.table().axis_bins(ieng)
  154 + if emin == None:
  155 + emin = psf.table().axis_lo(ieng, 0)
  156 + if emax == None:
  157 + emax = psf.table().axis_hi(ieng, neng-1)
  158 +
  159 + # Determine offset angle range
  160 + itheta = psf.table().axis('THETA')
  161 + ntheta = psf.table().axis_bins(itheta)
  162 + if tmin == None:
  163 + tmin = psf.table().axis_lo(itheta, 0)
  164 + if tmax == None:
  165 + tmax = psf.table().axis_hi(itheta, ntheta-1)
  166 +
  167 + # Use log energies
  168 + emin = math.log10(emin)
  169 + emax = math.log10(emax)
  170 +
  171 + # Set axes
  172 + denergy = (emax - emin)/(nengs-1)
  173 + dtheta = (tmax - tmin)/(nthetas-1)
  174 + logenergies = [emin+i*denergy for i in range(nengs)]
  175 + thetas = [tmax-i*dtheta for i in range(nthetas)]
  176 +
  177 + # Initialise image
  178 + image = []
  179 +
  180 + # Loop over offset angles
  181 + for theta in thetas:
  182 +
  183 + # Initialise row
  184 + row = []
  185 +
  186 + # Loop over energies
  187 + for logenergy in logenergies:
  188 +
  189 + # Get containment radius value
  190 + value = psf.containment_radius(0.68, logenergy, theta*gammalib.deg2rad) * \
  191 + gammalib.rad2deg
  192 +
  193 + # Limit containment radius values to [0.0,0.5]
  194 + if value < 0.0:
  195 + value = 0.0
  196 + if value > 0.5:
  197 + value = 0.5
  198 +
  199 + # Append value
  200 + row.append(value)
  201 +
  202 + # Append row
  203 + image.append(row)
  204 +
  205 + # Plot image
  206 + sub.imshow(image, extent=[emin,emax,tmin,tmax], aspect=0.5)
  207 +
  208 + # Show boundary contours
  209 + contours = sub.contour(logenergies, thetas, image, [0.0], colors=('white'))
  210 + sub.clabel(contours, inline=1, fontsize=8)
  211 +
  212 + # Plot title and axis
  213 + if psf.classname() == 'GCTAPsfKing':
  214 + sub.set_title('King function PSF 68% containment radius')
  215 + else:
  216 + sub.set_title('Gaussian PSF 68% containment radius')
  217 + sub.set_xlabel('log10(E/TeV)')
  218 + sub.set_ylabel('Offset angle (deg)')
  219 +
  220 + # Return
  221 + return
  222 +
  223 +
  224 +# ====================== #
  225 +# Plot energy dispersion #
  226 +# ====================== #
  227 +def plot_edisp(sub, edisp, emin=None, emax=None, tmin=None, tmax=None,
  228 + nengs=100, nthetas=100):
  229 + """
  230 + Plot Background template
  231 +
  232 + Parameters
  233 + ----------
  234 + sub : figure
  235 + Subplot
  236 + edisp : `~gammalib.GCTAEdisp2D`
  237 + Instrument Response Function
  238 + emin : float, optional
  239 + Minimum energy (TeV)
  240 + emax : float, optional
  241 + Maximum energy (TeV)
  242 + tmin : float, optional
  243 + Minimum offset angle (deg)
  244 + tmax : float, optional
  245 + Maximum offset angle (deg)
  246 + nengs : int, optional
  247 + Number of energies
  248 + nthetas : int, optional
  249 + Number of offset angles
  250 + """
  251 + # Determine energy range
  252 + ieng = edisp.table().axis('ETRUE')
  253 + neng = edisp.table().axis_bins(ieng)
  254 + if emin == None:
  255 + emin = edisp.table().axis_lo(ieng, 0)
  256 + if emax == None:
  257 + emax = edisp.table().axis_hi(ieng, neng-1)
  258 +
  259 + # Determine offset angle range
  260 + itheta = edisp.table().axis('THETA')
  261 + ntheta = edisp.table().axis_bins(itheta)
  262 + if tmin == None:
  263 + tmin = edisp.table().axis_lo(itheta, 0)
  264 + if tmax == None:
  265 + tmax = edisp.table().axis_hi(itheta, ntheta-1)
  266 +
  267 + # Use log energies
  268 + emin = math.log10(emin)
  269 + emax = math.log10(emax)
  270 +
  271 + # Set axes
  272 + denergy = (emax - emin)/(nengs-1)
  273 + dtheta = (tmax - tmin)/(nthetas-1)
  274 + logenergies = [emin+i*denergy for i in range(nengs)]
  275 + thetas = [tmax-i*dtheta for i in range(nthetas)]
  276 +
  277 + # Initialise image
  278 + image = []
  279 +
  280 + # Loop over offset angles
  281 + for theta in thetas:
  282 +
  283 + # Initialise row
  284 + row = []
  285 +
  286 + # Compute detx and dety
  287 + detx = theta*gammalib.deg2rad
  288 + dety = 0.0
  289 +
  290 + # Loop over energies
  291 + for logenergy in logenergies:
  292 +
  293 + # Get no migration energy dispersion value
  294 + value = edisp(logenergy, logenergy, theta*gammalib.deg2rad)
  295 +
  296 + # Append value
  297 + row.append(value)
  298 +
  299 + # Append row
  300 + image.append(row)
  301 +
  302 + # Plot image
  303 + sub.imshow(image, extent=[emin,emax,tmin,tmax], aspect=0.5)
  304 +
  305 + # Show boundary contours
  306 + contours = sub.contour(logenergies, thetas, image, [0.0], colors=('white'))
  307 + sub.clabel(contours, inline=1, fontsize=8)
  308 +
  309 + # Plot title and axis
  310 + sub.set_title('No migration energy dispersion')
  311 + sub.set_xlabel('log10(E/TeV)')
  312 + sub.set_ylabel('Offset angle (deg)')
  313 +
  314 + # Return
  315 + return
  316 +
  317 +
  318 +# =============== #
  319 +# Plot Background #
  320 +# =============== #
  321 +def plot_bkg(sub, bkg, emin=None, emax=None, tmin=None, tmax=None,
  322 + nengs=100, nthetas=100):
  323 + """
  324 + Plot Background template
  325 +
  326 + Parameters
  327 + ----------
  328 + sub : figure
  329 + Subplot
  330 + bkg : `~gammalib.GCTABackground3D`
  331 + Instrument Response Function
  332 + emin : float, optional
  333 + Minimum energy (TeV)
  334 + emax : float, optional
  335 + Maximum energy (TeV)
  336 + tmin : float, optional
  337 + Minimum offset angle (deg)
  338 + tmax : float, optional
  339 + Maximum offset angle (deg)
  340 + nengs : int, optional
  341 + Number of energies
  342 + nthetas : int, optional
  343 + Number of offset angles
  344 + """
  345 + # Determine energy range
  346 + ieng = bkg.table().axis('ENERG')
  347 + neng = bkg.table().axis_bins(ieng)
  348 + if emin == None:
  349 + emin = bkg.table().axis_lo(ieng, 0)
  350 + if emax == None:
  351 + emax = bkg.table().axis_hi(ieng, neng-1)
  352 +
  353 + # Determine offset angle range
  354 + if tmin == None:
  355 + tmin = 0.0
  356 + if tmax == None:
  357 + tmax = 6.0
  358 +
  359 + # Use log energies
  360 + emin = math.log10(emin)
  361 + emax = math.log10(emax)
  362 +
  363 + # Set axes
  364 + denergy = (emax - emin)/(nengs-1)
  365 + dtheta = (tmax - tmin)/(nthetas-1)
  366 + logenergies = [emin+i*denergy for i in range(nengs)]
  367 + thetas = [tmax-i*dtheta for i in range(nthetas)]
  368 +
  369 + # Initialise image
  370 + image = []
  371 +
  372 + # Loop over offset angles
  373 + for theta in thetas:
  374 +
  375 + # Initialise row
  376 + row = []
  377 +
  378 + # Compute detx and dety
  379 + detx = theta*gammalib.deg2rad
  380 + dety = 0.0
  381 +
  382 + # Loop over energies
  383 + for logenergy in logenergies:
  384 +
  385 + # Get containment radius value
  386 + value = bkg(logenergy, detx, dety)
  387 +
  388 + # Append value
  389 + row.append(value)
  390 +
  391 + # Append row
  392 + image.append(row)
  393 +
  394 + # Plot image
  395 + sub.imshow(image, extent=[emin,emax,tmin,tmax], aspect=0.5)
  396 +
  397 + # Show boundary contours
  398 + contours = sub.contour(logenergies, thetas, image, [0.0], colors=('white'))
  399 + sub.clabel(contours, inline=1, fontsize=8)
  400 +
  401 + # Plot title and axis
  402 + sub.set_title('Background acceptance')
  403 + sub.set_xlabel('log10(E/TeV)')
  404 + sub.set_ylabel('Offset angle (deg)')
  405 +
  406 + # Return
  407 + return
  408 +
  409 +
  410 +# ======== #
  411 +# Plot IRF #
  412 +# ======== #
  413 +def plot_irf(irf, emin, emax, tmin, tmax, plotfile):
  414 + """
  415 + Plot Instrument Response Function
  416 +
  417 + Parameters
  418 + ----------
  419 + irf : `~gammalib.GCTAResponseIrf`
  420 + Instrument Response Function
  421 + emin : float
  422 + Minimum energy (TeV)
  423 + emax : float
  424 + Maximum energy (TeV)
  425 + tmin : float
  426 + Minimum offset angle (deg)
  427 + tmax : float
  428 + Maximum offset angle (deg)
  429 + plotfile : str
  430 + Plot filename
  431 + """
  432 + # Build selection string
  433 + selection = ''
  434 + eselection = ''
  435 + tselection = ''
  436 + if emin != None and emax != None:
  437 + eselection += '%.3f-%.1f TeV' % (emin, emax)
  438 + elif emin != None:
  439 + eselection += ' >%.3f TeV' % (emin)
  440 + elif emax != None:
  441 + eselection += ' <%.1f TeV' % (emax)
  442 + if tmin != None and tmax != None:
  443 + tselection += '%.1f-%.1f deg' % (tmin, tmax)
  444 + elif tmin != None:
  445 + tselection += ' >%.1f deg' % (tmin)
  446 + elif tmax != None:
  447 + tselection += ' <%.1f deg' % (tmax)
  448 + if len(eselection) > 0 and len(tselection) > 0:
  449 + selection = ' (%s, %s)' % (eselection, tselection)
  450 + elif len(eselection) > 0:
  451 + selection = ' (%s)' % (eselection)
  452 + elif len(tselection) > 0:
  453 + selection = ' (%s)' % (tselection)
  454 +
  455 + # Build title
  456 + mission = irf.caldb().mission()
  457 + instrument = irf.caldb().instrument()
  458 + response = irf.rspname()
  459 + title = '%s "%s" Instrument Response Function "%s"%s' % \
  460 + (gammalib.toupper(mission), instrument, response, selection)
  461 +
  462 + # Create figure
  463 + fig = plt.figure(figsize=(12,8))
  464 +
  465 + # Add title
  466 + fig.suptitle(title, fontsize=16)
  467 +
  468 + # Plot Aeff
  469 + ax1 = fig.add_subplot(221)
  470 + plot_aeff(ax1, irf.aeff(), emin=emin, emax=emax, tmin=tmin, tmax=tmax)
  471 +
  472 + # Plot Psf
  473 + ax2 = fig.add_subplot(222)
  474 + plot_psf(ax2, irf.psf(), emin=emin, emax=emax, tmin=tmin, tmax=tmax)
  475 +
  476 + # Plot Edisp
  477 + ax3 = fig.add_subplot(223)
  478 + plot_edisp(ax3, irf.edisp(), emin=emin, emax=emax, tmin=tmin, tmax=tmax)
  479 +
  480 + # Plot Background
  481 + ax4 = fig.add_subplot(224)
  482 + plot_bkg(ax4, irf.background(), emin=emin, emax=emax, tmin=tmin, tmax=tmax)
  483 +
  484 + # Show plots or save it into file
  485 + if len(plotfile) > 0:
  486 + plt.savefig(plotfile)
  487 + else:
  488 + plt.show()
  489 +
  490 + # Return
  491 + return
  492 +
  493 +
  494 +# ======== #
  495 +# Show IRF #
  496 +# ======== #
  497 +def show_irf():
  498 + """
  499 + Show Instrument Response Function
  500 + """
  501 + # Set usage string
  502 + usage = 'show_irf.py [-p plotfile] caldb irf'
  503 +
  504 + # Set default options
  505 + options = [{'option': '-p', 'value': ''},
  506 + {'option': '-emin', 'value': None},
  507 + {'option': '-emax', 'value': None},
  508 + {'option': '-tmin', 'value': None},
  509 + {'option': '-tmax', 'value': None}]
  510 +
  511 + # Get arguments and options from command line arguments
  512 + args, options = cscripts.ioutils.get_args_options(options, usage)
  513 +
  514 + # Extract script parameters from options
  515 + plotfile = options[0]['value']
  516 + emin = options[1]['value']
  517 + emax = options[2]['value']
  518 + tmin = options[3]['value']
  519 + tmax = options[4]['value']
  520 +
  521 + # Convert limits to float
  522 + if emin != None:
  523 + emin = float(emin)
  524 + if emax != None:
  525 + emax = float(emax)
  526 + if tmin != None:
  527 + tmin = float(tmin)
  528 + if tmax != None:
  529 + tmax = float(tmax)
  530 +
  531 + # Get IRF
  532 + caldb = gammalib.GCaldb('cta', args[0])
  533 + irf = gammalib.GCTAResponseIrf(args[1], caldb)
  534 +
  535 + # Plot IRF
  536 + plot_irf(irf, emin, emax, tmin, tmax, plotfile)
  537 +
  538 + # Return
  539 + return
  540 +
  541 +
  542 +# ======================== #
  543 +# Main routine entry point #
  544 +# ======================== #
  545 +if __name__ == '__main__':
  546 +
  547 + # Show IRF
  548 + show_irf()
... ...