Source code for bxa.sherpa.background.xmm

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
from __future__ import print_function

"""
Empirical XMM background model.

Written by Zhu Liu, adapted by Torben Simm and Johannes Buchner
(C) 2013-2016

For example usage, see examples/sherpa/background/xmm/fit.py
"""	

import os
import numpy

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

print("""

Using XMM empirical background model originally by Richard Sturm.
Please reference Maggi P., et al., 2014, A&A, 561, AA76.

""")


[docs] def get_embedded_file(filename): """ Gets the path of a file in the same folder as this script """ return os.path.join(os.path.dirname(__file__), filename)
[docs] def get_pn_bkg_model(i, galabs, fit=False): #================================================================= # Parameters: # i = sherpa ID of data set # galabs = name of model for galactic absorption # pnbrsp = response model of bkg data set # fit: True = fit bkg spectrum; False = just set_bkg_full_model # returns unit response for bkg model #================================================================= #get instrument response model for bkg spectrum if get_rmf(i).energ_lo[0] == 0: get_rmf(i).energ_lo[0] = 0.001 if get_arf(i).energ_lo[0] == 0: get_arf(i).energ_lo[0] = 0.001 pnbrsp = get_response(i,bkg_id=1) pnscale = get_bkg_scale(i) #get background scaling factor #create unit response dia_pn_rmf=get_embedded_file('pn_dia.rmf') dia_pn_arf=get_embedded_file('pn_dia.arf') copy_data(i,1000+2) load_bkg_rmf(1000+2, dia_pn_rmf) #load diagonal bkg matrices load_bkg_arf(1000+2, dia_pn_arf) pnbunitrsp = get_response(1000+2, bkg_id=1) delete_data(1000+2) #gaussian line center energy, line width, and initial normalization; for *PN background* 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] pnbkgcons, pnbkgspline1, pnbkgexpdec, pnbkgsmedge1, pnbkgsmedge2, pnbkgspline2, pnbkginspl, pnbkgline1, pnbkgline2, pnbkgline3, pnbkgline4, pnbkgline5, pnbkgline6, pnbkgline7, pnbkgline8, pnbkgline9, pnbkgline10, pnbkgline11, pnbkgpl, pnbkgapec, pnbkglcapec, = (xsconstant.pncons, xsspline.pnspline1, xsexpdec.pnexpdec, xssmedge.pnsmedge1, xssmedge.pnsmedge2, xsspline.pnspline2, xspowerlaw.pnbkpl, xsgaussian.pngau1, xsgaussian.pngau2, xsgaussian.pngau3, xsgaussian.pngau4, xsgaussian.pngau5, xsgaussian.pngau6, xsgaussian.pngau7, xsgaussian.pngau8, xsgaussian.pngau9, xsgaussian.pngau10, xsgaussian.pngau11, xspowerlaw.pnpnexpl, xsapec.pnapec, xsapec.pnlcapec) pnlines = [pnbkgline1, pnbkgline2, pnbkgline3, pnbkgline4, pnbkgline5, pnbkgline6, pnbkgline7, pnbkgline8, pnbkgline9, pnbkgline10, pnbkgline11] pnfixwid = [pnbkgline2, pnbkgline3, pnbkgline4, pnbkgline5, pnbkgline6, pnbkgline8, pnbkgline11] pnfree = [pnbkgline1, pnbkgline7, pnbkgline9, pnbkgline10] #define PN background model pn_bkg = pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline3 + pnbkgline4 + pnbkgline5 + pnbkgline6 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10 + pnbkgline11))) + pnbrsp(galabs*(pnbkgpl+pnbkgapec)+pnbkglcapec) for l, c in zip(pnlines, pncenters): l.LineE = c l.LineE.min = c - 0.05 l.LineE.max = c + 0.05 pnbkgline2.LineE = pnbkgline1.LineE pnbkgline8.LineE = pnbkgline7.LineE for l, s in zip(pnlines, pnlinewidth): l.Sigma = s for l in pnfree: l.Sigma.min = 1E-5 l.Sigma.max = 0.2 for l in pnfixwid: l.Sigma.freeze() for l, n in zip(pnlines, pnlinenorm): l.norm = n l.norm.min = 1E-10 l.norm.max = 1E10 # Scaling constant pnbkgcons.factor = 1.0 pnbkgcons.factor.freeze() # thermal radiation pnbkgapec.kT = 0.286928 pnbkgapec.kT.min = 0.008 pnbkgapec.kT.max = 64 pnbkgapec.Abundanc = 1.0 pnbkgapec.Abundanc.freeze() pnbkgapec.Redshift = 0.0 pnbkgapec.Redshift.freeze() pnbkgapec.norm = 5.58410E-05 pnbkgapec.norm.min = 1e-10 pnbkgapec.norm.max = 1e10 # local thermal radiation pnbkglcapec.kT = 0.1 pnbkglcapec.kT.freeze() pnbkglcapec.Abundanc = 1.0 pnbkglcapec.Abundanc.freeze() pnbkglcapec.Redshift = 0.0 pnbkglcapec.Redshift.freeze() pnbkglcapec.norm = 3.89164E-05 pnbkglcapec.norm.min = 1e-10 pnbkglcapec.norm.max = 1e10 # exponential decay pnbkgexpdec.factor = 44.3418 pnbkgexpdec.factor.min = 0 pnbkgexpdec.factor.max = 100 pnbkgexpdec.norm = 6830.89 pnbkgexpdec.norm.freeze() # Smear function pnbkgsmedge1.edgeE = 0.538408 pnbkgsmedge1.edgeE.freeze() pnbkgsmedge1.MaxTau = 1.40238 pnbkgsmedge1.MaxTau.min = 0 pnbkgsmedge1.MaxTau.max = 10 pnbkgsmedge1.index = -2.67000 pnbkgsmedge1.index.freeze() pnbkgsmedge1.width = 0.313365 pnbkgsmedge1.width.min = 0.01 pnbkgsmedge1.width.max = 100 pnbkgsmedge2.edgeE = 1.38826 pnbkgsmedge2.edgeE.freeze() pnbkgsmedge2.MaxTau.min = 0 pnbkgsmedge2.MaxTau.max = 10 pnbkgsmedge2.MaxTau = 9.37167 pnbkgsmedge2.index = -2.67000 pnbkgsmedge2.index.freeze() pnbkgsmedge2.width = 5.7642 pnbkgsmedge2.width.min = 0.01 pnbkgsmedge2.width.max = 100 # Spline funtion pnbkgspline1.Estart = 0.200000 pnbkgspline1.Estart.freeze() pnbkgspline1.Ystart = -1.31506 pnbkgspline1.Ystart.min = -1e+6 pnbkgspline1.Ystart.max = 1e+6 pnbkgspline1.Yend = 1064.16 pnbkgspline1.Yend.min = -1e+6 pnbkgspline1.Yend.max = 1e+6 pnbkgspline1.YPstart = -106.183 pnbkgspline1.YPstart.min = -1e+6 pnbkgspline1.YPstart.max = 1e+6 pnbkgspline1.YPend = -366.092 pnbkgspline1.YPend.min = -1e+6 pnbkgspline1.YPend.max = 1e6 pnbkgspline1.Eend = 1.74715 pnbkgspline1.Eend.min = 0 pnbkgspline1.Eend.max = 100 pnbkgspline2.Estart = 3.29056 pnbkgspline2.Ystart = 1.00643 pnbkgspline2.Ystart.min = -1e+6 pnbkgspline2.Ystart.max = 1e+6 pnbkgspline2.Yend = 0.887026 pnbkgspline2.Yend.min = -1e+6 pnbkgspline2.Yend.max = 1e+6 pnbkgspline2.YPstart = -0.278401 pnbkgspline2.YPstart.min = -1e+6 pnbkgspline2.YPstart.max = 1e+6 pnbkgspline2.YPend = 4.84809E-03 pnbkgspline2.YPend.min = -1e+6 pnbkgspline2.YPend.max = 1e6 pnbkgspline2.Eend = 7.32701 pnbkgspline2.Eend.min = 0 pnbkgspline2.Eend.max = 100 # Background power law pnbkginspl.PhoIndex = 0.279 pnbkginspl.PhoIndex.min = -2 pnbkginspl.PhoIndex.max = 9 pnbkginspl.norm = 8.23614E-03 pnbkginspl.norm.min = 1E-10 pnbkginspl.norm.max = 1E6 # Extragalactic background pnbkgpl.PhoIndex = 1.46 pnbkgpl.PhoIndex.freeze() pnbkgpl.norm = 1.25288E-05 pnbkgpl.norm.min = 1e-10 pnbkgpl.norm.max = 1e3 if(fit): #fit bkg set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl)))) print('Fitting (1/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline2)))) print('Fitting (2/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2)))) print('Fitting (3/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline7 + pnbkgline8)))) print('Fitting (4/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10)))) print('Fitting (5/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline3 + pnbkgline4 + pnbkgline5 + pnbkgline6 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10 + pnbkgline11)))) print('Fitting (6/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline3 + pnbkgline4 + pnbkgline5 + pnbkgline6 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10 + pnbkgline11))) + pnbrsp(galabs*(pnbkgapec))) print('Fitting (7/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline3 + pnbkgline4 + pnbkgline5 + pnbkgline6 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10 + pnbkgline11))) + pnbrsp(galabs*(pnbkgapec+pnbkgpl))) print('Fitting (8/8)...') fit_bkg(i) set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline3 + pnbkgline4 + pnbkgline5 + pnbkgline6 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10 + pnbkgline11))) + pnbrsp(galabs*(pnbkgapec+pnbkgpl)+pnbkglcapec)) fit_bkg(i) freeze(pnbkgcons, pnbkgspline1, pnbkgexpdec, pnbkgsmedge1, pnbkgsmedge2, pnbkgspline2, pnbkginspl, pnbkgline1, pnbkgline2, pnbkgline3, pnbkgline4, pnbkgline5, pnbkgline6, pnbkgline7, pnbkgline8, pnbkgline9, pnbkgline10, pnbkgline11, galabs, pnbkgapec, pnbkgpl, pnbkglcapec) print(' ') print('PN background model set up and fitted') print('Please double-check that it is a good fit') print(' ') else: set_bkg_full_model(i,pnbunitrsp(pnbkgcons*(pnbkgspline1*pnbkgexpdec + pnbkgsmedge1*pnbkgsmedge2*(pnbkgspline2*pnbkginspl + pnbkgline1 + pnbkgline2 + pnbkgline3 + pnbkgline4 + pnbkgline5 + pnbkgline6 + pnbkgline7 + pnbkgline8 + pnbkgline9 + pnbkgline10 + pnbkgline11))) + pnbrsp(galabs*(pnbkgapec+pnbkgpl)+pnbkglcapec)) print(' ') print('PN background model set up') print(' ') return pnscale*(pnbunitrsp(pncons*(pnspline1*pnexpdec + pnsmedge1*pnsmedge2*(pnspline2*pnbkpl + pngau1 + pngau2 + pngau3 + pngau4 + pngau5 + pngau6 + pngau7 + pngau8 + pngau9 + pngau10 + pngau11))) + pnbrsp(galabs*(pnapec+pnpnexpl)+pnlcapec))
[docs] def get_mos_bkg_model(i, galabs, fit=False): #================================================================= # Parameters: # i = sherpa ID of data set # galabs = name of model for galactic absorption # fit: True = fit bkg spectrum; False = just set_bkg_full_model # returns unit response for bkg model #================================================================= #get instrument response model for bkg spectrum if get_rmf(i).energ_lo[0] == 0: get_rmf(i).energ_lo[0] = 0.001 if get_arf(i).energ_lo[0] == 0: get_arf(i).energ_lo[0] = 0.001 mosbrsp = get_response(i,bkg_id=1) mosscale = get_bkg_scale(i) #get background scaling factor #create unit response dia_mos_rmf=get_embedded_file('mos_dia.rmf') dia_mos_arf=get_embedded_file('mos_dia.arf') copy_data(i,1000+2) load_bkg_rmf(1000+2, dia_mos_rmf) #load diagonal bkg matrices load_bkg_arf(1000+2, dia_mos_arf) mosbunitrsp = get_response(1000+2, bkg_id=1) delete_data(1000+2) #gaussian line center energy, line width, and initial normalization; for *MOS background* moscenters = [1.48600, 1.48700, 1.74000, 5.41000, 5.89500, 6.42000, 9.71000] moslinewidth = [3.84602E-02, 0.165816, 3.54985E-02, 9.77018E-02, 7.45076E-02, 7.42365E-02, 9.04855E-02] moslinenorm = [9.93119E-03, 1.67028E-03, 1.75461E-03, 2.86358E-04, 2.07525E-04, 3.07555E-04, 4.58115E-04] # model component prefix mos = 'mos%s' % i mosbkgcons, mosbkgsmedge, mosbkgspline, mosbkgbknpl, mosbkgline1, mosbkgline2, mosbkgline3, mosbkgline4, mosbkgline5, mosbkgline6, mosbkgline7, mosbkgpl, mosbkgapec, mosbkglcapec = (xsconstant(mos+'cons'), xssmedge(mos+'smedge'), xsspline(mos+'spline'), xsbknpower(mos+'bknpl'), xsgaussian(mos+'gau1'), xsgaussian(mos+'gau2'), xsgaussian(mos+'gau3'), xsgaussian(mos+'gau4'), xsgaussian(mos+'gau5'), xsgaussian(mos+'gau6'), xsgaussian(mos+'gau7'), xspowerlaw(mos+'expl'), xsapec(mos+'apec'), xsapec(mos+'lcapec')) moslines = [mosbkgline1, mosbkgline2, mosbkgline3, mosbkgline4, mosbkgline5, mosbkgline6, mosbkgline7] mos_bkg = 'mos_bkg_model' #define MOS background model mos_bkg = mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7) + mosbrsp(galabs*(mosbkgpl + mosbkgapec) + mosbkglcapec) for l, c in zip(moslines, moscenters): l.LineE = c l.LineE.min = c - 0.05 l.LineE.max = c + 0.05 for l, s in zip(moslines, moslinewidth): l.Sigma = s l.Sigma.min = 1E-4 l.Sigma.max = 0.2 mosbkgline1.Sigma.min = 1E-4 mosbkgline1.Sigma.max = 1E-1 for l, n in zip(moslines, moslinenorm): l.norm = n l.norm.min = 1E-10 l.norm.max = 1E10 # Constant factor mosbkgcons.factor = 1.0 mosbkgcons.factor.freeze() # Smear function mosbkgsmedge.edgeE = 0.538408 mosbkgsmedge.edgeE.freeze() mosbkgsmedge.MaxTau = 0.246633 mosbkgsmedge.MaxTau.min = 0.0 mosbkgsmedge.MaxTau.max = 10.0 mosbkgsmedge.index = -2.67 mosbkgsmedge.index.freeze() mosbkgsmedge.width = 1E-02 mosbkgsmedge.width.min = 1E-02 mosbkgsmedge.width.max = 1E2 # Spline funtion mosbkgspline.Estart = 3.08175 mosbkgspline.Ystart = 1.00984 mosbkgspline.Yend = 1.99144 mosbkgspline.YPstart = -2.90195e-02 mosbkgspline.YPend = 5.49102e-02 mosbkgspline.Estart.freeze() mosbkgspline.Ystart.freeze() mosbkgspline.Yend.freeze() mosbkgspline.YPstart.freeze() mosbkgspline.YPend.freeze() mosbkgspline.Eend = 13.6492 mosbkgspline.Eend.min = 0 mosbkgspline.Eend.max = 100 # Broken power law mosbkgbknpl.PhoIndx1 = 1.48636 mosbkgbknpl.PhoIndx1.min = -2 mosbkgbknpl.PhoIndx1.max = 9 mosbkgbknpl.BreakE = 0.415173 mosbkgbknpl.PhoIndx2 = 0.315615 mosbkgbknpl.BreakE.freeze() mosbkgbknpl.PhoIndx2.freeze() mosbkgbknpl.norm = 2.90071E-03 mosbkgbknpl.norm.min = 1E-10 mosbkgbknpl.norm.max = 1E10 # Extragalactic background mosbkgpl.PhoIndex = 1.46 mosbkgpl.PhoIndex.freeze() mosbkgpl.norm = 4.58115E-04 mosbkgpl.norm.min = 1E-10 mosbkgpl.norm.max = 1E10 # thermal radiation mosbkgapec.kT = 0.286928 mosbkgapec.kT.min = 0.008 mosbkgapec.kT.max = 64 mosbkgapec.Abundanc = 1.0 mosbkgapec.Abundanc.freeze() mosbkgapec.Redshift = 0.0 mosbkgapec.Redshift.freeze() mosbkgapec.norm = 5.58410E-05 mosbkgapec.norm.min = 1e-10 mosbkgapec.norm.max = 1e10 # local thermal radiation mosbkglcapec.kT = 0.1 mosbkglcapec.kT.freeze() mosbkglcapec.Abundanc = 1.0 mosbkglcapec.Abundanc.freeze() mosbkglcapec.Redshift = 0.0 mosbkglcapec.Redshift.freeze() mosbkglcapec.norm = 3.89164E-05 mosbkglcapec.norm.min = 1e-10 mosbkglcapec.norm.max = 1e10 if(fit): #fit bkg set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl)) print('Fitting (1/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline2)) print('Fitting (2/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline2 + mosbkgline3)) print('Fitting (3/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3)) print('Fitting (4/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7)) print('Fitting (5/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7) + mosbrsp(galabs*(mosbkgapec))) print('Fitting (6/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7) + mosbrsp(galabs*(mosbkgapec+mosbkgpl))) print('Fitting (7/8)...') fit_bkg(i) set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7) + mosbrsp(galabs*(mosbkgapec+mosbkgpl)+mosbkglcapec)) print('Fitting (8/8)...') fit_bkg(i) freeze(mosbkgcons, mosbkgsmedge, mosbkgspline, mosbkgbknpl, mosbkgline1, mosbkgline2, mosbkgline3, mosbkgline4, mosbkgline5, mosbkgline6, mosbkgline7, mosbkgpl, mosbkgapec, mosbkglcapec) print(' ') print('MOS background model set up and fitted') print('Please double-check that it is a good fit') print(' ') else: set_bkg_full_model(i,mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7) + mosbrsp(galabs*(mosbkgapec+mosbkgpl)+mosbkglcapec)) print(' ') print('MOS background model set up') print(' ') return mosscale*(mosbunitrsp(mosbkgcons*mosbkgsmedge*mosbkgspline*mosbkgbknpl + mosbkgline1 + mosbkgline2 + mosbkgline3 + mosbkgline4 + mosbkgline5 + mosbkgline6 + mosbkgline7) + mosbrsp(galabs*(mosbkgapec+mosbkgpl)+mosbkglcapec))
#return mosscale*(mosbunitrsp(moscons*mossmedge*mosspline*mosbknpl + mosgau1 + mosgau2 + mosgau3 + # mosgau4 + mosgau5 + mosgau6 + mosgau7) + mosbrsp(galabs*(mosapec+mosexpl)+moslcapec))
[docs] def get_mos_bkg_model_cached(i, galabs): filename = get_bkg(i).name + '.bkgpars' if os.path.exists(filename): bkgmodel = get_mos_bkg_model(i, galabs, fit=False) for p, v in zip(bkgmodel.pars, numpy.loadtxt(filename)): p.val = v else: bkgmodel = get_mos_bkg_model(i, galabs, fit=True) numpy.savetxt(filename, [p.val for p in bkgmodel.pars]) for p in bkgmodel.pars: p.freeze() return bkgmodel
[docs] def get_pn_bkg_model_cached(i, galabs): filename = get_bkg(i).name + '.bkgpars' if os.path.exists(filename): bkgmodel = get_pn_bkg_model(i, galabs, fit=False) for p, v in zip(bkgmodel.pars, numpy.loadtxt(filename)): p.val = v else: bkgmodel = get_pn_bkg_model(i, galabs, fit=True) numpy.savetxt(filename, [p.val for p in bkgmodel.pars]) for p in bkgmodel.pars: p.freeze() return bkgmodel