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/examples/pipeline_binned_disk.py b/examples/pipeline_binned_disk.py index 4bbbca8..a6dffa7 100755 --- a/examples/pipeline_binned_disk.py +++ b/examples/pipeline_binned_disk.py @@ -155,4 +155,5 @@ if __name__ == '__main__': obs.models(gammalib.GModels(datadir+'/crab.xml')) # Run analysis pipeline - run_pipeline(obs, model=datadir+'/crab.xml') + run_pipeline(obs, model=datadir+'/crab.xml', enumbins=10, nxpix=40, + nypix=40, binsz=0.1) diff --git a/examples/pipeline_binned_mem.py b/examples/pipeline_binned_mem.py index 9abd992..30b0639 100755 --- a/examples/pipeline_binned_mem.py +++ b/examples/pipeline_binned_mem.py @@ -126,4 +126,4 @@ if __name__ == '__main__': obs.models(gammalib.GModels(datadir+'/crab.xml')) # Run analysis pipeline - run_pipeline(obs) + run_pipeline(obs, enumbins=10, nxpix=40, nypix=40, binsz=0.1) diff --git a/examples/pipeline_stacked_disk.py b/examples/pipeline_stacked_disk.py index 63b8951..008b3cf 100755 --- a/examples/pipeline_stacked_disk.py +++ b/examples/pipeline_stacked_disk.py @@ -205,4 +205,5 @@ if __name__ == '__main__': obs.models(gammalib.GModels(datadir+'/crab.xml')) # Run analysis pipeline - run_pipeline(obs, model=datadir+'/crab.xml') + run_pipeline(obs, model=datadir+'/crab.xml', enumbins=10, nxpix=40, + nypix=40, binsz=0.1) diff --git a/examples/pipeline_stacked_mem.py b/examples/pipeline_stacked_mem.py index 9e91f5b..083b64e 100755 --- a/examples/pipeline_stacked_mem.py +++ b/examples/pipeline_stacked_mem.py @@ -178,4 +178,4 @@ if __name__ == '__main__': obs.models(gammalib.GModels(datadir+'/crab.xml')) # Run analysis pipeline - run_pipeline(obs) + run_pipeline(obs, enumbins=10, nxpix=40, nypix=40, binsz=0.1) diff --git a/test/test_cspull.py b/test/test_cspull.py index 66e49b6..2a32961 100755 --- a/test/test_cspull.py +++ b/test/test_cspull.py @@ -146,7 +146,7 @@ class Test(test): pull = cscripts.cspull() pull['inmodel'] = self._model pull['outfile'] = 'cspull_py2.dat' - pull['ntrials'] = 2 + pull['ntrials'] = 1 pull['caldb'] = self._caldb pull['irf'] = self._irf pull['ra'] = 83.6331 @@ -157,8 +157,8 @@ class Test(test): pull['tmax'] = 100.0 pull['deadc'] = 0.95 pull['rad'] = 5.0 - pull['npix'] = 100 - pull['binsz'] = 0.02 + pull['npix'] = 20 + pull['binsz'] = 0.2 pull['coordsys'] = 'CEL' pull['proj'] = 'TAN' pull['logfile'] = 'cspull_py2.log' @@ -168,14 +168,14 @@ class Test(test): pull.execute() # Check pull distribution file - self._check_pull_file('cspull_py2.dat') + self._check_pull_file('cspull_py2.dat', rows=2) # Set-up cspull from event list pull = cscripts.cspull() pull['inobs'] = self._events pull['inmodel'] = self._model pull['outfile'] = 'cspull_py3.dat' - pull['ntrials'] = 2 + pull['ntrials'] = 1 pull['caldb'] = self._caldb pull['irf'] = self._irf pull['enumbins'] = 0 @@ -186,7 +186,7 @@ class Test(test): pull.execute() # Check pull distribution file - self._check_pull_file('cspull_py3.dat') + self._check_pull_file('cspull_py3.dat', rows=2) # Build observation container with unbinned observation cta = gammalib.GCTAObservation(self._events) @@ -197,7 +197,7 @@ class Test(test): pull = cscripts.cspull(obs) pull['inmodel'] = self._model pull['outfile'] = 'cspull_py4.dat' - pull['ntrials'] = 2 + pull['ntrials'] = 1 pull['caldb'] = self._caldb pull['irf'] = self._irf pull['enumbins'] = 0 @@ -208,7 +208,7 @@ class Test(test): pull.execute() # Check pull distribution file - self._check_pull_file('cspull_py4.dat') + self._check_pull_file('cspull_py4.dat', rows=2) # Set-up stacked cspull with one observation where response cubes # are specified in the input observation @@ -216,7 +216,7 @@ class Test(test): pull['inobs'] = self._inobs pull['inmodel'] = self._stacked_model pull['outfile'] = 'cspull_py5.dat' - pull['ntrials'] = 2 + pull['ntrials'] = 1 pull['enumbins'] = 0 pull['logfile'] = 'cspull_py5.log' pull['chatter'] = 4 @@ -225,7 +225,7 @@ class Test(test): pull.execute() # Check pull distribution file - self._check_pull_file('cspull_py5.dat') + self._check_pull_file('cspull_py5.dat', rows=2) # Set-up stacked cspull with two observations from which response # cubes will be computed internally. The IRFs are also specified @@ -235,7 +235,7 @@ class Test(test): pull['inobs'] = self._inobs_two pull['inmodel'] = self._model pull['outfile'] = 'cspull_py6.dat' - pull['ntrials'] = 2 + pull['ntrials'] = 1 pull['emin'] = 0.02 pull['emax'] = 100.0 pull['enumbins'] = 10 @@ -250,13 +250,13 @@ class Test(test): pull.execute() # Check pull distribution file - self._check_pull_file('cspull_py6.dat') + self._check_pull_file('cspull_py6.dat', rows=2) # Return return # Check pull file - def _check_pull_file(self, filename): + def _check_pull_file(self, filename, rows=3): """ Check pull file """ @@ -264,8 +264,8 @@ class Test(test): pulls = gammalib.GCsv(filename) # Check dimensions - self.test_value(pulls.nrows(), 3, - 'Check for 3 rows in pull file') + self.test_value(pulls.nrows(), rows, + 'Check for number of rows in pull file') # Return return 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 diff --git a/test/test_ctmapcube.py b/test/test_ctmapcube.py index 22137ee..e5604db 100644 --- a/test/test_ctmapcube.py +++ b/test/test_ctmapcube.py @@ -71,8 +71,8 @@ class Test(test): # Setup ctmapcube command cmd = ctmapcube+' inmodel="'+self._model+'"'+ \ ' outcube="ctmapcube_cmd1.fits"'+\ - ' emin=0.1 emax=100.0 enumbins=20 ebinalg=LOG'+ \ - ' nxpix=200 nypix=200 binsz=0.02 coordsys=CEL'+ \ + ' emin=0.1 emax=100.0 enumbins=10 ebinalg=LOG'+ \ + ' nxpix=100 nypix=100 binsz=0.04 coordsys=CEL'+ \ ' xref=83.63 yref=22.01 proj=CAR'+ \ ' logfile="ctmapcube_cmd1.log" chatter=1' @@ -90,8 +90,8 @@ class Test(test): # Setup ctmapcube command cmd = ctmapcube+' inmodel="events_that_do_not_exist.fits"'+ \ ' outcube="ctmapcube_cmd2.fits"'+\ - ' emin=0.1 emax=100.0 enumbins=20 ebinalg=LOG'+ \ - ' nxpix=200 nypix=200 binsz=0.02 coordsys=CEL'+ \ + ' emin=0.1 emax=100.0 enumbins=10 ebinalg=LOG'+ \ + ' nxpix=100 nypix=100 binsz=0.04 coordsys=CEL'+ \ ' xref=83.63 yref=22.01 proj=CAR'+ \ ' logfile="ctmapcube_cmd2.log"' @@ -114,10 +114,10 @@ class Test(test): mapcube['ebinalg'] = 'LOG' mapcube['emin'] = 0.1 mapcube['emax'] = 100.0 - mapcube['enumbins'] = 20 - mapcube['nxpix'] = 200 - mapcube['nypix'] = 200 - mapcube['binsz'] = 0.02 + mapcube['enumbins'] = 10 + mapcube['nxpix'] = 100 + mapcube['nypix'] = 100 + mapcube['binsz'] = 0.04 mapcube['coordsys'] = 'CEL' mapcube['proj'] = 'CAR' mapcube['xref'] = 83.63 @@ -152,12 +152,12 @@ class Test(test): cube = gammalib.GModelSpatialDiffuseCube(filename) # Test map cube - #self.test_value(cube.maps(), 20, "20 maps") - #self.test_value(cube.pixels(), 40000, "40000 map pixels") + #self.test_value(cube.maps(), 10, "10 maps") + #self.test_value(cube.pixels(), 10000, "10000 map pixels") # The map cube is not loaded by default !!!! We should add the method so that # the attributes are known !!!! - self.test_value(cube.maps(), 0, '20 maps') - self.test_value(cube.pixels(), 0, '40000 map pixels') + self.test_value(cube.maps(), 0, '10 maps') + self.test_value(cube.pixels(), 0, '10000 map pixels') # Return return diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index 925694c..e81ad49 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -69,7 +69,7 @@ class Test(test): ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ ' outmap="cttsmap_cmd1.fits"'+ \ - ' nxpix=5 nypix=5 binsz=0.02'+ \ + ' nxpix=3 nypix=3 binsz=0.02'+ \ ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ ' logfile="cttsmap_cmd1.log" chatter=1' @@ -89,7 +89,7 @@ class Test(test): ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ ' outmap="cttsmap_cmd2.fits"'+ \ - ' nxpix=5 nypix=5 binsz=0.02'+ \ + ' nxpix=3 nypix=3 binsz=0.02'+ \ ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ ' logfile="cttsmap_cmd2.log" chatter=2' @@ -113,8 +113,8 @@ class Test(test): tsmap['caldb'] = self._caldb tsmap['irf'] = self._irf tsmap['outmap'] = 'cttsmap_py1.fits' - tsmap['nxpix'] = 5 - tsmap['nypix'] = 5 + tsmap['nxpix'] = 3 + tsmap['nypix'] = 3 tsmap['binsz'] = 0.02 tsmap['coordsys'] = 'CEL' tsmap['proj'] = 'CAR' @@ -149,14 +149,14 @@ class Test(test): # Check dimensions self.test_value(ts.naxis(), 2, 'Check for 2 dimensions') - self.test_value(ts.naxes(0), 5, 'Check for 5 pixels in X') - self.test_value(ts.naxes(1), 5, 'Check for 5 pixels in Y') + self.test_value(ts.naxes(0), 3, 'Check for 3 pixels in X') + self.test_value(ts.naxes(1), 3, 'Check for 3 pixels in Y') self.test_value(prefactor.naxis(), 2, 'Check for 2 dimensions') - self.test_value(prefactor.naxes(0), 5, 'Check for 5 pixels in X') - self.test_value(prefactor.naxes(1), 5, 'Check for 5 pixels in Y') + self.test_value(prefactor.naxes(0), 3, 'Check for 3 pixels in X') + self.test_value(prefactor.naxes(1), 3, 'Check for 3 pixels in Y') self.test_value(index.naxis(), 2, 'Check for 2 dimensions') - self.test_value(index.naxes(0), 5, 'Check for 5 pixels in X') - self.test_value(index.naxes(1), 5, 'Check for 5 pixels in Y') + self.test_value(index.naxes(0), 3, 'Check for 3 pixels in X') + self.test_value(index.naxes(1), 3, 'Check for 3 pixels in Y') # Return return diff --git a/test/test_examples.py.in b/test/test_examples.py.in index 1c2142d..d39aff9 100755 --- a/test/test_examples.py.in +++ b/test/test_examples.py.in @@ -186,7 +186,7 @@ class Test(gammalib.GPythonTestSuite): Test make_ts_distributions """ # Set arguments - args = '-n 2 -e 0 -t 1800 -d %s' % self._datadir + args = '-n 1 -e 0 -t 100 -d %s' % self._datadir # Execute script self._execute_python('make_ts_distributions', args=args)