#! /usr/bin/env python # ========================================================================== # Shows one or several CTA response functions. # # Copyright (C) 2014-2017 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 sys import math import gammalib import cscripts try: import matplotlib.pyplot as plt plt.figure() plt.close() except (ImportError, RuntimeError): print('This script needs the "matplotlib" module') sys.exit() # ================================== # # Compute Eq. (17) of Li & Ma (1983) # # ================================== # def sigma_lima(Non, Noff, alpha=0.2): """ Compute Eq. (17) of Li & Ma (1983). Parameters ---------- Non : float Number of On counts Noff : float Number of Off counts alpha : float, optional Ratio of on-to-off exposure Returns ------- sigma : float Detection significance in Gaussian sigma """ # Compute sensitivity in Gaussian sigma alpha1 = alpha + 1.0 Ntotal = Non + Noff arg1 = Non / Ntotal arg2 = Noff / Ntotal term1 = Non * math.log((alpha1/alpha)*arg1) term2 = Noff * math.log(alpha1*arg2) sigma = math.sqrt(2.0 * (term1 + term2)) # Return sensitivity return sigma # ======================================== # # Solve Eq. (17) of Li & Ma (1983) for Non # # ======================================== # def Non_lima(sigma, Noff, alpha=0.2): """ Solve Eq. (17) of Li & Ma (1983) for Non Parameters ---------- sigma : float Required detection significance in Gaussian sigma Noff : float Number of Off counts alpha : float, optional Ratio of on-to-off exposure Returns ------- Non : float Number of On counts """ # Set initial guess for Non S = 25.0 Non = S + alpha*Noff # Compute initial sigma estimate s = sigma_lima(Non, Noff, alpha=alpha) delta = s - sigma # Iterate while abs(delta) > 1.0e-6: S *= sigma/s Non = S + alpha*Noff s = sigma_lima(Non, Noff, alpha=alpha) delta = s - sigma # Return return Non # ======================= # # Show one effective area # # ======================= # def show_one_effective_area(rsp, name, color='r'): """ Show one effective area Parameters ---------- rsp : `~gammalib.GCTAResponse` Response function name : str Name of the response function color : str, optional Color for plot """ # Set to figure 1 plt.figure(1) # Generate logE(TeV) vector nbins = 20 logE = [-1.7+0.2*i for i in range(nbins)] E = [math.pow(10.0, e) for e in logE] # Generate effective area vector aeff = [0.0 for i in range(nbins)] for i in range(nbins): aeff[i] = rsp.aeff()(logE[i]) # Plot data plt.loglog(E, aeff, color+'-', label=name) plt.loglog(E, aeff, color+'o') # Set axes plt.xlabel('Energy (TeV)') plt.ylabel(r'Effective area (cm$^{2}$)') # Set legend plt.legend(loc='lower right') # Return return # ======================== # # Show one background rate # # ======================== # def show_one_background_rate(rsp, name, color='r'): """ Show one background rate Parameters ---------- rsp : `~gammalib.GCTAResponse` Response function name : str Name of the response function color : str, optional Color for plot """ # Set to figure 2 plt.figure(2) # Generate logE(TeV) vector nbins = 20 logE = [-1.7+0.2*i for i in range(nbins)] E = [math.pow(10.0, e) for e in logE] # Generate on-axis background rate vector bgrate = [0.0 for i in range(nbins)] for i in range(nbins): bgrate[i] = rsp.background()(logE[i], 0.0, 0.0) if bgrate[i] < 1.0e-10: bgrate[i] = 0.0 # Plot data plt.loglog(E, bgrate, color+'-', label=name) plt.loglog(E, bgrate, color+'o') # Set axes plt.xlabel('Energy (TeV)') plt.ylabel(r'Background rate (counts s$^{-1}$ MeV$^{-1}$ sr$^{-1}$)') # Set legend plt.legend(loc='lower left') # Return return # ==================== # # Show one sensitivity # # ==================== # def show_one_sensitivity(rsp, name, color='r', duration=180000.0, alpha=0.2, sigma=5.0): """ Show one sensitivity Parameters ---------- rsp : `~gammalib.GCTAResponse` Response function name : str Name of the response function color : str, optional Color for plot duration : float, optional Duration of observation (s) alpha : float, optional Ratio of on-to-off exposure sigma : float, optional Required detection significance in Gaussian sigma """ # Set constants r68_to_sigma = 0.6624305 TeV2erg = 1.0e6 * gammalib.MeV2erg # Set to figure 3 plt.figure(3) # Generate logE(TeV) vector nbins = 20 logE = [-1.7+0.2*i for i in range(nbins)] E = [math.pow(10.0, e) for e in logE] # Initialise sensitivity vector flux = [0.0 for i in range(nbins)] # Loop over all energy bins for i in range(nbins): # Compute the energy width in MeV ewidth = (math.pow(10.0, logE[i]+0.1) - math.pow(10.0, logE[i]-0.1)) * \ 1.0e6 # Extract effective area in cm2 from Instrument Response Function aeff = rsp.aeff()(logE[i]) # Extract 68% containment radius in radians from Instrument Response # Function. # A kluge is used since this information is not directly available, # but the "delta_max()" method returns 5 times the "sigma" of the # Gaussian PSF. Hence dividing by 5 gives the "sigma", and dividing # by the 68% containment radius to sigma conversion factor gives # r68. r68 = rsp.psf().delta_max(logE[i]) / 5.0 / r68_to_sigma # Compute the solid angle of the 68% containment radius solidangle = gammalib.twopi * (1.0 - math.cos(r68)) # Compute the total number of background counts by multiplying the # background rate by the observation duration, the energy width and # the solid angle bgd_counts = rsp.background()(logE[i], 0.0, 0.0) * \ duration * ewidth * solidangle # The number of Off counts is the total number of background counts # divided by the "alpha" parameter Noff = bgd_counts/alpha # If there are Off counts then compute the corresponding On counts # for "sigma" detection significance if Noff > 0: # Compute On counts using Li & Ma Non = Non_lima(sigma, Noff) # Get the source counts by subtracting the background counts # under the source. Make sure that the source counts are at least # 5% of the background counts (low-energy background systematics) # and that there are at least 10 source events (high-energy limit) src_counts = Non - bgd_counts if src_counts < 0.05*bgd_counts: src_counts = 0.05*bgd_counts if src_counts < 10: src_counts = 10.0 # Compute flux under the assumption that the source spectrum is # a power law with index of -2.6. The flux is the source counts # divided by the duration and 68% of the effective area times # a conversion factor. The conversion factor is given by # # conv = 1.6021765 * E^2 / F(E) # # where E is the mean energy of the bin and F(E) is the flux in # the bin emin = gammalib.GEnergy(math.pow(10.0, logE[i]-0.1), 'TeV') emax = gammalib.GEnergy(math.pow(10.0, logE[i]+0.1), 'TeV') epivot = gammalib.GEnergy(math.pow(10.0, logE[i]), 'TeV') plaw = gammalib.GModelSpectralPlaw(1.0e-6, -2.6, epivot) conv = TeV2erg*E[i]*E[i]/plaw.flux(emin, emax) flux[i] = conv * src_counts / (duration * aeff*0.68) # Plot data plt.loglog(E, flux, color+'-', label=name) plt.loglog(E, flux, color+'o') # Set axes plt.xlabel('Energy (TeV)') plt.ylabel(r'Sensitivity (erg cm$^{-2}$ s$^{-1}$)') # Set legend plt.legend(loc='upper right') # Return return # ==================== # # Show one sensitivity # # ==================== # def show_one_response(rspname, dbname, name, rootdir=None, color='r'): """ Show one response. Parameters ---------- rspname : str Response name dbname : str Database name name : str Name of the response function rootdir : str, optional Response root directory color : str, optional Color for plot """ # Set-up calibration database caldb = gammalib.GCaldb() if rootdir != None: caldb.rootdir(rootdir) if gammalib.dir_exists(dbname): caldb.rootdir(dbname) else: caldb.open('cta', dbname) # Load response function rsp = gammalib.GCTAResponseIrf(rspname, caldb) # Show effective area show_one_effective_area(rsp, name, color=color) # Show background rate show_one_background_rate(rsp, name, color=color) # Show sensitivity show_one_sensitivity(rsp, name, color=color) # Return return # ==================== # # Show one sensitivity # # ==================== # def plot_response(plotfile): """ Show response Parameters ---------- plotfile : str Plot filename """ # Create figures plt.figure(1) plt.title('Effective area') plt.figure(2) plt.title('Background rate') plt.figure(3) plt.title('Sensitivity') # Set response dictionary rsps = [{'dbname': 'prod2', 'rspname': 'South_0.5h', 'name': 'Prod2 (30 min)', 'color': 'b'}, {'dbname': 'prod2', 'rspname': 'South_50h', 'name': 'Prod2 (50 h)', 'color': 'r'}] # Loop over all responses for rsp in rsps: # Get attributes dbname = rsp['dbname'] rspname = rsp['rspname'] color = rsp['color'] name = rsp['name'] rootdir = rsp.setdefault('rootdir', None) # Show response show_one_response(rspname, dbname, name, rootdir=rootdir, color=color) # Show figure if len(plotfile) > 0: plt.savefig(plotfile) else: plt.show() # Return return # ============= # # Show response # # ============= # def show_response(): """ Show response """ # Set usage string usage = 'show_response.py [-p plotfile]' # Set default options options = [{'option': '-p', 'value': ''}] # Get arguments and options from command line arguments _, options = cscripts.ioutils.get_args_options(options, usage) # Extract script parameters from options plotfile = options[0]['value'] # Show response plot_response(plotfile) # Return return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': # Show response show_response()