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