From f59f73ff10ef59393f84860ff6d88beb02884dc1 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 24 Jun 2016 22:37:11 +0200 Subject: [PATCH] Move some example scripts in test/dev and add Prod3 IRF generation to unit tests --- examples/benchmark_cta_analysis.py | 361 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- examples/check_models.py | 128 -------------------------------------------------------------------------------------------------------------------------------- examples/generate_prod3_irfs.py | 61 +++++++++++++++++++++++++++++++++++++++++++++---------------- examples/simulate_events.py | 267 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- test/dev/benchmark_cta_analysis.py | 361 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/dev/check_models.py | 128 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/dev/simulate_events.py | 267 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/test_examples.py | 14 ++++++++++++++ 8 files changed, 815 insertions(+), 772 deletions(-) delete mode 100755 examples/benchmark_cta_analysis.py delete mode 100755 examples/check_models.py delete mode 100755 examples/simulate_events.py create mode 100755 test/dev/benchmark_cta_analysis.py create mode 100755 test/dev/check_models.py create mode 100755 test/dev/simulate_events.py diff --git a/examples/benchmark_cta_analysis.py b/examples/benchmark_cta_analysis.py deleted file mode 100755 index 5fa01d4..0000000 --- a/examples/benchmark_cta_analysis.py +++ /dev/null @@ -1,361 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# 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. -# -# Copyright (C) 2014-2016 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 ctools -import time -try: - import matplotlib.pyplot as plt - has_matplotlib = True -except: - has_matplotlib = False - - -# ========================== # -# Unbinned analysis pipeline # -# ========================== # -def unbinned_pipeline(model_name, duration): - """ - Unbinned analysis pipeline. - """ - # Set script parameters - caldb = "prod2" - irf = "South_50h" - ra = 83.63 - dec = 22.01 - rad_sim = 10.0 - tstart = 0.0 - tstop = duration - emin = 0.1 - emax = 100.0 - rad_select = 3.0 - - # Get start CPU time - cpu_start = time.clock() - - # Simulate events - sim = ctools.ctobssim() - sim["inmodel"] = model_name - sim["caldb"] = caldb - sim["irf"] = irf - sim["ra"] = ra - sim["dec"] = dec - sim["rad"] = rad_sim - sim["tmin"] = tstart - sim["tmax"] = tstop - sim["emin"] = emin - sim["emax"] = emax - sim.run() - - # Select events - select = ctools.ctselect(sim.obs()) - select["ra"] = ra - select["dec"] = dec - select["rad"] = rad_select - select["tmin"] = tstart - select["tmax"] = tstop - select["emin"] = emin - select["emax"] = emax - select.run() - - # Get ctlike start CPU time - cpu_ctlike = time.clock() - - # Perform maximum likelihood fitting - like = ctools.ctlike(select.obs()) - like.run() - - # Get stop CPU time and compute elapsed times - cpu_stop = time.clock() - cpu_elapsed = cpu_stop - cpu_start - cpu_ctlike = cpu_stop - cpu_ctlike - - # Return - return cpu_elapsed, cpu_ctlike - - -# ======================== # -# Binned analysis pipeline # -# ======================== # -def binned_pipeline(model_name, duration): - """ - Binned analysis pipeline. - """ - # Set script parameters - caldb = "prod2" - irf = "South_50h" - ra = 83.63 - dec = 22.01 - rad_sim = 10.0 - tstart = 0.0 - tstop = duration - emin = 0.1 - emax = 100.0 - enumbins = 40 - nxpix = 200 - nypix = 200 - binsz = 0.02 - coordsys = "CEL" - proj = "CAR" - - # Get start CPU time - cpu_start = time.clock() - - # Simulate events - sim = ctools.ctobssim() - sim["inmodel"] = model_name - sim["caldb"] = caldb - sim["irf"] = irf - sim["ra"] = ra - sim["dec"] = dec - sim["rad"] = rad_sim - sim["tmin"] = tstart - sim["tmax"] = tstop - sim["emin"] = emin - sim["emax"] = emax - sim.run() - - # Bin events into counts map - bin = ctools.ctbin(sim.obs()) - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax - bin["enumbins"] = enumbins - bin["nxpix"] = nxpix - bin["nypix"] = nypix - bin["binsz"] = binsz - bin["coordsys"] = coordsys - bin["xref"] = ra - bin["yref"] = dec - bin["proj"] = proj - bin.run() - - # Get ctlike start CPU time - cpu_ctlike = time.clock() - - # Perform maximum likelihood fitting - like = ctools.ctlike(bin.obs()) - like.run() - - # Get stop CPU time and compute elapsed times - cpu_stop = time.clock() - cpu_elapsed = cpu_stop - cpu_start - cpu_ctlike = cpu_stop - cpu_ctlike - - # Return - return cpu_elapsed, cpu_ctlike - - -# ========================= # -# Stacked analysis pipeline # -# ========================= # -def stacked_pipeline(model_name, duration): - """ - Stacked analysis pipeline. - """ - # Set script parameters - caldb = "prod2" - irf = "South_50h" - ra = 83.63 - dec = 22.01 - rad_sim = 10.0 - tstart = 0.0 - tstop = duration - emin = 0.1 - emax = 100.0 - enumbins = 40 - nxpix = 200 - nypix = 200 - binsz = 0.02 - coordsys = "CEL" - proj = "CAR" - - # Get start CPU time - cpu_start = time.clock() - - # Simulate events - sim = ctools.ctobssim() - sim["inmodel"] = model_name - sim["caldb"] = caldb - sim["irf"] = irf - sim["ra"] = ra - sim["dec"] = dec - sim["rad"] = rad_sim - sim["tmin"] = tstart - sim["tmax"] = tstop - sim["emin"] = emin - sim["emax"] = emax - sim.run() - - # Bin events into counts map - bin = ctools.ctbin(sim.obs()) - bin["ebinalg"] = "LOG" - bin["emin"] = emin - bin["emax"] = emax - bin["enumbins"] = enumbins - bin["nxpix"] = nxpix - bin["nypix"] = nypix - bin["binsz"] = binsz - bin["coordsys"] = coordsys - bin["proj"] = proj - bin["xref"] = ra - bin["yref"] = dec - bin.run() - - # Create exposure cube - expcube = ctools.ctexpcube(sim.obs()) - expcube["incube"] = "NONE" - expcube["caldb"] = caldb - expcube["irf"] = irf - expcube["ebinalg"] = "LOG" - expcube["emin"] = emin - expcube["emax"] = emax - expcube["enumbins"] = enumbins - expcube["nxpix"] = nxpix - expcube["nypix"] = nypix - expcube["binsz"] = binsz - expcube["coordsys"] = coordsys - expcube["proj"] = proj - expcube["xref"] = ra - expcube["yref"] = dec - expcube.run() - - # Create PSF cube - psfcube = ctools.ctpsfcube(sim.obs()) - psfcube["incube"] = "NONE" - psfcube["caldb"] = caldb - psfcube["irf"] = irf - psfcube["ebinalg"] = "LOG" - psfcube["emin"] = emin - psfcube["emax"] = emax - psfcube["enumbins"] = enumbins - psfcube["nxpix"] = 10 - psfcube["nypix"] = 10 - psfcube["binsz"] = 1.0 - psfcube["coordsys"] = coordsys - psfcube["proj"] = proj - psfcube["xref"] = ra - psfcube["yref"] = dec - psfcube.run() - - # Create background cube - bkgcube = ctools.ctbkgcube(sim.obs()) - bkgcube["incube"] = "NONE" - bkgcube["ebinalg"] = "LOG" - bkgcube["emin"] = emin - bkgcube["emax"] = emax - bkgcube["enumbins"] = enumbins - bkgcube["nxpix"] = 10 - bkgcube["nypix"] = 10 - bkgcube["binsz"] = 1.0 - bkgcube["coordsys"] = coordsys - bkgcube["proj"] = proj - bkgcube["xref"] = ra - bkgcube["yref"] = dec - bkgcube.run() - - # Attach background model to observation container - bin.obs().models(bkgcube.models()) - - # Set Exposure and Psf cube for first CTA observation - # (ctbin will create an observation with a single container) - bin.obs()[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube()) - - # Get ctlike start CPU time - cpu_ctlike = time.clock() - - # Perform maximum likelihood fitting - like = ctools.ctlike(bin.obs()) - like.run() - - # Get stop CPU time and compute elapsed times - cpu_stop = time.clock() - cpu_elapsed = cpu_stop - cpu_start - cpu_ctlike = cpu_stop - cpu_ctlike - - # Return - return cpu_elapsed, cpu_ctlike - - -#==========================# -# Main routine entry point # -#==========================# -if __name__ == '__main__': - - # Initialise parameters - title = "CTA analysis benchmark" - model_name = "${CTOOLS}/share/models/crab.xml" - - # Log header - print("*************************************") - print("* CTA analysis benchmark *") - print("*************************************") - print(" Duration Unbinned Binned Stacked ctlike-UB ctlike-B ctlike-S") - - # Initialize arrays - a_duration = [] - a_unbinned = [] - a_binned = [] - a_stacked = [] - c_unbinned = [] - c_binned = [] - c_stacked = [] - - # Perform benchmarks - duration = 1800.0 - for i in range(10): - - # Perform analyses - t_unbinned, ct_unbinned = unbinned_pipeline(model_name, duration) - t_binned, ct_binned = binned_pipeline(model_name, duration) - t_stacked, ct_stacked = stacked_pipeline(model_name, duration) - - # Collect results - a_duration.append(duration/3600.0) - a_unbinned.append(t_unbinned) - a_binned.append(t_binned) - a_stacked.append(t_stacked) - c_unbinned.append(ct_unbinned) - c_binned.append(ct_binned) - c_stacked.append(ct_stacked) - - # Print results - print("%10.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f" % \ - (duration, t_unbinned, t_binned, t_stacked, \ - ct_unbinned, ct_binned, ct_stacked)) - - # Increment duration - duration *= 2.0 - - # Optionally plot results - if has_matplotlib: - plt.figure(1) - plt.title(title) - plt.loglog(a_duration, a_unbinned, 'ro-', label='unbinned') - plt.loglog(a_duration, a_binned, 'bo-', label='binned') - plt.loglog(a_duration, a_stacked, 'go-', label='stacked') - plt.loglog(a_duration, c_unbinned, 'ro--') - plt.loglog(a_duration, c_binned, 'bo--') - plt.loglog(a_duration, c_stacked, 'go--') - plt.xlabel("Observation duration (hours)") - plt.ylabel("CPU time (seconds)") - plt.legend(loc="lower right") - plt.show() diff --git a/examples/check_models.py b/examples/check_models.py deleted file mode 100755 index 2a4902b..0000000 --- a/examples/check_models.py +++ /dev/null @@ -1,128 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# Performs unbinned maximum likelihood analysis for various models -# -# Copyright (C) 2011-2016 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 sys -import ctools - - -# ================= # -# Analysis pipeline # -# ================= # -def pipeline(model_name): - """ - Unbinned in-memory analysis pipeline - - 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 - caldb = 'prod2' - irf = 'South_0.5h' - ra = 83.63 - dec = 22.01 - rad_sim = 5.0 - tstart = 0.0 - tstop = 180.0 - emin = 0.1 - emax = 100.0 - rad_select = 3.0 - - # Write model name - print('*** Model: '+model_name+' ************************************') - - # Initialise timing - wall_seconds = 0.0 - cpu_seconds = 0.0 - - # Simulate events - sim = ctools.ctobssim() - sim['inmodel'] = model_name - sim['caldb'] = caldb - sim['irf'] = irf - sim['ra'] = ra - sim['dec'] = dec - sim['rad'] = rad_sim - sim['tmin'] = tstart - sim['tmax'] = tstop - sim['emin'] = emin - sim['emax'] = emax - sim.run() - - # Select events - select = ctools.ctselect(sim.obs()) - select['ra'] = ra - select['dec'] = dec - select['rad'] = rad_select - select['tmin'] = tstart - select['tmax'] = tstop - select['emin'] = emin - select['emax'] = emax - select.run() - - # Perform maximum likelihood fitting - like = ctools.ctlike(select.obs()) - like.run() - - # Show model fitting results - print(like.obs().models()[0]) - - # Return - return - - -#==========================# -# Main routine entry point # -#==========================# -if __name__ == '__main__': - - # Initialise flags - need_help = False - - # Test for command line arguments - print(sys.argv[0]) - if (len(sys.argv) > 1): - if sys.argv[1] == '-h': - need_help = True - else: - need_help = True - - # Print help if needed and exit - if need_help: - print('Usage: check_models.py [OPTIONS]') - print(' -h Display this usage message') - sys.exit() - - # Dump header - print('*******************************************') - print('* Check models using an unbinned analysis *') - print('*******************************************') - - # Perform analysis for Crab model - pipeline('$CTOOLS/share/models/crab.xml') - - # Perform analysis for disk model - pipeline('$CTOOLS/share/models/disk.xml') - - # Perform analysis for Gaussian model - pipeline('$CTOOLS/share/models/gauss.xml') - - # Perform analysis for shell model - pipeline('$CTOOLS/share/models/shell.xml') diff --git a/examples/generate_prod3_irfs.py b/examples/generate_prod3_irfs.py index fbb263d..15e4a91 100755 --- a/examples/generate_prod3_irfs.py +++ b/examples/generate_prod3_irfs.py @@ -30,12 +30,21 @@ def set_prod3_desy(): Set Prod3 DESY database. """ # Set database attributes - path = "/project-data/cta/performance/prod3" - psftype = "Gauss" + #path = '/project-data/cta/performance/prod3' + #file = 'DESY.d20160304.V5.ID0NIM2.prod3-paranalp05-NN.S.3HB1-ND-3.180000s.root' + path = 'data' + file = 'irf.root' + psftype = 'Gauss' + oversample = 3 + norm1d = True + bgdinfill = False + bgdethres = 10.0 # Set database content - db = [{'inst': "prod3", 'id': "South_50h", 'psftype': psftype, 'path': path, - 'file': "DESY.d20160304.V5.ID0NIM2.prod3-paranalp05-NN.S.3HB1-ND-3.180000s.root"} + db = [{'inst': 'prod3', 'id': 'South_50h', 'psftype': psftype, + 'oversample': oversample, 'norm1d': norm1d, + 'bgdinfill': bgdinfill, 'bgdethres': bgdethres, + 'path': path, 'file': file} ] # Return database @@ -66,23 +75,43 @@ if __name__ == '__main__': bgdscale = entry['bgdscale'] else: bgdscale = 1.0 + if entry.has_key('bgdethres'): + bgdethres = entry['bgdethres'] + else: + bgdethres = 1000.0 + if entry.has_key('bgdinfill'): + bgdinfill = entry['bgdinfill'] + else: + bgdinfill = False + if entry.has_key('oversample'): + oversample = entry['oversample'] + else: + oversample = 1 + if entry.has_key('norm1d'): + norm1d = entry['norm1d'] + else: + norm1d = False # Set filename - filename = entry['path']+"/"+entry['file'] + filename = entry['path']+'/'+entry['file'] # Set-up csroot2caldb caldb = cscripts.csroot2caldb() - caldb["infile"] = filename - caldb["outdir"] = "prod3" - caldb["inst"] = entry['inst'] - caldb["id"] = entry['id'] - caldb["version"] = entry['inst'] - caldb["analysis"] = "DESY" - caldb["psftype"] = entry['psftype'] - caldb["rebin"] = rebin - caldb["eascale"] = eascale - caldb["bgdscale"] = bgdscale - caldb["logfile"] = "generate_prod3_irfs.log" + caldb['infile'] = filename + caldb['outdir'] = 'prod3' + caldb['inst'] = entry['inst'] + caldb['id'] = entry['id'] + caldb['version'] = entry['inst'] + caldb['analysis'] = 'DESY' + caldb['psftype'] = entry['psftype'] + caldb['norm1d'] = norm1d + caldb['rebin'] = rebin + caldb['eascale'] = eascale + caldb['bgdscale'] = bgdscale + caldb['bgdoversample'] = oversample + caldb['bgdethres'] = bgdethres + caldb['bgdinfill'] = bgdinfill + caldb['logfile'] = 'generate_prod3_irfs.log' # Add CALDB entry caldb.execute() diff --git a/examples/simulate_events.py b/examples/simulate_events.py deleted file mode 100755 index 01fcac6..0000000 --- a/examples/simulate_events.py +++ /dev/null @@ -1,267 +0,0 @@ -#! /usr/bin/env python -# ========================================================================== -# This script simulates events from a data-space model and compares them to -# the model. -# -# If matplotlib is installed, the event spectrum will be displayed. -# -# -------------------------------------------------------------------------- -# -# Copyright (C) 2014 Juergen Knoedlseder -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . -# -# ========================================================================== -import gammalib -import ctools -import math -try: - import matplotlib.pyplot as plt - has_matplotlib = True -except: - has_matplotlib = False - - -# ======================== # -# Create CTA observation # -# ======================== # -def createobs(ra=86.171648, dec=-1.4774586, rad=5.0, - emin=0.1, emax=100.0, duration=360000.0, deadc=0.95, - irf="South_50h", caldb="prod2"): - """ - Create CTA observation. - """ - # Allocate CTA observation - obs = gammalib.GCTAObservation() - - # Set calibration database - db = gammalib.GCaldb() - if (gammalib.dir_exists(caldb)): - db.rootdir(caldb) - else: - db.open("cta", caldb) - - # Set pointing direction - pntdir = gammalib.GSkyDir() - pntdir.radec_deg(ra, dec) - pnt = gammalib.GCTAPointing() - pnt.dir(pntdir) - obs.pointing(pnt) - - # Set ROI - roi = gammalib.GCTARoi() - instdir = gammalib.GCTAInstDir() - instdir.dir(pntdir) - roi.centre(instdir) - roi.radius(rad) - - # Set GTI - gti = gammalib.GGti() - start = gammalib.GTime(0.0) - stop = gammalib.GTime(duration) - gti.append(start, stop) - - # Set energy boundaries - ebounds = gammalib.GEbounds() - e_min = gammalib.GEnergy() - e_max = gammalib.GEnergy() - e_min.TeV(emin) - e_max.TeV(emax) - ebounds.append(e_min, e_max) - - # Allocate event list - events = gammalib.GCTAEventList() - events.roi(roi) - events.gti(gti) - events.ebounds(ebounds) - obs.events(events) - - # Set instrument response - obs.response(irf, db) - - # Set ontime, livetime, and deadtime correction factor - obs.ontime(duration) - obs.livetime(duration*deadc) - obs.deadc(deadc) - - # Return observation - return obs - - -# =============== # -# Simulate events # -# =============== # -def simulate(obs, xmlname): - """ - Simulate events. - """ - # Allocate random number generator - ran = gammalib.GRan() - - # Load models and extract first model - models = gammalib.GModels(xmlname) - model = models[0] - print(model) - - # Simulate events - events = model.mc(obs, ran) - - # Print event statistics - npred = obs.npred(models) - print(str(len(events)) + " events simulated.") - print(str(npred) + " events expected from Npred.") - - # Check distance - pntdir = obs.pointing().dir() - for event in events: - distance = gammalib.GCTAInstDir(event.dir()).dir().dist_deg(pntdir) - if (distance > 5.0): - print("Too far: "+str(distance)) - - # Return events - return events - - -# ============================== # -# Simulate events using ctobssim # -# ============================== # -def simulate_ctobssim(obs, xmlname, seed=0): - """ - Simulate events using ctobssim. - """ - # Create containers - observations = gammalib.GObservations() - observations.append(obs) - - # Append models - models = gammalib.GModels(xmlname) - observations.models(models) - print(models[0]) - - # Allocate ctobssim application and set parameters - sim = ctools.ctobssim(observations) - sim['seed'].integer(seed) - - # Run simulator - sim.run() - - # Retrieve events - events = sim.obs()[0].events().copy() - - # Print event statistics - npred = obs.npred(models) - print(str(len(events)) + " events simulated.") - print(str(npred) + " events expected from Npred.") - - # Delete the simulation - del sim - - # Return events - return events - - -# =========== # -# Show events # -# =========== # -def show_events(events, xmlname, duration, emin, emax, ebins=30): - """ - Show events using matplotlib. - """ - # Create figure - plt.figure(1) - plt.title("MC simulated event spectrum (" + str(emin) + '-' + str(emax) + " TeV)") - - # Setup energy range covered by data - ebds = gammalib.GEbounds(ebins, gammalib.GEnergy(emin, "TeV"), - gammalib.GEnergy(emax, "TeV")) - - # Create energy axis - energy = [] - for i in range(ebds.size()): - energy.append(ebds.elogmean(i).TeV()) - - # Fill histogram - counts = [0.0 for i in range(ebds.size())] - for event in events: - index = ebds.index(event.energy()) - counts[index] = counts[index] + 1.0 - - # Create error bars - error = [math.sqrt(c) for c in counts] - - # Get model values - sum = 0.0 - models = gammalib.GModels(xmlname) - m = models[0] - model = [] - t = gammalib.GTime() - for i in range(ebds.size()): - eval = ebds.elogmean(i) - ewidth = ebds.emax(i) - ebds.emin(i) - f = m.npred(eval, t, obs) * ewidth.MeV() * duration - sum += f - model.append(f) - print(str(sum) + " events expected from spectrum (integration).") - - # Plot data - plt.loglog(energy, counts, 'ro') - plt.errorbar(energy, counts, error, fmt=None, ecolor='r') - - # Plot model - plt.plot(energy, model, 'b-') - - # Set axes - plt.xlabel("Energy (TeV)") - plt.ylabel("Number of events") - - # Notify - print("PLEASE CLOSE WINDOW TO CONTINUE ...") - - # Show plot - plt.show() - - # Return - return - - -#==========================# -# Main routine entry point # -#==========================# -if __name__ == '__main__': - - # Dump header - print('') - print('*******************') - print('* Simulate events *') - print('*******************') - - # Set XML names - xmlname = 'model_bgd.xml' - - # Set simulation parameters - duration = 360000.0 # Observation duration (s) - emin = 0.1 # 0.1 TeV - emax = 100.0 # 100 TeV - ebins = 30 - - # Set observation - obs = createobs(emin=emin, emax=emax, duration=duration) - - # Perform simulation - #events = simulate(obs, xmlname) - events = simulate_ctobssim(obs, xmlname) - - # Show events - if has_matplotlib: - show_events(events, xmlname, duration, emin, emax, ebins=ebins) diff --git a/test/dev/benchmark_cta_analysis.py b/test/dev/benchmark_cta_analysis.py new file mode 100755 index 0000000..5fa01d4 --- /dev/null +++ b/test/dev/benchmark_cta_analysis.py @@ -0,0 +1,361 @@ +#! /usr/bin/env python +# ========================================================================== +# 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. +# +# Copyright (C) 2014-2016 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 ctools +import time +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except: + has_matplotlib = False + + +# ========================== # +# Unbinned analysis pipeline # +# ========================== # +def unbinned_pipeline(model_name, duration): + """ + Unbinned analysis pipeline. + """ + # Set script parameters + caldb = "prod2" + irf = "South_50h" + ra = 83.63 + dec = 22.01 + rad_sim = 10.0 + tstart = 0.0 + tstop = duration + emin = 0.1 + emax = 100.0 + rad_select = 3.0 + + # Get start CPU time + cpu_start = time.clock() + + # Simulate events + sim = ctools.ctobssim() + sim["inmodel"] = model_name + sim["caldb"] = caldb + sim["irf"] = irf + sim["ra"] = ra + sim["dec"] = dec + sim["rad"] = rad_sim + sim["tmin"] = tstart + sim["tmax"] = tstop + sim["emin"] = emin + sim["emax"] = emax + sim.run() + + # Select events + select = ctools.ctselect(sim.obs()) + select["ra"] = ra + select["dec"] = dec + select["rad"] = rad_select + select["tmin"] = tstart + select["tmax"] = tstop + select["emin"] = emin + select["emax"] = emax + select.run() + + # Get ctlike start CPU time + cpu_ctlike = time.clock() + + # Perform maximum likelihood fitting + like = ctools.ctlike(select.obs()) + like.run() + + # Get stop CPU time and compute elapsed times + cpu_stop = time.clock() + cpu_elapsed = cpu_stop - cpu_start + cpu_ctlike = cpu_stop - cpu_ctlike + + # Return + return cpu_elapsed, cpu_ctlike + + +# ======================== # +# Binned analysis pipeline # +# ======================== # +def binned_pipeline(model_name, duration): + """ + Binned analysis pipeline. + """ + # Set script parameters + caldb = "prod2" + irf = "South_50h" + ra = 83.63 + dec = 22.01 + rad_sim = 10.0 + tstart = 0.0 + tstop = duration + emin = 0.1 + emax = 100.0 + enumbins = 40 + nxpix = 200 + nypix = 200 + binsz = 0.02 + coordsys = "CEL" + proj = "CAR" + + # Get start CPU time + cpu_start = time.clock() + + # Simulate events + sim = ctools.ctobssim() + sim["inmodel"] = model_name + sim["caldb"] = caldb + sim["irf"] = irf + sim["ra"] = ra + sim["dec"] = dec + sim["rad"] = rad_sim + sim["tmin"] = tstart + sim["tmax"] = tstop + sim["emin"] = emin + sim["emax"] = emax + sim.run() + + # Bin events into counts map + bin = ctools.ctbin(sim.obs()) + bin["ebinalg"] = "LOG" + bin["emin"] = emin + bin["emax"] = emax + bin["enumbins"] = enumbins + bin["nxpix"] = nxpix + bin["nypix"] = nypix + bin["binsz"] = binsz + bin["coordsys"] = coordsys + bin["xref"] = ra + bin["yref"] = dec + bin["proj"] = proj + bin.run() + + # Get ctlike start CPU time + cpu_ctlike = time.clock() + + # Perform maximum likelihood fitting + like = ctools.ctlike(bin.obs()) + like.run() + + # Get stop CPU time and compute elapsed times + cpu_stop = time.clock() + cpu_elapsed = cpu_stop - cpu_start + cpu_ctlike = cpu_stop - cpu_ctlike + + # Return + return cpu_elapsed, cpu_ctlike + + +# ========================= # +# Stacked analysis pipeline # +# ========================= # +def stacked_pipeline(model_name, duration): + """ + Stacked analysis pipeline. + """ + # Set script parameters + caldb = "prod2" + irf = "South_50h" + ra = 83.63 + dec = 22.01 + rad_sim = 10.0 + tstart = 0.0 + tstop = duration + emin = 0.1 + emax = 100.0 + enumbins = 40 + nxpix = 200 + nypix = 200 + binsz = 0.02 + coordsys = "CEL" + proj = "CAR" + + # Get start CPU time + cpu_start = time.clock() + + # Simulate events + sim = ctools.ctobssim() + sim["inmodel"] = model_name + sim["caldb"] = caldb + sim["irf"] = irf + sim["ra"] = ra + sim["dec"] = dec + sim["rad"] = rad_sim + sim["tmin"] = tstart + sim["tmax"] = tstop + sim["emin"] = emin + sim["emax"] = emax + sim.run() + + # Bin events into counts map + bin = ctools.ctbin(sim.obs()) + bin["ebinalg"] = "LOG" + bin["emin"] = emin + bin["emax"] = emax + bin["enumbins"] = enumbins + bin["nxpix"] = nxpix + bin["nypix"] = nypix + bin["binsz"] = binsz + bin["coordsys"] = coordsys + bin["proj"] = proj + bin["xref"] = ra + bin["yref"] = dec + bin.run() + + # Create exposure cube + expcube = ctools.ctexpcube(sim.obs()) + expcube["incube"] = "NONE" + expcube["caldb"] = caldb + expcube["irf"] = irf + expcube["ebinalg"] = "LOG" + expcube["emin"] = emin + expcube["emax"] = emax + expcube["enumbins"] = enumbins + expcube["nxpix"] = nxpix + expcube["nypix"] = nypix + expcube["binsz"] = binsz + expcube["coordsys"] = coordsys + expcube["proj"] = proj + expcube["xref"] = ra + expcube["yref"] = dec + expcube.run() + + # Create PSF cube + psfcube = ctools.ctpsfcube(sim.obs()) + psfcube["incube"] = "NONE" + psfcube["caldb"] = caldb + psfcube["irf"] = irf + psfcube["ebinalg"] = "LOG" + psfcube["emin"] = emin + psfcube["emax"] = emax + psfcube["enumbins"] = enumbins + psfcube["nxpix"] = 10 + psfcube["nypix"] = 10 + psfcube["binsz"] = 1.0 + psfcube["coordsys"] = coordsys + psfcube["proj"] = proj + psfcube["xref"] = ra + psfcube["yref"] = dec + psfcube.run() + + # Create background cube + bkgcube = ctools.ctbkgcube(sim.obs()) + bkgcube["incube"] = "NONE" + bkgcube["ebinalg"] = "LOG" + bkgcube["emin"] = emin + bkgcube["emax"] = emax + bkgcube["enumbins"] = enumbins + bkgcube["nxpix"] = 10 + bkgcube["nypix"] = 10 + bkgcube["binsz"] = 1.0 + bkgcube["coordsys"] = coordsys + bkgcube["proj"] = proj + bkgcube["xref"] = ra + bkgcube["yref"] = dec + bkgcube.run() + + # Attach background model to observation container + bin.obs().models(bkgcube.models()) + + # Set Exposure and Psf cube for first CTA observation + # (ctbin will create an observation with a single container) + bin.obs()[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube()) + + # Get ctlike start CPU time + cpu_ctlike = time.clock() + + # Perform maximum likelihood fitting + like = ctools.ctlike(bin.obs()) + like.run() + + # Get stop CPU time and compute elapsed times + cpu_stop = time.clock() + cpu_elapsed = cpu_stop - cpu_start + cpu_ctlike = cpu_stop - cpu_ctlike + + # Return + return cpu_elapsed, cpu_ctlike + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + + # Initialise parameters + title = "CTA analysis benchmark" + model_name = "${CTOOLS}/share/models/crab.xml" + + # Log header + print("*************************************") + print("* CTA analysis benchmark *") + print("*************************************") + print(" Duration Unbinned Binned Stacked ctlike-UB ctlike-B ctlike-S") + + # Initialize arrays + a_duration = [] + a_unbinned = [] + a_binned = [] + a_stacked = [] + c_unbinned = [] + c_binned = [] + c_stacked = [] + + # Perform benchmarks + duration = 1800.0 + for i in range(10): + + # Perform analyses + t_unbinned, ct_unbinned = unbinned_pipeline(model_name, duration) + t_binned, ct_binned = binned_pipeline(model_name, duration) + t_stacked, ct_stacked = stacked_pipeline(model_name, duration) + + # Collect results + a_duration.append(duration/3600.0) + a_unbinned.append(t_unbinned) + a_binned.append(t_binned) + a_stacked.append(t_stacked) + c_unbinned.append(ct_unbinned) + c_binned.append(ct_binned) + c_stacked.append(ct_stacked) + + # Print results + print("%10.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f" % \ + (duration, t_unbinned, t_binned, t_stacked, \ + ct_unbinned, ct_binned, ct_stacked)) + + # Increment duration + duration *= 2.0 + + # Optionally plot results + if has_matplotlib: + plt.figure(1) + plt.title(title) + plt.loglog(a_duration, a_unbinned, 'ro-', label='unbinned') + plt.loglog(a_duration, a_binned, 'bo-', label='binned') + plt.loglog(a_duration, a_stacked, 'go-', label='stacked') + plt.loglog(a_duration, c_unbinned, 'ro--') + plt.loglog(a_duration, c_binned, 'bo--') + plt.loglog(a_duration, c_stacked, 'go--') + plt.xlabel("Observation duration (hours)") + plt.ylabel("CPU time (seconds)") + plt.legend(loc="lower right") + plt.show() diff --git a/test/dev/check_models.py b/test/dev/check_models.py new file mode 100755 index 0000000..2a4902b --- /dev/null +++ b/test/dev/check_models.py @@ -0,0 +1,128 @@ +#! /usr/bin/env python +# ========================================================================== +# Performs unbinned maximum likelihood analysis for various models +# +# Copyright (C) 2011-2016 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 sys +import ctools + + +# ================= # +# Analysis pipeline # +# ================= # +def pipeline(model_name): + """ + Unbinned in-memory analysis pipeline + + 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 + caldb = 'prod2' + irf = 'South_0.5h' + ra = 83.63 + dec = 22.01 + rad_sim = 5.0 + tstart = 0.0 + tstop = 180.0 + emin = 0.1 + emax = 100.0 + rad_select = 3.0 + + # Write model name + print('*** Model: '+model_name+' ************************************') + + # Initialise timing + wall_seconds = 0.0 + cpu_seconds = 0.0 + + # Simulate events + sim = ctools.ctobssim() + sim['inmodel'] = model_name + sim['caldb'] = caldb + sim['irf'] = irf + sim['ra'] = ra + sim['dec'] = dec + sim['rad'] = rad_sim + sim['tmin'] = tstart + sim['tmax'] = tstop + sim['emin'] = emin + sim['emax'] = emax + sim.run() + + # Select events + select = ctools.ctselect(sim.obs()) + select['ra'] = ra + select['dec'] = dec + select['rad'] = rad_select + select['tmin'] = tstart + select['tmax'] = tstop + select['emin'] = emin + select['emax'] = emax + select.run() + + # Perform maximum likelihood fitting + like = ctools.ctlike(select.obs()) + like.run() + + # Show model fitting results + print(like.obs().models()[0]) + + # Return + return + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + + # Initialise flags + need_help = False + + # Test for command line arguments + print(sys.argv[0]) + if (len(sys.argv) > 1): + if sys.argv[1] == '-h': + need_help = True + else: + need_help = True + + # Print help if needed and exit + if need_help: + print('Usage: check_models.py [OPTIONS]') + print(' -h Display this usage message') + sys.exit() + + # Dump header + print('*******************************************') + print('* Check models using an unbinned analysis *') + print('*******************************************') + + # Perform analysis for Crab model + pipeline('$CTOOLS/share/models/crab.xml') + + # Perform analysis for disk model + pipeline('$CTOOLS/share/models/disk.xml') + + # Perform analysis for Gaussian model + pipeline('$CTOOLS/share/models/gauss.xml') + + # Perform analysis for shell model + pipeline('$CTOOLS/share/models/shell.xml') diff --git a/test/dev/simulate_events.py b/test/dev/simulate_events.py new file mode 100755 index 0000000..01fcac6 --- /dev/null +++ b/test/dev/simulate_events.py @@ -0,0 +1,267 @@ +#! /usr/bin/env python +# ========================================================================== +# This script simulates events from a data-space model and compares them to +# the model. +# +# If matplotlib is installed, the event spectrum will be displayed. +# +# -------------------------------------------------------------------------- +# +# Copyright (C) 2014 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools +import math +try: + import matplotlib.pyplot as plt + has_matplotlib = True +except: + has_matplotlib = False + + +# ======================== # +# Create CTA observation # +# ======================== # +def createobs(ra=86.171648, dec=-1.4774586, rad=5.0, + emin=0.1, emax=100.0, duration=360000.0, deadc=0.95, + irf="South_50h", caldb="prod2"): + """ + Create CTA observation. + """ + # Allocate CTA observation + obs = gammalib.GCTAObservation() + + # Set calibration database + db = gammalib.GCaldb() + if (gammalib.dir_exists(caldb)): + db.rootdir(caldb) + else: + db.open("cta", caldb) + + # Set pointing direction + pntdir = gammalib.GSkyDir() + pntdir.radec_deg(ra, dec) + pnt = gammalib.GCTAPointing() + pnt.dir(pntdir) + obs.pointing(pnt) + + # Set ROI + roi = gammalib.GCTARoi() + instdir = gammalib.GCTAInstDir() + instdir.dir(pntdir) + roi.centre(instdir) + roi.radius(rad) + + # Set GTI + gti = gammalib.GGti() + start = gammalib.GTime(0.0) + stop = gammalib.GTime(duration) + gti.append(start, stop) + + # Set energy boundaries + ebounds = gammalib.GEbounds() + e_min = gammalib.GEnergy() + e_max = gammalib.GEnergy() + e_min.TeV(emin) + e_max.TeV(emax) + ebounds.append(e_min, e_max) + + # Allocate event list + events = gammalib.GCTAEventList() + events.roi(roi) + events.gti(gti) + events.ebounds(ebounds) + obs.events(events) + + # Set instrument response + obs.response(irf, db) + + # Set ontime, livetime, and deadtime correction factor + obs.ontime(duration) + obs.livetime(duration*deadc) + obs.deadc(deadc) + + # Return observation + return obs + + +# =============== # +# Simulate events # +# =============== # +def simulate(obs, xmlname): + """ + Simulate events. + """ + # Allocate random number generator + ran = gammalib.GRan() + + # Load models and extract first model + models = gammalib.GModels(xmlname) + model = models[0] + print(model) + + # Simulate events + events = model.mc(obs, ran) + + # Print event statistics + npred = obs.npred(models) + print(str(len(events)) + " events simulated.") + print(str(npred) + " events expected from Npred.") + + # Check distance + pntdir = obs.pointing().dir() + for event in events: + distance = gammalib.GCTAInstDir(event.dir()).dir().dist_deg(pntdir) + if (distance > 5.0): + print("Too far: "+str(distance)) + + # Return events + return events + + +# ============================== # +# Simulate events using ctobssim # +# ============================== # +def simulate_ctobssim(obs, xmlname, seed=0): + """ + Simulate events using ctobssim. + """ + # Create containers + observations = gammalib.GObservations() + observations.append(obs) + + # Append models + models = gammalib.GModels(xmlname) + observations.models(models) + print(models[0]) + + # Allocate ctobssim application and set parameters + sim = ctools.ctobssim(observations) + sim['seed'].integer(seed) + + # Run simulator + sim.run() + + # Retrieve events + events = sim.obs()[0].events().copy() + + # Print event statistics + npred = obs.npred(models) + print(str(len(events)) + " events simulated.") + print(str(npred) + " events expected from Npred.") + + # Delete the simulation + del sim + + # Return events + return events + + +# =========== # +# Show events # +# =========== # +def show_events(events, xmlname, duration, emin, emax, ebins=30): + """ + Show events using matplotlib. + """ + # Create figure + plt.figure(1) + plt.title("MC simulated event spectrum (" + str(emin) + '-' + str(emax) + " TeV)") + + # Setup energy range covered by data + ebds = gammalib.GEbounds(ebins, gammalib.GEnergy(emin, "TeV"), + gammalib.GEnergy(emax, "TeV")) + + # Create energy axis + energy = [] + for i in range(ebds.size()): + energy.append(ebds.elogmean(i).TeV()) + + # Fill histogram + counts = [0.0 for i in range(ebds.size())] + for event in events: + index = ebds.index(event.energy()) + counts[index] = counts[index] + 1.0 + + # Create error bars + error = [math.sqrt(c) for c in counts] + + # Get model values + sum = 0.0 + models = gammalib.GModels(xmlname) + m = models[0] + model = [] + t = gammalib.GTime() + for i in range(ebds.size()): + eval = ebds.elogmean(i) + ewidth = ebds.emax(i) - ebds.emin(i) + f = m.npred(eval, t, obs) * ewidth.MeV() * duration + sum += f + model.append(f) + print(str(sum) + " events expected from spectrum (integration).") + + # Plot data + plt.loglog(energy, counts, 'ro') + plt.errorbar(energy, counts, error, fmt=None, ecolor='r') + + # Plot model + plt.plot(energy, model, 'b-') + + # Set axes + plt.xlabel("Energy (TeV)") + plt.ylabel("Number of events") + + # Notify + print("PLEASE CLOSE WINDOW TO CONTINUE ...") + + # Show plot + plt.show() + + # Return + return + + +#==========================# +# Main routine entry point # +#==========================# +if __name__ == '__main__': + + # Dump header + print('') + print('*******************') + print('* Simulate events *') + print('*******************') + + # Set XML names + xmlname = 'model_bgd.xml' + + # Set simulation parameters + duration = 360000.0 # Observation duration (s) + emin = 0.1 # 0.1 TeV + emax = 100.0 # 100 TeV + ebins = 30 + + # Set observation + obs = createobs(emin=emin, emax=emax, duration=duration) + + # Perform simulation + #events = simulate(obs, xmlname) + events = simulate_ctobssim(obs, xmlname) + + # Show events + if has_matplotlib: + show_events(events, xmlname, duration, emin, emax, ebins=ebins) diff --git a/test/test_examples.py b/test/test_examples.py index 78b174c..b8ceb0e 100755 --- a/test/test_examples.py +++ b/test/test_examples.py @@ -77,6 +77,7 @@ class Test(gammalib.GPythonTestSuite): self.name('examples') # Append tests + self.append(self.test_generate_prod3_irfs, 'Test generate_prod3_irfs') self.append(self.test_pipeline_binned_disk, 'Test pipeline_binned_disk') self.append(self.test_pipeline_binned_mem, 'Test pipeline_binned_mem') self.append(self.test_pipeline_stacked_disk, 'Test pipeline_stacked_disk') @@ -87,6 +88,19 @@ class Test(gammalib.GPythonTestSuite): # Return return + # Test generate_prod3_irfs + def test_generate_prod3_irfs(self): + """ + Test generate_prod3_irfs + """ + # Execute script + self._execute_python('generate_prod3_irfs') + + #TODO: Do any testing + + # Return + return + # Test pipeline_binned_disk def test_pipeline_binned_disk(self): """ -- libgit2 0.21.4