rasmus/nafuma
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Implement Lorentzian, PV and double peak fit

Browse source
This commit is contained in:
rasmusthog 2022-10-10 20:32:37 +02:00
parent f30426e95d
commit 4496c34fd2
1 changed files with 132 additions and 23 deletions
Download patch file Download diff file Expand all files Collapse all files

133
nafuma/xanes/calib.py
View file

@ -1044,7 +1044,7 @@ def fit_pre_edge_feautre(data: dict, options={}) -> pd.DataFrame:
default_options = { default_options = {
'remove_background': True, 'remove_background': True,
'background_model': 'exponential', 'background_model': 'exponential',
'fit_model': 'gaussian', 'peak_model': 'gaussian',
'extract_range': None, 'extract_range': None,
'extract_range_increments': [0, 0], 'extract_range_increments': [0, 0],
'background_limits': [None, None], 'background_limits': [None, None],
@ -1142,22 +1142,73 @@ def fit_pre_edge_feautre(data: dict, options={}) -> pd.DataFrame:
# Fit Gaussian # Fit Gaussian
# FIXME Should have options for Lorentzian and Pseudo-Voigt here as well. # FIXME Should have options for Lorentzian and Pseudo-Voigt here as well.
if options['fit_model'] == 'gaussian':
# FIXME Understand what this deprecation warning means and what changes should be made to make it future proof # FIXME Understand what this deprecation warning means and what changes should be made to make it future proof
warnings.filterwarnings(action='ignore', category=np.VisibleDeprecationWarning) warnings.filterwarnings(action='ignore', category=np.VisibleDeprecationWarning)
sigma_init = removed_background_df['ZapEnergy'].loc[removed_background_df[filename] == removed_background_df[filename].max()] mu_init = float(removed_background_df['ZapEnergy'].loc[removed_background_df[filename] == removed_background_df[filename].max()])
popt, pcov = curve_fit( gauss,
if options['peak_model'] == 'gaussian':
popt, pcov = curve_fit(gauss,
removed_background_df['ZapEnergy'], removed_background_df['ZapEnergy'],
removed_background_df[filename], removed_background_df[filename],
p0=[0.0005, sigma_init, 0.001], p0=[0.0005, mu_init, 0.001],
bounds=[ bounds=[
(0, sigma_init-0.002, -np.inf), (0, mu_init-0.002, -np.inf),
(0.5, sigma_init+0.002, np.inf) (0.5, mu_init+0.002, np.inf)
] ]
) )
elif options['peak_model'] == '2gaussian':
popt, pcov = curve_fit(_2gauss,
removed_background_df['ZapEnergy'],
removed_background_df[filename],
p0=[0.0005, mu_init, 0.001, 0.0005, mu_init+0.003, 0.0005],
#bounds=[
# (0, mu_init-0.002, -np.inf),
# (0.5, mu_init+0.002, np.inf)
# ]
#)
)
elif options['peak_model'] == 'lorentzian':
popt, pcov = curve_fit(lorentz,
removed_background_df['ZapEnergy'],
removed_background_df[filename],
p0=[1, mu_init, 0.001],
# bounds=[
# (mu_init-0.001, 0.00001),
# (mu_init+0.001, 0.01)
# ]
)
elif options['peak_model'] == '2lorentzian':
popt, pcov = curve_fit(_2lorentz,
removed_background_df['ZapEnergy'],
removed_background_df[filename],
p0=[1, mu_init, 0.001, 1, mu_init+0.003, 0.001],
# bounds=[
# (mu_init-0.001, 0.00001),
# (mu_init+0.001, 0.01)
# ]
)
elif options['peak_model'] == 'pseudo-voigt':
popt, pcov = curve_fit(pseudo_voigt,
removed_background_df['ZapEnergy'],
removed_background_df[filename],
p0=[1, mu_init, 0.001, 0.5]
)
centroids.append(popt) centroids.append(popt)
errors.append(pcov) errors.append(pcov)
@ -1175,24 +1226,63 @@ def fit_pre_edge_feautre(data: dict, options={}) -> pd.DataFrame:
axes[0].axvline(x=options['background_limits'][0][1], color='black', ls='--') axes[0].axvline(x=options['background_limits'][0][1], color='black', ls='--')
axes[0].axvline(x=options['background_limits'][1][0], color='black', ls='--') axes[0].axvline(x=options['background_limits'][1][0], color='black', ls='--')
axes[0].axvline(x=options['background_limits'][1][1], color='black', ls='--') axes[0].axvline(x=options['background_limits'][1][1], color='black', ls='--')
axes[0].set_title(f'{os.path.basename(filename)} - Background removal', size=20) axes[0].set_title(f'{os.path.basename(filename)} - Background removal', size=10)
axes[0].set_ylabel('Normalised x$\mu$(E)', size=20)
axes[0].set_xlabel('Energy (keV)', size=20)
peak_background.plot(x='ZapEnergy', y=filename, ax=axes[0], color='green', kind='scatter') peak_background.plot(x='ZapEnergy', y=filename, ax=axes[0], color='green', kind='scatter')
if options['xlim']: if options['xlim']:
axes[0].set_xlim(options['xlim']) axes[0].set_xlim(options['xlim'])
if options['ylim']: if options['ylim']:
axes[0].set_ylim(options['ylim']) axes[0].set_ylim(options['ylim'])
axes[0].set_ylabel('Normalised x$\mu$(E)', size=20)
axes[0].set_xlabel('Energy (keV)', size=20)
axes[0].axhline(y=0, ls='--', color='black')
# Fitted curve # Fitted curve
if options['peak_model'] == 'gaussian':
y_fit = gauss(removed_background_df['ZapEnergy'], *popt) y_fit = gauss(removed_background_df['ZapEnergy'], *popt)
components = [y_fit]
elif options['peak_model'] == 'lorentzian':
y_fit = lorentz(removed_background_df['ZapEnergy'], *popt)
components = [y_fit]
elif options['peak_model'] == 'pseudo-voigt':
y_fit = pseudo_voigt(removed_background_df['ZapEnergy'], *popt)
components = [y_fit]
elif options['peak_model'] == '2gaussian':
y_fit = _2gauss(removed_background_df['ZapEnergy'], *popt)
y_fit1 = gauss(removed_background_df['ZapEnergy'], *popt[0:3])
y_fit2 = gauss(removed_background_df['ZapEnergy'], *popt[3:6])
components = [y_fit1, y_fit2]
elif options['peak_model'] == '2lorentzian':
y_fit = _2lorentz(removed_background_df['ZapEnergy'], *popt)
y_fit1 = lorentz(removed_background_df['ZapEnergy'], *popt[0:3])
y_fit2 = lorentz(removed_background_df['ZapEnergy'], *popt[3:6])
components = [y_fit1, y_fit2]
removed_background_df.plot(x='ZapEnergy', y=filename, ax=axes[1], color='black', label='Background subtracted', kind='scatter') removed_background_df.plot(x='ZapEnergy', y=filename, ax=axes[1], color='black', label='Background subtracted', kind='scatter')
axes[1].plot(removed_background_df['ZapEnergy'], y_fit, color='red', label=f'Fit data ({options["fit_model"]})')
axes[1].set_title(f'{os.path.basename(filename)} - Pre-edge feature fit', size=20) axes[1].plot(removed_background_df['ZapEnergy'], y_fit, color='red', label=f'Fit data ({options["peak_model"]})')
for comp in components:
axes[1].fill_between(x=removed_background_df['ZapEnergy'], y1=comp, y2=0, alpha=0.2)
residuals = (removed_background_df[filename] - y_fit) - 0.1*removed_background_df[filename].max()
axes[1].scatter(x=removed_background_df['ZapEnergy'], y=residuals)
axes[1].axhline(y=-0.1*removed_background_df[filename].max(), ls='--', color='black')
axes[1].set_title(f'{os.path.basename(filename)} - Pre-edge feature fit', size=10)
axes[1].set_ylabel('Normalised x$\mu$(E)', size=20) axes[1].set_ylabel('Normalised x$\mu$(E)', size=20)
axes[1].set_xlabel('Energy (keV)', size=20) axes[1].set_xlabel('Energy (keV)', size=20)
axes[1].legend()
axes[1].axhline(y=0, ls='--', color='black')
if options['xlim']: if options['xlim']:
axes[1].set_xlim(options['xlim']) axes[1].set_xlim(options['xlim'])
@ -1220,12 +1310,31 @@ def fit_pre_edge_feautre(data: dict, options={}) -> pd.DataFrame:
def gauss(x, A, mu, sigma): def gauss(x, A, mu, sigma):
return (A/(sigma*np.sqrt(np.pi)))*np.exp(-(x-mu)**2/(2*sigma**2)) return (A/(sigma*np.sqrt(np.pi)))*np.exp(-(x-mu)**2/(2*sigma**2))
def _2gauss(x, A1, mu1, sigma1, A2, mu2, sigma2):
return (A1/(sigma1*np.sqrt(np.pi)))*np.exp(-(x-mu1)**2/(2*sigma1**2))+(A2/(sigma2*np.sqrt(np.pi)))*np.exp(-(x-mu2)**2/(2*sigma2**2))
def lorentz(x, A, mu, sigma):
return (A/np.pi * ((sigma)/(((x-mu)**2) + (sigma)**2)))
def _2lorentz(x, A1, mu1, sigma1, A2, mu2, sigma2):
return (A1/np.pi * ((sigma1)/(((x-mu1)**2) + (sigma1)**2))) + (A2/np.pi * ((sigma2)/(((x-mu2)**2) + (sigma2)**2)))
def pseudo_voigt(x, A, mu, sigma, eta):
G = gauss(x, A, mu, sigma)
L = lorentz(x, A, mu, sigma)
return eta*G + (1-eta)*L
def arctan(x,a,b,c,d): def arctan(x,a,b,c,d):
return a*np.arctan(x*b+c) + d return a*np.arctan(x*b+c) + d
def save_data(data, options={}): def save_data(data, options={}):
default_options = { default_options = {