From 1661676c505eb9389c997489972b0c80526c91db Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 1 Aug 2016 22:53:42 +0200 Subject: [PATCH] Fix regression bug introduces in cssens.py during refactoring --- cscripts/cssens.py | 65 +++++++++++++++++++++++++++++++++++++++++++++-------------------- test/test_cssens.py | 6 +++--- 2 files changed, 48 insertions(+), 23 deletions(-) diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 03e9a86..1d6cdf7 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -373,7 +373,17 @@ class cssens(ctools.cscript): logL = fit.opt().value() npred = fit.obs().npred() models = fit.obs().models() - ts = models[self._srcname].ts() + source = models[self._srcname] + ts = source.ts() + + # Get fitted Crab, photon and energy fluxes + crab_flux = source['Prefactor'].value() / crab_prefactor + photon_flux = source.spectral().flux(emin, emax) + energy_flux = source.spectral().eflux(emin, emax) + + # Compute differential sensitivity in unit erg/cm2/s + energy = gammalib.GEnergy(e_mean, 'TeV') + sensitivity = source.spectral().eval(energy) * e_mean*erg_mean * 1.0e6 # Assess quality based on a comparison between Npred and Nevents quality = npred - nevents @@ -400,27 +410,30 @@ class cssens(ctools.cscript): for par in model: self._log(str(par)+'\n') - # Start over if TS was non-positive + # If TS was non-positive then increase the test flux and start + # over if ts <= 0.0: - if self._logExplicit(): - self._log('Non positive TS. Start over.\n') - continue - # Get fitted Crab, photon and energy fluxes - crab_flux = models[self._srcname]['Prefactor'].value() / \ - crab_prefactor - photon_flux = models[self._srcname].spectral().flux(emin, emax) - energy_flux = models[self._srcname].spectral().eflux(emin, emax) + # If the number of iterations was exceeded then stop + if (iterations >= self._max_iter): + if self._logTerse(): + self._log(' Test ended after %d iterations.\n' % + self._max_iter) + break - # Compute differential sensitivity in unit erg/cm2/s - energy = gammalib.GEnergy(e_mean, 'TeV') - sensitivity = models[self._srcname].spectral().eval(energy) * \ - e_mean*erg_mean * 1.0e6 + # Increase test flux by a factor of 2 + test_crab_flux = test_crab_flux * 2.0 + + # Signal start we start over + if self._logExplicit(): + self._log('Non positive TS, increase test flux and start ' + 'over.\n') + + # ... and start over + continue # Compute flux correction factor based on average TS - correct = 1.0 - if ts > 0: - correct = math.sqrt(self._ts_thres/ts) + correct = math.sqrt(self._ts_thres/ts) # Compute extrapolated fluxes crab_flux = correct * crab_flux @@ -432,7 +445,7 @@ class cssens(ctools.cscript): energy_flux_value.append(energy_flux) sensitivity_value.append(sensitivity) - # Write background and test source fit results + # Write fit results if self._logExplicit(): self._log.parformat('Photon flux') self._log(photon_flux) @@ -453,7 +466,7 @@ class cssens(ctools.cscript): for par in model: self._log(str(par)+'\n') elif self._logTerse(): - self._log.parformat('Iteration '+str(iterations)) + self._log.parformat('Iteration %d' % iterations) self._log('TS=') self._log(ts) self._log(' ') @@ -589,7 +602,19 @@ class cssens(ctools.cscript): results = [] # Set test source model for this observation - models = modutils.test_source(self._obs.models(), self._srcname) + models = modutils.test_source(self._obs.models(), self._srcname, + ra=self._ra, dec=self._dec) + + # Write observations into logger + if self._logTerse(): + self._log('\n') + self._log.header1(gammalib.number('Observation',len(self._obs))) + self._log(str(self._obs)) + self._log('\n') + if self._logExplicit(): + for obs in self._obs: + self._log(str(obs)) + self._log('\n') # Write models into logger if self._logTerse(): diff --git a/test/test_cssens.py b/test/test_cssens.py index c7b5570..2f2bd84 100755 --- a/test/test_cssens.py +++ b/test/test_cssens.py @@ -72,7 +72,7 @@ class Test(test): cmd = cssens+' inmodel="'+self._model+'"'+ \ ' srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ - ' duration=500.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \ + ' duration=1800.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \ ' outfile="cssens_cmd1.dat"'+ \ ' logfile="cssens_cmd1.log" chatter=1' @@ -91,7 +91,7 @@ class Test(test): cmd = cssens+' inmodel="model_that_does_not_exist.xml"'+ \ ' srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ - ' duration=500.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \ + ' duration=1800.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \ ' outfile="cssens_cmd2.dat"'+ \ ' logfile="cssens_cmd2.log" chatter=2' @@ -115,7 +115,7 @@ class Test(test): sens['caldb'] = self._caldb sens['irf'] = self._irf sens['outfile'] = 'cssens_py1.dat' - sens['duration'] = 500.0 + sens['duration'] = 1800.0 sens['rad'] = 3.0 sens['emin'] = 1.0 sens['emax'] = 10.0 -- libgit2 0.21.4