Leonardo Di Venere / ctools

Commit 9048805f303b11d51355feb28d9314d42df3cb70

Authored by knodlseder
1 parent 7e01dd12

initial version

Inline Side-by-side
Showing 1 changed file with 132 additions and 0 deletions
test/test_coverage.py 0 → 100755
  View file @9048805
  1 +#! /usr/bin/env python
  2 +# ===========================================================================================#
  3 +# This script illustrates how to build a binned analysis pipeline using CTAtools.
  4 +#
  5 +# If matplotlib is installed, a counts spectrum and an offset histogram will be displayed.
  6 +# ===========================================================================================#
  7 +from ctatools import *
  8 +from gammalib import *
  9 +from math import *
  10 +import os
  11 +import glob
  12 +import sys
  13 +import csv
  14 +
  15 +# Parameters
  16 +events_name = "events.fits"
  17 +cntmap_name = "cntmap.fits"
  18 +result_name = "results.xml"
  19 +caldb = "irf"
  20 +irf = "kb_E_50h_v3"
  21 +ra = 83.63
  22 +dec = 22.01
  23 +rad_sim = 10.0
  24 +tstart = 0.0
  25 +tstop = 1e4
  26 +emin = 1
  27 +emax = 10
  28 +enumbins = 5
  29 +nxpix = 200
  30 +nypix = 200
  31 +binsz = 0.02
  32 +coordsys = "CEL"
  33 +proj = "CAR"
  34 +
  35 +def run_sim(seed, model_name):
  36 + sim = ctobssim()
  37 + sim.logFileOpen() # We need this to explicitely open the log file in Python mode
  38 + sim["infile"].filename(model_name)
  39 + sim["outfile"].filename(events_name)
  40 + sim["caldb"].string(caldb)
  41 + sim["irf"].string(irf)
  42 + sim["seed"].integer(seed)
  43 + sim["ra"].real(ra)
  44 + sim["dec"].real(dec)
  45 + sim["rad"].real(rad_sim)
  46 + sim["tmin"].real(tstart)
  47 + sim["tmax"].real(tstop)
  48 + sim["emin"].real(emin)
  49 + sim["emax"].real(emax)
  50 + sim.execute()
  51 + return sim
  52 +
  53 +def run_bin():
  54 + bin = ctbin()
  55 + bin.logFileOpen() # We need this to explicitely open the log file in Python mode
  56 + bin["evfile"].filename(events_name)
  57 + bin["outfile"].filename(cntmap_name)
  58 + bin["emin"].real(emin)
  59 + bin["emax"].real(emax)
  60 + bin["enumbins"].integer(enumbins)
  61 + bin["nxpix"].integer(nxpix)
  62 + bin["nypix"].integer(nypix)
  63 + bin["binsz"].real(binsz)
  64 + bin["coordsys"].string(coordsys)
  65 + bin["xref"].real(ra)
  66 + bin["yref"].real(dec)
  67 + bin["proj"].string(proj)
  68 + bin.execute()
  69 + return bin
  70 +
  71 +def run_like(model_name):
  72 + like = ctlike()
  73 + like.logFileOpen() # We need this to explicitely open the log file in Python mode
  74 + like["cntmap"].filename(cntmap_name)
  75 + like["srcmdl"].filename(model_name)
  76 + like["outmdl"].filename(result_name)
  77 + like["method"].string("BINNED")
  78 + like["caldb"].string(caldb)
  79 + like["irf"].string(irf)
  80 + like.execute()
  81 + return like
  82 +
  83 +def test_coverage(n_trials=3, model_name="data/crab.xml"):
  84 + """Simulate and fit ntrials time and look at the
  85 + pull distribution for all free parameters.
  86 + It should be Gaussian with mean 0 and rms 1"""
  87 +
  88 + results = []
  89 + colnames = [] # to preserve the order!
  90 + for seed in range(n_trials):
  91 + print 'seed:', seed
  92 + sim = run_sim(seed, model_name)
  93 + run_bin()
  94 + like = run_like(model_name)
  95 + models = like.obs().models()
  96 + result = {}
  97 + colnames = []
  98 + for par in range(models.npars()):
  99 + p = models.par(par)
  100 + if p.isfree():
  101 + name = p.name()
  102 + colnames.append(name)
  103 + fitted_value = p.value()
  104 + real_value = sim.obs().models().par(par).value()
  105 + error = p.error()
  106 + pull = (fitted_value - real_value) / error
  107 + result[name] = pull
  108 + results.append(result)
  109 + return results, colnames
  110 +
  111 +if __name__ == '__main__':
  112 + usage = "test_coverage xml_model_file n_trials results_file"
  113 + if len(sys.argv) < 4:
  114 + print usage
  115 + sys.exit()
  116 +
  117 + xml_model_file = sys.argv[1]
  118 + n_trials = int(sys.argv[2])
  119 + results_file = sys.argv[3]
  120 +
  121 + print("Testing coverage for file {0} with {1} trials"
  122 + "".format(xml_model_file, n_trials))
  123 + results, colnames = test_coverage(n_trials, xml_model_file)
  124 +
  125 + print("Writing pull values for all free parameters to file {0}"
  126 + "".format(results_file))
  127 +
  128 + f = open(results_file, 'w')
  129 + writer = csv.DictWriter(f, colnames)
  130 + writer.writerow(dict((_,_) for _ in colnames))
  131 + writer.writerows(results)
  132 + f.close()
0 133 \ No newline at end of file
... ...