From bbe311d49fb48a70f9a3222e0026f966a05d1ac7 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 24 Jun 2016 16:52:37 +0200 Subject: [PATCH] Adapt example scripts --- examples/Makefile.am | 2 +- examples/README.md | 139 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------- examples/check_models.py | 80 ++++++++++++++++++++++++++++++++++++-------------------------------------------- examples/pipeline_binned_disk.py | 3 +-- examples/pipeline_binned_mem.py | 3 +-- examples/pipeline_stacked_disk.py | 3 +-- examples/pipeline_stacked_mem.py | 3 +-- examples/pipeline_unbinned_disk.py | 3 +-- examples/pipeline_unbinned_mem.py | 3 +-- src/support/ctool.cpp | 72 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------- 10 files changed, 178 insertions(+), 133 deletions(-) diff --git a/examples/Makefile.am b/examples/Makefile.am index e867429..c37ee57 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,7 +1,7 @@ ## Process this file with automake to produce Makefile.in # Scripts to be distributed -scripts = README \ +scripts = README.md \ benchmark_cta_analysis.py \ check_models.py \ make_pull_at_sensitivity_limit.py \ diff --git a/examples/README.md b/examples/README.md index a7fc211..291918c 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1,86 +1,99 @@ +ctools examples +=============== + 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. +`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. + +`check_models.py` +> Checks a variety of models using an unbinned analysis. + +`generate_prod3_irfs/py` +> Generate Prod3 Instrument Response Functions from FITS files. + +`make_pointings.py` +> This script creates pointings for CTA surveys. For the moment the +> Galactic Plane Survey (gps), the Extragalactic Survey (extgal), +> the Galactic Centre Survey (gc) and the LMC survey (lmc) are +> implemented. The script outputs and ASCII file that can be used by +> csobsdef. -check_models.py - Checks a variety of models using an unbinned analysis. +`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_pointings.py - This script creates pointings for CTA surveys. For the moment the - Galactic Plane Survey (gps), the Extragalactic Survey (extgal), - the Galactic Centre Survey (gc) and the LMC survey (lmc) are - implemented. The script outputs and ASCII file that can be used by - csobsdef. +`make_spectrum.py` +> This script shows how to make a spectrum using obsutils. -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_survey.py` +> This script simulates a CTA galactic plane survey. -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. -make_survey.py - This script simulates a CTA galactic plane survey. +`pipeline_binned_disk.py` +> Example pipeline for simulating events and performing a binned +> analysis by saving intermediate results to disk. -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. +`pipeline_binned_mem.py` +> Example pipeline for simulating events and performing a binned +> analysis in memory. No files will be written to disk. -benchmark_cta_analysis.py - Performs performance benchmarks of the various CTA analysis types. +`pipeline_stacked_disk.py` +> Example pipeline for simulating events and performing a stacked +> analysis by saving intermediate results to disk. -pipeline_binned_disk.py - Example pipeline for simulating events and performing a binned - analysis by saving intermediate results to disk. +`pipeline_stacked_mem.py` +> Example pipeline for simulating events and performing a stacked +> analysis in memory. No files will be written to disk. -pipeline_binned_mem.py - Example pipeline for simulating events and performing a binned - analysis in memory. No files will be written to disk. +`pipeline_unbinned_disk.py` +> Example pipeline for simulating events and performing an unbinned +> analysis by saving intermediate results to disk. -pipeline_stacked_disk.py - Example pipeline for simulating events and performing a stacked - analysis by saving intermediate results to disk. +`pipeline_unbinned_mem.py` +> Example pipeline for simulating events and performing an unbinned +> analysis in memory. No files will be written to disk. -pipeline_stacked_mem.py - Example pipeline for simulating events and performing a stacked - analysis in memory. No files will be written to disk. +`show_butterfly.py` +> Displays a butterfly generated using the ctbutterfly tool. -pipeline_unbinned_disk.py - Example pipeline for simulating events and performing an unbinned - analysis by saving intermediate results to disk. +`show_model.py` +> Displays a spectral model. -pipeline_unbinned_mem.py - Example pipeline for simulating events and performing an unbinned - analysis in memory. No files will be written to disk. +`show_pha.py` +> Displays a spectrum in PHA format. -show_butterfly.py - Displays a butterfly generated using the ctbutterfly tool. +`show_pull_evolution.py` +> Displays the evolution of the mean and rms pull as function of number of +> trials. -show_model.py - Displays a spectral model. +`show_pull_histogram.py` +> Displays pull histograms generated by the cspull script. -show_pha.py - Displays a spectrum in PHA format. +`show_response.py` +> Displays one or several CTA response functions. -show_pull_histogram.py - Displays pull histograms generated by the cspull script. +`show_sensitivity.py` +> Displays sensitivity information generated by the cssens script. -show_response.py - Displays one or several CTA response functions. +`show_spectrum.py` +> Displays spectrum generated by the csspec script. -show_sensitivity.py - Displays sensitivity information generated by the cssens script. +`show_ts_distribution.py` +> Displays Test Statistic histograms generated using the cstsdist +> script. -show_ts_distribution.py - Displays Test Statistic histograms generated using the cstsdist - script. +`simulate_events.py` +> Simulates events and compares the simulation to the model. diff --git a/examples/check_models.py b/examples/check_models.py index 936d61d..2a4902b 100755 --- a/examples/check_models.py +++ b/examples/check_models.py @@ -1,14 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This script performs an unbinned maximum likelihood analysis for a -# variety of models. This allows checking the models. +# Performs unbinned maximum likelihood analysis for various models # -# Usage: -# ./check_models.py -# -# ========================================================================== -# -# Copyright (C) 2011-2015 Juergen Knoedlseder +# 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 @@ -33,26 +27,26 @@ import ctools # ================= # def pipeline(model_name): """ - Unbinned analysis pipeline - keep intermediate results in memory. + 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_50h" + caldb = 'prod2' + irf = 'South_0.5h' ra = 83.63 dec = 22.01 - rad_sim = 10.0 + rad_sim = 5.0 tstart = 0.0 - tstop = 1800.0 + tstop = 180.0 emin = 0.1 emax = 100.0 rad_select = 3.0 # Write model name - print("*** Model: "+model_name+" ************************************") + print('*** Model: '+model_name+' ************************************') # Initialise timing wall_seconds = 0.0 @@ -60,27 +54,27 @@ def pipeline(model_name): # 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['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['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 @@ -98,39 +92,37 @@ def pipeline(model_name): # Main routine entry point # #==========================# if __name__ == '__main__': - """ - Perform unbinned analyses for various models. - """ + # Initialise flags need_help = False # Test for command line arguments print(sys.argv[0]) if (len(sys.argv) > 1): - if sys.argv[1] == "-h": + 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") + 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("*******************************************") + print('*******************************************') + print('* Check models using an unbinned analysis *') + print('*******************************************') # Perform analysis for Crab model - pipeline("$CTOOLS/share/models/crab.xml") + pipeline('$CTOOLS/share/models/crab.xml') # Perform analysis for disk model - pipeline("$CTOOLS/share/models/disk.xml") + pipeline('$CTOOLS/share/models/disk.xml') # Perform analysis for Gaussian model - pipeline("$CTOOLS/share/models/gauss.xml") + pipeline('$CTOOLS/share/models/gauss.xml') # Perform analysis for shell model - pipeline("$CTOOLS/share/models/shell.xml") + pipeline('$CTOOLS/share/models/shell.xml') diff --git a/examples/pipeline_binned_disk.py b/examples/pipeline_binned_disk.py index 89c3e4d..38b42f7 100755 --- a/examples/pipeline_binned_disk.py +++ b/examples/pipeline_binned_disk.py @@ -179,8 +179,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header print("********************************************") print("* CTA binned analysis pipeline *") diff --git a/examples/pipeline_binned_mem.py b/examples/pipeline_binned_mem.py index f85a5b0..88724da 100755 --- a/examples/pipeline_binned_mem.py +++ b/examples/pipeline_binned_mem.py @@ -153,8 +153,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header print("********************************************") print("* CTA binned analysis pipeline *") diff --git a/examples/pipeline_stacked_disk.py b/examples/pipeline_stacked_disk.py index 3277caf..1c35626 100755 --- a/examples/pipeline_stacked_disk.py +++ b/examples/pipeline_stacked_disk.py @@ -224,8 +224,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header print("********************************************") print("* CTA stacked analysis pipeline *") diff --git a/examples/pipeline_stacked_mem.py b/examples/pipeline_stacked_mem.py index 555055a..36ab850 100755 --- a/examples/pipeline_stacked_mem.py +++ b/examples/pipeline_stacked_mem.py @@ -203,8 +203,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header print("********************************************") print("* CTA stacked analysis pipeline *") diff --git a/examples/pipeline_unbinned_disk.py b/examples/pipeline_unbinned_disk.py index 55fae66..44c7ea6 100755 --- a/examples/pipeline_unbinned_disk.py +++ b/examples/pipeline_unbinned_disk.py @@ -154,8 +154,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header print("********************************************") print("* CTA unbinned analysis pipeline *") diff --git a/examples/pipeline_unbinned_mem.py b/examples/pipeline_unbinned_mem.py index 392c21c..b2deb08 100755 --- a/examples/pipeline_unbinned_mem.py +++ b/examples/pipeline_unbinned_mem.py @@ -139,8 +139,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - """ + # Dump header print("********************************************") print("* CTA unbinned analysis pipeline *") diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index ec19f93..2c78676 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -52,6 +52,7 @@ #define G_SETUP_OBSERVATION "ctool::setup_observations(GObservations&)" #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)" #define G_CREATE_EBOUNDS "ctool::create_ebounds()" +#define G_SET_OBS_RESPONSE "ctool::set_obs_response(GCTAObservation*)" #define G_PROVIDE_HELP "ctool::provide_help()" /* __ Debug definitions __________________________________________________ */ @@ -931,6 +932,9 @@ void ctool::set_response(GObservations& obs) * * @param[in,out] obs CTA observation * + * @exception GException::invalid_value + * Energy dispersion requested but not energy dispersion cube given + * * Set the response for one CTA observation. If the CTA observation contains * a counts cube the method first attempts to set the response information * using the following user parameters @@ -964,13 +968,23 @@ void ctool::set_obs_response(GCTAObservation* obs) // response information for a stacked analysis is provided, and if so, // set that response information for the CTA observation if (dynamic_cast(obs->events()) != NULL) { - if (has_par("expcube") && has_par("psfcube") && - has_par("edispcube") && has_par("bkgcube")) { + + // Check if response cube parameters are provided. We need all four + // response cube parameters. + if (has_par("expcube") && has_par("psfcube") && + has_par("bkgcube") && has_par("edispcube")) { + + // Check if we need to query the energy dispersion cube. We do + // not need to query the cube is an "edisp" parameter exists and + // if this parameter is set to "no". + bool query_edisp = true; + if (has_par("edisp")) { + query_edisp = (*this)["edisp"].boolean(); + } // Get filenames GFilename expcube((*this)["expcube"].filename()); GFilename psfcube((*this)["psfcube"].filename()); - GFilename edispcube((*this)["edispcube"].filename()); GFilename bkgcube((*this)["bkgcube"].filename()); // Extract stacked response information if available @@ -981,21 +995,53 @@ void ctool::set_obs_response(GCTAObservation* obs) !psfcube.is_empty() && !bkgcube.is_empty()) { - // Get exposure, PSF and background cubes + // Load exposure, PSF and background cubes GCTACubeExposure exposure(expcube); GCTACubePsf psf(psfcube); GCTACubeBackground background(bkgcube); - // If an energy dispersion cube is available, get also - // this cube, and set all information as stacked response - // for this CTA observation - if ((gammalib::toupper(edispcube.url()) != "NONE") && - !edispcube.is_empty()) { - GCTACubeEdisp edisp(edispcube); - obs->response(exposure, psf, edisp, background); - } + // If querying of energy dispersion cube is requested then + // query it now + if (query_edisp) { + + // Get filename + GFilename edispcube((*this)["edispcube"].filename()); + + // If filename is valid then use energy dispersion ... + if ((gammalib::toupper(edispcube.url()) != "NONE") && + !edispcube.is_empty()) { + + // Load energy dispersion cube + GCTACubeEdisp edisp(edispcube); + + // Set response with all four cubes + obs->response(exposure, psf, edisp, background); + + } // endif: energy dispersion cube was provided + + // ... otherwise work without energy dispersion + else { + + // If energy dispersion was requested but no cube + // was provided then throw an exception + if (has_par("edisp") && (*this)["edisp"].boolean()) { + std::string msg = "Energy dispersion requested but " + "no energy dispersion cube was " + "specified."; + throw GException::invalid_value(G_SET_OBS_RESPONSE, + msg); + } + + // Set response without energy dispersion + obs->response(exposure, psf, background); + + } // endelse: no energy dispersion was available + + } // endif: energy dispersion needed + + // ... otherwise work without energy dispersion else { - obs->response(exposure, psf, background); + obs->response(exposure, psf, background); } // Signal that response is available -- libgit2 0.21.4