diff --git a/cscripts/cssens.par b/cscripts/cssens.par index fbf139f..6966fe4 100755 --- a/cscripts/cssens.par +++ b/cscripts/cssens.par @@ -28,19 +28,22 @@ outfile, f, a, sensitivity.fits,,, "Output sensitivity FITS file" # # Script parameters #================== -offset, r, h, 0.0,,, "Source offset angle (deg)" -duration, r, a, 180000.0,,, "Effective exposure time (s)" -rad, r, a, 5.0,,, "Radius of ROI (deg)" -emin, r, a, 0.03163,,, "Lower energy limit (TeV)" -emax, r, a, 125.8925,,, "Upper energy limit (TeV)" -bins, i, a, 18,,, "Number of energy bins for sensitivity computation" -enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)" -npix, i, h, 200,,, "Number of pixels for binned analysis" -binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)" -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" +offset, r, h, 0.0,,, "Source offset angle (deg)" +duration, r, a, 180000.0,,, "Effective exposure time (s)" +rad, r, a, 5.0,,, "Radius of ROI (deg)" +emin, r, a, 0.03163,,, "Lower energy limit (TeV)" +emax, r, a, 125.8925,,, "Upper energy limit (TeV)" +bins, i, a, 18,,, "Number of energy bins for sensitivity computation" +enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)" +npix, i, h, 200,,, "Number of pixels for binned analysis" +binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)" +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" +min_src_events, i, h, -1,,, "minimal number of source counts" +non_source_uncert_frac, r, h, -1,,, "fraction of no-source counts serving as minimal threshold for source counts" +non_source_uncert_rad, r, h, 0.33,,, "radius for circular region in which to integrate no-source counts for threshold" # # Standard parameters diff --git a/cscripts/cssens.py b/cscripts/cssens.py index f057967..c3e9521 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -20,6 +20,7 @@ # ========================================================================== import sys import math +import statistics import gammalib import ctools from cscripts import obsutils @@ -418,15 +419,7 @@ class cssens(ctools.csobservation): # Write header for iteration into logger self._log_header2(gammalib.EXPLICIT, 'Iteration '+str(iterations)) - # Create a copy of the test models, set the normalisation parameter - # of the test source in the models, and append the models to the - # observation. "crab_prefactor" is the prefactor that corresponds - # to a flux of 1 Crab. - models = test_model.copy() - prefactor = modutils.normalisation_parameter(models[self._srcname]) - crab_prefactor = prefactor.value() * crab_unit - prefactor.value(crab_prefactor * test_crab_flux) - self.obs().models(models) + crab_prefactor = self._set_src_prefactor(test_model, crab_unit, test_crab_flux) # Simulate events for the models. "sim" holds an observation # container with observations containing the simulated events. @@ -560,6 +553,31 @@ class cssens(ctools.csobservation): ' Test ended after %d iterations.' % max_iter) break + n_bck_evts = None + for iterations in range(max_iter): + pass_evt_cut, n_bck_evts = self._sim_evt_excess( + enumbins=enumbins, + binsz=binsz, + npix=npix, + n_bck_evts=n_bck_evts, + ) + if pass_evt_cut: + break + + # increase the test flux if the cuts failed + correct = 1 + ratio_precision + crab_flux = correct * crab_flux + photon_flux = correct * photon_flux + energy_flux = correct * energy_flux + sensitivity = correct * sensitivity + + _ = self._set_src_prefactor(test_model, crab_unit, crab_flux) + + # Write results into logger + name = 'Iteration %d of source count check' % iterations + value = 'Fit=%9.4f mCrab Sens=%e erg/cm2/s' % (crab_flux*1000.0, sensitivity) + self._log_value(gammalib.TERSE, name, value) + # Write fit results into logger self._log_header3(gammalib.TERSE, 'Fit results') self._log_value(gammalib.TERSE, 'Photon flux', @@ -593,6 +611,181 @@ class cssens(ctools.csobservation): # Return result return result + # + def _set_src_prefactor(self, test_model, crab_unit, test_crab_flux): + """ + Create a copy of the test models, set the normalisation parameter + of the test source in the models, and append the models to the observation. + + Parameters + ---------- + test_model : `~gammalib.GModels` + Test source model + crab_unit : float + Crab unit factor + test_crab_flux : float + Test flux in Crab units (100 mCrab) + + Returns + ------- + crab_prefactor : float + the prefactor that corresponds to a flux of 1 Crab. + """ + + models = test_model.copy() + prefactor = modutils.normalisation_parameter(models[self._srcname]) + crab_prefactor = prefactor.value() * crab_unit + + val_margin = 0.01 + min_pref = prefactor.min() + max_pref = prefactor.max() + + val_now = crab_prefactor * test_crab_flux + val_now = max(val_now, min_pref * (1 + val_margin)) + val_now = min(val_now, max_pref * (1 - val_margin)) + crab_prefactor = val_now / test_crab_flux + + prefactor.value(val_now) + + self.obs().models(models) + + return crab_prefactor + + def _sim_evt_excess(self, enumbins, binsz, npix, n_bck_evts): + """ + Return the number of excess events for the source model, compared + to all other models + + Parameters + ---------- + enumbins : int + number of bins for the observation simulation + binsz : float + bin size for the observation simulation + npix : int + pixel size for the observation simulation + n_bck_evts : int + number of non-source counts (given to avoid computing the same quantity multiple times) + + Returns + ------- + pass_evt_cut : bool + does the source pass the minimal excess criteria + n_bck_evts : int + number of non-source counts + """ + + # number of times to simulate the observation + n_sims = 15 + + # radius of integration for hte non-source background (not necessarily the + # same as self['rad'], to match the current standard for non-ctools analyses) + rad_deg = self['non_source_uncert_rad'].real() + + # minimal number of events for the source + min_src_events = self['min_src_events'].integer() + # threshold by fraction of non-source counts (5 times the assumed uncertainty + # on the background, which is taken as 1% of the non-source counts) + non_source_uncert_frac = self['non_source_uncert_frac'].real() + + # end here if no cuts have been spesified + if min_src_events <= 0 and non_source_uncert_frac <= 0: + return True + + models = self.obs().models() + + src_model = gammalib.GModels() + bck_models = gammalib.GModels() + for s in range(models.size()): + model = models[s] + if model.name() == self._srcname: + src_model.append(model) + else: + bck_models.append(model) + + src_obs = self.obs().copy() + src_obs.models(src_model) + + bck_obs = self.obs().copy() + bck_obs.models(bck_models) + + n_src_evts = 0 + for is_src in [True, False]: + if not is_src: + # if the non-source fraction cut is not applied, no need + # to calculate the background rate. + # n_bck_evts will be None onlu in the firs iteration, where + # we only need to compute the non-source counts once, since the + # models dont change with the itterations. + if non_source_uncert_frac <= 0 or n_bck_evts is not None: + continue + + obs_now = src_obs if is_src else bck_obs + + n_sim_evts = [] + for n_sim in range(n_sims): + seed = n_sim + + sim_now = obsutils.sim( + obs_now, + nbins=enumbins, + seed=seed, + binsz=binsz, + npix=npix, + log=self._log_clients, + debug=self['debug'].boolean(), + edisp=self['edisp'].boolean(), + nthreads=1 + ) + + # consider all events in the entire RoI for the minimal cource count cut + if is_src: + n_sim_evts += [len(sim_now[0].events())] + # consider only background within the thetaSqr cut for the systematic + # check, in accordance with the current CTA convention for Li&Ma calculations + else: + select = ctools.ctselect(sim_now) + select['rad'] = rad_deg + select['emin'] = 'INDEF' + select['tmin'] = 'INDEF' + select.run() + + obs = select.obs() + n_sim_evts += [len(obs[0].events())] + + if is_src: + n_src_evts = statistics.median(n_sim_evts) + else: + n_bck_evts = statistics.median(n_sim_evts) + + name = 'number of ' + ('source' if is_src else 'non-source') + ' events' + value = ( + ','.join([str(i) for i in n_sim_evts]) + + ' --> median is: ' + str(n_src_evts if is_src else n_bck_evts) + ) + self._log_value(gammalib.TERSE, name, value) + + if n_bck_evts is None: + n_bck_evts = 0 + + has_min_evts = n_src_evts >= min_src_events + is_above_bck = n_src_evts >= n_bck_evts * non_source_uncert_frac + pass_evt_cut = has_min_evts and is_above_bck + + name = 'minimal excess threshold' + value = ( + 'has counts above \"min source counts(' + + str(min_src_events) + ' events)\": ' + str(has_min_evts) + + ' ; pass background uncertainty threshold(' + + str(max(non_source_uncert_frac, 0) * 100) + '% == ' + + str(n_bck_evts * non_source_uncert_frac) + + ' events): ' + str(is_above_bck) + + ' --> pass selection: ' + str(pass_evt_cut) + ) + self._log_value(gammalib.TERSE, name, value) + + return pass_evt_cut, n_bck_evts + def _predict_flux(self, results, ts): """ Predict Crab flux for a given TS value @@ -635,9 +828,15 @@ class cssens(ctools.csobservation): mean_xy *= norm mean_xx *= norm mean_yy *= norm - rxy = (mean_xy - mean_x * mean_y) / \ - math.sqrt((mean_xx - mean_x*mean_x) * - (mean_yy - mean_y*mean_y)) + + rxy_norm = (mean_xx - mean_x * mean_x) * (mean_yy - mean_y * mean_y) + if rxy_norm < 1e-10: + rxy_norm = 1 + else: + rxy_norm = math.sqrt(rxy_norm) + + rxy = (mean_xy - mean_x * mean_y) / rxy_norm + regcoeff = rxy*rxy # Compute regression line slope @@ -884,4 +1083,4 @@ if __name__ == '__main__': app = cssens(sys.argv) # Run application - app.execute() + app.execute() \ No newline at end of file