diff --git a/examples/make_unbinned_analysis.py b/examples/make_unbinned_analysis.py
new file mode 100755
index 0000000..cd389d7
--- /dev/null
+++ b/examples/make_unbinned_analysis.py
@@ -0,0 +1,427 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This script illustrates how to build an unbinned analysis pipeline with
+# the ctatools. If matplotlib is installed, the spectrum will be displayed
+# on the screen.
+#
+# Required 3rd party modules:
+# None
+#
+# Copyright (C) 2011 Jurgen Knodlseder
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+from ctatools import *
+from gammalib import *
+from math import *
+import os
+import glob
+import sys
+
+
+# ============================= #
+# Analysis pipeline - version 1 #
+# ============================= #
+def pipeline_v1():
+	"""
+	Unbinned analysis pipeline - save intermediate results.
+    
+    This function implements an analysis pipeline that successively calls
+    ctobssim, ctselect and ctlike by saving the intermediate results as FITS
+    files on disk. This replicates the way how the analysis would be done
+    in a ftools-like approach.
+	"""
+	# Set script parameters
+	model_name           = "${GAMMALIB}/share/models/crab.xml"
+	events_name          = "events.fits"
+	selected_events_name = "selected_events.fits"
+	result_name          = "results.xml"
+	caldb                = "${GAMMALIB}/share/caldb/cta"
+	irf                  = "kb_E_50h_v3"
+	ra                   =   83.63
+	dec                  =   22.01
+	rad_sim              =   10.0
+	tstart               =    0.0
+	tstop                = 1800.0
+	emin                 =    0.1
+	emax                 =  100.0
+	rad_select           =    3.0
+
+	# Initialise timing
+	wall_seconds = 0.0
+	cpu_seconds  = 0.0
+
+	# Simulate events
+	sim = ctobssim()
+	sim.logFileOpen()
+	sim["infile"].filename(model_name)
+	sim["outfile"].filename(events_name)
+	sim["caldb"].string(caldb)
+	sim["irf"].string(irf)
+	sim["ra"].real(ra)
+	sim["dec"].real(dec)
+	sim["rad"].real(rad_sim)
+	sim["tmin"].real(tstart)
+	sim["tmax"].real(tstop)
+	sim["emin"].real(emin)
+	sim["emax"].real(emax)
+	sim.execute()
+	print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
+
+	# Update timing
+	wall_seconds += sim.telapse()
+	cpu_seconds  += sim.celapse()
+	
+	# Select events
+	select = ctselect()
+	select.logFileOpen()
+	select["infile"].filename(events_name)
+	select["outfile"].filename(selected_events_name)
+	select["ra"].real(ra)
+	select["dec"].real(dec)
+	select["rad"].real(rad_select)
+	select["tmin"].real(tstart)
+	select["tmax"].real(tstop)
+	select["emin"].real(emin)
+	select["emax"].real(emax)
+	select.execute()
+	print "Select event ("+str(select.celapse())+" CPU seconds)"
+
+	# Update timing
+	wall_seconds += select.telapse()
+	cpu_seconds  += select.celapse()
+	
+	# Perform maximum likelihood fitting
+	like = ctlike()
+	like.logFileOpen()
+	like["evfile"].filename(selected_events_name)
+	like["srcmdl"].filename(model_name)
+	like["outmdl"].filename(result_name)
+	like["method"].string("UNBINNED")
+	like["caldb"].string(caldb)
+	like["irf"].string(irf)
+	like.execute()
+	print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
+
+	# Update timing
+	wall_seconds += like.telapse()
+	cpu_seconds  += like.celapse()
+	
+	# Show total times
+	print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
+	print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
+
+	# Return
+	return
+
+
+# ============================= #
+# Analysis pipeline - version 2 #
+# ============================= #
+def pipeline_v2(show_data):
+	"""
+	Unbinned analysis pipeline - keep intermediate results in memory.
+    
+    This function implements an analysis pipeline that successively calls
+    ctobssim, ctselect and ctlike without saving the intermediate results as
+	FITS files on disk. All data is only hold in memory.
+	"""
+	# Set script parameters
+	model_name  = "${GAMMALIB}/share/models/crab.xml"
+	caldb       = "${GAMMALIB}/share/caldb/cta"
+	irf         = "kb_E_50h_v3"
+	ra          =   83.63
+	dec         =   22.01
+	rad_sim     =   10.0
+	tstart      =    0.0
+	tstop       = 1800.0
+	emin        =    0.1
+	emax        =  100.0
+	rad_select  =    3.0
+
+	# Initialise timing
+	wall_seconds = 0.0
+	cpu_seconds  = 0.0
+
+	# Simulate events
+	sim = ctobssim()
+	sim["infile"].filename(model_name)
+	sim["caldb"].string(caldb)
+	sim["irf"].string(irf)
+	sim["ra"].real(ra)
+	sim["dec"].real(dec)
+	sim["rad"].real(rad_sim)
+	sim["tmin"].real(tstart)
+	sim["tmax"].real(tstop)
+	sim["emin"].real(emin)
+	sim["emax"].real(emax)
+	sim.run()
+	print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
+
+	# Update timing
+	wall_seconds += sim.telapse()
+	cpu_seconds  += sim.celapse()
+
+	# Select events
+	select = ctselect(sim.obs())
+	select["ra"].real(ra)
+	select["dec"].real(dec)
+	select["rad"].real(rad_select)
+	select["tmin"].real(tstart)
+	select["tmax"].real(tstop)
+	select["emin"].real(emin)
+	select["emax"].real(emax)
+	select.run()
+	print "Select event ("+str(select.celapse())+" CPU seconds)"
+
+	# Update timing
+	wall_seconds += select.telapse()
+	cpu_seconds  += select.celapse()
+
+	# Perform maximum likelihood fitting
+	like = ctlike(select.obs())
+	like.run()
+	print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
+
+	# Update timing
+	wall_seconds += like.telapse()
+	cpu_seconds  += like.celapse()
+		
+	# Show total times
+	print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
+	print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
+
+	# Show model fitting results
+	#print like.obs().models()
+	
+	# Optionally plot counts
+	if show_data:
+		plot_counts(like.obs())
+	
+	# Return
+	return
+
+
+# =========== #
+# Plot counts #
+# =========== #
+def plot_counts(observations):
+	"""
+	Plot counts.
+	"""
+	# Only proceed if matplotlib is available
+	try:
+		# Import matplotlib
+		import matplotlib.pyplot as plt
+
+		# Set legend fontsize
+		params = {'legend.fontsize': 10}
+		plt.rcParams.update(params)
+
+		# Set plot styles
+		styles = ['b-', 'g-', 'y-', 'n-']
+
+		# Dump header
+		print ""
+		print "Make plots (using matplotlib):"
+		print "=============================="
+		
+		# Create figure 1
+		plt.figure(1,figsize=(12,6))
+		plt.subplots_adjust(hspace=.7)
+
+		# Create subplot 1
+		plt.subplot(121)
+		plt.title("Spectrum (summed over all pixels)")
+		
+		# Loop over observations
+		for obs in observations:
+		
+			# Get event list
+			list = cast_GCTAEventList(obs.events())
+			
+			# Create energy axis
+			ebounds = GEbounds()
+			emin    = GEnergy()
+			emax    = GEnergy()
+			emin.TeV(0.1)
+			emax.TeV(100.0)
+			ebounds.setlog(emin, emax, 10)
+			energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
+
+			# Create spectrum
+			print "Extract data:"
+			counts = [0.0 for i in range(ebounds.size())]
+			print list.size()
+			for atom in list:
+				index         = ebounds.index(atom.energy())
+				counts[index] = counts[index] + 1.0
+
+			# Create error bars
+			error = [sqrt(c) for c in counts]
+
+			# Plot spectrum
+			plt.loglog(energy, counts, 'ro', label='data')
+			#plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
+
+			# Extract models
+			print "Extract models:"
+			sum_model = [0.0 for i in range(ebounds.size())]
+			for k, m in enumerate(observations.models()):
+				print "- "+m.name()
+				model = [0.0 for i in range(ebounds.size())]
+				for atom in list:
+					index        = ebounds.index(atom.energy())
+					prob         = m.eval(atom, obs)
+					model[index] = model[index] + prob * atom.size()
+				for i in range(ebounds.size()):
+					sum_model[i] = sum_model[i] + model[i]
+				#plt.loglog(energy, model, styles[k], label=m.name())
+			#plt.loglog(energy, sum_model, 'r-', label='total')
+
+		# Put labels
+		plt.xlabel("Energy (TeV)")
+		plt.ylabel("Counts")
+		plt.legend(loc="lower left")
+
+		# Create subplot 2
+		plt.subplot(122)
+		plt.title("Offset (summed over all energies)")
+
+		# Set Crab direction
+		crab = GSkyDir()
+		crab.radec_deg(83.63, 22.01)
+		
+		# Loop over observations
+		for obs in observations:
+		
+			# Get event cube
+			list = cast_GCTAEventList(obs.events())
+			
+			# Create offset histogram
+			print "Extract data:"
+			nx       = 30
+			doffset2 = 0.01
+			offset2  = [(i+0.5)*doffset2 for i in range(nx)]
+			counts   = [0.0  for i in range(nx)]
+			for atom in list:
+				off   = atom.dir().dist_deg(crab)
+				off2  = off*off
+				index = int(off2/doffset2)
+				if index < nx:
+					counts[index] = counts[index] + 1.0
+
+			# Create error bars
+			error = [sqrt(c) for c in counts]
+
+			# Plot distribution
+			plt.semilogy(offset2, counts, 'ro', label='data')
+
+			# Extract models
+			print "Extract models:"
+			sum_model = [0.0 for i in range(nx)]
+			for k, m in enumerate(observations.models()):
+				print "- "+m.name()
+				model = [0.0 for i in range(nx)]
+				for atom in list:
+					off   = atom.dir().dist_deg(crab)
+					off2  = off*off
+					index = int(off2/doffset2)
+					if index < nx:
+						prob         = m.eval(atom, obs)
+						model[index] = model[index] + prob * atom.size()
+				for i in range(nx):
+					sum_model[i] = sum_model[i] + model[i]
+				#plt.plot(offset2, model, styles[k], label=m.name())
+			#plt.plot(offset2, sum_model, 'r-', label='total')
+			#plt.ylim(ymin=0.1)
+
+		# Put labels
+		plt.xlabel("Offset (deg^2)")
+		plt.ylabel("Counts")
+		plt.legend(loc="upper right")
+
+		# Show counts spectra
+		plt.show()
+
+	except ImportError:
+		print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
+
+	# Return
+	return
+
+
+#==========================#
+# Main routine entry point #
+#==========================#
+if __name__ == '__main__':
+	"""
+	CTA unbinned analysis script.
+	
+	This script implements the analysis steps to perform an unbinned analysis.
+	The steps include:
+	- Event simulation
+	- Event selection
+	- Binned maximum likelihood fitting
+	
+	Two variations of the script are implemented. One that calls the classes
+	like executables, storing the intermediate results on disk, and a second
+	that uses the results of the precedent class directly without storing
+	any data on disk.
+	"""
+	# Initialise flags
+	need_help = False
+	show_data = True
+	
+	# Test for command line arguments
+	print sys.argv[0]
+	if (len(sys.argv) > 1):
+		if sys.argv[1] == "-h":
+			need_help = True
+		elif sys.argv[1] == "-noshow":
+			show_data = False
+		else:
+			need_help = True			
+
+	# Print help if needed and exit
+	if need_help:
+		print "Usage: example_binned.py [OPTIONS]"
+		print "     -h       Display this usage message"
+		print "     -noshow  Do not show data using matplotlib"
+		sys.exit()
+
+	# Dump header
+	print "*********************************"
+	print "* CTA unbinned analysis scripts *"
+	print "*********************************"
+	
+	# Remove any existing result files
+	list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
+	for files in list:
+		for file in files:
+			os.remove(file)
+	
+	# Save intermediate results on disk
+	print ""
+	print "Executable analysis pipeline:"
+	print "============================="
+	pipeline_v1()
+	
+	# Do analysis in memory
+	print ""
+	print "In memory analysis pipeline:"
+	print "============================"
+	pipeline_v2(show_data)
+	print ""