Leonardo Di Venere / ctools

Commit 437ae5426b619787275c918751cbac7afbc34d19

Authored by knodlseder
1 parent 531864e0

initial version

Inline Side-by-side
Showing 1 changed file with 203 additions and 0 deletions
test/example_binned.py 0 → 100755
  View file @437ae54
  1 +#! /usr/bin/env python
  2 +# ===========================================================================================#
  3 +# This script illustrates how to build a binned analysis chain using CTAtools.
  4 +#
  5 +# If matplotlib is installed, the spectrum will be displayed on the screen.
  6 +# ===========================================================================================#
  7 +from ctatools import *
  8 +from math import *
  9 +
  10 +
  11 +# =================== #
  12 +# Processing pipeline #
  13 +# =================== #
  14 +def process():
  15 + """
  16 + Binned analysis processing pipeline.
  17 + """
  18 + # Set script parameters
  19 + model_name = "data/crab.xml"
  20 + events_name = "events.fits"
  21 + cntmap_name = "cntmap.fits"
  22 + result_name = "results.xml"
  23 + caldb = "irf"
  24 + irf = "kb_E_50h_v3"
  25 + ra = 83.63
  26 + dec = 22.01
  27 + rad_sim = 10.0
  28 + tstart = 0.0
  29 + tstop = 1800.0
  30 + emin = 0.1
  31 + emax = 100.0
  32 + enumbins = 20
  33 + nxpix = 200
  34 + nypix = 200
  35 + binsz = 0.02
  36 + coordsys = "CEL"
  37 + proj = "CAR"
  38 +
  39 + # Initialise timing
  40 + wall_seconds = 0.0
  41 + cpu_seconds = 0.0
  42 +
  43 + # Simulate events
  44 + sim = ctobssim()
  45 + sim["infile"].filename(model_name)
  46 + sim["outfile"].filename(events_name)
  47 + sim["caldb"].string(caldb)
  48 + sim["irf"].string(irf)
  49 + sim["ra"].real(ra)
  50 + sim["dec"].real(dec)
  51 + sim["rad"].real(rad_sim)
  52 + sim["tmin"].real(tstart)
  53 + sim["tmax"].real(tstop)
  54 + sim["emin"].real(emin)
  55 + sim["emax"].real(emax)
  56 + sim.run()
  57 + print "Simulated events ("+str(sim.celapse())+" CPU seconds)"
  58 +
  59 + # Update timing
  60 + wall_seconds += sim.telapse()
  61 + cpu_seconds += sim.celapse()
  62 +
  63 + # Bin events into counts map
  64 + bin = ctbin()
  65 + bin["evfile"].filename(events_name)
  66 + bin["outfile"].filename(cntmap_name)
  67 + bin["emin"].real(emin)
  68 + bin["emax"].real(emax)
  69 + bin["enumbins"].integer(enumbins)
  70 + bin["nxpix"].integer(nxpix)
  71 + bin["nypix"].integer(nypix)
  72 + bin["binsz"].real(binsz)
  73 + bin["coordsys"].string(coordsys)
  74 + bin["xref"].real(ra)
  75 + bin["yref"].real(dec)
  76 + bin["proj"].string(proj)
  77 + bin.run()
  78 + print "Binned events into counts map ("+str(bin.celapse())+" CPU seconds)"
  79 +
  80 + # Update timing
  81 + wall_seconds += bin.telapse()
  82 + cpu_seconds += bin.celapse()
  83 +
  84 + # Perform maximum likelihood fitting
  85 + like = ctlike()
  86 + like["cntmap"].filename(cntmap_name)
  87 + like["srcmdl"].filename(model_name)
  88 + like["outmdl"].filename(result_name)
  89 + like["method"].string("BINNED")
  90 + like["stat"].string("POISSON")
  91 + like["caldb"].string(caldb)
  92 + like["irf"].string(irf)
  93 + like.run()
  94 + print "Maximum likelihood fitting ("+str(like.celapse())+" CPU seconds)"
  95 +
  96 + # Update timing
  97 + wall_seconds += like.telapse()
  98 + cpu_seconds += like.celapse()
  99 +
  100 + # Show total times
  101 + print "Total wall time elapsed: "+str(wall_seconds)+" seconds"
  102 + print "Total CPU time used ...: "+str(cpu_seconds)+" seconds"
  103 +
  104 + # Return
  105 + return
  106 +
  107 +
  108 +# ============ #
  109 +# Show photons #
  110 +# ============ #
  111 +def show_photons(photons, xmlname, e_min, e_max, area, duration, ebins=30):
  112 + """
  113 + Show photons using matplotlib (if available).
  114 + """
  115 + # Only proceed if matplotlib is available
  116 + try:
  117 + # Import matplotlib
  118 + import matplotlib.pyplot as plt
  119 +
  120 + # Create figure
  121 + plt.figure(1)
  122 + plt.title("MC simulated photon spectrum ("+str(e_min)+'-'+str(e_max)+" TeV)")
  123 +
  124 + # Setup energy range covered by data
  125 + emin = GEnergy()
  126 + emax = GEnergy()
  127 + emin.TeV(e_min)
  128 + emax.TeV(e_max)
  129 + ebds = GEbounds()
  130 + ebds.setlog(emin, emax, ebins)
  131 +
  132 + # Create energy axis
  133 + energy = []
  134 + for i in range(ebds.size()):
  135 + energy.append(ebds.elogmean(i).TeV())
  136 +
  137 + # Fill histogram
  138 + counts = [0.0 for i in range(ebds.size())]
  139 + for photon in photons:
  140 + index = ebds.index(photon.energy())
  141 + counts[index] = counts[index] + 1.0
  142 +
  143 + # Create error bars
  144 + error = [sqrt(c) for c in counts]
  145 +
  146 + # Get model values
  147 + models = GModels(xmlname)
  148 + crab = cast_GModelSky(models[0])
  149 + model = []
  150 + d = GSkyDir()
  151 + d.radec_deg(83.6331, 22.0145)
  152 + t = GTime()
  153 + for i in range(ebds.size()):
  154 + eval = ebds.elogmean(i)
  155 + ewidth = ebds.emax(i) - ebds.emin(i)
  156 + f = crab.value(d, eval, t) * area * duration * ewidth.MeV()
  157 + model.append(f)
  158 +
  159 + # Plot data
  160 + plt.loglog(energy, counts, 'ro')
  161 + plt.errorbar(energy, counts, error, fmt=None, ecolor='r')
  162 +
  163 + # Plot model
  164 + plt.plot(energy, model, 'b-')
  165 +
  166 + # Set axes
  167 + plt.xlabel("Energy (TeV)")
  168 + plt.ylabel("Number of incident photons")
  169 +
  170 + # Allocate histogram
  171 + # Show plot
  172 + plt.show()
  173 +
  174 + except:
  175 + print "Matplotlib is not (correctly) installed on your system. No data are shown."
  176 +
  177 + # Return
  178 + return
  179 +
  180 +
  181 +#==========================#
  182 +# Main routine entry point #
  183 +#==========================#
  184 +if __name__ == '__main__':
  185 + """
  186 + CTA binned analysis script.
  187 +
  188 + This script implements the analysis steps to perform a binned analysis.
  189 + The steps include:
  190 + - Event simulation
  191 + - Event binning
  192 + - Binned maximum likelihood fitting
  193 + """
  194 + # Dump header
  195 + print "*****************************"
  196 + print "* CTA binned analysis script *"
  197 + print "*****************************"
  198 +
  199 + # Process data
  200 + process()
  201 +
  202 + # Show photons
  203 + #show_photons(photons, xmlname, e_min, e_max, area, duration)
... ...