Luigi Tibaldo / ctools

Commit 1a1d50d1dd17a99043c9c6d15e25ed5356f230a7

Authored by knodlseder
1 parent 8e6b49b8

renamed from unbinned_analysis.py

Inline Side-by-side
Showing 1 changed file with 427 additions and 0 deletions
examples/make_unbinned_analysis.py 0 → 100755
  View file @1a1d50d
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script illustrates how to build an unbinned analysis pipeline with
  4 +# the ctatools. If matplotlib is installed, the spectrum will be displayed
  5 +# on the screen.
  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 +import sys
  32 +
  33 +
  34 +# ============================= #
  35 +# Analysis pipeline - version 1 #
  36 +# ============================= #
  37 +def pipeline_v1():
  38 + """
  39 + Unbinned analysis pipeline - save intermediate results.
  40 +
  41 + This function implements an analysis pipeline that successively calls
  42 + ctobssim, ctselect and ctlike by saving the intermediate results as FITS
  43 + files on disk. This replicates the way how the analysis would be done
  44 + in a ftools-like approach.
  45 + """
  46 + # Set script parameters
  47 + model_name = "${GAMMALIB}/share/models/crab.xml"
  48 + events_name = "events.fits"
  49 + selected_events_name = "selected_events.fits"
  50 + result_name = "results.xml"
  51 + caldb = "${GAMMALIB}/share/caldb/cta"
  52 + irf = "kb_E_50h_v3"
  53 + ra = 83.63
  54 + dec = 22.01
  55 + rad_sim = 10.0
  56 + tstart = 0.0
  57 + tstop = 1800.0
  58 + emin = 0.1
  59 + emax = 100.0
  60 + rad_select = 3.0
  61 +
  62 + # Initialise timing
  63 + wall_seconds = 0.0
  64 + cpu_seconds = 0.0
  65 +
  66 + # Simulate events
  67 + sim = ctobssim()
  68 + sim.logFileOpen()
  69 + sim["infile"].filename(model_name)
  70 + sim["outfile"].filename(events_name)
  71 + sim["caldb"].string(caldb)
  72 + sim["irf"].string(irf)
  73 + sim["ra"].real(ra)
  74 + sim["dec"].real(dec)
  75 + sim["rad"].real(rad_sim)
  76 + sim["tmin"].real(tstart)
  77 + sim["tmax"].real(tstop)
  78 + sim["emin"].real(emin)
  79 + sim["emax"].real(emax)
  80 + sim.execute()
  81 + print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
  82 +
  83 + # Update timing
  84 + wall_seconds += sim.telapse()
  85 + cpu_seconds += sim.celapse()
  86 +
  87 + # Select events
  88 + select = ctselect()
  89 + select.logFileOpen()
  90 + select["infile"].filename(events_name)
  91 + select["outfile"].filename(selected_events_name)
  92 + select["ra"].real(ra)
  93 + select["dec"].real(dec)
  94 + select["rad"].real(rad_select)
  95 + select["tmin"].real(tstart)
  96 + select["tmax"].real(tstop)
  97 + select["emin"].real(emin)
  98 + select["emax"].real(emax)
  99 + select.execute()
  100 + print "Select event ("+str(select.celapse())+" CPU seconds)"
  101 +
  102 + # Update timing
  103 + wall_seconds += select.telapse()
  104 + cpu_seconds += select.celapse()
  105 +
  106 + # Perform maximum likelihood fitting
  107 + like = ctlike()
  108 + like.logFileOpen()
  109 + like["evfile"].filename(selected_events_name)
  110 + like["srcmdl"].filename(model_name)
  111 + like["outmdl"].filename(result_name)
  112 + like["method"].string("UNBINNED")
  113 + like["caldb"].string(caldb)
  114 + like["irf"].string(irf)
  115 + like.execute()
  116 + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  117 +
  118 + # Update timing
  119 + wall_seconds += like.telapse()
  120 + cpu_seconds += like.celapse()
  121 +
  122 + # Show total times
  123 + print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
  124 + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  125 +
  126 + # Return
  127 + return
  128 +
  129 +
  130 +# ============================= #
  131 +# Analysis pipeline - version 2 #
  132 +# ============================= #
  133 +def pipeline_v2(show_data):
  134 + """
  135 + Unbinned analysis pipeline - keep intermediate results in memory.
  136 +
  137 + This function implements an analysis pipeline that successively calls
  138 + ctobssim, ctselect and ctlike without saving the intermediate results as
  139 + FITS files on disk. All data is only hold in memory.
  140 + """
  141 + # Set script parameters
  142 + model_name = "${GAMMALIB}/share/models/crab.xml"
  143 + caldb = "${GAMMALIB}/share/caldb/cta"
  144 + irf = "kb_E_50h_v3"
  145 + ra = 83.63
  146 + dec = 22.01
  147 + rad_sim = 10.0
  148 + tstart = 0.0
  149 + tstop = 1800.0
  150 + emin = 0.1
  151 + emax = 100.0
  152 + rad_select = 3.0
  153 +
  154 + # Initialise timing
  155 + wall_seconds = 0.0
  156 + cpu_seconds = 0.0
  157 +
  158 + # Simulate events
  159 + sim = ctobssim()
  160 + sim["infile"].filename(model_name)
  161 + sim["caldb"].string(caldb)
  162 + sim["irf"].string(irf)
  163 + sim["ra"].real(ra)
  164 + sim["dec"].real(dec)
  165 + sim["rad"].real(rad_sim)
  166 + sim["tmin"].real(tstart)
  167 + sim["tmax"].real(tstop)
  168 + sim["emin"].real(emin)
  169 + sim["emax"].real(emax)
  170 + sim.run()
  171 + print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
  172 +
  173 + # Update timing
  174 + wall_seconds += sim.telapse()
  175 + cpu_seconds += sim.celapse()
  176 +
  177 + # Select events
  178 + select = ctselect(sim.obs())
  179 + select["ra"].real(ra)
  180 + select["dec"].real(dec)
  181 + select["rad"].real(rad_select)
  182 + select["tmin"].real(tstart)
  183 + select["tmax"].real(tstop)
  184 + select["emin"].real(emin)
  185 + select["emax"].real(emax)
  186 + select.run()
  187 + print "Select event ("+str(select.celapse())+" CPU seconds)"
  188 +
  189 + # Update timing
  190 + wall_seconds += select.telapse()
  191 + cpu_seconds += select.celapse()
  192 +
  193 + # Perform maximum likelihood fitting
  194 + like = ctlike(select.obs())
  195 + like.run()
  196 + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  197 +
  198 + # Update timing
  199 + wall_seconds += like.telapse()
  200 + cpu_seconds += like.celapse()
  201 +
  202 + # Show total times
  203 + print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
  204 + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  205 +
  206 + # Show model fitting results
  207 + #print like.obs().models()
  208 +
  209 + # Optionally plot counts
  210 + if show_data:
  211 + plot_counts(like.obs())
  212 +
  213 + # Return
  214 + return
  215 +
  216 +
  217 +# =========== #
  218 +# Plot counts #
  219 +# =========== #
  220 +def plot_counts(observations):
  221 + """
  222 + Plot counts.
  223 + """
  224 + # Only proceed if matplotlib is available
  225 + try:
  226 + # Import matplotlib
  227 + import matplotlib.pyplot as plt
  228 +
  229 + # Set legend fontsize
  230 + params = {'legend.fontsize': 10}
  231 + plt.rcParams.update(params)
  232 +
  233 + # Set plot styles
  234 + styles = ['b-', 'g-', 'y-', 'n-']
  235 +
  236 + # Dump header
  237 + print ""
  238 + print "Make plots (using matplotlib):"
  239 + print "=============================="
  240 +
  241 + # Create figure 1
  242 + plt.figure(1,figsize=(12,6))
  243 + plt.subplots_adjust(hspace=.7)
  244 +
  245 + # Create subplot 1
  246 + plt.subplot(121)
  247 + plt.title("Spectrum (summed over all pixels)")
  248 +
  249 + # Loop over observations
  250 + for obs in observations:
  251 +
  252 + # Get event list
  253 + list = cast_GCTAEventList(obs.events())
  254 +
  255 + # Create energy axis
  256 + ebounds = GEbounds()
  257 + emin = GEnergy()
  258 + emax = GEnergy()
  259 + emin.TeV(0.1)
  260 + emax.TeV(100.0)
  261 + ebounds.setlog(emin, emax, 10)
  262 + energy = [ebounds.elogmean(i).TeV() for i in range(ebounds.size())]
  263 +
  264 + # Create spectrum
  265 + print "Extract data:"
  266 + counts = [0.0 for i in range(ebounds.size())]
  267 + print list.size()
  268 + for atom in list:
  269 + index = ebounds.index(atom.energy())
  270 + counts[index] = counts[index] + 1.0
  271 +
  272 + # Create error bars
  273 + error = [sqrt(c) for c in counts]
  274 +
  275 + # Plot spectrum
  276 + plt.loglog(energy, counts, 'ro', label='data')
  277 + #plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  278 +
  279 + # Extract models
  280 + print "Extract models:"
  281 + sum_model = [0.0 for i in range(ebounds.size())]
  282 + for k, m in enumerate(observations.models()):
  283 + print "- "+m.name()
  284 + model = [0.0 for i in range(ebounds.size())]
  285 + for atom in list:
  286 + index = ebounds.index(atom.energy())
  287 + prob = m.eval(atom, obs)
  288 + model[index] = model[index] + prob * atom.size()
  289 + for i in range(ebounds.size()):
  290 + sum_model[i] = sum_model[i] + model[i]
  291 + #plt.loglog(energy, model, styles[k], label=m.name())
  292 + #plt.loglog(energy, sum_model, 'r-', label='total')
  293 +
  294 + # Put labels
  295 + plt.xlabel("Energy (TeV)")
  296 + plt.ylabel("Counts")
  297 + plt.legend(loc="lower left")
  298 +
  299 + # Create subplot 2
  300 + plt.subplot(122)
  301 + plt.title("Offset (summed over all energies)")
  302 +
  303 + # Set Crab direction
  304 + crab = GSkyDir()
  305 + crab.radec_deg(83.63, 22.01)
  306 +
  307 + # Loop over observations
  308 + for obs in observations:
  309 +
  310 + # Get event cube
  311 + list = cast_GCTAEventList(obs.events())
  312 +
  313 + # Create offset histogram
  314 + print "Extract data:"
  315 + nx = 30
  316 + doffset2 = 0.01
  317 + offset2 = [(i+0.5)*doffset2 for i in range(nx)]
  318 + counts = [0.0 for i in range(nx)]
  319 + for atom in list:
  320 + off = atom.dir().dist_deg(crab)
  321 + off2 = off*off
  322 + index = int(off2/doffset2)
  323 + if index < nx:
  324 + counts[index] = counts[index] + 1.0
  325 +
  326 + # Create error bars
  327 + error = [sqrt(c) for c in counts]
  328 +
  329 + # Plot distribution
  330 + plt.semilogy(offset2, counts, 'ro', label='data')
  331 +
  332 + # Extract models
  333 + print "Extract models:"
  334 + sum_model = [0.0 for i in range(nx)]
  335 + for k, m in enumerate(observations.models()):
  336 + print "- "+m.name()
  337 + model = [0.0 for i in range(nx)]
  338 + for atom in list:
  339 + off = atom.dir().dist_deg(crab)
  340 + off2 = off*off
  341 + index = int(off2/doffset2)
  342 + if index < nx:
  343 + prob = m.eval(atom, obs)
  344 + model[index] = model[index] + prob * atom.size()
  345 + for i in range(nx):
  346 + sum_model[i] = sum_model[i] + model[i]
  347 + #plt.plot(offset2, model, styles[k], label=m.name())
  348 + #plt.plot(offset2, sum_model, 'r-', label='total')
  349 + #plt.ylim(ymin=0.1)
  350 +
  351 + # Put labels
  352 + plt.xlabel("Offset (deg^2)")
  353 + plt.ylabel("Counts")
  354 + plt.legend(loc="upper right")
  355 +
  356 + # Show counts spectra
  357 + plt.show()
  358 +
  359 + except ImportError:
  360 + print "Matplotlib is not (correctly) installed on your system. No counts spectra are shown."
  361 +
  362 + # Return
  363 + return
  364 +
  365 +
  366 +#==========================#
  367 +# Main routine entry point #
  368 +#==========================#
  369 +if __name__ == '__main__':
  370 + """
  371 + CTA unbinned analysis script.
  372 +
  373 + This script implements the analysis steps to perform an unbinned analysis.
  374 + The steps include:
  375 + - Event simulation
  376 + - Event selection
  377 + - Binned maximum likelihood fitting
  378 +
  379 + Two variations of the script are implemented. One that calls the classes
  380 + like executables, storing the intermediate results on disk, and a second
  381 + that uses the results of the precedent class directly without storing
  382 + any data on disk.
  383 + """
  384 + # Initialise flags
  385 + need_help = False
  386 + show_data = True
  387 +
  388 + # Test for command line arguments
  389 + print sys.argv[0]
  390 + if (len(sys.argv) > 1):
  391 + if sys.argv[1] == "-h":
  392 + need_help = True
  393 + elif sys.argv[1] == "-noshow":
  394 + show_data = False
  395 + else:
  396 + need_help = True
  397 +
  398 + # Print help if needed and exit
  399 + if need_help:
  400 + print "Usage: example_binned.py [OPTIONS]"
  401 + print " -h Display this usage message"
  402 + print " -noshow Do not show data using matplotlib"
  403 + sys.exit()
  404 +
  405 + # Dump header
  406 + print "*********************************"
  407 + print "* CTA unbinned analysis scripts *"
  408 + print "*********************************"
  409 +
  410 + # Remove any existing result files
  411 + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
  412 + for files in list:
  413 + for file in files:
  414 + os.remove(file)
  415 +
  416 + # Save intermediate results on disk
  417 + print ""
  418 + print "Executable analysis pipeline:"
  419 + print "============================="
  420 + pipeline_v1()
  421 +
  422 + # Do analysis in memory
  423 + print ""
  424 + print "In memory analysis pipeline:"
  425 + print "============================"
  426 + pipeline_v2(show_data)
  427 + print ""
... ...