Fabio Iocco / ctools

Commit 910cb4f3baa4faa0ce0f4933e293f6f8b48c6852

Authored by knodlseder
1 parent 1a1d50d1

renamed from binned_analysis.py

Inline Side-by-side
Showing 1 changed file with 424 additions and 0 deletions
examples/make_binned_analysis.py 0 → 100755
  View file @910cb4f
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script illustrates how to build a binned analysis pipeline using the
  4 +# ctatools. If matplotlib is installed, a counts spectrum and an offset
  5 +# histogram will be displayed.
  6 +#
  7 +# Required 3rd party modules:
  8 +# None
  9 +#
  10 +# Copyright (C) 2011 Jurgen Knodlseder
  11 +#
  12 +# This program is free software: you can redistribute it and/or modify
  13 +# it under the terms of the GNU General Public License as published by
  14 +# the Free Software Foundation, either version 3 of the License, or
  15 +# (at your option) any later version.
  16 +#
  17 +# This program is distributed in the hope that it will be useful,
  18 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  20 +# GNU General Public License for more details.
  21 +#
  22 +# You should have received a copy of the GNU General Public License
  23 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  24 +#
  25 +# ==========================================================================
  26 +from ctatools import *
  27 +from gammalib import *
  28 +from math import *
  29 +import os
  30 +import glob
  31 +
  32 +
  33 +# ============================= #
  34 +# Analysis pipeline - version 1 #
  35 +# ============================= #
  36 +def pipeline_v1():
  37 + """
  38 + Binned analysis pipeline - save intermediate results.
  39 +
  40 + This function implements an analysis pipeline that successively calls
  41 + ctobssim, ctbin and ctlike by saving the intermediate results as FITS
  42 + files on disk. This replicates the way how the analysis would be done
  43 + in a ftools-like approach.
  44 + """
  45 + # Set script parameters
  46 + model_name = "${GAMMALIB}/share/models/crab.xml"
  47 + events_name = "events.fits"
  48 + cntmap_name = "cntmap.fits"
  49 + result_name = "results.xml"
  50 + caldb = "${GAMMALIB}/share/caldb/cta"
  51 + irf = "kb_E_50h_v3"
  52 + ra = 83.63
  53 + dec = 22.01
  54 + rad_sim = 10.0
  55 + tstart = 0.0
  56 + tstop = 1800.0
  57 + emin = 0.1
  58 + emax = 100.0
  59 + enumbins = 20
  60 + nxpix = 200
  61 + nypix = 200
  62 + binsz = 0.02
  63 + coordsys = "CEL"
  64 + proj = "CAR"
  65 +
  66 + # Initialise timing
  67 + wall_seconds = 0.0
  68 + cpu_seconds = 0.0
  69 +
  70 + # Simulate events
  71 + sim = ctobssim()
  72 + sim.logFileOpen() # We need this to explicitely open the log file in Python mode
  73 + sim["infile"].filename(model_name)
  74 + sim["outfile"].filename(events_name)
  75 + sim["caldb"].string(caldb)
  76 + sim["irf"].string(irf)
  77 + sim["ra"].real(ra)
  78 + sim["dec"].real(dec)
  79 + sim["rad"].real(rad_sim)
  80 + sim["tmin"].real(tstart)
  81 + sim["tmax"].real(tstop)
  82 + sim["emin"].real(emin)
  83 + sim["emax"].real(emax)
  84 + sim.execute()
  85 + print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
  86 +
  87 + # Update timing
  88 + wall_seconds += sim.telapse()
  89 + cpu_seconds += sim.celapse()
  90 +
  91 + # Bin events into counts map
  92 + bin = ctbin()
  93 + bin.logFileOpen() # We need this to explicitely open the log file in Python mode
  94 + bin["evfile"].filename(events_name)
  95 + bin["outfile"].filename(cntmap_name)
  96 + bin["emin"].real(emin)
  97 + bin["emax"].real(emax)
  98 + bin["enumbins"].integer(enumbins)
  99 + bin["nxpix"].integer(nxpix)
  100 + bin["nypix"].integer(nypix)
  101 + bin["binsz"].real(binsz)
  102 + bin["coordsys"].string(coordsys)
  103 + bin["xref"].real(ra)
  104 + bin["yref"].real(dec)
  105 + bin["proj"].string(proj)
  106 + bin.execute()
  107 + print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
  108 +
  109 + # Update timing
  110 + wall_seconds += bin.telapse()
  111 + cpu_seconds += bin.celapse()
  112 +
  113 + # Perform maximum likelihood fitting
  114 + like = ctlike()
  115 + like.logFileOpen() # We need this to explicitely open the log file in Python mode
  116 + like["cntmap"].filename(cntmap_name)
  117 + like["srcmdl"].filename(model_name)
  118 + like["outmdl"].filename(result_name)
  119 + like["method"].string("BINNED")
  120 + like["caldb"].string(caldb)
  121 + like["irf"].string(irf)
  122 + like.execute()
  123 + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  124 +
  125 + # Update timing
  126 + wall_seconds += like.telapse()
  127 + cpu_seconds += like.celapse()
  128 +
  129 + # Show total times
  130 + print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
  131 + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  132 +
  133 + # Return
  134 + return
  135 +
  136 +
  137 +# ============================= #
  138 +# Analysis pipeline - version 2 #
  139 +# ============================= #
  140 +def pipeline_v2():
  141 + """
  142 + Binned analysis pipeline - keep intermediate results in memory.
  143 +
  144 + This function implements an analysis pipeline that successively calls
  145 + ctobssim, ctbin and ctlike without saving the intermediate results as
  146 + FITS files on disk. All data are hold in memory.
  147 +
  148 + At the end, results are plotted (if matplotlib is installed)
  149 + """
  150 + # Set script parameters
  151 + model_name = "${GAMMALIB}/share/models/crab.xml"
  152 + caldb = "${GAMMALIB}/share/caldb/cta"
  153 + irf = "kb_E_50h_v3"
  154 + ra = 83.63
  155 + dec = 22.01
  156 + rad_sim = 10.0
  157 + tstart = 0.0
  158 + tstop = 1800.0
  159 + emin = 0.1
  160 + emax = 100.0
  161 + enumbins = 20
  162 + nxpix = 200
  163 + nypix = 200
  164 + binsz = 0.02
  165 + coordsys = "CEL"
  166 + proj = "CAR"
  167 +
  168 + # Initialise timing
  169 + wall_seconds = 0.0
  170 + cpu_seconds = 0.0
  171 +
  172 + # Simulate events
  173 + sim = ctobssim()
  174 + sim["infile"].filename(model_name)
  175 + sim["caldb"].string(caldb)
  176 + sim["irf"].string(irf)
  177 + sim["ra"].real(ra)
  178 + sim["dec"].real(dec)
  179 + sim["rad"].real(rad_sim)
  180 + sim["tmin"].real(tstart)
  181 + sim["tmax"].real(tstop)
  182 + sim["emin"].real(emin)
  183 + sim["emax"].real(emax)
  184 + sim.run()
  185 + print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
  186 +
  187 + # Update timing
  188 + wall_seconds += sim.telapse()
  189 + cpu_seconds += sim.celapse()
  190 +
  191 + # Bin events into counts map
  192 + bin = ctbin(sim.obs())
  193 + bin["emin"].real(emin)
  194 + bin["emax"].real(emax)
  195 + bin["enumbins"].integer(enumbins)
  196 + bin["nxpix"].integer(nxpix)
  197 + bin["nypix"].integer(nypix)
  198 + bin["binsz"].real(binsz)
  199 + bin["coordsys"].string(coordsys)
  200 + bin["xref"].real(ra)
  201 + bin["yref"].real(dec)
  202 + bin["proj"].string(proj)
  203 + bin.run()
  204 + print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
  205 +
  206 + # Update timing
  207 + wall_seconds += bin.telapse()
  208 + cpu_seconds += bin.celapse()
  209 +
  210 + # Perform maximum likelihood fitting
  211 + like = ctlike(bin.obs())
  212 + like.run()
  213 + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  214 +
  215 + # Update timing
  216 + wall_seconds += like.telapse()
  217 + cpu_seconds += like.celapse()
  218 +
  219 + # Show total times
  220 + print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
  221 + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  222 +
  223 + # Show model fitting results
  224 + #print like.obs().models()
  225 +
  226 + # Plot counts
  227 + plot_counts(bin.obs())
  228 +
  229 + # Return
  230 + return
  231 +
  232 +
  233 +# =========== #
  234 +# Plot counts #
  235 +# =========== #
  236 +def plot_counts(observations):
  237 + """
  238 + Plot counts.
  239 + """
  240 + # Only proceed if matplotlib is available
  241 + try:
  242 + # Import matplotlib
  243 + import matplotlib.pyplot as plt
  244 +
  245 + # Set legend fontsize
  246 + params = {'legend.fontsize': 10}
  247 + plt.rcParams.update(params)
  248 +
  249 + # Set plot styles
  250 + styles = ['b-', 'g-', 'y-', 'n-']
  251 +
  252 + # Dump header
  253 + print ""
  254 + print "Make plots (using matplotlib):"
  255 + print "=============================="
  256 +
  257 + # Create figure 1
  258 + plt.figure(1,figsize=(12,6))
  259 + plt.subplots_adjust(hspace=.7)
  260 +
  261 + # Create subplot 1
  262 + plt.subplot(121)
  263 + plt.title("Spectrum (summed over all pixels)")
  264 +
  265 + # Loop over observations
  266 + for obs in observations:
  267 +
  268 + # Get event cube
  269 + cube = cast_GCTAEventCube(obs.events())
  270 +
  271 + # Create energy axis
  272 + energy = []
  273 + ebounds = cube.ebounds()
  274 + for i in range(ebounds.size()):
  275 + energy.append(ebounds.elogmean(i).TeV())
  276 +
  277 + # Create spectrum
  278 + print "Extract data:"
  279 + counts = [0.0 for i in range(ebounds.size())]
  280 + for bin in cube:
  281 + index = ebounds.index(bin.energy())
  282 + counts[index] = counts[index] + bin.counts()
  283 +
  284 + # Create error bars
  285 + error = [sqrt(c) for c in counts]
  286 +
  287 + # Plot spectrum
  288 + plt.loglog(energy, counts, 'ro', label='data')
  289 + plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  290 +
  291 + # Extract models
  292 + print "Extract models:"
  293 + sum_model = [0.0 for i in range(ebounds.size())]
  294 + for k, m in enumerate(observations.models()):
  295 + print "- "+m.name()
  296 + model = [0.0 for i in range(ebounds.size())]
  297 + for bin in cube:
  298 + index = ebounds.index(bin.energy())
  299 + prob = m.eval(bin, obs)
  300 + model[index] = model[index] + prob * bin.size()
  301 + for i in range(ebounds.size()):
  302 + sum_model[i] = sum_model[i] + model[i]
  303 + plt.loglog(energy, model, styles[k], label=m.name())
  304 + plt.loglog(energy, sum_model, 'r-', label='total')
  305 +
  306 + # Put labels
  307 + plt.xlabel("Energy (TeV)")
  308 + plt.ylabel("Counts")
  309 + plt.legend(loc="lower left")
  310 +
  311 + # Create subplot 2
  312 + plt.subplot(122)
  313 + plt.title("Offset (summed over all energies)")
  314 +
  315 + # Set Crab direction
  316 + crab = GSkyDir()
  317 + crab.radec_deg(83.63, 22.01)
  318 +
  319 + # Loop over observations
  320 + for obs in observations:
  321 +
  322 + # Get event cube
  323 + cube = cast_GCTAEventCube(obs.events())
  324 +
  325 + # Create offset histogram
  326 + print "Extract data:"
  327 + nx = 30
  328 + doffset2 = 0.01
  329 + offset2 = [(i+0.5)*doffset2 for i in range(nx)]
  330 + counts = [0.0 for i in range(nx)]
  331 + for bin in cube:
  332 + off = bin.dir().dist_deg(crab)
  333 + off2 = off*off
  334 + index = int(off2/doffset2)
  335 + if index < nx:
  336 + counts[index] = counts[index] + bin.counts()
  337 +
  338 + # Create error bars
  339 + error = [sqrt(c) for c in counts]
  340 +
  341 + # Plot distribution
  342 + plt.semilogy(offset2, counts, 'ro', label='data')
  343 +
  344 + # Extract models
  345 + print "Extract models:"
  346 + sum_model = [0.0 for i in range(nx)]
  347 + for k, m in enumerate(observations.models()):
  348 + print "- "+m.name()
  349 + model = [0.0 for i in range(nx)]
  350 + for bin in cube:
  351 + off = bin.dir().dist_deg(crab)
  352 + off2 = off*off
  353 + index = int(off2/doffset2)
  354 + if index < nx:
  355 + prob = m.eval(bin, obs)
  356 + model[index] = model[index] + prob * bin.size()
  357 + for i in range(nx):
  358 + sum_model[i] = sum_model[i] + model[i]
  359 + plt.plot(offset2, model, styles[k], label=m.name())
  360 + plt.plot(offset2, sum_model, 'r-', label='total')
  361 + plt.ylim(ymin=0.1)
  362 +
  363 + # Put labels
  364 + plt.xlabel("Offset (deg^2)")
  365 + plt.ylabel("Counts")
  366 + plt.legend(loc="upper right")
  367 +
  368 + # Show counts spectra
  369 + plt.show()
  370 +
  371 + except:
  372 + print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
  373 +
  374 + # Return
  375 + return
  376 +
  377 +
  378 +#==========================#
  379 +# Main routine entry point #
  380 +#==========================#
  381 +if __name__ == '__main__':
  382 + """
  383 + CTA binned analysis script.
  384 +
  385 + This script implements the analysis steps to perform a binned analysis.
  386 + The steps include:
  387 + - Event simulation
  388 + - Event binning
  389 + - Binned maximum likelihood fitting
  390 +
  391 + Two variations of the script are implemented. One that calls the classes
  392 + like executables, storing the intermediate results on disk, and a second
  393 + that uses the results of the precedent class directly without storing
  394 + any data on disk.
  395 +
  396 + The second variation will also show a counts spectrum of the data after
  397 + fitting (only if matplotlib is available on your system). Note that the
  398 + counts spectrum is summed over the entire data space, i.e. no event selection
  399 + around the Crab has been performed.
  400 + """
  401 + # Dump header
  402 + print "*************************************"
  403 + print "* CTA binned analysis scripts *"
  404 + print "*************************************"
  405 + print "... this script will take 1-2 minutes"
  406 +
  407 + # Remove any existing result files
  408 + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
  409 + for files in list:
  410 + for file in files:
  411 + os.remove(file)
  412 +
  413 + # Save intermediate results on disk
  414 + print ""
  415 + print "Executable analysis pipeline:"
  416 + print "============================="
  417 + pipeline_v1()
  418 +
  419 + # Do analysis in memory
  420 + print ""
  421 + print "In memory analysis pipeline:"
  422 + print "============================"
  423 + pipeline_v2()
  424 + print ""