From 75b1cd4290a1e949993dd7e9a17083d27534b4e6 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 2 May 2019 16:40:48 +0200 Subject: [PATCH] Add H.E.S.S. background lookup to models and enhance csbkgmodel unit test --- cscripts/Makefile.am | 3 ++- cscripts/csbkgmodel.py | 190 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------ models/Makefile.am | 3 ++- models/hess_bkg_lookup.fits | Bin 0 -> 8640 bytes test/data/hess_bkg_lookup.fits | Bin 0 -> 8640 bytes test/test_csbkgmodel.py | 109 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 177 insertions(+), 128 deletions(-) create mode 100644 models/hess_bkg_lookup.fits create mode 100644 test/data/hess_bkg_lookup.fits diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 651d69d..f408c3a 100755 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -235,7 +235,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/onoff_rmf2.fits.gz \ $(top_srcdir)/test/data/model_crab_radialacceptance.xml \ $(top_srcdir)/test/data/sigdist_map.fits.gz \ - $(top_srcdir)/test/data/variability.fits.gz + $(top_srcdir)/test/data/variability.fits.gz \ + $(top_srcdir)/test/data/hess_bkg_lookup.fits # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \ diff --git a/cscripts/csbkgmodel.py b/cscripts/csbkgmodel.py index 1f466ab..4b8637e 100755 --- a/cscripts/csbkgmodel.py +++ b/cscripts/csbkgmodel.py @@ -30,7 +30,6 @@ class csbkgmodel(ctools.csobservation): """ Generates background model """ - # Constructor def __init__(self, *argv): """ @@ -46,7 +45,7 @@ class csbkgmodel(ctools.csobservation): # Return return - # State methods por pickling + # __getstate__ method for pickling def __getstate__(self): """ Extend ctools.csobservation getstate method to include some members @@ -64,6 +63,7 @@ class csbkgmodel(ctools.csobservation): # Return pickled dictionary return state + # __setstate__ method for pickling def __setstate__(self, state): """ Extend ctools.csobservation setstate method to include some members @@ -102,8 +102,9 @@ class csbkgmodel(ctools.csobservation): self['smax'].real() if spatial == 'LOOKUP': self['slufile'].filename() - if spatial == 'LOOKUP' or spatial == 'GAUSS' or \ - spatial == 'GAUSS(E)' or spatial == 'PROFILE': + if spatial == 'LOOKUP' or spatial == 'GAUSS' or \ + spatial == 'GAUSS(E)' or spatial == 'PROFILE' or \ + spatial == 'POLYNOM': self['gradient'].boolean() spectral = self['spectral'].string() if spectral == 'NODES': @@ -178,133 +179,74 @@ class csbkgmodel(ctools.csobservation): spectral = gammalib.GModelSpectralPlaw(1.0e-13, -2.5, epivot) model = gammalib.GCTAModelAeffBackground(spectral) - # Handle LOOKUP model - elif self['spatial'].string() == 'LOOKUP': - - # Set spatial model - factor1 = gammalib.GCTAModelSpatialLookup(self['slufile'].filename()) - - # Optionally add gradient - if self['gradient'].boolean(): - spatial = gammalib.GCTAModelSpatialMultiplicative() - factor2 = gammalib.GCTAModelSpatialGradient() - spatial.append(factor1) - spatial.append(factor2) - else: - spatial = factor1 - - # Set spectral model - epivot = gammalib.GEnergy(1.0, 'TeV') - spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot) - spectral['Prefactor'].min(1.0e-8) - - # Set background model - model = gammalib.GCTAModelBackground(spatial, spectral) - - # Handle GAUSS model - elif self['spatial'].string() == 'GAUSS': - - # Set spatial model - factor1 = gammalib.GCTAModelRadialGauss(3.0) - factor1['Sigma'].min(sigma_min) - factor1['Sigma'].max(sigma_max) - - # Optionally add gradient - if self['gradient'].boolean(): - spatial = gammalib.GCTAModelSpatialMultiplicative() - factor2 = gammalib.GCTAModelSpatialGradient() - spatial.append(factor1) - spatial.append(factor2) - else: - spatial = factor1 + # Handle other models + else: - # Set spectral model - epivot = gammalib.GEnergy(1.0, 'TeV') - spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot) - spectral['Prefactor'].min(1.0e-8) + # Handle LOOKUP model + if self['spatial'].string() == 'LOOKUP': - # Set background model - model = gammalib.GCTAModelBackground(spatial, spectral) + # Set spatial model + factor1 = gammalib.GCTAModelSpatialLookup(self['slufile'].filename()) - # Handle GAUSS(E) model - elif self['spatial'].string() == 'GAUSS(E)': + # Handle GAUSS model + elif self['spatial'].string() == 'GAUSS': - # Set spatial model - if self['snumbins'].integer() == 1: + # Set spatial model factor1 = gammalib.GCTAModelRadialGauss(3.0) factor1['Sigma'].min(sigma_min) factor1['Sigma'].max(sigma_max) - else: - emin = gammalib.GEnergy(self['smin'].real(), 'TeV') - emax = gammalib.GEnergy(self['smax'].real(), 'TeV') - energies = gammalib.GEnergies(self['snumbins'].integer(), emin, emax) - spectrum = gammalib.GModelSpectralConst(3.0) - nodes = gammalib.GModelSpectralNodes(spectrum, energies) - for i in range(nodes.nodes()): - nodes[i*2+1].min(sigma_min) - nodes[i*2+1].max(sigma_max) - nodes.autoscale() - factor1 = gammalib.GCTAModelSpatialGaussSpectrum(nodes) - # Optionally add gradient - if self['gradient'].boolean(): - spatial = gammalib.GCTAModelSpatialMultiplicative() - factor2 = gammalib.GCTAModelSpatialGradient() - spatial.append(factor1) - spatial.append(factor2) - else: - spatial = factor1 - - # Set spectral model - epivot = gammalib.GEnergy(1.0, 'TeV') - spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot) - spectral['Prefactor'].min(1.0e-8) - - # Set background model - model = gammalib.GCTAModelBackground(spatial, spectral) - - # Handle PROFILE model - elif self['spatial'].string() == 'PROFILE': - - # Set spatial model - factor1 = gammalib.GCTAModelRadialProfile(2.0, 4.0, 5.0) - factor1['Width'].min(1.0) - factor1['Width'].max(10.0) - factor1['Core'].min(1.0) - factor1['Core'].max(10.0) - factor1['Tail'].min(-10.0) - factor1['Tail'].max(10.0) - factor1['Tail'].free() - - # Optionally add gradient - if self['gradient'].boolean(): - spatial = gammalib.GCTAModelSpatialMultiplicative() - factor2 = gammalib.GCTAModelSpatialGradient() - spatial.append(factor1) - spatial.append(factor2) + # Handle GAUSS(E) model + elif self['spatial'].string() == 'GAUSS(E)': + + # Set spatial model. If a single energy bin is specified then + # use the GAUSS model. Otherwise allocate a node spectrum for + # the energy nodes of GAUSS(E). + if self['snumbins'].integer() == 1: + factor1 = gammalib.GCTAModelRadialGauss(3.0) + factor1['Sigma'].min(sigma_min) + factor1['Sigma'].max(sigma_max) + else: + emin = gammalib.GEnergy(self['smin'].real(), 'TeV') + emax = gammalib.GEnergy(self['smax'].real(), 'TeV') + energies = gammalib.GEnergies(self['snumbins'].integer(), emin, emax) + spectrum = gammalib.GModelSpectralConst(3.0) + nodes = gammalib.GModelSpectralNodes(spectrum, energies) + for i in range(nodes.nodes()): + nodes[i*2+1].min(sigma_min) + nodes[i*2+1].max(sigma_max) + nodes.autoscale() + factor1 = gammalib.GCTAModelSpatialGaussSpectrum(nodes) + + # Handle PROFILE model + elif self['spatial'].string() == 'PROFILE': + + # Set spatial model + factor1 = gammalib.GCTAModelRadialProfile(2.0, 4.0, 5.0) + factor1['Width'].min(1.0) + factor1['Width'].max(10.0) + factor1['Core'].min(1.0) + factor1['Core'].max(10.0) + factor1['Tail'].min(1.0) + factor1['Tail'].max(10.0) + factor1['Tail'].free() + + # Handle POLYNOM model + elif self['spatial'].string() == 'POLYNOM': + + # Set spatial model + factor1 = gammalib.GCTAModelRadialPolynom([1.0, -0.1, +0.1]) + factor1['Coeff0'].min(0.1) + factor1['Coeff0'].max(10.0) + factor1['Coeff0'].fix() + factor1['Coeff1'].min(-10.0) + factor1['Coeff1'].max(10.0) + factor1['Coeff2'].min(-10.0) + factor1['Coeff2'].max(10.0) + + # Any other strings (should never occur) else: - spatial = factor1 - - # Set spectral model - epivot = gammalib.GEnergy(1.0, 'TeV') - spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot) - spectral['Prefactor'].min(1.0e-8) - - # Set background model - model = gammalib.GCTAModelBackground(spatial, spectral) - - # Handle POLYNOM model - elif self['spatial'].string() == 'POLYNOM': - - # Set spatial model - factor1 = gammalib.GCTAModelRadialPolynom([1.0, -0.1, +0.1]) - factor1['Coeff0'].min(0.1) - factor1['Coeff0'].max(10.0) - factor1['Coeff0'].fix() - factor1['Coeff1'].min(-10.0) - factor1['Coeff1'].max(10.0) - factor1['Coeff2'].min(-10.0) - factor1['Coeff2'].max(10.0) + model = None # Optionally add gradient if self['gradient'].boolean(): @@ -323,10 +265,6 @@ class csbkgmodel(ctools.csobservation): # Set background model model = gammalib.GCTAModelBackground(spatial, spectral) - # Any other strings (should never occur) - else: - model = None - # Set background name and instrument if model is not None: model.name('Background') diff --git a/models/Makefile.am b/models/Makefile.am index 859f199..ca18b19 100755 --- a/models/Makefile.am +++ b/models/Makefile.am @@ -37,7 +37,8 @@ models = 1dc_howto.xml \ crab_profile.xml \ disk.xml \ gauss.xml \ - shell.xml + shell.xml \ + hess_bkg_lookup.fits # Install models modelsdir = $(datadir)/models diff --git a/models/hess_bkg_lookup.fits b/models/hess_bkg_lookup.fits new file mode 100644 index 0000000..9f47d80 Binary files /dev/null and b/models/hess_bkg_lookup.fits differ diff --git a/test/data/hess_bkg_lookup.fits b/test/data/hess_bkg_lookup.fits new file mode 100644 index 0000000..9f47d80 Binary files /dev/null and b/test/data/hess_bkg_lookup.fits differ diff --git a/test/test_csbkgmodel.py b/test/test_csbkgmodel.py index 1ab4242..d5d96bb 100755 --- a/test/test_csbkgmodel.py +++ b/test/test_csbkgmodel.py @@ -45,6 +45,7 @@ class Test(test): # Set test datasets and parameters self._myevents1 = self._datadir + '/crab_offaxis1.fits' self._myevents2 = self._datadir + '/crab_offaxis2.fits' + self._lookup = self._datadir + '/hess_bkg_lookup.fits' # Return return @@ -250,6 +251,114 @@ class Test(test): # Check background model self._check_bkg_model('csbkgmodel_py5.xml', nmodels=2) + # Test GAUSS(E) spatial model + bkgmodel = cscripts.csbkgmodel() + bkgmodel['inobs'] = self._events + bkgmodel['caldb'] = self._caldb + bkgmodel['irf'] = self._irf + bkgmodel['instrument'] = 'CTA' + bkgmodel['spatial'] = 'GAUSS(E)' + bkgmodel['snumbins'] = 2 + bkgmodel['smin'] = 1.0 + bkgmodel['smax'] = 10.0 + bkgmodel['gradient'] = True + bkgmodel['spectral'] = 'NODES' + bkgmodel['ebinalg'] = 'LOG' + bkgmodel['emin'] = 1.0 + bkgmodel['emax'] = 100.0 + bkgmodel['enumbins'] = 8 + bkgmodel['runwise'] = True + bkgmodel['rad'] = 2.0 + bkgmodel['chatter'] = 2 + bkgmodel['outmodel'] = 'csbkgmodel_py6.xml' + bkgmodel['logfile'] = 'csbkgmodel_py6.log' + + # Run csbkgmodel script and save background model + bkgmodel.logFileOpen() # Make sure we get a log file + bkgmodel.execute() + + # Check background model + self._check_bkg_model('csbkgmodel_py6.xml') + + # Test LOOKUP spatial model + bkgmodel = cscripts.csbkgmodel() + bkgmodel['inobs'] = self._events + bkgmodel['caldb'] = self._caldb + bkgmodel['irf'] = self._irf + bkgmodel['instrument'] = 'CTA' + bkgmodel['spatial'] = 'LOOKUP' + bkgmodel['slufile'] = self._lookup + bkgmodel['gradient'] = True + bkgmodel['spectral'] = 'NODES' + bkgmodel['ebinalg'] = 'LOG' + bkgmodel['emin'] = 1.0 + bkgmodel['emax'] = 100.0 + bkgmodel['enumbins'] = 8 + bkgmodel['runwise'] = True + bkgmodel['rad'] = 2.0 + bkgmodel['chatter'] = 2 + bkgmodel['outmodel'] = 'csbkgmodel_py7.xml' + bkgmodel['logfile'] = 'csbkgmodel_py7.log' + + # Run csbkgmodel script and save background model + bkgmodel.logFileOpen() # Make sure we get a log file + bkgmodel.execute() + + # Check background model + self._check_bkg_model('csbkgmodel_py7.xml') + + # Test PROFILE spatial model + bkgmodel = cscripts.csbkgmodel() + bkgmodel['inobs'] = self._events + bkgmodel['caldb'] = self._caldb + bkgmodel['irf'] = self._irf + bkgmodel['instrument'] = 'CTA' + bkgmodel['spatial'] = 'PROFILE' + bkgmodel['gradient'] = True + bkgmodel['spectral'] = 'NODES' + bkgmodel['ebinalg'] = 'LOG' + bkgmodel['emin'] = 1.0 + bkgmodel['emax'] = 100.0 + bkgmodel['enumbins'] = 8 + bkgmodel['runwise'] = True + bkgmodel['rad'] = 2.0 + bkgmodel['chatter'] = 2 + bkgmodel['outmodel'] = 'csbkgmodel_py8.xml' + bkgmodel['logfile'] = 'csbkgmodel_py8.log' + + # Run csbkgmodel script and save background model + bkgmodel.logFileOpen() # Make sure we get a log file + bkgmodel.execute() + + # Check background model + self._check_bkg_model('csbkgmodel_py8.xml') + + # Test POLYNOM spatial model + bkgmodel = cscripts.csbkgmodel() + bkgmodel['inobs'] = self._events + bkgmodel['caldb'] = self._caldb + bkgmodel['irf'] = self._irf + bkgmodel['instrument'] = 'CTA' + bkgmodel['spatial'] = 'POLYNOM' + bkgmodel['gradient'] = True + bkgmodel['spectral'] = 'NODES' + bkgmodel['ebinalg'] = 'LOG' + bkgmodel['emin'] = 1.0 + bkgmodel['emax'] = 100.0 + bkgmodel['enumbins'] = 8 + bkgmodel['runwise'] = True + bkgmodel['rad'] = 2.0 + bkgmodel['chatter'] = 2 + bkgmodel['outmodel'] = 'csbkgmodel_py9.xml' + bkgmodel['logfile'] = 'csbkgmodel_py9.log' + + # Run csbkgmodel script and save background model + bkgmodel.logFileOpen() # Make sure we get a log file + bkgmodel.execute() + + # Check background model + self._check_bkg_model('csbkgmodel_py9.xml') + # Return return -- libgit2 0.21.4