#!/usr/bin/env python3
#
'''
.. module:: spectrumfit
:synopsis: fit the Holzer alpha and beta Mn54 lines to the spectrum
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import matplotlib
if __name__ == "__main__":
matplotlib.use('GTK3Agg')
import matplotlib.pyplot as plt
import numpy
from tesfdmtools.methods.holzer import holzer,holzerfit
[docs]def spectrumfit(spectrum,bins,debug=False,afitrange=[5.860,5.920]):
'''
fit the Holzer alpha and beta Mn54 lines to the spectrum
Args:
* `spectrum` = spectrum to be fitted
* `bins` = centers of the spectrum bins
* `debug` = if True, plot intermediate fit solution
* `afitrange` = range in keV of the fit window for the Mn-55 K-alpha line [default: 5.860,5.920]
Returns:
* `nppp` = best fit parameters for Mn-55 K-alpha complex
* nppp[0] : shift (eV)
* nppp[1] : compression factor (gain)
* nppp[2] : normalization
* nppp[3] : gaussian convolution width (keV)
* nppp[4] : C-stat
* nppp[5] : degrees of freedom
* nppp[6,7,8,9] : 1-sigma errors of nppp[0,1,2,3]
* nppp[10,11] : conversion parameters: energy = nppp[10]+nppp[11]*bins
* `spectrum` = input spectrum within fit window (K-Alpha lines)
* `xxx` = fitted energy axis for spectrum
* `holy` = holzer curve fit for xxx energy axis
* `lratio` = linearity ratio: (alpha-bin/beta-bin)/(E-alpha/E-beta)
'''
print("fitspectrum")
palpha=spectrum.argmax() # pos of alpha line
if (spectrum.size-3) <= palpha:
return (None,None,None,None,None)
for i in numpy.arange(palpha,(spectrum.size-2)): # search from position of alpha max for end of alpha line energies
if ( spectrum[i] == 0 ) and ( spectrum[i+1] == 0 ) and ( spectrum[i+2] == 0 ) and ( spectrum[i+3] == 0 ):
break
pz=i
pbeta=spectrum[pz:].argmax()+pz # pos of beta line
print("positions of alpha, beta lines (bins): ",bins[palpha],bins[pbeta])
if ( bins[pbeta] != bins[pbeta] ) or ( bins[palpha] != bins[palpha] ): # check if valid data (check for NaN)
return (None,None,None,None,None)
binw=(6.49045-5.897867)/(pbeta-palpha)
aee=(numpy.arange(spectrum.size)-palpha)*binw+5.897867 # approximate energy scale
bindx,=numpy.where( (aee > 6.4200) & (aee < 6.5200)) # beta line fit positions
if bindx.size <= 1:
return (None,None,None,None,None)
ppp=holzerfit(aee[bindx],spectrum[bindx],line='Beta')
energy=aee+ppp[0]/1000.0
print("beta line shift (prel. energy scale): ",ppp[0])
pmin=afitrange[0] # fixed alpha line interval
pmax=afitrange[1]
print("alpha lines interval: ",pmin,pmax)
indx,=numpy.where((energy >= pmin) & (energy <= pmax) )
if indx.size <= 1:
return (None,None,None,None,None)
ppp=numpy.zeros(10,dtype=float)
gindx=False
if indx.size > 5:
ppp=holzerfit(energy[indx],spectrum[indx])
gindx=True
else:
indx=numpy.arange(energy.size)
print("alpha line shift (prel. energy scale): ",ppp[0])
avapos=(5897.867+5887.743)/2.0/1000.0 # average alpha line position (keV)
bpos=6.49045 # position of beta line
avashift=ppp[0] # average alpha line shift (eV)
avashiev=avashift/1000.0 # shift in keV
xxx=(energy[indx]-bpos)/((avapos-avashiev)-bpos)*avashiev+energy[indx] # proper energy axis
if debug:
fxxx=(energy-bpos)/((avapos-avashiev)-bpos)*avashiev+energy
f, axarr = plt.subplots(1)
axarr.set_xlabel('Energy (keV)')
axarr.set_ylabel('Counts')
axarr.set_title('Calibrated spectrum')
axarr.plot(fxxx,spectrum)
rr=axarr.axis()
yr=[rr[2],rr[3]]
axarr.plot([energy[palpha],energy[palpha]],yr,'r-')
axarr.plot([energy[pbeta],energy[pbeta]],yr,'r-')
axarr.plot([pmin,pmin],yr,'g-')
axarr.plot([pmax,pmax],yr,'g-')
plt.show()
plt.close('all')
nppp=numpy.zeros(12,dtype=float)
if gindx:
nppp[0:10]=holzerfit(xxx,spectrum[indx])
holy=holzer(nppp,xxx)
else:
nppp[0:10]=ppp
holy=numpy.zeros(xxx.size)
e2b=numpy.polyfit(xxx,bins[indx],1) # conversion from energy to bins
abin=e2b[1]+e2b[0]*avapos # bin for alpha position
bbin=e2b[1]+e2b[0]*bpos # bin for beta position
lratio=(abin/bbin)/(avapos/bpos) # liniarity ratio
nppp[10]=-e2b[1]/e2b[0] # pass on conversion bin -> eenergy
nppp[11]=1.0/e2b[0]
print("alpha,beta bin positions: %6.3f %6.3f :: linearity ratio: %7.4f" % (abin,bbin,lratio) )
return (nppp,spectrum[indx],xxx,holy,lratio)
#===========================================================================================
if __name__ == "__main__":
'''
Test program on fixed dataset
'''
from tesfdmtools.hdf.HMUX import HMUX
from tesfdmtools.utils.reventlist import reventlist
from tesfdmtools.methods.HDF_Eventutils import makespectrum,fitspectrum
from tesfdmtools.utils.fiterrplot import fiterrplot
indx,optfit,basel,peaks,hdffile,rtim,ftim,fenrg=reventlist("/stage/xmmdat10/Athena/new_data/20170624T51mK_test/"+\
"eventlist_Run62_201706241244_xray_100_100_100_px0_px000.lst")
spectrum,bins=makespectrum(optfit)
pltdata={}
iplt=1
hdf=HMUX(filename=hdffile)
fitspectrum(pltdata,iplt,hdf,0,0,spectrum,bins)
fp=pltdata[iplt]
f, ax = plt.subplots(1)
fiterrplot(fp['xax'],fp['spectrum'],fp['yerr'],fp['holy'],\
xtitle="Energy (keV)",ytitle="Counts",ptitle=fp['ptitle'],ptxt=fp['ptxt'],pltax=ax)
plt.show()
plt.close('all')
# nppp,spectrum,xxx,holy=fitspectrum(spectrum,bins,debug=True)