#!/usr/bin/env python3
#
'''
.. module:: specfit
:synopsis: fit lines in (MXS) eventlist
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
Fit a multiline spectrum. Fit energy resolutions and liniarity.
'''
import numpy
from tesfdmtools.methods.holzer import holzerfit,ehe1,ehe2,ehi1,ehi2,ebe1,ebi1,holzer,hbeta,holzerelements
from tesfdmtools.utils.fiterrplot import fiterrplot
from tesfdmtools.methods.MCeventlist import plotspec
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
# compute average line positions
alphapos={}
betapos={}
for el in ehe1:
alphapos[el]=(numpy.sum(ehe1[el]*ehi1[el])+numpy.sum(ehe2[el]*ehi2[el]))/\
(numpy.sum(ehi1[el])+numpy.sum(ehi2[el]))
if ebi1[el][0] != 0.0:
betapos[el] =numpy.sum(ebe1[el]*ebi1[el])/numpy.sum(ebi1[el])
#print(alphapos)
#print(betapos)
def erptxt(fpar):
txt=' '.join((' shift: %8.2f (eV)\n' % fpar[0]).split())+'\n'
txt=txt+' '.join((' fwhm: %8.1f (eV)\n' % (fpar[3]*1000.0)).split())+'\n'
txt=txt+' '.join(('C-stat: %8.1f\n' % fpar[4]).split())+'\n'
txt=txt+' '.join(('deg.frdm: %8.1f' % fpar[5]).split())
return txt
[docs]def specfit(elist,elements=['cr','cu'],fitrange=[-30.0,30.0],debug=False,pltax=None):
'''
Fit multi-line spectrum. Specify the elements for which lines have to be fitted.
The first two elements in the list must specify the two strongest alpha lines.
At this moment this routine assumes fluorescent lines only. No Bremsstrahlung continuum.
Args:
* `elist` = float array with the list of event energies
**Kwargs:
* `elements` = list of elements which have lines in the spectrum. The first two elements must have the strongest alpha lines.
* `fitrange` = range around the line centers to use for the fit.
* `debug` = logical, when **True**, show some debug information, including plot output on screen.
* `pltax` = plot instance (**Figure**) for plot output. If **None** show output on screen.
'''
# define elements of interest for fit
holzerelements(elements)
# first make course spectrum to determine approximate line positions
lpos=[] # energy range to consider
for el in elements:
lpos.append(alphapos[el])
lpos.append(betapos[el])
lrange=[min(lpos),max(lpos)]
delta=0.1*(lrange[1]-lrange[0])
erange=numpy.array([lrange[0]-delta,lrange[1]+delta])
sbins=int((erange[1]-erange[0])/30.0+1.0) # binning for course spectrum
lrange=[numpy.amin(elist),numpy.amax(elist)]
dl=0.1*(lrange[1]-lrange[0])
hrange=[lrange[0]-dl,lrange[1]+dl]
cspectrum,binedges=numpy.histogram(elist,bins=sbins,range=hrange)
cbinpos=(binedges[:-1]+binedges[1:])/2.0
if debug:
pltax=None
plotspec(cbinpos,cspectrum)
# find the peaks
peakpos=numpy.zeros(100)
peakint=numpy.zeros(100)
cspectrum[0]=5
cspectrum[-1]=5
pindx,=numpy.where( cspectrum > (2.0*numpy.average(cspectrum)) ) # find peaks
gpos=pindx[1:]-pindx[0:-1] #
gpos[0]=3
gpos[-1]=3
ppos,=numpy.where( gpos > 2 ) # gaps
k=0
for j,ii in enumerate(ppos[0:-1]):
i1=pindx[ii+1]-1
i2=pindx[ppos[j+1]]+1
peakint[k]=numpy.sum(cspectrum[i1:i2])
peakpos[k]=numpy.sum(numpy.arange(i1,i2,dtype=float)*cspectrum[i1:i2])/\
float(numpy.sum(cspectrum[i1:i2]))
k=k+1
sindx=numpy.flip(numpy.argsort(peakint[:k]),axis=0)
peakpos=cbinpos[(peakpos[sindx]+0.5).astype(int)] # peak positions in elist
peakint=peakint[sindx] # peak intensities
if debug:
print(peakpos,'\n',peakint)
# compute approximate energy scale, based on two highest peaks corresponding to alpha lines of first
# two elements.
minl=9E9
maxl=-9E9
for el in elements[0:2]:
minl=min([minl,alphapos[el]])
maxl=max([maxl,alphapos[el]])
mine=peakpos[0:2].min()
maxe=peakpos[0:2].max()
apreng=(cbinpos-mine)/(maxe-mine)*(maxl-minl)+minl # approximate energy scale
bescale=(maxe-mine)/(maxl-minl)
if debug:
plotspec(apreng,cspectrum)
#
# make fine spectrum binning and fit each line separately within its fitranges
#
mins=9E9
maxs=-9E9
for el in elements:
mins=min([mins,alphapos[el],betapos[el]])
maxs=max([maxs,alphapos[el],betapos[el]])
mins=mins+2.0*fitrange[0] # energy range for spectrum binning
maxs=maxs+2.0*fitrange[1]
if debug:
print('fine spectrum range: ',mins,maxs)
minb=(mins-minl)*bescale+mine
maxb=(maxs-minl)*bescale+mine
fnbins=int((maxs-mins)/0.2+0.5)
if debug:
print('histogram: ',minb,maxb,fnbins)
fspectrum,binedges=numpy.histogram(elist,bins=fnbins,range=[minb,maxb]) # fine binning spectrum
fbins=(binedges[:-1]+binedges[1:])/2.0
fenergies=(fbins-mine)/(maxe-mine)*(maxl-minl)+minl
if debug:
plotspec(fbins,fspectrum)
nfits=len(elements)*2
shifts=numpy.zeros(nfits,dtype=float) # fitted line shifts
efits=numpy.zeros(nfits,dtype=float) # average line positions
res=numpy.zeros(nfits,dtype=float) # resolutions
reserr=numpy.zeros(nfits,dtype=float) # resolution errors
gain=numpy.zeros(nfits,dtype=float) # gain
gainerr=numpy.zeros(nfits,dtype=float) # gain errors
elmnts=numpy.repeat(' ',nfits) # elements
ab=numpy.zeros(nfits,dtype=int)
albe=['$_\\alpha$','$_\\beta$']
fsize=None
if pltax is not None:
gs = gridspec.GridSpec(nrows=(len(elements)+1), ncols=2)
fsize=6
for niter in [1,2]: # iterate 2 times
if debug:
print("=============== iter %d =============================" % niter)
n=0
for je,el in enumerate(elements): # fit each line
mina=alphapos[el]+fitrange[0]-shifts[je*2]
maxa=alphapos[el]+fitrange[1]-shifts[je*2]
minb=(mina-minl)*bescale+mine
maxb=(maxa-minl)*bescale+mine
if debug:
print("line plot: ",minb,maxb)
indx,=numpy.where( ( fbins > minb ) & ( fbins < maxb ) )
fpar=holzerfit(fenergies[indx]/1000.0,fspectrum[indx],line='Alpha',elements=[el])
yfit=holzer(fpar,fenergies[indx]/1000.0)
if debug or ( niter == 2 ):
if pltax is not None:
eax=pltax.add_subplot(gs[je, 0])
else:
eax=None
fiterrplot(fenergies[indx],fspectrum[indx],numpy.sqrt(fspectrum[indx]),yfit,\
xtitle='Energy (eV)',ytitle='Counts',ptitle=(el+" (Alpha)"),\
pltax=eax,ptxt=erptxt(fpar),ptxtsize=fsize)
# plotspec(fenergies[indx],fspectrum[indx])
shifts[n]=fpar[0]
res[n]=fpar[3]*1000.0
reserr[n]=fpar[9]*1000.0
gain[n]=fpar[1]
gainerr[n]=fpar[7]
efits[n]=alphapos[el]
elmnts[n]=el
ab[n]=0
n=n+1
mina=betapos[el]+fitrange[0]-shifts[je*2+1]
maxa=betapos[el]+fitrange[1]-shifts[je*2+1]
minb=(mina-minl)*bescale+mine
maxb=(maxa-minl)*bescale+mine
if debug:
print("line plot: ",minb,maxb)
indx,=numpy.where( ( fbins > minb ) & ( fbins < maxb ) )
fpar=holzerfit(fenergies[indx]/1000.0,fspectrum[indx],line='Beta',elements=[el])
yfit=hbeta(fpar,fenergies[indx]/1000.0)
if debug or ( niter == 2 ):
if pltax is not None:
eax=pltax.add_subplot(gs[je, 1])
else:
eax=None
fiterrplot(fenergies[indx],fspectrum[indx],numpy.sqrt(fspectrum[indx]),yfit,\
xtitle='Energy (eV)',ytitle='Counts',ptitle=(el+" (Beta)"),\
pltax=eax,ptxt=erptxt(fpar),ptxtsize=fsize)
# plotspec(fenergies[indx],fspectrum[indx])
shifts[n]=fpar[0]
res[n]=fpar[3]*1000.0
reserr[n]=fpar[9]*1000.0
gain[n]=fpar[1]
gainerr[n]=fpar[7]
efits[n]=betapos[el]
elmnts[n]=el
ab[n]=1
n=n+1
fline=numpy.polyfit(efits,shifts,1) # fit energy shift systematics
ss=fline[0]*fenergies+fline[1]
fenergies=fenergies+ss
if debug: # show shift systematics
f,ax=plt.subplots(1)
ax.plot(efits,shifts,'bo')
ax.plot(fenergies,ss,'r-')
for k,eei in enumerate(efits):
ax.text(eei,shifts[k],' '+elmnts[k]+albe[ab[k]]+'\n ')
ax.set_xlabel('Energy')
ax.set_ylabel('shift')
ax.set_title('fit systematics')
plt.show()
plt.close('all')
if pltax is None:
f,ax=plt.subplots(2,sharex=True) # show fitted resolutions and gains
ax0=ax[0]
ax1=ax[1]
else:
ax0=pltax.add_subplot(gs[len(elements), 0])
ax1=pltax.add_subplot(gs[len(elements), 1])
ss=numpy.argsort(efits)
ax0.plot(efits[ss],res[ss],'ro')
ax0.errorbar(efits[ss],res[ss],yerr=reserr,color='b')
ax1.plot(efits[ss],gain[ss],'ro')
ax1.errorbar(efits[ss],gain[ss],yerr=gainerr,color='b')
for sss in ss:
ax0.text(efits[sss],res[sss],' '+elmnts[sss]+albe[ab[sss]]+'\n ')
ax1.text(efits[sss],gain[sss],' '+elmnts[sss]+albe[ab[sss]]+'\n ')
ax0.set_xlabel('Energy (eV)')
ax0.set_ylabel('Resolution (eV)')
ax1.set_ylabel('gain')
ax0.set_title('Specfits results')
if pltax is None:
plt.show()
plt.close('all')
else:
plt.tight_layout()
#=============================================================================
if __name__ == "__main__":
# test program
import sys
from tesfdmtools.methods.MCeventlist import MCeventlist
# from matplotlib.backends.backend_pdf import PdfPages
elements=['mn','cu','cr']
lines=['Alpha','Beta']
bspec=None
if len(sys.argv) > 1:
elements=[]
lines=[]
for pp in sys.argv[1:]:
if pp in ['Alpha','Beta']:
lines.append(pp)
elif pp in ['mn','cr','cu','ti']:
elements.append(pp)
elif pp == 'bspec':
bspec=True
else:
sys.exit("Use: %s element1 [ element2 ....] line1 [ line2 ]" % sys.argv[0] )
elist,ss,ee=MCeventlist(counts=5000,elements=elements,lines=lines,res=3.0,bspec=bspec)
av=numpy.average(elist)
elist=elist/av
# pdf = PdfPages('specfit.pdf')
fig=plt.figure(figsize=tuple(numpy.array([8.27,11.69])*0.9))
specfit(elist,elements=elements,fitrange=[-30.0,30.0],debug=False,pltax=fig)
# plt.savefig(pdf,format='pdf',papertype='a4')
# pdf.close()
plt.show()
plt.close('all')