Johan Bregeon / ctools

Commit aa19272d713a9e838b1438fe8b0216dc619c0e6b

Authored by Jürgen Knödlseder
1 parent 210ceab1

Add disk analysis pipelines.

Inline Side-by-side
Showing 4 changed files with 618 additions and 0 deletions
examples/README
  View file @aa19272
@@ -4,14 +4,26 @@ usage of the ctools. @@ -4,14 +4,26 @@ usage of the ctools.
4 benchmark_cta_analysis.py 4 benchmark_cta_analysis.py
5 Performs performance benchmarks of the various CTA analysis types. 5 Performs performance benchmarks of the various CTA analysis types.
6 6
  7 +pipeline_binned_disk.py
  8 + Example pipeline for simulating events and performing a binned
  9 + analysis by saving intermediate results to disk.
  10 +
7 pipeline_binned_mem.py 11 pipeline_binned_mem.py
8 Example pipeline for simulating events and performing a binned 12 Example pipeline for simulating events and performing a binned
9 analysis in memory. No files will be written to disk. 13 analysis in memory. No files will be written to disk.
10 14
  15 +pipeline_stacked_disk.py
  16 + Example pipeline for simulating events and performing a stacked
  17 + analysis by saving intermediate results to disk.
  18 +
11 pipeline_stacked_mem.py 19 pipeline_stacked_mem.py
12 Example pipeline for simulating events and performing a stacked 20 Example pipeline for simulating events and performing a stacked
13 analysis in memory. No files will be written to disk. 21 analysis in memory. No files will be written to disk.
14 22
  23 +pipeline_unbinned_disk.py
  24 + Example pipeline for simulating events and performing an unbinned
  25 + analysis by saving intermediate results to disk.
  26 +
15 pipeline_unbinned_mem.py 27 pipeline_unbinned_mem.py
16 Example pipeline for simulating events and performing an unbinned 28 Example pipeline for simulating events and performing an unbinned
17 analysis in memory. No files will be written to disk. 29 analysis in memory. No files will be written to disk.
examples/pipeline_binned_disk.py 0 → 100755
  View file @aa19272
  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_disk.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", \
  98 + model="${CTOOLS}/share/models/crab.xml", \
  99 + caldb="dummy", irf="cta_dummy_irf", \
  100 + debug=False):
  101 + """
  102 + Simulation and binned analysis pipeline.
  103 +
  104 + Keywords:
  105 + emin - Minimum energy of cube [TeV] (default: 0.1)
  106 + emax - Maximum energy of cube [TeV] (default: 100.0)
  107 + enumbins - Number of energy bins in cube (default: 20)
  108 + nxpix - Number of RA pixels in cube (default: 200)
  109 + nypix - Number of DEC pixels in cube (default: 200)
  110 + binsz - Spatial cube bin size [deg] (default: 0.02)
  111 + coordsys - Cube coordinate system (CEL or GAL)
  112 + proj - Cube World Coordinate System (WCS) projection
  113 + model - Model Xml file
  114 + caldb - Calibration database path (default: "dummy")
  115 + irf - Instrument response function (default: cta_dummy_irf)
  116 + debug - Enable debugging (default: False)
  117 + """
  118 + # Simulate events
  119 + sim = ctools.ctobssim(obs)
  120 + sim["debug"].boolean(debug)
  121 + sim["outevents"].filename("obs.xml")
  122 + sim.execute()
  123 +
  124 + # Bin events by looping over all observations in the container
  125 + sim_obs = gammalib.GObservations("obs.xml")
  126 + obs = gammalib.GObservations()
  127 + for run in sim_obs:
  128 +
  129 + # Get event filename and set counts cube filename
  130 + eventfile = run.eventfile()
  131 + cubefile = "cube_"+eventfile
  132 +
  133 + # Bin events for that observation
  134 + bin = ctools.ctbin()
  135 + bin["inobs"].filename(eventfile)
  136 + bin["outcube"].filename(cubefile)
  137 + bin["ebinalg"].string("LOG")
  138 + bin["emin"].real(emin)
  139 + bin["emax"].real(emax)
  140 + bin["enumbins"].integer(enumbins)
  141 + bin["nxpix"].integer(nxpix)
  142 + bin["nypix"].integer(nypix)
  143 + bin["binsz"].real(binsz)
  144 + bin["coordsys"].string(coordsys)
  145 + bin["usepnt"].boolean(True)
  146 + bin["proj"].string(proj)
  147 + bin.execute()
  148 +
  149 + # Set observation ID
  150 + bin.obs()[0].id(cubefile)
  151 + bin.obs()[0].eventfile(cubefile)
  152 +
  153 + # Append result to observations
  154 + obs.extend(bin.obs())
  155 +
  156 + # Save XML file
  157 + xml = gammalib.GXml()
  158 + obs.write(xml)
  159 + xml.save("obs_cube.xml")
  160 +
  161 + # Perform maximum likelihood fitting
  162 + like = ctools.ctlike()
  163 + like["inobs"].filename("obs_cube.xml")
  164 + like["inmodel"].filename(model)
  165 + like["outmodel"].filename("fit_results.xml")
  166 + like["expcube"].filename("NONE")
  167 + like["psfcube"].filename("NONE")
  168 + like["caldb"].string(caldb)
  169 + like["irf"].string(irf)
  170 + like["debug"].boolean(True) # Switch this always on for results in console
  171 + like.execute()
  172 +
  173 + # Return
  174 + return
  175 +
  176 +
  177 +# ======================== #
  178 +# Main routine entry point #
  179 +# ======================== #
  180 +if __name__ == '__main__':
  181 + """
  182 + """
  183 + # Dump header
  184 + print("********************************************")
  185 + print("* CTA binned analysis pipeline *")
  186 + print("********************************************")
  187 +
  188 + # Setup observations
  189 + obs = setup_observations()
  190 +
  191 + # Setup model
  192 + obs = setup_model(obs)
  193 +
  194 + # Run analysis pipeline
  195 + run_pipeline(obs)
examples/pipeline_stacked_disk.py 0 → 100755
  View file @aa19272
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  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 +#
  7 +# Usage:
  8 +# ./pipeline_stacked_disk.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, emin=0.1, emax=100.0, \
  96 + enumbins=20, nxpix=200, nypix=200, binsz=0.02, \
  97 + coordsys="CEL", proj="CAR", \
  98 + model="${CTOOLS}/share/models/crab.xml", \
  99 + caldb="dummy", irf="cta_dummy_irf", \
  100 + debug=False):
  101 + """
  102 + Simulation and stacked analysis pipeline.
  103 +
  104 + Keywords:
  105 + ra - RA of cube centre [deg] (default: 83.6331)
  106 + dec - DEC of cube centre [deg] (default: 22.0145)
  107 + emin - Minimum energy of cube [TeV] (default: 0.1)
  108 + emax - Maximum energy of cube [TeV] (default: 100.0)
  109 + enumbins - Number of energy bins in cube (default: 20)
  110 + nxpix - Number of RA pixels in cube (default: 200)
  111 + nypix - Number of DEC pixels in cube (default: 200)
  112 + binsz - Spatial cube bin size [deg] (default: 0.02)
  113 + coordsys - Cube coordinate system (CEL or GAL)
  114 + proj - Cube World Coordinate System (WCS) projection
  115 + debug - Enable debugging (default: False)
  116 + """
  117 + # Simulate events
  118 + sim = ctools.ctobssim(obs)
  119 + sim["debug"].boolean(debug)
  120 + sim["outevents"].filename("obs.xml")
  121 + sim.execute()
  122 +
  123 + # Bin events into counts map
  124 + bin = ctools.ctbin()
  125 + bin["inobs"].filename("obs.xml")
  126 + bin["outcube"].filename("cntcube.fits")
  127 + bin["ebinalg"].string("LOG")
  128 + bin["emin"].real(emin)
  129 + bin["emax"].real(emax)
  130 + bin["enumbins"].integer(enumbins)
  131 + bin["nxpix"].integer(nxpix)
  132 + bin["nypix"].integer(nypix)
  133 + bin["binsz"].real(binsz)
  134 + bin["coordsys"].string(coordsys)
  135 + bin["proj"].string(proj)
  136 + bin["xref"].real(ra)
  137 + bin["yref"].real(dec)
  138 + bin["debug"].boolean(debug)
  139 + bin.execute()
  140 +
  141 + # Create exposure cube
  142 + expcube = ctools.ctexpcube()
  143 + expcube["inobs"].filename("obs.xml")
  144 + expcube["incube"].filename("cntcube.fits")
  145 + expcube["outcube"].filename("expcube.fits")
  146 + expcube["caldb"].string(caldb)
  147 + expcube["irf"].string(irf)
  148 + expcube["ebinalg"].string("LOG")
  149 + expcube["emin"].real(emin)
  150 + expcube["emax"].real(emax)
  151 + expcube["enumbins"].integer(enumbins)
  152 + expcube["nxpix"].integer(nxpix)
  153 + expcube["nypix"].integer(nypix)
  154 + expcube["binsz"].real(binsz)
  155 + expcube["coordsys"].string(coordsys)
  156 + expcube["proj"].string(proj)
  157 + expcube["xref"].real(ra)
  158 + expcube["yref"].real(dec)
  159 + expcube["debug"].boolean(debug)
  160 + expcube.execute()
  161 +
  162 + # Create PSF cube
  163 + psfcube = ctools.ctpsfcube()
  164 + psfcube["inobs"].filename("obs.xml")
  165 + psfcube["incube"].filename("NONE")
  166 + psfcube["outcube"].filename("psfcube.fits")
  167 + psfcube["caldb"].string(caldb)
  168 + psfcube["irf"].string(irf)
  169 + psfcube["ebinalg"].string("LOG")
  170 + psfcube["emin"].real(emin)
  171 + psfcube["emax"].real(emax)
  172 + psfcube["enumbins"].integer(enumbins)
  173 + psfcube["nxpix"].integer(10)
  174 + psfcube["nypix"].integer(10)
  175 + psfcube["binsz"].real(1.0)
  176 + psfcube["coordsys"].string(coordsys)
  177 + psfcube["proj"].string(proj)
  178 + psfcube["xref"].real(ra)
  179 + psfcube["yref"].real(dec)
  180 + psfcube["debug"].boolean(debug)
  181 + psfcube.execute()
  182 +
  183 + # Create background cube
  184 + bkgcube = ctools.ctbkgcube()
  185 + bkgcube["inobs"].filename("obs.xml")
  186 + bkgcube["inmodel"].filename(model)
  187 + bkgcube["incube"].filename("cntcube.fits")
  188 + bkgcube["outcube"].filename("bkgcube.fits")
  189 + bkgcube["outmodel"].filename("model_bkg.xml")
  190 + bkgcube["caldb"].string(caldb)
  191 + bkgcube["irf"].string(irf)
  192 + bkgcube["ebinalg"].string("LOG")
  193 + bkgcube["emin"].real(emin)
  194 + bkgcube["emax"].real(emax)
  195 + bkgcube["enumbins"].integer(enumbins)
  196 + bkgcube["nxpix"].integer(10)
  197 + bkgcube["nypix"].integer(10)
  198 + bkgcube["binsz"].real(1.0)
  199 + bkgcube["coordsys"].string(coordsys)
  200 + bkgcube["proj"].string(proj)
  201 + bkgcube["xref"].real(ra)
  202 + bkgcube["yref"].real(dec)
  203 + bkgcube["debug"].boolean(debug)
  204 + bkgcube.execute()
  205 +
  206 + # Perform maximum likelihood fitting
  207 + like = ctools.ctlike()
  208 + like["inobs"].filename("cntcube.fits")
  209 + like["inmodel"].filename("model_bkg.xml")
  210 + like["outmodel"].filename("fit_results.xml")
  211 + like["expcube"].filename("expcube.fits")
  212 + like["psfcube"].filename("psfcube.fits")
  213 + like["caldb"].string(caldb)
  214 + like["irf"].string(irf)
  215 + like["debug"].boolean(True) # Switch this always on for results in console
  216 + like.execute()
  217 +
  218 + # Return
  219 + return
  220 +
  221 +
  222 +# ======================== #
  223 +# Main routine entry point #
  224 +# ======================== #
  225 +if __name__ == '__main__':
  226 + """
  227 + """
  228 + # Dump header
  229 + print("********************************************")
  230 + print("* CTA stacked analysis pipeline *")
  231 + print("********************************************")
  232 +
  233 + # Setup observations
  234 + obs = setup_observations()
  235 +
  236 + # Setup model
  237 + obs = setup_model(obs)
  238 +
  239 + # Run analysis pipeline
  240 + run_pipeline(obs)
examples/pipeline_unbinned_disk.py 0 → 100755
  View file @aa19272
  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_disk.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 + model="${CTOOLS}/share/models/crab.xml", \
  99 + caldb="dummy", irf="cta_dummy_irf", \
  100 + debug=False):
  101 + """
  102 + Simulation and unbinned analysis pipeline.
  103 +
  104 + Keywords:
  105 + ra - RA of cube centre [deg] (default: 83.63)
  106 + dec - DEC of cube centre [deg] (default: 22.01)
  107 + rad - Selection radius [deg] (default: 3.0)
  108 + emin - Minimum energy of cube [TeV] (default: 0.1)
  109 + emax - Maximum energy of cube [TeV] (default: 100.0)
  110 + tmin - Start time [MET] (default: 0.0)
  111 + tmax - Stop time [MET] (default: 0.0)
  112 + model - Model Xml file
  113 + caldb - Calibration database path (default: "dummy")
  114 + irf - Instrument response function (default: cta_dummy_irf)
  115 + debug - Enable debugging (default: False)
  116 + """
  117 + # Get model
  118 +
  119 + # Simulate events
  120 + sim = ctools.ctobssim(obs)
  121 + sim["debug"].boolean(debug)
  122 + sim["outevents"].filename("obs.xml")
  123 + sim.execute()
  124 +
  125 + # Select events
  126 + select = ctools.ctselect()
  127 + select["inobs"].filename("obs.xml")
  128 + select["outobs"].filename("obs_selected.xml")
  129 + select["ra"].real(ra)
  130 + select["dec"].real(dec)
  131 + select["rad"].real(rad)
  132 + select["emin"].real(emin)
  133 + select["emax"].real(emax)
  134 + select["tmin"].real(tmin)
  135 + select["tmax"].real(tmax)
  136 + select["debug"].boolean(debug)
  137 + select.execute()
  138 +
  139 + # Perform maximum likelihood fitting
  140 + like = ctools.ctlike()
  141 + like["inobs"].filename("obs_selected.xml")
  142 + like["inmodel"].filename(model)
  143 + like["outmodel"].filename("fit_results.xml")
  144 + like["caldb"].string(caldb)
  145 + like["irf"].string(irf)
  146 + like["debug"].boolean(True) # Switch this always on for results in console
  147 + like.execute()
  148 +
  149 + # Return
  150 + return
  151 +
  152 +
  153 +# ======================== #
  154 +# Main routine entry point #
  155 +# ======================== #
  156 +if __name__ == '__main__':
  157 + """
  158 + """
  159 + # Dump header
  160 + print("********************************************")
  161 + print("* CTA unbinned analysis pipeline *")
  162 + print("********************************************")
  163 +
  164 + # Setup observations
  165 + obs = setup_observations()
  166 +
  167 + # Setup model
  168 + obs = setup_model(obs)
  169 +
  170 + # Run analysis pipeline
  171 + run_pipeline(obs)