From a0bc9bdf4397173fbe301dacec03224575318355 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 24 Nov 2014 09:05:15 +0100 Subject: [PATCH] Correct analysis scripts --- ChangeLog | 8 ++++++++ NEWS | 184 +++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------- examples/make_binned_analysis.py | 730 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- examples/make_unbinned_analysis.py | 723 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 4 files changed, 755 insertions(+), 890 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1ba9210..f49ef45 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2014-11-24 Juergen Knoedlseder + + * Version 1.0.0 released + ======================== + + Review and adapt example scripts + + 2014-11-08 Juergen Knoedlseder * Version 0.8.0 released diff --git a/NEWS b/NEWS index 75df195..e012065 100644 --- a/NEWS +++ b/NEWS @@ -1,8 +1,8 @@ -New Features and Important Changes in ctools 00-08-00 +New Features and Important Changes in ctools 1.0.0 Author(s) : Juergen Knoedlseder -8 November 2014 +24 November 2014 1. Introduction @@ -10,16 +10,6 @@ Author(s) : Juergen Knoedlseder This document briefly summarizes executable by executable the most important changes made since the last release of ctools. -For all ctools, the public method interface was cleaned up so that -only methods that are really useful for an outside client are -exposed. Low-level methods that are used by the application are now -protected, as improper usage of these methods may lead to -unpredictable results. - -The ctool base class has been introduced and all tools now derive -from ctool instead of GApplication. This allows to avoid code -duplications and makes the code generally slimmer and more readable. - 2. Configuration ----------------- @@ -33,108 +23,32 @@ None 4. Support model ---------------- -A support module has been introduced that will contain common classes -and functions. For the moment the only class in this module is the -ctool base class. +None 5. ctobssim - CTA observation simulation ---------------------------------------- -ctobssim was not fully thread save as the models were shared among the -threads which led to concurrent memory access when evaluating the Monte -Carlo cache. The models are now copied within the parallel section, -avoiding concurrent memory access (#1258). - -ctobssim now also supports the handling of energy dispersion. Energy -dispersion can be switched on or off using the hidden edisp parameter. - -The memory needs for ctobssim have been reduced by first working on a -local copy of the CTA observation (this de-allocates properly and memory -cache associated to the response functions), and by saving and disposing -the event lists immediately if called using the execute() method (this -is the method called when the tool is run from the command line). +None 6. ctselect - CTA event selection --------------------------------- -The temporary file is now properly close to avoid file opening errors -when iterating over ctselect (#1149). - -Non-CTA time references are now correctly handled by ctselect. While the -tmin and tmax parameters are given in CTA reference time, they are -internally converted to the reference time of the GTI before event -selection (#1329). - -ctselect now supports save and user energy thresholds (#1039). There -are three different possibilities for specifying a threshold. They -are toggled using the hidden usethres parameter: - -NONE: Use the emin and emax parameters to select an energy range. - This is the default mode which is identical to the previous - workings. The only difference is now that if the emin and/or - emax parameter is set to 0.0, the threshold will be ignored. - This allows to select all events above or below a given - threshold. - -DEFAULT: Use the save threshold that is found in the effective - area FITS header keywords LO_THRES and HI_THRES. If the - emin and emax parameters are non-zero, they will provide - an additional energy constraint. - -USER: Use the user specified threshold that are found in the - XML observation definition file. If the emin and emax - parameters are non-zero, they will provide an additional - energy constraint. - -Note that any existing energy boundaries will be used as envelope, -i.e. a wider energy threshold specification will have no effect. +None 7. ctbin - CTA event binning ---------------------------- -WARNING: The philosophy behind the ctbin tool changed fundamentally! - -ctbin now always produces a single counts cube on output. If a single -event list is specified on input, the functionnality is basically -unchanged. If multiple event lists are specified, however, events are -now merged into a single counts cube file. No observation definition -XML file is written anymore! - -So far, ctbin still writes pointing information in the counts cube, -although this is not necessarily meaningful if merging event lists -with different pointings. In the long run, pointing information will -be summarized in a pointings definition file, and will be removed -from the counts cube header. The logic is that the pointing definition -is set to the map centre. This is only identical to the pointing in -case of binning a single event list and using the usepnt=yes option -(or specifying the pointing manually as the xref and yref parameters). - -Memory usage has been reduced by disposing the event lists after their -usage. +None 8. ctlike - Maximum likelihood fitting -------------------------------------- -ctlike now has full access to the calibration database, allowing -the use for 2D response functions in FITS files. - -ctlike now also supports the handling of energy dispersion. Energy -dispersion can be switched on or off using the hidden edisp parameter. +None 9. cttsmap - Generation of Test Statistics map ---------------------------------------------- -The cttsmap tool has been added that provides a test statistics map -for a given set of observations. This tool is not restricted to CTA -data and will work for all types of instruments. - -cttsmap will move a test sources over a specified grid of positions -and compute the likelihood difference between having the source at -a given position and having not the source. Twice the likelihood -difference is called the Test Statistic. The computation is done for -a grid of positions, and the results is saved into a FITS file as -TS map. Also the value of the fitted parameters are stored as -extension in the FITS file. +None 10. ctskymap - CTA sky mapping @@ -144,105 +58,45 @@ None 11. ctmodel - CTA counts cube modelling --------------------------------------- -WARNING: The philosophy behind the ctmodel tool changed fundamentally! - -ctmodel now always produces a single model cube on output. There are -two important input parameters: "infile" and "obsfile". "infile" -specifies a counts or model cube that is used to define the binning -of the model cube. The functionnality for this parameter is as before. -If "infile=NONE", the cube will be defined according to task parameters -that specify the spatial and the energy binning. The "obsfile" -parameter is new and allows specifying the observation definition -for which the model cube should be computed. If an event list or a -counts cube is presented, the pointing and Good Time Intervals will -be extracted from the FITS file header. If an observation definition -XML file is presented, ctmodel will now loop over all CTA observations -in the file and add the model values to the model cube. If no -response information is associated to the "obsfile" (which is typically -the case for event lists or counts cubes), the "caldb" and "irf" -parameters will be used to identify the response. - -ctmodel now has full access to the calibration database, allowing -the use for 2D response functions in FITS files. - -ctmodel now also supports the handling of energy dispersion. Energy -dispersion can be switched on or off using the hidden edisp parameter. - -Memory usage has been reduced by disposing the event lists after their -usage. +None 12. ctexpcube - CTA exposure cube generation -------------------------------------------- -The first version of the ctexpcube tool is included in the release. -ctexpcube computes the exposure cube for an observation definition -file. The cube definition may be taken from a counts cube or may be -specified by user parameters. +None 13. ctpsfcube - CTA point spread function cube generation --------------------------------------------------------- -The first version of the ctpsfcube tool is included in the release. -ctpsfcube computes the point spread function cube for an observation -definition file. The cube definition may be taken from a counts cube -or may be specified by user parameters. +None 14. ctbkgcube - CTA background cube generation ---------------------------------------------- -The first version of the ctbkgcube tool is included in the release. -ctbkgcube computes the background cube for an observation definition -file. The cube definition may be taken from a counts cube or may be -specified by user parameters. The tools also updates the input model -by replacing all IACT background model components by a single -background cube background model with a power law as spectral -parameter. The model container is available to Python using the -models() method or can be written into an XML file on disk. +None 15. ctcubemask - CTA cube filtering ----------------------------------- -The first version of the ctcubemask tool is included in the release. -ctcubemask allows filtering of bins in a counts cube by setting the -bin content to -1. Filtering can be done by energy, for a circular -Region of Interest (ROI) or using ds9 exclusion regions. +None 16. ctbutterfly - Generate butterfly diagrams ============================================= -The first version of the ctbutterfly tool is included in the release. -This tools derives a butterfly diagram from the fitted model parameters -for one specified source. The result will be written into an ASCII -file. +None 17. Examples ------------ -Added the following example scripts: -* pipeline_stacked.py (simulation and stacked analysis) +The following scripts have been adapted to fit the actual +interface definition +* make_binned_analysis.py +* make_unbinned_analysis.py 18. Scripts ----------- -Add logL and event statistics to cspull output file. - -Convert all tabs into spaces. - -The obsutils script has now been moved into the ctools module namespace. - -The following functions have been added to obsutils: -* set_obs() replaces the now obsolete set() function -* set_obs_list() returns an observation container filled with a list of - observations -* set_obs_patterns() defines standard patterns of CTA observations. - Current options are: - - "single": single pointing - - "four": four pointings offset by specific angle around centre -* spectrum() generates spectral points by fitting a source model in - specified energy bins. - -Add the csresmap script that uses ctbin and ctmodel to generate a -residual counts cube for a given model. +None 19. caldb diff --git a/examples/make_binned_analysis.py b/examples/make_binned_analysis.py index 86cacad..6cc0b3c 100755 --- a/examples/make_binned_analysis.py +++ b/examples/make_binned_analysis.py @@ -4,7 +4,7 @@ # ctools. If matplotlib is installed, a counts spectrum and an offset # histogram will be displayed. # -# Copyright (C) 2011-2012 Jurgen Knodlseder +# Copyright (C) 2011-2014 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 @@ -20,404 +20,404 @@ # along with this program. If not, see . # # ========================================================================== -from ctools import * -from gammalib import * -from math import * +import gammalib +import ctools import os import glob import sys +import math try: - import matplotlib.pyplot as plt - has_matplotlib = True + import matplotlib.pyplot as plt + has_matplotlib = True except: - has_matplotlib = False + has_matplotlib = False # ============================= # # Analysis pipeline - version 1 # # ============================= # def pipeline_v1(): - """ - Binned analysis pipeline - save intermediate results. + """ + Binned analysis pipeline - save intermediate results. This function implements an analysis pipeline that successively calls ctobssim, ctbin 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 = "${CTOOLS}/share/models/crab.xml" - events_name = "events.fits" - cntmap_name = "cntmap.fits" - result_name = "results.xml" - caldb = "${CALDB}/data/cta/dummy/bcf" - irf = "cta_dummy_irf" - ra = 83.63 - dec = 22.01 - rad_sim = 10.0 - tstart = 0.0 - tstop = 1800.0 - emin = 0.1 - emax = 100.0 - enumbins = 20 - nxpix = 200 - nypix = 200 - binsz = 0.02 - coordsys = "CEL" - proj = "CAR" - - # Initialise timing - wall_seconds = 0.0 - cpu_seconds = 0.0 - - # Simulate events - sim = ctobssim() - sim.logFileOpen() # We need this to explicitely open the log file in Python mode - 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() - sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n") - - # Update timing - wall_seconds += sim.telapse() - cpu_seconds += sim.celapse() - - # Bin events into counts map - bin = ctbin() - bin.logFileOpen() # We need this to explicitely open the log file in Python mode - bin["evfile"].filename(events_name) - bin["outfile"].filename(cntmap_name) - bin["ebinalg"].string("LOG") - bin["emin"].real(emin) - bin["emax"].real(emax) - bin["enumbins"].integer(enumbins) - bin["nxpix"].integer(nxpix) - bin["nypix"].integer(nypix) - bin["binsz"].real(binsz) - bin["coordsys"].string(coordsys) - bin["xref"].real(ra) - bin["yref"].real(dec) - bin["proj"].string(proj) - bin.execute() - sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n") - - # Update timing - wall_seconds += bin.telapse() - cpu_seconds += bin.celapse() - - # Perform maximum likelihood fitting - like = ctlike() - like.logFileOpen() # We need this to explicitely open the log file in Python mode - like["infile"].filename(cntmap_name) - like["srcmdl"].filename(model_name) - like["outmdl"].filename(result_name) - like["caldb"].string(caldb) - like["irf"].string(irf) - like.execute() - sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n") - - # Update timing - wall_seconds += like.telapse() - cpu_seconds += like.celapse() - - # Show total times - sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n") - sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n") - - # Return - return + """ + # Set script parameters + model_name = "${CTOOLS}/share/models/crab.xml" + events_name = "events.fits" + cntmap_name = "cntmap.fits" + result_name = "results.xml" + caldb = "${CALDB}/data/cta/dummy/bcf" + irf = "cta_dummy_irf" + ra = 83.63 + dec = 22.01 + rad_sim = 10.0 + tstart = 0.0 + tstop = 1800.0 + emin = 0.1 + emax = 100.0 + enumbins = 20 + nxpix = 200 + nypix = 200 + binsz = 0.02 + coordsys = "CEL" + proj = "CAR" + + # Initialise timing + wall_seconds = 0.0 + cpu_seconds = 0.0 + + # Simulate events + sim = ctools.ctobssim() + sim.logFileOpen() # We need this to explicitely open the log file in Python mode + 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() + sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n") + + # Update timing + wall_seconds += sim.telapse() + cpu_seconds += sim.celapse() + + # Bin events into counts map + bin = ctools.ctbin() + bin.logFileOpen() # We need this to explicitely open the log file in Python mode + bin["evfile"].filename(events_name) + bin["outfile"].filename(cntmap_name) + bin["ebinalg"].string("LOG") + bin["emin"].real(emin) + bin["emax"].real(emax) + bin["enumbins"].integer(enumbins) + bin["nxpix"].integer(nxpix) + bin["nypix"].integer(nypix) + bin["binsz"].real(binsz) + bin["coordsys"].string(coordsys) + bin["xref"].real(ra) + bin["yref"].real(dec) + bin["proj"].string(proj) + bin.execute() + sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n") + + # Update timing + wall_seconds += bin.telapse() + cpu_seconds += bin.celapse() + + # Perform maximum likelihood fitting + like = ctools.ctlike() + like.logFileOpen() # We need this to explicitely open the log file in Python mode + like["infile"].filename(cntmap_name) + like["srcmdl"].filename(model_name) + like["outmdl"].filename(result_name) + like["caldb"].string(caldb) + like["irf"].string(irf) + like.execute() + sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n") + + # Update timing + wall_seconds += like.telapse() + cpu_seconds += like.celapse() + + # Show total times + sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n") + sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n") + + # Return + return # ============================= # # Analysis pipeline - version 2 # # ============================= # def pipeline_v2(): - """ - Binned analysis pipeline - keep intermediate results in memory. + """ + Binned analysis pipeline - keep intermediate results in memory. This function implements an analysis pipeline that successively calls ctobssim, ctbin and ctlike without saving the intermediate results as - FITS files on disk. All data are hold in memory. - - At the end, results are plotted (if matplotlib is installed) - """ - # Set script parameters - model_name = "${CTOOLS}/share/models/crab.xml" - caldb = "${CALDB}/data/cta/dummy/bcf" - irf = "cta_dummy_irf" - ra = 83.63 - dec = 22.01 - rad_sim = 10.0 - tstart = 0.0 - tstop = 1800.0 - emin = 0.1 - emax = 100.0 - enumbins = 20 - nxpix = 200 - nypix = 200 - binsz = 0.02 - coordsys = "CEL" - proj = "CAR" - - # 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() - sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n") - - # Update timing - wall_seconds += sim.telapse() - cpu_seconds += sim.celapse() - - # Bin events into counts map - bin = ctbin(sim.obs()) - bin["ebinalg"].string("LOG") - bin["emin"].real(emin) - bin["emax"].real(emax) - bin["enumbins"].integer(enumbins) - bin["nxpix"].integer(nxpix) - bin["nypix"].integer(nypix) - bin["binsz"].real(binsz) - bin["coordsys"].string(coordsys) - bin["xref"].real(ra) - bin["yref"].real(dec) - bin["proj"].string(proj) - bin.run() - sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n") - - # Update timing - wall_seconds += bin.telapse() - cpu_seconds += bin.celapse() - - # Perform maximum likelihood fitting - like = ctlike(bin.obs()) - like.run() - sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n") - - # Update timing - wall_seconds += like.telapse() - cpu_seconds += like.celapse() - - # Show total times - sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n") - sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n") - - # Show model fitting results - #sys.stdout.write(like.obs().models()+"\n") - - # Plot counts - if has_matplotlib: - plot_counts(bin.obs()) - else: - sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n") - - # Return - return + FITS files on disk. All data are hold in memory. + + At the end, results are plotted (if matplotlib is installed) + """ + # Set script parameters + model_name = "${CTOOLS}/share/models/crab.xml" + caldb = "${CALDB}/data/cta/dummy/bcf" + irf = "cta_dummy_irf" + ra = 83.63 + dec = 22.01 + rad_sim = 10.0 + tstart = 0.0 + tstop = 1800.0 + emin = 0.1 + emax = 100.0 + enumbins = 20 + nxpix = 200 + nypix = 200 + binsz = 0.02 + coordsys = "CEL" + proj = "CAR" + + # Initialise timing + wall_seconds = 0.0 + cpu_seconds = 0.0 + + # Simulate events + sim = ctools.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() + sys.stdout.write("Simulated events ("+str(sim.celapse())+" CPU seconds)\n") + + # Update timing + wall_seconds += sim.telapse() + cpu_seconds += sim.celapse() + + # Bin events into counts map + bin = ctools.ctbin(sim.obs()) + bin["ebinalg"].string("LOG") + bin["emin"].real(emin) + bin["emax"].real(emax) + bin["enumbins"].integer(enumbins) + bin["nxpix"].integer(nxpix) + bin["nypix"].integer(nypix) + bin["binsz"].real(binsz) + bin["coordsys"].string(coordsys) + bin["xref"].real(ra) + bin["yref"].real(dec) + bin["proj"].string(proj) + bin.run() + sys.stdout.write("Binned events into counts map ("+str(bin.celapse())+" CPU seconds)\n") + + # Update timing + wall_seconds += bin.telapse() + cpu_seconds += bin.celapse() + + # Perform maximum likelihood fitting + like = ctools.ctlike(bin.obs()) + like.run() + sys.stdout.write("Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)\n") + + # Update timing + wall_seconds += like.telapse() + cpu_seconds += like.celapse() + + # Show total times + sys.stdout.write("Total wall time elapsed: "+str(wall_seconds)+" seconds\n") + sys.stdout.write("Total CPU time used ...: "+str(cpu_seconds)+" seconds\n") + + # Show model fitting results + #sys.stdout.write(like.obs().models()+"\n") + + # Plot counts + if has_matplotlib: + plot_counts(bin.obs()) + else: + sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n") + + # Return + return # =========== # # Plot counts # # =========== # def plot_counts(observations): - """ - Plot counts. - """ - # Set legend fontsize - params = {'legend.fontsize': 10} - plt.rcParams.update(params) - - # Set plot styles - styles = ['b-', 'g-', 'y-', 'n-'] - - # Dump header - sys.stdout.write("\n") - sys.stdout.write("Make plots (using matplotlib):\n") - sys.stdout.write("==============================\n") - - # 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 cube - cube = obs.events() - - # Create energy axis - energy = [] - ebounds = cube.ebounds() - for i in range(ebounds.size()): - energy.append(ebounds.elogmean(i).TeV()) - - # Create spectrum - sys.stdout.write("Extract data:\n") - counts = [0.0 for i in range(ebounds.size())] - for bin in cube: - index = ebounds.index(bin.energy()) - counts[index] = counts[index] + bin.counts() - - # 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 - sys.stdout.write("Extract models:\n") - sum_model = [0.0 for i in range(ebounds.size())] - for k, m in enumerate(observations.models()): - sys.stdout.write("- "+m.name()+"\n") - model = [0.0 for i in range(ebounds.size())] - for bin in cube: - index = ebounds.index(bin.energy()) - prob = m.eval(bin, obs) - model[index] = model[index] + prob * bin.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 - cube = obs.events() - - # Create offset histogram - sys.stdout.write("Extract data:\n") - nx = 30 - doffset2 = 0.01 - offset2 = [(i+0.5)*doffset2 for i in range(nx)] - counts = [0.0 for i in range(nx)] - for bin in cube: - off = GCTAInstDir(bin.dir()).dir().dist_deg(crab) - off2 = off*off - index = int(off2/doffset2) - if index < nx: - counts[index] = counts[index] + bin.counts() - - # Create error bars - error = [sqrt(c) for c in counts] - - # Plot distribution - plt.semilogy(offset2, counts, 'ro', label='data') - - # Extract models - sys.stdout.write("Extract models:\n") - sum_model = [0.0 for i in range(nx)] - for k, m in enumerate(observations.models()): - sys.stdout.write("- "+m.name()+"\n") - model = [0.0 for i in range(nx)] - for bin in cube: - off = GCTAInstDir(bin.dir()).dir().dist_deg(crab) - off2 = off*off - index = int(off2/doffset2) - if index < nx: - prob = m.eval(bin, obs) - model[index] = model[index] + prob * bin.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() - - # Return - return + """ + Plot counts. + """ + # Set legend fontsize + params = {'legend.fontsize': 10} + plt.rcParams.update(params) + + # Set plot styles + styles = ['b-', 'g-', 'y-', 'n-'] + + # Dump header + sys.stdout.write("\n") + sys.stdout.write("Make plots (using matplotlib):\n") + sys.stdout.write("==============================\n") + + # 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 cube + cube = obs.events() + + # Create energy axis + energy = [] + ebounds = cube.ebounds() + for i in range(ebounds.size()): + energy.append(ebounds.elogmean(i).TeV()) + + # Create spectrum + sys.stdout.write("Extract data:\n") + counts = [0.0 for i in range(ebounds.size())] + for bin in cube: + index = ebounds.index(bin.energy()) + counts[index] = counts[index] + bin.counts() + + # Create error bars + error = [math.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 + sys.stdout.write("Extract models:\n") + sum_model = [0.0 for i in range(ebounds.size())] + for k, m in enumerate(observations.models()): + sys.stdout.write("- "+m.name()+"\n") + model = [0.0 for i in range(ebounds.size())] + for bin in cube: + index = ebounds.index(bin.energy()) + prob = m.eval(bin, obs) + model[index] = model[index] + prob * bin.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 = gammalib.GSkyDir() + crab.radec_deg(83.63, 22.01) + + # Loop over observations + for obs in observations: + + # Get event cube + cube = obs.events() + + # Create offset histogram + sys.stdout.write("Extract data:\n") + nx = 30 + doffset2 = 0.01 + offset2 = [(i+0.5)*doffset2 for i in range(nx)] + counts = [0.0 for i in range(nx)] + for bin in cube: + off = gammalib.GCTAInstDir(bin.dir()).dir().dist_deg(crab) + off2 = off*off + index = int(off2/doffset2) + if index < nx: + counts[index] = counts[index] + bin.counts() + + # Create error bars + error = [math.sqrt(c) for c in counts] + + # Plot distribution + plt.semilogy(offset2, counts, 'ro', label='data') + + # Extract models + sys.stdout.write("Extract models:\n") + sum_model = [0.0 for i in range(nx)] + for k, m in enumerate(observations.models()): + sys.stdout.write("- "+m.name()+"\n") + model = [0.0 for i in range(nx)] + for bin in cube: + off = gammalib.GCTAInstDir(bin.dir()).dir().dist_deg(crab) + off2 = off*off + index = int(off2/doffset2) + if index < nx: + prob = m.eval(bin, obs) + model[index] = model[index] + prob * bin.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() + + # Return + return #==========================# # Main routine entry point # #==========================# if __name__ == '__main__': - """ - CTA binned analysis script. - - This script implements the analysis steps to perform a binned analysis. - The steps include: - - Event simulation - - Event binning - - 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. - - The second variation will also show a counts spectrum of the data after - fitting (only if matplotlib is available on your system). Note that the - counts spectrum is summed over the entire data space, i.e. no event selection - around the Crab has been performed. - """ - # Dump header - sys.stdout.write("*************************************\n") - sys.stdout.write("* CTA binned analysis scripts *\n") - sys.stdout.write("*************************************\n") - sys.stdout.write("... this script will take 1-2 minutes\n") - - # 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 - sys.stdout.write("\n") - sys.stdout.write("Executable analysis pipeline:\n") - sys.stdout.write("=============================\n") - pipeline_v1() - - # Do analysis in memory - sys.stdout.write("\n") - sys.stdout.write("In memory analysis pipeline:\n") - sys.stdout.write("============================\n") - pipeline_v2() - sys.stdout.write("\n") + """ + CTA binned analysis script. + + This script implements the analysis steps to perform a binned analysis. + The steps include: + - Event simulation + - Event binning + - 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. + + The second variation will also show a counts spectrum of the data after + fitting (only if matplotlib is available on your system). Note that the + counts spectrum is summed over the entire data space, i.e. no event selection + around the Crab has been performed. + """ + # Dump header + sys.stdout.write("*************************************\n") + sys.stdout.write("* CTA binned analysis scripts *\n") + sys.stdout.write("*************************************\n") + sys.stdout.write("... this script will take 1-2 minutes\n") + + # 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 + sys.stdout.write("\n") + sys.stdout.write("Executable analysis pipeline:\n") + sys.stdout.write("=============================\n") + pipeline_v1() + + # Do analysis in memory + sys.stdout.write("\n") + sys.stdout.write("In memory analysis pipeline:\n") + sys.stdout.write("============================\n") + pipeline_v2() + sys.stdout.write("\n") diff --git a/examples/make_unbinned_analysis.py b/examples/make_unbinned_analysis.py index ea723f7..5011035 100755 --- a/examples/make_unbinned_analysis.py +++ b/examples/make_unbinned_analysis.py @@ -7,7 +7,7 @@ # Required 3rd party modules: # None # -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2014 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 @@ -23,399 +23,402 @@ # along with this program. If not, see . # # ========================================================================== -from ctools import * -from gammalib import * -from math import * +#from ctools import * +#from gammalib import * +#from math import * +import gammalib +import ctools +import math import os import glob import sys try: - import matplotlib.pyplot as plt - has_matplotlib = True + import matplotlib.pyplot as plt + has_matplotlib = True except: - has_matplotlib = False + has_matplotlib = False # ============================= # # Analysis pipeline - version 1 # # ============================= # def pipeline_v1(): - """ - Unbinned analysis pipeline - save intermediate results. + """ + 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 = "${CTOOLS}/share/models/crab.xml" - events_name = "events.fits" - selected_events_name = "selected_events.fits" - result_name = "results.xml" - caldb = "${CALDB}/data/cta/dummy/bcf" - irf = "cta_dummy_irf" - 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["infile"].filename(selected_events_name) - like["srcmdl"].filename(model_name) - like["outmdl"].filename(result_name) - 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 + """ + # Set script parameters + model_name = "${CTOOLS}/share/models/crab.xml" + events_name = "events.fits" + selected_events_name = "selected_events.fits" + result_name = "results.xml" + caldb = "${CALDB}/data/cta/dummy/bcf" + irf = "cta_dummy_irf" + 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 = ctools.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 = ctools.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 = ctools.ctlike() + like.logFileOpen() + like["infile"].filename(selected_events_name) + like["srcmdl"].filename(model_name) + like["outmdl"].filename(result_name) + 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. + """ + 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 = "${CTOOLS}/share/models/crab.xml" - caldb = "${CALDB}/data/cta/dummy/bcf" - irf = "cta_dummy_irf" - 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: - if has_matplotlib: - plot_counts(like.obs()) - else: - sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n") - - # Return - return + FITS files on disk. All data is only hold in memory. + """ + # Set script parameters + model_name = "${CTOOLS}/share/models/crab.xml" + caldb = "${CALDB}/data/cta/dummy/bcf" + irf = "cta_dummy_irf" + 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 = ctools.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 = ctools.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 = ctools.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: + if has_matplotlib: + plot_counts(like.obs()) + else: + sys.stdout.write("Matplotlib is not (correctly) installed on your system. No counts spectra are shown.\n") + + # Return + return # =========== # # Plot counts # # =========== # def plot_counts(observations): - """ - Plot counts. - """ - # 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 = obs.events() - - # Create energy axis - ebounds = GEbounds(10, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV")) - 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 = 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 = GCTAInstDir(atom.dir()).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 = GCTAInstDir(atom.dir()).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() - - # Return - return + """ + Plot counts. + """ + # 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 = obs.events() + + # Create energy axis + ebounds = gammalib.GEbounds(10, gammalib.GEnergy(0.1, "TeV"), gammalib.GEnergy(100.0, "TeV")) + 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 = [math.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 = gammalib.GSkyDir() + crab.radec_deg(83.63, 22.01) + + # Loop over observations + for obs in observations: + + # Get event cube + list = 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 = gammalib.GCTAInstDir(atom.dir()).dir().dist_deg(crab) + off2 = off*off + index = int(off2/doffset2) + if index < nx: + counts[index] = counts[index] + 1.0 + + # Create error bars + error = [math.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 = gammalib.GCTAInstDir(atom.dir()).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() + + # 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 "" + """ + 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 "" -- libgit2 0.21.4