Konstancja Satalecka / ctools

Commit 210ceab145af8c6598a2ca84802095a9d9d4a8b1

Authored by Jürgen Knödlseder
1 parent 0de79511

Add in-memory analysis pipelines

Inline Side-by-side
Showing 4 changed files with 358 additions and 17 deletions
examples/README
  View file @210ceab
... ... @@ -4,10 +4,21 @@ usage of the ctools.
4 4 benchmark_cta_analysis.py
5 5 Performs performance benchmarks of the various CTA analysis types.
6 6  
  7 +pipeline_binned_mem.py
  8 + Example pipeline for simulating events and performing a binned
  9 + analysis in memory. No files will be written to disk.
  10 +
7 11 pipeline_stacked_mem.py
8 12 Example pipeline for simulating events and performing a stacked
9 13 analysis in memory. No files will be written to disk.
10 14  
  15 +pipeline_unbinned_mem.py
  16 + Example pipeline for simulating events and performing an unbinned
  17 + analysis in memory. No files will be written to disk.
  18 +
  19 +
  20 +
  21 +
11 22 make_binned_analysis.py
12 23 Example Python pipelines for binned analysis.
13 24 The first pipeline writes intermediate results to the disk while
... ...
examples/pipeline_binned_mem.py 0 → 100755
  View file @210ceab
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script illustrates how to perform a binned CTA analysis based on
  4 +# simulated CTA data. You may use and adapt this script to implement your
  5 +# own pipeline.
  6 +#
  7 +# Usage:
  8 +# ./pipeline_binned_mem.py
  9 +#
  10 +# ==========================================================================
  11 +#
  12 +# Copyright (C) 2015 Juergen Knoedlseder
  13 +#
  14 +# This program is free software: you can redistribute it and/or modify
  15 +# it under the terms of the GNU General Public License as published by
  16 +# the Free Software Foundation, either version 3 of the License, or
  17 +# (at your option) any later version.
  18 +#
  19 +# This program is distributed in the hope that it will be useful,
  20 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  22 +# GNU General Public License for more details.
  23 +#
  24 +# You should have received a copy of the GNU General Public License
  25 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  26 +#
  27 +# ==========================================================================
  28 +import gammalib
  29 +import ctools
  30 +from ctools import obsutils
  31 +
  32 +
  33 +# ================== #
  34 +# Setup observations #
  35 +# ================== #
  36 +def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, \
  37 + emin=0.1, emax=100.0, rad=5.0, duration=1800.0, \
  38 + deadc=0.95, \
  39 + caldb="dummy", irf="cta_dummy_irf"):
  40 + """
  41 + Returns an observation container.
  42 +
  43 + Keywords:
  44 + pattern - Pointing pattern, either "single" or "four"
  45 + ra - RA of pattern centre [deg] (default: 83.6331)
  46 + dec - DEC of pattern centre [deg] (default: 22.0145)
  47 + offset - Offset between observations of pattern [deg] (default: 1.5)
  48 + emin - Minimum energy [TeV] (default: 0.1)
  49 + emax - Maximum energy [TeV] (default: 100.0)
  50 + rad - ROI radius used for analysis [deg] (default: 5.0)
  51 + duration - Duration of one CTA observation [seconds] (default: 1800.0)
  52 + deadc - Deadtime correction factor (default: 0.95)
  53 + caldb - Calibration database path (default: "dummy")
  54 + irf - Instrument response function (default: cta_dummy_irf)
  55 + """
  56 + # Set list of observations
  57 + obs_def_list = obsutils.set_obs_patterns(pattern, \
  58 + ra=ra, \
  59 + dec=dec, \
  60 + offset=offset)
  61 +
  62 + # Get observation container
  63 + obs = obsutils.set_obs_list(obs_def_list, \
  64 + duration=duration, \
  65 + emin=emin, \
  66 + emax=emax, \
  67 + rad=rad, \
  68 + caldb=caldb, \
  69 + irf=irf)
  70 +
  71 + # Return observation container
  72 + return obs
  73 +
  74 +
  75 +# =========== #
  76 +# Setup model #
  77 +# =========== #
  78 +def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
  79 + """
  80 + Setup model for analysis.
  81 +
  82 + Keywords:
  83 + model - Model Xml file
  84 + """
  85 + # Append model from file to observation container
  86 + obs.models(gammalib.GModels(model))
  87 +
  88 + # Return observation container
  89 + return obs
  90 +
  91 +
  92 +# ================================ #
  93 +# Simulation and analysis pipeline #
  94 +# ================================ #
  95 +def run_pipeline(obs, emin=0.1, emax=100.0, \
  96 + enumbins=20, nxpix=200, nypix=200, binsz=0.02, \
  97 + coordsys="CEL", proj="CAR", debug=False):
  98 + """
  99 + Simulation and binned analysis pipeline.
  100 +
  101 + Keywords:
  102 + emin - Minimum energy of cube [TeV] (default: 0.1)
  103 + emax - Maximum energy of cube [TeV] (default: 100.0)
  104 + enumbins - Number of energy bins in cube (default: 20)
  105 + nxpix - Number of RA pixels in cube (default: 200)
  106 + nypix - Number of DEC pixels in cube (default: 200)
  107 + binsz - Spatial cube bin size [deg] (default: 0.02)
  108 + coordsys - Cube coordinate system (CEL or GAL)
  109 + proj - Cube World Coordinate System (WCS) projection
  110 + debug - Enable debugging (default: False)
  111 + """
  112 + # Simulate events
  113 + sim = ctools.ctobssim(obs)
  114 + sim["debug"].boolean(debug)
  115 + sim.run()
  116 +
  117 + # Bin events by looping over all observations in the container
  118 + obs = gammalib.GObservations()
  119 + obs.models(sim.obs().models())
  120 + for run in sim.obs():
  121 +
  122 + # Create container with a single observation
  123 + container = gammalib.GObservations()
  124 + container.append(run)
  125 +
  126 + # Bin events for that observation
  127 + bin = ctools.ctbin(container)
  128 + bin["ebinalg"].string("LOG")
  129 + bin["emin"].real(emin)
  130 + bin["emax"].real(emax)
  131 + bin["enumbins"].integer(enumbins)
  132 + bin["nxpix"].integer(nxpix)
  133 + bin["nypix"].integer(nypix)
  134 + bin["binsz"].real(binsz)
  135 + bin["coordsys"].string(coordsys)
  136 + bin["usepnt"].boolean(True)
  137 + bin["proj"].string(proj)
  138 + bin.run()
  139 +
  140 + # Append result to observations
  141 + obs.extend(bin.obs())
  142 +
  143 + # Perform maximum likelihood fitting
  144 + like = ctools.ctlike(obs)
  145 + like["debug"].boolean(True) # Switch this always on for results in console
  146 + like.run()
  147 +
  148 + # Return
  149 + return
  150 +
  151 +
  152 +# ======================== #
  153 +# Main routine entry point #
  154 +# ======================== #
  155 +if __name__ == '__main__':
  156 + """
  157 + """
  158 + # Dump header
  159 + print("********************************************")
  160 + print("* CTA binned analysis pipeline *")
  161 + print("********************************************")
  162 +
  163 + # Setup observations
  164 + obs = setup_observations()
  165 +
  166 + # Setup model
  167 + obs = setup_model(obs)
  168 +
  169 + # Run analysis pipeline
  170 + run_pipeline(obs)
... ...
examples/pipeline_stacked_mem.py
  View file @210ceab
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script illustrates how to perform a stacked CTA analysis from the
4   -# simulation to the analysis. You may use and adapt this script to
5   -# implement your own pipeline. You can run the example by typing
  3 +# This script illustrates how to perform a stacked CTA analysis based on
  4 +# simulated CTA data. You may use and adapt this script to implement your
  5 +# own pipeline.
6 6 #
  7 +# Usage:
7 8 # ./pipeline_stacked_mem.py
8 9 #
9   -# Copyright (C) 2014 Juergen Knoedlseder
  10 +# ==========================================================================
  11 +#
  12 +# Copyright (C) 2014-2015 Juergen Knoedlseder
10 13 #
11 14 # This program is free software: you can redistribute it and/or modify
12 15 # it under the terms of the GNU General Public License as published by
... ... @@ -30,7 +33,7 @@ from ctools import obsutils
30 33 # ================== #
31 34 # Setup observations #
32 35 # ================== #
33   -def setup_observations(pattern="four", ra=83.6331, dec=22.0145, offset=1.5, \
  36 +def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, \
34 37 emin=0.1, emax=100.0, rad=5.0, duration=1800.0, \
35 38 deadc=0.95, \
36 39 caldb="dummy", irf="cta_dummy_irf"):
... ... @@ -77,7 +80,7 @@ def setup_model(obs, model=&quot;${CTOOLS}/share/models/crab.xml&quot;):
77 80 Setup model for analysis.
78 81  
79 82 Keywords:
80   - model - Model Xml file
  83 + model - Model Xml file
81 84 """
82 85 # Append model from file to observation container
83 86 obs.models(gammalib.GModels(model))
... ... @@ -89,23 +92,24 @@ def setup_model(obs, model=&quot;${CTOOLS}/share/models/crab.xml&quot;):
89 92 # ================================ #
90 93 # Simulation and analysis pipeline #
91 94 # ================================ #
92   -def run_pipeline(obs, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0, \
  95 +def run_pipeline(obs, ra=83.63, dec=22.01, emin=0.1, emax=100.0, \
93 96 enumbins=20, nxpix=200, nypix=200, binsz=0.02, \
94 97 coordsys="CEL", proj="CAR", debug=False):
95 98 """
96 99 Simulation and stacked analysis pipeline.
97 100  
98 101 Keywords:
99   - ra - RA of cube centre [deg] (default: 83.6331)
100   - dec - DEC of cube centre [deg] (default: 22.0145)
101   - emin - Minimum energy of cube [TeV] (default: 0.1)
102   - emax - Maximum energy of cube [TeV] (default: 100.0)
103   - enumbins - Number of energy bins in cube (default: 20)
104   - nxpix - Number of RA pixels in cube (default: 200)
105   - nypix - Number of DEC pixels in cube (default: 200)
106   - binsz - Spatial cube bin size [deg] (default: 0.02)
107   - coordsys - Cube coordinate system (CEL or GAL)
108   - proj - Cube World Coordinate System (WCS) projection
  102 + ra - RA of cube centre [deg] (default: 83.6331)
  103 + dec - DEC of cube centre [deg] (default: 22.0145)
  104 + emin - Minimum energy of cube [TeV] (default: 0.1)
  105 + emax - Maximum energy of cube [TeV] (default: 100.0)
  106 + enumbins - Number of energy bins in cube (default: 20)
  107 + nxpix - Number of RA pixels in cube (default: 200)
  108 + nypix - Number of DEC pixels in cube (default: 200)
  109 + binsz - Spatial cube bin size [deg] (default: 0.02)
  110 + coordsys - Cube coordinate system (CEL or GAL)
  111 + proj - Cube World Coordinate System (WCS) projection
  112 + debug - Enable debugging (default: False)
109 113 """
110 114 # Simulate events
111 115 sim = ctools.ctobssim(obs)
... ...
examples/pipeline_unbinned_mem.py 0 → 100755
  View file @210ceab
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script illustrates how to perform an unbinned CTA analysis based on
  4 +# simulated CTA data. You may use and adapt this script to implement your
  5 +# own pipeline.
  6 +#
  7 +# Usage:
  8 +# ./pipeline_unbinned_mem.py
  9 +#
  10 +# ==========================================================================
  11 +#
  12 +# Copyright (C) 2015 Juergen Knoedlseder
  13 +#
  14 +# This program is free software: you can redistribute it and/or modify
  15 +# it under the terms of the GNU General Public License as published by
  16 +# the Free Software Foundation, either version 3 of the License, or
  17 +# (at your option) any later version.
  18 +#
  19 +# This program is distributed in the hope that it will be useful,
  20 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  22 +# GNU General Public License for more details.
  23 +#
  24 +# You should have received a copy of the GNU General Public License
  25 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  26 +#
  27 +# ==========================================================================
  28 +import gammalib
  29 +import ctools
  30 +from ctools import obsutils
  31 +
  32 +
  33 +# ================== #
  34 +# Setup observations #
  35 +# ================== #
  36 +def setup_observations(pattern="four", ra=83.63, dec=22.01, offset=1.5, \
  37 + emin=0.1, emax=100.0, rad=5.0, duration=1800.0, \
  38 + deadc=0.95, \
  39 + caldb="dummy", irf="cta_dummy_irf"):
  40 + """
  41 + Returns an observation container.
  42 +
  43 + Keywords:
  44 + pattern - Pointing pattern, either "single" or "four"
  45 + ra - RA of pattern centre [deg] (default: 83.6331)
  46 + dec - DEC of pattern centre [deg] (default: 22.0145)
  47 + offset - Offset between observations of pattern [deg] (default: 1.5)
  48 + emin - Minimum energy [TeV] (default: 0.1)
  49 + emax - Maximum energy [TeV] (default: 100.0)
  50 + rad - ROI radius used for analysis [deg] (default: 5.0)
  51 + duration - Duration of one CTA observation [seconds] (default: 1800.0)
  52 + deadc - Deadtime correction factor (default: 0.95)
  53 + caldb - Calibration database path (default: "dummy")
  54 + irf - Instrument response function (default: cta_dummy_irf)
  55 + """
  56 + # Set list of observations
  57 + obs_def_list = obsutils.set_obs_patterns(pattern, \
  58 + ra=ra, \
  59 + dec=dec, \
  60 + offset=offset)
  61 +
  62 + # Get observation container
  63 + obs = obsutils.set_obs_list(obs_def_list, \
  64 + duration=duration, \
  65 + emin=emin, \
  66 + emax=emax, \
  67 + rad=rad, \
  68 + caldb=caldb, \
  69 + irf=irf)
  70 +
  71 + # Return observation container
  72 + return obs
  73 +
  74 +
  75 +# =========== #
  76 +# Setup model #
  77 +# =========== #
  78 +def setup_model(obs, model="${CTOOLS}/share/models/crab.xml"):
  79 + """
  80 + Setup model for analysis.
  81 +
  82 + Keywords:
  83 + model - Model Xml file
  84 + """
  85 + # Append model from file to observation container
  86 + obs.models(gammalib.GModels(model))
  87 +
  88 + # Return observation container
  89 + return obs
  90 +
  91 +
  92 +# ================================ #
  93 +# Simulation and analysis pipeline #
  94 +# ================================ #
  95 +def run_pipeline(obs, ra=83.63, dec=22.01, rad=3.0, \
  96 + emin=0.1, emax=100.0, \
  97 + tmin=0.0, tmax=0.0, \
  98 + debug=False):
  99 + """
  100 + Simulation and unbinned analysis pipeline.
  101 +
  102 + Keywords:
  103 + ra - RA of cube centre [deg] (default: 83.63)
  104 + dec - DEC of cube centre [deg] (default: 22.01)
  105 + rad - Selection radius [deg] (default: 3.0)
  106 + emin - Minimum energy of cube [TeV] (default: 0.1)
  107 + emax - Maximum energy of cube [TeV] (default: 100.0)
  108 + tmin - Start time [MET] (default: 0.0)
  109 + tmax - Stop time [MET] (default: 0.0)
  110 + debug - Enable debugging (default: False)
  111 + """
  112 + # Simulate events
  113 + sim = ctools.ctobssim(obs)
  114 + sim["debug"].boolean(debug)
  115 + sim.run()
  116 +
  117 + # Select events
  118 + select = ctools.ctselect(sim.obs())
  119 + select["ra"].real(ra)
  120 + select["dec"].real(dec)
  121 + select["rad"].real(rad)
  122 + select["emin"].real(emin)
  123 + select["emax"].real(emax)
  124 + select["tmin"].real(tmin)
  125 + select["tmax"].real(tmax)
  126 + select["debug"].boolean(debug)
  127 + select.run()
  128 +
  129 + # Perform maximum likelihood fitting
  130 + like = ctools.ctlike(select.obs())
  131 + like["debug"].boolean(True) # Switch this always on for results in console
  132 + like.run()
  133 +
  134 + # Return
  135 + return
  136 +
  137 +
  138 +# ======================== #
  139 +# Main routine entry point #
  140 +# ======================== #
  141 +if __name__ == '__main__':
  142 + """
  143 + """
  144 + # Dump header
  145 + print("********************************************")
  146 + print("* CTA unbinned analysis pipeline *")
  147 + print("********************************************")
  148 +
  149 + # Setup observations
  150 + obs = setup_observations()
  151 +
  152 + # Setup model
  153 + obs = setup_model(obs)
  154 +
  155 + # Run analysis pipeline
  156 + run_pipeline(obs)
... ...