import numpy as np
from pydiana.tools import pyroot_fitting as prf
from pydiana.tools import dataframe_manip as dfman
from pydiana.tools import graph_tools as gt
from pydiana import ROOT,Diana
import pandas as pd
[docs]
def get_iron_data(num=2,decimals:int=0):
"""
Function for getting the iron peaks of Ka and Kb
Parameters:
-----------
num : integer specifying how many lines we want to return. If 3 it returns the two Ka lines and the Kb line, if 2 it averages the two Ka lines and returns one Ka and one Kb, if 1 it only returns the averaged Ka line.
decimals : integer saying how many decimal places to use. The approximation is done to the order of magnitude specified by 10^(-decimals) (if decimals is negative it will approximate to tens, hundreds, etc...)
Returns
-------
2D array: First row containing the energy of the iron lines in eV. The second row contains the intensities of the iron lines.
"""
intensities = [8.5,16.9,1.01+1.98]
energies = [5888,5899,6490]
if num==2:
energies = [(energies[0]*intensities[0]+energies[1]*intensities[1])/(intensities[0]+intensities[1]),energies[2]]
intensities = [intensities[0]+intensities[1],intensities[2]]
elif num==1:
energies = [(energies[0]*intensities[0]+energies[1]*intensities[1])/(intensities[0]+intensities[1])]
intensities = [intensities[0]+intensities[1]]
for i in range(len(energies)):
energies[i]=round(energies[i],decimals)
intensities[i]=round(intensities[i],decimals)
if decimals<=0:
energies[i]=int(energies[i])
intensities[i]=int(intensities[i])
return energies,intensities
[docs]
def iron_fit(histogram,params,
parfix:list = None,
unit='mV',
xrange:list=None,
advanced:bool=False,
constrained:bool=False,
degree:int=None,
polType:str='Chebyshev',
fit_options:str='QSL',
):
"""
Utility for fitting a double gaussian to the iron peak
Parameters
----------
histogram : dictionary containing "counts" with the heights of the histogram bars, and "bin_edges" containing the edges of the histogram bins (left and right)
params : array-like containing the parameters of the gaussians
parfix : array-like contatining the index of the parameters to fix
advanced : Boolean to switch off and on the advanced 3 gaussian fitting. It can be used to fit two separate gaussians in case the two different K_a lines can be resolved.
constrained : Boolean used to toggle a constained gaussian fitting (apart from the parameters of the first gaussian, the mean values and intensities of the remaining ones are physically determined).
Returns
-------
fit_results : pandas dataframe containing the fit results
func : root fitted function
"""
if xrange is None:
xrange = (min(histogram['bin_edges']),max(histogram['bin_edges']))
else:
xrange = (min(xrange),max(xrange))
expr=''
if advanced:
expr = 'gaus(0)+gaus(3)+gaus(6)+[9]'
parnames={
0:f"Constant_Ka1[counts]",
1:f"Mean_Ka1[{unit}]",
2:f"StdDev_Ka1[{unit}]",
3:f"Constant_Ka2[counts]",
4:f"Mean_Ka2[{unit}]",
5:f"StdDev_Ka2[{unit}]",
6:f"Constant_Kb[counts]",
7:f"Mean_Kb[{unit}]",
8:f"StdDev_Kb[{unit}]",
}
parlims={
0:(0,max(histogram['counts'])),
1:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
2:(0,max(histogram['bin_edges'])),
3:(0,max(histogram['counts'])),
4:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
5:(0,max(histogram['bin_edges'])),
6:(0,max(histogram['counts'])),
7:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
8:(0,max(histogram['bin_edges'])),
}
if constrained:
ironinfo =get_iron_data(3)
mean_ratio1 = ironinfo[0][1]/ironinfo[0][0]
int_ratio1 = ironinfo[1][1]/ironinfo[1][0]
mean_ratio2 = ironinfo[0][2]/ironinfo[0][0]
int_ratio2 = ironinfo[1][2]/ironinfo[1][0]
expr = f"gaus(0)+[0]*{int_ratio1}*exp(-0.5*((x-([1]*{mean_ratio1}))/[3])*((x-([1]*{mean_ratio1}))/[3]))+[0]*{int_ratio2}*exp(-0.5*((x-([1]*{mean_ratio2}))/[4])*((x-([1]*{mean_ratio2}))/[4]))"
parnames={
0:f"Constant_Ka1[counts]",
1:f"Mean_Ka1[{unit}]",
2:f"StdDev_Ka1[{unit}]",
3:f"StdDev_Ka2[{unit}]",
4:f"StdDev_Kb[{unit}]",
}
parlims={
0:(0,max(histogram['counts'])),
1:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
2:(0,max(histogram['bin_edges'])),
3:(0,max(histogram['bin_edges'])),
4:(0,max(histogram['bin_edges'])),
5:(0,max(histogram['counts'])),
}
else:
expr = 'gaus(0)+gaus(3)'
parnames={0:f"Constant_Ka[counts]",
1:f"Mean_Ka[{unit}]",
2:f"StdDev_Ka[{unit}]",
3:f"Constant_Kb[counts]",
4:f"Mean_Kb[{unit}]",
5:f"StdDev_Kb[{unit}]",
}
parlims={
0:(0,max(histogram['counts'])),
1:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
2:(0,max(histogram['bin_edges'])),
3:(0,max(histogram['counts'])),
4:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
5:(0,max(histogram['bin_edges'])),
}
if constrained:
ironinfo =get_iron_data(2)
mean_ratio = ironinfo[0][1]/ironinfo[0][0]
int_ratio = ironinfo[1][1]/ironinfo[1][0]
expr = f"gaus(0)+[0]*{int_ratio}*exp(-0.5*((x-([1]*{mean_ratio}))/[3])*((x-([1]*{mean_ratio}))/[3]))"
parnames={0:f"Constant_Ka[counts]",
1:f"Mean_Ka[{unit}]",
2:f"StdDev_Ka[{unit}]",
3:f"StdDev_Kb[{unit}]",
}
parlims={
0:(0,max(histogram['counts'])),
1:(min(histogram['bin_edges']),max(histogram['bin_edges'])),
2:(0,max(histogram['bin_edges'])),
3:(0,max(histogram['bin_edges'])),
}
expr += prf.make_background_function(degree=degree,start_idx=len(parnames),type=polType)
if degree is not None:
params = np.concatenate((params,[0]*degree))
func = prf.root_func(expr=expr,
name='Iron Peak Fit',
parnames=parnames,
parlims=parlims,
parvals={i:pp for i,pp in enumerate(params)},
parfix=parfix,
xrange=xrange)
fit_result,func = prf.root_hist_fit(heights=histogram['counts'],
bin_edges=histogram['bin_edges'],
func=func,
get_detail_matrices=False,
options=fit_options)
fit_result = prf.extract_fit_values(func)
fit_result['Units']['MinX']=unit
fit_result['Units']['MaxX']=unit
if not constrained:
tmp=pd.DataFrame()
if advanced:
tmp['Variables']=('Mean Ratio Ka2/Ka1',"Mean Ratio Kb/Ka1",'Intensity Ratio Ka2/Ka1',"Intensity Ratio Kb/Ka1")
mkaratio = fit_result['Values']['Mean_Ka2']/fit_result['Values']['Mean_Ka1']
mkaratioerr = np.sqrt((fit_result['Errors']['Mean_Ka2']/fit_result['Values']['Mean_Ka1'])**2 + (fit_result['Values']['Mean_Ka2']*fit_result['Errors']['Mean_Ka1']/fit_result['Values']['Mean_Ka1']**2)**2)
mkabratio = fit_result['Values']['Mean_Kb']/fit_result['Values']['Mean_Ka1']
mkabratioerr = np.sqrt((fit_result['Errors']['Mean_Kb']/fit_result['Values']['Mean_Ka1'])**2 + (fit_result['Values']['Mean_Kb']*fit_result['Errors']['Mean_Ka1']/fit_result['Values']['Mean_Ka1']**2)**2)
ckaratio = fit_result['Values']['Constant_Ka2']/fit_result['Values']['Constant_Ka1']
ckaratioerr = np.sqrt((fit_result['Errors']['Constant_Ka2']/fit_result['Values']['Constant_Ka1'])**2 + (fit_result['Values']['Constant_Ka2']*fit_result['Errors']['Constant_Ka1']/fit_result['Values']['Constant_Ka1']**2)**2)
ckabratio = fit_result['Values']['Constant_Ka1']/fit_result['Values']['Constant_Kb']
ckabratioerr = np.sqrt((fit_result['Errors']['Constant_Ka1']/fit_result['Values']['Constant_Kb'])**2 + (fit_result['Values']['Constant_Ka1']*fit_result['Errors']['Constant_Kb']/fit_result['Values']['Constant_Kb']**2)**2)
tmp['Values']=(mkaratio,mkabratio,ckaratio,ckabratio)
tmp['Errors']=(mkaratioerr,mkabratioerr,ckaratioerr,ckabratioerr)
else:
tmp['Variables']=("Mean Ratio Kb/Ka","Intensity Ratio Kb/Ka")
mkratio = fit_result['Values']['Mean_Kb']/fit_result['Values']['Mean_Ka']
mkratioerr = np.sqrt((fit_result['Errors']['Mean_Kb']/fit_result['Values']['Mean_Ka'])**2 + (fit_result['Values']['Mean_Kb']*fit_result['Errors']['Mean_Ka']/fit_result['Values']['Mean_Ka']**2)**2)
ckaratio = fit_result['Values']['Constant_Ka']/fit_result['Values']['Constant_Kb'] * fit_result['Values']['StdDev_Ka']/fit_result['Values']['StdDev_Kb']
ckaratioerr = np.sqrt((fit_result['Errors']['Constant_Ka']/fit_result['Values']['Constant_Ka'] * ckaratio)**2 +(fit_result['Errors']['Constant_Kb']*ckaratio/fit_result['Values']['Constant_Kb'])**2 + (fit_result['Errors']['StdDev_Ka']/fit_result['Values']['StdDev_Ka'] * ckaratio)**2 +(fit_result['Errors']['StdDev_Kb']*ckaratio/fit_result['Values']['StdDev_Kb'])**2 )
tmp['Values']=(mkratio,ckaratio)
tmp['Errors']=(mkratioerr,ckaratioerr)
tmp['Units']=["-"]*len(tmp)
tmp.set_index('Variables',inplace=True)
fit_result = dfman.dataframe_insert(fit_result,tmp,'Chi2')
return fit_result,func
[docs]
def iron_calibrate(fit_peaks,parfix=None,iron_peak_energies=None,iron_peak_errors=None,unit_iron='eV',get_detail_matrices:bool=False,showplot:bool=True):
"""
Function for fitting the iron calibration constant
Parameters
----------
fit_peaks : pandas dataframe with the fit result of the histogram fitting of the iron lines
iron_peak_energies : array-like with the energies of the 2 iron lines
iron_peak_errors : array-like with the error on the energies of the 2 iron lines
unit_iron : string with the unit of the iron energies
Returns
-------
fit_results : pandas dataframe containing the fit results
func : root fitted function
"""
mviron = [fit_peaks['Values'][v] for v in fit_peaks.index if 'Mean_K' in v]
mvironerr = [fit_peaks['Errors'][v] for v in fit_peaks.index if 'Mean_K' in v]
if iron_peak_energies is None:
iron_peak_energies=get_iron_data(len(mviron))[0]
if iron_peak_errors is None:
iron_peak_errors=np.zeros(len(iron_peak_energies))
fig = gt.make_scatter(y=mviron,
x=iron_peak_energies,
err_y=mvironerr,
err_x=iron_peak_errors,
title='Iron Peak Calibration',
showplot=False,
xlab=f'Iron Peaks Energy [{unit_iron}]',
ylab=f"Iron Peaks Amplitude [{fit_peaks['Units']['MaxX']}]",
)
root_scatter = prf.root_scatter(y=mviron,
x = iron_peak_energies,
yerr =mvironerr,
xerr=iron_peak_errors,
)
func = prf.root_func(expr='[0]*x',
name='iron_cal',
parnames={0:f"Responsivity[{fit_peaks['Units']['MinX']}/{unit_iron}]"},
parfix=parfix,
xrange=(0,max(iron_peak_energies)),
)
fit_results,func = prf.root_fitting(root_scatter,func,get_detail_matrices=get_detail_matrices)
fit_result = prf.extract_fit_values(func)
fit_result['Units']['MinX']=fit_peaks['Units']['MinX']
fit_result['Units']['MaxX']=fit_peaks['Units']['MaxX']
fig = gt.overlay_fit_func(figin=fig,fit=fit_result,func=func,funlab='Iron Calibration',showplot=False)
fig.update_yaxes(range=(0,1.05*func.Eval(iron_peak_energies[-1])))
fig = gt.make_residuals(fig,x=iron_peak_energies,y=mviron,func=func,fit_result=fit_result,err_x=iron_peak_errors,err_y=mvironerr,showplot=False)
tmp = pd.DataFrame()
tmp['Variables']=["Calibration Constant"]
tmp['Values']=[1/fit_result['Values']['Responsivity']]
tmp['Errors']=[fit_result['Errors']['Responsivity']/fit_result['Values']['Responsivity']**2]
tmp['Units']=[f"{unit_iron}/{fit_peaks['Units']['MinX']}"]
tmp.set_index('Variables',inplace=True)
fit_result = dfman.dataframe_insert(fit_result,tmp,'Chi2')
fig = gt.add_summary_legend(fig,fit_result,x=0.3,y=0.7,showlegend=False,showplot=False,printchi2=False)
if showplot:
fig.show()
return fit_result,func,root_scatter,fig