"""
Background models for various telescopes and instruments
"""
import numpy
import json
import logging
import warnings
import os
if 'MAKESPHINXDOC' not in os.environ:
import sherpa.astro.ui as ui
from sherpa.stats import Cash, CStat
from sherpa.models.parameter import Parameter
from sherpa.models import ArithmeticModel, CompositeModel
from sherpa.astro.ui import *
from sherpa.astro.instrument import RSPModelNoPHA, RMFModelNoPHA
else:
CompositeModel, ArithmeticModel = object, object
class BackgroundModel(object):
pass
logbm = logging.getLogger('bxa.BackgroundModel')
logbm.setLevel(logging.INFO)
"""
Background model for Chandra, as for the CDFS.
Uses a flat continuum, two broad gaussians and 8 narrow instrumental lines.
"""
[docs]
class ChandraBackground(BackgroundModel):
def __init__(self, storage):
self.storage = storage
i = self.storage.i
self.centers = [1.486, 1.739, 2.142, 7.478, 8.012, 8.265, 8.4939, 9.7133]
continuum, softsoftend, softend = box1d('continuum_%s' % i), gauss1d('softsoftend_%s' % i), gauss1d('softend_%s' % i)
line1, line2, line3, line4, line5, line6, line7, line8 = [
gauss1d('line%d_%s' % (j, i)) for j in range(1, 9)]
self.lines = [line1, line2, line3, line4, line5, line6, line7, line8]
self.abslines = [line5, line6, line7]
self.linelines = [line1, line2, line3, line4, line8]
#lines = linelines + abslines
self.boxes = [continuum]
self.softboxes = [softend, softsoftend]
self.plboxes = self.boxes + self.softboxes
logbm.info('creating Chandra background model')
for l, c in zip(self.lines, self.centers):
l.pos = c
l.pos.min = c - 0.1
l.pos.max = c + 0.1
for l in self.lines:
l.ampl = 2000
l.ampl.min = 1e-8
l.ampl.max = 1e12
l.fwhm = 0.02
l.fwhm.min = 0.02
l.pos.freeze()
for l in self.abslines:
l.fwhm.max = 0.4
for l in self.linelines:
l.fwhm.max = 0.1
# narrow soft end
softsoftend.pos = 0.3
softsoftend.pos.min = 0
softsoftend.pos.max = 0.6
softsoftend.pos.freeze()
softsoftend.fwhm = 0.5
softsoftend.fwhm.min = 0.2
softsoftend.fwhm.max = 0.7
softsoftend.fwhm.freeze()
# wide gaussian
softend.pos = 0
softend.pos.max = 1
softend.pos.min = -1
softend.pos.freeze()
softend.fwhm = 3.8
softend.fwhm.min = 2
softend.fwhm.max = 7
softend.fwhm.freeze()
for b in self.plboxes:
b.ampl.val = 1
b.ampl.freeze()
changepoints = [0.2, 0.8, 1.3, 2.5, 8.2, 8.4, 8.75, 12]
changepoints = [0., 12]
for b,clo,chi in zip(self.boxes, changepoints[:-1], changepoints[1:]):
b.xlow = clo
b.xhi = chi
b.xlow.freeze()
b.xhi.freeze()
b.ampl.min = 1e-8
b.ampl.max = 1e8
b.ampl.val = 1
contlevel = const1d("contlevel_%s" % i)
contlevel.integrate = False
softlevel = const1d("softlevel_%s" % i)
softlevel.integrate = False
softsoftlevel = const1d("softsoftlevel_%s" % i)
softsoftlevel.integrate = False
contlevel.c0.min = 1e-6
softlevel.c0.min = 1e-4
softsoftlevel.c0.min = 1e-4
contlevel.c0.max = 1e10
softlevel.c0.max = 1e10
softsoftlevel.c0.max = 1e10
# init
contlevel.c0.val = 1e-3
softlevel.c0.val = 1e1
softsoftlevel.c0.val = 1e1
self.stages = ['continuum', 'softfeatures'] + ['line%d' % i for i, line in enumerate(self.linelines + self.abslines)] + ['full']
self.contlevel = contlevel
self.softlevel = softlevel
self.softsoftlevel = softsoftlevel
[docs]
def set_zero(self):
for l in self.lines + self.plboxes:
l.ampl.min = 0
l.ampl.val = 0
"""
range over which this model is valid
"""
[docs]
def set_filter(self):
ignore(None, 0.4)
ignore(9.8, None)
notice(0.4, 9.8)
[docs]
def set_model(self, stage):
i = self.storage.i
withsoft = stage not in ('continuum')
withlines = stage not in ('continuum', 'softfeatures')
logbm.info('stage "%s" for %s ...' % (stage, i))
continuum = self.boxes[0]
contlevel = self.contlevel
softlevel = self.softlevel
softsoftlevel = self.softsoftlevel
softend, softsoftend = self.softboxes
bg_model = continuum * contlevel
bunitrsp = self.storage.bunitrsp
set_bkg_model(i, bg_model)
set_bkg_full_model(i, bunitrsp(bg_model))
logbm.debug('zooming to %.1f %.1f' % (2.5, 7))
ignore(None, 2.5)
ignore(7, None)
notice(2.5, 7)
self.stagepars = [contlevel.c0]
self.pars = list(self.stagepars)
if stage == 'continuum':
return
logbm.debug('adding soft end...')
bg_model += (softend * softlevel + softsoftend * softsoftlevel) * contlevel
delete_bkg_model(i)
set_bkg_model(i, bg_model)
set_bkg_full_model(i, bunitrsp(bg_model))
ignore(None, 0.4)
ignore(2., None)
notice(0.4, 2)
logbm.debug('zooming to %.1f %.1f' % (0.4, 2.))
self.stagepars = [softsoftlevel.c0, softlevel.c0, softsoftlevel.c0]
self.pars += self.stagepars
if stage == 'softfeatures':
return
logbm.debug('adding lines...')
for j, l in enumerate(self.linelines + self.abslines):
bg_model += l * contlevel
set_bkg_model(i, bg_model)
set_bkg_full_model(i, bunitrsp(bg_model))
if l.ampl.frozen: continue
lo = l.pos.val - max(3*l.fwhm.val, 0.2)
hi = l.pos.val + max(3*l.fwhm.val, 0.2)
logbm.debug('zooming to %.1f %.1f' % (lo, hi))
ignore(None, lo)
ignore(hi, None)
notice(lo, hi)
self.stagepars = [l.ampl]
self.pars += self.stagepars + [l.fwhm]
if stage == 'line%d' % j:
return
self.set_filter()
# finally, full fit
self.stagepars = self.pars
"""
Background model for XMM/PN, as for the XMM-XXL.
Developed by Richard Sturm at MPE
Written for Sherpa by Zhu Liu
Adapted into this framework by Johannes Buchner
"""
class XMMPNBackground(BackgroundModel):
pncenters = [1.49165, 1.49165, 4.53177, 5.42516, 6.38155, 7.48675, 8.04087, 8.04087, 8.60924, 8.89395, 9.56160]
pnlinewidth = [5.73813E-02, 3.63469E-02, 6.10487E-02, 7.08380E-02, 9.59053E-02, 6.52422E-02, 9.48594E-02, 6.26174E-05, 0.120893, 0.114254, 0.108717]
pnlinenorm = [7.81356E-03, 3.96601E-03, 7.30727E-04, 4.96413E-04, 5.31295E-9, 6.84796E-04, 3.01564E-02, 1.41847E-04, 8.87887E-03, 5.75592E-03, 1.71367E-03]
def __init__(self, storage, galactic_NH = 1e20):
self.storage = storage
i = self.storage.i
self.pnrsp = get_response(i)
self.pnbrsp = get_response(i,bkg_id=1)
self.pnscale = get_bkg_scale(i)
copy_data(i,1000+2)
self.pnbunitrsp = get_response(1000+2, bkg_id=1)
delete_data(1000+2)
self.pnbkgcons = xsconstant('pncons_%s' % i)
self.pnbkgspline1 = xsspline('pnspline1_%s' % i)
self.pnbkgexpdec = xsexpdec('pnexpdec_%s' % i)
self.pnbkgsmedge1 = xssmedge('pnsmedge1_%s' % i)
self.pnbkgspline2 = xsspline('pnspline2_%s' % i)
self.pnbkgsmedge2 = xssmedge('pnsmedge2_%s' % i)
self.pnbkginspl = xspowerlaw('pnbkpl_%s' % i)
self.pnlines = []
for j in range(1, 11+1):
comp = xsgaussian('pngau%d_%s' % (j,i))
self.pnlines.append(comp)
setattr(self, 'pnbkgline%d' % j, comp)
self.pnbkgpl = xspowerlaw('pnpnexpl_%s' % i)
self.pnbkgapec = xsapec('pnapec_%s' % i)
self.pnbkglcapec = xsapec('pnlcapec_%s' % i)
self.galabs = xsphabs('absgal_%s' % i)
self.galabs.nH = galactic_NH / 1e22
self.galabs.nH.freeze()
self.pnfixwid = [self.pnbkgline2, self.pnbkgline3, self.pnbkgline4, self.pnbkgline5, self.pnbkgline6, self.pnbkgline8, self.pnbkgline11]
self.pnfree = [self.pnbkgline1, self.pnbkgline7, self.pnbkgline9, self.pnbkgline10]
#define PN background model
#pn_bkg = self.pnbunitrsp(self.pnbkgcons*(self.pnbkgspline1*self.pnbkgexpdec + self.pnbkgsmedge1*self.pnbkgsmedge2*(self.pnbkgspline2*self.pnbkginspl +
# self.pnbkgline1 + self.pnbkgline2 + self.pnbkgline3 + self.pnbkgline4 + self.pnbkgline5 + self.pnbkgline6 + self.pnbkgline7 + self.pnbkgline8 +
# self.pnbkgline9 + self.pnbkgline10 + self.pnbkgline11))) + self.pnbrsp(self.galabs*(self.pnbkgpl+self.pnbkgapec)+self.pnbkglcapec)
# Define source model: galabs*(intgalabs*(srcpl+torus+pnsoftpl))
#set_model(i,galabs*(intgalabs*(xspowerlaw.srcpl)))
#set_full_model(i,(pnrsp(galabs*(intgalabs*(xspowerlaw.srcpl))) + pnscale*(pn_bkg)))
for l, c in zip(self.pnlines, self.pncenters):
l.LineE = c
l.LineE.min = c - 0.05
l.LineE.max = c + 0.05
self.pnbkgline2.LineE = self.pnbkgline1.LineE
self.pnbkgline8.LineE = self.pnbkgline7.LineE
for l, s in zip(self.pnlines, self.pnlinewidth):
l.Sigma = s
for l in self.pnfree:
l.Sigma.min = 1E-5
l.Sigma.max = 0.2
for l in self.pnfixwid:
l.Sigma.freeze()
for l, n in zip(self.pnlines, self.pnlinenorm):
l.norm = n
l.norm.min = 1E-10
l.norm.max = 1E10
# Scaling constant
self.pnbkgcons.factor = 1.0
self.pnbkgcons.factor.freeze()
# thermal radiation
self.pnbkgapec.kT = 0.286928
self.pnbkgapec.kT.min = 0.008
self.pnbkgapec.kT.max = 64
self.pnbkgapec.Abundanc = 1.0
self.pnbkgapec.Abundanc.freeze()
self.pnbkgapec.Redshift = 0.0
self.pnbkgapec.Redshift.freeze()
self.pnbkgapec.norm = 5.58410E-05
self.pnbkgapec.norm.min = 1e-10
self.pnbkgapec.norm.max = 1e10
# local thermal radiation
self.pnbkglcapec.kT = 0.1
self.pnbkglcapec.kT.freeze()
self.pnbkglcapec.Abundanc = 1.0
self.pnbkglcapec.Abundanc.freeze()
self.pnbkglcapec.Redshift = 0.0
self.pnbkglcapec.Redshift.freeze()
self.pnbkglcapec.norm = 3.89164E-05
self.pnbkglcapec.norm.min = 1e-10
self.pnbkglcapec.norm.max = 1e10
# exponential decay
self.pnbkgexpdec.factor = 44.3418
self.pnbkgexpdec.factor.min = 0
self.pnbkgexpdec.factor.max = 100
self.pnbkgexpdec.norm = 6830.89
self.pnbkgexpdec.norm.freeze()
# Smear function
self.pnbkgsmedge1.edgeE = 0.538408
self.pnbkgsmedge1.edgeE.freeze()
self.pnbkgsmedge1.MaxTau = 1.40238
self.pnbkgsmedge1.MaxTau.min = 0
self.pnbkgsmedge1.MaxTau.max = 10
self.pnbkgsmedge1.index = -2.67000
self.pnbkgsmedge1.index.freeze()
self.pnbkgsmedge1.width = 0.313365
self.pnbkgsmedge1.width.min = 0.01
self.pnbkgsmedge1.width.max = 100
self.pnbkgsmedge2.edgeE = 1.38826
self.pnbkgsmedge2.edgeE.freeze()
self.pnbkgsmedge2.MaxTau.min = 0
self.pnbkgsmedge2.MaxTau.max = 10
self.pnbkgsmedge2.MaxTau = 9.37167
self.pnbkgsmedge2.index = -2.67000
self.pnbkgsmedge2.index.freeze()
self.pnbkgsmedge2.width = 5.7642
self.pnbkgsmedge2.width.min = 0.01
self.pnbkgsmedge2.width.max = 100
# Spline funtion
self.pnbkgspline1.Estart = 0.200000
self.pnbkgspline1.Estart.freeze()
self.pnbkgspline1.Ystart = -1.31506
self.pnbkgspline1.Ystart.min = -1e+6
self.pnbkgspline1.Ystart.max = 1e+6
self.pnbkgspline1.Yend = 1064.16
self.pnbkgspline1.Yend.min = -1e+6
self.pnbkgspline1.Yend.max = 1e+6
self.pnbkgspline1.YPstart = -106.183
self.pnbkgspline1.YPstart.min = -1e+6
self.pnbkgspline1.YPstart.max = 1e+6
self.pnbkgspline1.YPend = -366.092
self.pnbkgspline1.YPend.min = -1e+6
self.pnbkgspline1.YPend.max = 1e6
self.pnbkgspline1.Eend = 1.74715
self.pnbkgspline1.Eend.min = 0
self.pnbkgspline1.Eend.max = 100
self.pnbkgspline2.Estart = 3.29056
self.pnbkgspline2.Ystart = 1.00643
self.pnbkgspline2.Ystart.min = -1e+6
self.pnbkgspline2.Ystart.max = 1e+6
self.pnbkgspline2.Yend = 0.887026
self.pnbkgspline2.Yend.min = -1e+6
self.pnbkgspline2.Yend.max = 1e+6
self.pnbkgspline2.YPstart = -0.278401
self.pnbkgspline2.YPstart.min = -1e+6
self.pnbkgspline2.YPstart.max = 1e+6
self.pnbkgspline2.YPend = 4.84809E-03
self.pnbkgspline2.YPend.min = -1e+6
self.pnbkgspline2.YPend.max = 1e6
self.pnbkgspline2.Eend = 7.32701
self.pnbkgspline2.Eend.min = 0
self.pnbkgspline2.Eend.max = 100
# Background power law
self.pnbkginspl.PhoIndex = 0.279
self.pnbkginspl.PhoIndex.min = -2
self.pnbkginspl.PhoIndex.max = 9
self.pnbkginspl.norm = 8.23614E-03
self.pnbkginspl.norm.min = 1E-10
self.pnbkginspl.norm.max = 1E6
# Extragalactic background
self.pnbkgpl.PhoIndex = 1.46
self.pnbkgpl.PhoIndex.freeze()
self.pnbkgpl.norm = 1.25288E-05
self.pnbkgpl.norm.min = 1e-10
self.pnbkgpl.norm.max = 1e3
self.stages = []
self.stage_models = {}
self.stage_models_full = {}
base = self.pnbkgspline2*self.pnbkginspl
logbm.info( ' XMMPNBackground: creating stages for %s' % i)
self.addStage('1', base)
self.addStage('2', base + self.pnbkgline2)
self.addStage('3', base + self.pnbkgline1 + self.pnbkgline2)
self.addStage('4', base + self.pnbkgline1 + self.pnbkgline2 + self.pnbkgline7 + self.pnbkgline8)
self.addStage('5', base + self.pnbkgline1 + self.pnbkgline2 + self.pnbkgline7 + self.pnbkgline8 + self.pnbkgline9 + self.pnbkgline10)
lines = self.pnbkgline1 + self.pnbkgline2 + self.pnbkgline3 + self.pnbkgline4 + self.pnbkgline5 + self.pnbkgline6 + self.pnbkgline7 + self.pnbkgline8 + self.pnbkgline9 + self.pnbkgline10 + self.pnbkgline11
self.addStage('6', base + lines)
self.addStage('7', base + lines, self.galabs*(self.pnbkgapec))
self.addStage('8', base + lines, self.galabs*(self.pnbkgapec+self.pnbkgpl))
self.addStage('9', base + lines, self.galabs*(self.pnbkgapec+self.pnbkgpl)+self.pnbkglcapec)
logbm.info( ' XMMPNBackground: setup for %s completed' % i)
def addStage(self, name, instrument_bg, external_bg=None):
model = self.pnbkgcons*(self.pnbkgspline1*self.pnbkgexpdec + self.pnbkgsmedge1*self.pnbkgsmedge2*(instrument_bg))
full_model = self.pnbunitrsp(self.pnbkgcons*(self.pnbkgspline1*self.pnbkgexpdec + self.pnbkgsmedge1*self.pnbkgsmedge2*(instrument_bg)))
if external_bg is not None:
model = model + external_bg
full_model = full_model + self.pnrsp(external_bg)
self.stage_models[name] = model
self.stage_models_full[name] = full_model
self.stages.append(name)
"""
range over which this model is valid
"""
def set_filter(self):
i = self.storage.i
notice(None, None)
ignore(None, 0.3)
ignore(10.0, None)
def set_model(self, stage):
i = self.storage.i
logbm.info( ' XMMPNBackground: going to stage %s on %s' % (stage, i))
bg_model = self.stage_models[stage]
bg_model_full = self.stage_models_full[stage]
self.stagepars = list(bg_model.pars)
self.pars = list(self.stagepars)
delete_bkg_model(i)
logbm.info( ' XMMPNBackground: setting model')
set_bkg_model(i, bg_model)
set_bkg_full_model(i, bg_model_full)
self.set_filter()
"""
Background model for Swift/XRT.
Uses a flat continuum, two broad gaussians and 8 narrow instrumental lines.
"""
[docs]
class SwiftXRTBackground(BackgroundModel):
def __init__(self, storage):
self.storage = storage
i = self.storage.i
logbm.info( ' SwiftXRTBackground: setting up for %s' % i)
dip = box1d('dip_%s' % i)
pbknpl = xsbknpower('pbknpl_%s' % i)
g1, g2, g3, g4 = [gauss1d('gauss%d_%s' % (j,i)) for j in [1, 2, 3, 4]]
pbknpl.BreakE.min = 0.2
pbknpl.BreakE.max = 5
pbknpl.BreakE.val = 2
pbknpl.PhoIndx1.max = 4
pbknpl.PhoIndx2.max = 4
pbknpl.PhoIndx1.min = 0.8
pbknpl.PhoIndx2.min = 0.8
pbknpl.PhoIndx1.val = 2
pbknpl.PhoIndx2.val = 1.5
pbknpl.norm.min = 1e-10
pbknpl.norm.max = 1
pbknpl.norm.val = 0.004
lines = [(0.1, 0.7, 1.1), (2, 2.15, 2.5), (1, 1.2, 1.4), (0, 0.4, 0.5)]
for g, (lo, mid, hi) in zip([g1, g2, g3, g4], lines):
g.pos.val = mid
g.pos.min = lo
g.pos.max = hi
g.fwhm.val = 0.2
g.fwhm.max = 1
g.fwhm.min = 0.01
g.ampl.val = 0.01
g.ampl.max = 1
g.ampl.min = 1e-6
dip.xlow = 2
dip.xhi = 3
dip.xlow.min = 1.75
dip.xlow.max = 2.25
dip.xhi.min = 2.75
dip.xhi.max = 3.25
dip.ampl.val = 0.5
dip.ampl.min = 1e-3
dip.ampl.max = 1 - 1e-3
self._pars = [pbknpl, dip, g1, g2, g3, g4]
self.stages = [2, 3, 4, 5, 6, 7]
[docs]
def set_filter(self):
logbm.debug('SwiftXRTBackground: setting filter 0.3-5keV')
ignore(None, 0.3)
ignore(5, None)
notice(0.3, 5)
[docs]
def set_model(self, stage):
i = self.storage.i
[pbknpl, dip, g1, g2, g3, g4] = self._pars
logbm.info('stage "%s" for ID=%s ...' % (stage, i))
bunitrsp = self.storage.bunitrsp
model = stage
self.bg_model = pbknpl
self.stagepars = list(pbknpl.pars)
if model > 1:
self.bg_model = self.bg_model + g1
self.stagepars += list(g1.pars)
if model > 2:
self.bg_model = self.bg_model + g2
self.stagepars = list(g2.pars)
if model > 3:
self.bg_model = self.bg_model + g3
self.stagepars = list(g3.pars)
if model > 4:
self.bg_model = (1 - dip) * self.bg_model
self.stagepars = list(dip.pars)
if model > 5:
self.bg_model = self.bg_model + g4
self.stagepars = list(g4.pars)
self.pars = [p for p in self.bg_model.pars if not p.link]
if model > 6: # finally, full fit
self.stagepars = self.pars
set_bkg_model(i, self.bg_model)
set_bkg_full_model(i, bunitrsp(self.bg_model))
self.set_filter()
logbm.debug('background model set for stage "%s" for ID=%s' % (stage, i))
[docs]
class SwiftXRTWTBackground(SwiftXRTBackground):
[docs]
def set_filter(self):
logbm.debug('SwiftXRTBackground: setting filter 0.4-5keV')
ignore(None, 0.4)
ignore(5, None)
notice(0.4, 5)
__all__ = ['SwiftXRTBackground', 'SwiftXRTWTBackground', 'ChandraBackground']