Peter Deiml / ctools

Commit 337bfcbc994d25f30013974d5d0590095e9e8e99

Authored by knodlseder
1 parent 9be69c6c

initial version

Inline Side-by-side
Showing 1 changed file with 185 additions and 0 deletions
test/test_python.py 0 → 100755
  View file @337bfcb
  1 +#! /usr/bin/env python
  2 +# ===========================================================================================#
  3 +# This script illustrates how to build a unbinned 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 gammalib import *
  9 +from math import *
  10 +import os
  11 +import glob
  12 +import sys
  13 +
  14 +# ============================= #
  15 +# Analysis pipeline - version 1 #
  16 +# ============================= #
  17 +def pipeline_v1():
  18 + """
  19 + Unbinned analysis pipeline - save intermediate results.
  20 +
  21 + This function implements an analysis pipeline that successively calls
  22 + ctobssim, ctselect and ctlike by saving the intermediate results as FITS
  23 + files on disk. This replicates the way how the analysis would be done
  24 + in a ftools-like approach.
  25 + """
  26 + # Set script parameters
  27 + model_name = "data/crab.xml"
  28 + events_name = "events.fits"
  29 + selected_events_name = "selected_events.fits"
  30 + result_name = "results.xml"
  31 + caldb = "irf"
  32 + irf = "kb_E_50h_v3"
  33 + ra = 83.63
  34 + dec = 22.01
  35 + rad_sim = 10.0
  36 + tstart = 0.0
  37 + tstop = 1800.0
  38 + emin = 0.1
  39 + emax = 100.0
  40 + rad_select = 3.0
  41 +
  42 + # Initialise timing
  43 + wall_seconds = 0.0
  44 + cpu_seconds = 0.0
  45 +
  46 + # Simulate events
  47 + sim = ctobssim()
  48 + sim.logFileOpen()
  49 + sim["infile"].filename(model_name)
  50 + sim["outfile"].filename(events_name)
  51 + sim["caldb"].string(caldb)
  52 + sim["irf"].string(irf)
  53 + sim["ra"].real(ra)
  54 + sim["dec"].real(dec)
  55 + sim["rad"].real(rad_sim)
  56 + sim["tmin"].real(tstart)
  57 + sim["tmax"].real(tstop)
  58 + sim["emin"].real(emin)
  59 + sim["emax"].real(emax)
  60 + sim.execute()
  61 + sys.stdout.write(".")
  62 +
  63 + # Select events
  64 + select = ctselect()
  65 + select.logFileOpen()
  66 + select["infile"].filename(events_name)
  67 + select["outfile"].filename(selected_events_name)
  68 + select["ra"].real(ra)
  69 + select["dec"].real(dec)
  70 + select["rad"].real(rad_select)
  71 + select["tmin"].real(tstart)
  72 + select["tmax"].real(tstop)
  73 + select["emin"].real(emin)
  74 + select["emax"].real(emax)
  75 + select.execute()
  76 + sys.stdout.write(".")
  77 +
  78 + # Perform maximum likelihood fitting
  79 + like = ctlike()
  80 + like.logFileOpen()
  81 + like["evfile"].filename(selected_events_name)
  82 + like["srcmdl"].filename(model_name)
  83 + like["outmdl"].filename(result_name)
  84 + like["method"].string("UNBINNED")
  85 + like["caldb"].string(caldb)
  86 + like["irf"].string(irf)
  87 + like.execute()
  88 + sys.stdout.write(".")
  89 +
  90 + # Return
  91 + return
  92 +
  93 +
  94 +# ============================= #
  95 +# Analysis pipeline - version 2 #
  96 +# ============================= #
  97 +def pipeline_v2():
  98 + """
  99 + Unbinned analysis pipeline - keep intermediate results in memory.
  100 +
  101 + This function implements an analysis pipeline that successively calls
  102 + ctobssim, ctselect and ctlike without saving the intermediate results as
  103 + FITS files on disk. All data is only hold in memory.
  104 + """
  105 + # Set script parameters
  106 + model_name = "data/crab.xml"
  107 + caldb = "irf"
  108 + irf = "kb_E_50h_v3"
  109 + ra = 83.63
  110 + dec = 22.01
  111 + rad_sim = 10.0
  112 + tstart = 0.0
  113 + tstop = 1800.0
  114 + emin = 0.1
  115 + emax = 100.0
  116 + rad_select = 3.0
  117 +
  118 + # Initialise timing
  119 + wall_seconds = 0.0
  120 + cpu_seconds = 0.0
  121 +
  122 + # Simulate events
  123 + sim = ctobssim()
  124 + sim["infile"].filename(model_name)
  125 + sim["caldb"].string(caldb)
  126 + sim["irf"].string(irf)
  127 + sim["ra"].real(ra)
  128 + sim["dec"].real(dec)
  129 + sim["rad"].real(rad_sim)
  130 + sim["tmin"].real(tstart)
  131 + sim["tmax"].real(tstop)
  132 + sim["emin"].real(emin)
  133 + sim["emax"].real(emax)
  134 + sim.run()
  135 + sys.stdout.write(".")
  136 +
  137 + # Select events
  138 + select = ctselect(sim.obs())
  139 + select["ra"].real(ra)
  140 + select["dec"].real(dec)
  141 + select["rad"].real(rad_select)
  142 + select["tmin"].real(tstart)
  143 + select["tmax"].real(tstop)
  144 + select["emin"].real(emin)
  145 + select["emax"].real(emax)
  146 + select.run()
  147 + sys.stdout.write(".")
  148 +
  149 + # Perform maximum likelihood fitting
  150 + like = ctlike(select.obs())
  151 + like.run()
  152 + sys.stdout.write(".")
  153 +
  154 + # Return
  155 + return
  156 +
  157 +
  158 +#==========================#
  159 +# Main routine entry point #
  160 +#==========================#
  161 +if __name__ == '__main__':
  162 + """
  163 + Test ctatools Python modules.
  164 + """
  165 + # Dump header
  166 + print ""
  167 + print "*************************************"
  168 + print "* ctatools Python interface testing *"
  169 + print "************************************"
  170 +
  171 + # Remove any existing result files
  172 + list = [glob.glob("*.fits"), glob.glob("*.log"), glob.glob("*.xml")]
  173 + for files in list:
  174 + for file in files:
  175 + os.remove(file)
  176 +
  177 + # Save intermediate results on disk
  178 + print "Test executable analysis: ",
  179 + pipeline_v1()
  180 + print " ok"
  181 +
  182 + # Keep results in memory
  183 + print "Test in-memory analysis: ",
  184 + pipeline_v2()
  185 + print " ok"