From e67369a1fa26fff55ed00f30db2b06483c43ef72 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 14 Aug 2016 03:21:41 +0200 Subject: [PATCH] Increase comment density (#1828) --- examples/show_obs.py | 17 +++++++++-------- examples/show_response.py | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 62 insertions(+), 17 deletions(-) diff --git a/examples/show_obs.py b/examples/show_obs.py index f674a9b..1cb6de7 100755 --- a/examples/show_obs.py +++ b/examples/show_obs.py @@ -61,7 +61,7 @@ def run_csobsinfo(filename, ra, dec, debug=True): info['inobs'] = obsdef.url() info['debug'] = debug - # Set offset + # Set offset if the "ra" and "dec" arguments are not empty strings if ra != '' and dec != '': info['ra'] = float(ra) info['dec'] = float(dec) @@ -94,14 +94,15 @@ def plot_information(info, ra, dec, plotfile): plotfile : str Plot filename """ - # Retrieve observation info + # Retrieve observation information zeniths = info.zeniths() azimuths = info.azimuths() offsets = info.offsets() ebounds = info.ebounds() gti = info.gti() - # Create figure with subplots + # Create figure with subplots. In case that offset and energy boundary + # information is available a column is added to the subplots nrows = 2 ncols = 2 if info['offset'].boolean() and ebounds.size() > 0: @@ -110,7 +111,7 @@ def plot_information(info, ra, dec, plotfile): plt.figure() plt.subplot(nrows, ncols, iplot) - # Plot zenith angle distribution + # Plot zenith angle histogram zmin = min(zeniths) zmax = max(zeniths) plt.hist(zeniths, bins=30, range=(zmin, zmax), fc='blue') @@ -118,7 +119,7 @@ def plot_information(info, ra, dec, plotfile): plt.ylabel('Frequency') plt.title('Zenith angle distribution') - # Plot azimuth angle distribution + # Plot azimuth angle histogram iplot += 1 plt.subplot(nrows, ncols, iplot) amin = min(azimuths) @@ -128,7 +129,7 @@ def plot_information(info, ra, dec, plotfile): plt.ylabel('Frequency') plt.title('Azimuth distribution') - # Plot offset if possible + # If available, plot offset from target histogram if info['offset'].boolean(): iplot += 1 plt.subplot(nrows, ncols, iplot) @@ -139,7 +140,7 @@ def plot_information(info, ra, dec, plotfile): plt.ylabel('Frequency') plt.title('Offset distribution') - # Plot energy thresholds if possible + # If available, plot energy thresholds histogram if ebounds.size() > 0: emin = [] emax = [] @@ -155,7 +156,7 @@ def plot_information(info, ra, dec, plotfile): plt.legend(loc='upper left') plt.title('Energy threshold') - # Plot observation point in time wrt zenith angle + # Plot zenith angle versus mean observing time for all observations iplot += 1 plt.subplot(nrows, ncols, iplot) times = [] diff --git a/examples/show_response.py b/examples/show_response.py index ae62b97..0386108 100755 --- a/examples/show_response.py +++ b/examples/show_response.py @@ -52,7 +52,7 @@ def sigma_lima(Non, Noff, alpha=0.2): sigma : float Detection significance in Gaussian sigma """ - # Compute sensitivity + # Compute sensitivity in Gaussian sigma alpha1 = alpha + 1.0 sum = Non + Noff arg1 = Non / sum @@ -230,23 +230,67 @@ def show_one_sensitivity(rsp, name, color='r', duration=180000.0, alpha=0.2, logE = [-1.7+0.2*i for i in range(nbins)] E = [math.pow(10.0, e) for e in logE] - # Generate sensitivity vector + # Initialise sensitivity vector flux = [0.0 for i in range(nbins)] + + # Loop over all energy bins for i in range(nbins): - ewidth = (math.pow(10.0, logE[i]+0.1) - \ - math.pow(10.0, logE[i]-0.1))*1.0e6 - aeff = rsp.aeff()(logE[i]) - r68 = rsp.psf().delta_max(logE[i]) / 5.0 / r68_to_sigma + + # 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)) - bgd_counts = rsp.background()(logE[i], 0.0, 0.0) * duration * ewidth * solidangle - Noff = bgd_counts/alpha + + # 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: - Non = Non_lima(sigma, Noff) + + # 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') -- libgit2 0.21.4