#!/usr/bin/env python3
#
'''
.. module:: holzer
:synopsis: Compute and fit K-alpha/K-beta lines according to Holzer
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
Compute and fit K-alpha/K-beta lines.
Original line data according to the profiles specified by
Holzer et al, Physical Review A, Vol. 56, No. 6, p. 4554, December 1997
Ti line from: Y.Ito et al, Physical Review A Vol. 94, p. 042506 (2016)
Mn alpha lines from GSFC document "Instrument Calibration report Natural Line Shapes"
authors M.Eckart,C. Kilbourne, Scott Porter
https://heasarc.gsfc.nasa.gov/docs/hitomi/calib/caldb_doc/asth_sxs_caldb_linefit_v20161223.pdf
'''
import numpy,scipy.optimize
import sys
# parameter dictionaries for the different elements
ehe1={}
ehe2={}
ehw1={}
ehw2={}
ehi1={}
ehi2={}
ebe1={}
ebw1={}
ebi1={}
# manganese alpha
ehe1['mn']=numpy.array([5898.882,5897.898,5894.864,5896.566,5899.444,5902.712])
ehe2['mn']=numpy.array([5887.772,5886.528])
ehw1['mn']=numpy.array([1.7145,2.0442,4.4985,2.6616,0.97669,1.5528])/2.
ehw2['mn']=numpy.array([2.3604,4.2168])/2.
ehi1['mn']=numpy.array([.784,.263,.067,.095,.071,0.011])
ehi2['mn']=numpy.array([.369,.10])
# manganese beta
ebe1['mn']=numpy.array([6490.89,6486.31,6477.73,6490.06,6488.83])
ebw1['mn']=numpy.array([1.83,9.40,13.22,1.81,2.81])/2.
ebi1['mn']=numpy.array([0.608,0.109,0.077,0.397,0.176])
# chrome alpha
ehe1['cr']=numpy.array([5414.874,5414.099,5412.745,5410.583,5418.304])
ehe2['cr']=numpy.array([5405.551,5403.986])
ehw1['cr']=numpy.array([1.457,1.760,3.138,5.149,1.988])/2.
ehw2['cr']=numpy.array([2.224,4.740])/2.
ehi1['cr']=numpy.array([0.822,0.237,0.085,0.045,0.015])
ehi2['cr']=numpy.array([0.386,0.036])
# chrome beta
ebe1['cr']=numpy.array([5947.00,5935.31,5946.24,5942.04,5944.93])
ebw1['cr']=numpy.array([1.70,15.98,1.90,6.69,3.37])/2.
ebi1['cr']=numpy.array([0.670,0.055,0.337,0.082,0.151])
# copper alpha
ehe1['cu']=numpy.array([8047.837,8045.367])
ehe2['cu']=numpy.array([8027.993,8026.504])
ehw1['cu']=numpy.array([2.285,3.358])/2.
ehw2['cu']=numpy.array([2.666,3.571])/2.
ehi1['cu']=numpy.array([0.957,0.090])
ehi2['cu']=numpy.array([0.334,0.111])
# chrome beta
ebe1['cu']=numpy.array([8905.532,8903.109,8908.462,8897.387,8911.393])
ebw1['cu']=numpy.array([3.52,3.52,3.55,8.08,5.31])/2.
ebi1['cu']=numpy.array([0.757,0.388,0.171,0.068,0.055])
# titanium alpha
ehe1['ti']=numpy.array([4510.832,4509.604,4514.029])
ehe2['ti']=numpy.array([4504.701,4505.66 ])
ehw1['ti']=numpy.array([1.393,2.78,1.538])/2.
ehw2['ti']=numpy.array([1.757,5.12])/2.
ehi1['ti']=numpy.array([100.0,25.39,3.364])/150.0
ehi2['ti']=numpy.array([49.7,25.27])/150.0
# tittanium beta (no data)
ebe1['ti']=numpy.array([4510.0])
ebw1['ti']=numpy.array([0.01])/2.
ebi1['ti']=numpy.array([0.0])
class Hpar:
def __init__(self,elements=['mn']):
for j,com in enumerate(elements):
if com not in ehe1:
sys.exit(("Error, element: '%s' is not a valid Holzer element\n" % com)+
"Valid elements are: '"+"' '".join(ehe1.keys())+"'")
if j == 0:
self.he1=ehe1[com]
self.he2=ehe2[com]
self.hw1=ehw1[com]
self.hw2=ehw2[com]
self.hi1=ehi1[com]
self.hi2=ehi2[com]
self.be1=ebe1[com]
self.bw1=ebw1[com]
self.bi1=ebi1[com]
else:
self.he1=numpy.append(self.he1,ehe1[com])
self.he2=numpy.append(self.he2,ehe2[com])
self.hw1=numpy.append(self.hw1,ehw1[com])
self.hw2=numpy.append(self.hw2,ehw2[com])
self.hi1=numpy.append(self.hi1,ehi1[com])
self.hi2=numpy.append(self.hi2,ehi2[com])
self.be1=numpy.append(self.be1,ebe1[com])
self.bw1=numpy.append(self.bw1,ebw1[com])
self.bi1=numpy.append(self.bi1,ebi1[com])
self.hapos=(numpy.append(self.he1,self.he2)
*numpy.append(self.hi1,self.hi2)).sum()/\
numpy.append(self.hi1,self.hi2).sum()
self.hd1=self.hapos-self.he1
self.hd2=self.hapos-self.he2
if self.be1.size == 1:
self.bapos=self.be1[0]
else:
self.bapos=(self.be1*self.bi1).sum()/self.bi1.sum()
self.bd1=self.bapos-self.be1
self.elements=elements
self.range={'Alpha':[self.he2.min(),self.he1.max()],
'Beta':[self.be1.min(),self.be1.max()]}
hels=Hpar()
# print "Average position Alpha, Beta: ",hapos,bapos
ptxt=[" shift"," gain"," amplitude","resolution"]
def extendx(x):
dx=x[1]-x[0] # delta x
ix=x.size//2 # number of bins to shift new x-axis
jx=ix+x.size
nx=x.size*2
x0=x[0]-ix*dx # start of new x-axisxx
xx=x0+numpy.arange(nx)*dx # new x-axis
return (xx,ix,jx) # return new x with pointers where old x-axis is included
def convolve(ff,xx,w):
gconv=numpy.exp(-(xx-xx[0])**2/(2.0*((w*1000.0)/2.35482)**2)) # gaussian for convolution
gconv[1:]=gconv[1:]+gconv[-1:0:-1] # make symmetric
gconv=gconv/gconv.sum() # normalize
cfunc=abs(numpy.fft.ifft(numpy.fft.fft(ff)*numpy.fft.fft(gconv))) # perform convolution
return cfunc
boffset=0.0 # bremsstrahlungs level offset for holzer fits
[docs]def holzer(p,x,offset=None):
'''
Generate Holzer profile for Holzer-elements Alpha lines.
Args:
* `p` = parameters for the profile to generate
* p[0]: shift
* p[1]: compression factor (gain)
* p[2]: normalization
* p[3]: gaussian convolution width (keV)
* `x` = list of energies (keV) for which to generate the profile
Kwargs:
* `offset` = brehmsstrahlungs offset
Example of use::
profile = holzer(p,x)
'''
global boffset # bremsstrahlungs level offset for holzer fits
if offset is None:
offset=boffset
xk=x*1000.0 # keV ==> eV
xx,x1,x2=extendx(xk) # extend x-axis
tot=numpy.zeros(len(xx),dtype='float')
for j in range(len(hels.hd1)):
tot=tot+hels.hi1[j]/(1.0+
((xx-(hels.hapos-p[0]-p[1]*hels.hd1[j]))/(p[1]*hels.hw1[j]))**2) # sum of lorentzians for alpha-1
for j in range(len(hels.hd2)):
tot=tot+hels.hi2[j]/(1.0+
((xx-(hels.hapos-p[0]-p[1]*hels.hd2[j]))/(p[1]*hels.hw2[j]))**2) # sum of lorentzians for alpha-2
tot=tot/tot.sum()*p[2] # intensity
if p[3] <= 0.0:
return tot[x1:x2]+offset
cc=convolve(tot,xx,p[3])
return cc[x1:x2]+offset
[docs]def hbeta(p,x,offset=None):
'''
Generate Holzer profile for Holzer-elements Beta lines.
Args:
* `p` = parameters for the profile to generate
* p[0]: shift
* p[1]: compression factor (gain)
* p[2]: normalization
* p[3]: gaussian convolution width (keV)
* `x` = list of energies (keV) for which to generate the profile
Kwargs:
* `offset` = bremsstrahlungs offset
Example of use::
profile = hbeta(p,x)
'''
global boffset
if offset is None:
offset=boffset
xk=x*1000.0 # keV ==> eV
xx,x1,x2=extendx(xk) # extend x-axis
tot=numpy.zeros(len(xx),dtype='float')
for j in range(len(hels.bd1)):
tot=tot+hels.bi1[j]/(1.0+
((xx-(hels.bapos-p[0]-p[1]*hels.bd1[j]))/(p[1]*hels.bw1[j]))**2) # sum of lorentzians for beta
tot=tot/tot.sum()*p[2] # intensity
# if len(p) > 4: # add baseline level as free parameter ?
# tot=tot+p[4]
if p[3] <= 0.0:
return tot[x1:x2]+offset
cc=convolve(tot,xx,p[3])
return cc[x1:x2]+offset
[docs]def holzerelements(elements=['mn']):
'''
Define Holzer elements to be used for holzer spectrum generation.
Valid elements are: 'mn','cu','cr','ti'
Kwargs:
* `elements` = tuple of elements generating lines
Returns:
* range = range in energy of all input elements (for both Alpha and Beta lines)
'''
global hels
hels=Hpar(elements=elements)
return hels.range
[docs]def holzerfit(x,y,line='Alpha',elements=['mn'],offset=0.0,bgain=None):
'''
Fit holzer profile to given count spectrum, using C-stat statistics
Args:
* `x` = list of energies (keV)
* `y` = list of counts (counts)
Kwargs:
* `line` = Type of line to fit 'Alpha' or 'Beta'
* `elements` = tuple of elements generating lines
* `offset` = Bremsstrahlungs level offset
* `bgain` = Fix Beta line gain parameter. When **None**, the gain is a free parameter.
Returns:
* `fitpar` = list of fit parameters
* fitpar[0]: shift (eV)
* fitpar[1]: compression factor (gain)
* fitpar[2]: normalization
* fitpar[3]: Gaussian resolution width (eV)
* fitpar[4]: C-stat
* fitpar[5]: degrees of freedom
* fitpar[6]: 1-sigma error in fitpar[0]
* fitpar[7]: 1-sigma error in fitpar[1]
* fitpar[8]: 1-sigma error in fitpar[2]
* fitpar[9]: 1-sigma error in fitpar[3]
Example of use::
fitpar = holzerfit(x,y)
'''
global ptxt,hels,boffset
hels=Hpar(elements=elements)
boffset=offset
if line == 'Beta':
if bgain is None:
hfunc=lambda x, p1, p2, p3, p4: hbeta((p1,p2,p3,p4),x)
else:
hfunc=lambda x, p1, p3, p4: hbeta((p1,bgain,p3,p4),x)
shift=hels.be1[0]-x[y.argmax()]*1000.0
# shift=be1[0]-(x[0]+x[-1])/2.0*1000.0
else:
hfunc=lambda x, p1, p2, p3, p4: holzer((p1,p2,p3,p4),x)
shift=hels.he1[1]-x[y.argmax()]*1000.0
# shift=(he1[1]+he2[0])/2.0-(x[0]+x[-1])/2.0*1000.0
norm=numpy.sum(y)
maxy=max(y)
top,=numpy.where(y > (0.5*float(maxy)))
fwhm=top.size*(x[1]-x[0]) # first est. fwhm (0.003 prev)
p=numpy.array([shift,1.0,norm,fwhm])
covar=numpy.zeros(p.size,dtype=float) # use in case of failed fit
# p=numpy.append(p,0.0) # baseline level 0.0
# print "p: ",p
fitok=True
bounds=([-numpy.inf,0.,0.,0.],numpy.inf) # parameter bounds
yy=numpy.where(y <= 0.0,1.0,y)
if ( bgain is not None ) and ( line == 'Beta' ):
p=numpy.array([p[0],p[2],p[3]],dtype=float) # ignore gain parameter for beta line
bounds=([-numpy.inf,0.,0.],numpy.inf)
try:
p1, covar = scipy.optimize.curve_fit(hfunc, x, y, p.copy(), numpy.sqrt(yy) ,\
bounds=bounds)
except:
p1=p.copy()
fitok=False
ps=p1.copy()
#
# recompute weights, based on fit and iterate (Wheaton et. al. 1995)
# needed for low count rate statistics (see e.g. SPEX manual 2.12 "spectra fitting" )
#
if line == 'Beta':
if bgain is not None:
ye=numpy.sqrt(hbeta(numpy.array([p1[0],bgain,p1[1],p1[2]]),x))
else:
ye=numpy.sqrt(hbeta(p1,x))
else:
ye=numpy.sqrt(holzer(p1,x))
try:
p1, covar = scipy.optimize.curve_fit(hfunc, x, y, p1.copy(), ye ,\
bounds=bounds,absolute_sigma=True)
except:
p1=ps.copy()
fitok=False
ps=p1.copy()
if line == 'Beta':
if bgain is not None:
ye=numpy.sqrt(hbeta(numpy.array([p1[0],bgain,p1[1],p1[2]]),x))
else:
ye=numpy.sqrt(hbeta(p1,x))
else:
ye=numpy.sqrt(holzer(p1,x))
try:
p1, covar = scipy.optimize.curve_fit(hfunc, x, y, p1.copy(), ye ,\
bounds=bounds,absolute_sigma=True)
except:
p1=ps.copy()
fitok=False
ps=p1.copy()
perr = numpy.sqrt(numpy.diag(covar)) # parameter error estimates
if ( bgain is not None ) and ( line == 'Beta'):
p1=numpy.array([p1[0],bgain,p1[1],p1[2]],dtype=float)
perr=numpy.array([perr[0],0.0,perr[1],perr[2]],dtype=float)
#
if not fitok:
p1[0]=0.0
p1[1]=1.00
p1[3]=0.0
print("Error, holzer fit failed to converge")
perr=numpy.zeros(p1.size)
print("holzerfit result: ")
for i in range(len(ptxt)):
print(" %s %10.3g +/- %10.3g" % (ptxt[i],p1[i],perr[i]))
if line == 'Beta':
yf=hbeta(p1,x)
else:
yf=holzer(p1,x)
c=2.0*numpy.sum(yf-y+y*numpy.log(yy/yf))
deg=x.size+p1.size
print(" C-stat: %7.1f ( %d degrees of freedom )" % (c,deg))
#
boffset=0.0
#
return numpy.append(p1,numpy.append([c,deg],perr))
#==============================================================
if __name__ == "__main__":
'''
Test program
'''
import sys
import matplotlib.pyplot as plt
if len(sys.argv) > 1:
elements=sys.argv[1:]
else:
elements=['mn']
tpar={}
for el in ehe1:
da=(ehe1[el].max()-ehe2[el].min())*0.8
db=(ebe1[el].max()-ebe1[el].min())*0.8
tpar[el]=[(ehe2[el].min()-da)/1000.0,(ehe1[el].max()+da)/1000.0,
numpy.average(ehe2[el])/1000.0,numpy.average(ehe1[el])/1000.0,0.0022,
(ebe1[el].min()-db)/1000.0,(ebe1[el].max()+db)/1000.0,
numpy.average(ebe1[el])/1000.0,100.0,220.0,100.0]
p=tpar[elements[0]]
print(p)
xx=numpy.arange(p[0],p[1],0.0002) # energy in keV !!
yy=p[8]*numpy.exp(-(xx-p[2])**2/(2.0*p[4]**2))+\
p[9]*numpy.exp(-(xx-p[3])**2/(2.0*p[4]**2))+\
20.0*numpy.random.rand(len(xx))
pp=holzerfit(xx,yy,elements=[elements[0]])
xb=numpy.arange(p[5],p[6],0.0002) # energy in keV !!
yb=p[10]*numpy.exp(-(xb-p[7])**2/(2.0*p[4]**2))+\
20.0*numpy.random.rand(len(xb))
if xb.size > 0:
bb=holzerfit(xb,yb,line='Beta',elements=[elements[0]])
f,ax=plt.subplots(2)
ax[0].set_xlabel("X")
ax[0].set_ylabel("Y")
ax[0].set_title("testing Holzer fit: '"+elements[0]+"'")
ax[0].plot(xx,yy)
ax[0].plot(xx,holzer(pp,xx))
print("Fitted parameters: ",pp)
print("FWHM: ",pp[3])
if xb.size > 0:
ax[1].set_xlabel("X")
ax[1].set_ylabel("Y")
ax[1].plot(xb,yb)
ax[1].plot(xb,hbeta(bb,xb))
print("Fitted parameters: ",bb)
print("FWHM: ",bb[3])
plt.tight_layout()
plt.show()
plt.close('all')