new_cssens.py
| 1 |
#! /usr/bin/env python
|
|---|---|
| 2 |
# ==========================================================================
|
| 3 |
# This script computes the CTA sensitivity.
|
| 4 |
#
|
| 5 |
# Copyright (C) 2011-2016 Juergen Knoedlseder
|
| 6 |
#
|
| 7 |
# This program is free software: you can redistribute it and/or modify
|
| 8 |
# it under the terms of the GNU General Public License as published by
|
| 9 |
# the Free Software Foundation, either version 3 of the License, or
|
| 10 |
# (at your option) any later version.
|
| 11 |
#
|
| 12 |
# This program is distributed in the hope that it will be useful,
|
| 13 |
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
| 14 |
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
| 15 |
# GNU General Public License for more details.
|
| 16 |
#
|
| 17 |
# You should have received a copy of the GNU General Public License
|
| 18 |
# along with this program. If not, see <http://www.gnu.org/licenses/>.
|
| 19 |
#
|
| 20 |
# ==========================================================================
|
| 21 |
import sys |
| 22 |
import csv |
| 23 |
import math |
| 24 |
import gammalib |
| 25 |
import ctools |
| 26 |
from cscripts import obsutils |
| 27 |
|
| 28 |
|
| 29 |
# ============ #
|
| 30 |
# cssens class #
|
| 31 |
# ============ #
|
| 32 |
class cssens(ctools.cscript): |
| 33 |
"""
|
| 34 |
Computes the CTA sensitivity.
|
| 35 |
|
| 36 |
This class computes the CTA sensitivity for a number of energy bins using
|
| 37 |
ctlike. Spectra are fitted in narrow energy bins to simulated data,
|
| 38 |
and the flux level is determined that leads to a particular significance.
|
| 39 |
|
| 40 |
The significance is determined using the Test Statistic, defined as twice
|
| 41 |
the likelihood difference between fitting with and without the test source.
|
| 42 |
"""
|
| 43 |
|
| 44 |
# Constructor
|
| 45 |
def __init__(self, *argv): |
| 46 |
"""
|
| 47 |
Constructor.
|
| 48 |
"""
|
| 49 |
# Set name
|
| 50 |
self._name = "cssens" |
| 51 |
self._version = "1.1.0" |
| 52 |
|
| 53 |
# Initialise class members
|
| 54 |
self._obs = gammalib.GObservations()
|
| 55 |
self._ebounds = gammalib.GEbounds()
|
| 56 |
self._obs_ebounds = []
|
| 57 |
self._srcname = "" |
| 58 |
self._outfile = gammalib.GFilename()
|
| 59 |
self._ra = None |
| 60 |
self._dec = None |
| 61 |
self._edisp = False |
| 62 |
self._emin = 0.020 |
| 63 |
self._emax = 200.0 |
| 64 |
self._bins = 21 |
| 65 |
self._enumbins = 0 |
| 66 |
self._npix = 200 |
| 67 |
self._binsz = 0.05 |
| 68 |
self._type = "Differential" |
| 69 |
self._ts_thres = 25.0 |
| 70 |
self._max_iter = 50 |
| 71 |
self._num_avg = 3 |
| 72 |
self._log_clients = False |
| 73 |
self._debug = False |
| 74 |
|
| 75 |
# Initialise application by calling the appropriate class
|
| 76 |
# constructor.
|
| 77 |
self._init_cscript(argv)
|
| 78 |
|
| 79 |
# Return
|
| 80 |
return
|
| 81 |
|
| 82 |
|
| 83 |
# Private methods
|
| 84 |
def _get_parameters(self): |
| 85 |
"""
|
| 86 |
Get user parameters from parfile.
|
| 87 |
"""
|
| 88 |
# Set observation if not done before
|
| 89 |
if self._obs == None or self._obs.size() == 0: |
| 90 |
self._obs = self._set_obs() |
| 91 |
|
| 92 |
# Set models if we have none
|
| 93 |
if self._obs.models().size() == 0: |
| 94 |
self._obs.models(self["inmodel"].filename()) |
| 95 |
|
| 96 |
# Get source name
|
| 97 |
self._srcname = self["srcname"].string() |
| 98 |
|
| 99 |
# Read further parameters
|
| 100 |
self._outfile = self["outfile"].filename() |
| 101 |
self._emin = self["emin"].real() |
| 102 |
self._emax = self["emax"].real() |
| 103 |
self._bins = self["bins"].integer() |
| 104 |
|
| 105 |
# Read parameters for binned if requested
|
| 106 |
self._enumbins = self["enumbins"].integer() |
| 107 |
if not self._enumbins == 0: |
| 108 |
self._npix = self["npix"].integer() |
| 109 |
self._binsz = self["binsz"].real() |
| 110 |
|
| 111 |
# Read remaining parameters
|
| 112 |
self._edisp = self["edisp"].boolean() |
| 113 |
self._ts_thres = self["sigma"].real()*self["sigma"].real() |
| 114 |
self._max_iter = self["max_iter"].integer() |
| 115 |
self._num_avg = self["num_avg"].integer() |
| 116 |
self._type = self["type"].string() |
| 117 |
|
| 118 |
# Set some fixed parameters
|
| 119 |
self._debug = self["debug"].boolean() # Debugging in client tools |
| 120 |
|
| 121 |
# Derive some parameters
|
| 122 |
self._ebounds = gammalib.GEbounds(self._bins, |
| 123 |
gammalib.GEnergy(self._emin, "TeV"), |
| 124 |
gammalib.GEnergy(self._emax, "TeV")) |
| 125 |
|
| 126 |
# Write input parameters into logger
|
| 127 |
if self._logTerse(): |
| 128 |
self._log_parameters()
|
| 129 |
self._log("\n") |
| 130 |
|
| 131 |
# Return
|
| 132 |
return
|
| 133 |
|
| 134 |
def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): |
| 135 |
"""
|
| 136 |
Set an observation container.
|
| 137 |
|
| 138 |
Kwargs:
|
| 139 |
lpnt: Galactic longitude of pointing [deg] (default: 0.0)
|
| 140 |
bpnt: Galactic latitude of pointing [deg] (default: 0.0)
|
| 141 |
emin: Minimum energy [TeV] (default: 0.1)
|
| 142 |
emax: Maximum energy [TeV] (default: 100.0)
|
| 143 |
|
| 144 |
Returns:
|
| 145 |
Observation container.
|
| 146 |
"""
|
| 147 |
# If an observation was provided on input then load it from XML
|
| 148 |
# file
|
| 149 |
filename = self["inobs"].filename() |
| 150 |
if filename != "NONE" and filename != "": |
| 151 |
obs = self._get_observations()
|
| 152 |
|
| 153 |
# ... otherwise allocate a single observation
|
| 154 |
else:
|
| 155 |
|
| 156 |
# Read relevant user parameters
|
| 157 |
caldb = self["caldb"].string() |
| 158 |
irf = self["irf"].string() |
| 159 |
deadc = self["deadc"].real() |
| 160 |
duration = self["duration"].real() |
| 161 |
rad = self["rad"].real() |
| 162 |
|
| 163 |
# Allocate observation container
|
| 164 |
obs = gammalib.GObservations() |
| 165 |
|
| 166 |
# Set single pointing
|
| 167 |
pntdir = gammalib.GSkyDir() |
| 168 |
pntdir.lb_deg(lpnt, bpnt) |
| 169 |
|
| 170 |
# Create CTA observation
|
| 171 |
run = obsutils.set_obs(pntdir, caldb=caldb, irf=irf, |
| 172 |
duration=duration, deadc=deadc, |
| 173 |
emin=emin, emax=emax, rad=rad) |
| 174 |
|
| 175 |
# Append observation to container
|
| 176 |
obs.append(run) |
| 177 |
|
| 178 |
# Set source position
|
| 179 |
offset = self["offset"].real() |
| 180 |
pntdir.lb_deg(lpnt, bpnt+offset) |
| 181 |
self._ra = pntdir.ra_deg()
|
| 182 |
self._dec = pntdir.dec_deg()
|
| 183 |
|
| 184 |
# Return observation container
|
| 185 |
return obs
|
| 186 |
|
| 187 |
def _set_models(self, fitspat=False, fitspec=False): |
| 188 |
"""
|
| 189 |
Set full model and background model.
|
| 190 |
|
| 191 |
Kwargs:
|
| 192 |
fitspat: Fit spatial parameter (default: False).
|
| 193 |
fitspec: Fit spectral parameters (default: False).
|
| 194 |
|
| 195 |
Returns:
|
| 196 |
Tuple containing full model and background model.
|
| 197 |
"""
|
| 198 |
# Retrieve full model from observation container
|
| 199 |
full_model = self._obs.models().copy()
|
| 200 |
|
| 201 |
# Get source model
|
| 202 |
model = full_model[self._srcname]
|
| 203 |
|
| 204 |
# Check that model has a Prefactor
|
| 205 |
if not model.has_par("Prefactor"): |
| 206 |
msg = "Model \""+self._srcname+"\" has no parameter "+\ |
| 207 |
"\"Prefactor\". Only spectral models with a "+\
|
| 208 |
"\"Prefactor\" parameter are supported."
|
| 209 |
raise RuntimeError(msg) |
| 210 |
|
| 211 |
# Set source position
|
| 212 |
if self._ra != None and self._dec != None: |
| 213 |
if model.has_par("RA") and model.has_par("DEC"): |
| 214 |
model["RA"].value(self._ra) |
| 215 |
model["DEC"].value(self._dec) |
| 216 |
|
| 217 |
# Fit or fix spatial parameters
|
| 218 |
if fitspat:
|
| 219 |
if model.has_par("RA"): |
| 220 |
model["RA"].free()
|
| 221 |
if model.has_par("DEC"): |
| 222 |
model["DEC"].free()
|
| 223 |
if model.has_par("Sigma"): |
| 224 |
model["Sigma"].free()
|
| 225 |
if model.has_par("Radius"): |
| 226 |
model["Radius"].free()
|
| 227 |
if model.has_par("Width"): |
| 228 |
model["Width"].free()
|
| 229 |
else:
|
| 230 |
if model.has_par("RA"): |
| 231 |
model["RA"].fix()
|
| 232 |
if model.has_par("DEC"): |
| 233 |
model["DEC"].fix()
|
| 234 |
if model.has_par("Sigma"): |
| 235 |
model["Sigma"].fix()
|
| 236 |
if model.has_par("Radius"): |
| 237 |
model["Radius"].fix()
|
| 238 |
if model.has_par("Width"): |
| 239 |
model["Width"].fix()
|
| 240 |
|
| 241 |
# Fit or fix spectral parameters
|
| 242 |
if fitspec:
|
| 243 |
if model.has_par("Index"): |
| 244 |
model["Index"].free()
|
| 245 |
if model.has_par("Cutoff"): |
| 246 |
model["Cutoff"].free()
|
| 247 |
else:
|
| 248 |
if model.has_par("Index"): |
| 249 |
model["Index"].fix()
|
| 250 |
if model.has_par("Cutoff"): |
| 251 |
model["Cutoff"].fix()
|
| 252 |
|
| 253 |
# Create background model
|
| 254 |
bkg_model = full_model.copy() |
| 255 |
bkg_model.remove(self._srcname)
|
| 256 |
|
| 257 |
# Return models
|
| 258 |
return full_model, bkg_model
|
| 259 |
|
| 260 |
def _set_obs_ebounds(self, emin, emax): |
| 261 |
"""
|
| 262 |
Set energy boundaries for observation container.
|
| 263 |
|
| 264 |
Sets the energy boundaries for all observations in the observation
|
| 265 |
container.
|
| 266 |
|
| 267 |
Args:
|
| 268 |
emin: Minimum energy
|
| 269 |
emax: Maximum energy
|
| 270 |
"""
|
| 271 |
# Loop over all observations in container
|
| 272 |
for obs in self._obs: |
| 273 |
|
| 274 |
# Get observation energy boundaries
|
| 275 |
obs_ebounds = obs.events().ebounds() |
| 276 |
|
| 277 |
# Get minimum and maximum energy
|
| 278 |
obs_emin = obs_ebounds.emin() |
| 279 |
obs_emax = obs_ebounds.emax() |
| 280 |
|
| 281 |
# Case A: bin fully contained in observation ebounds
|
| 282 |
if obs_emin <= emin and obs_emax >= emax: |
| 283 |
ebounds = gammalib.GEbounds(emin, emax) |
| 284 |
|
| 285 |
# Case B: bin fully outside of obs ebounds
|
| 286 |
elif emax < obs_emin or emin > obs_emax: |
| 287 |
|
| 288 |
# Set zero range (inspired by ctselect)
|
| 289 |
e0 = gammalib.GEnergy(0.0, "TeV") |
| 290 |
ebounds = gammalib.GEbounds(e0, e0) |
| 291 |
|
| 292 |
# Case C: bin partly overlapping with observation ebounds
|
| 293 |
else:
|
| 294 |
|
| 295 |
# Set energy range as obs ebounds were fully contained inside energy bin
|
| 296 |
set_emin = emin |
| 297 |
set_emax = emax |
| 298 |
|
| 299 |
# Adjust energy bin to respect observation energy boundary
|
| 300 |
if emin < obs_emin:
|
| 301 |
set_emin = obs_emin |
| 302 |
if emax > obs_emax:
|
| 303 |
set_emax = obs_emax |
| 304 |
|
| 305 |
#Set energy boundaries
|
| 306 |
ebounds = gammalib.GEbounds(set_emin, set_emax) |
| 307 |
|
| 308 |
# Set energy boundaries
|
| 309 |
obs.events().ebounds(ebounds) |
| 310 |
|
| 311 |
# Return
|
| 312 |
return
|
| 313 |
|
| 314 |
def _get_crab_flux(self, emin, emax): |
| 315 |
"""
|
| 316 |
Return Crab photon flux in a given energy interval.
|
| 317 |
|
| 318 |
Args:
|
| 319 |
emin: Minimum energy
|
| 320 |
emax: Maximum energy
|
| 321 |
|
| 322 |
Returns:
|
| 323 |
Crab photon flux in specified energy interval.
|
| 324 |
"""
|
| 325 |
# Set Crab TeV spectral model based on a power law
|
| 326 |
crab = gammalib.GModelSpectralPlaw(5.7e-16, -2.48, |
| 327 |
gammalib.GEnergy(0.3, "TeV")) |
| 328 |
|
| 329 |
# Determine flux
|
| 330 |
flux = crab.flux(emin, emax) |
| 331 |
|
| 332 |
# Return flux
|
| 333 |
return flux
|
| 334 |
|
| 335 |
def _get_sensitivity(self, emin, emax, bkg_model, full_model): |
| 336 |
"""
|
| 337 |
Determine sensitivity for given observations.
|
| 338 |
|
| 339 |
Args:
|
| 340 |
emin: Minimum energy for fitting and flux computation
|
| 341 |
emax: Maximum energy for fitting and flux computation
|
| 342 |
bkg_model: Background model
|
| 343 |
full_model: Source model
|
| 344 |
|
| 345 |
Returns:
|
| 346 |
Result dictionary
|
| 347 |
"""
|
| 348 |
# Set TeV->erg conversion factor
|
| 349 |
tev2erg = 1.6021764
|
| 350 |
|
| 351 |
# Set energy boundaries
|
| 352 |
self._set_obs_ebounds(emin, emax)
|
| 353 |
|
| 354 |
# Determine energy boundaries from first observation in the container
|
| 355 |
loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) |
| 356 |
e_mean = math.pow(10.0, loge)
|
| 357 |
erg_mean = e_mean * tev2erg |
| 358 |
|
| 359 |
# Compute Crab unit (this is the factor with which the Prefactor needs
|
| 360 |
# to be multiplied to get 1 Crab
|
| 361 |
crab_flux = self._get_crab_flux(emin, emax)
|
| 362 |
src_flux = full_model[self._srcname].spectral().flux(emin, emax)
|
| 363 |
crab_unit = crab_flux/src_flux |
| 364 |
|
| 365 |
# Write header
|
| 366 |
if self._logTerse(): |
| 367 |
self._log("\n") |
| 368 |
self._log.header2("Energies: "+str(emin)+" - "+str(emax)) |
| 369 |
self._log.parformat("Crab flux") |
| 370 |
self._log(crab_flux)
|
| 371 |
self._log(" ph/cm2/s\n") |
| 372 |
self._log.parformat("Source model flux") |
| 373 |
self._log(src_flux)
|
| 374 |
self._log(" ph/cm2/s\n") |
| 375 |
self._log.parformat("Crab unit factor") |
| 376 |
self._log(crab_unit)
|
| 377 |
self._log("\n") |
| 378 |
|
| 379 |
# Initialise loop
|
| 380 |
crab_flux_value = [] |
| 381 |
photon_flux_value = [] |
| 382 |
energy_flux_value = [] |
| 383 |
sensitivity_value = [] |
| 384 |
iterations = 0
|
| 385 |
test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab) |
| 386 |
|
| 387 |
# Loop until we break
|
| 388 |
while True: |
| 389 |
|
| 390 |
# Update iteration counter
|
| 391 |
iterations += 1
|
| 392 |
|
| 393 |
# Write header
|
| 394 |
if self._logExplicit(): |
| 395 |
self._log.header2("Iteration "+str(iterations)) |
| 396 |
|
| 397 |
# Set source model. crab_prefactor is the Prefactor that
|
| 398 |
# corresponds to 1 Crab
|
| 399 |
src_model = full_model.copy() |
| 400 |
crab_prefactor = src_model[self._srcname]['Prefactor'].value() * \ |
| 401 |
crab_unit |
| 402 |
src_model[self._srcname]['Prefactor'].value(crab_prefactor * \ |
| 403 |
test_crab_flux) |
| 404 |
self._obs.models(src_model)
|
| 405 |
|
| 406 |
# Simulate events
|
| 407 |
sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations, |
| 408 |
binsz=self._binsz, npix=self._npix, |
| 409 |
log=self._log_clients, debug=self._debug, |
| 410 |
edisp=self._edisp)
|
| 411 |
|
| 412 |
# Determine number of events in simulation
|
| 413 |
nevents = 0.0
|
| 414 |
for run in sim: |
| 415 |
nevents += run.events().number() |
| 416 |
|
| 417 |
# Write simulation results
|
| 418 |
if self._logExplicit(): |
| 419 |
self._log.header3("Simulation") |
| 420 |
self._log.parformat("Number of simulated events") |
| 421 |
self._log(nevents)
|
| 422 |
self._log("\n") |
| 423 |
|
| 424 |
# Fit background only
|
| 425 |
sim.models(bkg_model) |
| 426 |
like = obsutils.fit(sim, log=self._log_clients,
|
| 427 |
debug=self._debug, edisp=self._edisp) |
| 428 |
result_bgm = like.obs().models().copy() |
| 429 |
LogL_bgm = like.opt().value() |
| 430 |
npred_bgm = like.obs().npred() |
| 431 |
|
| 432 |
# Assess quality based on a comparison between Npred and Nevents
|
| 433 |
quality_bgm = npred_bgm-nevents |
| 434 |
|
| 435 |
# Write background fit results
|
| 436 |
if self._logExplicit(): |
| 437 |
self._log.header3("Background model fit") |
| 438 |
self._log.parformat("log likelihood") |
| 439 |
self._log(LogL_bgm)
|
| 440 |
self._log("\n") |
| 441 |
self._log.parformat("Number of predicted events") |
| 442 |
self._log(npred_bgm)
|
| 443 |
self._log("\n") |
| 444 |
self._log.parformat("Fit quality") |
| 445 |
self._log(quality_bgm)
|
| 446 |
self._log("\n") |
| 447 |
|
| 448 |
# Write model fit results
|
| 449 |
if self._logExplicit(): |
| 450 |
for model in result_bgm: |
| 451 |
self._log.parformat("Model") |
| 452 |
self._log(model.name())
|
| 453 |
self._log("\n") |
| 454 |
for par in model: |
| 455 |
self._log(str(par)+"\n") |
| 456 |
|
| 457 |
# Fit background and test source
|
| 458 |
sim.models(src_model) |
| 459 |
like = obsutils.fit(sim, log=self._log_clients,
|
| 460 |
debug=self._debug, edisp=self._edisp) |
| 461 |
result_all = like.obs().models().copy() |
| 462 |
LogL_all = like.opt().value() |
| 463 |
npred_all = like.obs().npred() |
| 464 |
ts = 2.0*(LogL_bgm-LogL_all)
|
| 465 |
|
| 466 |
# Assess quality based on a comparison between Npred and Nevents
|
| 467 |
quality_all = npred_all-nevents |
| 468 |
|
| 469 |
# Write background and test source fit results
|
| 470 |
if self._logExplicit(): |
| 471 |
self._log.header3("Background and test source model fit") |
| 472 |
self._log.parformat("Test statistics") |
| 473 |
self._log(ts)
|
| 474 |
self._log("\n") |
| 475 |
self._log.parformat("log likelihood") |
| 476 |
self._log(LogL_all)
|
| 477 |
self._log("\n") |
| 478 |
self._log.parformat("Number of predicted events") |
| 479 |
self._log(npred_all)
|
| 480 |
self._log("\n") |
| 481 |
self._log.parformat("Fit quality") |
| 482 |
self._log(quality_all)
|
| 483 |
self._log("\n") |
| 484 |
#
|
| 485 |
for model in result_all: |
| 486 |
self._log.parformat("Model") |
| 487 |
self._log(model.name())
|
| 488 |
self._log("\n") |
| 489 |
for par in model: |
| 490 |
self._log(str(par)+"\n") |
| 491 |
|
| 492 |
# Start over if TS was non-positive
|
| 493 |
if ts <= 0.0: |
| 494 |
if self._logExplicit(): |
| 495 |
self._log("Non positive TS. Start over.\n") |
| 496 |
continue
|
| 497 |
|
| 498 |
# Get fitted Crab, photon and energy fluxes
|
| 499 |
crab_flux = result_all[self._srcname]['Prefactor'].value() / \ |
| 500 |
crab_prefactor |
| 501 |
#crab_flux_err = result_all[self._srcname]['Prefactor'].error() / \
|
| 502 |
# crab_prefactor
|
| 503 |
photon_flux = result_all[self._srcname].spectral().flux(emin, emax)
|
| 504 |
energy_flux = result_all[self._srcname].spectral().eflux(emin, emax)
|
| 505 |
|
| 506 |
# Compute differential sensitivity in unit erg/cm2/s
|
| 507 |
energy = gammalib.GEnergy(e_mean, "TeV")
|
| 508 |
time = gammalib.GTime() |
| 509 |
sensitivity = result_all[self._srcname].spectral().eval(energy, time) * \
|
| 510 |
e_mean*erg_mean * 1.0e6
|
| 511 |
|
| 512 |
# Compute flux correction factor based on average TS
|
| 513 |
correct = 1.0
|
| 514 |
if ts > 0: |
| 515 |
correct = math.sqrt(self._ts_thres/ts)
|
| 516 |
|
| 517 |
# Compute extrapolated fluxes
|
| 518 |
crab_flux = correct * crab_flux |
| 519 |
photon_flux = correct * photon_flux |
| 520 |
energy_flux = correct * energy_flux |
| 521 |
sensitivity = correct * sensitivity |
| 522 |
crab_flux_value.append(crab_flux) |
| 523 |
photon_flux_value.append(photon_flux) |
| 524 |
energy_flux_value.append(energy_flux) |
| 525 |
sensitivity_value.append(sensitivity) |
| 526 |
|
| 527 |
# Write background and test source fit results
|
| 528 |
if self._logExplicit(): |
| 529 |
self._log.parformat("Photon flux") |
| 530 |
self._log(photon_flux)
|
| 531 |
self._log(" ph/cm2/s\n") |
| 532 |
self._log.parformat("Energy flux") |
| 533 |
self._log(energy_flux)
|
| 534 |
self._log(" erg/cm2/s\n") |
| 535 |
self._log.parformat("Crab flux") |
| 536 |
self._log(crab_flux*1000.0) |
| 537 |
self._log(" mCrab\n") |
| 538 |
self._log.parformat("Differential sensitivity") |
| 539 |
self._log(sensitivity)
|
| 540 |
self._log(" erg/cm2/s\n") |
| 541 |
for model in result_all: |
| 542 |
self._log.parformat("Model") |
| 543 |
self._log(model.name())
|
| 544 |
self._log("\n") |
| 545 |
for par in model: |
| 546 |
self._log(str(par)+"\n") |
| 547 |
elif self._logTerse(): |
| 548 |
self._log.parformat("Iteration "+str(iterations)) |
| 549 |
self._log("TS=") |
| 550 |
self._log(ts)
|
| 551 |
self._log(" ") |
| 552 |
self._log("corr=") |
| 553 |
self._log(correct)
|
| 554 |
self._log(" ") |
| 555 |
self._log(photon_flux)
|
| 556 |
self._log(" ph/cm2/s = ") |
| 557 |
self._log(energy_flux)
|
| 558 |
self._log(" erg/cm2/s = ") |
| 559 |
self._log(crab_flux*1000.0) |
| 560 |
self._log(" mCrab = ") |
| 561 |
self._log(sensitivity)
|
| 562 |
self._log(" erg/cm2/s\n") |
| 563 |
|
| 564 |
# Compute sliding average of extrapolated fitted prefactor,
|
| 565 |
# photon and energy flux. This damps out fluctuations and
|
| 566 |
# improves convergence
|
| 567 |
crab_flux = 0.0
|
| 568 |
photon_flux = 0.0
|
| 569 |
energy_flux = 0.0
|
| 570 |
sensitivity = 0.0
|
| 571 |
num = 0.0
|
| 572 |
for k in range(self._num_avg): |
| 573 |
inx = len(crab_flux_value) - k - 1 |
| 574 |
if inx >= 0: |
| 575 |
crab_flux += crab_flux_value[inx] |
| 576 |
photon_flux += photon_flux_value[inx] |
| 577 |
energy_flux += energy_flux_value[inx] |
| 578 |
sensitivity += sensitivity_value[inx] |
| 579 |
num += 1.0
|
| 580 |
crab_flux /= num |
| 581 |
photon_flux /= num |
| 582 |
energy_flux /= num |
| 583 |
sensitivity /= num |
| 584 |
|
| 585 |
# Compare average flux to last average
|
| 586 |
if iterations > self._num_avg: |
| 587 |
if test_crab_flux > 0: |
| 588 |
ratio = crab_flux/test_crab_flux |
| 589 |
|
| 590 |
# We have 2 convergence criteria:
|
| 591 |
# 1. The average flux does not change
|
| 592 |
# 2. The flux correction factor is small
|
| 593 |
if ratio >= 0.99 and ratio <= 1.01 and \ |
| 594 |
correct >= 0.9 and correct <= 1.1: |
| 595 |
if self._logTerse(): |
| 596 |
self._log(" Converged ("+str(ratio)+")\n") |
| 597 |
break
|
| 598 |
else:
|
| 599 |
if self._logTerse(): |
| 600 |
self._log(" Flux is zero.\n") |
| 601 |
break
|
| 602 |
|
| 603 |
# Use average for next iteration
|
| 604 |
test_crab_flux = crab_flux |
| 605 |
|
| 606 |
# Exit loop if number of trials exhausted
|
| 607 |
if (iterations >= self._max_iter): |
| 608 |
if self._logTerse(): |
| 609 |
self._log(" Test ended after "+str(self._max_iter)+ |
| 610 |
" iterations.\n")
|
| 611 |
break
|
| 612 |
|
| 613 |
# Write fit results
|
| 614 |
if self._logTerse(): |
| 615 |
self._log.header3("Fit results") |
| 616 |
self._log.parformat("Test statistics") |
| 617 |
self._log(ts)
|
| 618 |
self._log("\n") |
| 619 |
self._log.parformat("Photon flux") |
| 620 |
self._log(photon_flux)
|
| 621 |
self._log(" ph/cm2/s\n") |
| 622 |
self._log.parformat("Energy flux") |
| 623 |
self._log(energy_flux)
|
| 624 |
self._log(" erg/cm2/s\n") |
| 625 |
self._log.parformat("Crab flux") |
| 626 |
self._log(crab_flux*1000.0) |
| 627 |
self._log(" mCrab\n") |
| 628 |
self._log.parformat("Differential sensitivity") |
| 629 |
self._log(sensitivity)
|
| 630 |
self._log(" erg/cm2/s\n") |
| 631 |
self._log.parformat("Number of simulated events") |
| 632 |
self._log(nevents)
|
| 633 |
self._log("\n") |
| 634 |
self._log.header3("Background and test source model fitting") |
| 635 |
self._log.parformat("log likelihood") |
| 636 |
self._log(LogL_all)
|
| 637 |
self._log("\n") |
| 638 |
self._log.parformat("Number of predicted events") |
| 639 |
self._log(npred_all)
|
| 640 |
self._log("\n") |
| 641 |
for model in result_all: |
| 642 |
self._log.parformat("Model") |
| 643 |
self._log(model.name())
|
| 644 |
self._log("\n") |
| 645 |
for par in model: |
| 646 |
self._log(str(par)+"\n") |
| 647 |
self._log.header3("Background model fit") |
| 648 |
self._log.parformat("log likelihood") |
| 649 |
self._log(LogL_bgm)
|
| 650 |
self._log("\n") |
| 651 |
self._log.parformat("Number of predicted events") |
| 652 |
self._log(npred_bgm)
|
| 653 |
self._log("\n") |
| 654 |
for model in result_bgm: |
| 655 |
self._log.parformat("Model") |
| 656 |
self._log(model.name())
|
| 657 |
self._log("\n") |
| 658 |
for par in model: |
| 659 |
self._log(str(par)+"\n") |
| 660 |
|
| 661 |
# Restore energy boundaries of observation container
|
| 662 |
for i, obs in enumerate(self._obs): |
| 663 |
obs.events().ebounds(self._obs_ebounds[i])
|
| 664 |
|
| 665 |
# Store result
|
| 666 |
result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
|
| 667 |
'crab_flux': crab_flux, 'photon_flux': photon_flux, \ |
| 668 |
'energy_flux': energy_flux, \
|
| 669 |
'sensitivity': sensitivity}
|
| 670 |
|
| 671 |
# Return result
|
| 672 |
return result
|
| 673 |
|
| 674 |
|
| 675 |
# Public methods
|
| 676 |
def run(self): |
| 677 |
"""
|
| 678 |
Run the script.
|
| 679 |
"""
|
| 680 |
# Switch screen logging on in debug mode
|
| 681 |
if self._logDebug(): |
| 682 |
self._log.cout(True) |
| 683 |
|
| 684 |
# Get parameters
|
| 685 |
self._get_parameters()
|
| 686 |
|
| 687 |
# Loop over observations and store ebounds
|
| 688 |
for obs in self._obs: |
| 689 |
self._obs_ebounds.append(obs.events().ebounds())
|
| 690 |
|
| 691 |
# Initialise script
|
| 692 |
colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', |
| 693 |
'energy_flux', 'sensitivity'] |
| 694 |
results = [] |
| 695 |
|
| 696 |
# Initialise models
|
| 697 |
full_model, bkg_model = self._set_models()
|
| 698 |
|
| 699 |
# Write models into logger
|
| 700 |
if self._logTerse(): |
| 701 |
self._log("\n") |
| 702 |
self._log.header1("Models") |
| 703 |
self._log.header2("Background model") |
| 704 |
self._log(str(bkg_model)) |
| 705 |
self._log("\n\n") |
| 706 |
self._log.header2("Full model") |
| 707 |
self._log(str(full_model)) |
| 708 |
self._log("\n") |
| 709 |
|
| 710 |
# Write header
|
| 711 |
if self._logTerse(): |
| 712 |
self._log("\n") |
| 713 |
self._log.header1("Sensitivity determination") |
| 714 |
self._log.parformat("Type") |
| 715 |
self._log(self._type) |
| 716 |
self._log("\n") |
| 717 |
|
| 718 |
# Loop over energy bins
|
| 719 |
for ieng in range(self._ebounds.size()): |
| 720 |
|
| 721 |
# Set energies
|
| 722 |
if self._type == "Differential": |
| 723 |
emin = self._ebounds.emin(ieng)
|
| 724 |
emax = self._ebounds.emax(ieng)
|
| 725 |
elif self._type == "Integral": |
| 726 |
emin = self._ebounds.emin(ieng)
|
| 727 |
emax = self._ebounds.emax()
|
| 728 |
else:
|
| 729 |
msg = "Invalid sensitivity type \""+self._type+"\" "+\ |
| 730 |
"encountered. Either specify \"Differential\" or "+\
|
| 731 |
"\"Integral\"."
|
| 732 |
raise RuntimeError(msg)
|
| 733 |
|
| 734 |
# Determine sensitivity
|
| 735 |
result = self._get_sensitivity(emin, emax,
|
| 736 |
bkg_model, full_model) |
| 737 |
|
| 738 |
# Write results
|
| 739 |
if ieng == 0: |
| 740 |
f = open(self._outfile.url(), 'w') |
| 741 |
writer = csv.DictWriter(f, colnames) |
| 742 |
headers = {}
|
| 743 |
for n in colnames: |
| 744 |
headers[n] = n |
| 745 |
writer.writerow(headers) |
| 746 |
writer.writerow(result) |
| 747 |
f.close() |
| 748 |
else:
|
| 749 |
f = open(self._outfile.url(), 'a') |
| 750 |
writer = csv.DictWriter(f, colnames) |
| 751 |
writer.writerow(result) |
| 752 |
f.close() |
| 753 |
|
| 754 |
# Append results
|
| 755 |
results.append(result) |
| 756 |
|
| 757 |
# Return
|
| 758 |
return
|
| 759 |
|
| 760 |
def execute(self): |
| 761 |
"""
|
| 762 |
Execute the script.
|
| 763 |
"""
|
| 764 |
# Open logfile
|
| 765 |
self.logFileOpen()
|
| 766 |
|
| 767 |
# Run the script
|
| 768 |
self.run()
|
| 769 |
|
| 770 |
# Return
|
| 771 |
return
|
| 772 |
|
| 773 |
|
| 774 |
|
| 775 |
# ======================== #
|
| 776 |
# Main routine entry point #
|
| 777 |
# ======================== #
|
| 778 |
if __name__ == '__main__': |
| 779 |
|
| 780 |
# Create instance of application
|
| 781 |
app = cssens(sys.argv) |
| 782 |
|
| 783 |
# Run application
|
| 784 |
app.execute() |