From 6de368966bca1957a75cfa89c9a938bf364b8df6 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 6 Feb 2015 18:31:02 +0100 Subject: [PATCH] Update example scripts. --- examples/Makefile.am | 17 ++++++++++++----- examples/README | 55 ++++++++++++++++++++++--------------------------------- examples/benchmark_cta_analysis.py | 4 ++-- examples/make_binned_analysis.py | 423 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- examples/make_cube_analysis.py | 145 ------------------------------------------------------------------------------------------------------------------------------------------------- examples/make_unbinned_analysis.py | 418 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 6 files changed, 36 insertions(+), 1026 deletions(-) delete mode 100755 examples/make_binned_analysis.py delete mode 100755 examples/make_cube_analysis.py delete mode 100755 examples/make_unbinned_analysis.py diff --git a/examples/Makefile.am b/examples/Makefile.am index 0e65ac1..813f68a 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -2,14 +2,21 @@ # Scripts to be distributed scripts = README \ - pipeline_stacked_mem.py \ - make_binned_analysis.py \ + benchmark_cta_analysis.py \ make_pull_at_sensitivity_limit.py \ + make_spectrum.py \ make_ts_distributions.py \ - make_unbinned_analysis.py \ + pipeline_binned_disk.py \ + pipeline_binned_mem.py \ + pipeline_stacked_disk.py \ + pipeline_stacked_mem.py \ + pipeline_unbinned_disk.py \ + pipeline_unbinned_mem.py \ + show_butterfly.py \ + show_pha.py \ + show_response.py \ show_pull_histogram.py \ - show_ts_distribution.py \ - show_butterfly.py + show_ts_distribution.py # Set directories scriptsdir = $(datadir)/examples/python diff --git a/examples/README b/examples/README index 5fae2de..47bde63 100644 --- a/examples/README +++ b/examples/README @@ -2,6 +2,28 @@ This folder contains Python example scripts that illustrate the usage of the ctools. benchmark_cta_analysis.py + This script peforms a benchmark of the various CTA analysis types + (unbinned, binned, stacked), and if matplotlib is installed, creates + a plot of the benchmark results. + +make_pull_at_sensitivity_limit.py + This script creates pull distributions as function of energy + for a test source that is at the 5 sigma sensitivity limit for + 50h of effective observing time. It takes as input file the + output file generated using cspull. + If the Python processing module is installed the script can make + use of multiple cores and runs in parallel. + +make_spectrum.py + This script shows how to make a spectrum using obsutils. + +make_ts_distributions.py + This script creates TS distributions as function of energy for + an off observation of 50h effective observing time. + If the Python processing module is installed the script can make + use of multiple cores and runs in parallel. + +benchmark_cta_analysis.py Performs performance benchmarks of the various CTA analysis types. pipeline_binned_disk.py @@ -28,39 +50,6 @@ pipeline_unbinned_mem.py Example pipeline for simulating events and performing an unbinned analysis in memory. No files will be written to disk. - - - -make_binned_analysis.py - Example Python pipelines for binned analysis. - The first pipeline writes intermediate results to the disk while - the second pipeline keeps the results in memory. - If matplotlib is installed on the system the fitted spectrum - and the observed offset angle distribution is displayed in - counts space. - -make_pull_at_sensitivity_limit.py - This script creates pull distributions as function of energy - for a test source that is at the 5 sigma sensitivity limit for - 50h of effective observing time. It takes as input file the - output file generated using cspull. - If the Python processing module is installed the script can make - use of multiple cores and runs in parallel. - -make_ts_distributions.py - This script creates TS distributions as function of energy for - an off observation of 50h effective observing time. - If the Python processing module is installed the script can make - use of multiple cores and runs in parallel. - -make_unbinned_analysis.py - Example Python pipelines for unbinned analysis. - The first pipeline writes intermediate results to the disk while - the second pipeline keeps the results in memory. - If matplotlib is installed on the system the fitted spectrum - and the observed offset angle distribution is displayed in - counts space. - show_butterfly.py Displays a butterfly generated using the ctbutterfly tool. diff --git a/examples/benchmark_cta_analysis.py b/examples/benchmark_cta_analysis.py index 4f3eca9..c073a84 100755 --- a/examples/benchmark_cta_analysis.py +++ b/examples/benchmark_cta_analysis.py @@ -1,10 +1,10 @@ #! /usr/bin/env python # ========================================================================== # This script peforms a benchmark of the various CTA analysis types -# (unbinned, binned, cube-style), and if matplotlib is installed, creates +# (unbinned, binned, stacked), and if matplotlib is installed, creates # a plot of the benchmark results. # -# Copyright (C) 2014 Jurgen Knodlseder +# Copyright (C) 2014-2015 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 diff --git a/examples/make_binned_analysis.py b/examples/make_binned_analysis.py deleted file mode 100755 index cd1d174..0000000 --- a/examples/make_binned_analysis.py +++ /dev/null @@ -1,423 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script illustrates how to build a binned analysis pipeline using the -# ctools. If matplotlib is installed, a counts spectrum and an offset -# histogram will be displayed. -# -# 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 -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -import gammalib -import ctools -import os -import glob -import sys -import math -try: - import matplotlib.pyplot as plt - has_matplotlib = True -except: - has_matplotlib = False - - -# ============================= # -# Analysis pipeline - version 1 # -# ============================= # -def pipeline_v1(): - """ - 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 = ctools.ctobssim() - sim.logFileOpen() # We need this to explicitely open the log file in Python mode - sim["inmodel"].filename(model_name) - sim["outevents"].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["inobs"].filename(events_name) - bin["outcube"].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["inobs"].filename(cntmap_name) - like["inmodel"].filename(model_name) - like["outmodel"].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. - - 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 = ctools.ctobssim() - sim["inmodel"].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 = [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") diff --git a/examples/make_cube_analysis.py b/examples/make_cube_analysis.py deleted file mode 100755 index f5d104c..0000000 --- a/examples/make_cube_analysis.py +++ /dev/null @@ -1,145 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script illustrates how to perform a stacked cube analysis. -# -# Copyright (C) 2014 Juergen Knoedlseder -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -import ctools -import gammalib -from ctools import obsutils - - -# ================= # -# Analysis pipeline # -# ================= # -def pipeline(caldb = "dummy", \ - irf = "cta_dummy_irf", \ - model = "$CTOOLS/share/models/crab.xml"): - """ - Analysis pipeline for a stacked cube analysis. The pipeline writes - intermediate results to FITS files on disk. Also log files of all - individual ctools will be written. - - Keywords: - caldb - Calibration database repository - irf - Name of instrument response - model - Source model - """ - # Set pointing pattern (4 pointings) - pointings = obsutils.set_obs_patterns("four", offset=2.0) - - # Set CTA observations (default: 1800 sec per pointing) - obs = obsutils.set_obs_list(pointings, caldb=caldb, irf=irf) - - # Load models - models = gammalib.GModels(model) - obs.models(models) - - # Run ctobssim (simulates 4 event lists) - sim = ctools.ctobssim(obs) - sim.logFileOpen() - sim.run() - - # Run ctbin (builds stacked event cube) - bin = ctools.ctbin(sim.obs()) - bin.logFileOpen() - bin["outcube"].filename("cntmap.fits") - bin["ebinalg"].string("LOG") - bin["emin"].real(0.1) - bin["emax"].real(100.0) - bin["enumbins"].integer(20) - bin["usepnt"].boolean(False) - bin["nxpix"].integer(300) - bin["nypix"].integer(300) - bin["binsz"].real(0.02) - bin["coordsys"].string("CEL") - bin["proj"].string("CAR") - bin["xref"].real(83.6331) - bin["yref"].real(22.0145) - bin.execute() - - # Run ctexpcube (builds exposure cube) - expcube = ctools.ctexpcube(sim.obs()) - expcube.logFileOpen() - expcube["incube"].filename("cntmap.fits") - expcube["outcube"].filename("expcube.fits") - expcube["caldb"].string(caldb) - expcube["irf"].string(irf) - expcube.execute() - - # Run ctpsfcube (builds PSF cube) - psfcube = ctools.ctpsfcube(sim.obs()) - psfcube.logFileOpen() - psfcube["incube"].filename("NONE") - psfcube["outcube"].filename("psfcube.fits") - psfcube["caldb"].string(caldb) - psfcube["irf"].string(irf) - psfcube["ebinalg"].string("LOG") - psfcube["emin"].real(0.1) - psfcube["emax"].real(100.0) - psfcube["enumbins"].integer(20) - psfcube["nxpix"].integer(5) - psfcube["nypix"].integer(5) - psfcube["binsz"].real(2.0) - psfcube["coordsys"].string("CEL") - psfcube["proj"].string("CAR") - psfcube["xref"].real(83.6331) - psfcube["yref"].real(22.0145) - psfcube["amax"].real(0.3) - psfcube["anumbins"].integer(60) - psfcube.execute() - - # Run ctbkgcube (builds background cube) - bkgcube = ctools.ctbkgcube(sim.obs()) - bkgcube.logFileOpen() - bkgcube["incube"].filename("cntmap.fits") - bkgcube["inmodel"].filename(model) - bkgcube["outcube"].filename("bkgcube.fits") - bkgcube["outmodel"].filename("bkgcube.xml") - bkgcube.execute() - - # Create observation container for stacked cube analysis - obs.clear() - cta = gammalib.GCTAObservation("cntmap.fits", "expcube.fits", "psfcube.fits") - obs.append(cta) - obs.models(gammalib.GModels("bkgcube.xml")) - - # Run ctlike - like = ctools.ctlike(obs) - like.logFileOpen() - like["debug"].boolean(True) - like["outmodel"].filename("cube_results.xml") - like.execute() - - # Return - return - - -#==========================# -# Main routine entry point # -#==========================# -if __name__ == '__main__': - """ - """ - # Dump header - print("*************************************") - print("* CTA stacked cube analysis *") - print("*************************************") - - # Run pipeline - pipeline() - \ No newline at end of file diff --git a/examples/make_unbinned_analysis.py b/examples/make_unbinned_analysis.py deleted file mode 100755 index 3859571..0000000 --- a/examples/make_unbinned_analysis.py +++ /dev/null @@ -1,418 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script illustrates how to build an unbinned analysis pipeline with -# the ctools. If matplotlib is installed, the spectrum will be displayed -# on the screen. -# -# Required 3rd party modules: -# None -# -# Copyright (C) 2011-2015 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 . -# -# ========================================================================== -import gammalib -import ctools -import math -import os -import glob -import sys -try: - import matplotlib.pyplot as plt - has_matplotlib = True -except: - has_matplotlib = False - - -# ============================= # -# 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 = "${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["inmodel"].filename(model_name) - sim["outevents"].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["inobs"].filename(events_name) - select["outobs"].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["inobs"].filename(selected_events_name) - like["inmodel"].filename(model_name) - like["outmodel"].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. - - 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 = ctools.ctobssim() - sim["inmodel"].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") - - # 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 = 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 "" -- libgit2 0.21.4