Justine Devin / gammalib

Commit c65243a6ac1e9f9bd60d7f2a954827209123f1eb

Authored by Jürgen Knödlseder
1 parent 7aee2824

Add PSF simulation test script.

Inline Side-by-side
Showing 1 changed file with 167 additions and 0 deletions
inst/cta/test/example_sim_psf.py 0 → 100755
  View file @c65243a
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script simulated the PSF distribution.
  4 +#
  5 +# Based on the MAGIC spectrum of the Crab nebula, and by assuming a powerlaw,
  6 +# it will create a Monte Carlo sample of photons.
  7 +#
  8 +# If matplotlib is installed, the spectrum will be displayed on the screen.
  9 +#
  10 +# --------------------------------------------------------------------------
  11 +#
  12 +# Copyright (C) 2013 Juergen Knoedlseder
  13 +#
  14 +# This program is free software: you can redistribute it and/or modify
  15 +# it under the terms of the GNU General Public License as published by
  16 +# the Free Software Foundation, either version 3 of the License, or
  17 +# (at your option) any later version.
  18 +#
  19 +# This program is distributed in the hope that it will be useful,
  20 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  22 +# GNU General Public License for more details.
  23 +#
  24 +# You should have received a copy of the GNU General Public License
  25 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  26 +#
  27 +# ==========================================================================
  28 +from gammalib import *
  29 +from math import *
  30 +
  31 +
  32 +# ======================== #
  33 +# Simulate CTA observation #
  34 +# ======================== #
  35 +def simulate(xmlname, e_min, e_max, area, duration):
  36 + """
  37 + Simulate CTA observation.
  38 + """
  39 + # Allocate MC parameters
  40 + dir = GSkyDir()
  41 + emin = GEnergy()
  42 + emax = GEnergy()
  43 + tmin = GTime(0.0)
  44 + tmax = GTime(duration)
  45 +
  46 + # Define MC parameters
  47 + dir.radec_deg(83.6331, 22.0145)
  48 + radius = 10.0
  49 + emin.TeV(e_min)
  50 + emax.TeV(e_max)
  51 +
  52 + # Allocate random number generator
  53 + ran = GRan()
  54 +
  55 + # Load models and extract first model
  56 + models = GModels(xmlname)
  57 + model = models[0]
  58 + print(model)
  59 +
  60 + # Simulate photons
  61 + photons = model.mc(area, dir, radius, emin, emax, tmin, tmax, ran)
  62 +
  63 + # Print photons
  64 + print(str(len(photons)) + " photons simulated.")
  65 +
  66 + # Return photons
  67 + return photons
  68 +
  69 +
  70 +# ======== #
  71 +# Show PSF #
  72 +# ======== #
  73 +def sim_psf(response, energy, r_max=0.8, rbins=1000, nmc=1000000):
  74 + """
  75 + Simulate PSF and show results using matplotlib (if available).
  76 + """
  77 + # Check if matplotlib is available
  78 + try:
  79 + import matplotlib.pyplot as plt
  80 + has_matplotlib = True
  81 + except ImportError:
  82 + print("Matplotlib is not (correctly) installed on your system.")
  83 + has_matplotlib = False
  84 +
  85 + # Continue only if matplotlib is available
  86 + if has_matplotlib:
  87 +
  88 + # Create figure
  89 + plt.figure(1)
  90 + plt.title("MC simulated PSF ("+str(energy)+" TeV)")
  91 +
  92 + # Set log10(energy)
  93 + logE = log10(energy)
  94 +
  95 + # Create delta axis
  96 + delta = []
  97 + dr = r_max/rbins
  98 + for i in range(rbins):
  99 + delta.append((i+0.5)*dr)
  100 +
  101 + # Reset histogram
  102 + counts = [0.0 for i in range(rbins)]
  103 +
  104 + # Allocate random number generator
  105 + ran = GRan()
  106 +
  107 + # Simulate offsets
  108 + for i in range(nmc):
  109 + offset = rsp.psf().mc(ran, logE) * 180.0/pi
  110 + #print(offset)
  111 + index = int(offset/dr)
  112 + if (index < rbins):
  113 + counts[index] += 1.0
  114 +
  115 + # Create error bars
  116 + error = [sqrt(c) for c in counts]
  117 +
  118 + # Get expected PSF
  119 + sum = 0.0
  120 + psf = []
  121 + for i in range(rbins):
  122 + r = delta[i]*pi/180.0
  123 + value = rsp.psf()(r, logE) * 2.0*pi * sin(r) * dr * pi/180.0 * nmc
  124 + sum += value
  125 + psf.append(value)
  126 + print(sum)
  127 +
  128 + # Plot simulated data
  129 + plt.plot(delta, counts, 'ro')
  130 + plt.errorbar(delta, counts, error, fmt=None, ecolor='r')
  131 +
  132 + # Plot PSF
  133 + plt.plot(delta, psf, 'b-')
  134 +
  135 + # Set axes
  136 + plt.xlabel("Offset angle (degrees)")
  137 + plt.ylabel("Number of counts")
  138 +
  139 + # Notify
  140 + print("PLEASE CLOSE WINDOW TO CONTINUE ...")
  141 +
  142 + # Show plot
  143 + plt.show()
  144 +
  145 + # Return
  146 + return
  147 +
  148 +
  149 +#==========================#
  150 +# Main routine entry point #
  151 +#==========================#
  152 +if __name__ == '__main__':
  153 + """
  154 + Simulate PSF.
  155 + """
  156 + # Dump header
  157 + print("")
  158 + print("****************")
  159 + print("* Simulate PSF *")
  160 + print("****************")
  161 +
  162 + # Load response
  163 + rsp = GCTAResponse("irf_test.fits",
  164 + "../caldb/data/cta/e/bcf/000002")
  165 +
  166 + # Simulate PSF
  167 + sim_psf(rsp, 1.0)
... ...