Stefan Eschbach / ctools

Commit b9778a813a83396098bf54cc7bbb25d3083c34f6

Authored by Jürgen Knödlseder
1 parent e77d88fa

Replace averaging by regression method in cssens (#1989)

Inline Side-by-side
Showing 6 changed files with 223 additions and 175 deletions
ChangeLog
  View file @b9778a8
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.3.0 released
4 4 ========================
5 5  
  6 + Implement regression scheme in cssens and remove num_avg parameter (#1989)
6 7 Copy cscripts into bin directory instead of creating a symbolic links (#1987)
7 8 If code originates from Git, always compile the Python bindings (#1985)
8 9 Raise exception instead of exiting with return code in Python tests (#1983)
... ...
... ... @@ -133,6 +133,10 @@ extension that can be used for a stacked analysis. The energy boundaries
133 133 are set as function of the variation in the effective area and background
134 134 rates of the Instrument Response Function (#1976).
135 135  
  136 +A regression scheme now replaces the averaging scheme in cssens as this
  137 +allows more accurately the prediction of the flux value at the threshold
  138 +Test Statistic. The num_avg parameter has been removed (#1989).
  139 +
136 140  
137 141 Examples
138 142 --------
... ...
cscripts/cssens.par
  View file @b9778a8
... ... @@ -20,7 +20,7 @@ inobs, f, h, NONE,,, "Input event list, counts cube, or observation definition
20 20 inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
21 21 srcname, s, a, Crab,,, "Source name"
22 22 caldb, s, a, prod2,,, "Calibration database"
23   -irf, s, a, South_0.5h,,, "Instrument response function"
  23 +irf, s, a, South_50h,,, "Instrument response function"
24 24 edisp, b, h, no,,, "Apply energy dispersion?"
25 25 deadc, r, h, 0.98,0,1, "Deadtime correction factor"
26 26 outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file"
... ... @@ -31,16 +31,15 @@ outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file"
31 31 offset, r, h, 0.0,,, "Source offset angle (deg)"
32 32 duration, r, a, 180000.0,,, "Effective exposure time (s)"
33 33 rad, r, a, 5.0,,, "Radius of ROI (deg)"
34   -emin, r, a, 0.020,,, "Lower energy limit (TeV)"
35   -emax, r, a, 200.0,,, "Upper energy limit (TeV)"
36   -bins, i, a, 21,,, "Number of energy bins for differential sensitivity computation"
  34 +emin, r, a, 0.030,,, "Lower energy limit (TeV)"
  35 +emax, r, a, 180.0,,, "Upper energy limit (TeV)"
  36 +bins, i, a, 20,,, "Number of energy bins for differential sensitivity computation"
37 37 enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)"
38 38 npix, i, h, 200,,, "Number of pixels for binned analysis"
39 39 binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)"
40 40 type, s, h, Differential,Differential|Integral,, "Sensitivity type"
41 41 sigma, r, h, 5.0,,, "Significance threshold (Gaussian sigma)"
42 42 max_iter, i, h, 50,,, "Maximum number of iterations"
43   -num_avg, i, h, 3,,, "Number of iterations for sliding average"
44 43  
45 44 #
46 45 # Standard parameters
... ... @@ -49,6 +48,6 @@ chatter, i, h, 2,0,4, "Chattiness of output"
49 48 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
50 49 debug, b, h, no,,, "Debugging mode activated"
51 50 mode, s, h, ql,,, "Mode of automatic parameters"
52   -logfile, f, h, cspull.log,,, "Log filename"
  51 +logfile, f, h, cssens.log,,, "Log filename"
53 52  
54 53  
... ...
cscripts/cssens.py
  View file @b9778a8
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# Computes the array sensitivity using a test source
  3 +# Computes the array sensitivity using the Test Statistic for a test source
4 4 #
5   -# Copyright (C) 2011-2016 Juergen Knoedlseder
  5 +# Copyright (C) 2011-2017 Juergen Knoedlseder
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -40,7 +40,7 @@ class cssens(ctools.cscript):
40 40 and the flux level is determined that leads to a particular significance.
41 41  
42 42 The significance is determined using the Test Statistic, defined as twice
43   - the likelihood difference between fitting with and without the test source.
  43 + the likelihood difference between fitting with and without the test source.
44 44 """
45 45  
46 46 # Constructor
... ... @@ -50,31 +50,20 @@ class cssens(ctools.cscript):
50 50 """
51 51 # Set name
52 52 self._name = 'cssens'
53   - self._version = '1.2.0'
  53 + self._version = '1.3.0'
54 54  
55 55 # Initialise class members
56   - self._obs = gammalib.GObservations()
57 56 self._ebounds = gammalib.GEbounds()
58 57 self._obs_ebounds = []
59 58 self._srcname = ''
60   - self._outfile = gammalib.GFilename()
61 59 self._ra = None
62 60 self._dec = None
63   - self._edisp = False
64   - self._emin = 0.020
65   - self._emax = 200.0
66   - self._bins = 21
67   - self._enumbins = 0
68   - self._npix = 200
69   - self._binsz = 0.05
70   - self._type = 'Differential'
71   - self._ts_thres = 25.0
72   - self._max_iter = 50
73   - self._num_avg = 3
74 61 self._log_clients = False
75   -
76   - # Initialise application by calling the appropriate base class
77   - # constructor
  62 +
  63 + # Initialise observation container from constructor arguments
  64 + self._obs, argv = self._set_input_obs(argv)
  65 +
  66 + # Initialise application by calling the appropriate class constructor
78 67 self._init_cscript(argv)
79 68  
80 69 # Return
... ... @@ -87,8 +76,8 @@ class cssens(ctools.cscript):
87 76 Get user parameters from parfile
88 77 """
89 78 # Set observation if not done before
90   - if self._obs == None or self._obs.size() == 0:
91   - self._obs = self._set_obs()
  79 + if self._obs.size() == 0:
  80 + self._obs = self._set_obs(self['emin'].real(), self['emax'].real())
92 81  
93 82 # Set models if we have none
94 83 if self._obs.models().size() == 0:
... ... @@ -98,31 +87,28 @@ class cssens(ctools.cscript):
98 87 self._srcname = self['srcname'].string()
99 88  
100 89 # Read further parameters
101   - self._outfile = self['outfile'].filename()
102   - self._emin = self['emin'].real()
103   - self._emax = self['emax'].real()
104   - self._bins = self['bins'].integer()
105   -
106   - # Read parameters for binned if requested
107   - self._enumbins = self['enumbins'].integer()
108   - if not self._enumbins == 0:
109   - self._npix = self['npix'].integer()
110   - self._binsz = self['binsz'].real()
111   -
112   - # Read remaining parameters
113   - self._edisp = self['edisp'].boolean()
114   - self._ts_thres = self['sigma'].real() * self['sigma'].real()
115   - self._max_iter = self['max_iter'].integer()
116   - self._num_avg = self['num_avg'].integer()
117   - self._type = self['type'].string()
118   -
119   - # Query remaining parameters
  90 + emin = self['emin'].real()
  91 + emax = self['emax'].real()
  92 + bins = self['bins'].integer()
  93 +
  94 + # Query parameters for binned if requested
  95 + enumbins = self['enumbins'].integer()
  96 + if not enumbins == 0:
  97 + self['npix'].integer()
  98 + self['binsz'].real()
  99 +
  100 + # Query input parameters
  101 + self['sigma'].real()
  102 + self['max_iter'].integer()
  103 + self['type'].string()
  104 + self['outfile'].filename()
  105 + self['edisp'].boolean()
120 106 self['debug'].boolean()
121 107  
122 108 # Derive some parameters
123   - self._ebounds = gammalib.GEbounds(self._bins,
124   - gammalib.GEnergy(self._emin, 'TeV'),
125   - gammalib.GEnergy(self._emax, 'TeV'))
  109 + self._ebounds = gammalib.GEbounds(bins,
  110 + gammalib.GEnergy(emin, 'TeV'),
  111 + gammalib.GEnergy(emax, 'TeV'))
126 112  
127 113 # Write input parameters into logger
128 114 self._log_parameters(gammalib.TERSE)
... ... @@ -130,20 +116,20 @@ class cssens(ctools.cscript):
130 116 # Return
131 117 return
132 118  
133   - def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  119 + def _set_obs(self, emin, emax, lpnt=0.0, bpnt=0.0):
134 120 """
135 121 Set an observation container
136 122  
137 123 Parameters
138 124 ----------
  125 + emin : float
  126 + Minimum energy (TeV)
  127 + emax : float
  128 + Maximum energy (TeV)
139 129 lpnt : float, optional
140 130 Galactic longitude of pointing (deg)
141 131 bpnt : float, optional
142 132 Galactic latitude of pointing (deg)
143   - emin : float, optional
144   - Minimum energy (TeV)
145   - emax : float, optional
146   - Maximum energy (TeV)
147 133  
148 134 Returns
149 135 -------
... ... @@ -205,23 +191,23 @@ class cssens(ctools.cscript):
205 191  
206 192 # Get energy boundaries of the observation
207 193 obs_ebounds = self._obs_ebounds[i]
208   -
  194 +
209 195 # Get minimum and maximum energy of the observation
210 196 obs_emin = obs_ebounds.emin()
211 197 obs_emax = obs_ebounds.emax()
212   -
  198 +
213 199 # If [emin,emax] is fully contained in the observation energy range
214 200 # the use [emin,emax] as energy boundaries
215 201 if obs_emin <= emin and obs_emax >= emax:
216 202 ebounds = gammalib.GEbounds(emin, emax)
217   -
  203 +
218 204 # ... otherwise, if [emin,emax] is completely outside the
219 205 # observation energy range then set the energy boundaries to the
220 206 # zero-width interval [0,0]
221 207 elif emax < obs_emin or emin > obs_emax:
222 208 e0 = gammalib.GEnergy(0.0, 'TeV')
223 209 ebounds = gammalib.GEbounds(e0, e0)
224   -
  210 +
225 211 # ... otherwise, if [emin,emax] overlaps partially with the
226 212 # observation energy range then set the energy boundaries to the
227 213 # overlapping part
... ... @@ -229,10 +215,10 @@ class cssens(ctools.cscript):
229 215 # Set overlapping energy range
230 216 set_emin = max(emin, obs_emin)
231 217 set_emax = min(emax, obs_emax)
232   -
  218 +
233 219 # Set energy boundaries
234 220 ebounds = gammalib.GEbounds(set_emin, set_emax)
235   -
  221 +
236 222 # Set the energy boundaries as the boundaries of the observation
237 223 obs.events().ebounds(ebounds)
238 224  
... ... @@ -286,6 +272,20 @@ class cssens(ctools.cscript):
286 272 # Set TeV->erg conversion factor
287 273 tev2erg = 1.6021764
288 274  
  275 + # Set parameters
  276 + ts_thres = self['sigma'].real() * self['sigma'].real()
  277 + max_iter = self['max_iter'].integer()
  278 + enumbins = self['enumbins'].integer()
  279 + if not enumbins == 0:
  280 + npix = self['npix'].integer()
  281 + binsz = self['binsz'].real()
  282 + else:
  283 + npix = 200
  284 + binsz = 0.05
  285 +
  286 + # Set flux ratio precision required for convergence to 5%
  287 + ratio_precision = 0.05
  288 +
289 289 # Set energy boundaries
290 290 self._set_obs_ebounds(emin, emax)
291 291  
... ... @@ -300,20 +300,27 @@ class cssens(ctools.cscript):
300 300 src_flux = test_model[self._srcname].spectral().flux(emin, emax)
301 301 crab_unit = crab_flux/src_flux
302 302  
303   - # Write header and some introductory information into logger
  303 + # Initialise regression coefficient
  304 + regcoeff = 0.0
  305 +
  306 + # Write header for energy bin
  307 + self._log_string(gammalib.TERSE, '')
304 308 self._log_header2(gammalib.TERSE, 'Energies: '+str(emin)+' - '+str(emax))
  309 +
  310 + # Write initial parameters
  311 + self._log_header3(gammalib.TERSE, 'Initial parameters')
305 312 self._log_value(gammalib.TERSE, 'Crab flux', str(crab_flux)+' ph/cm2/s')
306   - self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+
307   - ' ph/cm2/s')
  313 + self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+' ph/cm2/s')
308 314 self._log_value(gammalib.TERSE, 'Crab unit factor', crab_unit)
309 315  
310 316 # Initialise loop
311   - crab_flux_value = []
312   - photon_flux_value = []
313   - energy_flux_value = []
314   - sensitivity_value = []
315   - iterations = 0
316   - test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab)
  317 + results = []
  318 + iterations = 0
  319 + test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab)
  320 +
  321 + # Write header for iterations for terse chatter level
  322 + if self._logTerse():
  323 + self._log_header3(gammalib.TERSE, 'Iterations')
317 324  
318 325 # Loop until we break
319 326 while True:
... ... @@ -335,11 +342,11 @@ class cssens(ctools.cscript):
335 342  
336 343 # Simulate events for the models. "sim" holds an observation
337 344 # container with observations containing the simulated events.
338   - sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations,
339   - binsz=self._binsz, npix=self._npix,
  345 + sim = obsutils.sim(self._obs, nbins=enumbins, seed=iterations,
  346 + binsz=binsz, npix=npix,
340 347 log=self._log_clients,
341 348 debug=self['debug'].boolean(),
342   - edisp=self._edisp)
  349 + edisp=self['edisp'].boolean())
343 350  
344 351 # Determine number of events in simulation by summing the events
345 352 # over all observations in the observation container
... ... @@ -354,7 +361,7 @@ class cssens(ctools.cscript):
354 361 # Fit test source to the simulated events in the observation
355 362 # container
356 363 fit = ctools.ctlike(sim)
357   - fit['edisp'] = edisp=self._edisp
  364 + fit['edisp'] = self['edisp'].boolean()
358 365 fit['debug'] = self['debug'].boolean()
359 366 fit['chatter'] = self['chatter'].integer()
360 367 fit.run()
... ... @@ -383,29 +390,24 @@ class cssens(ctools.cscript):
383 390 # Assess quality based on a comparison between Npred and Nevents
384 391 quality = npred - nevents
385 392  
386   - # Write test source fit results into logger
387   - self._log_header3(gammalib.EXPLICIT, 'Test source model fit')
388   - self._log_value(gammalib.EXPLICIT, 'Test statistics', ts)
389   - self._log_value(gammalib.EXPLICIT, 'log likelihood', logL)
390   - self._log_value(gammalib.EXPLICIT, 'Number of predicted events', quality)
391   - self._log_value(gammalib.EXPLICIT, 'Fit quality', npred)
392   - for model in models:
393   - self._log_value(gammalib.EXPLICIT, 'Model', model.name())
394   - for par in model:
395   - self._log_string(gammalib.EXPLICIT, str(par))
396   -
397   - # If TS was non-positive then increase the test flux and start
398   - # over
  393 + # Write fit results into logger
  394 + name = 'Iteration %d' % iterations
  395 + value = ('TS=%10.4f Sim=%9.4f mCrab Fit=%9.4f mCrab '
  396 + 'Sens=%e erg/cm2/s' %
  397 + (ts, test_crab_flux*1000.0, crab_flux*1000.0, sensitivity))
  398 + self._log_value(gammalib.TERSE, name, value)
  399 +
  400 + # If TS was non-positive then increase the test flux and start over
399 401 if ts <= 0.0:
400 402  
401 403 # If the number of iterations was exceeded then stop
402   - if (iterations >= self._max_iter):
  404 + if (iterations >= max_iter):
403 405 self._log_string(gammalib.TERSE,
404   - ' Test ended after %d iterations.' % self._max_iter)
  406 + ' Test ended after %d iterations.' % max_iter)
405 407 break
406 408  
407   - # Increase test flux by a factor of 2
408   - test_crab_flux = test_crab_flux * 2.0
  409 + # Increase test flux
  410 + test_crab_flux *= 3.0
409 411  
410 412 # Signal start we start over
411 413 self._log_string(gammalib.EXPLICIT,
... ... @@ -414,95 +416,64 @@ class cssens(ctools.cscript):
414 416 # ... and start over
415 417 continue
416 418  
417   - # Compute flux correction factor based on average TS. We
418   - # assume here that the significance is proportional to the source
419   - # flux, and that the significance scales as the square root of
420   - # the Test Statistic
421   - correct = math.sqrt(self._ts_thres/ts)
  419 + # Append result entry to result list
  420 + result = {'ts': ts, 'crab_flux': crab_flux,
  421 + 'photon_flux': photon_flux,
  422 + 'energy_flux': energy_flux}
  423 + results.append(result)
  424 +
  425 + # Predict Crab flux at threshold TS using a linear regression of
  426 + # the log(TS) and log(crab_flux) values that have so far been
  427 + # computed. If not enough results are available than use a simple
  428 + # TS scaling relation.
  429 + if len(results) > 1:
  430 + pred_crab_flux, regcoeff = self._predict_flux(results, ts_thres)
  431 + correct = pred_crab_flux / crab_flux
  432 + else:
  433 + correct = math.sqrt(ts_thres/ts)
422 434  
423 435 # Compute extrapolated fluxes based on the flux correction factor
424 436 crab_flux = correct * crab_flux
425 437 photon_flux = correct * photon_flux
426 438 energy_flux = correct * energy_flux
427 439 sensitivity = correct * sensitivity
428   - crab_flux_value.append(crab_flux)
429   - photon_flux_value.append(photon_flux)
430   - energy_flux_value.append(energy_flux)
431   - sensitivity_value.append(sensitivity)
432   -
433   - # Write fit results
434   - if self._logExplicit():
435   - self._log_value(gammalib.EXPLICIT, 'Photon flux',
436   - str(photon_flux)+' ph/cm2/s')
437   - self._log_value(gammalib.EXPLICIT, 'Energy flux',
438   - str(energy_flux)+' erg/cm2/s')
439   - self._log_value(gammalib.EXPLICIT, 'Crab flux',
440   - str(crab_flux*1000.0)+' mCrab')
441   - self._log_value(gammalib.EXPLICIT, 'Differential sensitivity',
442   - str(sensitivity)+' erg/cm2/s')
443   - for model in models:
444   - self._log_value(gammalib.EXPLICIT, 'Model', model.name())
445   - for par in model:
446   - self._log_string(gammalib.EXPLICIT, str(par))
447   - elif self._logTerse():
448   - name = 'Iteration %d' % iterations
449   - value = ('TS=%.4f corr=%e %e ph/cm2/s = %e erg/cm2/s = %.2f '
450   - 'mCrab = %e erg/cm2/s' % (ts, correct, photon_flux,
451   - energy_flux, crab_flux*1000.0, sensitivity))
452   - self._log_value(gammalib.TERSE, name, value)
453   -
454   - # Compute sliding average of extrapolated fitted prefactor,
455   - # photon and energy flux. This damps out fluctuations and
456   - # improves convergence.
457   - crab_flux = 0.0
458   - photon_flux = 0.0
459   - energy_flux = 0.0
460   - sensitivity = 0.0
461   - num = 0.0
462   - for k in range(self._num_avg):
463   - inx = len(crab_flux_value) - k - 1
464   - if inx >= 0:
465   - crab_flux += crab_flux_value[inx]
466   - photon_flux += photon_flux_value[inx]
467   - energy_flux += energy_flux_value[inx]
468   - sensitivity += sensitivity_value[inx]
469   - num += 1.0
470   - crab_flux /= num
471   - photon_flux /= num
472   - energy_flux /= num
473   - sensitivity /= num
474   -
475   - # Compare average flux to last average
476   - if iterations > self._num_avg:
  440 +
  441 + # If we have at least 3 results then check if the flux determination
  442 + # at the TS threshold has converged
  443 + if len(results) > 3:
477 444 if test_crab_flux > 0:
  445 +
  446 + # Compute fractional change in the Crab flux between two
  447 + # iterations
478 448 ratio = crab_flux/test_crab_flux
479 449  
480   - # We have 2 convergence criteria:
481   - # 1. The average flux does not change
482   - # 2. The flux correction factor is small
483   - if ratio >= 0.99 and ratio <= 1.01 and \
484   - correct >= 0.9 and correct <= 1.1:
485   - if self._logTerse():
486   - self._log(' Converged ('+str(ratio)+')\n')
  450 + # If fractional change is smaller than the required position
  451 + # the iterations are stopped
  452 + if ratio > 1.0-ratio_precision and \
  453 + ratio < 1.0+ratio_precision:
  454 + value = ('TS=%10.4f Sim=%9.4f mCrab '
  455 + ' Sens=%e erg/cm2/s' %
  456 + (ts, crab_flux*1000.0, sensitivity))
  457 + self._log_value(gammalib.TERSE, 'Converged result', value)
  458 + self._log_value(gammalib.TERSE, 'Converged flux ratio', ratio)
  459 + self._log_value(gammalib.TERSE, 'Regression coefficient',
  460 + regcoeff)
487 461 break
488 462 else:
489   - if self._logTerse():
490   - self._log(' Flux is zero.\n')
  463 + self._log_value(gammalib.TERSE, 'Not converged', 'Flux is zero')
491 464 break
492 465  
493   - # Use average for next iteration
  466 + # Set test flux for next iteration
494 467 test_crab_flux = crab_flux
495 468  
496 469 # Exit loop if number of trials exhausted
497   - if (iterations >= self._max_iter):
498   - if self._logTerse():
499   - self._log(' Test ended after '+str(self._max_iter)+
500   - ' iterations.\n')
  470 + if (iterations >= max_iter):
  471 + self._log_string(gammalib.TERSE,
  472 + ' Test ended after %d iterations.' % max_iter)
501 473 break
502 474  
503 475 # Write fit results into logger
504 476 self._log_header3(gammalib.TERSE, 'Fit results')
505   - self._log_value(gammalib.TERSE, 'Test statistics', ts)
506 477 self._log_value(gammalib.TERSE, 'Photon flux',
507 478 str(photon_flux)+' ph/cm2/s')
508 479 self._log_value(gammalib.TERSE, 'Energy flux',
... ... @@ -523,16 +494,84 @@ class cssens(ctools.cscript):
523 494 # Restore energy boundaries of observation container
524 495 for i, obs in enumerate(self._obs):
525 496 obs.events().ebounds(self._obs_ebounds[i])
526   -
  497 +
527 498 # Store result
528 499 result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
529 500 'crab_flux': crab_flux, 'photon_flux': photon_flux, \
530 501 'energy_flux': energy_flux, \
531   - 'sensitivity': sensitivity}
  502 + 'sensitivity': sensitivity, 'regcoeff': regcoeff, \
  503 + 'nevents': nevents, 'npred': npred}
532 504  
533 505 # Return result
534 506 return result
535 507  
  508 + def _predict_flux(self, results, ts):
  509 + """
  510 + Predict Crab flux for a given TS value
  511 +
  512 + The Crab flux at a given Test Statistic value is predicted by doing a
  513 + linear regression of the log(TS) and log(crab_flux) values in a results
  514 + list.
  515 +
  516 + See https://en.wikipedia.org/wiki/Simple_linear_regression
  517 +
  518 + Parameters
  519 + ----------
  520 + results : list of dict
  521 + List of results
  522 + ts : float
  523 + Test Statistic value
  524 +
  525 + Returns
  526 + -------
  527 + crab_flux_prediction, regcoeff : tuple of float
  528 + Predicted Crab flux and regression coefficient
  529 + """
  530 + # Compute means and regression coefficient
  531 + mean_x = 0.0
  532 + mean_y = 0.0
  533 + mean_xy = 0.0
  534 + mean_xx = 0.0
  535 + mean_yy = 0.0
  536 + for result in results:
  537 + x = math.log(float(result['ts']))
  538 + y = math.log(float(result['crab_flux']))
  539 + mean_x += x
  540 + mean_y += y
  541 + mean_xy += x * y
  542 + mean_xx += x * x
  543 + mean_yy += y * y
  544 + norm = 1.0 / float(len(results))
  545 + mean_x *= norm
  546 + mean_y *= norm
  547 + mean_xy *= norm
  548 + mean_xx *= norm
  549 + mean_yy *= norm
  550 + rxy = (mean_xy - mean_x * mean_y) / \
  551 + math.sqrt((mean_xx - mean_x*mean_x) *
  552 + (mean_yy - mean_y*mean_y))
  553 + regcoeff = rxy*rxy
  554 +
  555 + # Compute regression line slope
  556 + beta_nom = 0.0
  557 + beta_denom = 0.0
  558 + for result in results:
  559 + x = math.log(float(result['ts']))
  560 + y = math.log(float(result['crab_flux']))
  561 + beta_nom += (x - mean_x) * (y - mean_y)
  562 + beta_denom += (x - mean_x) * (x - mean_x)
  563 + beta = beta_nom / beta_denom
  564 +
  565 + # Compute regression line offset
  566 + alpha = mean_y - beta * mean_x
  567 +
  568 + # Predict Crab flux at TS threshold
  569 + log_ts_thres = math.log(ts)
  570 + crab_flux_prediction = math.exp(alpha + beta * log_ts_thres)
  571 +
  572 + # Return
  573 + return (crab_flux_prediction, regcoeff)
  574 +
536 575  
537 576 # Public methods
538 577 def run(self):
... ... @@ -545,15 +584,16 @@ class cssens(ctools.cscript):
545 584  
546 585 # Get parameters
547 586 self._get_parameters()
548   -
  587 +
549 588 # Loop over observations and store a deep copy of the energy
550 589 # boundaries for later use
551 590 for obs in self._obs:
552 591 self._obs_ebounds.append(obs.events().ebounds().copy())
553   -
  592 +
554 593 # Initialise script
555 594 colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux',
556   - 'energy_flux', 'sensitivity']
  595 + 'energy_flux', 'sensitivity', 'regcoeff', 'nevents',
  596 + 'npred']
557 597 results = []
558 598  
559 599 # Set test source model for this observation
... ... @@ -566,30 +606,35 @@ class cssens(ctools.cscript):
566 606 # Write models into logger
567 607 self._log_models(gammalib.NORMAL, models, 'Model')
568 608  
  609 + # Get sensitivity type
  610 + sensitivity_type = self['type'].string()
  611 +
569 612 # Write header
570 613 self._log_header1(gammalib.TERSE, 'Sensitivity determination')
571   - self._log_value(gammalib.TERSE, 'Type', str(self._type))
  614 + self._log_value(gammalib.TERSE, 'Type', sensitivity_type)
572 615  
573 616 # Loop over energy bins
574 617 for ieng in range(self._ebounds.size()):
575 618  
576 619 # Set energies
577   - if self._type == 'Differential':
  620 + if sensitivity_type == 'Differential':
578 621 emin = self._ebounds.emin(ieng)
579 622 emax = self._ebounds.emax(ieng)
580   - elif self._type == 'Integral':
  623 + elif sensitivity_type == 'Integral':
581 624 emin = self._ebounds.emin(ieng)
582 625 emax = self._ebounds.emax()
583 626 else:
584 627 msg = ('Invalid sensitivity type "%s" encountered. Either '
585   - 'specify "Differential" or "Integral".' % self._type)
  628 + 'specify "Differential" or "Integral".' %
  629 + sensitivity_type)
586 630 raise RuntimeError(msg)
587 631  
588 632 # Determine sensitivity
589 633 result = self._get_sensitivity(emin, emax, models)
590 634  
591 635 # Write out trial result
592   - ioutils.write_csv_row(self._outfile.url(), ieng, colnames, result)
  636 + ioutils.write_csv_row(self['outfile'].filename().url(), ieng,
  637 + colnames, result)
593 638  
594 639 # Append results
595 640 results.append(result)
... ...
doc/source/users/reference_manual/cssens.rst
  View file @b9778a8
... ... @@ -100,9 +100,6 @@ General parameters
100 100 ``(max_iter = 50) [integer]``
101 101 Maximum number of iterations.
102 102  
103   -``(num_avg = 3) [integer]``
104   - Number of iterations for sliding average.
105   -
106 103  
107 104 Standard parameters
108 105 -------------------
... ...
src/ctobssim/ctobssim.cpp
  View file @b9778a8
... ... @@ -1471,7 +1471,8 @@ void ctobssim::simulate_background(GCTAObservation* obs,
1471 1471  
1472 1472 // Debug code: signal that we step into the response MC method
1473 1473 #if defined(G_BACKGROUND_DEBUG)
1474   - std::cout << "ctobssim::simulate_background: model->mc in" << std::endl;
  1474 + std::cout << "ctobssim::simulate_background: model->mc in"
  1475 + << std::endl;
1475 1476 #endif
1476 1477  
1477 1478 // Get simulated CTA event list. Note that this method
... ... @@ -1481,7 +1482,8 @@ void ctobssim::simulate_background(GCTAObservation* obs,
1481 1482  
1482 1483 // Debug code: signal that we came back
1483 1484 #if defined(G_BACKGROUND_DEBUG)
1484   - std::cout << "ctobssim::simulate_background: model->mc out" << std::endl;
  1485 + std::cout << "ctobssim::simulate_background: model->mc out"
  1486 + << std::endl;
1485 1487 #endif
1486 1488  
1487 1489 // Continue only if we got a CTA event list
... ...