From 6fc75aff4e295eec0a856c0980a73a43d21e44e0 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 20 May 2021 12:33:53 +0200 Subject: [PATCH] Rebase code and update documentation (#3668) --- ChangeLog | 5 +++-- NEWS | 8 +++++++- README.md | 2 +- cscripts/cssens.par | 4 ++-- cscripts/cssens.py | 57 +++++++++++++++++++++++++-------------------------------- doc/source/admin/release_history/2.0.rst | 8 +++++--- doc/source/users/reference_manual/cssens.rst | 17 +++++++++++++++++ 7 files changed, 60 insertions(+), 41 deletions(-) diff --git a/ChangeLog b/ChangeLog index d5e91a9..ca4e487 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,13 +1,14 @@ -2021-05-19 +2021-05-20 * Version 2.0.0 released ======================== - Use test source position in cssens (#3555) + Add source counts constraints to cssens (#3668) Add ASCL shields to "About" page (#3630) Use file function constructor in csspec (#3626) Add BINS method to csspec (#3606) Add hidden "fix_spat" parameter to cttsmap (#3562) + Use test source position in cssens (#3555) Add optimizer setter to ctlikelihood base class Add Sphinx tutorials test script (#3510) Add "refit_if_failed" parameter to ctlike diff --git a/NEWS b/NEWS index 8f6d1c0..e847bdd 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 2.0.0 -19 May 2021 +20 May 2021 Introduction @@ -283,6 +283,12 @@ the width, height and position angle of the rectangular sky regions (#3248). cssens - Computes sensitivity curve ----------------------------------- +Added constraints on the minimum number of source counts in the sensitivity +computations. The constraints are specified by the new parameters "mincounts", +"bkgexcess" and "bkgrad". In addition, a "seed" parameter was added allowing +to specify the seed value for the random number generators. Changes were also +made to improve the stability of the sensitivity computation (#3668). + The test source position is now used as pointing direction, eventually offset in galactic latitude by a number of degrees specified using the "offset" parameter (#3555). diff --git a/README.md b/README.md index 8e9f654..7e2b042 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ctools information ================== -* Version: 2.0.0.dev (19 May 2021) +* Version: 2.0.0.dev (20 May 2021) * GammaLib dependency: 2.0.0.dev [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/) diff --git a/cscripts/cssens.par b/cscripts/cssens.par index 3141230..38f3650 100755 --- a/cscripts/cssens.par +++ b/cscripts/cssens.par @@ -41,8 +41,8 @@ type, s, h, Differential,Differential|Integral,, "Sensitivity type" sigma, r, h, 5.0,,, "Significance threshold (Gaussian sigma)" max_iter, i, h, 50,,, "Maximum number of iterations" statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic" -mincounts, i, h, 0,,, "Minimal number of source counts" -bkgcut, r, h, 0.0,,, "Background uncertainty fraction" +mincounts, i, h, 0,,, "Minimum number of source counts" +bkgexcess, r, h, 0.0,,, "Background uncertainty fraction" bkgrad, r, h, 0.33,,, "Background radius (deg)" seed, i, h, 1,,, "Initial random number generator seed" diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 7d3b594..650ed3a 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -24,7 +24,6 @@ import gammalib import ctools from cscripts import obsutils from cscripts import modutils -from cscripts import ioutils from cscripts import mputils @@ -46,7 +45,7 @@ class cssens(ctools.csobservation): # Constructor def __init__(self, *argv): """ - Constructor. + Constructor """ # Initialise application by calling the appropriate class constructor self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) @@ -154,7 +153,7 @@ class cssens(ctools.csobservation): self['edisp'].boolean() self['debug'].boolean() self['mincounts'].integer() - + # Read seed self._seed = self['seed'].integer() @@ -417,9 +416,6 @@ class cssens(ctools.csobservation): edisp_orig = self['edisp'].boolean() self['edisp'] = False - # Get seed parameter for generating events - seed = self['seed'].integer() - # Set flux ratio precision required for convergence to 5% ratio_precision = 0.05 @@ -465,6 +461,9 @@ class cssens(ctools.csobservation): # Update iteration counter iterations += 1 + # Increment the seed value + self._seed += 1 + # Write header for iteration into logger self._log_header2(gammalib.EXPLICIT, 'Iteration '+str(iterations)) @@ -472,7 +471,7 @@ class cssens(ctools.csobservation): # Simulate events for the models. "sim" holds an observation # container with observations containing the simulated events. - sim = obsutils.sim(self.obs(), nbins=enumbins, seed=(iterations+seed), + sim = obsutils.sim(self.obs(), nbins=enumbins, seed=self._seed, binsz=binsz, npix=npix, log=self._log_clients, debug=self['debug'].boolean(), @@ -589,7 +588,7 @@ class cssens(ctools.csobservation): self._log_value(gammalib.TERSE, 'Converged flux ratio', ratio) self._log_value(gammalib.TERSE, 'Regression coefficient', regcoeff) - + # If the flux has converged then check if the original # value of edisp was set, and is not what has fo far # been used @@ -648,7 +647,7 @@ class cssens(ctools.csobservation): break # ... otherwise increase the test flux - correct = 1 + ratio_precision + correct = 1.0 + ratio_precision crab_flux = correct * crab_flux photon_flux = correct * photon_flux energy_flux = correct * energy_flux @@ -744,7 +743,7 @@ class cssens(ctools.csobservation): # Return Prefactor return crab_prefactor - def _simulate_events(self, obs, rad, seed, enumbins, binsz, npix): + def _simulate_events(self, obs, rad, enumbins, binsz, npix): """ Simulate events @@ -754,8 +753,6 @@ class cssens(ctools.csobservation): Observation container rad : float Simulation selection radius (degrees) - seed : integer - Random number generator seed enumbins : integer Number of energy bins binsz : float @@ -769,19 +766,14 @@ class cssens(ctools.csobservation): # Simulate events for n_sim in range(15): - # Increment the seed - seed += n_sim + # Increment the seed value + self._seed += 1 # Simulate observations - sim = obsutils.sim(obs, - nbins=enumbins, - seed=seed, - binsz=binsz, - npix=npix, - log=self._log_clients, + sim = obsutils.sim(obs, nbins=enumbins, seed=self._seed, binsz=binsz, + npix=npix, log=self._log_clients, debug=self['debug'].boolean(), - edisp=self['edisp'].boolean(), - nthreads=1) + edisp=self['edisp'].boolean(), nthreads=1) # If a selection radius is given then select only the events within # that selection radius @@ -843,14 +835,11 @@ class cssens(ctools.csobservation): # Get user parameters mincounts = self['mincounts'].integer() - bkgcut = self['bkgcut'].real() + bkgexcess = self['bkgexcess'].real() bkgrad = self['bkgrad'].real() # Continue only if cuts were specified - if mincounts > 0 or bkgcut > 0.0: - - # Get seed parameter for generating events - seed = self['seed'].integer() + if mincounts > 0 or bkgexcess > 0.0: # Get models and split them into source and background models = self.obs().models() @@ -869,16 +858,16 @@ class cssens(ctools.csobservation): bck_obs.models(bck_models) # Simulate source events - n_src_evts = self._simulate_events(src_obs, 0.0, seed, enumbins, binsz, npix) + n_src_evts = self._simulate_events(src_obs, 0.0, enumbins, binsz, npix) # Simulate background events in case that a background fraction cut # is applied and if the number of background events was not yet # estimated - if bkgcut > 0.0 and n_bck_evts == None: - n_bck_evts = self._simulate_events(bck_obs, bkgrad, seed, enumbins, binsz, npix) + if bkgexcess > 0.0 and n_bck_evts == None: + n_bck_evts = self._simulate_events(bck_obs, bkgrad, enumbins, binsz, npix) # Set cut results - min_bkg_events = n_bck_evts * bkgcut + min_bkg_events = n_bck_evts * bkgexcess has_min_evts = n_src_evts >= mincounts is_above_bck = n_src_evts >= min_bkg_events pass_evt_cut = has_min_evts and is_above_bck @@ -937,7 +926,7 @@ class cssens(ctools.csobservation): mean_xy *= norm mean_xx *= norm mean_yy *= norm - + # Compute regression coefficient rxy_norm = (mean_xx - mean_x * mean_x) * (mean_yy - mean_y * mean_y) if rxy_norm < 1e-10: @@ -1065,6 +1054,10 @@ class cssens(ctools.csobservation): table.card('SIGMA', self['sigma'].real(), '[sigma] Sensitivity threshold') table.card('MAX_ITER', self['max_iter'].integer(), 'Maximum number of iterations') table.card('STAT', self['statistic'].string(), 'Optimization statistic') + table.card('MINCOUNT', self['mincounts'].integer(), 'Minimum number of source counts') + table.card('BKGEXCES', self['bkgexcess'].real(), 'Background uncertainty fraction') + table.card('BKGRAD', self['bkgrad'].real(), '[deg] Background radius') + table.card('SEED', self['seed'].integer(), 'Seed value for random numbers') # Append filled columns to fits table table.append(e_mean) diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index c85d03b..ee35000 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -ctools 2.0 release (19 May 2021) +ctools 2.0 release (20 May 2021) ================================ Overview @@ -24,8 +24,8 @@ Bug fixes Improvements ------------ -* [`3555 `_] - - Use test source position in :ref:`cssens` +* [`3668 `_] - + Add source counts constraints to :ref:`cssens` * [`3630 `_] - Add ASCL shields to **About** page * [`3626 `_] - @@ -34,6 +34,8 @@ Improvements Add ``BINS`` method to :ref:`csspec` * [`3562 `_] - Add hidden ``fix_spat`` parameter to :ref:`cttsmap` +* [`3555 `_] - + Use test source position in :ref:`cssens` * Add optimizer setter to ``ctlikelihood`` base class. * Add ``refit_if_failed`` parameter to :ref:`ctlike` * [`3441 `_] - diff --git a/doc/source/users/reference_manual/cssens.rst b/doc/source/users/reference_manual/cssens.rst index 9443837..c6da43e 100755 --- a/doc/source/users/reference_manual/cssens.rst +++ b/doc/source/users/reference_manual/cssens.rst @@ -121,6 +121,23 @@ General parameters observations, and ``CHI2`` uses the Chi squared statistic for all binned or stacked observations. +``(mincounts = 0) [integer]`` + Constraint on the minimum number of required source counts. Conventionally, + a constraint for a minimum number of 10 source counts is applied for CTA + sensitivity estimates. If ``0`` is specified then no source counts limit + will be applied. + +``(bkgexcess = 0.0) [real]`` + Constraint on the minimum number of required source counts with respect to + the number of background counts. This value is a fraction that is + conventionally set to ``0.05`` (or 5%) which means that the number of source + counts needs to exceed 5% of the number of background counts. If ``0.0``is + specifie then no constraint will be applied. + +``(bkgrad = 0.33) [real]`` + Radius in degrees used to estimate the number of background counts + underlying the source. This radius is only used if ``bkgexcess > 0.0``. + Standard parameters ------------------- -- libgit2 0.21.4