Iftach Sadeh / ctools

Commit 955940f1ca82363890ca1052dea055ff17021612

Authored by Jürgen Knödlseder
1 parent e8f4b286

Refactor code, change names of new parameters and adapt code to coding conventions (#3668)

Inline Side-by-side
Showing 2 changed files with 261 additions and 192 deletions
cscripts/cssens.par
  View file @955940f
@@ -28,23 +28,23 @@ outfile, f, a, sensitivity.fits,,, "Output sensitivity FITS file" @@ -28,23 +28,23 @@ outfile, f, a, sensitivity.fits,,, "Output sensitivity FITS file"
28 # 28 #
29 # Script parameters 29 # Script parameters
30 #================== 30 #==================
31 -offset, r, h, 0.0,,, "Source offset angle (deg)"  
32 -duration, r, a, 180000.0,,, "Effective exposure time (s)"  
33 -rad, r, a, 5.0,,, "Radius of ROI (deg)"  
34 -emin, r, a, 0.03163,,, "Lower energy limit (TeV)"  
35 -emax, r, a, 125.8925,,, "Upper energy limit (TeV)"  
36 -bins, i, a, 18,,, "Number of energy bins for sensitivity computation"  
37 -enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)"  
38 -npix, i, h, 200,,, "Number of pixels for binned analysis"  
39 -binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)"  
40 -type, s, h, Differential,Differential|Integral,, "Sensitivity type"  
41 -sigma, r, h, 5.0,,, "Significance threshold (Gaussian sigma)"  
42 -max_iter, i, h, 50,,, "Maximum number of iterations"  
43 -statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"  
44 -min_src_events, i, h, -1,,, "minimal number of source counts"  
45 -non_source_uncert_frac, r, h, -1,,, "fraction of no-source counts serving as minimal threshold for source counts"  
46 -non_source_uncert_rad, r, h, 0.33,,, "radius for circular region in which to integrate no-source counts for threshold"  
47 -seed, i, h, 0,,, "optional seed parameter for generating events" 31 +offset, r, h, 0.0,,, "Source offset angle (deg)"
  32 +duration, r, a, 180000.0,,, "Effective exposure time (s)"
  33 +rad, r, a, 5.0,,, "Radius of ROI (deg)"
  34 +emin, r, a, 0.03163,,, "Lower energy limit (TeV)"
  35 +emax, r, a, 125.8925,,, "Upper energy limit (TeV)"
  36 +bins, i, a, 18,,, "Number of energy bins for sensitivity computation"
  37 +enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)"
  38 +npix, i, h, 200,,, "Number of pixels for binned analysis"
  39 +binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)"
  40 +type, s, h, Differential,Differential|Integral,, "Sensitivity type"
  41 +sigma, r, h, 5.0,,, "Significance threshold (Gaussian sigma)"
  42 +max_iter, i, h, 50,,, "Maximum number of iterations"
  43 +statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
  44 +mincounts, i, h, 0,,, "Minimal number of source counts"
  45 +bkgcut, r, h, 0.0,,, "Background uncertainty fraction"
  46 +bkgrad, r, h, 0.33,,, "Background radius (deg)"
  47 +seed, i, h, 1,,, "Initial random number generator seed"
48 48
49 # 49 #
50 # Standard parameters 50 # Standard parameters
cscripts/cssens.py
  View file @955940f
@@ -20,7 +20,6 @@ @@ -20,7 +20,6 @@
20 # ========================================================================== 20 # ==========================================================================
21 import sys 21 import sys
22 import math 22 import math
23 -import statistics  
24 import gammalib 23 import gammalib
25 import ctools 24 import ctools
26 from cscripts import obsutils 25 from cscripts import obsutils
@@ -61,6 +60,7 @@ class cssens(ctools.csobservation): @@ -61,6 +60,7 @@ class cssens(ctools.csobservation):
61 self._dec = None 60 self._dec = None
62 self._log_clients = False 61 self._log_clients = False
63 self._models = gammalib.GModels() 62 self._models = gammalib.GModels()
  63 + self._seed = 1
64 self._nthreads = 0 64 self._nthreads = 0
65 65
66 # Return 66 # Return
@@ -86,6 +86,7 @@ class cssens(ctools.csobservation): @@ -86,6 +86,7 @@ class cssens(ctools.csobservation):
86 'dec' : self._dec, 86 'dec' : self._dec,
87 'log_clients' : self._log_clients, 87 'log_clients' : self._log_clients,
88 'models' : self._models, 88 'models' : self._models,
  89 + 'seed' : self._seed,
89 'nthreads' : self._nthreads} 90 'nthreads' : self._nthreads}
90 91
91 # Return pickled dictionary 92 # Return pickled dictionary
@@ -109,6 +110,7 @@ class cssens(ctools.csobservation): @@ -109,6 +110,7 @@ class cssens(ctools.csobservation):
109 self._dec = state['dec'] 110 self._dec = state['dec']
110 self._log_clients = state['log_clients'] 111 self._log_clients = state['log_clients']
111 self._models = state['models'] 112 self._models = state['models']
  113 + self._seed = state['seed']
112 self._nthreads = state['nthreads'] 114 self._nthreads = state['nthreads']
113 115
114 # Return 116 # Return
@@ -151,6 +153,10 @@ class cssens(ctools.csobservation): @@ -151,6 +153,10 @@ class cssens(ctools.csobservation):
151 self['type'].string() 153 self['type'].string()
152 self['edisp'].boolean() 154 self['edisp'].boolean()
153 self['debug'].boolean() 155 self['debug'].boolean()
  156 + self['mincounts'].integer()
  157 +
  158 + # Read seed
  159 + self._seed = self['seed'].integer()
154 160
155 # Derive some parameters 161 # Derive some parameters
156 self._ebounds = gammalib.GEbounds(bins, 162 self._ebounds = gammalib.GEbounds(bins,
@@ -170,6 +176,41 @@ class cssens(ctools.csobservation): @@ -170,6 +176,41 @@ class cssens(ctools.csobservation):
170 # Return 176 # Return
171 return 177 return
172 178
  179 + def _median(self, array):
  180 + """
  181 + Compute median value for an array
  182 +
  183 + Parameters
  184 + ----------
  185 + array : list of floats
  186 + Array
  187 +
  188 + Returns
  189 + -------
  190 + median : float
  191 + Median value of array
  192 + """
  193 + # Initialise median
  194 + median = 0.0
  195 +
  196 + # Get number of elements in array
  197 + n = len(array)
  198 +
  199 + # Continue only if there are elements in the array
  200 + if n > 0:
  201 +
  202 + # Sort the array
  203 + wrk_array = sorted(array)
  204 +
  205 + # Get median
  206 + if n % 2 != 0:
  207 + median = wrk_array[int(n/2)]
  208 + else:
  209 + median = (wrk_array[int((n-1)/2)] + wrk_array[int(n/2)]) / 2.0
  210 +
  211 + # Return median
  212 + return median
  213 +
173 def _set_source(self): 214 def _set_source(self):
174 """ 215 """
175 Set source name and position 216 Set source name and position
@@ -371,12 +412,12 @@ class cssens(ctools.csobservation): @@ -371,12 +412,12 @@ class cssens(ctools.csobservation):
371 npix = 200 412 npix = 200
372 binsz = 0.05 413 binsz = 0.05
373 414
374 - # keep track of the original value of edisp  
375 - edisp_orig = self['edisp'].boolean()  
376 - # set edisp off for the first iterations 415 + # Keep track of the original value of edisp and switch-off energy
  416 + # dispersion for the first iterations to speed-up the computations
  417 + edisp_orig = self['edisp'].boolean()
377 self['edisp'] = False 418 self['edisp'] = False
378 419
379 - # optional seed parameter for generating events 420 + # Get seed parameter for generating events
380 seed = self['seed'].integer() 421 seed = self['seed'].integer()
381 422
382 # Set flux ratio precision required for convergence to 5% 423 # Set flux ratio precision required for convergence to 5%
@@ -431,18 +472,15 @@ class cssens(ctools.csobservation): @@ -431,18 +472,15 @@ class cssens(ctools.csobservation):
431 472
432 # Simulate events for the models. "sim" holds an observation 473 # Simulate events for the models. "sim" holds an observation
433 # container with observations containing the simulated events. 474 # container with observations containing the simulated events.
434 - sim = obsutils.sim(self.obs(), nbins=enumbins, seed=(iterations + seed), 475 + sim = obsutils.sim(self.obs(), nbins=enumbins, seed=(iterations+seed),
435 binsz=binsz, npix=npix, 476 binsz=binsz, npix=npix,
436 log=self._log_clients, 477 log=self._log_clients,
437 debug=self['debug'].boolean(), 478 debug=self['debug'].boolean(),
438 edisp=self['edisp'].boolean(), 479 edisp=self['edisp'].boolean(),
439 nthreads=1) 480 nthreads=1)
440 481
441 - # Determine number of events in simulation by summing the events  
442 - # over all observations in the observation container  
443 - nevents = 0.0  
444 - for run in sim:  
445 - nevents += run.events().number() 482 + # Determine number of events in simulation
  483 + nevents = sim.nobserved()
446 484
447 # Write simulation results into logger 485 # Write simulation results into logger
448 self._log_header3(gammalib.EXPLICIT, 'Simulation') 486 self._log_header3(gammalib.EXPLICIT, 'Simulation')
@@ -513,7 +551,7 @@ class cssens(ctools.csobservation): @@ -513,7 +551,7 @@ class cssens(ctools.csobservation):
513 551
514 # Predict Crab flux at threshold TS using a linear regression of 552 # Predict Crab flux at threshold TS using a linear regression of
515 # the log(TS) and log(crab_flux) values that have so far been 553 # the log(TS) and log(crab_flux) values that have so far been
516 - # computed. If not enough results are available than use a simple 554 + # computed. If not enough results are available then use a simple
517 # TS scaling relation. 555 # TS scaling relation.
518 correct = math.sqrt(ts_thres / ts) 556 correct = math.sqrt(ts_thres / ts)
519 if len(results) > 1: 557 if len(results) > 1:
@@ -521,8 +559,9 @@ class cssens(ctools.csobservation): @@ -521,8 +559,9 @@ class cssens(ctools.csobservation):
521 pred_crab_flux, regcoeff = self._predict_flux(results, ts_thres) 559 pred_crab_flux, regcoeff = self._predict_flux(results, ts_thres)
522 correct = pred_crab_flux / crab_flux 560 correct = pred_crab_flux / crab_flux
523 except Exception as e: 561 except Exception as e:
524 - # if the regression fails, retain the simple scaling with ts_thres  
525 - self._log_value(gammalib.TERSE, 'skipping failed regression with error:', e) 562 + # If the regression fails, retain the simple scaling with
  563 + # ts_thres
  564 + self._log_value(gammalib.TERSE, 'Skipping failed regression', e)
526 565
527 # Compute extrapolated fluxes based on the flux correction factor 566 # Compute extrapolated fluxes based on the flux correction factor
528 crab_flux = correct * crab_flux 567 crab_flux = correct * crab_flux
@@ -551,17 +590,24 @@ class cssens(ctools.csobservation): @@ -551,17 +590,24 @@ class cssens(ctools.csobservation):
551 self._log_value(gammalib.TERSE, 'Regression coefficient', 590 self._log_value(gammalib.TERSE, 'Regression coefficient',
552 regcoeff) 591 regcoeff)
553 592
554 - # if the flux has converged, check if the original value of  
555 - # edisp was set, and is not what has fo far been used 593 + # If the flux has converged then check if the original
  594 + # value of edisp was set, and is not what has fo far
  595 + # been used
556 if edisp_orig and not self['edisp'].boolean(): 596 if edisp_orig and not self['edisp'].boolean():
557 - # set edisp on and continue the calculation 597 +
  598 + # Set edisp on and continue the calculation
558 self['edisp'] = True 599 self['edisp'] = True
559 600
560 - value = 'set edisp on after initial convergence withouth it'  
561 - self._log_value(gammalib.TERSE, 'Converged result', value) 601 + # Log action
  602 + self._log_value(gammalib.TERSE, 'Converged result',
  603 + 'Now use energy dispersion after '
  604 + 'initial convergence without it')
  605 +
  606 + # ... otherwise finish the calculation after convergence
562 else: 607 else:
563 - # finish the calculation after convergence  
564 break 608 break
  609 +
  610 + # ... otherwise break with a zero flux
565 else: 611 else:
566 self._log_value(gammalib.TERSE, 'Not converged', 'Flux is zero') 612 self._log_value(gammalib.TERSE, 'Not converged', 'Flux is zero')
567 break 613 break
@@ -575,35 +621,42 @@ class cssens(ctools.csobservation): @@ -575,35 +621,42 @@ class cssens(ctools.csobservation):
575 ' Test ended after %d iterations.' % max_iter) 621 ' Test ended after %d iterations.' % max_iter)
576 break 622 break
577 623
578 - # make sure edisp is set correctly before the count cuts  
579 - if edisp_orig and not self['edisp'].boolean():  
580 - self['edisp'] = True 624 + # Write header for iterations for terse chatter level
  625 + if self._logTerse():
  626 + self._log_header3(gammalib.TERSE, 'Iterations for source counts cuts')
  627 +
  628 + # Recover original energy dispersion setting
  629 + self['edisp'] = edisp_orig
581 630
  631 + # Take into account the excess event cuts
582 n_bck_evts = None 632 n_bck_evts = None
583 for iterations in range(max_iter): 633 for iterations in range(max_iter):
584 - pass_evt_cut, n_bck_evts = self._sim_evt_excess(  
585 - enumbins=enumbins,  
586 - binsz=binsz,  
587 - npix=npix,  
588 - n_bck_evts=n_bck_evts,  
589 - ) 634 +
  635 + # Simulate event excess
  636 + pass_evt_cut, n_bck_evts = self._sim_evt_excess(enumbins=enumbins,
  637 + binsz=binsz,
  638 + npix=npix,
  639 + n_bck_evts=n_bck_evts)
  640 +
  641 + # Write results into logger
  642 + name = 'Iteration %d' % iterations
  643 + value = 'Fit=%9.4f mCrab Sens=%e erg/cm2/s' % (crab_flux*1000.0, sensitivity)
  644 + self._log_value(gammalib.TERSE, name, value)
  645 +
  646 + # If event excess cuts are passed then break now
590 if pass_evt_cut: 647 if pass_evt_cut:
591 break 648 break
592 649
593 - # increase the test flux if the cuts failed 650 + # ... otherwise increase the test flux
594 correct = 1 + ratio_precision 651 correct = 1 + ratio_precision
595 crab_flux = correct * crab_flux 652 crab_flux = correct * crab_flux
596 photon_flux = correct * photon_flux 653 photon_flux = correct * photon_flux
597 energy_flux = correct * energy_flux 654 energy_flux = correct * energy_flux
598 sensitivity = correct * sensitivity 655 sensitivity = correct * sensitivity
599 656
  657 + # ...
600 _ = self._set_src_prefactor(test_model, crab_unit, crab_flux) 658 _ = self._set_src_prefactor(test_model, crab_unit, crab_flux)
601 659
602 - # Write results into logger  
603 - name = 'Iteration %d of source count check' % iterations  
604 - value = 'Fit=%9.4f mCrab Sens=%e erg/cm2/s' % (crab_flux*1000.0, sensitivity)  
605 - self._log_value(gammalib.TERSE, name, value)  
606 -  
607 # Write fit results into logger 660 # Write fit results into logger
608 self._log_header3(gammalib.TERSE, 'Fit results') 661 self._log_header3(gammalib.TERSE, 'Fit results')
609 self._log_value(gammalib.TERSE, 'Photon flux', 662 self._log_value(gammalib.TERSE, 'Photon flux',
@@ -656,37 +709,111 @@ class cssens(ctools.csobservation): @@ -656,37 +709,111 @@ class cssens(ctools.csobservation):
656 crab_prefactor : float 709 crab_prefactor : float
657 the prefactor that corresponds to a flux of 1 Crab. 710 the prefactor that corresponds to a flux of 1 Crab.
658 """ 711 """
659 - 712 + # Initialise variables
660 models = test_model.copy() 713 models = test_model.copy()
661 prefactor = modutils.normalisation_parameter(models[self._srcname]) 714 prefactor = modutils.normalisation_parameter(models[self._srcname])
662 crab_prefactor = prefactor.value() * crab_unit 715 crab_prefactor = prefactor.value() * crab_unit
  716 + val_margin = 0.01
  717 + min_pref = prefactor.min() * (1.0 + val_margin)
  718 + max_pref = prefactor.max() * (1.0 - val_margin)
  719 + val_now = crab_prefactor * test_crab_flux
663 720
664 - val_margin = 0.01  
665 - min_pref = prefactor.min() * (1 + val_margin)  
666 - max_pref = prefactor.max() * (1 - val_margin)  
667 -  
668 - val_now = crab_prefactor * test_crab_flux  
669 - 721 + # Check whether prefactor is in valid range
670 if val_now < min_pref or val_now > max_pref: 722 if val_now < min_pref or val_now > max_pref:
671 - log_name = 'test prefactor beyond limits'  
672 - log_value = (  
673 - 'test prefactor is too low/high, having value '  
674 - + str(val_now) + ' compared to limits: ['  
675 - + str(prefactor.min()) + ' , ' + str(prefactor.max()) + ']'  
676 - )  
677 - self._log_value(gammalib.TERSE, log_name, log_value)  
678 -  
679 - val_now = max(val_now, min_pref)  
680 - val_now = min(val_now, max_pref)  
681 - 723 +
  724 + # Store old prefactor value for logging
  725 + val_old = val_now
  726 +
  727 + # Set new prefactor value
  728 + val_now = max(val_now, min_pref)
  729 + val_now = min(val_now, max_pref)
682 crab_prefactor = val_now / test_crab_flux 730 crab_prefactor = val_now / test_crab_flux
683 731
  732 + # Log prefactor modification
  733 + self._log_value(gammalib.EXPLICIT, 'Prefactor range',
  734 + '['+str(min_pref)+','+str(max_pref)+']')
  735 + self._log_value(gammalib.EXPLICIT, 'Initial Prefactor', val_old)
  736 + self._log_value(gammalib.EXPLICIT, 'Updated Prefactor', val_now)
  737 +
  738 + # Store prefactor value
684 prefactor.value(val_now) 739 prefactor.value(val_now)
685 740
  741 + # Update the models
686 self.obs().models(models) 742 self.obs().models(models)
687 743
  744 + # Return Prefactor
688 return crab_prefactor 745 return crab_prefactor
689 746
  747 + def _simulate_events(self, obs, rad, seed, enumbins, binsz, npix):
  748 + """
  749 + Simulate events
  750 +
  751 + Parameters
  752 + ----------
  753 + obs : `~gammalib.GObservations`
  754 + Observation container
  755 + rad : float
  756 + Simulation selection radius (degrees)
  757 + seed : integer
  758 + Random number generator seed
  759 + enumbins : integer
  760 + Number of energy bins
  761 + binsz : float
  762 + Pixel size
  763 + npix : integer
  764 + Number of pixels
  765 + """
  766 + # Initialise list of number of simulated events
  767 + n_sim_evts = []
  768 +
  769 + # Simulate events
  770 + for n_sim in range(15):
  771 +
  772 + # Increment the seed
  773 + seed += n_sim
  774 +
  775 + # Simulate observations
  776 + sim = obsutils.sim(obs,
  777 + nbins=enumbins,
  778 + seed=seed,
  779 + binsz=binsz,
  780 + npix=npix,
  781 + log=self._log_clients,
  782 + debug=self['debug'].boolean(),
  783 + edisp=self['edisp'].boolean(),
  784 + nthreads=1)
  785 +
  786 + # If a selection radius is given then select only the events within
  787 + # that selection radius
  788 + if rad > 0.0:
  789 +
  790 + # Run event selection tool
  791 + select = ctools.ctselect(sim)
  792 + select['rad'] = rad
  793 + select['emin'] = 'INDEF'
  794 + select['tmin'] = 'INDEF'
  795 + select.run()
  796 +
  797 + # Store number of events
  798 + n_sim_evts.append(select.obs().nobserved())
  799 +
  800 + # ... otherwise we simply store the number of simulated events
  801 + else:
  802 + n_sim_evts.append(sim.nobserved())
  803 +
  804 + # Determine median number of events
  805 + median = self._median(n_sim_evts)
  806 +
  807 + # Log results
  808 + if rad > 0.0:
  809 + name = 'Median source events'
  810 + else:
  811 + name = 'Median background events'
  812 + self._log_value(gammalib.EXPLICIT, name, median)
  813 +
  814 + # Return
  815 + return median
  816 +
690 def _sim_evt_excess(self, enumbins, binsz, npix, n_bck_evts): 817 def _sim_evt_excess(self, enumbins, binsz, npix, n_bck_evts):
691 """ 818 """
692 Return the number of excess events for the source model, compared 819 Return the number of excess events for the source model, compared
@@ -695,134 +822,77 @@ class cssens(ctools.csobservation): @@ -695,134 +822,77 @@ class cssens(ctools.csobservation):
695 Parameters 822 Parameters
696 ---------- 823 ----------
697 enumbins : int 824 enumbins : int
698 - number of bins for the observation simulation 825 + Number of bins for the observation simulation
699 binsz : float 826 binsz : float
700 - bin size for the observation simulation 827 + Bin size for the observation simulation
701 npix : int 828 npix : int
702 - pixel size for the observation simulation 829 + Pixel size for the observation simulation
703 n_bck_evts : int 830 n_bck_evts : int
704 - number of non-source counts (given to avoid computing the same quantity multiple times) 831 + Number of background counts from previous call
705 832
706 Returns 833 Returns
707 ------- 834 -------
708 pass_evt_cut : bool 835 pass_evt_cut : bool
709 - does the source pass the minimal excess criteria 836 + Signals that the source passes the minimal excess criteria
710 n_bck_evts : int 837 n_bck_evts : int
711 - number of non-source counts 838 + Number of background counts
712 """ 839 """
713 -  
714 - # number of times to simulate the observation  
715 - n_sims = 15  
716 -  
717 - # radius of integration for hte non-source background (not necessarily the  
718 - # same as self['rad'], to match the current standard for non-ctools analyses)  
719 - rad_deg = self['non_source_uncert_rad'].real()  
720 -  
721 - # minimal number of events for the source  
722 - min_src_events = self['min_src_events'].integer()  
723 - # threshold by fraction of non-source counts (5 times the assumed uncertainty  
724 - # on the background, which is taken as 1% of the non-source counts)  
725 - non_source_uncert_frac = self['non_source_uncert_frac'].real()  
726 -  
727 - # optional seed parameter for generating events  
728 - seed = self['seed'].integer()  
729 -  
730 - # end here if no cuts have been spesified  
731 - if min_src_events <= 0 and non_source_uncert_frac <= 0:  
732 - return True  
733 -  
734 - models = self.obs().models()  
735 -  
736 - src_model = gammalib.GModels()  
737 - bck_models = gammalib.GModels()  
738 - for s in range(models.size()):  
739 - model = models[s]  
740 - if model.name() == self._srcname:  
741 - src_model.append(model)  
742 - else:  
743 - bck_models.append(model)  
744 -  
745 - src_obs = self.obs().copy()  
746 - src_obs.models(src_model)  
747 -  
748 - bck_obs = self.obs().copy()  
749 - bck_obs.models(bck_models)  
750 -  
751 - n_src_evts = 0  
752 - for is_src in [True, False]:  
753 - if not is_src:  
754 - # if the non-source fraction cut is not applied, no need  
755 - # to calculate the background rate.  
756 - # n_bck_evts will be None onlu in the firs iteration, where  
757 - # we only need to compute the non-source counts once, since the  
758 - # models dont change with the itterations.  
759 - if non_source_uncert_frac <= 0 or n_bck_evts is not None:  
760 - continue  
761 -  
762 - obs_now = src_obs if is_src else bck_obs  
763 -  
764 - n_sim_evts = []  
765 - for n_sim in range(n_sims):  
766 - seed = n_sim + seed  
767 -  
768 - sim_now = obsutils.sim(  
769 - obs_now,  
770 - nbins=enumbins,  
771 - seed=seed,  
772 - binsz=binsz,  
773 - npix=npix,  
774 - log=self._log_clients,  
775 - debug=self['debug'].boolean(),  
776 - edisp=self['edisp'].boolean(),  
777 - nthreads=1  
778 - )  
779 -  
780 - # consider all events in the entire RoI for the minimal cource count cut  
781 - if is_src:  
782 - n_sim_evts += [len(sim_now[0].events())]  
783 - # consider only background within the thetaSqr cut for the systematic  
784 - # check, in accordance with the current CTA convention for Li&Ma calculations 840 + # Initialise results
  841 + pass_evt_cut = True
  842 + n_bck_evts = 0
  843 +
  844 + # Get user parameters
  845 + mincounts = self['mincounts'].integer()
  846 + bkgcut = self['bkgcut'].real()
  847 + bkgrad = self['bkgrad'].real()
  848 +
  849 + # Continue only if cuts were specified
  850 + if mincounts > 0 or bkgcut > 0.0:
  851 +
  852 + # Get seed parameter for generating events
  853 + seed = self['seed'].integer()
  854 +
  855 + # Get models and split them into source and background
  856 + models = self.obs().models()
  857 + src_model = gammalib.GModels()
  858 + bck_models = gammalib.GModels()
  859 + for model in models:
  860 + if model.name() == self._srcname:
  861 + src_model.append(model)
785 else: 862 else:
786 - select = ctools.ctselect(sim_now)  
787 - select['rad'] = rad_deg  
788 - select['emin'] = 'INDEF'  
789 - select['tmin'] = 'INDEF'  
790 - select.run()  
791 -  
792 - obs = select.obs()  
793 - n_sim_evts += [len(obs[0].events())]  
794 -  
795 - if is_src:  
796 - n_src_evts = statistics.median(n_sim_evts)  
797 - else:  
798 - n_bck_evts = statistics.median(n_sim_evts)  
799 -  
800 - name = 'number of ' + ('source' if is_src else 'non-source') + ' events'  
801 - value = (  
802 - ','.join([str(i) for i in n_sim_evts])  
803 - + ' --> median is: ' + str(n_src_evts if is_src else n_bck_evts)  
804 - )  
805 - self._log_value(gammalib.TERSE, name, value)  
806 -  
807 - if n_bck_evts is None:  
808 - n_bck_evts = 0  
809 -  
810 - has_min_evts = n_src_evts >= min_src_events  
811 - is_above_bck = n_src_evts >= n_bck_evts * non_source_uncert_frac  
812 - pass_evt_cut = has_min_evts and is_above_bck  
813 -  
814 - name = 'minimal excess threshold'  
815 - value = (  
816 - 'has counts above \"min source counts('  
817 - + str(min_src_events) + ' events)\": ' + str(has_min_evts)  
818 - + ' ; pass background uncertainty threshold('  
819 - + str(max(non_source_uncert_frac, 0) * 100) + '% == '  
820 - + str(n_bck_evts * non_source_uncert_frac)  
821 - + ' events): ' + str(is_above_bck)  
822 - + ' --> pass selection: ' + str(pass_evt_cut)  
823 - )  
824 - self._log_value(gammalib.TERSE, name, value)  
825 - 863 + bck_models.append(model)
  864 +
  865 + # Create observations for source and background
  866 + src_obs = self.obs().copy()
  867 + bck_obs = self.obs().copy()
  868 + src_obs.models(src_model)
  869 + bck_obs.models(bck_models)
  870 +
  871 + # Simulate source events
  872 + n_src_evts = self._simulate_events(src_obs, 0.0, seed, enumbins, binsz, npix)
  873 +
  874 + # Simulate background events in case that a background fraction cut
  875 + # is applied and if the number of background events was not yet
  876 + # estimated
  877 + if bkgcut > 0.0 and n_bck_evts == None:
  878 + n_bck_evts = self._simulate_events(bck_obs, bkgrad, seed, enumbins, binsz, npix)
  879 +
  880 + # Set cut results
  881 + min_bkg_events = n_bck_evts * bkgcut
  882 + has_min_evts = n_src_evts >= mincounts
  883 + is_above_bck = n_src_evts >= min_bkg_events
  884 + pass_evt_cut = has_min_evts and is_above_bck
  885 +
  886 + # Log results
  887 + self._log_value(gammalib.EXPLICIT, 'Pass minimum counts cut', has_min_evts)
  888 + self._log_value(gammalib.EXPLICIT, 'Pass background cut', is_above_bck)
  889 + self._log_value(gammalib.EXPLICIT, 'Pass event cut', pass_evt_cut)
  890 + self._log_value(gammalib.EXPLICIT, 'Minimum counts threshold', mincounts)
  891 + self._log_value(gammalib.EXPLICIT, 'Background threshold', min_bkg_events)
  892 + self._log_value(gammalib.EXPLICIT, 'Source events', n_src_evts)
  893 + self._log_value(gammalib.EXPLICIT, 'Background events', n_bck_evts)
  894 +
  895 + # Return passing flag and number of background events
826 return pass_evt_cut, n_bck_evts 896 return pass_evt_cut, n_bck_evts
827 897
828 def _predict_flux(self, results, ts): 898 def _predict_flux(self, results, ts):
@@ -868,14 +938,13 @@ class cssens(ctools.csobservation): @@ -868,14 +938,13 @@ class cssens(ctools.csobservation):
868 mean_xx *= norm 938 mean_xx *= norm
869 mean_yy *= norm 939 mean_yy *= norm
870 940
  941 + # Compute regression coefficient
871 rxy_norm = (mean_xx - mean_x * mean_x) * (mean_yy - mean_y * mean_y) 942 rxy_norm = (mean_xx - mean_x * mean_x) * (mean_yy - mean_y * mean_y)
872 if rxy_norm < 1e-10: 943 if rxy_norm < 1e-10:
873 rxy_norm = 1 944 rxy_norm = 1
874 else: 945 else:
875 rxy_norm = math.sqrt(rxy_norm) 946 rxy_norm = math.sqrt(rxy_norm)
876 -  
877 - rxy = (mean_xy - mean_x * mean_y) / rxy_norm  
878 - 947 + rxy = (mean_xy - mean_x * mean_y) / rxy_norm
879 regcoeff = rxy*rxy 948 regcoeff = rxy*rxy
880 949
881 # Compute regression line slope 950 # Compute regression line slope
@@ -1122,4 +1191,4 @@ if __name__ == &#39;__main__&#39;: @@ -1122,4 +1191,4 @@ if __name__ == &#39;__main__&#39;:
1122 app = cssens(sys.argv) 1191 app = cssens(sys.argv)
1123 1192
1124 # Run application 1193 # Run application
1125 - app.execute()  
1126 \ No newline at end of file 1194 \ No newline at end of file
  1195 + app.execute()