diff --git a/examples/model_bgd.xml b/examples/model_bgd.xml
new file mode 100644
index 0000000..b927387
--- /dev/null
+++ b/examples/model_bgd.xml
@@ -0,0 +1,12 @@
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/simulate_events.py b/examples/simulate_events.py
new file mode 100755
index 0000000..2f5b919
--- /dev/null
+++ b/examples/simulate_events.py
@@ -0,0 +1,269 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This script simulates events from a data-space model and compares them to
+# the model.
+#
+# If matplotlib is installed, the event spectrum will be displayed.
+#
+# --------------------------------------------------------------------------
+#
+# Copyright (C) 2014 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 gammalib
+import ctools
+import math
+try:
+ import matplotlib.pyplot as plt
+ has_matplotlib = True
+except:
+ has_matplotlib = False
+
+
+# ======================== #
+# Create CTA observation #
+# ======================== #
+def createobs(ra=86.171648, dec=-1.4774586, rad=5.0, \
+ emin=0.1, emax=100.0, duration=360000.0, deadc=0.95, \
+ irf="cta_dummy_irf", caldb="$CTOOLS/share/caldb/data/cta/dummy"):
+ """
+ Create CTA observation.
+ """
+ # Allocate CTA observation
+ obs = gammalib.GCTAObservation()
+
+ # Set calibration database
+ db = gammalib.GCaldb()
+ if (gammalib.dir_exists(caldb)):
+ db.rootdir(caldb)
+ else:
+ db.open("cta", caldb)
+
+ # Set pointing direction
+ pntdir = gammalib.GSkyDir()
+ pntdir.radec_deg(ra, dec)
+ pnt = gammalib.GCTAPointing()
+ pnt.dir(pntdir)
+ obs.pointing(pnt)
+
+ # Set ROI
+ roi = gammalib.GCTARoi()
+ instdir = gammalib.GCTAInstDir()
+ instdir.dir(pntdir)
+ roi.centre(instdir)
+ roi.radius(rad)
+
+ # Set GTI
+ gti = gammalib.GGti()
+ start = gammalib.GTime(0.0)
+ stop = gammalib.GTime(duration)
+ gti.append(start, stop)
+
+ # Set energy boundaries
+ ebounds = gammalib.GEbounds()
+ e_min = gammalib.GEnergy()
+ e_max = gammalib.GEnergy()
+ e_min.TeV(emin)
+ e_max.TeV(emax)
+ ebounds.append(e_min, e_max)
+
+ # Allocate event list
+ events = gammalib.GCTAEventList()
+ events.roi(roi)
+ events.gti(gti)
+ events.ebounds(ebounds)
+ obs.events(events)
+
+ # Set instrument response
+ obs.response(irf, db)
+
+ # Set ontime, livetime, and deadtime correction factor
+ obs.ontime(duration)
+ obs.livetime(duration*deadc)
+ obs.deadc(deadc)
+
+ # Return observation
+ return obs
+
+
+# =============== #
+# Simulate events #
+# =============== #
+def simulate(obs, xmlname):
+ """
+ Simulate events.
+ """
+ # Allocate random number generator
+ ran = gammalib.GRan()
+
+ # Load models and extract first model
+ models = gammalib.GModels(xmlname)
+ model = models[0]
+ print(model)
+
+ # Simulate events
+ events = model.mc(obs, ran)
+
+ # Print event statistics
+ npred = obs.npred(models)
+ print(str(len(events)) + " events simulated.")
+ print(str(npred) + " events expected from Npred.")
+
+ # Check distance
+ pntdir = obs.pointing().dir()
+ for event in events:
+ distance = gammalib.GCTAInstDir(event.dir()).dir().dist_deg(pntdir)
+ if (distance > 5.0):
+ print("Too far: "+str(distance))
+
+ # Return events
+ return events
+
+
+# ============================== #
+# Simulate events using ctobssim #
+# ============================== #
+def simulate_ctobssim(obs, xmlname, seed=0):
+ """
+ Simulate events using ctobssim.
+ """
+ # Create containers
+ observations = gammalib.GObservations()
+ observations.append(obs)
+
+ # Append models
+ models = gammalib.GModels(xmlname)
+ observations.models(models)
+ print(models[0])
+
+ # Allocate ctobssim application and set parameters
+ sim = ctools.ctobssim(observations)
+ sim['seed'].integer(seed)
+
+ # Run simulator
+ sim.run()
+
+ # Retrieve events
+ events = sim.obs()[0].events().copy()
+
+ # Print event statistics
+ npred = obs.npred(models)
+ print(str(len(events)) + " events simulated.")
+ print(str(npred) + " events expected from Npred.")
+
+ # Delete the simulation
+ del sim
+
+ # Return events
+ return events
+
+
+# =========== #
+# Show events #
+# =========== #
+def show_events(events, xmlname, duration, emin, emax, ebins=30):
+ """
+ Show events using matplotlib.
+ """
+ # Create figure
+ plt.figure(1)
+ plt.title("MC simulated event spectrum (" + str(emin) + '-' + str(emax) + " TeV)")
+
+ # Setup energy range covered by data
+ ebds = gammalib.GEbounds(ebins, gammalib.GEnergy(emin, "TeV"),
+ gammalib.GEnergy(emax, "TeV"))
+
+ # Create energy axis
+ energy = []
+ for i in range(ebds.size()):
+ energy.append(ebds.elogmean(i).TeV())
+
+ # Fill histogram
+ counts = [0.0 for i in range(ebds.size())]
+ for event in events:
+ index = ebds.index(event.energy())
+ counts[index] = counts[index] + 1.0
+
+ # Create error bars
+ error = [math.sqrt(c) for c in counts]
+
+ # Get model values
+ sum = 0.0
+ models = gammalib.GModels(xmlname)
+ m = models[0]
+ model = []
+ t = gammalib.GTime()
+ for i in range(ebds.size()):
+ eval = ebds.elogmean(i)
+ ewidth = ebds.emax(i) - ebds.emin(i)
+ f = m.npred(eval, t, obs) * ewidth.MeV() * duration
+ sum += f
+ model.append(f)
+ print(str(sum) + " events expected from spectrum (integration).")
+
+ # Plot data
+ plt.loglog(energy, counts, 'ro')
+ plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
+
+ # Plot model
+ plt.plot(energy, model, 'b-')
+
+ # Set axes
+ plt.xlabel("Energy (TeV)")
+ plt.ylabel("Number of events")
+
+ # Notify
+ print("PLEASE CLOSE WINDOW TO CONTINUE ...")
+
+ # Show plot
+ plt.show()
+
+ # Return
+ return
+
+
+#==========================#
+# Main routine entry point #
+#==========================#
+if __name__ == '__main__':
+ """
+ Simulate events.
+ """
+ # Dump header
+ print("")
+ print("*******************")
+ print("* Simulate events *")
+ print("*******************")
+
+ # Set XML names
+ xmlname = "model_bgd.xml"
+
+ # Set simulation parameters
+ duration = 360000.0 # Observation duration (s)
+ emin = 0.1 # 0.1 TeV
+ emax = 100.0 # 100 TeV
+ ebins = 30
+
+ # Set observation
+ obs = createobs(emin=emin, emax=emax, duration=duration)
+
+ # Perform simulation
+ #events = simulate(obs, xmlname)
+ events = simulate_ctobssim(obs, xmlname)
+
+ # Show events
+ if has_matplotlib:
+ show_events(events, xmlname, duration, emin, emax, ebins=ebins)