Maria Haupt / ctools

Commit 60cc0fd36bf91d78decf5bf072708452d449584d

Authored by knodlseder
1 parent af530c33

initial version

Inline Side-by-side
Showing 1 changed file with 649 additions and 0 deletions
scripts/cssens 0 → 100755
  View file @60cc0fd
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script computes the CTA sensitivity for a number of energy bins using
  4 +# ctlike. Crab spectra are fitted in narrow energy bins to simulated data,
  5 +# and the flux level is determined that leads to a particular significance.
  6 +# The significance can be either determined using the Test Statistic (which
  7 +# is defined as twice the likelihood difference between fitting with and
  8 +# without the test source) or using the source model scaling factor (or
  9 +# Prefactor) divided by its statistical error. The latter has been turned
  10 +# out to be more reliable, and is used as default parameter.
  11 +#
  12 +# Copyright (C) 2011 Jurgen Knodlseder
  13 +#
  14 +# This program is free software: you can redistribute it and/or modify
  15 +# it under the terms of the GNU General Public License as published by
  16 +# the Free Software Foundation, either version 3 of the License, or
  17 +# (at your option) any later version.
  18 +#
  19 +# This program is distributed in the hope that it will be useful,
  20 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  22 +# GNU General Public License for more details.
  23 +#
  24 +# You should have received a copy of the GNU General Public License
  25 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  26 +#
  27 +# ==========================================================================
  28 +from ctatools import *
  29 +from gammalib import *
  30 +import obsutils
  31 +import sys
  32 +import csv
  33 +import math
  34 +
  35 +
  36 +# ============ #
  37 +# cssens class #
  38 +# ============ #
  39 +class cssens(GApplication):
  40 + """
  41 + This class implements the sensitivity computation script. It derives
  42 + from the GammaLib::GApplication class which provides support for parameter
  43 + files, command line arguments, and logging. In that way the Python
  44 + script behaves just as a regular ctool.
  45 + """
  46 + def __init__(self):
  47 + """
  48 + Void constructor.
  49 + """
  50 + # Set name
  51 + self.name = "cssens"
  52 + self.version = "0.1.0"
  53 +
  54 + # Make sure that parfile exists
  55 + file = self.parfile()
  56 +
  57 + # Initialise application
  58 + GApplication.__init__(self, self.name, self.version)
  59 +
  60 + # Set logger properties
  61 + self.logFileOpen()
  62 + self.log_header()
  63 + self.log.date(True)
  64 +
  65 + # Return
  66 + return
  67 +
  68 + def __init__(self, argv):
  69 + """
  70 + Command line constructor.
  71 + """
  72 + # Set name
  73 + self.name = "cssens"
  74 + self.version = "0.1.0"
  75 +
  76 + # Make sure that parfile exists
  77 + file = self.parfile()
  78 +
  79 + # Initialise application
  80 + GApplication.__init__(self, self.name, self.version, argv)
  81 +
  82 + # Set logger properties
  83 + self.logFileOpen()
  84 + self.log_header()
  85 + self.log.date(True)
  86 +
  87 + # Return
  88 + return
  89 +
  90 + def __del__(self):
  91 + """
  92 + Destructor.
  93 + """
  94 + # Write separator into logger
  95 + if self.logTerse():
  96 + self.log("\n")
  97 +
  98 + # Return
  99 + return
  100 +
  101 + def parfile(self):
  102 + """
  103 + Check if parfile exists. If parfile does not exist then create a
  104 + default parfile. This kluge avoids shipping the cscript with a parfile.
  105 + """
  106 + # Set parfile name
  107 + parfile = self.name+".par"
  108 +
  109 + try:
  110 + pars = GPars(parfile)
  111 + except:
  112 + # Signal if parfile was not found
  113 + print "Parfile "+parfile+" not found. Create default parfile."
  114 +
  115 + # Create default parfile
  116 + pars = GPars()
  117 + pars.append(GPar("outfile","f","h","sensitivity.dat","","","Output file name"))
  118 + pars.append(GPar("caldb","s","h","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  119 + pars.append(GPar("irf","s","h","kb_E_50h_v3","","","Instrument response function"))
  120 + pars.append(GPar("duration","r","h","180000.0","","","Effective exposure time (s)"))
  121 + pars.append(GPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  122 + pars.append(GPar("enumbins","i","h","0","","","Number of energy bins"))
  123 + pars.append(GPar("sigma","r","h","5.0","","","Significance threshold"))
  124 + pars.append(GPar("ts_use","b","h","no","","","Use TS instead of Sigma?"))
  125 + pars.append(GPar("index","r","h","-2.48","","","Assumed spectral index"))
  126 + pars.append(GPar("max_iter","i","h","50","","","Maximum number of iterations"))
  127 + pars.append(GPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
  128 + pars.append_standard()
  129 + pars.save(parfile)
  130 +
  131 + # Return
  132 + return
  133 +
  134 + def get_parameters(self):
  135 + """
  136 + Get parameters from parfile.
  137 + """
  138 + # Get parameters
  139 + self.m_outfile = self["outfile"].filename()
  140 + self.m_caldb = self["caldb"].string()
  141 + self.m_irf = self["irf"].string()
  142 + self.m_duration = self["duration"].real()
  143 + self.m_roi = self["rad"].real()
  144 + self.m_nbins = self["enumbins"].integer()
  145 + self.m_ts_thres = self["sigma"].real()*self["sigma"].real()
  146 + self.m_use_ts = self["ts_use"].boolean()
  147 + self.m_index = self["index"].real()
  148 + self.m_max_iter = self["max_iter"].integer()
  149 + self.m_num_avg = self["num_avg"].integer()
  150 +
  151 + # Set some fixed parameters
  152 + self.m_log = False # Logging in client tools
  153 + self.m_debug = False # Debugging in client tools
  154 +
  155 + # Return
  156 + return
  157 +
  158 + def execute(self):
  159 + """
  160 + Execute the script.
  161 + """
  162 + # Run the script
  163 + self.run()
  164 +
  165 + # Return
  166 + return
  167 +
  168 + def run(self):
  169 + """
  170 + Run the script.
  171 + """
  172 + # Switch screen logging on in debug mode
  173 + if self.logDebug():
  174 + self.log.cout(True)
  175 +
  176 + # Get parameters
  177 + self.get_parameters()
  178 +
  179 + # Write input parameters into logger
  180 + if self.logTerse():
  181 + self.log_parameters()
  182 + self.log("\n")
  183 +
  184 + # Initialise script
  185 + colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux']
  186 + results = []
  187 +
  188 + # Initialise models. Note that we centre the point source at the Galactic
  189 + # center as our observation is also centred at the Galactic centre, so
  190 + # we're onaxis.
  191 + bkg_model = GModels()
  192 + full_model = GModels()
  193 + bkg_model.append(self.set_bkg_model())
  194 + full_model.append(self.set_bkg_model())
  195 + full_model.append(self.set_src_model(0.0, 0.0, index=self.m_index))
  196 +
  197 + # Write models into logger
  198 + if self.logTerse():
  199 + self.log("\n")
  200 + self.log.header1("Models")
  201 + self.log.header2("Background model")
  202 + self.log(str(bkg_model))
  203 + self.log("\n\n")
  204 + self.log.header2("Full model")
  205 + self.log(str(full_model))
  206 + self.log("\n")
  207 +
  208 + # Write heder
  209 + if self.logTerse():
  210 + self.log("\n")
  211 + self.log.header1("Sensitivity determination")
  212 +
  213 + # Loop over energy bands. The number of energy band is still
  214 + # fixed ...
  215 + for ieng in range(21):
  216 +
  217 + # Set energies
  218 + loge = -1.7 + ieng * 0.2
  219 + emean = pow(10.0, loge)
  220 + emin = pow(10.0, loge-0.1)
  221 + emax = pow(10.0, loge+0.1)
  222 +
  223 + # Setup observation(s)
  224 + obs = self.set_obs(emin=emin, emax=emax)
  225 +
  226 + # Determine sensitivity
  227 + result = self.get_sensitivity(obs, bkg_model, full_model)
  228 +
  229 + # Write results
  230 + if ieng == 0:
  231 + f = open(self.m_outfile, 'w')
  232 + writer = csv.DictWriter(f, colnames)
  233 + writer.writerow(dict((_,_) for _ in colnames))
  234 + writer.writerow(result)
  235 + f.close()
  236 + else:
  237 + f = open(self.m_outfile, 'a')
  238 + writer = csv.DictWriter(f, colnames)
  239 + writer.writerow(result)
  240 + f.close()
  241 +
  242 + # Append results
  243 + results.append(result)
  244 +
  245 + # Return
  246 + return
  247 +
  248 + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  249 + """
  250 + Returns an observation container with a single CTA observation.
  251 +
  252 + Keywords:
  253 + lpnt - Galactic longitude of pointing [deg] (default: 0.0)
  254 + bpnt - Galactic latitude of pointing [deg] (default: 0.0)
  255 + emin - Minimum energy [TeV] (default: 0.1)
  256 + emax - Maximum energy [TeV] (default: 100.0)
  257 + """
  258 + # Allocate observation container
  259 + obs = GObservations()
  260 +
  261 + # Set single pointing
  262 + pntdir = GSkyDir()
  263 + pntdir.lb_deg(lpnt, bpnt)
  264 +
  265 + # Create CTA observation
  266 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  267 + duration=self.m_duration, \
  268 + emin=emin, emax=emax, rad=self.m_roi)
  269 +
  270 + # Append observation to container
  271 + obs.append(run)
  272 +
  273 + # Return observation container
  274 + return obs
  275 +
  276 + def set_bkg_model(self, fitidx=False):
  277 + """
  278 + Setup standard CTA background model. We use a simple power law here,
  279 + scaled to Konrad's E configuration performance table.
  280 +
  281 + Keywords:
  282 + fitidx - Fit spectral index (default: False)
  283 + """
  284 + # Define radial component
  285 + radial = GCTAModelRadialGauss(3.0)
  286 +
  287 + # Define spectral component
  288 + spectrum = GModelSpectralPlaw(61.8, -1.85)
  289 + spectrum["Prefactor"].scale(1.0e-6)
  290 + spectrum["PivotEnergy"].value(1.0)
  291 + spectrum["PivotEnergy"].scale(1.0e6)
  292 + if fitidx:
  293 + spectrum["Index"].free()
  294 + else:
  295 + spectrum["Index"].fix()
  296 +
  297 + # Create background model
  298 + model = GCTAModelRadialAcceptance(radial, spectrum)
  299 + model.name("Background")
  300 + model.instruments("CTA")
  301 +
  302 + # Return background model
  303 + return model
  304 +
  305 + def set_src_model(self, l, b, flux=1.0, index=-2.48, \
  306 + type="point", sigma=1.0, radius=1.0, width=0.1, \
  307 + fitpos=False, fitidx=False):
  308 + """
  309 + Returns a single source with Crab-like spectrum. The source flux
  310 + can be scaled in Crab units. The Crab spectrum is based on MAGIC
  311 + observations (Albert et al. 2008, ApJ, 674, 1037).
  312 +
  313 + Parameters:
  314 + l - Galactic longitude of source location [deg]
  315 + b - Galactic latitude of source location [deg]
  316 + Keywords:
  317 + flux - Source flux [Crabs]
  318 + index - Spectral index
  319 + type - Source type ("point", "gauss", "disk", "shell")
  320 + sigma - Gaussian sigma (for type="gauss")
  321 + radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
  322 + width - Shell width [deg] (for type="shell")
  323 + fitpos - Fit position? (default: True)
  324 + fitidx - Fit index? (default: True)
  325 + """
  326 + # Set source location
  327 + location = GSkyDir()
  328 + location.lb_deg(l, b)
  329 +
  330 + # Set source spectrum
  331 + spectrum = GModelSpectralPlaw(flux, index)
  332 + spectrum["Prefactor"].scale(5.7e-16)
  333 + spectrum["PivotEnergy"].value(0.3)
  334 + spectrum["PivotEnergy"].scale(1.0e6)
  335 + if fitidx:
  336 + spectrum["Index"].free()
  337 + else:
  338 + spectrum["Index"].fix()
  339 +
  340 + # Set source
  341 + if type == "point":
  342 + spatial = GModelSpatialPtsrc(location)
  343 + if fitpos:
  344 + spatial[0].free()
  345 + spatial[1].free()
  346 + source = GModelPointSource(spatial, spectrum)
  347 + elif type == "gauss":
  348 + radial = GModelRadialGauss(location, sigma)
  349 + if fitpos:
  350 + radial[0].free()
  351 + radial[1].free()
  352 + source = GModelExtendedSource(radial, spectrum)
  353 + elif type == "disk":
  354 + radial = GModelRadialDisk(location, radius)
  355 + if fitpos:
  356 + radial[0].free()
  357 + radial[1].free()
  358 + source = GModelExtendedSource(radial, spectrum)
  359 + elif type == "shell":
  360 + radial = GModelRadialShell(location, radius, width)
  361 + if fitpos:
  362 + radial[0].free()
  363 + radial[1].free()
  364 + source = GModelExtendedSource(radial, spectrum)
  365 + else:
  366 + self.log("ERROR: Unknown source type '"+type+"'.\n")
  367 + return None
  368 +
  369 + # Set source name
  370 + source.name("Test")
  371 +
  372 + # Return source
  373 + return source
  374 +
  375 + def get_sensitivity(self, obs, bkg_model, full_model):
  376 + """
  377 + Determine sensitivity for a given observations.
  378 +
  379 + Parameters:
  380 + obs - Observation container
  381 + bkg_model - Background model
  382 + full_model - Source model
  383 + """
  384 + # Determine energy boundaries from first observation in
  385 + # the container
  386 + for run in obs:
  387 + emin = run.events().ebounds().emin()
  388 + emax = run.events().ebounds().emax()
  389 + loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
  390 + break
  391 +
  392 + # Write header
  393 + if self.logTerse():
  394 + self.log("\n")
  395 + self.log.header2("Energies: "+str(emin)+" - "+str(emax))
  396 +
  397 + # Initialise loop
  398 + flux_value = []
  399 + pflux_value = []
  400 + eflux_value = []
  401 + iter = 0
  402 + test_flux = 0.01 # This is the initial test flux in Crab units
  403 +
  404 + # Loop until we break
  405 + while True:
  406 +
  407 + # Update iteration counter
  408 + iter += 1
  409 +
  410 + # Write header
  411 + if self.logExplicit():
  412 + self.log.header2("Iteation "+str(iter))
  413 +
  414 + # Set test flux
  415 + full_model[1]['Prefactor'].value(test_flux)
  416 + obs.models(full_model)
  417 +
  418 + # Simulate events
  419 + sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
  420 + log=self.m_log, debug=self.m_debug)
  421 +
  422 + # Determine number of events in simulation
  423 + nevents = 0.0
  424 + for run in sim:
  425 + nevents += run.events().number()
  426 +
  427 + # Write simulation results
  428 + if self.logExplicit():
  429 + self.log.header3("Simulation")
  430 + self.log.parformat("Number of simulated events")
  431 + self.log(nevents)
  432 + self.log("\n")
  433 +
  434 + # Fit background only
  435 + sim.models(bkg_model)
  436 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  437 + result_bgm = like.obs().models().copy()
  438 + LogL_bgm = like.opt().value()
  439 + npred_bgm = like.obs().npred()
  440 +
  441 + # Store results
  442 + Pre_bgm = result_bgm[0]['Prefactor'].value()
  443 + ePre_bgm = result_bgm[0]['Prefactor'].error()
  444 + Inx_bgm = result_bgm[0]['Index'].value()
  445 + eInx_bgm = result_bgm[0]['Index'].error()
  446 + Sig_bgm = result_bgm[0]['Sigma'].value()
  447 + eSig_bgm = result_bgm[0]['Sigma'].error()
  448 +
  449 + # Write background fit results
  450 + if self.logExplicit():
  451 + self.log.header3("Background model fit")
  452 + self.log.parformat("log likelihood")
  453 + self.log(LogL_bgm)
  454 + self.log("\n")
  455 + self.log.parformat("Number of predicted events")
  456 + self.log(npred_bgm)
  457 + self.log("\n")
  458 + for model in result_bgm:
  459 + self.log.parformat("Model")
  460 + self.log(model.name())
  461 + self.log("\n")
  462 + for par in model:
  463 + self.log(str(par)+"\n")
  464 +
  465 + # Fit background and test source
  466 + sim.models(full_model)
  467 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  468 + result_all = like.obs().models().copy()
  469 + LogL_all = like.opt().value()
  470 + npred_all = like.obs().npred()
  471 + ts = 2.0*(LogL_bgm-LogL_all)
  472 +
  473 + # Get fitted Crab, photon and energy fluxes
  474 + cflux = result_all[0]['Prefactor'].value()
  475 + pflux = cast_GModelSky(result_all["Test"]).spectral().flux(emin, emax)
  476 + eflux = cast_GModelSky(result_all["Test"]).spectral().eflux(emin, emax)
  477 +
  478 + # Store results
  479 + Pre_all = result_all[0]['Prefactor'].value()
  480 + ePre_all = result_all[0]['Prefactor'].error()
  481 + Inx_all = result_all[0]['Index'].value()
  482 + eInx_all = result_all[0]['Index'].error()
  483 + Sig_all = result_all[0]['Sigma'].value()
  484 + eSig_all = result_all[0]['Sigma'].error()
  485 + SPre_all = result_all[1]['Prefactor'].value()
  486 + eSPre_all = result_all[1]['Prefactor'].error()
  487 + SInx_all = result_all[1]['Index'].value()
  488 + eSInx_all = result_all[1]['Index'].error()
  489 +
  490 + # Write background and test source fit results
  491 + if self.logExplicit():
  492 + self.log.header3("Background and test source model fit")
  493 + self.log.parformat("Test statistics")
  494 + self.log(ts)
  495 + self.log("\n")
  496 + self.log.parformat("log likelihood")
  497 + self.log(LogL_all)
  498 + self.log("\n")
  499 + self.log.parformat("Number of predicted events")
  500 + self.log(npred_all)
  501 + self.log("\n")
  502 + self.log.parformat("Photon flux")
  503 + self.log(pflux)
  504 + self.log(" ph/cm2/s\n")
  505 + self.log.parformat("Energy flux")
  506 + self.log(eflux)
  507 + self.log(" erg/cm2/s\n")
  508 + self.log.parformat("Crab flux")
  509 + self.log(cflux*1000.0)
  510 + self.log(" mCrab\n")
  511 + for model in result_all:
  512 + self.log.parformat("Model")
  513 + self.log(model.name())
  514 + self.log("\n")
  515 + for par in model:
  516 + self.log(str(par)+"\n")
  517 + elif self.logTerse():
  518 + self.log.parformat("Iteration "+str(iter))
  519 + self.log("TS=")
  520 + self.log(ts)
  521 + self.log(" ")
  522 + self.log(pflux)
  523 + self.log(" ph/cm2/s = ")
  524 + self.log(eflux)
  525 + self.log(" erg/cm2/s = ")
  526 + self.log(cflux*1000.0)
  527 + self.log(" mCrab\n")
  528 +
  529 + # Compute extrapolated fitted prefactor, photon and energy fluxes
  530 + correct = 0.0
  531 + if self.m_use_ts:
  532 + if ts > 0:
  533 + correct = math.sqrt(self.m_ts_thres/ts)
  534 + else:
  535 + if eSPre_all > 0:
  536 + correct = math.sqrt(self.m_ts_thres)/(SPre_all/eSPre_all)
  537 + flux = correct * SPre_all
  538 + pflux = correct * pflux
  539 + eflux = correct * eflux
  540 + flux_value.append(flux)
  541 + pflux_value.append(pflux)
  542 + eflux_value.append(eflux)
  543 +
  544 + # Compute sliding average of extrapolated fitted prefactor,
  545 + # photon and energy flux. This damps out fluctuations and
  546 + # improves convergence
  547 + flux = 0.0
  548 + pflux = 0.0
  549 + eflux = 0.0
  550 + num = 0.0
  551 + for k in range(self.m_num_avg):
  552 + inx = iter - k - 1
  553 + if inx >= 0:
  554 + flux += flux_value[inx]
  555 + pflux += pflux_value[inx]
  556 + eflux += eflux_value[inx]
  557 + num += 1.0
  558 + flux /= num
  559 + pflux /= num
  560 + eflux /= num
  561 +
  562 + # Compare average flux to last average
  563 + if iter > self.m_num_avg:
  564 + if test_flux > 0:
  565 + ratio = flux/test_flux
  566 + if ratio >= 0.99 and ratio <= 1.01:
  567 + if self.logTerse():
  568 + self.log(" Converged ("+str(ratio)+")\n")
  569 + break
  570 + else:
  571 + if self.logTerse():
  572 + self.log(" Flux is zero.\n")
  573 + break
  574 +
  575 + # Use average for next iteration
  576 + test_flux = flux
  577 +
  578 + # Exit loop if number of trials exhausted
  579 + if (iter >= self.m_max_iter):
  580 + if self.logTerse():
  581 + self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n")
  582 + break
  583 +
  584 + # Write fit results
  585 + if self.logTerse():
  586 + self.log.header3("Fit results")
  587 + self.log.parformat("Test statistics")
  588 + self.log(ts)
  589 + self.log("\n")
  590 + self.log.parformat("Photon flux")
  591 + self.log(pflux)
  592 + self.log(" ph/cm2/s\n")
  593 + self.log.parformat("Energy flux")
  594 + self.log(eflux)
  595 + self.log(" erg/cm2/s\n")
  596 + self.log.parformat("Crab flux")
  597 + self.log(cflux*1000.0)
  598 + self.log(" mCrab\n")
  599 + self.log.parformat("Number of simulated events")
  600 + self.log(nevents)
  601 + self.log("\n")
  602 + self.log.header3("Background and test source model fitting")
  603 + self.log.parformat("log likelihood")
  604 + self.log(LogL_all)
  605 + self.log("\n")
  606 + self.log.parformat("Number of predicted events")
  607 + self.log(npred_all)
  608 + self.log("\n")
  609 + for model in result_all:
  610 + self.log.parformat("Model")
  611 + self.log(model.name())
  612 + self.log("\n")
  613 + for par in model:
  614 + self.log(str(par)+"\n")
  615 + self.log.header3("Background model fit")
  616 + self.log.parformat("log likelihood")
  617 + self.log(LogL_bgm)
  618 + self.log("\n")
  619 + self.log.parformat("Number of predicted events")
  620 + self.log(npred_bgm)
  621 + self.log("\n")
  622 + for model in result_bgm:
  623 + self.log.parformat("Model")
  624 + self.log(model.name())
  625 + self.log("\n")
  626 + for par in model:
  627 + self.log(str(par)+"\n")
  628 +
  629 + # Store result
  630 + result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
  631 + 'crab': flux, 'pflux': pflux, 'eflux': eflux}
  632 +
  633 + # Return result
  634 + return result
  635 +
  636 +
  637 +# ======================== #
  638 +# Main routine entry point #
  639 +# ======================== #
  640 +if __name__ == '__main__':
  641 + """
  642 + CTA ctlike-based sensitivity estimator.
  643 + """
  644 + # Create instance of application
  645 + app = cssens(sys.argv)
  646 +
  647 + # Run application
  648 + app.execute()
  649 +
0 650 \ No newline at end of file
... ...