Konstancja Satalecka / ctools

Commit 557ef70ec55523369fd1759a77e76ae4ffd61ed4

Authored by knodlseder
1 parent d02652da

initial version

Inline Side-by-side
Showing 1 changed file with 128 additions and 0 deletions
test/example_models.py 0 → 100755
  View file @557ef70
  1 +#! /usr/bin/env python
  2 +# ===========================================================================================#
  3 +# This script illustrates how to build a unbinned analysis chain using CTAtools.
  4 +# ===========================================================================================#
  5 +from ctatools import *
  6 +from gammalib import *
  7 +from math import *
  8 +import os
  9 +import glob
  10 +import sys
  11 +
  12 +
  13 +# ================= #
  14 +# Analysis pipeline #
  15 +# ================= #
  16 +def pipeline(model_name):
  17 + """
  18 + Unbinned analysis pipeline - keep intermediate results in memory.
  19 +
  20 + This function implements an analysis pipeline that successively calls
  21 + ctobssim, ctselect and ctlike without saving the intermediate results as
  22 + FITS files on disk. All data is only hold in memory.
  23 + """
  24 + # Set script parameters
  25 + caldb = "irf"
  26 + irf = "kb_E_50h_v3"
  27 + ra = 83.63
  28 + dec = 22.01
  29 + rad_sim = 10.0
  30 + tstart = 0.0
  31 + tstop = 1800.0
  32 + emin = 0.1
  33 + emax = 100.0
  34 + rad_select = 3.0
  35 +
  36 + # Write model name
  37 + print "Model: "+model_name
  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["caldb"].string(caldb)
  47 + sim["irf"].string(irf)
  48 + sim["ra"].real(ra)
  49 + sim["dec"].real(dec)
  50 + sim["rad"].real(rad_sim)
  51 + sim["tmin"].real(tstart)
  52 + sim["tmax"].real(tstop)
  53 + sim["emin"].real(emin)
  54 + sim["emax"].real(emax)
  55 + sim.run()
  56 +
  57 + # Select events
  58 + select = ctselect(sim.obs())
  59 + select["ra"].real(ra)
  60 + select["dec"].real(dec)
  61 + select["rad"].real(rad_select)
  62 + select["tmin"].real(tstart)
  63 + select["tmax"].real(tstop)
  64 + select["emin"].real(emin)
  65 + select["emax"].real(emax)
  66 + select.run()
  67 +
  68 + # Perform maximum likelihood fitting
  69 + like = ctlike(select.obs())
  70 + like.run()
  71 +
  72 + # Show model fitting results
  73 + print like.obs().models()[0]
  74 +
  75 + # Return
  76 + return
  77 +
  78 +
  79 +#==========================#
  80 +# Main routine entry point #
  81 +#==========================#
  82 +if __name__ == '__main__':
  83 + """
  84 + Perform unbinned analyses for various models.
  85 +
  86 + This script implements the analysis steps to perform an unbinned analysis.
  87 + The steps include:
  88 + - Event simulation
  89 + - Event selection
  90 + - Binned maximum likelihood fitting
  91 + """
  92 + # Initialise flags
  93 + need_help = False
  94 + show_data = True
  95 +
  96 + # Test for command line arguments
  97 + print sys.argv[0]
  98 + if (len(sys.argv) > 1):
  99 + if sys.argv[1] == "-h":
  100 + need_help = True
  101 + elif sys.argv[1] == "-noshow":
  102 + show_data = False
  103 + else:
  104 + need_help = True
  105 +
  106 + # Print help if needed and exit
  107 + if need_help:
  108 + print "Usage: example_models.py [OPTIONS]"
  109 + print " -h Display this usage message"
  110 + print " -noshow Do not show data using matplotlib"
  111 + sys.exit()
  112 +
  113 + # Dump header
  114 + print "********************************"
  115 + print "* CTA unbinned analysis script *"
  116 + print "********************************"
  117 +
  118 + # Perform analysis for Crab model
  119 + pipeline("data/crab.xml")
  120 +
  121 + # Perform analysis for disk model
  122 + pipeline("data/disk.xml")
  123 +
  124 + # Perform analysis for Gaussian model
  125 + pipeline("data/gauss.xml")
  126 +
  127 + # Perform analysis for shell model
  128 + pipeline("data/shell.xml")
... ...