0 
import numpy as np
import os
import matplotlib.pyplot as plt
from lmfit.models import GaussianModel, ExponentialModel, LorentzianModel, VoigtModel
import scipy
from scipy.optimize import leastsq
###############################################################################
###############################################################################

os.chdir('D:/reference')

def readfio(filename):
"""
READ *.fio file into dictionary of arrays
"""
    cols = []
# Open file
    f = open(filename +'.fio')
    counter = 0
    for n, line in enumerate(f):
        if line.strip().startswith('Col'):
            linestoskip = n +  1
            cols.append(line.split()[2]) #To have only the column headers in a list
            counter = counter + 1 #gives the number of columns, in principle it is not necessary
    f.close() # close the file 
    # Read numerical data without header
    #print (cols)
    data = np.genfromtxt(filename+'.fio',skip_header=linestoskip, skip_footer =1)
    return data, cols
###############################################################################
###############################################################################
def fitting():
    data, cols = readfio(filename)
    x = (2*np.pi/wavelength(energy))*(np.sin(np.radians(data[:,cols.index('om')]))+
    np.sin(np.radians(data[:,cols.index('tt')])-np.radians(data[:,cols.index('om')])))  
    y = data[:,cols.index('signalcounter_atten')]/data[:,cols.index('petra_beamcurrent')]

    exp_mod = ExponentialModel(prefix='exp_')
    pars = exp_mod.guess(y, x=x)

    gauss1  = LorentzianModel(prefix='g1_')
    pars.update( gauss1.make_params())

    pars['g1_center'].set(3.33, min=3.2, max=3.45)
    pars['g1_sigma'].set(0.02, min=0.01)
    pars['g1_amplitude'].set(2000, min=1)

    gauss2  = LorentzianModel(prefix='g2_')

    pars.update(gauss2.make_params())

    pars['g2_center'].set(3.59, min=3.45, max=3.7)
    pars['g2_sigma'].set(0.01, min=0.001)
    pars['g2_amplitude'].set(10000, min=1)

    mod = gauss1 + gauss2 + exp_mod


    init = mod.eval(pars, x=x)
    plt.semilogy(x, y)
    plt.semilogy(x, init, 'k--')

    out = mod.fit(y, pars, x=x)

    print(out.fit_report(min_correl=0.01))

    plt.semilogy(x, out.best_fit, 'r-')
    plt.show()
###############################################################################    ###############################################################################

def wavelength(energy):
    h = 6.626e-34  #joules
    c = 2.998e8    #m/sec
    eV = 1.602e-19 #Joules
    wavelength = h*c/(energy*1000*eV)*1e10       
    return wavelength

##############################################################################    
energy = 10.000 # in KeV    
filename = 'alignment_00241'
fitting()

我有这个程序来使用洛伦兹函数来拟合两个峰值。但是，我没有得到最合适的。谁能帮我这个。这是我用来拟合的文件。数据文件可以从此链接下载

4

1 回答 1

1

我相信您的 y 数据中包含 NaN，从而阻止了拟合。您可能会看到参数值的 NaN。如果是这样，请执行以下操作::

for ix in np.where(np.isnan(y)):
    y[ix]  = (y[ix-1] + y[ix+1]) / 2.0

在进行合身之前可能会有所帮助。

于 2016-05-20T02:35:20.783 回答