Michael Mayer / ctools

Commit bbe311d49fb48a70f9a3222e0026f966a05d1ac7

Authored by Jürgen Knödlseder
1 parent 80c639dc

Adapt example scripts

Inline Side-by-side
Showing 10 changed files with 178 additions and 133 deletions
examples/Makefile.am
  View file @bbe311d
1 1 ## Process this file with automake to produce Makefile.in
2 2  
3 3 # Scripts to be distributed
4   -scripts = README \
  4 +scripts = README.md \
5 5 benchmark_cta_analysis.py \
6 6 check_models.py \
7 7 make_pull_at_sensitivity_limit.py \
... ...
examples/README.md
  View file @bbe311d
  1 +ctools examples
  2 +===============
  3 +
1 4 This folder contains Python example scripts that illustrate the
2 5 usage of the ctools.
3 6  
4   -benchmark_cta_analysis.py
5   - This script peforms a benchmark of the various CTA analysis types
6   - (unbinned, binned, stacked), and if matplotlib is installed, creates
7   - a plot of the benchmark results.
  7 +`benchmark_cta_analysis.py`
  8 +> This script peforms a benchmark of the various CTA analysis types
  9 +> (unbinned, binned, stacked), and if matplotlib is installed, creates
  10 +> a plot of the benchmark results.
  11 +
  12 +`check_models.py`
  13 +> Checks a variety of models using an unbinned analysis.
  14 +
  15 +`generate_prod3_irfs/py`
  16 +> Generate Prod3 Instrument Response Functions from FITS files.
  17 +
  18 +`make_pointings.py`
  19 +> This script creates pointings for CTA surveys. For the moment the
  20 +> Galactic Plane Survey (gps), the Extragalactic Survey (extgal),
  21 +> the Galactic Centre Survey (gc) and the LMC survey (lmc) are
  22 +> implemented. The script outputs and ASCII file that can be used by
  23 +> csobsdef.
8 24  
9   -check_models.py
10   - Checks a variety of models using an unbinned analysis.
  25 +`make_pull_at_sensitivity_limit.py`
  26 +> This script creates pull distributions as function of energy
  27 +> for a test source that is at the 5 sigma sensitivity limit for
  28 +> 50h of effective observing time. It takes as input file the
  29 +> output file generated using cspull.
  30 +> If the Python processing module is installed the script can make
  31 +> use of multiple cores and runs in parallel.
11 32  
12   -make_pointings.py
13   - This script creates pointings for CTA surveys. For the moment the
14   - Galactic Plane Survey (gps), the Extragalactic Survey (extgal),
15   - the Galactic Centre Survey (gc) and the LMC survey (lmc) are
16   - implemented. The script outputs and ASCII file that can be used by
17   - csobsdef.
  33 +`make_spectrum.py`
  34 +> This script shows how to make a spectrum using obsutils.
18 35  
19   -make_pull_at_sensitivity_limit.py
20   - This script creates pull distributions as function of energy
21   - for a test source that is at the 5 sigma sensitivity limit for
22   - 50h of effective observing time. It takes as input file the
23   - output file generated using cspull.
24   - If the Python processing module is installed the script can make
25   - use of multiple cores and runs in parallel.
  36 +`make_survey.py`
  37 +> This script simulates a CTA galactic plane survey.
26 38  
27   -make_spectrum.py
28   - This script shows how to make a spectrum using obsutils.
  39 +`make_ts_distributions.py`
  40 +> This script creates TS distributions as function of energy for
  41 +> an off observation of 50h effective observing time.
  42 +> If the Python processing module is installed the script can make
  43 +> use of multiple cores and runs in parallel.
29 44  
30   -make_survey.py
31   - This script simulates a CTA galactic plane survey.
  45 +`pipeline_binned_disk.py`
  46 +> Example pipeline for simulating events and performing a binned
  47 +> analysis by saving intermediate results to disk.
32 48  
33   -make_ts_distributions.py
34   - This script creates TS distributions as function of energy for
35   - an off observation of 50h effective observing time.
36   - If the Python processing module is installed the script can make
37   - use of multiple cores and runs in parallel.
  49 +`pipeline_binned_mem.py`
  50 +> Example pipeline for simulating events and performing a binned
  51 +> analysis in memory. No files will be written to disk.
38 52  
39   -benchmark_cta_analysis.py
40   - Performs performance benchmarks of the various CTA analysis types.
  53 +`pipeline_stacked_disk.py`
  54 +> Example pipeline for simulating events and performing a stacked
  55 +> analysis by saving intermediate results to disk.
41 56  
42   -pipeline_binned_disk.py
43   - Example pipeline for simulating events and performing a binned
44   - analysis by saving intermediate results to disk.
  57 +`pipeline_stacked_mem.py`
  58 +> Example pipeline for simulating events and performing a stacked
  59 +> analysis in memory. No files will be written to disk.
45 60  
46   -pipeline_binned_mem.py
47   - Example pipeline for simulating events and performing a binned
48   - analysis in memory. No files will be written to disk.
  61 +`pipeline_unbinned_disk.py`
  62 +> Example pipeline for simulating events and performing an unbinned
  63 +> analysis by saving intermediate results to disk.
49 64  
50   -pipeline_stacked_disk.py
51   - Example pipeline for simulating events and performing a stacked
52   - analysis by saving intermediate results to disk.
  65 +`pipeline_unbinned_mem.py`
  66 +> Example pipeline for simulating events and performing an unbinned
  67 +> analysis in memory. No files will be written to disk.
53 68  
54   -pipeline_stacked_mem.py
55   - Example pipeline for simulating events and performing a stacked
56   - analysis in memory. No files will be written to disk.
  69 +`show_butterfly.py`
  70 +> Displays a butterfly generated using the ctbutterfly tool.
57 71  
58   -pipeline_unbinned_disk.py
59   - Example pipeline for simulating events and performing an unbinned
60   - analysis by saving intermediate results to disk.
  72 +`show_model.py`
  73 +> Displays a spectral model.
61 74  
62   -pipeline_unbinned_mem.py
63   - Example pipeline for simulating events and performing an unbinned
64   - analysis in memory. No files will be written to disk.
  75 +`show_pha.py`
  76 +> Displays a spectrum in PHA format.
65 77  
66   -show_butterfly.py
67   - Displays a butterfly generated using the ctbutterfly tool.
  78 +`show_pull_evolution.py`
  79 +> Displays the evolution of the mean and rms pull as function of number of
  80 +> trials.
68 81  
69   -show_model.py
70   - Displays a spectral model.
  82 +`show_pull_histogram.py`
  83 +> Displays pull histograms generated by the cspull script.
71 84  
72   -show_pha.py
73   - Displays a spectrum in PHA format.
  85 +`show_response.py`
  86 +> Displays one or several CTA response functions.
74 87  
75   -show_pull_histogram.py
76   - Displays pull histograms generated by the cspull script.
  88 +`show_sensitivity.py`
  89 +> Displays sensitivity information generated by the cssens script.
77 90  
78   -show_response.py
79   - Displays one or several CTA response functions.
  91 +`show_spectrum.py`
  92 +> Displays spectrum generated by the csspec script.
80 93  
81   -show_sensitivity.py
82   - Displays sensitivity information generated by the cssens script.
  94 +`show_ts_distribution.py`
  95 +> Displays Test Statistic histograms generated using the cstsdist
  96 +> script.
83 97  
84   -show_ts_distribution.py
85   - Displays Test Statistic histograms generated using the cstsdist
86   - script.
  98 +`simulate_events.py`
  99 +> Simulates events and compares the simulation to the model.
... ...
examples/check_models.py
  View file @bbe311d
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script performs an unbinned maximum likelihood analysis for a
4   -# variety of models. This allows checking the models.
  3 +# Performs unbinned maximum likelihood analysis for various models
5 4 #
6   -# Usage:
7   -# ./check_models.py
8   -#
9   -# ==========================================================================
10   -#
11   -# Copyright (C) 2011-2015 Juergen Knoedlseder
  5 +# Copyright (C) 2011-2016 Juergen Knoedlseder
12 6 #
13 7 # This program is free software: you can redistribute it and/or modify
14 8 # it under the terms of the GNU General Public License as published by
... ... @@ -33,26 +27,26 @@ import ctools
33 27 # ================= #
34 28 def pipeline(model_name):
35 29 """
36   - Unbinned analysis pipeline - keep intermediate results in memory.
  30 + Unbinned in-memory analysis pipeline
37 31  
38 32 This function implements an analysis pipeline that successively calls
39 33 ctobssim, ctselect and ctlike without saving the intermediate results as
40 34 FITS files on disk. All data is only hold in memory.
41 35 """
42 36 # Set script parameters
43   - caldb = "prod2"
44   - irf = "South_50h"
  37 + caldb = 'prod2'
  38 + irf = 'South_0.5h'
45 39 ra = 83.63
46 40 dec = 22.01
47   - rad_sim = 10.0
  41 + rad_sim = 5.0
48 42 tstart = 0.0
49   - tstop = 1800.0
  43 + tstop = 180.0
50 44 emin = 0.1
51 45 emax = 100.0
52 46 rad_select = 3.0
53 47  
54 48 # Write model name
55   - print("*** Model: "+model_name+" ************************************")
  49 + print('*** Model: '+model_name+' ************************************')
56 50  
57 51 # Initialise timing
58 52 wall_seconds = 0.0
... ... @@ -60,27 +54,27 @@ def pipeline(model_name):
60 54  
61 55 # Simulate events
62 56 sim = ctools.ctobssim()
63   - sim["inmodel"] = model_name
64   - sim["caldb"] = caldb
65   - sim["irf"] = irf
66   - sim["ra"] = ra
67   - sim["dec"] = dec
68   - sim["rad"] = rad_sim
69   - sim["tmin"] = tstart
70   - sim["tmax"] = tstop
71   - sim["emin"] = emin
72   - sim["emax"] = emax
  57 + sim['inmodel'] = model_name
  58 + sim['caldb'] = caldb
  59 + sim['irf'] = irf
  60 + sim['ra'] = ra
  61 + sim['dec'] = dec
  62 + sim['rad'] = rad_sim
  63 + sim['tmin'] = tstart
  64 + sim['tmax'] = tstop
  65 + sim['emin'] = emin
  66 + sim['emax'] = emax
73 67 sim.run()
74 68  
75 69 # Select events
76 70 select = ctools.ctselect(sim.obs())
77   - select["ra"] = ra
78   - select["dec"] = dec
79   - select["rad"] = rad_select
80   - select["tmin"] = tstart
81   - select["tmax"] = tstop
82   - select["emin"] = emin
83   - select["emax"] = emax
  71 + select['ra'] = ra
  72 + select['dec'] = dec
  73 + select['rad'] = rad_select
  74 + select['tmin'] = tstart
  75 + select['tmax'] = tstop
  76 + select['emin'] = emin
  77 + select['emax'] = emax
84 78 select.run()
85 79  
86 80 # Perform maximum likelihood fitting
... ... @@ -98,39 +92,37 @@ def pipeline(model_name):
98 92 # Main routine entry point #
99 93 #==========================#
100 94 if __name__ == '__main__':
101   - """
102   - Perform unbinned analyses for various models.
103   - """
  95 +
104 96 # Initialise flags
105 97 need_help = False
106 98  
107 99 # Test for command line arguments
108 100 print(sys.argv[0])
109 101 if (len(sys.argv) > 1):
110   - if sys.argv[1] == "-h":
  102 + if sys.argv[1] == '-h':
111 103 need_help = True
112 104 else:
113 105 need_help = True
114 106  
115 107 # Print help if needed and exit
116 108 if need_help:
117   - print("Usage: check_models.py [OPTIONS]")
118   - print(" -h Display this usage message")
  109 + print('Usage: check_models.py [OPTIONS]')
  110 + print(' -h Display this usage message')
119 111 sys.exit()
120 112  
121 113 # Dump header
122   - print("*******************************************")
123   - print("* Check models using an unbinned analysis *")
124   - print("*******************************************")
  114 + print('*******************************************')
  115 + print('* Check models using an unbinned analysis *')
  116 + print('*******************************************')
125 117  
126 118 # Perform analysis for Crab model
127   - pipeline("$CTOOLS/share/models/crab.xml")
  119 + pipeline('$CTOOLS/share/models/crab.xml')
128 120  
129 121 # Perform analysis for disk model
130   - pipeline("$CTOOLS/share/models/disk.xml")
  122 + pipeline('$CTOOLS/share/models/disk.xml')
131 123  
132 124 # Perform analysis for Gaussian model
133   - pipeline("$CTOOLS/share/models/gauss.xml")
  125 + pipeline('$CTOOLS/share/models/gauss.xml')
134 126  
135 127 # Perform analysis for shell model
136   - pipeline("$CTOOLS/share/models/shell.xml")
  128 + pipeline('$CTOOLS/share/models/shell.xml')
... ...
examples/pipeline_binned_disk.py
  View file @bbe311d
... ... @@ -179,8 +179,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0,
179 179 # Main routine entry point #
180 180 # ======================== #
181 181 if __name__ == '__main__':
182   - """
183   - """
  182 +
184 183 # Dump header
185 184 print("********************************************")
186 185 print("* CTA binned analysis pipeline *")
... ...
examples/pipeline_binned_mem.py
  View file @bbe311d
... ... @@ -153,8 +153,7 @@ def run_pipeline(obs, emin=0.1, emax=100.0,
153 153 # Main routine entry point #
154 154 # ======================== #
155 155 if __name__ == '__main__':
156   - """
157   - """
  156 +
158 157 # Dump header
159 158 print("********************************************")
160 159 print("* CTA binned analysis pipeline *")
... ...
examples/pipeline_stacked_disk.py
  View file @bbe311d
... ... @@ -224,8 +224,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0,
224 224 # Main routine entry point #
225 225 # ======================== #
226 226 if __name__ == '__main__':
227   - """
228   - """
  227 +
229 228 # Dump header
230 229 print("********************************************")
231 230 print("* CTA stacked analysis pipeline *")
... ...
examples/pipeline_stacked_mem.py
  View file @bbe311d
... ... @@ -203,8 +203,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0,
203 203 # Main routine entry point #
204 204 # ======================== #
205 205 if __name__ == '__main__':
206   - """
207   - """
  206 +
208 207 # Dump header
209 208 print("********************************************")
210 209 print("* CTA stacked analysis pipeline *")
... ...
examples/pipeline_unbinned_disk.py
  View file @bbe311d
... ... @@ -154,8 +154,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0,
154 154 # Main routine entry point #
155 155 # ======================== #
156 156 if __name__ == '__main__':
157   - """
158   - """
  157 +
159 158 # Dump header
160 159 print("********************************************")
161 160 print("* CTA unbinned analysis pipeline *")
... ...
examples/pipeline_unbinned_mem.py
  View file @bbe311d
... ... @@ -139,8 +139,7 @@ def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0,
139 139 # Main routine entry point #
140 140 # ======================== #
141 141 if __name__ == '__main__':
142   - """
143   - """
  142 +
144 143 # Dump header
145 144 print("********************************************")
146 145 print("* CTA unbinned analysis pipeline *")
... ...
src/support/ctool.cpp
  View file @bbe311d
... ... @@ -52,6 +52,7 @@
52 52 #define G_SETUP_OBSERVATION "ctool::setup_observations(GObservations&)"
53 53 #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)"
54 54 #define G_CREATE_EBOUNDS "ctool::create_ebounds()"
  55 +#define G_SET_OBS_RESPONSE "ctool::set_obs_response(GCTAObservation*)"
55 56 #define G_PROVIDE_HELP "ctool::provide_help()"
56 57  
57 58 /* __ Debug definitions __________________________________________________ */
... ... @@ -931,6 +932,9 @@ void ctool::set_response(GObservations& obs)
931 932 *
932 933 * @param[in,out] obs CTA observation
933 934 *
  935 + * @exception GException::invalid_value
  936 + * Energy dispersion requested but not energy dispersion cube given
  937 + *
934 938 * Set the response for one CTA observation. If the CTA observation contains
935 939 * a counts cube the method first attempts to set the response information
936 940 * using the following user parameters
... ... @@ -964,13 +968,23 @@ void ctool::set_obs_response(GCTAObservation* obs)
964 968 // response information for a stacked analysis is provided, and if so,
965 969 // set that response information for the CTA observation
966 970 if (dynamic_cast<const GCTAEventCube*>(obs->events()) != NULL) {
967   - if (has_par("expcube") && has_par("psfcube") &&
968   - has_par("edispcube") && has_par("bkgcube")) {
  971 +
  972 + // Check if response cube parameters are provided. We need all four
  973 + // response cube parameters.
  974 + if (has_par("expcube") && has_par("psfcube") &&
  975 + has_par("bkgcube") && has_par("edispcube")) {
  976 +
  977 + // Check if we need to query the energy dispersion cube. We do
  978 + // not need to query the cube is an "edisp" parameter exists and
  979 + // if this parameter is set to "no".
  980 + bool query_edisp = true;
  981 + if (has_par("edisp")) {
  982 + query_edisp = (*this)["edisp"].boolean();
  983 + }
969 984  
970 985 // Get filenames
971 986 GFilename expcube((*this)["expcube"].filename());
972 987 GFilename psfcube((*this)["psfcube"].filename());
973   - GFilename edispcube((*this)["edispcube"].filename());
974 988 GFilename bkgcube((*this)["bkgcube"].filename());
975 989  
976 990 // Extract stacked response information if available
... ... @@ -981,21 +995,53 @@ void ctool::set_obs_response(GCTAObservation* obs)
981 995 !psfcube.is_empty() &&
982 996 !bkgcube.is_empty()) {
983 997  
984   - // Get exposure, PSF and background cubes
  998 + // Load exposure, PSF and background cubes
985 999 GCTACubeExposure exposure(expcube);
986 1000 GCTACubePsf psf(psfcube);
987 1001 GCTACubeBackground background(bkgcube);
988 1002  
989   - // If an energy dispersion cube is available, get also
990   - // this cube, and set all information as stacked response
991   - // for this CTA observation
992   - if ((gammalib::toupper(edispcube.url()) != "NONE") &&
993   - !edispcube.is_empty()) {
994   - GCTACubeEdisp edisp(edispcube);
995   - obs->response(exposure, psf, edisp, background);
996   - }
  1003 + // If querying of energy dispersion cube is requested then
  1004 + // query it now
  1005 + if (query_edisp) {
  1006 +
  1007 + // Get filename
  1008 + GFilename edispcube((*this)["edispcube"].filename());
  1009 +
  1010 + // If filename is valid then use energy dispersion ...
  1011 + if ((gammalib::toupper(edispcube.url()) != "NONE") &&
  1012 + !edispcube.is_empty()) {
  1013 +
  1014 + // Load energy dispersion cube
  1015 + GCTACubeEdisp edisp(edispcube);
  1016 +
  1017 + // Set response with all four cubes
  1018 + obs->response(exposure, psf, edisp, background);
  1019 +
  1020 + } // endif: energy dispersion cube was provided
  1021 +
  1022 + // ... otherwise work without energy dispersion
  1023 + else {
  1024 +
  1025 + // If energy dispersion was requested but no cube
  1026 + // was provided then throw an exception
  1027 + if (has_par("edisp") && (*this)["edisp"].boolean()) {
  1028 + std::string msg = "Energy dispersion requested but "
  1029 + "no energy dispersion cube was "
  1030 + "specified.";
  1031 + throw GException::invalid_value(G_SET_OBS_RESPONSE,
  1032 + msg);
  1033 + }
  1034 +
  1035 + // Set response without energy dispersion
  1036 + obs->response(exposure, psf, background);
  1037 +
  1038 + } // endelse: no energy dispersion was available
  1039 +
  1040 + } // endif: energy dispersion needed
  1041 +
  1042 + // ... otherwise work without energy dispersion
997 1043 else {
998   - obs->response(exposure, psf, background);
  1044 + obs->response(exposure, psf, background);
999 1045 }
1000 1046  
1001 1047 // Signal that response is available
... ...