Michael Mayer / gammalib

Commit d34270a7382e2b02b68969f232fda2d2ae5bc90f

Authored by Jürgen Knödlseder
1 parent 33d830a2

Add test script for diffuse map cubes.

Inline Side-by-side
Showing 1 changed file with 188 additions and 0 deletions
test/test_diffuse_cube.py 0 → 100755
  View file @d34270a
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script tests the simulation of diffuse map cubes.
  4 +#
  5 +# If matplotlib is installed, the simulation results will be displayed.
  6 +#
  7 +# --------------------------------------------------------------------------
  8 +#
  9 +# Copyright (C) 2013 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 +from gammalib import *
  26 +from math import *
  27 +
  28 +
  29 +# ======================== #
  30 +# Simulate CTA observation #
  31 +# ======================== #
  32 +def simulate(xmlname, e_min, e_max, area, duration, dir, radius):
  33 + """
  34 + Simulate CTA observation.
  35 + """
  36 + # Allocate MC parameters
  37 + emin = GEnergy()
  38 + emax = GEnergy()
  39 + tmin = GTime(0.0)
  40 + tmax = GTime(duration)
  41 +
  42 + # Define MC parameters
  43 + emin.TeV(e_min)
  44 + emax.TeV(e_max)
  45 +
  46 + # Allocate random number generator
  47 + ran = GRan()
  48 +
  49 + # Load models and extract first model
  50 + models = GModels(xmlname)
  51 + model = models[0]
  52 + #print(model)
  53 + #print(model.spatial())
  54 + #print(model.spatial().spectrum())
  55 +
  56 + # Simulate photons
  57 + photons = model.mc(area, dir, radius, emin, emax, tmin, tmax, ran)
  58 +
  59 + # Print photons
  60 + print(str(len(photons)) + " photons simulated.")
  61 +
  62 + # Return photons
  63 + return photons
  64 +
  65 +
  66 +# ============ #
  67 +# Show photons #
  68 +# ============ #
  69 +def show_photons(photons, xmlname, e_min, e_max, area, duration, dir, radius, ebins=30):
  70 + """
  71 + Show photons using matplotlib (if available).
  72 + """
  73 + # Only proceed if matplotlib is available
  74 + try:
  75 + # Import matplotlib
  76 + import matplotlib.pyplot as plt
  77 +
  78 + # Create figure
  79 + plt.figure(1)
  80 + plt.title("MC simulated photon spectrum (" + str(e_min) + '-' + str(e_max) + " TeV)")
  81 +
  82 + # Setup energy range covered by data
  83 + ebds = GEbounds(ebins, GEnergy(e_min, "TeV"), GEnergy(e_max, "TeV"))
  84 +
  85 + # Create energy axis
  86 + energy = []
  87 + for i in range(ebds.size()):
  88 + energy.append(ebds.elogmean(i).TeV())
  89 +
  90 + # Fill histogram
  91 + counts = [0.0 for i in range(ebds.size())]
  92 + for photon in photons:
  93 + index = ebds.index(photon.energy())
  94 + counts[index] = counts[index] + 1.0
  95 +
  96 + # Create error bars
  97 + error = [sqrt(c) for c in counts]
  98 +
  99 + # Get model values
  100 + models = GModels(xmlname)
  101 + mod = models[0]
  102 + model = []
  103 + t = GTime()
  104 + mod.spatial().set_mc_cone(dir, radius)
  105 + for i in range(ebds.size()):
  106 + eng = ebds.elogmean(i)
  107 + ewidth = ebds.ewidth(i)
  108 + f = mod.spatial().spectrum().eval(eng, t) * \
  109 + area * duration * ewidth.MeV()
  110 + model.append(f)
  111 +
  112 + # Plot data
  113 + plt.loglog(energy, counts, 'ro')
  114 + plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  115 +
  116 + # Plot model
  117 + plt.plot(energy, model, 'b-')
  118 +
  119 + # Set axes
  120 + plt.xlabel("Energy (TeV)")
  121 + plt.ylabel("Number of incident photons")
  122 +
  123 + # Create figure
  124 + plt.figure(2)
  125 + plt.title("MC simulated photon map")
  126 +
  127 + # Create RA and DEC arrays
  128 + ra = []
  129 + dec = []
  130 + for photon in photons:
  131 + ra.append(photon.dir().ra_deg())
  132 + dec.append(photon.dir().dec_deg())
  133 +
  134 + # Make scatter plot
  135 + plt.scatter(ra, dec, marker=".")
  136 +
  137 + # Set axes
  138 + plt.xlabel("Right Ascension (deg)")
  139 + plt.ylabel("Declination (deg)")
  140 +
  141 + # Notify
  142 + print("PLEASE CLOSE WINDOW TO CONTINUE ...")
  143 +
  144 + # Allocate histogram
  145 + # Show plot
  146 + plt.show()
  147 +
  148 + except ImportError:
  149 + print("Matplotlib is not (correctly) installed on your system.")
  150 +
  151 + # Return
  152 + return
  153 +
  154 +
  155 +#==========================#
  156 +# Main routine entry point #
  157 +#==========================#
  158 +if __name__ == '__main__':
  159 + """
  160 + Simulate photons.
  161 + """
  162 + # Dump header
  163 + print("")
  164 + print("********************")
  165 + print("* Simulate photons *")
  166 + print("********************")
  167 +
  168 + # Set model filename
  169 + xmlnames = ["data/model_diffuse_cube.xml"]
  170 +
  171 + # Set simulation parameters
  172 + e_min = 0.1 # 0.1 TeV
  173 + e_max = 100.0 # 100 TeV
  174 + area = 1.0 # 1 cm2
  175 + duration = 1.0 # 1 sec
  176 + dir = GSkyDir()
  177 + dir.radec_deg(84.17263, 22.01444)
  178 + radius = 1.0
  179 +
  180 +
  181 + # Loop over models
  182 + for xmlname in xmlnames:
  183 +
  184 + # Perform simulation
  185 + photons = simulate(xmlname, e_min, e_max, area, duration, dir, radius)
  186 +
  187 + # Show photons
  188 + show_photons(photons, xmlname, e_min, e_max, area, duration, dir, radius)