#!/usr/bin/env python3
'''
.. module:: mxsspec
:synopsis: utilities to interpret the mxs(Cr/Cu)+Mn55 spectrum
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import os.path
import sys
import numpy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.widgets import Cursor
from mpl_toolkits.axes_grid1 import make_axes_locatable
from tesfdmtools.utils.fiterrplot import fiterrplot
from tesfdmtools.methods.expfit import expfit,expfunc
from tesfdmtools.methods.gaussfit import gaussfit,gauss
from tesfdmtools.methods.fitcor import fitcor
from tesfdmtools.methods.xcorplot import xcorplot
from tesfdmtools.methods.HDF_Eventutils import makeplot
# MXS and Mn55 line positions
from tesfdmtools.methods.holzer import ehe1,ehe2,ehi1,ehi2,ebe1,ebi1,holzerfit,holzer,hbeta
epos=numpy.zeros(6,dtype=float)
ltype=numpy.chararray(6) # line type 'alpha' 'beta'
etype=numpy.empty(6,dtype='object') # element type (Mn,Cr,Cu)
abpos={}
for j,el in enumerate(['mn','cu','cr']):
jj=j*2
epos[jj]=(numpy.sum(ehe1[el]*ehi1[el])+numpy.sum(ehe2[el]*ehi2[el]))/\
(numpy.sum(ehi1[el])+numpy.sum(ehi2[el])) # average alpha line position
ltype[jj]='a'
etype[jj]=el
jj=jj+1
epos[jj]=numpy.sum(ebe1[el]*ebi1[el])/numpy.sum(ebi1[el]) # average beta line position
ltype[jj]='b'
etype[jj]=el
abpos[el]=[epos[jj-1],epos[jj]]
epindx=epos.argsort() # list of sorted line positions
epos=epos[epindx] # list of sorted line positions
ltype=ltype[epindx] # sorted list of line types
etype=etype[epindx]
alpha,=numpy.where( ltype == b'a' ) # index to sorted alpha lines
beta,=numpy.where( ltype == b'b' ) # index to sorted beta lines
eparatio=numpy.zeros(alpha.size,dtype=float) # ratio's within alpha line positions
epalids=[(1,0,2),(2,0,1),(2,1,0)] # two line-id's for ratios and other line
for j in numpy.arange(alpha.size):
eparatio[j]=epos[alpha[epalids[j][0]]]/epos[alpha[epalids[j][1]]]
#print("eparatio: ",eparatio)
# ================================================================
[docs]def getpos(event):
'''
Get interactive mouse position
'''
global pfig,px,py,cid
px=event.xdata
py=event.ydata
pfig.canvas.mpl_disconnect(cid)
pfig.canvas.manager.window.destroy()
# ================================================================
[docs]def findpeaks(x,y,debug=False,fitas=False,manfit=False,pltax=None):
'''
Find peaks in spectrum, and compute approximate energy scale
Args:
* `x` = the x-values of the spectrum (arbitrary scale)
* `y` = the y-values(counts) of the spectrum
Kwargs:
* `debug` = when **True** show extra info for debugging
* `fitas` = fit the Al and Si lines in the spectrum
* `manfit` = manually fit the Al and Si lines
* `pltax` = plot instance for the plot. If **None**, ouput to screen.
Returns:
* `escale` = approximate energy scale
* `pindx` = (sorted) index positions of lines used
* `epos` = (sorted) energies of lines used
* `lfit` = the lines used for the energy calibration
'''
global epos,pfig,px,py,cid,eparatio,epalids,alpha,beta
# auto-correlation of the spectrum to find the peak widths
n=50 # number of bins for autocorrelation
corr=numpy.zeros(2*n,dtype=float)
xcor=numpy.arange(n*2)
rspec=numpy.roll(y,-n)
for i in numpy.arange(n*2):
corr[i]=corr[i]+numpy.sum(y*rspec)
rspec=numpy.roll(rspec,1)
# fit gauss function to top 1/3 of curve
coff=0.66*corr.max() # cutoff at 2/3 of max
indx,=numpy.where( corr > coff ) # data to use
if indx.size < 10:
print("mxsspec.findpeaks: insufficient width for gaussian fit\n",
" increase number of spectral bins")
return (None,None,None,None,None)
p,perr,yfit=gaussfit(xcor[indx],corr[indx])
if debug:
gg=gauss(p,xcor)
f,ax=plt.subplots(1)
ax.plot(xcor-n,corr,xcor-n,gg)
ax.set_xlabel('X-scale (arb. units)')
ax.set_ylabel('Cross corr.')
ax.set_title('Gauss fit to crosscorr.')
plt.show()
plt.close('all')
lwidth=p[2]/numpy.sqrt(2.0)*(x[1]-x[0])
yy=y
yclean=numpy.zeros(y.size,dtype=float)
frac=0.75
# frac=1.0
level=numpy.mean(yy)
if debug:
print("spectrum overall level: ",level)
xlwidth=p[2]/numpy.sqrt(2.0)
for i in numpy.arange(50): # 50 iteretions of 'clean' to find major peaks
mx=numpy.argmax(yy)
pos=x[mx]
top=numpy.mean(yy[(mx-1):(mx+2)])
# top=yy[mx]
if top < (2*level):
if debug:
print("top, level: ",top,level)
yy[mx]=0.0
continue
gg=gauss(numpy.array([top,pos,lwidth]),x)
if debug and ( i < 8 ):
f,ax=plt.subplots(1)
ax.plot(x,yy,x,gg,x,yclean)
ax.set_xlabel('X-scale (arb. units)')
ax.set_ylabel('Spectrum (counts)')
ax.set_title('line fit')
plt.show()
plt.close('all')
yy=yy-frac*gg
zindx,=numpy.where( yy < 0 )
yy[zindx]=0.0
yclean[mx]=yclean[mx]+frac*top
if debug:
f,ax=plt.subplots(1)
ax.plot(x,y,'b-',x,yy,'g-',x,yclean,'r-')
ax.set_xlabel('X-scale (arb. units)')
ax.set_ylabel('Spectrum (counts)')
ax.set_title('line fit')
plt.show()
plt.close('all')
merge=True # merge peaks which are less than 8 sigma apart
mdist=8.0*lwidth/(x[1]-x[0])
if debug:
print("Max peak distance merge:",mdist)
if debug:
pindx,=numpy.where( yclean > 0.0 ) # peaks
dists=pindx[1:]-pindx[:-1]
print('dists: ',dists)
nbins=int(max((10,(dists.max()/(2*lwidth/(x[1]-x[0])))))+0.5)
print('peak dist nbins: ',nbins)
f,ax=plt.subplots(1)
ax.hist(dists,bins=nbins)
ax.set_xlabel('Peak distance')
ax.set_ylabel('N')
ax.set_title('Peak distances')
plt.show()
plt.close('all')
hh,bb=numpy.histogram(dists,bins=nbins)
zz,=numpy.where( hh == 0 )
# mdist=bb[zz[0]]
print("Max peak distance merge:",mdist)
while merge:
pindx,=numpy.where( yclean > 0.0 ) # peaks found
dist=numpy.concatenate((pindx[1:],[y.size-1]))-pindx
dindx=numpy.argsort(dist)
if dist[dindx[0]] < mdist:
p1=pindx[dindx[0]] # first peak to merge
p2=pindx[dindx[0]+1] # second peak to merge
i1=yclean[p1] # first peak intensity
i2=yclean[p2] # second peak intensity
itot=i1+i2 # merged peak intensity
pnew=(p1*i1+p2*i2)/(i1+i2) # merged peak position
# print("Merge peak dist: ",dist[dindx[0]],p1,p2,i1,i2,pnew)
yclean[p1]=0.0
yclean[p2]=0.0
yclean[int(pnew+0.5)]=itot
else:
merge=False
# find proper line identifications.
# start with two brightest (alpha) lines.
xx=numpy.arange(x.size)
# if debug:
# print("found line index,positions: ",pindx,x[pindx])
lx=numpy.argsort(yclean[pindx]) # sort peaks on brightness
if debug:
print("3 brightest lines (positions,intensity): ",x[pindx[lx[-3:]]],yclean[pindx[lx[-3:]]])
l1=pindx[lx[-1]] # brightest line position
l2=pindx[lx[-2]] # second brightest line position
l3=pindx[lx[-3]] # second brightest line position
if l2 > l1: # make l1 the line at highest energy
l0=l1
l1=l2
l2=l0
# lminb=yclean[pindx[lx[-7]]] # minimum brightness for a valid beta peak
ratio=float(x[l1])/float(x[l2])
pparatio=numpy.array([1.07416140,1.366346,1.2720121],dtype=float) # best estimate of positions for
# nonlinear energy scale; ratios of Mn/Cr, Cu/Cr, Cu/Mn
if debug:
print("l1,x[l1],l2,x[l2], ratio: ",l1,x[l1],l2,x[l2],ratio)
print("eparatio:",eparatio)
print("pparatio:",pparatio)
print("epos: ",epos)
drat=0.05
for k,paratio in enumerate(pparatio):
if ( ratio > ((1.00-drat)*paratio)) and ( ratio < ((1.00+drat)*paratio)): # do ratios match within 5%?
p1=epos[alpha[epalids[k][0]]] # line position of l1
p2=epos[alpha[epalids[k][1]]] # line position of l2
lother=epalids[k][2] # other alpha line
break
if debug:
print("lines,energies:",l1,p1,l2,p2)
pfit=numpy.polyfit(numpy.array([l1,l2],dtype=float),numpy.array([p1,p2],dtype=float),1)
escale=pfit[0]*xx+pfit[1] # energy scale based on linear fit of two lines
lfit={'mn':[False,False],'cu':[True,False],'cr':[True,False]}
if k == 1:
if yclean[l3] > (0.5*yclean[l2]): # is the 3rd line a more than half the Cr-alpha intensity?
Mnline=True # yes, this must be the Mn_alpha
else: # if not, this must be the Cr-beta line
Mnline=False # prelimenary scale based on two lines only (Cr,Cu) no Mn
pindx=numpy.array([l1,l2],dtype=float)
fepos=numpy.array([p1,p2],dtype=float)
eposfit=pfit[0]*pindx+pfit[1]
if debug:
print("No Mn line -------------------------------------------------")
print("1st-2nd lines position, intensity: ",x[l1],x[l2],yclean[l1],yclean[l2])
print("3rd line position, intensity, ratio: ",x[l3],yclean[l3],(x[l3]/x[l2]))
print("line index: %8d %8d" % (l1,l2))
print("line position: %8.6f %8.6f" % (x[l1],x[l2]))
print("line energy: %8.1f %8.1f" % (p1,p2))
else:
Mnline=True
if Mnline:
pother=epos[alpha[lother]]
if(debug):
print("pother: ",pother, "minscale:",escale.min()," maxscale",escale.max())
sw=0.05
l3indx,=numpy.where( ( escale > ((1.0-sw)*pother) ) & \
( escale < ((1.0+sw)*pother) ) ) # approximate search area for 3rd line
yclean[l1]=0.0 # clear brightest lines, since they have been found
yclean[l2]=0.0
l3=yclean[l3indx].argmax()+l3indx[0] # 3rd alpha line is brightest in search area
p3=epos[alpha[lother]]
yclean[l3]=0.0
if debug:
f,ax=plt.subplots(1)
xplt=numpy.arange(y.size)
ax.plot(xplt,y,'b-',xplt[l3indx],y[l3indx],'r-')
yline=numpy.array([0,y.max()])
ax.plot([l1,l1],yline,'g-',[l2,l2],yline,'g-')
ax.plot(xplt,yclean,'c-')
ax.plot([l3,l3],yline,'m-')
plt.show()
plt.close('all')
print("line index: %8d %8d %8d" % (l1,l2,l3))
print("line position: %8.6f %8.6f %8.6f" % (x[l1],x[l2],x[l3]))
print("line energy: %8.1f %8.1f %8.1f" % (p1,p2,p3))
pindx=numpy.array([l1,l2,l3],dtype=float)
fepos=numpy.array([p1,p2,p3],dtype=float)
pfit,eposfit=expfit(pindx,fepos,debug=debug) # first estimate of energy scale based on alpha lines
if pfit is None: # only update energy scale when fit is ok, otherwise return
return (escale,pindx,fepos,xx,lfit)
lfit['mn'][0]=True
escale=expfunc(pfit,xx)
if debug:
f,ax=plt.subplots(1)
ax.plot(pindx,eposfit,'bo',xx,escale,'r-')
ax.set_title('First estimate energy fit')
ax.set_xlabel('spectrum bin')
ax.set_ylabel('Energy (eV)')
plt.show()
plt.close('all')
bs=0.05 # relative search width for beta lines
for j,bpos in enumerate(epos[beta]): # find the associated beta lines
if ( not Mnline ) and ( etype[beta[j]] == 'mn' ):
continue
bindx,=numpy.where( ( escale > ( (1.0-bs)*bpos ) ) & ( escale < ( (1.0+bs)*bpos ) ) )
lb=yclean[bindx].argmax()+bindx[0] # position of highest peak in interval
if yclean[lb] > (3.0*level): # only valid if brightness > minimum level
pindx=numpy.concatenate((pindx,[lb]))
fepos=numpy.concatenate((fepos,[bpos]))
lfit[etype[beta[j]]][1]=True
yclean[lb]=0.0
if debug:
print('pindx: ',pindx)
print('fepos: ',fepos)
f,ax=plt.subplots(2)
ax[0].plot(escale,y,'b-')
ax[0].plot(escale,yclean,'g-')
ax[0].set_xlabel('Energy (keV)')
ax[1].plot(xx,y,'b-')
yline=numpy.array([0,y.max()])
for ppos in pindx:
ax[1].plot([ppos,ppos],yline,'r-')
plt.show()
plt.close('all')
pfit,eposfit=expfit(pindx,fepos) # second estimate of energy scale based on alpha and beta lines
escale=expfunc(pfit,xx)
if debug:
print("eposfit: ",eposfit)
print(" fepos: ",fepos)
if numpy.absolute(eposfit-fepos).max() > 100.0: # more than 100 eV mismatch in energy for a line?
pindx=pindx[:3] # discard beta lines for energy fit
fepos=fepos[:3]
pfit,eposfit=expfit(pindx,fepos) # estimate of energy scale without beta lines
escale=expfunc(pfit,xx)
print("Prelimenary energy scale based on alpha lines only")
else:
print("Prelimenary energy scale based on both alpha and beta lines")
# fit Al and Si alpha lines ??
if fitas:
asindx,=numpy.where( ( escale > 1000.0 ) & ( escale < 2200.0 ) ) # region (eV) of Al and Si lines
elpos=numpy.zeros(2,dtype=float)
eltxt=['Al','Si']
cnv=numpy.array([1,3,6,3,1],dtype=float)
cnv=cnv/cnv.sum()
yy=numpy.convolve(y[asindx].astype(float),cnv,'same')
sw=2 # asumed width of peak for center of mass computation
for i in [0,1]:
if manfit: # fit manually ?
pfig,ax=plt.subplots(1)
ax.plot(xx[asindx],yy)
ax.set_xlabel('Optimal fit scale')
ax.set_ylabel('Counts')
ax.set_title('Select '+eltxt[i]+'-alpha line')
cursor=Cursor(ax, useblit=True, color='red', linewidth=1 ) # show cursor with crosshair
cid = pfig.canvas.mpl_connect('button_press_event', getpos) # activate 'onclick' routine
plt.show()
elpos[i]=px
print('Manual line position for '+eltxt[i]+':',px)
else:
mp=yy.argmax() # position of highest peak
elpos[i]=numpy.sum(xx[asindx[(mp-sw):(mp+sw+1)]]*yy[(mp-sw):(mp+sw+1)])/\
numpy.sum(yy[(mp-sw):(mp+sw+1)])
yy[(mp-sw-20):(mp+sw+21)]=0.0
elpos.sort()
if debug:
fig,ax=plt.subplots(1)
ax.plot(xx[asindx],y[asindx],'b-')
ax.plot([elpos[0],elpos[0]],[0.0,y[asindx].max()],'r-')
ax.plot([elpos[1],elpos[1]],[0.0,y[asindx].max()],'r-')
ax.set_xlabel('Optimal fit scale')
ax.set_ylabel('Counts')
ax.set_title('Al and Si line positions')
plt.show()
plt.close('all')
pindx=numpy.concatenate((elpos,pindx))
fepos=numpy.concatenate((numpy.array([1486.,1739.]),fepos))
if debug:
print("epos: ",fepos)
print("pindx: ",pindx)
# pfit=numpy.polyfit(pindx,epos,2)
# xx=numpy.arange(x.size)
sxpindx=pindx.argsort() # sort lines found
pindx=pindx[sxpindx]
fepos=fepos[sxpindx]
pfit,eposfit=expfit(pindx,fepos)
if pfit is None:
pfit=numpy.polyfit(pindx,fepos,1)
escale=pfit[0]*xx+pfit[1]
else:
escale=expfunc(pfit,xx)
if debug:
f,ax=plt.subplots(1)
ax.plot(pindx,fepos,'bo',xx,escale,'r-')
ax.set_xlabel('spectrum position')
ax.set_ylabel('line energies (eV)')
ax.set_title('approximate energy calibration')
plt.show()
plt.close('all')
f,ax=plt.subplots(1)
ax.plot(escale,y,'b-')
ax.set_xlabel('Approximate energy (eV)')
ax.set_ylabel('Counts')
ax.set_title('Spectrum')
plt.show()
plt.close('all')
axsize=1.2
if pltax is None:
fig = plt.figure()
ax1 = fig.add_subplot(111)
else:
if 'pecal' in pltax:
ax1 = pltax['pecal']
axsize=0.6
else:
return (escale,pindx,fepos,xx,lfit)
ax1.plot(pindx,fepos,'bo',xx,escale,'r-')
ax1.set_ylabel('line energies (eV)')
ax1.set_title('approximate energy calibration')
div = make_axes_locatable(ax1)
ax2 = div.append_axes("bottom", axsize, pad=0.0, sharex=ax1)
ax2.plot([xx[0],xx[-1]],[0.0,0.0],'r-')
ax2.plot(pindx,(fepos-escale[pindx.astype(int)]),'bo')
ax2.set_ylabel(r'$\Delta$(eV)')
ax2.set_xlabel('spectrum bin')
if pltax is None:
plt.tight_layout()
plt.show()
plt.close('all')
return (escale,pindx,fepos,xx,lfit)
# ============================================================================
[docs]def fitlines(energy,spectrum,debug=False,\
lfit={'mn':[True,True],'cu':[True,True],'cr':[True,True]},pltaxs=None,
lindx=None,lenergy=None,fixbgain=False,optbins=None,npoly=None):
'''
Fit the separate lines in the spectrum, using a limited spectral width around each line.
Args:
* `energy` = approximate energy scale
* `spectrum` = the spectrum (counts)
Kwargs:
* `debug` = print debug information when **True**
* `lfit` = dictionary of lines to fit. **[True,True]** will fit both alpha and beta lines
* `pltaxs` = dictionary of plot instances for line fit plots
* `lindx` = (sorted) spectrum index of lines used for approximate energy scale
* `lenergy` = (sorted) energies of lines used for approximate energy scale
* `fixbgain` = when **True** set gain for beta line holzer fit to alpha gain fitted value
* `optbins` = optimal fit amplitude scale (corresponds to **energy** )
* `npoly` = degree of polynomial for energy scale. When **None** use exponential fit
Returns:
* `escale` = fitted energy scale
* `afitpar` = fitted parameters for alpha lines
'''
global abpos,epos
fwidth=30.0 # width (eV) around line position taken for fitting the line
ltypes=['Alpha','Beta']
eekev=energy/1000.0 # energy in keV
nmax=10
resolutions=numpy.zeros((2,nmax),dtype=float)
reserr=numpy.zeros((2,nmax),dtype=float)
shifts=numpy.zeros((2,nmax),dtype=float)
shifterr=numpy.zeros((2,nmax),dtype=float)
shiftpos=numpy.zeros((2,nmax),dtype=float)
gains=numpy.zeros((2,nmax),dtype=float)
gainerr=numpy.zeros((2,nmax),dtype=float)
energies=numpy.zeros(nmax)
afitpar={} # alpha lines fit results
j=0
for element in lfit: # go through all elements
for ltype in [0,1]: # go through both type of lines (Alpha,Beta)
if lfit[element][ltype]: # do we want to process this line?
lpos=abpos[element][ltype] # line position in eV
e1=lpos-fwidth
e2=lpos+fwidth
indx,=numpy.where( ( energy >= e1 ) & ( energy <= e2 ) )
lwindx,=numpy.where( ( energy >= (e1-fwidth) ) & ( energy <= e1 ) )
hwindx,=numpy.where( ( energy >= e2 ) & ( energy <= (e2+fwidth) ) )
br1=spectrum[lwindx].sum()/float(lwindx.size) # (brems) offset level low wing
br2=spectrum[hwindx].sum()/float(hwindx.size) # (brems) offset level high wing
boffset=min([br1,br2]) # estmated bremsstrahlung level from minimum of winglevels
if debug:
print('Bremsstrahlung level for ',element,ltypes[ltype],' : ',boffset)
if ( ltype == 0 ) or ( not fixbgain ):
bgain=None # no fixed beta gain for alpha line
fitpar=holzerfit(eekev[indx],spectrum[indx],line=ltypes[ltype],elements=[element],\
offset=boffset,bgain=bgain)
windx,=numpy.where( ( energy >= (e1-fwidth) ) & ( energy <= (e2+fwidth) ) )
if ltype == 0:
bgain=fitpar[1] # use alpha gain for next fit of beta line
afitpar[element]=fitpar # return these fit parmeters
fit=holzer(fitpar,eekev[windx],offset=boffset)
else:
fit=hbeta(fitpar,eekev[windx],offset=boffset)
if optbins is not None: # compute lratio ?
ealpha=abpos[element][0]
ebeta=abpos[element][1]
abin=optbins[numpy.absolute(energy-(lpos+afitpar[element][0])).argmin()]
bbin=optbins[numpy.absolute(energy-(lpos+fitpar[0])).argmin()]
lratio=(abin/bbin)/(ealpha/ebeta)
afitpar[element]=numpy.concatenate((afitpar[element],[lratio]))
specerr=numpy.sqrt(spectrum[windx])
tt=['Shift:','Gain:','Counts:','FWHM:','C-stat:','Degr. of frdm:']
ptxt= '%7s %7.4f\n' % (tt[0],fitpar[0])
ptxt=ptxt+'%7s %7.4f\n' % (tt[1],fitpar[1])
ptxt=ptxt+'%7s %7d\n' % (tt[2],fitpar[2])
ptxt=ptxt+'%7s %5.2f (eV)\n' % (tt[3],fitpar[3]*1000.0)
ptxt=ptxt+'%7s %7.1f\n' % (tt[4],fitpar[4])
ptxt=ptxt+'%s %d\n' % (tt[5],fitpar[5])
if debug:
fiterrplot(energy[windx],spectrum[windx],specerr,fit,\
ptxt=ptxt,ptitle=(element+" "+ltypes[ltype]),\
xtitle='Energy (eV)',ytitle='Counts')
if ( ltype == 0 ) and ( pltaxs is not None ):
if element in pltaxs:
fiterrplot(energy[windx],spectrum[windx],specerr,fit,\
ptxt=ptxt,ptitle=(element+" "+ltypes[ltype]),\
xtitle='Energy (eV)',ytitle='Counts',pltax=pltaxs[element],
ptxtsize=7)
resolutions[ltype,j//2]=fitpar[3]*1000.0
reserr[ltype,j//2]=fitpar[9]*1000.0
shifts[ltype,j//2]=fitpar[0]
shifterr[ltype,j//2]=fitpar[6]
shiftpos[ltype,j//2]=lpos
gains[ltype,j//2]=fitpar[1]
gainerr[ltype,j//2]=fitpar[7]
energies[j]=lpos
j=j+1
if debug:
arange=numpy.arange(0,j,2)
brange=arange+1
fig,ax=plt.subplots(3,sharex=True)
ax[0].errorbar(energies[brange],resolutions[1,:j//2],yerr=reserr[1,:j//2],fmt='bo')
ax[0].errorbar(energies[arange],resolutions[0,:j//2],yerr=reserr[0,:j//2],fmt='ro')
ax[0].set_ylabel('FWHM resolultion')
ax[1].errorbar(energies[brange],shifts[1,:j//2],shifterr[1,:j//2],fmt='bo')
ax[1].errorbar(energies[arange],shifts[0,:j//2],shifterr[0,:j//2],fmt='ro')
ax[1].set_ylabel('Shift (eV)')
ax[2].errorbar(energies[brange],gains[1,:j//2],gainerr[1,:j//2],fmt='bo')
ax[2].errorbar(energies[arange],gains[0,:j//2],gainerr[0,:j//2],fmt='ro')
ax[2].set_xlabel('Energy (eV)')
ax[2].set_ylabel('Gain')
plt.show()
plt.close('all')
if pltaxs is not None:
if 'fwhm' in pltaxs:
pltaxs['fwhm'].errorbar(energies[brange],resolutions[1,:j//2],yerr=reserr[1,:j//2],fmt='bo')
pltaxs['fwhm'].errorbar(energies[arange],resolutions[0,:j//2],yerr=reserr[0,:j//2],fmt='ro')
pltaxs['fwhm'].set_xlabel('Energy (eV)')
pltaxs['fwhm'].set_ylabel('FWHM resolultion')
if 'shift' in pltaxs:
pltaxs['shift'].errorbar(energies[brange],shifts[1,:j//2],shifterr[1,:j//2],fmt='bo')
pltaxs['shift'].errorbar(energies[arange],shifts[0,:j//2],shifterr[0,:j//2],fmt='ro')
pltaxs['shift'].set_xlabel('Energy (eV)')
pltaxs['shift'].set_ylabel('Shift (eV)')
if 'gain' in pltaxs:
pltaxs['gain'].errorbar(energies[brange],gains[1,:j//2],gainerr[1,:j//2],fmt='bo')
pltaxs['gain'].errorbar(energies[arange],gains[0,:j//2],gainerr[0,:j//2],fmt='ro')
pltaxs['gain'].set_xlabel('Energy (eV)')
pltaxs['gain'].set_ylabel('Gain')
# correct energy scale for fitted shifts
ashifts=shifts[0,:]
ashiftpos=shiftpos[0,:]
ashifterr=shifterr[0,:]
jb=j
j=j//2
# include beta lines
ashifts[j:jb]=shifts[1,:j]
ashiftpos[j:jb]=shiftpos[1,:j]
ashifterr[j:jb]=shifterr[1,:j]+2.0 # assume extra 2eV systematic error for Beta lines
jp=j
j=jb
# print("lenergy: ",lenergy)
# print("epos: ",epos)
if ( lenergy[1] < epos.min() ) and ( ( j == 3 ) or ( j == 6 ) ): # add Al and Si line positions
asindx,=numpy.where( ( energy > 1000.0 ) & ( energy < 2200.0 ) ) # region (eV) of Al and Si lines
elpos=numpy.zeros(2,dtype=float)
cnv=numpy.array([1,3,6,3,1],dtype=float)
cnv=cnv/cnv.sum()
yy=numpy.convolve(spectrum[asindx].astype(float),cnv,'same')
sw=2 # asumed width of peak for center of mass computation
for i in [0,1]:
mp=yy.argmax() # position of highest peak
elpos[i]=numpy.sum(energy[asindx[(mp-sw):(mp+sw+1)]]*yy[(mp-sw):(mp+sw+1)])/\
numpy.sum(yy[(mp-sw):(mp+sw+1)])
yy[(mp-sw-20):(mp+sw+21)]=0.0
elpos.sort()
# print("elpos: ",elpos)
# print("energy[lindx[:2]]:",energy[lindx[:2].astype(int)])
mxerr=ashifterr[:jp].max()*1.5
ashifts[j:(j+2)]=-(energy[lindx[:2].astype(int)]-elpos)
ashiftpos[j:(j+2)]=elpos
ashifterr[j:(j+2)]=[mxerr,mxerr] # manually set Al and Si line uncertainties
j=j+2
if debug:
print("shifts: ",ashifts)
print("ashiftpos: ",ashiftpos)
print("ashifterr: ",ashifterr)
f,ax=plt.subplots(1)
ax.plot(energy[asindx],yy,'g-')
ax.plot(energy[asindx],spectrum[asindx],'b-')
yline=ax.axis()[2:]
ax.plot([elpos[0],elpos[0]],yline,'r-')
ax.plot([elpos[1],elpos[1]],yline,'r-')
ax.set_title('Refit Al/Si lines')
plt.show()
plt.close('all')
ashiftw=1./ashifterr[:j] # weights for the shifts
posindx=numpy.zeros(ashifts.size,dtype=int)
for i in numpy.arange(j):
dd=numpy.absolute(energy-ashiftpos[i])
posindx[i]=dd.argmin() # index for given line energy
if debug:
print("Energies: ",energy[posindx[:j]])
print(" shifts: ",ashifts[:j])
ee=energy[posindx[:j]]+ashifts[:j] # correct energy scale for shifts
xx=numpy.arange(energy.size)
if npoly is not None: # use polynomial energy scale for final fit, when npoly is not None
pnew=numpy.polyfit(posindx[:j],ee,npoly,w=ashiftw)
newscale=numpy.polyval(pnew,xx)
else: # otherwise use exponential fit
pnew,eefit=expfit(posindx[:j],ee,yerr=ashifterr[:j])
if pnew is None:
newscale=energy # keep old energy scale when fit cannot be made
else:
newscale=expfunc(pnew,xx)
if debug:
fig,ax=plt.subplots(2)
ax[0].plot(xx,newscale-energy,'r-')
ax[0].plot(posindx[:j],ashifts[:j],'bo')
ax[0].set_xlabel('indx')
ax[0].set_ylabel('energy-newscale')
ax[0].set_title('energy scale update')
ax[1].plot(newscale,spectrum,'b-')
ax[1].set_xlabel('New energy scale (eV)')
ax[1].set_ylabel('Counts')
plt.show()
plt.close('all')
else:
if pltaxs is not None:
if 'spectrum' in pltaxs:
pltaxs['spectrum'].plot(newscale,spectrum,'b-')
pltaxs['spectrum'].set_xlabel('Energy (eV)')
pltaxs['spectrum'].set_ylabel('Counts')
pltaxs['spectrum'].set_title('Calibrated spectrum')
if 'ecal' in pltaxs:
pltaxs['ecal'].plot(xx,newscale-energy,'r-')
pltaxs['ecal'].plot(posindx[:j],ashifts[:j],'bo')
pltaxs['ecal'].set_xlabel('spectrum bin')
pltaxs['ecal'].set_ylabel('Delta update (eV)')
pltaxs['ecal'].set_title('Line fit differences')
return (newscale,afitpar)
# ===================================================================
[docs]def mbasecorr(optfit,bins,escale,basls,debug=False,pltax=None,pdf=None,ptitle=""):
'''
Compute baseline corrections as function of energy
Args:
* `optfit` = optimal fit parameters of the events
* `bins` = optimal fit bins of the spectrum
* `escale` = energies for the bins
* `basls` = eventlist containing the baselines of the events
Kwargs:
* `debug` = when **True**, provide debug information
* `pltax` = plotinstance for result plot
* `pdf` = pdf file for pdf plot
* `ptitle` = title for pdf plot
Returns:
* `tfit` = linear fit of baseline-dependence as function of energy
* `avbaseline` = average baseline value
'''
global abpos
fwidth=30.0 # width around alpha lines for correlation
elements=['cr','mn','cu']
if debug or ( pdf is not None ):
f,ax=plt.subplots(len(elements))
plt.suptitle(ptitle,size=12)
centers=numpy.zeros(len(elements)+1)
tilts=numpy.zeros(len(elements)+1)
avbaseline=0.0
for k,element in enumerate(elements): # go through all elements
lpos=abpos[element][0] # position of alpha line
e1=lpos-fwidth
e2=lpos+fwidth
eindx,=numpy.where( ( escale > e1 ) & ( escale < e2 ) )
lindx,=numpy.where( ( optfit > bins[eindx[0]] ) & ( optfit < bins[eindx[-1]] ) )
newoptfit,center,ytilt=fitcor(basls[lindx],optfit[lindx],\
nolim=True,debug=debug) # find linear correlation
centers[k]=center[1]
tilts[k]=ytilt
avbaseline=avbaseline+center[0]
if debug or ( pdf is not None ):
# ax[k].plot(basls[lindx],optfit[lindx],'b.')
xcorplot(ax[k],basls[lindx],optfit[lindx])
rr=ax[k].axis()
y1=center[1]+(rr[0]-center[0])*ytilt
y2=center[1]+(rr[1]-center[0])*ytilt
ax[k].plot([rr[0],rr[1]],[y1,y2],'r-',label=('%5.3f' % (ytilt*100000.0)))
ax[k].legend()
ax[k].set_ylabel('Opt.fit')
avbaseline=avbaseline/float(len(elements))
# tfit=numpy.polyfit(centers,tilts,1)
tfit=numpy.zeros(2,dtype=float)
tfit[0]=tilts[:-1].sum()/centers[:-1].sum() # force fit through (0,0)
if debug or ( pdf is not None):
ax[(len(elements)-1)].set_xlabel('Baseline')
if debug:
plt.show()
plt.close('all')
f,ax=plt.subplots(1)
ax.plot(centers,tilts,'bo')
rr=numpy.array(ax.axis())
lfit=tfit[0]*rr[:2]+tfit[1]
ax.plot(rr[:2],lfit,'r-')
plt.show()
print('tilt fit: ',tfit,' average baseline: ',avbaseline)
else:
makeplot(pdf,ps=True)
plt.close('all')
if pltax is not None:
if 'tilts' in pltax:
pltax['tilts'].plot(centers,tilts,'bo')
rr=numpy.array(pltax['tilts'].axis())
lfit=tfit[0]*rr[:2]+tfit[1]
pltax['tilts'].plot(rr[:2],lfit,'r-')
pltax['tilts'].set_xlabel('Optimal fit par.')
pltax['tilts'].set_ylabel('Correction')
pltax['tilts'].set_title('Baseline correction')
coptfit=optfit-(tfit[0]*optfit+tfit[1])*(basls-avbaseline) # correct energies for baseline
cspectrum,bins=numpy.histogram(coptfit,range=[bins[0],bins[-1]],\
bins=bins.size) # make corrected spectrum
return (cspectrum,(centers,tilts,tfit))
# ==========================================================================
[docs]def pixpars(hdf,channel,freq,spectrum):
'''
Make dictionary with pixel info
Args:
* `hdf` = HDF instrance of the hdf data
* `channel` = selected FDM channel
* `freq` = select frequency (pixel)
* `spectrum` = the events spectrum
Returns:
* `fitpars` = dictionary with pixel parameters
'''
fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq )
freqindx=fffx[0]
khz=hdf.channel[channel].freq_conf['freq'][freqindx]
fname=os.path.basename(hdf.hdffile).rsplit('.',1)[0]
fid="%s ch=%d px=%d %7.2f (kHz)" % (fname,channel,freq,khz) # plot title
fitpars={'khz':khz,'counts':numpy.sum(spectrum),'fid':fid}
return fitpars
[docs]def pixinfo(hdf,channel,freq):
'''
Return pixel info
Args:
* `hdf` = HDF instrance of the hdf data
* `channel` = selected FDM channel
* `freq` = select frequency (pixel)
Returns:
* `fid` = pixel info
'''
fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq )
freqindx=fffx[0]
khz=hdf.channel[channel].freq_conf['freq'][freqindx]
fid="ch=%d, px=%d %7.2f (kHz)" % (channel,freq,khz) # plot title
return fid
# =========================================================================
[docs]def showescale(pltax,escale,pindx,epos,xx):
'''
Show plot of energy scale fit
Args:
* `pltax` = plot instance for plot
* `escale` = energy scale
* `pindx` = fitted line index
* `epos` = line positions
* `xx` = bins for energy scale
'''
if pltax is None:
return
else:
if 'pecal' in pltax:
ax1 = pltax['pecal']
axsize=0.6
else:
return
ax1.plot(pindx,epos,'bo',xx,escale,'r-')
ax1.set_ylabel('line energies (eV)')
ax1.set_title('approximate energy calibration')
div = make_axes_locatable(ax1)
ax2 = div.append_axes("bottom", axsize, pad=0.0, sharex=ax1)
ax2.plot([xx[0],xx[-1]],[0.0,0.0],'r-')
ax2.plot(pindx,(epos-escale[pindx.astype(int)]),'bo')
ax2.set_ylabel(r'$\Delta$(eV)')
ax2.set_xlabel('spectrum bin')
# ===========================================================================
[docs]def showtilts(pltax,tpar):
'''
Show plot of baseline correlaion results
Args:
* `tpar` = tilt correlation parameters
* `pltax` = plotinstance for plot
'''
centers=tpar[0]
tilts=tpar[1]
tfit=tpar[2]
if pltax is not None:
if 'tilts' in pltax:
pltax['tilts'].plot(centers,tilts,'bo')
rr=numpy.array(pltax['tilts'].axis())
lfit=tfit[0]*rr[:2]+tfit[1]
pltax['tilts'].plot(rr[:2],lfit,'r-')
pltax['tilts'].set_xlabel('Optimal fit par.')
pltax['tilts'].set_ylabel('Correction')
pltax['tilts'].set_title('Baseline correction')