forked from optics/eib700
1040 lines
31 KiB
Python
1040 lines
31 KiB
Python
#!/usr/bin/env python3
|
|
# -*- coding: utf-8 -*-
|
|
"""
|
|
19.2.2025 gaussfit_const nun mit drei Rückgabeparametern.
|
|
|
|
|
|
13.2.2025 VOIGT_AVAILABLE eingef"uhrt
|
|
class ReadAsciiData
|
|
|
|
|
|
def Fitten( fitfun, x, y, par, fix=None)
|
|
|
|
def gauss (x, Y0, Y1, Y2, sigma, x0, A)
|
|
def lorentz(x, Y0, Y1, Y2, gamma, x0, A)
|
|
def voigt (x, Y0, Y1, Y2, sigma, gamma, x0, A)
|
|
|
|
def Gauss(x, *par ) # Y0, Y1, Y2, gamma, x0, A jede gausskurve mit eigener sigma
|
|
def gauss_single(x, Y0, Y1, Y2, sigma, x0, A): # obsolet, ersetzt durch gauss bzw Gauss
|
|
|
|
def gaussfit_const(x, y) # Gausskurve mit konstantem Untergrund
|
|
def gaussfit_lin(x, y) # Gausskurve mit linearem Untergrund
|
|
|
|
def PMOS_lorentz ( x, *par)
|
|
def PMOS_lorentz_falt( x, *par)
|
|
def VOIGT_NUM ( x, *par)
|
|
def lorentzNfree ( x, *par)
|
|
|
|
def voigtFWHM(gamma,sigma, plot=None) # Bestimmt die FWHM eines Voigt-profils
|
|
def is_number(string)
|
|
|
|
def ReadPara(name):
|
|
def WritePar(name,par,fix):
|
|
|
|
def plot(x,y):
|
|
def plot2(y1,y2):
|
|
|
|
|
|
14.1.2025: Prozedur Fitten erlaubt es Parameter festzuhalten.
|
|
15.1.2025: lorentz-profile nun mit quadratischem Untergrund
|
|
|
|
@author: rolf
|
|
"""
|
|
import sys
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from scipy.optimize import curve_fit
|
|
|
|
VOIGT_AVAILABLE=False;
|
|
|
|
def test(x, Y0, Y1, sigma, x0, A):
|
|
print("x=",x, " Y0=",Y0, " Y1=",Y1, "sigma = ",sigma, " x0 = ",x0," A=",A)
|
|
|
|
return Y0+ Y1*x+ A*np.exp(-(x-x0)**2/(2*sigma**2))
|
|
|
|
|
|
# 19.2.2025 Nun mit drei Rückgabeparametern
|
|
def gaussfit_const(x, y):
|
|
# gauss (x, Y0, Y1, Y2, sigma, x0, A)
|
|
|
|
sum_p = y.sum() # Berechne Mittelwert und rms
|
|
dum = x*y
|
|
sum_ep = dum.sum()
|
|
dum = dum*x
|
|
sum_epp = dum.sum()
|
|
par=np.zeros(6)
|
|
|
|
par[0] = 0 # Untergrund
|
|
par[1] = 0 # steigung
|
|
par[2] = 0 # Krümmung
|
|
|
|
par[4] = sum_ep/sum_p # Schwerpunkt
|
|
par[3] = np.sqrt( sum_epp/sum_p - par[4]*par[4]) # sigma
|
|
par[3] = par[3] / np.sqrt(2) # 12.2.2023 geht besser
|
|
par[5] = y.max()
|
|
|
|
fix=np.array([0,1,1,0,0,0])
|
|
#print(par)
|
|
Popt, Done=Fitten(gauss, x, y, par, fix)
|
|
|
|
print("gaussfit : Popt =",Popt)
|
|
# Popt, pcov = curve_fit(gauss, x, y, p0=par)
|
|
# print("gaussfit : Popt =",Popt)
|
|
|
|
return Popt, par, Done
|
|
|
|
|
|
def gaussfit_lin(x, y):
|
|
# gauss (x, Y0, Y1, Y2, sigma, x0, A)
|
|
|
|
sum_p = y.sum() # Berechne Mittelwert und rms
|
|
dum = x*y
|
|
sum_ep = dum.sum()
|
|
dum = dum*x
|
|
sum_epp = dum.sum()
|
|
par=np.zeros(6)
|
|
|
|
par[0] = 0 # Untergrund
|
|
par[1] = 0 # steigung
|
|
par[2] = 0 # Krümmung
|
|
|
|
par[4] = sum_ep/sum_p # Schwerpunkt ist bei 0.
|
|
par[3] = np.sqrt( sum_epp/sum_p - par[4]*par[4]) # sigma
|
|
par[3] = par[3] / np.sqrt(2) # 12.2.2023 geht besser
|
|
par[5] = y.max()
|
|
|
|
fix=np.array([0,0,1,0,0,0])
|
|
#print(par)
|
|
Popt, Done=Fitten(gauss, x, y, par, fix)
|
|
|
|
print("gaussfit : Popt =",Popt)
|
|
# Popt, pcov = curve_fit(gauss, x, y, p0=par)
|
|
# print("gaussfit : Popt =",Popt)
|
|
|
|
return Popt, par
|
|
|
|
|
|
if (VOIGT_AVAILABLE):
|
|
from scipy.special import voigt_profile
|
|
|
|
def voigt_single(x, Y0, Y1, Y2, sigma, gamma, x0, A):
|
|
sigma = np.abs(sigma)
|
|
gamma = np.abs(gamma)
|
|
Norm= voigt_profile( 0 ,sigma,gamma)
|
|
Profile= A*voigt_profile( (x-x0),sigma,gamma)
|
|
ReturnValue=Y0+ Y1*x + Y2*x*x + Profile / Norm
|
|
return( ReturnValue)
|
|
|
|
########### N Voigt-profile mit gleichen Werten für gamma und sigm
|
|
|
|
def voigt(x,*par):
|
|
dim = len(par)
|
|
# print("len=",dim);
|
|
N = int((dim-5) / 2)
|
|
Y0 = par[0]
|
|
Y1 = par[1]
|
|
Y2 = par[2]
|
|
sigma = par[3]
|
|
gamma = par[4]
|
|
|
|
sigma=np.abs(sigma)
|
|
|
|
NN= np.arange(N) # 0 .. N-1
|
|
x0= np.zeros(N)
|
|
A = np.zeros(N)
|
|
# print("sigma = {0:9.2f}".format(sigma))
|
|
for i in NN:
|
|
x0[i] = par[5 + 2*i ]
|
|
A [i] = par[5 + 2*i+1]
|
|
# print("i= " ,i," x[{0:1d}] = {1:9.2f}".format(i,x0[i]))
|
|
|
|
V = Y0 + Y1*x + Y2*x*x
|
|
for i in NN:
|
|
V = V + voigt_single(x,0,0,0,sigma,gamma,x0[i],A[i])
|
|
return(V)
|
|
|
|
################### Bestimmt die FWHM einer Voigt-funktion ####################3
|
|
# Die Linei muss symmetreisch um 0 liegen.
|
|
def voigtFWHM(gamma,sigma, plot=None):
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
width = gamma + sigma
|
|
x=np.arange(0,5*width,0.01*width)
|
|
y=voigt(x,0,0,0,sigma,gamma,0,1)
|
|
i1=np.where (y<0.5)[0][0]
|
|
i2=np.where (y>0.5)[0][-1]
|
|
FWHM = x[i1]+x[i2]
|
|
if (plot==True):
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
ax1.plot(x,y)
|
|
ax1.set_xlim( left=0)
|
|
ax1.set_ylim( top=1, bottom=0)
|
|
ax1.grid()
|
|
ax1.hlines(y=0.5 , xmin=0, xmax=FWHM/2, linewidth=2, color='grey')
|
|
ax1.vlines(x=FWHM/2 , ymin=0, ymax=1/2, linewidth=2, color='grey')
|
|
ax1.text(FWHM/2, 0.5,"HWHM = {0:5.2f}".format(FWHM/2) ,fontsize=12, color='grey' , transform=ax1.transData, va='bottom',ha='left')
|
|
plt.show()
|
|
return (FWHM)
|
|
|
|
#
|
|
|
|
#def gauss(x, *pars):
|
|
# Y0 =pars[0]
|
|
# x0 =pars[1]
|
|
# sigma=pars[2]
|
|
# A =pars[3]
|
|
# return Y0+ A*np.exp(-(x-x0)**2/(2*sigma**2))
|
|
def gauss_single(x, Y0, Y1, Y2, sigma, x0, A):
|
|
return Y0+ Y1*x+ Y2*x*x + A*np.exp(-(x-x0)**2/(2*sigma**2))
|
|
|
|
|
|
def gauss(x, *par ): # p[0] + p[1]*x + p[2]*x*x + p[5] * exp[ - ( x - p[4])^2 /( 2 p[3])^2 ]
|
|
dim = len(par)
|
|
# print("len=",dim);
|
|
N = int((dim-4) / 2)
|
|
if ((dim-4) % 2 !=0):
|
|
print("\n gauss: wrong parameter number ",dim)
|
|
print("\n gauss: p= ",par, " => exit program\n")
|
|
sys.exit()
|
|
|
|
#print(N, " Gauss lines");
|
|
NN= np.arange(N) # 0 .. N-1
|
|
x0= np.zeros(N)
|
|
A = np.zeros(N)
|
|
|
|
Y0 = par[0]
|
|
Y1 = par[1]
|
|
Y2 = par[2]
|
|
sigma = par[3]
|
|
#print("sigma = {0:9.2f}".format(sigma))
|
|
for i in NN:
|
|
x0[i] = par[4 + 2*i ]
|
|
A [i] = par[4+ 2*i+1]
|
|
# print("i= " ,i," x[{0:1d}] = {1:9.2f}".format(i,x0[i]))
|
|
|
|
G = Y0 + Y1*x + Y2*x*x
|
|
for i in NN:
|
|
G = G + gauss_single(x,0,0,0,sigma,x0[i],A[i])
|
|
return G
|
|
|
|
|
|
def Gauss(x, *par ): # p[0] + p[1]*x + + p[2]*x*x + p[5] * exp[ - ( x - p[4])^2 /( 2 p[3])^2 ]
|
|
|
|
# jede gausskurve mit unterschiedlichem sigma
|
|
dim = len(par)
|
|
print("len=",dim);
|
|
N = int((dim-3) / 3)
|
|
if ((dim-3) % 3 !=0):
|
|
print("\n gauss: wrong parameter number ",dim)
|
|
print("\n gauss: p= ",par, " => exit program\n")
|
|
sys.exit()
|
|
|
|
#print(N, " Gauss lines");
|
|
NN= np.arange(N) # 0 .. N-1
|
|
sigma= np.zeros(N)
|
|
x0 = np.zeros(N)
|
|
A = np.zeros(N)
|
|
|
|
Y0 = par[0]
|
|
Y1 = par[1]
|
|
Y2 = par[2]
|
|
#print("sigma = {0:9.2f}".format(sigma))
|
|
for i in NN:
|
|
sigma[i] = par[3 + 3*i ]
|
|
x0[i] = par[3 + 3*i+1]
|
|
A [i] = par[3 + 2*i+2]
|
|
|
|
print("i= " ,i," sigma= {1:9.2f} x= {2:9.2f} A = {3:9.2f} ".format(i, sigma[i] , x0[i], A [i] ))
|
|
|
|
G = Y0 + Y1*x + Y2*x*x
|
|
for i in NN:
|
|
G = G + gauss_single(x,0,0,0,sigma[i],x0[i],A[i])
|
|
return G
|
|
|
|
|
|
def lorentz_single(x, Y0, Y1, Y2, gamma, x0, A):
|
|
return Y0+ Y1*x+ Y2*x*x+ A* gamma**2/( (x-x0)**2 + gamma**2)
|
|
|
|
def lorentz_old(x, *par ):
|
|
dim = len(par)
|
|
print("lorentz len =",dim);
|
|
N = int((dim-4) / 2)
|
|
print("lorentz N =",N)
|
|
|
|
NN= np.arange(N) # 0 .. N-1
|
|
x0= np.zeros(N)
|
|
A = np.zeros(N)
|
|
|
|
Y0 = par[0]
|
|
Y1 = par[1]
|
|
Y2 = par[2]
|
|
|
|
gamma = par[3]
|
|
#print(N, " Lorentz lines gamma=",gamma);
|
|
|
|
for i in NN:
|
|
x0[i] = par[4 + 2*i ]
|
|
A [i] = par[4+ 2*i+1]
|
|
#print("x[{0:1d}] = {1:9.2f}, A[{0:1f}] ={2:9.2f}".format(i,x0[i],A[i]))
|
|
|
|
G = Y0 + Y1*x + Y2*x*x
|
|
|
|
for i in NN:
|
|
G = G + lorentz_single(x,0,0,0,gamma, x0[i], A[i])
|
|
|
|
return G
|
|
|
|
|
|
def lorentz(x, *par ):
|
|
Y0 = 0; Y1=0; Y2=0;
|
|
|
|
dim = len(par); # print("lorentz len =",dim);
|
|
|
|
if (dim==0): return 0;
|
|
if (dim>0): Y0=par[0]
|
|
if (dim>1): Y1=par[1]
|
|
if (dim>2): Y2=par[2]
|
|
|
|
G = Y0 + Y1*x + Y2*x*x
|
|
|
|
if (dim < 4):
|
|
return G
|
|
else:
|
|
N = int((dim-4) / 2)
|
|
# print("lorentz N =",N)
|
|
|
|
NN= np.arange(N) # 0 .. N-1
|
|
x0= np.zeros(N)
|
|
A = np.zeros(N)
|
|
|
|
gamma = par[3]
|
|
#print(N, " Lorentz lines gamma=",gamma);
|
|
|
|
for i in NN:
|
|
x0[i] = par[4 + 2*i ]
|
|
A [i] = par[4+ 2*i+1]
|
|
#print("x[{0:1d}] = {1:9.2f}, A[{0:1f}] ={2:9.2f}".format(i,x0[i],A[i]))
|
|
G = G + lorentz_single(x,0,0,0,gamma, x0[i], A[i])
|
|
|
|
return G
|
|
|
|
# Alle Linien mit unterschiedlicher Breite, quadratischer Untergrund
|
|
def lorentzNfree(x, *par ):
|
|
dim = len(par)
|
|
# print("len=",dim);
|
|
N = int((dim-3) / 3)
|
|
|
|
|
|
NN= np.arange(N) # 0 .. N-1
|
|
x0= np.zeros(N)
|
|
A = np.zeros(N)
|
|
gamma = np.zeros(N)
|
|
|
|
Y0 = par[0]
|
|
Y1 = par[1]
|
|
Y2 = par[2]
|
|
Y2=0
|
|
G = Y0 + Y1*x + Y2*x*x
|
|
|
|
for i in NN:
|
|
x0[i] = par[3 + 3*i ]
|
|
gamma[i]= par[3 + 3*i+1]
|
|
A [i] = par[3+ 3*i+2]
|
|
G = G + lorentz(x,0,0, gamma[i], x0[i], A[i])
|
|
#print("x[{0:1d}] = {1:9.2f}, A[{0:1f}] ={2:9.2f}".format(i,x0[i],A[i]))
|
|
|
|
|
|
return G
|
|
|
|
# wozu habe ich das?
|
|
# Aufruf z.B. : singlefit(roi1)
|
|
#Line =[]
|
|
#Parabel=[]
|
|
def singlefit(r):
|
|
i= int ( (r.ny)/2 )
|
|
x = r.yy[:,0] # y-koordinaten des vertikalen Schnittes,
|
|
xx = r.xx
|
|
|
|
xs= []; ys=[];
|
|
|
|
for i in np.arange(0,r.nx,10) :
|
|
print("i=",i)
|
|
y = r.roi[:,i] # vertikaler Schnitt bei x = i
|
|
p1= gaussfit(x, y)
|
|
|
|
ys = np.append(ys, p1[2] ) # Schwerpunkt
|
|
xs = np.append(xs, xx[0,i] ) # Schwerpunkt
|
|
|
|
print("\nResult of gaussfit: ", p1)
|
|
LinePar=np.polyfit(xs,ys,1)
|
|
ParaPar=np.polyfit(xs,ys,2)
|
|
Line= np.polyval(LinePar,xs)
|
|
Para= np.polyval(ParaPar,xs)
|
|
|
|
|
|
|
|
def WrappedFitFunc(x, *par):
|
|
global Ind; # Index der anzeigt welcher Parameter das ist
|
|
global FITFUNC
|
|
global PAR #
|
|
|
|
# print("WrappedFitFunc: *par= ",*par)
|
|
# print("WrappedFitFunc: called with ", len(par));
|
|
# for i in np.arange(0,len(par) ,1):
|
|
# print("par {0} = {1} ".format(i,par[i]) );
|
|
|
|
# for i in np.arange(0,PAR.size,1):
|
|
# print("WrappedFitFunc: PAR {0} = {1} ".format(i,PAR[i]));
|
|
|
|
fitpar = 1.0*PAR; # Parameter mit denen die tatsächliche Function aufgerufen wird
|
|
for i in np.arange(0,Ind.size,1):
|
|
fitpar[Ind[i]] = par[i]
|
|
|
|
# print(" Parameter to call the real fitfunc")
|
|
# for i in np.arange(0,fitpar.size,1):
|
|
# print("WrappedFitFunc: fitpar {0} = {1} ".format(i,fitpar[i]));
|
|
|
|
# print("call fitfun with ",fitpar)
|
|
res=FITFUNC(x,*fitpar)
|
|
|
|
# print("WrappedFitFunc: Done")
|
|
|
|
return res
|
|
|
|
|
|
|
|
def TestFit():
|
|
p=2.5*np.arange(0,5,1)
|
|
fix=np.zeros(p.size)
|
|
fix[3]=1
|
|
x=np.array([4]); y=[9]
|
|
Fitten(test,x,y, p,fix)
|
|
|
|
|
|
|
|
class ReadAsciiData:
|
|
def __init__(self, name):
|
|
print, name
|
|
self.name=name
|
|
if (name == ""):
|
|
self.filename=""
|
|
self.x =[0.0,1.0]
|
|
self.z =[0.0,1.0]
|
|
else:
|
|
self.filename=name
|
|
print('ReadAsciiProfile: Filename=', self.filename )
|
|
f = open(self.filename,'r')
|
|
self.x = []
|
|
self.z = []
|
|
self.dim=1
|
|
line = f.readline()
|
|
line = line.strip()
|
|
while line:
|
|
line =line.replace("D", "E")
|
|
columns = line.split()
|
|
l0 = columns[0].strip()
|
|
l1 = columns[1].strip()
|
|
# print(line," ",l0, " ", l1 )
|
|
if ( is_number(l0) and is_number(l1) ) :
|
|
x = float(l0)
|
|
z = float(l1)
|
|
self.x.append(x)
|
|
self.z.append(z)
|
|
# print("append",u," ",h,"\n")
|
|
|
|
#else:
|
|
# print("Zeile enthaelt keine Zahlen")
|
|
line = f.readline()
|
|
line = line.rstrip()
|
|
line = line.strip()
|
|
f.close() # close file
|
|
#---- Mache Felder aus den Listen
|
|
self.x=np.array(self.x)
|
|
self.z=np.array(self.z)
|
|
# print("exit")
|
|
|
|
def is_number(string):
|
|
try:
|
|
float(string)
|
|
return True
|
|
except ValueError:
|
|
return False
|
|
|
|
|
|
def ReadPara(name):
|
|
print('ReadParameter: Filename=',name )
|
|
par = []
|
|
fix = []
|
|
try:
|
|
f = open(name,'r')
|
|
except:
|
|
print("\nFile not found : ",name)
|
|
sys.exit()
|
|
|
|
while (True):
|
|
line = f.readline()
|
|
|
|
if len(line) == 0:
|
|
# print("End of file")
|
|
f.close() # close file
|
|
par=np.array(par)
|
|
fix=np.array(fix)
|
|
|
|
# print ("done ", par)
|
|
# print("fix = ",fix)
|
|
return par, fix
|
|
|
|
line = line.strip()
|
|
# print("line = <",line,">, len= ",len(line))
|
|
if ( len(line)!=0): # falls leere Zeile
|
|
if (line[0]!="#") :
|
|
|
|
columns = line.split()
|
|
|
|
# print( "columns = ",columns)
|
|
|
|
if (line[0]=="F"):
|
|
fix.append(1)
|
|
lp = columns[1].strip()
|
|
else:
|
|
fix.append(0)
|
|
lf="f"
|
|
lp = columns[0].strip()
|
|
|
|
if is_number(lp):
|
|
par.append(float(lp))
|
|
|
|
def WritePar(name,par,fix):
|
|
print("WritePar: Open ",name)
|
|
f = open(name,'w')
|
|
index=0
|
|
for i in fix:
|
|
if (i==1):
|
|
f.write("F {}\n".format(par[index]))
|
|
else:
|
|
f.write(" {}\n".format(par[index]))
|
|
index+=1
|
|
|
|
f.close()
|
|
|
|
|
|
#par, fix = ReadPara("test.par")
|
|
#for i in np.arange(0,len(par)-1,1):
|
|
# print( "par[",i,"] = " , par[i]," fix = ",fix[i])
|
|
|
|
def plot(x,y):
|
|
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
ax1.plot(x,y)
|
|
plt.show()
|
|
|
|
def plot2(y1,y2):
|
|
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
ax1.plot(y1, color="blue", label="y1")
|
|
ax1.plot(y2, color="red", label="y2")
|
|
ax1.legend()
|
|
ax1.grid()
|
|
plt.show()
|
|
|
|
def plot_2(x, y1,y2):
|
|
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
ax1.scatter(x,y1, color="blue", label="y1")
|
|
ax1.plot(x,y2, color="red", label="y2")
|
|
ax1.legend()
|
|
ax1.grid()
|
|
plt.show()
|
|
|
|
|
|
|
|
def Beispiel1(): # plot voigt functions from Plot1_54 and compare it with the calculation here.
|
|
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
|
|
|
|
par=[ 9.28550841702634E-0002, # Fit done with Plot1_54
|
|
-4.78363903962303E-0004,
|
|
8.32490259463868E-0007,
|
|
3.44260812426823E+0000 ,
|
|
5.13720842095014E+0000 ,
|
|
3.61689342344058E+0002,
|
|
9.30078922943656E-0001 ,
|
|
3.92041504453909E+0002 ,
|
|
8.21359979370255E-0001 ,
|
|
4.22429183438332E+0002 ,
|
|
3.51056223438879E-0001 ,
|
|
4.52093268651923E+0002 ,
|
|
1.29974026072931E-0001 ]
|
|
|
|
# par = [0,0,2,4,30,10,70,5]
|
|
|
|
voigt = ReadAsciiData("data/PSRD132V.THE") # Fit done with Plot1_54
|
|
|
|
x= np.arange(voigt.x.min(),voigt.x.max(),1)
|
|
y= voigt(x,*par)
|
|
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
ax1.scatter(voigt.x, voigt.z,color="blue", label="plot1_54")
|
|
ax1.plot(x,y, color="red", label="python")
|
|
plt.show()
|
|
|
|
|
|
|
|
############ Berechne die Funktion die gefaltet werden soll #######
|
|
# hier eine lorentzfunktion ohne Untergrund
|
|
def faltfun_lorentz(x,*par):
|
|
dim = len(par)
|
|
N = int((dim-4) / 2)
|
|
gamma = par[3]
|
|
NN = np.arange(N) # 0 .. N-1
|
|
x0 = np.zeros(N)
|
|
A = np.zeros(N)
|
|
|
|
L = 0
|
|
for i in NN:
|
|
x0[i] = par[4 + 2*i]
|
|
A [i] = par[5 + 2*i]
|
|
# print("x(",i,"=" ,x0[i], "A=",A[i])
|
|
L = L + lorentz(x,0,0,0,gamma,x0[i],A[i])
|
|
|
|
# print("faltfun {0:6.2f}= {1}".format(x,L))
|
|
return L
|
|
|
|
#
|
|
|
|
# Falte eine faltfun_lorentz mit einer Gaussfunktion
|
|
# - wird in FaltFun.py mit der eingebauten voigt-funktion verglichen
|
|
#- wird später obsolet
|
|
# addiere danach einen quadratischen Untergrund
|
|
# par = Y0,Y1,Y2, sigma, FaltFunPar)
|
|
# FaltFunPar sind die Parameter für die zu faltenede Funktion
|
|
#
|
|
# V0, Y1, Y2, sigma, gamma, x0 , A
|
|
def VOIGT_NUM(x,*par):
|
|
|
|
dim = len(par)
|
|
N = int((dim-5) / 2)
|
|
Y0 = par[0]
|
|
Y1 = par[1]
|
|
Y2 = par[2]
|
|
sigma = par[3]
|
|
sigma = np.abs(sigma)
|
|
|
|
FaltFunPar = np.delete(par, [3])
|
|
|
|
Nsig = 6 # Bereich über den integriert wird: +/- Nsig * sigma
|
|
Ni = 100+1 # Anzahl der Stützstellen
|
|
|
|
Xlim = Nsig * sigma; # Integrationsgrenze
|
|
XI = np.linspace(-Xlim,Xlim, Ni)
|
|
# print("XI = ",XI)
|
|
|
|
dxi = XI[1] - XI[0]
|
|
# print("Funpar = ", FaltFunPar," dxi = ",dxi)
|
|
|
|
Int = 0.0;
|
|
for xi in XI: # Faltungsintegral
|
|
g = gauss(xi,0,0,0,sigma,0,1)
|
|
f = faltfun_lorentz(x - xi, *FaltFunPar)
|
|
Int = Int + g*f
|
|
|
|
Norm = 1 # Normiere gaussfunktion auf Amplitude=1 : 1 oder Fläche=1 : 1/ np.sqrt( (2*np.pi) * sigma)
|
|
Norm = 1/ ( np.sqrt(2*np.pi) * sigma)
|
|
|
|
Int = Int * dxi * Norm
|
|
Untergrund = Y0 + Y1*x + Y2*x*x
|
|
|
|
value = Untergrund + Int
|
|
# print("VOIGT(",x," )= ", value)
|
|
return value
|
|
|
|
|
|
# erzeuge parameter-array das nur die freien Parameter enthält
|
|
# der array Ind zeigt an welcher Parameter das ursprünglich war
|
|
def Fitten( fitfun, x, y, par, FIX=None):
|
|
global Ind; # Index der anzeigt welcher Parameter das ist
|
|
global FITFUNC
|
|
|
|
global PAR # Speichere alle Parameter für WrappedFixFunc. Dort werden die freien Parameter hinkopiert
|
|
|
|
debug=1
|
|
if (debug): print("---------------- Fitten 0 -----------------------")
|
|
|
|
FITFUNC= fitfun
|
|
Q = []
|
|
fitpar = []
|
|
PAR = 1.0*par
|
|
|
|
try:
|
|
if (FIX==None):
|
|
print("fitten is called without parameter for fix")
|
|
FIX=0.0*par
|
|
except:
|
|
print("fitten is called with parametervalues for fix")
|
|
|
|
if (len(FIX)==0): FIX=0.0*par
|
|
|
|
|
|
if (debug):
|
|
print("---------------- Fitten 1 -----------------------")
|
|
print("fix: len= ", len(FIX),": ",FIX)
|
|
print("par: len= ",len(par)," : ",par)
|
|
|
|
for i in np.arange(0,par.size,1):
|
|
print("par {0:2} = {1:10.6f} FIX = {2}".format(i,par[i], FIX[i] ));
|
|
|
|
if (debug): print("---------Fitten: shuffle varible parameter -to front ----------------------")
|
|
|
|
for i in np.arange(0,par.size,1):
|
|
if (FIX[i] == 0):
|
|
Q.append(i)
|
|
fitpar.append(par[i])
|
|
|
|
Ind = np.array(Q)
|
|
FreeParam=Ind.size
|
|
fitpar = np.array(fitpar)
|
|
|
|
print("parameter for call of curve_fit, FreeParam= ",FreeParam)
|
|
for i in np.arange(0,fitpar.size,1):
|
|
print("fitpar {0:2} = {1:10.5f} index={2} ".format(i,fitpar[i],Ind[i]));
|
|
|
|
if (debug): print("---------------- Fitten: call curve_fit -----------------------")
|
|
|
|
done=1
|
|
|
|
|
|
try:
|
|
FITPAR, pcov = curve_fit(WrappedFitFunc, x, y, fitpar)
|
|
except:
|
|
print("\n!!!!!!!!!!!!! Fit failed, !!!!!!!!!!!!!")
|
|
print ("!!!!! come back with startparameter!!!!! \n")
|
|
done=0
|
|
return par, 0
|
|
|
|
if (debug): print("---------------- Fitten: reorganize parameter-----------------------")
|
|
|
|
ReturnPar=par
|
|
for i in np.arange(0,FreeParam,1):
|
|
ReturnPar[Ind[i]] = FITPAR[i]
|
|
|
|
for i in np.arange(0,fitpar.size,1):
|
|
print("Returnpar {0:2} = {1:16f} ".format(i,ReturnPar[i]));
|
|
|
|
if (debug): print("---------------------- Fitten exit----------------")
|
|
|
|
return ReturnPar , 1
|
|
################################################################
|
|
# Fitte gesamtes PMOS signal nicht gefaltet
|
|
def FitPMOS():
|
|
from MyPMOS import MyShots
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
|
|
|
|
D=MyShots(source='file', fname='data/run0020_520_580.txt') ; MonoPar= '150 l/mm, 3rd order, PSRC132, #20' #Data
|
|
ParStart, fix = ReadPara('data/run0020_520_580_PMOS_lorentz.par')
|
|
print(ParStart)
|
|
|
|
|
|
xr= D.e; yr = D.fmean # xr,yr: eingelesene Daten
|
|
delete = np.where(xr>600)
|
|
x= np.delete(xr, delete) # Dieser Bereich wird für den Fit verwendet
|
|
y= np.delete(yr, delete)
|
|
|
|
|
|
ys=PMOS_lorentz(x,*ParStart) # Kurve für die Startparameter
|
|
|
|
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
ax1.scatter(x,y)
|
|
ax1.plot(x,ys, color="green")
|
|
|
|
|
|
Par , Done = Fitten(PMOS_lorentz, x, y, ParStart, FIX=fix)
|
|
|
|
yt=PMOS_lorentz(x,*Par) # Kurve für die Fitparameter
|
|
|
|
ax1.plot(x,yt, color="red")
|
|
|
|
plt.show()
|
|
|
|
########### Berechne die Funktion die gefaltet werden soll #######
|
|
# hier eine TransmissionsMessung nach dem GAT oihne Faltung
|
|
# par 0-3: d.. quadratischer Untergrund
|
|
# par 3-5: Gaussförmiges I0 , sigma, x0, Amplitude
|
|
# dann Extinction mit
|
|
# par 6-8: quadratischer Untergrund
|
|
# 9 : Gamma für alle Lorentzprofile
|
|
# 10,11 : Lorentzprofil .
|
|
# ...
|
|
#
|
|
def PMOS_lorentz(x,*par):
|
|
dim = len(par)
|
|
#print("\nPMOS_lorentz: len=", dim)
|
|
|
|
I0par = par[0:6]
|
|
# print("PMOS_lorentz: I0par=", I0par)
|
|
I0 = gauss(x,*I0par)
|
|
if (dim>6):
|
|
ExtPar = par[6:dim]
|
|
#print("PMOS_lorentz: extPar=", ExtPar)
|
|
|
|
extinction=lorentz(x,*ExtPar)
|
|
|
|
fun = I0 * np.exp(-1.0 * extinction)
|
|
else:
|
|
fun= I0
|
|
|
|
return fun
|
|
|
|
# Fitfunktion zum Fitten einer PMOS-Messung mit Faltung
|
|
# par 0: sigma für die Faltung
|
|
# par 1-3: d.. quadratischer Untergrund
|
|
# par 4-6: Gaussförmiges I0 , sigma, x0, Amplitude
|
|
# dann Extinction mit
|
|
# par 7-9: quadratischer Untergrund
|
|
# 10 : Gamma für alle Lorentzprofile
|
|
# 11,12 : Lorentzprofil .
|
|
# ...
|
|
# Gaussförmiges I0-Signal
|
|
# Beliebig viele Lorentzprofile in der Extinction
|
|
#
|
|
# Kopie von VOIGT_NUM
|
|
# Faltet eine Funktion f = func(x , *FaltFunPar)
|
|
# mit einer Gaussfunktion
|
|
# sigma der Gaussfunktion ist der erste Parameter.
|
|
# die Parameter für func folgen danach.
|
|
def PMOS_lorentz_falt( x,*par):
|
|
|
|
dim = len(par)
|
|
sigma = np.abs(par[0])
|
|
FunPar = par[1:dim]
|
|
|
|
|
|
Nsig = 6 # Bereich über den integriert wird: +/- Nsig * sigma
|
|
Ni = 100+1 # Anzahl der Stützstellen
|
|
|
|
Xlim = Nsig * sigma; # Integrationsgrenze
|
|
XI = np.linspace(-Xlim,Xlim, Ni) # print("XI = ",XI)
|
|
dxi = XI[1] - XI[0]
|
|
print("PMOS_lorentz_falt: Funpar len=",len(FunPar)," : ", FunPar," dxi = ",dxi)
|
|
|
|
Int = 0.0;
|
|
for xi in XI: # Faltungsintegral
|
|
g = gauss(xi,0,0,0,sigma,0,1)
|
|
f = PMOS_lorentz(x - xi, *FunPar)
|
|
Int = Int + g*f
|
|
|
|
Norm = 1 # Normiere gaussfunktion auf Amplitude=1 : 1 oder Fläche=1 : 1/ np.sqrt( (2*np.pi) * sigma)
|
|
Norm = 1/ ( np.sqrt(2*np.pi) * sigma)
|
|
|
|
Int = Int * dxi * Norm
|
|
|
|
# sys.exit()
|
|
|
|
# print("PMOS_lorentz_falt(",x," )= ", Int)
|
|
return Int
|
|
##################################################################################################################
|
|
def FitPMOS_falt():
|
|
from MyPMOS import MyShots
|
|
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator)
|
|
from pathlib import Path
|
|
from MyImage import now
|
|
|
|
from Plotstyles import SRI2024_Poster_style, PowerPoint_style, myfigure, Init, myAxis, SRI2024_Poster_style_two
|
|
Init()
|
|
SRI2024_Poster_style()
|
|
|
|
FitLineFunc=PMOS_lorentz_falt
|
|
|
|
if (FitLineFunc==gauss) : sFitLineFunc="Gauss"
|
|
elif (FitLineFunc==lorentz): sFitLineFunc="Lorentz"
|
|
elif (FitLineFunc==voigt) : sFitLineFunc="Voigt"
|
|
elif (FitLineFunc==PMOS_lorentz_falt) : sFitLineFunc="PMOS_lorentz_falt"
|
|
else : sFitLineFunc="undef"
|
|
|
|
D=MyShots(source='file', fname='data/run0020_520_580.txt') ; MonoPar= '150 l/mm, 3rd order, PSRC132, #20' #Data
|
|
runs= Path(D.filename).stem
|
|
|
|
ParStart, fix = ReadPara('data/PMOS_lorentz_falt_run20.par')
|
|
|
|
print(" \nParameter read: ",ParStart)
|
|
|
|
############################ Fit the transmission ################
|
|
|
|
xr= D.e; yr = D.fmean # xr,yr: eingelesene Daten
|
|
delete = np.where(xr>800)
|
|
|
|
x= np.delete(xr, delete) # Dieser Bereich wird für den Fit verwendet
|
|
y= np.delete(yr, delete)
|
|
|
|
|
|
ys=PMOS_lorentz_falt(x,*ParStart) # Kurve für die Startparameter
|
|
|
|
|
|
Par , Done = Fitten(PMOS_lorentz_falt, x, y, ParStart, FIX=fix)
|
|
|
|
WritePar('output/done.par',Par,fix)
|
|
|
|
yt=PMOS_lorentz_falt(x,*Par) # Kurve für die Fitparameter
|
|
|
|
|
|
###################### Evaluate results ###########
|
|
|
|
if (FitLineFunc==PMOS_lorentz_falt):
|
|
|
|
I0Par = Par[0:9]
|
|
|
|
Nlines = (Par.size-11)//2
|
|
xi = np.zeros(Nlines)
|
|
Ai = np.zeros(Nlines)
|
|
Linie = np.zeros((Nlines, len(x)))
|
|
|
|
for i in np.arange(0, Nlines):
|
|
xi[i]= Par[11+2*i]
|
|
Ai[i]= Par[12+2*i]
|
|
gamma = Par[10]
|
|
Sigma = Par[4]
|
|
x1= Par[11]
|
|
x2= Par[13]
|
|
A1= Par[14]
|
|
dEdx= 230 / (x2-x1)
|
|
e= (x-x1) *dEdx/1000 + 400.8
|
|
FWHM=dEdx*2*gamma
|
|
Monores = np.sqrt(FWHM**2 - 113**2 )
|
|
|
|
|
|
I0 = FitLineFunc(x,*I0Par)
|
|
for i in np.arange(0, Nlines):
|
|
pLine= Par[0:11]
|
|
pLine=np.append(pLine, xi[i])
|
|
pLine=np.append(pLine, Ai[i])
|
|
print("calc Line ",i," par = ",pLine)
|
|
Linie[i]=FitLineFunc(x,*pLine)
|
|
|
|
else:
|
|
print("\nExtract FWHM: FitLineFunc ", FitLineFunc," not handled -> stop")
|
|
sys.exit()
|
|
|
|
|
|
|
|
###################### Plot Transmission ###########
|
|
|
|
fig1 = myfigure()
|
|
ax1 = myAxis(fig1)
|
|
#x1.scatter(e,y)
|
|
ax1.scatter(e,y,color='blue', s=6, label='N2-data')
|
|
|
|
ax1.plot(e,ys, color="green")
|
|
ax1.plot(e,yt, color="red")
|
|
|
|
for i in np.arange(0, Nlines):
|
|
ax1.plot( e, Linie[i], color="green",lw=1,label='N2-fit')
|
|
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator)
|
|
|
|
ax1.plot (e,I0, color="black")
|
|
|
|
ax1.set_title(MonoPar)
|
|
ax1.set_ylim( bottom=-0.1,top=round(I0.max() ) )
|
|
#ax1.xaxis.set_major_formatter(FormatStrFormatter('%.0f'))
|
|
#ax1.xaxis.set_ticks(np.arange(400, 403.1, 0.5))
|
|
#ax1.xaxis.set_minor_locator(MultipleLocator(4))
|
|
ax1.minorticks_on()
|
|
ax1.tick_params(bottom=True, top=True, left=True, right=True)
|
|
ax1.tick_params(top=True, labeltop=False, bottom=True, labelbottom=True, direction="in")
|
|
plt.gca().tick_params(axis='x', which='minor', top=True)
|
|
plt.gca().tick_params(axis='y', which='minor', right=True)
|
|
ax1.set_xlabel("Photon Energy (eV)")
|
|
ax1.set_ylabel("Intensity (a.u)")
|
|
|
|
outputname="output/run20_lorentzfalt_1.png"
|
|
|
|
fig1.text(1.0,0.0,now + " "+outputname,fontsize=8 , va='bottom',ha='right')
|
|
|
|
plt.savefig(outputname)
|
|
plt.show()
|
|
print("Save to " ,outputname)
|
|
|
|
|
|
###################### Plot Extinction ###########
|
|
|
|
fig2 = myfigure()
|
|
ax2 = myAxis(fig2)
|
|
|
|
t = -1 * np.log(y/I0)
|
|
tyt = -1 * np.log(yt/I0)
|
|
ax2.scatter(e,t,color='blue', s=6, label='N2-data')
|
|
#ax1.scatter(np.arange(0,y.size,1) ,y,color='blue', s=6, label='N2-data')
|
|
|
|
ax2.plot (e,tyt, color="red")
|
|
|
|
for i in np.arange(0, Nlines):
|
|
ax2.plot( e, -1* np.log(Linie[i]/I0), color="green",lw=1,label='N2-fit')
|
|
|
|
|
|
ax2.set_xlim([399.8,402])
|
|
ax1.set_xlim( left=399.8, right=403.1)
|
|
|
|
|
|
|
|
ax2.set_title(MonoPar)
|
|
|
|
ax2.set_ylim( bottom=-0.1,top=round(2*tyt.max()+0.5)/2 )
|
|
ax2.xaxis.set_major_formatter(FormatStrFormatter('%.1f'))
|
|
ax2.xaxis.set_ticks(np.arange(400, 402.6, 0.5))
|
|
ax2.xaxis.set_minor_locator(MultipleLocator(4))
|
|
ax2.minorticks_on()
|
|
ax2.tick_params(bottom=True, top=True, left=True, right=True)
|
|
ax2.tick_params(top=True, labeltop=False, bottom=True, labelbottom=True, direction="in")
|
|
plt.gca().tick_params(axis='x', which='minor', top=True)
|
|
plt.gca().tick_params(axis='y', which='minor', right=True)
|
|
ax2.set_xlabel("Photon Energy (eV)")
|
|
ax2.set_ylabel("Extinction (a.u)")
|
|
|
|
|
|
xpos=0.55
|
|
ypos=0.93
|
|
fig2.text(xpos,ypos ,'Profile {}'.format(sFitLineFunc) , fontsize=14, transform=ax1.transAxes , va='bottom',ha='left', backgroundcolor='white')
|
|
fig2.text(xpos,ypos-0.05,'FWHM = {0:3.0f} meV '.format(FWHM) , fontsize=14, transform=ax1.transAxes , va='bottom',ha='left', backgroundcolor='white')
|
|
fig2.text(xpos,ypos-0.1, 'Nat. Linewidth = {0:3.0f} meV '.format(113) , fontsize=14, transform=ax1.transAxes , va='bottom',ha='left', backgroundcolor='white')
|
|
fig2.text(xpos,ypos-0.15,'Instrum. Contrib = {0:3.0f} meV '.format(Monores) , fontsize=14, transform=ax1.transAxes , va='bottom',ha='left', backgroundcolor='white')
|
|
fig2.text(0,1.015, 'Furka' , fontsize=20, transform=ax1.transAxes , va='bottom',ha='left', backgroundcolor='white')
|
|
|
|
|
|
|
|
outputname="output/run20_lorentzfalt_2.png"
|
|
|
|
fig2.text(1.0,0.0,now + " "+outputname,fontsize=8 , va='bottom',ha='right')
|
|
|
|
plt.savefig(outputname)
|
|
plt.show()
|
|
print("Save to " ,outputname)
|
|
|
|
|
|
#FitPMOS_falt()
|
|
#FitPMOS()
|