jean perrot / ctools

Commit 39a5aa700397234982b228c1ed353e80986e50f5

Authored by Jürgen Knödlseder
1 parent 9c5ea323

Add event simulation example script and test XML file

Inline Side-by-side
Showing 2 changed files with 281 additions and 0 deletions
examples/model_bgd.xml 0 → 100644
  View file @39a5aa7
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Background" type="RadialAcceptance" instrument="CTA">
  4 + <spectrum type="FileFunction" file="$CTOOLS/share/models/bkg_dummy.txt">
  5 + <parameter scale="1.0" name="Normalization" min="0.0" max="1000.0" value="1.0" free="1"/>
  6 + </spectrum>
  7 + <radialModel type="Gaussian">
  8 + <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
  9 + </radialModel>
  10 + </source>
  11 +</source_library>
  12 +
... ...
examples/simulate_events.py 0 → 100755
  View file @39a5aa7
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script simulates events from a data-space model and compares them to
  4 +# the model.
  5 +#
  6 +# If matplotlib is installed, the event spectrum will be displayed.
  7 +#
  8 +# --------------------------------------------------------------------------
  9 +#
  10 +# Copyright (C) 2014 Juergen Knoedlseder
  11 +#
  12 +# This program is free software: you can redistribute it and/or modify
  13 +# it under the terms of the GNU General Public License as published by
  14 +# the Free Software Foundation, either version 3 of the License, or
  15 +# (at your option) any later version.
  16 +#
  17 +# This program is distributed in the hope that it will be useful,
  18 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  20 +# GNU General Public License for more details.
  21 +#
  22 +# You should have received a copy of the GNU General Public License
  23 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  24 +#
  25 +# ==========================================================================
  26 +import gammalib
  27 +import ctools
  28 +import math
  29 +try:
  30 + import matplotlib.pyplot as plt
  31 + has_matplotlib = True
  32 +except:
  33 + has_matplotlib = False
  34 +
  35 +
  36 +# ======================== #
  37 +# Create CTA observation #
  38 +# ======================== #
  39 +def createobs(ra=86.171648, dec=-1.4774586, rad=5.0, \
  40 + emin=0.1, emax=100.0, duration=360000.0, deadc=0.95, \
  41 + irf="cta_dummy_irf", caldb="$CTOOLS/share/caldb/data/cta/dummy"):
  42 + """
  43 + Create CTA observation.
  44 + """
  45 + # Allocate CTA observation
  46 + obs = gammalib.GCTAObservation()
  47 +
  48 + # Set calibration database
  49 + db = gammalib.GCaldb()
  50 + if (gammalib.dir_exists(caldb)):
  51 + db.rootdir(caldb)
  52 + else:
  53 + db.open("cta", caldb)
  54 +
  55 + # Set pointing direction
  56 + pntdir = gammalib.GSkyDir()
  57 + pntdir.radec_deg(ra, dec)
  58 + pnt = gammalib.GCTAPointing()
  59 + pnt.dir(pntdir)
  60 + obs.pointing(pnt)
  61 +
  62 + # Set ROI
  63 + roi = gammalib.GCTARoi()
  64 + instdir = gammalib.GCTAInstDir()
  65 + instdir.dir(pntdir)
  66 + roi.centre(instdir)
  67 + roi.radius(rad)
  68 +
  69 + # Set GTI
  70 + gti = gammalib.GGti()
  71 + start = gammalib.GTime(0.0)
  72 + stop = gammalib.GTime(duration)
  73 + gti.append(start, stop)
  74 +
  75 + # Set energy boundaries
  76 + ebounds = gammalib.GEbounds()
  77 + e_min = gammalib.GEnergy()
  78 + e_max = gammalib.GEnergy()
  79 + e_min.TeV(emin)
  80 + e_max.TeV(emax)
  81 + ebounds.append(e_min, e_max)
  82 +
  83 + # Allocate event list
  84 + events = gammalib.GCTAEventList()
  85 + events.roi(roi)
  86 + events.gti(gti)
  87 + events.ebounds(ebounds)
  88 + obs.events(events)
  89 +
  90 + # Set instrument response
  91 + obs.response(irf, db)
  92 +
  93 + # Set ontime, livetime, and deadtime correction factor
  94 + obs.ontime(duration)
  95 + obs.livetime(duration*deadc)
  96 + obs.deadc(deadc)
  97 +
  98 + # Return observation
  99 + return obs
  100 +
  101 +
  102 +# =============== #
  103 +# Simulate events #
  104 +# =============== #
  105 +def simulate(obs, xmlname):
  106 + """
  107 + Simulate events.
  108 + """
  109 + # Allocate random number generator
  110 + ran = gammalib.GRan()
  111 +
  112 + # Load models and extract first model
  113 + models = gammalib.GModels(xmlname)
  114 + model = models[0]
  115 + print(model)
  116 +
  117 + # Simulate events
  118 + events = model.mc(obs, ran)
  119 +
  120 + # Print event statistics
  121 + npred = obs.npred(models)
  122 + print(str(len(events)) + " events simulated.")
  123 + print(str(npred) + " events expected from Npred.")
  124 +
  125 + # Check distance
  126 + pntdir = obs.pointing().dir()
  127 + for event in events:
  128 + distance = gammalib.GCTAInstDir(event.dir()).dir().dist_deg(pntdir)
  129 + if (distance > 5.0):
  130 + print("Too far: "+str(distance))
  131 +
  132 + # Return events
  133 + return events
  134 +
  135 +
  136 +# ============================== #
  137 +# Simulate events using ctobssim #
  138 +# ============================== #
  139 +def simulate_ctobssim(obs, xmlname, seed=0):
  140 + """
  141 + Simulate events using ctobssim.
  142 + """
  143 + # Create containers
  144 + observations = gammalib.GObservations()
  145 + observations.append(obs)
  146 +
  147 + # Append models
  148 + models = gammalib.GModels(xmlname)
  149 + observations.models(models)
  150 + print(models[0])
  151 +
  152 + # Allocate ctobssim application and set parameters
  153 + sim = ctools.ctobssim(observations)
  154 + sim['seed'].integer(seed)
  155 +
  156 + # Run simulator
  157 + sim.run()
  158 +
  159 + # Retrieve events
  160 + events = sim.obs()[0].events().copy()
  161 +
  162 + # Print event statistics
  163 + npred = obs.npred(models)
  164 + print(str(len(events)) + " events simulated.")
  165 + print(str(npred) + " events expected from Npred.")
  166 +
  167 + # Delete the simulation
  168 + del sim
  169 +
  170 + # Return events
  171 + return events
  172 +
  173 +
  174 +# =========== #
  175 +# Show events #
  176 +# =========== #
  177 +def show_events(events, xmlname, duration, emin, emax, ebins=30):
  178 + """
  179 + Show events using matplotlib.
  180 + """
  181 + # Create figure
  182 + plt.figure(1)
  183 + plt.title("MC simulated event spectrum (" + str(emin) + '-' + str(emax) + " TeV)")
  184 +
  185 + # Setup energy range covered by data
  186 + ebds = gammalib.GEbounds(ebins, gammalib.GEnergy(emin, "TeV"),
  187 + gammalib.GEnergy(emax, "TeV"))
  188 +
  189 + # Create energy axis
  190 + energy = []
  191 + for i in range(ebds.size()):
  192 + energy.append(ebds.elogmean(i).TeV())
  193 +
  194 + # Fill histogram
  195 + counts = [0.0 for i in range(ebds.size())]
  196 + for event in events:
  197 + index = ebds.index(event.energy())
  198 + counts[index] = counts[index] + 1.0
  199 +
  200 + # Create error bars
  201 + error = [math.sqrt(c) for c in counts]
  202 +
  203 + # Get model values
  204 + sum = 0.0
  205 + models = gammalib.GModels(xmlname)
  206 + m = models[0]
  207 + model = []
  208 + t = gammalib.GTime()
  209 + for i in range(ebds.size()):
  210 + eval = ebds.elogmean(i)
  211 + ewidth = ebds.emax(i) - ebds.emin(i)
  212 + f = m.npred(eval, t, obs) * ewidth.MeV() * duration
  213 + sum += f
  214 + model.append(f)
  215 + print(str(sum) + " events expected from spectrum (integration).")
  216 +
  217 + # Plot data
  218 + plt.loglog(energy, counts, 'ro')
  219 + plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  220 +
  221 + # Plot model
  222 + plt.plot(energy, model, 'b-')
  223 +
  224 + # Set axes
  225 + plt.xlabel("Energy (TeV)")
  226 + plt.ylabel("Number of events")
  227 +
  228 + # Notify
  229 + print("PLEASE CLOSE WINDOW TO CONTINUE ...")
  230 +
  231 + # Show plot
  232 + plt.show()
  233 +
  234 + # Return
  235 + return
  236 +
  237 +
  238 +#==========================#
  239 +# Main routine entry point #
  240 +#==========================#
  241 +if __name__ == '__main__':
  242 + """
  243 + Simulate events.
  244 + """
  245 + # Dump header
  246 + print("")
  247 + print("*******************")
  248 + print("* Simulate events *")
  249 + print("*******************")
  250 +
  251 + # Set XML names
  252 + xmlname = "model_bgd.xml"
  253 +
  254 + # Set simulation parameters
  255 + duration = 360000.0 # Observation duration (s)
  256 + emin = 0.1 # 0.1 TeV
  257 + emax = 100.0 # 100 TeV
  258 + ebins = 30
  259 +
  260 + # Set observation
  261 + obs = createobs(emin=emin, emax=emax, duration=duration)
  262 +
  263 + # Perform simulation
  264 + #events = simulate(obs, xmlname)
  265 + events = simulate_ctobssim(obs, xmlname)
  266 +
  267 + # Show events
  268 + if has_matplotlib:
  269 + show_events(events, xmlname, duration, emin, emax, ebins=ebins)
... ...