diff --git a/ChangeLog b/ChangeLog index ebe243c..26d0fb4 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.3.0 released ======================== + Implement regression scheme in cssens and remove num_avg parameter (#1989) Copy cscripts into bin directory instead of creating a symbolic links (#1987) If code originates from Git, always compile the Python bindings (#1985) Raise exception instead of exiting with return code in Python tests (#1983) diff --git a/NEWS b/NEWS index 60abd03..649abee 100644 --- a/NEWS +++ b/NEWS @@ -133,6 +133,10 @@ extension that can be used for a stacked analysis. The energy boundaries are set as function of the variation in the effective area and background rates of the Instrument Response Function (#1976). +A regression scheme now replaces the averaging scheme in cssens as this +allows more accurately the prediction of the flux value at the threshold +Test Statistic. The num_avg parameter has been removed (#1989). + Examples -------- diff --git a/cscripts/cssens.par b/cscripts/cssens.par index e79fbb4..5cda99c 100644 --- a/cscripts/cssens.par +++ b/cscripts/cssens.par @@ -20,7 +20,7 @@ inobs, f, h, NONE,,, "Input event list, counts cube, or observation definition inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file" srcname, s, a, Crab,,, "Source name" caldb, s, a, prod2,,, "Calibration database" -irf, s, a, South_0.5h,,, "Instrument response function" +irf, s, a, South_50h,,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" deadc, r, h, 0.98,0,1, "Deadtime correction factor" outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file" @@ -31,16 +31,15 @@ outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file" 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.020,,, "Lower energy limit (TeV)" -emax, r, a, 200.0,,, "Upper energy limit (TeV)" -bins, i, a, 21,,, "Number of energy bins for differential sensitivity computation" +emin, r, a, 0.030,,, "Lower energy limit (TeV)" +emax, r, a, 180.0,,, "Upper energy limit (TeV)" +bins, i, a, 20,,, "Number of energy bins for differential 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" -num_avg, i, h, 3,,, "Number of iterations for sliding average" # # Standard parameters @@ -49,6 +48,6 @@ chatter, i, h, 2,0,4, "Chattiness of output" clobber, b, h, yes,,, "Overwrite existing output files with new output files?" debug, b, h, no,,, "Debugging mode activated" mode, s, h, ql,,, "Mode of automatic parameters" -logfile, f, h, cspull.log,,, "Log filename" +logfile, f, h, cssens.log,,, "Log filename" diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 5d3a8b0..9136d1a 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# Computes the array sensitivity using a test source +# Computes the array sensitivity using the Test Statistic for a test source # -# Copyright (C) 2011-2016 Juergen Knoedlseder +# Copyright (C) 2011-2017 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -40,7 +40,7 @@ class cssens(ctools.cscript): and the flux level is determined that leads to a particular significance. The significance is determined using the Test Statistic, defined as twice - the likelihood difference between fitting with and without the test source. + the likelihood difference between fitting with and without the test source. """ # Constructor @@ -50,31 +50,20 @@ class cssens(ctools.cscript): """ # Set name self._name = 'cssens' - self._version = '1.2.0' + self._version = '1.3.0' # Initialise class members - self._obs = gammalib.GObservations() self._ebounds = gammalib.GEbounds() self._obs_ebounds = [] self._srcname = '' - self._outfile = gammalib.GFilename() self._ra = None self._dec = None - self._edisp = False - self._emin = 0.020 - self._emax = 200.0 - self._bins = 21 - self._enumbins = 0 - self._npix = 200 - self._binsz = 0.05 - self._type = 'Differential' - self._ts_thres = 25.0 - self._max_iter = 50 - self._num_avg = 3 self._log_clients = False - - # Initialise application by calling the appropriate base class - # constructor + + # Initialise observation container from constructor arguments + self._obs, argv = self._set_input_obs(argv) + + # Initialise application by calling the appropriate class constructor self._init_cscript(argv) # Return @@ -87,8 +76,8 @@ class cssens(ctools.cscript): Get user parameters from parfile """ # Set observation if not done before - if self._obs == None or self._obs.size() == 0: - self._obs = self._set_obs() + if self._obs.size() == 0: + self._obs = self._set_obs(self['emin'].real(), self['emax'].real()) # Set models if we have none if self._obs.models().size() == 0: @@ -98,31 +87,28 @@ class cssens(ctools.cscript): self._srcname = self['srcname'].string() # Read further parameters - self._outfile = self['outfile'].filename() - self._emin = self['emin'].real() - self._emax = self['emax'].real() - self._bins = self['bins'].integer() - - # Read parameters for binned if requested - self._enumbins = self['enumbins'].integer() - if not self._enumbins == 0: - self._npix = self['npix'].integer() - self._binsz = self['binsz'].real() - - # Read remaining parameters - self._edisp = self['edisp'].boolean() - self._ts_thres = self['sigma'].real() * self['sigma'].real() - self._max_iter = self['max_iter'].integer() - self._num_avg = self['num_avg'].integer() - self._type = self['type'].string() - - # Query remaining parameters + emin = self['emin'].real() + emax = self['emax'].real() + bins = self['bins'].integer() + + # Query parameters for binned if requested + enumbins = self['enumbins'].integer() + if not enumbins == 0: + self['npix'].integer() + self['binsz'].real() + + # Query input parameters + self['sigma'].real() + self['max_iter'].integer() + self['type'].string() + self['outfile'].filename() + self['edisp'].boolean() self['debug'].boolean() # Derive some parameters - self._ebounds = gammalib.GEbounds(self._bins, - gammalib.GEnergy(self._emin, 'TeV'), - gammalib.GEnergy(self._emax, 'TeV')) + self._ebounds = gammalib.GEbounds(bins, + gammalib.GEnergy(emin, 'TeV'), + gammalib.GEnergy(emax, 'TeV')) # Write input parameters into logger self._log_parameters(gammalib.TERSE) @@ -130,20 +116,20 @@ class cssens(ctools.cscript): # Return return - def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): + def _set_obs(self, emin, emax, lpnt=0.0, bpnt=0.0): """ Set an observation container Parameters ---------- + emin : float + Minimum energy (TeV) + emax : float + Maximum energy (TeV) lpnt : float, optional Galactic longitude of pointing (deg) bpnt : float, optional Galactic latitude of pointing (deg) - emin : float, optional - Minimum energy (TeV) - emax : float, optional - Maximum energy (TeV) Returns ------- @@ -205,23 +191,23 @@ class cssens(ctools.cscript): # Get energy boundaries of the observation obs_ebounds = self._obs_ebounds[i] - + # Get minimum and maximum energy of the observation obs_emin = obs_ebounds.emin() obs_emax = obs_ebounds.emax() - + # If [emin,emax] is fully contained in the observation energy range # the use [emin,emax] as energy boundaries if obs_emin <= emin and obs_emax >= emax: ebounds = gammalib.GEbounds(emin, emax) - + # ... otherwise, if [emin,emax] is completely outside the # observation energy range then set the energy boundaries to the # zero-width interval [0,0] elif emax < obs_emin or emin > obs_emax: e0 = gammalib.GEnergy(0.0, 'TeV') ebounds = gammalib.GEbounds(e0, e0) - + # ... otherwise, if [emin,emax] overlaps partially with the # observation energy range then set the energy boundaries to the # overlapping part @@ -229,10 +215,10 @@ class cssens(ctools.cscript): # Set overlapping energy range set_emin = max(emin, obs_emin) set_emax = min(emax, obs_emax) - + # Set energy boundaries ebounds = gammalib.GEbounds(set_emin, set_emax) - + # Set the energy boundaries as the boundaries of the observation obs.events().ebounds(ebounds) @@ -286,6 +272,20 @@ class cssens(ctools.cscript): # Set TeV->erg conversion factor tev2erg = 1.6021764 + # Set parameters + ts_thres = self['sigma'].real() * self['sigma'].real() + max_iter = self['max_iter'].integer() + enumbins = self['enumbins'].integer() + if not enumbins == 0: + npix = self['npix'].integer() + binsz = self['binsz'].real() + else: + npix = 200 + binsz = 0.05 + + # Set flux ratio precision required for convergence to 5% + ratio_precision = 0.05 + # Set energy boundaries self._set_obs_ebounds(emin, emax) @@ -300,20 +300,27 @@ class cssens(ctools.cscript): src_flux = test_model[self._srcname].spectral().flux(emin, emax) crab_unit = crab_flux/src_flux - # Write header and some introductory information into logger + # Initialise regression coefficient + regcoeff = 0.0 + + # Write header for energy bin + self._log_string(gammalib.TERSE, '') self._log_header2(gammalib.TERSE, 'Energies: '+str(emin)+' - '+str(emax)) + + # Write initial parameters + self._log_header3(gammalib.TERSE, 'Initial parameters') self._log_value(gammalib.TERSE, 'Crab flux', str(crab_flux)+' ph/cm2/s') - self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+ - ' ph/cm2/s') + self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+' ph/cm2/s') self._log_value(gammalib.TERSE, 'Crab unit factor', crab_unit) # Initialise loop - crab_flux_value = [] - photon_flux_value = [] - energy_flux_value = [] - sensitivity_value = [] - iterations = 0 - test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab) + results = [] + iterations = 0 + test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab) + + # Write header for iterations for terse chatter level + if self._logTerse(): + self._log_header3(gammalib.TERSE, 'Iterations') # Loop until we break while True: @@ -335,11 +342,11 @@ class cssens(ctools.cscript): # Simulate events for the models. "sim" holds an observation # container with observations containing the simulated events. - sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations, - binsz=self._binsz, npix=self._npix, + sim = obsutils.sim(self._obs, nbins=enumbins, seed=iterations, + binsz=binsz, npix=npix, log=self._log_clients, debug=self['debug'].boolean(), - edisp=self._edisp) + edisp=self['edisp'].boolean()) # Determine number of events in simulation by summing the events # over all observations in the observation container @@ -354,7 +361,7 @@ class cssens(ctools.cscript): # Fit test source to the simulated events in the observation # container fit = ctools.ctlike(sim) - fit['edisp'] = edisp=self._edisp + fit['edisp'] = self['edisp'].boolean() fit['debug'] = self['debug'].boolean() fit['chatter'] = self['chatter'].integer() fit.run() @@ -383,29 +390,24 @@ class cssens(ctools.cscript): # Assess quality based on a comparison between Npred and Nevents quality = npred - nevents - # Write test source fit results into logger - self._log_header3(gammalib.EXPLICIT, 'Test source model fit') - self._log_value(gammalib.EXPLICIT, 'Test statistics', ts) - self._log_value(gammalib.EXPLICIT, 'log likelihood', logL) - self._log_value(gammalib.EXPLICIT, 'Number of predicted events', quality) - self._log_value(gammalib.EXPLICIT, 'Fit quality', npred) - for model in models: - self._log_value(gammalib.EXPLICIT, 'Model', model.name()) - for par in model: - self._log_string(gammalib.EXPLICIT, str(par)) - - # If TS was non-positive then increase the test flux and start - # over + # Write fit results into logger + name = 'Iteration %d' % iterations + value = ('TS=%10.4f Sim=%9.4f mCrab Fit=%9.4f mCrab ' + 'Sens=%e erg/cm2/s' % + (ts, test_crab_flux*1000.0, crab_flux*1000.0, sensitivity)) + self._log_value(gammalib.TERSE, name, value) + + # If TS was non-positive then increase the test flux and start over if ts <= 0.0: # If the number of iterations was exceeded then stop - if (iterations >= self._max_iter): + if (iterations >= max_iter): self._log_string(gammalib.TERSE, - ' Test ended after %d iterations.' % self._max_iter) + ' Test ended after %d iterations.' % max_iter) break - # Increase test flux by a factor of 2 - test_crab_flux = test_crab_flux * 2.0 + # Increase test flux + test_crab_flux *= 3.0 # Signal start we start over self._log_string(gammalib.EXPLICIT, @@ -414,95 +416,64 @@ class cssens(ctools.cscript): # ... and start over continue - # Compute flux correction factor based on average TS. We - # assume here that the significance is proportional to the source - # flux, and that the significance scales as the square root of - # the Test Statistic - correct = math.sqrt(self._ts_thres/ts) + # Append result entry to result list + result = {'ts': ts, 'crab_flux': crab_flux, + 'photon_flux': photon_flux, + 'energy_flux': energy_flux} + results.append(result) + + # Predict Crab flux at threshold TS using a linear regression of + # the log(TS) and log(crab_flux) values that have so far been + # computed. If not enough results are available than use a simple + # TS scaling relation. + if len(results) > 1: + pred_crab_flux, regcoeff = self._predict_flux(results, ts_thres) + correct = pred_crab_flux / crab_flux + else: + correct = math.sqrt(ts_thres/ts) # Compute extrapolated fluxes based on the flux correction factor crab_flux = correct * crab_flux photon_flux = correct * photon_flux energy_flux = correct * energy_flux sensitivity = correct * sensitivity - crab_flux_value.append(crab_flux) - photon_flux_value.append(photon_flux) - energy_flux_value.append(energy_flux) - sensitivity_value.append(sensitivity) - - # Write fit results - if self._logExplicit(): - self._log_value(gammalib.EXPLICIT, 'Photon flux', - str(photon_flux)+' ph/cm2/s') - self._log_value(gammalib.EXPLICIT, 'Energy flux', - str(energy_flux)+' erg/cm2/s') - self._log_value(gammalib.EXPLICIT, 'Crab flux', - str(crab_flux*1000.0)+' mCrab') - self._log_value(gammalib.EXPLICIT, 'Differential sensitivity', - str(sensitivity)+' erg/cm2/s') - for model in models: - self._log_value(gammalib.EXPLICIT, 'Model', model.name()) - for par in model: - self._log_string(gammalib.EXPLICIT, str(par)) - elif self._logTerse(): - name = 'Iteration %d' % iterations - value = ('TS=%.4f corr=%e %e ph/cm2/s = %e erg/cm2/s = %.2f ' - 'mCrab = %e erg/cm2/s' % (ts, correct, photon_flux, - energy_flux, crab_flux*1000.0, sensitivity)) - self._log_value(gammalib.TERSE, name, value) - - # Compute sliding average of extrapolated fitted prefactor, - # photon and energy flux. This damps out fluctuations and - # improves convergence. - crab_flux = 0.0 - photon_flux = 0.0 - energy_flux = 0.0 - sensitivity = 0.0 - num = 0.0 - for k in range(self._num_avg): - inx = len(crab_flux_value) - k - 1 - if inx >= 0: - crab_flux += crab_flux_value[inx] - photon_flux += photon_flux_value[inx] - energy_flux += energy_flux_value[inx] - sensitivity += sensitivity_value[inx] - num += 1.0 - crab_flux /= num - photon_flux /= num - energy_flux /= num - sensitivity /= num - - # Compare average flux to last average - if iterations > self._num_avg: + + # If we have at least 3 results then check if the flux determination + # at the TS threshold has converged + if len(results) > 3: if test_crab_flux > 0: + + # Compute fractional change in the Crab flux between two + # iterations ratio = crab_flux/test_crab_flux - # We have 2 convergence criteria: - # 1. The average flux does not change - # 2. The flux correction factor is small - if ratio >= 0.99 and ratio <= 1.01 and \ - correct >= 0.9 and correct <= 1.1: - if self._logTerse(): - self._log(' Converged ('+str(ratio)+')\n') + # If fractional change is smaller than the required position + # the iterations are stopped + if ratio > 1.0-ratio_precision and \ + ratio < 1.0+ratio_precision: + value = ('TS=%10.4f Sim=%9.4f mCrab ' + ' Sens=%e erg/cm2/s' % + (ts, crab_flux*1000.0, sensitivity)) + self._log_value(gammalib.TERSE, 'Converged result', value) + self._log_value(gammalib.TERSE, 'Converged flux ratio', ratio) + self._log_value(gammalib.TERSE, 'Regression coefficient', + regcoeff) break else: - if self._logTerse(): - self._log(' Flux is zero.\n') + self._log_value(gammalib.TERSE, 'Not converged', 'Flux is zero') break - # Use average for next iteration + # Set test flux for next iteration test_crab_flux = crab_flux # Exit loop if number of trials exhausted - if (iterations >= self._max_iter): - if self._logTerse(): - self._log(' Test ended after '+str(self._max_iter)+ - ' iterations.\n') + if (iterations >= max_iter): + self._log_string(gammalib.TERSE, + ' Test ended after %d iterations.' % max_iter) break # Write fit results into logger self._log_header3(gammalib.TERSE, 'Fit results') - self._log_value(gammalib.TERSE, 'Test statistics', ts) self._log_value(gammalib.TERSE, 'Photon flux', str(photon_flux)+' ph/cm2/s') self._log_value(gammalib.TERSE, 'Energy flux', @@ -523,16 +494,84 @@ class cssens(ctools.cscript): # Restore energy boundaries of observation container for i, obs in enumerate(self._obs): obs.events().ebounds(self._obs_ebounds[i]) - + # Store result result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ 'crab_flux': crab_flux, 'photon_flux': photon_flux, \ 'energy_flux': energy_flux, \ - 'sensitivity': sensitivity} + 'sensitivity': sensitivity, 'regcoeff': regcoeff, \ + 'nevents': nevents, 'npred': npred} # Return result return result + def _predict_flux(self, results, ts): + """ + Predict Crab flux for a given TS value + + The Crab flux at a given Test Statistic value is predicted by doing a + linear regression of the log(TS) and log(crab_flux) values in a results + list. + + See https://en.wikipedia.org/wiki/Simple_linear_regression + + Parameters + ---------- + results : list of dict + List of results + ts : float + Test Statistic value + + Returns + ------- + crab_flux_prediction, regcoeff : tuple of float + Predicted Crab flux and regression coefficient + """ + # Compute means and regression coefficient + mean_x = 0.0 + mean_y = 0.0 + mean_xy = 0.0 + mean_xx = 0.0 + mean_yy = 0.0 + for result in results: + x = math.log(float(result['ts'])) + y = math.log(float(result['crab_flux'])) + mean_x += x + mean_y += y + mean_xy += x * y + mean_xx += x * x + mean_yy += y * y + norm = 1.0 / float(len(results)) + mean_x *= norm + mean_y *= norm + 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)) + regcoeff = rxy*rxy + + # Compute regression line slope + beta_nom = 0.0 + beta_denom = 0.0 + for result in results: + x = math.log(float(result['ts'])) + y = math.log(float(result['crab_flux'])) + beta_nom += (x - mean_x) * (y - mean_y) + beta_denom += (x - mean_x) * (x - mean_x) + beta = beta_nom / beta_denom + + # Compute regression line offset + alpha = mean_y - beta * mean_x + + # Predict Crab flux at TS threshold + log_ts_thres = math.log(ts) + crab_flux_prediction = math.exp(alpha + beta * log_ts_thres) + + # Return + return (crab_flux_prediction, regcoeff) + # Public methods def run(self): @@ -545,15 +584,16 @@ class cssens(ctools.cscript): # Get parameters self._get_parameters() - + # Loop over observations and store a deep copy of the energy # boundaries for later use for obs in self._obs: self._obs_ebounds.append(obs.events().ebounds().copy()) - + # Initialise script colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', - 'energy_flux', 'sensitivity'] + 'energy_flux', 'sensitivity', 'regcoeff', 'nevents', + 'npred'] results = [] # Set test source model for this observation @@ -566,30 +606,35 @@ class cssens(ctools.cscript): # Write models into logger self._log_models(gammalib.NORMAL, models, 'Model') + # Get sensitivity type + sensitivity_type = self['type'].string() + # Write header self._log_header1(gammalib.TERSE, 'Sensitivity determination') - self._log_value(gammalib.TERSE, 'Type', str(self._type)) + self._log_value(gammalib.TERSE, 'Type', sensitivity_type) # Loop over energy bins for ieng in range(self._ebounds.size()): # Set energies - if self._type == 'Differential': + if sensitivity_type == 'Differential': emin = self._ebounds.emin(ieng) emax = self._ebounds.emax(ieng) - elif self._type == 'Integral': + elif sensitivity_type == 'Integral': emin = self._ebounds.emin(ieng) emax = self._ebounds.emax() else: msg = ('Invalid sensitivity type "%s" encountered. Either ' - 'specify "Differential" or "Integral".' % self._type) + 'specify "Differential" or "Integral".' % + sensitivity_type) raise RuntimeError(msg) # Determine sensitivity result = self._get_sensitivity(emin, emax, models) # Write out trial result - ioutils.write_csv_row(self._outfile.url(), ieng, colnames, result) + ioutils.write_csv_row(self['outfile'].filename().url(), ieng, + colnames, result) # Append results results.append(result) diff --git a/doc/source/users/reference_manual/cssens.rst b/doc/source/users/reference_manual/cssens.rst index bf22f90..2cec9c5 100644 --- a/doc/source/users/reference_manual/cssens.rst +++ b/doc/source/users/reference_manual/cssens.rst @@ -100,9 +100,6 @@ General parameters ``(max_iter = 50) [integer]`` Maximum number of iterations. -``(num_avg = 3) [integer]`` - Number of iterations for sliding average. - Standard parameters ------------------- diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 7e232a6..1618df5 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -1471,7 +1471,8 @@ void ctobssim::simulate_background(GCTAObservation* obs, // Debug code: signal that we step into the response MC method #if defined(G_BACKGROUND_DEBUG) - std::cout << "ctobssim::simulate_background: model->mc in" << std::endl; + std::cout << "ctobssim::simulate_background: model->mc in" + << std::endl; #endif // Get simulated CTA event list. Note that this method @@ -1481,7 +1482,8 @@ void ctobssim::simulate_background(GCTAObservation* obs, // Debug code: signal that we came back #if defined(G_BACKGROUND_DEBUG) - std::cout << "ctobssim::simulate_background: model->mc out" << std::endl; + std::cout << "ctobssim::simulate_background: model->mc out" + << std::endl; #endif // Continue only if we got a CTA event list