Francois Brun / ctools

Browse Code »

Commit d8a272885fd5a2018d87df351ab03a74c46e0822

Authored by knodlseder 2011-08-03 06:15:15 +0000

1 parent cb25c212

move here from test directory

Inline Side-by-side

Showing 2 changed files with 845 additions and 0 deletions

examples/binned_analysis.py 0 → 100755

View file @d8a2728

	1	+#! /usr/bin/env python
	2	+# ==========================================================================
	3	+# This script illustrates how to build a binned analysis pipeline using the
	4	+# ctatools. If matplotlib is installed, a counts spectrum and an offset
	5	+# histogram will be displayed.
	6	+#
	7	+# Copyright (C) 2011 Jurgen Knodlseder
	8	+#
	9	+# This program is free software: you can redistribute it and/or modify
	10	+# it under the terms of the GNU General Public License as published by
	11	+# the Free Software Foundation, either version 3 of the License, or
	12	+# (at your option) any later version.
	13	+#
	14	+# This program is distributed in the hope that it will be useful,
	15	+# but WITHOUT ANY WARRANTY; without even the implied warranty of
	16	+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	17	+# GNU General Public License for more details.
	18	+#
	19	+# You should have received a copy of the GNU General Public License
	20	+# along with this program. If not, see <http://www.gnu.org/licenses/>.
	21	+#
	22	+# ==========================================================================
	23	+from ctatools import *
	24	+from gammalib import *
	25	+from math import *
	26	+import os
	27	+import glob
	28	+
	29	+
	30	+# ============================= #
	31	+# Analysis pipeline - version 1 #
	32	+# ============================= #
	33	+def pipeline_v1():
	34	+ """
	35	+ Binned analysis pipeline - save intermediate results.
	36	+
	37	+ This function implements an analysis pipeline that successively calls
	38	+ ctobssim, ctbin and ctlike by saving the intermediate results as FITS
	39	+ files on disk. This replicates the way how the analysis would be done
	40	+ in a ftools-like approach.
	41	+ """
	42	+ # Set script parameters
	43	+ model_name = "${GAMMALIB}/share/models/crab.xml"
	44	+ events_name = "events.fits"
	45	+ cntmap_name = "cntmap.fits"
	46	+ result_name = "results.xml"
	47	+ caldb = "${GAMMALIB}/share/caldb/cta"
	48	+ irf = "kb_E_50h_v3"
	49	+ ra = 83.63
	50	+ dec = 22.01
	51	+ rad_sim = 10.0
	52	+ tstart = 0.0
	53	+ tstop = 1800.0
	54	+ emin = 0.1
	55	+ emax = 100.0
	56	+ enumbins = 20
	57	+ nxpix = 200
	58	+ nypix = 200
	59	+ binsz = 0.02
	60	+ coordsys = "CEL"
	61	+ proj = "CAR"
	62	+
	63	+ # Initialise timing
	64	+ wall_seconds = 0.0
	65	+ cpu_seconds = 0.0
	66	+
	67	+ # Simulate events
	68	+ sim = ctobssim()
	69	+ sim.logFileOpen() # We need this to explicitely open the log file in Python mode
	70	+ sim["infile"].filename(model_name)
	71	+ sim["outfile"].filename(events_name)
	72	+ sim["caldb"].string(caldb)
	73	+ sim["irf"].string(irf)
	74	+ sim["ra"].real(ra)
	75	+ sim["dec"].real(dec)
	76	+ sim["rad"].real(rad_sim)
	77	+ sim["tmin"].real(tstart)
	78	+ sim["tmax"].real(tstop)
	79	+ sim["emin"].real(emin)
	80	+ sim["emax"].real(emax)
	81	+ sim.execute()
	82	+ print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
	83	+
	84	+ # Update timing
	85	+ wall_seconds += sim.telapse()
	86	+ cpu_seconds += sim.celapse()
	87	+
	88	+ # Bin events into counts map
	89	+ bin = ctbin()
	90	+ bin.logFileOpen() # We need this to explicitely open the log file in Python mode
	91	+ bin["evfile"].filename(events_name)
	92	+ bin["outfile"].filename(cntmap_name)
	93	+ bin["emin"].real(emin)
	94	+ bin["emax"].real(emax)
	95	+ bin["enumbins"].integer(enumbins)
	96	+ bin["nxpix"].integer(nxpix)
	97	+ bin["nypix"].integer(nypix)
	98	+ bin["binsz"].real(binsz)
	99	+ bin["coordsys"].string(coordsys)
	100	+ bin["xref"].real(ra)
	101	+ bin["yref"].real(dec)
	102	+ bin["proj"].string(proj)
	103	+ bin.execute()
	104	+ print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
	105	+
	106	+ # Update timing
	107	+ wall_seconds += bin.telapse()
	108	+ cpu_seconds += bin.celapse()
	109	+
	110	+ # Perform maximum likelihood fitting
	111	+ like = ctlike()
	112	+ like.logFileOpen() # We need this to explicitely open the log file in Python mode
	113	+ like["cntmap"].filename(cntmap_name)
	114	+ like["srcmdl"].filename(model_name)
	115	+ like["outmdl"].filename(result_name)
	116	+ like["method"].string("BINNED")
	117	+ like["caldb"].string(caldb)
	118	+ like["irf"].string(irf)
	119	+ like.execute()
	120	+ print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
	121	+
	122	+ # Update timing
	123	+ wall_seconds += like.telapse()
	124	+ cpu_seconds += like.celapse()
	125	+
	126	+ # Show total times
	127	+ print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
	128	+ print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
	129	+
	130	+ # Return
	131	+ return
	132	+
	133	+
	134	+# ============================= #
	135	+# Analysis pipeline - version 2 #
	136	+# ============================= #
	137	+def pipeline_v2():
	138	+ """
	139	+ Binned analysis pipeline - keep intermediate results in memory.
	140	+
	141	+ This function implements an analysis pipeline that successively calls
	142	+ ctobssim, ctbin and ctlike without saving the intermediate results as
	143	+ FITS files on disk. All data are hold in memory.
	144	+
	145	+ At the end, results are plotted (if matplotlib is installed)
	146	+ """
	147	+ # Set script parameters
	148	+ model_name = "${GAMMALIB}/share/models/crab.xml"
	149	+ caldb = "${GAMMALIB}/share/caldb/cta"
	150	+ irf = "kb_E_50h_v3"
	151	+ ra = 83.63
	152	+ dec = 22.01
	153	+ rad_sim = 10.0
	154	+ tstart = 0.0
	155	+ tstop = 1800.0
	156	+ emin = 0.1
	157	+ emax = 100.0
	158	+ enumbins = 20
	159	+ nxpix = 200
	160	+ nypix = 200
	161	+ binsz = 0.02
	162	+ coordsys = "CEL"
	163	+ proj = "CAR"
	164	+
	165	+ # Initialise timing
	166	+ wall_seconds = 0.0
	167	+ cpu_seconds = 0.0
	168	+
	169	+ # Simulate events
	170	+ sim = ctobssim()
	171	+ sim["infile"].filename(model_name)
	172	+ sim["caldb"].string(caldb)
	173	+ sim["irf"].string(irf)
	174	+ sim["ra"].real(ra)
	175	+ sim["dec"].real(dec)
	176	+ sim["rad"].real(rad_sim)
	177	+ sim["tmin"].real(tstart)
	178	+ sim["tmax"].real(tstop)
	179	+ sim["emin"].real(emin)
	180	+ sim["emax"].real(emax)
	181	+ sim.run()
	182	+ print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
	183	+
	184	+ # Update timing
	185	+ wall_seconds += sim.telapse()
	186	+ cpu_seconds += sim.celapse()
	187	+
	188	+ # Bin events into counts map
	189	+ bin = ctbin(sim.obs())
	190	+ bin["emin"].real(emin)
	191	+ bin["emax"].real(emax)
	192	+ bin["enumbins"].integer(enumbins)
	193	+ bin["nxpix"].integer(nxpix)
	194	+ bin["nypix"].integer(nypix)
	195	+ bin["binsz"].real(binsz)
	196	+ bin["coordsys"].string(coordsys)
	197	+ bin["xref"].real(ra)
	198	+ bin["yref"].real(dec)
	199	+ bin["proj"].string(proj)
	200	+ bin.run()
	201	+ print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
	202	+
	203	+ # Update timing
	204	+ wall_seconds += bin.telapse()
	205	+ cpu_seconds += bin.celapse()
	206	+
	207	+ # Perform maximum likelihood fitting
	208	+ like = ctlike(bin.obs())
	209	+ like.run()
	210	+ print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
	211	+
	212	+ # Update timing
	213	+ wall_seconds += like.telapse()
	214	+ cpu_seconds += like.celapse()
	215	+
	216	+ # Show total times
	217	+ print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
	218	+ print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
	219	+
	220	+ # Show model fitting results
	221	+ #print like.obs().models()
	222	+
	223	+ # Plot counts
	224	+ plot_counts(bin.obs())
	225	+
	226	+ # Return
	227	+ return
	228	+
	229	+
	230	+# =========== #
	231	+# Plot counts #
	232	+# =========== #
	233	+def plot_counts(observations):
	234	+ """
	235	+ Plot counts.
	236	+ """
	237	+ # Only proceed if matplotlib is available
	238	+ try:
	239	+ # Import matplotlib
	240	+ import matplotlib.pyplot as plt
	241	+
	242	+ # Set legend fontsize
	243	+ params = {'legend.fontsize': 10}
	244	+ plt.rcParams.update(params)
	245	+
	246	+ # Set plot styles
	247	+ styles = ['b-', 'g-', 'y-', 'n-']
	248	+
	249	+ # Dump header
	250	+ print ""
	251	+ print "Make plots (using matplotlib):"
	252	+ print "=============================="
	253	+
	254	+ # Create figure 1
	255	+ plt.figure(1,figsize=(12,6))
	256	+ plt.subplots_adjust(hspace=.7)
	257	+
	258	+ # Create subplot 1
	259	+ plt.subplot(121)
	260	+ plt.title("Spectrum (summed over all pixels)")
	261	+
	262	+ # Loop over observations
	263	+ for obs in observations:
	264	+
	265	+ # Get event cube
	266	+ cube = cast_GCTAEventCube(obs.events())
	267	+
	268	+ # Create energy axis
	269	+ energy = []
	270	+ ebounds = cube.ebounds()
	271	+ for i in range(ebounds.size()):
	272	+ energy.append(ebounds.elogmean(i).TeV())
	273	+
	274	+ # Create spectrum
	275	+ print "Extract data:"
	276	+ counts = [0.0 for i in range(ebounds.size())]
	277	+ for bin in cube:
	278	+ index = ebounds.index(bin.energy())
	279	+ counts[index] = counts[index] + bin.counts()
	280	+
	281	+ # Create error bars
	282	+ error = [sqrt(c) for c in counts]
	283	+
	284	+ # Plot spectrum
	285	+ plt.loglog(energy, counts, 'ro', label='data')
	286	+ plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
	287	+
	288	+ # Extract models
	289	+ print "Extract models:"
	290	+ sum_model = [0.0 for i in range(ebounds.size())]
	291	+ for k, m in enumerate(observations.models()):
	292	+ print "- "+m.name()
	293	+ model = [0.0 for i in range(ebounds.size())]
	294	+ for bin in cube:
	295	+ index = ebounds.index(bin.energy())
	296	+ prob = m.eval(bin, obs)
	297	+ model[index] = model[index] + prob * bin.size()
	298	+ for i in range(ebounds.size()):
	299	+ sum_model[i] = sum_model[i] + model[i]
	300	+ plt.loglog(energy, model, styles[k], label=m.name())
	301	+ plt.loglog(energy, sum_model, 'r-', label='total')
	302	+
	303	+ # Put labels
	304	+ plt.xlabel("Energy (TeV)")
	305	+ plt.ylabel("Counts")
	306	+ plt.legend(loc="lower left")
	307	+
	308	+ # Create subplot 2
	309	+ plt.subplot(122)
	310	+ plt.title("Offset (summed over all energies)")
	311	+
	312	+ # Set Crab direction
	313	+ crab = GSkyDir()
	314	+ crab.radec_deg(83.63, 22.01)
	315	+
	316	+ # Loop over observations
	317	+ for obs in observations:
	318	+
	319	+ # Get event cube
	320	+ cube = cast_GCTAEventCube(obs.events())
	321	+
	322	+ # Create offset histogram
	323	+ print "Extract data:"
	324	+ nx = 30
	325	+ doffset2 = 0.01
	326	+ offset2 = [(i+0.5)*doffset2 for i in range(nx)]
	327	+ counts = [0.0 for i in range(nx)]
	328	+ for bin in cube:
	329	+ off = bin.dir().dist_deg(crab)
	330	+ off2 = off*off
	331	+ index = int(off2/doffset2)
	332	+ if index < nx:
	333	+ counts[index] = counts[index] + bin.counts()
	334	+
	335	+ # Create error bars
	336	+ error = [sqrt(c) for c in counts]
	337	+
	338	+ # Plot distribution
	339	+ plt.semilogy(offset2, counts, 'ro', label='data')
	340	+
	341	+ # Extract models
	342	+ print "Extract models:"
	343	+ sum_model = [0.0 for i in range(nx)]
	344	+ for k, m in enumerate(observations.models()):
	345	+ print "- "+m.name()
	346	+ model = [0.0 for i in range(nx)]
	347	+ for bin in cube:
	348	+ off = bin.dir().dist_deg(crab)
	349	+ off2 = off*off
	350	+ index = int(off2/doffset2)
	351	+ if index < nx:
	352	+ prob = m.eval(bin, obs)
	353	+ model[index] = model[index] + prob * bin.size()
	354	+ for i in range(nx):
	355	+ sum_model[i] = sum_model[i] + model[i]
	356	+ plt.plot(offset2, model, styles[k], label=m.name())
	357	+ plt.plot(offset2, sum_model, 'r-', label='total')
	358	+ plt.ylim(ymin=0.1)
	359	+
	360	+ # Put labels
	361	+ plt.xlabel("Offset (deg^2)")
	362	+ plt.ylabel("Counts")
	363	+ plt.legend(loc="upper right")
	364	+
	365	+ # Show counts spectra
	366	+ plt.show()
	367	+
	368	+ except:
	369	+ print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
	370	+
	371	+ # Return
	372	+ return
	373	+
	374	+
	375	+#==========================#
	376	+# Main routine entry point #
	377	+#==========================#
	378	+if __name__ == '__main__':
	379	+ """
	380	+ CTA binned analysis script.
	381	+
	382	+ This script implements the analysis steps to perform a binned analysis.
	383	+ The steps include:
	384	+ - Event simulation
	385	+ - Event binning
	386	+ - Binned maximum likelihood fitting
	387	+
	388	+ Two variations of the script are implemented. One that calls the classes
	389	+ like executables, storing the intermediate results on disk, and a second
	390	+ that uses the results of the precedent class directly without storing
	391	+ any data on disk.
	392	+
	393	+ The second variation will also show a counts spectrum of the data after
	394	+ fitting (only if matplotlib is available on your system). Note that the
	395	+ counts spectrum is summed over the entire data space, i.e. no event selection
	396	+ around the Crab has been performed.
	397	+ """
	398	+ # Dump header
	399	+ print "*************************************"
	400	+ print "* CTA binned analysis scripts *"
	401	+ print "*************************************"
	402	+ print "... this script will take 1-2 minutes"
	403	+
	404	+ # Remove any existing result files
	405	+ list = [glob.glob(".fits"), glob.glob(".log"), glob.glob("*.xml")]
	406	+ for files in list:
	407	+ for file in files:
	408	+ os.remove(file)
	409	+
	410	+ # Save intermediate results on disk
	411	+ print ""
	412	+ print "Executable analysis pipeline:"
	413	+ print "============================="
	414	+ pipeline_v1()
	415	+
	416	+ # Do analysis in memory
	417	+ print ""
	418	+ print "In memory analysis pipeline:"
	419	+ print "============================"
	420	+ pipeline_v2()
	421	+ print ""
...	...

examples/unbinned_analysis.py 0 → 100755

View file @d8a2728

	1	+#! /usr/bin/env python
	2	+# ==========================================================================
	3	+# This script illustrates how to build an unbinned analysis pipeline with
	4	+# the ctatools. If matplotlib is installed, the spectrum will be displayed
	5	+# on the screen.
	6	+#
	7	+# Copyright (C) 2011 Jurgen Knodlseder
	8	+#
	9	+# This program is free software: you can redistribute it and/or modify
	10	+# it under the terms of the GNU General Public License as published by
	11	+# the Free Software Foundation, either version 3 of the License, or
	12	+# (at your option) any later version.
	13	+#
	14	+# This program is distributed in the hope that it will be useful,
	15	+# but WITHOUT ANY WARRANTY; without even the implied warranty of
	16	+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	17	+# GNU General Public License for more details.
	18	+#
	19	+# You should have received a copy of the GNU General Public License
	20	+# along with this program. If not, see <http://www.gnu.org/licenses/>.
	21	+#
	22	+# ==========================================================================
	23	+from ctatools import *
	24	+from gammalib import *
	25	+from math import *
	26	+import os
	27	+import glob
	28	+import sys
	29	+
	30	+
	31	+# ============================= #
	32	+# Analysis pipeline - version 1 #
	33	+# ============================= #
	34	+def pipeline_v1():
	35	+ """
	36	+ Unbinned analysis pipeline - save intermediate results.
	37	+
	38	+ This function implements an analysis pipeline that successively calls
	39	+ ctobssim, ctselect and ctlike by saving the intermediate results as FITS
	40	+ files on disk. This replicates the way how the analysis would be done
	41	+ in a ftools-like approach.
	42	+ """
	43	+ # Set script parameters
	44	+ model_name = "${GAMMALIB}/share/models/crab.xml"
	45	+ events_name = "events.fits"
	46	+ selected_events_name = "selected_events.fits"
	47	+ result_name = "results.xml"
	48	+ caldb = "${GAMMALIB}/share/caldb/cta"
	49	+ irf = "kb_E_50h_v3"
	50	+ ra = 83.63
	51	+ dec = 22.01
	52	+ rad_sim = 10.0
	53	+ tstart = 0.0
	54	+ tstop = 1800.0
	55	+ emin = 0.1
	56	+ emax = 100.0
	57	+ rad_select = 3.0
	58	+
	59	+ # Initialise timing
	60	+ wall_seconds = 0.0
	61	+ cpu_seconds = 0.0
	62	+
	63	+ # Simulate events
	64	+ sim = ctobssim()
	65	+ sim.logFileOpen()
	66	+ sim["infile"].filename(model_name)
	67	+ sim["outfile"].filename(events_name)
	68	+ sim["caldb"].string(caldb)
	69	+ sim["irf"].string(irf)
	70	+ sim["ra"].real(ra)
	71	+ sim["dec"].real(dec)
	72	+ sim["rad"].real(rad_sim)
	73	+ sim["tmin"].real(tstart)
	74	+ sim["tmax"].real(tstop)
	75	+ sim["emin"].real(emin)
	76	+ sim["emax"].real(emax)
	77	+ sim.execute()
	78	+ print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
	79	+
	80	+ # Update timing
	81	+ wall_seconds += sim.telapse()
	82	+ cpu_seconds += sim.celapse()
	83	+
	84	+ # Select events
	85	+ select = ctselect()
	86	+ select.logFileOpen()
	87	+ select["infile"].filename(events_name)
	88	+ select["outfile"].filename(selected_events_name)
	89	+ select["ra"].real(ra)
	90	+ select["dec"].real(dec)
	91	+ select["rad"].real(rad_select)
	92	+ select["tmin"].real(tstart)
	93	+ select["tmax"].real(tstop)
	94	+ select["emin"].real(emin)
	95	+ select["emax"].real(emax)
	96	+ select.execute()
	97	+ print "Select event ("+str(select.celapse())+" CPU seconds)"
	98	+
	99	+ # Update timing
	100	+ wall_seconds += select.telapse()
	101	+ cpu_seconds += select.celapse()
	102	+
	103	+ # Perform maximum likelihood fitting
	104	+ like = ctlike()
	105	+ like.logFileOpen()
	106	+ like["evfile"].filename(selected_events_name)
	107	+ like["srcmdl"].filename(model_name)
	108	+ like["outmdl"].filename(result_name)
	109	+ like["method"].string("UNBINNED")
	110	+ like["caldb"].string(caldb)
	111	+ like["irf"].string(irf)
	112	+ like.execute()
	113	+ print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
	114	+
	115	+ # Update timing
	116	+ wall_seconds += like.telapse()
	117	+ cpu_seconds += like.celapse()
	118	+
	119	+ # Show total times
	120	+ print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
	121	+ print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
	122	+
	123	+ # Return
	124	+ return
	125	+
	126	+
	127	+# ============================= #
	128	+# Analysis pipeline - version 2 #
	129	+# ============================= #
	130	+def pipeline_v2(show_data):
	131	+ """
	132	+ Unbinned analysis pipeline - keep intermediate results in memory.
	133	+
	134	+ This function implements an analysis pipeline that successively calls
	135	+ ctobssim, ctselect and ctlike without saving the intermediate results as
	136	+ FITS files on disk. All data is only hold in memory.
	137	+ """
	138	+ # Set script parameters
	139	+ model_name = "${GAMMALIB}/share/models/crab.xml"
	140	+ caldb = "${GAMMALIB}/share/caldb/cta"
	141	+ irf = "kb_E_50h_v3"
	142	+ ra = 83.63
	143	+ dec = 22.01
	144	+ rad_sim = 10.0
	145	+ tstart = 0.0
	146	+ tstop = 1800.0
	147	+ emin = 0.1
	148	+ emax = 100.0
	149	+ rad_select = 3.0
	150	+
	151	+ # Initialise timing
	152	+ wall_seconds = 0.0
	153	+ cpu_seconds = 0.0
	154	+
	155	+ # Simulate events
	156	+ sim = ctobssim()
	157	+ sim["infile"].filename(model_name)
	158	+ sim["caldb"].string(caldb)
	159	+ sim["irf"].string(irf)
	160	+ sim["ra"].real(ra)
	161	+ sim["dec"].real(dec)
	162	+ sim["rad"].real(rad_sim)
	163	+ sim["tmin"].real(tstart)
	164	+ sim["tmax"].real(tstop)
	165	+ sim["emin"].real(emin)
	166	+ sim["emax"].real(emax)
	167	+ sim.run()
	168	+ print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
	169	+
	170	+ # Update timing
	171	+ wall_seconds += sim.telapse()
	172	+ cpu_seconds += sim.celapse()
	173	+
	174	+ # Select events
	175	+ select = ctselect(sim.obs())
	176	+ select["ra"].real(ra)
	177	+ select["dec"].real(dec)
	178	+ select["rad"].real(rad_select)
	179	+ select["tmin"].real(tstart)
	180	+ select["tmax"].real(tstop)
	181	+ select["emin"].real(emin)
	182	+ select["emax"].real(emax)
	183	+ select.run()
	184	+ print "Select event ("+str(select.celapse())+" CPU seconds)"
	185	+
	186	+ # Update timing
	187	+ wall_seconds += select.telapse()
	188	+ cpu_seconds += select.celapse()
	189	+
	190	+ # Perform maximum likelihood fitting
	191	+ like = ctlike(select.obs())
	192	+ like.run()
	193	+ print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
	194	+
	195	+ # Update timing
	196	+ wall_seconds += like.telapse()
	197	+ cpu_seconds += like.celapse()
	198	+
	199	+ # Show total times
	200	+ print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
	201	+ print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
	202	+
	203	+ # Show model fitting results
	204	+ #print like.obs().models()
	205	+
	206	+ # Optionally plot counts
	207	+ if show_data:
	208	+ plot_counts(like.obs())
	209	+
	210	+ # Return
	211	+ return
	212	+
	213	+
	214	+# =========== #
	215	+# Plot counts #
	216	+# =========== #
	217	+def plot_counts(observations):
	218	+ """
	219	+ Plot counts.
	220	+ """
	221	+ # Only proceed if matplotlib is available
	222	+ try:
	223	+ # Import matplotlib
	224	+ import matplotlib.pyplot as plt
	225	+
	226	+ # Set legend fontsize
	227	+ params = {'legend.fontsize': 10}
	228	+ plt.rcParams.update(params)
	229	+
	230	+ # Set plot styles
	231	+ styles = ['b-', 'g-', 'y-', 'n-']
	232	+
	233	+ # Dump header
	234	+ print ""
	235	+ print "Make plots (using matplotlib):"
	236	+ print "=============================="
	237	+
	238	+ # Create figure 1
	239	+ plt.figure(1,figsize=(12,6))
	240	+ plt.subplots_adjust(hspace=.7)
	241	+
	242	+ # Create subplot 1
	243	+ plt.subplot(121)
	244	+ plt.title("Spectrum (summed over all pixels)")
	245	+
	246	+ # Loop over observations
	247	+ for obs in observations:
	248	+
	249	+ # Get event list
	250	+ list = cast_GCTAEventList(obs.events())
	251	+
	252	+ # Create energy axis
	253	+ ebounds = GEbounds()
	254	+ emin = GEnergy()
	255	+ emax = GEnergy()
	256	+ emin.TeV(0.1)
	257	+ emax.TeV(100.0)
	258	+ ebounds.setlog(emin, emax, 10)
	259	+ energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
	260	+
	261	+ # Create spectrum
	262	+ print "Extract data:"
	263	+ counts = [0.0 for i in range(ebounds.size())]
	264	+ print list.size()
	265	+ for atom in list:
	266	+ index = ebounds.index(atom.energy())
	267	+ counts[index] = counts[index] + 1.0
	268	+
	269	+ # Create error bars
	270	+ error = [sqrt(c) for c in counts]
	271	+
	272	+ # Plot spectrum
	273	+ plt.loglog(energy, counts, 'ro', label='data')
	274	+ #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
	275	+
	276	+ # Extract models
	277	+ print "Extract models:"
	278	+ sum_model = [0.0 for i in range(ebounds.size())]
	279	+ for k, m in enumerate(observations.models()):
	280	+ print "- "+m.name()
	281	+ model = [0.0 for i in range(ebounds.size())]
	282	+ for atom in list:
	283	+ index = ebounds.index(atom.energy())
	284	+ prob = m.eval(atom, obs)
	285	+ model[index] = model[index] + prob * atom.size()
	286	+ for i in range(ebounds.size()):
	287	+ sum_model[i] = sum_model[i] + model[i]
	288	+ #plt.loglog(energy, model, styles[k], label=m.name())
	289	+ #plt.loglog(energy, sum_model, 'r-', label='total')
	290	+
	291	+ # Put labels
	292	+ plt.xlabel("Energy (TeV)")
	293	+ plt.ylabel("Counts")
	294	+ plt.legend(loc="lower left")
	295	+
	296	+ # Create subplot 2
	297	+ plt.subplot(122)
	298	+ plt.title("Offset (summed over all energies)")
	299	+
	300	+ # Set Crab direction
	301	+ crab = GSkyDir()
	302	+ crab.radec_deg(83.63, 22.01)
	303	+
	304	+ # Loop over observations
	305	+ for obs in observations:
	306	+
	307	+ # Get event cube
	308	+ list = cast_GCTAEventList(obs.events())
	309	+
	310	+ # Create offset histogram
	311	+ print "Extract data:"
	312	+ nx = 30
	313	+ doffset2 = 0.01
	314	+ offset2 = [(i+0.5)*doffset2 for i in range(nx)]
	315	+ counts = [0.0 for i in range(nx)]
	316	+ for atom in list:
	317	+ off = atom.dir().dist_deg(crab)
	318	+ off2 = off*off
	319	+ index = int(off2/doffset2)
	320	+ if index < nx:
	321	+ counts[index] = counts[index] + 1.0
	322	+
	323	+ # Create error bars
	324	+ error = [sqrt(c) for c in counts]
	325	+
	326	+ # Plot distribution
	327	+ plt.semilogy(offset2, counts, 'ro', label='data')
	328	+
	329	+ # Extract models
	330	+ print "Extract models:"
	331	+ sum_model = [0.0 for i in range(nx)]
	332	+ for k, m in enumerate(observations.models()):
	333	+ print "- "+m.name()
	334	+ model = [0.0 for i in range(nx)]
	335	+ for atom in list:
	336	+ off = atom.dir().dist_deg(crab)
	337	+ off2 = off*off
	338	+ index = int(off2/doffset2)
	339	+ if index < nx:
	340	+ prob = m.eval(atom, obs)
	341	+ model[index] = model[index] + prob * atom.size()
	342	+ for i in range(nx):
	343	+ sum_model[i] = sum_model[i] + model[i]
	344	+ #plt.plot(offset2, model, styles[k], label=m.name())
	345	+ #plt.plot(offset2, sum_model, 'r-', label='total')
	346	+ #plt.ylim(ymin=0.1)
	347	+
	348	+ # Put labels
	349	+ plt.xlabel("Offset (deg^2)")
	350	+ plt.ylabel("Counts")
	351	+ plt.legend(loc="upper right")
	352	+
	353	+ # Show counts spectra
	354	+ plt.show()
	355	+
	356	+ except ImportError:
	357	+ print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
	358	+
	359	+ # Return
	360	+ return
	361	+
	362	+
	363	+#==========================#
	364	+# Main routine entry point #
	365	+#==========================#
	366	+if __name__ == '__main__':
	367	+ """
	368	+ CTA unbinned analysis script.
	369	+
	370	+ This script implements the analysis steps to perform an unbinned analysis.
	371	+ The steps include:
	372	+ - Event simulation
	373	+ - Event selection
	374	+ - Binned maximum likelihood fitting
	375	+
	376	+ Two variations of the script are implemented. One that calls the classes
	377	+ like executables, storing the intermediate results on disk, and a second
	378	+ that uses the results of the precedent class directly without storing
	379	+ any data on disk.
	380	+ """
	381	+ # Initialise flags
	382	+ need_help = False
	383	+ show_data = True
	384	+
	385	+ # Test for command line arguments
	386	+ print sys.argv[0]
	387	+ if (len(sys.argv) > 1):
	388	+ if sys.argv[1] == "-h":
	389	+ need_help = True
	390	+ elif sys.argv[1] == "-noshow":
	391	+ show_data = False
	392	+ else:
	393	+ need_help = True
	394	+
	395	+ # Print help if needed and exit
	396	+ if need_help:
	397	+ print "Usage: example_binned.py [OPTIONS]"
	398	+ print " -h Display this usage message"
	399	+ print " -noshow Do not show data using matplotlib"
	400	+ sys.exit()
	401	+
	402	+ # Dump header
	403	+ print "*********************************"
	404	+ print "* CTA unbinned analysis scripts *"
	405	+ print "*********************************"
	406	+
	407	+ # Remove any existing result files
	408	+ list = [glob.glob(".fits"), glob.glob(".log"), glob.glob("*.xml")]
	409	+ for files in list:
	410	+ for file in files:
	411	+ os.remove(file)
	412	+
	413	+ # Save intermediate results on disk
	414	+ print ""
	415	+ print "Executable analysis pipeline:"
	416	+ print "============================="
	417	+ pipeline_v1()
	418	+
	419	+ # Do analysis in memory
	420	+ print ""
	421	+ print "In memory analysis pipeline:"
	422	+ print "============================"
	423	+ pipeline_v2(show_data)
	424	+ print ""
...	...