用python拟合费米函数与高斯函数的卷积

Question

我最近在处理光电发射数据，想通过拟合费米函数与高斯函数的卷积来估计仪器的分辨率。现在，我忽略了我也必须加入背景。 我试着根据我在Link1中发现的东西来做一个苏尔特

不幸的是，我现在卡住了，请你帮忙。我提供了一个小例子，模拟数据（通过眼睛判断）实际上是缺失背景的真实数据的良好表示：

import array as ar
import math
from math import e, sqrt, pi
import pylab as plt
from numpy import array, linspace, arange, zeros, ceil, amax, amin, argmax, argmin, abs
from numpy import polyfit, polyval, seterr, trunc, mean
from numpy.linalg import norm
from scipy.interpolate import interp1d
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt

kb = 8.6173e-5 # Boltzmann k in eV/K


def f(e,scale,T=300,muf=0):
    return scale*(1.0 / (np.exp((-e - muf)/(kb*T)) + 1)) # -e due to binding energy scale

def gaussian(x, mu, sig,scale):
    return scale*(np.exp(-np.power(-x - mu, 2.) / (2 * np.power(sig, 2.))))  # -x due to binding energy scale




Points=10000 # If i change the points the "modelt data" doesn not look good anymore?
espan = np.linspace(-10, 10, Points)


# simulate the fermi edge data such that is close to the read real data
# note that it actually is a convoluiton of a fermi and a gaussian function.
Sim_data = np.random.normal((8.07+np.convolve(f(espan,scale=1.05), gaussian(espan, mu=0, sig=0.3,scale=1.05),mode="same")),10)

# min max norm not ( not used for now )
# Sim_data_norm = ((Sim_data - Sim_data.min()) / (Sim_data.max() - Sim_data.min()))


# convolution
# I know the actual data was meassured at T=300K so i preset this value in the fermi function
# The position of the fermi function and the gaussian should be 0 so mu and muf are set to 0
# i guess x_values would have to be replaced for the x axis of my real data later?
def Gaus_fermi_convolution(size,sigC,scaleC,xC=espan):
    return size* np.convolve(gaussian(x=xC,sig=sigC,scale=scaleC,mu=0), f(e=espan,scale=scaleC,T=300,muf=0), mode="same") 

# try to fit
lowerBounds = [0.0, 0.0,0.0] # A, B, C ,D lower bounds
upperBounds = [440.0, 0.6,440.0] # A, B, C, D upper bounds

popt, pcov = scipy.optimize.curve_fit(Gaus_fermi_convolution, espan, Sim_data,bounds=[lowerBounds, upperBounds])
Opt_size,Opt_sigC,Opt_scaleC = popt

# Print out the coefficients determined by the optimization
print(Opt_size,Opt_sigC,Opt_scaleC)

# Plotting 
plt.plot(espan,Sim_data)
#The result is clearly wrong
plt.plot(espan,Gaus_fermi_convolution(size=Opt_size,sigC=Opt_sigC,scaleC=Opt_scaleC,xC=espan))

plt.xlim([2, -2])
plt.xlabel( r'$Binding \: Energy  \: [eV]$')
plt.ylabel(r'$Intensity\: [arb. u.] $')
plt.show()