Commit d8b41f9b authored by Matheus Müller's avatar Matheus Müller

Convolui com resolucao experimental

parent 18c98e68
......@@ -157,6 +157,9 @@ def ewin_build(window, xmax, xstep):
Entrysymbol.pack()
EntryLS.pack()
##############################################################################
# Buttons
but = tk.Button(EPframel2, command=create, text='Add element', bd=1, height=1,width=15)
but.pack()#side='left')
......
......@@ -6,13 +6,13 @@ from numerical import *
import sys
##############################################################################
# Initial parameters
# Arbitrary parameters to declare dictionaries
param = dict(dedx=192. ,Theta_out=70., dW2dx=20740., E0= 100000., Theta_in=0., FWHM0=180.)
ionb = dict(Z=1, mass=1.0079)
init=0
##############################################################################
# Report file, writes parameters
def config_write():
config = open("config.cfg", 'w+')
config.write( (Entrydedx.get())+'\n' )
......@@ -50,7 +50,7 @@ def init_calc():
for i in range(len(Edict)):
Edict[i]['dist'] = np.loadtxt(Edict[i]['symbol']+".prof")
spectro(param, str(modelvar.get()), ionb, float(EntryEmin.get()),float(EntryEmax.get()), float(EntryEstep.get()), float(EntryDstep.get()), int(LSvar.get()) )
spectro(param, str(modelvar.get()), ionb, float(EntryEmin.get()),float(EntryEmax.get()), float(EntryEstep.get()), float(EntryDstep.get()), int(LSvar.get()), float(EntryFWHM0.get()) )
##############################################################################
# Main window
......@@ -207,8 +207,9 @@ ExitButton.pack(side='bottom')
ewin_build(Label1, 15, float(EntryDstep.get()))
##############################################################################
#spectro(param, 'gaussiana', ionb, float(EntryEmin.get()), float(EntryEmax.get()), float(EntryEstep.get()), float(EntryDstep.get()), int(LSvar.get()) )
if init == 0:
CalcButton.invoke()
init = 1
root.mainloop()
Label1.mainloop()
......
......@@ -39,29 +39,51 @@ def convoluiF0(espectro, alpha,passoE):
#########################################################
######## Função de convolução com Gaussiana #######
def convoluiGauss(espectro, Sigma,passoE):
areaantes=espectro.sum()
temp=zeros(len(espectro),float)
wg4_step=4.0*Sigma/passoE;
wg2=Sigma**2
eixo=arange(-6*Sigma,6*Sigma,passoE)
eixo=arange(0,len(temp)*passoE,passoE)
print len(eixo)
for i in range(0,len(eixo)):
for j in range(0,i):
temp[i]+=espectro[j]*exp(-1.0*pow(passoE*(i-j-wg4_step),2.0)/(wg2*2.0));
return temp*areaantes/temp.sum(), wg4_step*passoE
# return temp/temp.max(), wg4_step*passoE # Incluido para teste com Gaussiana -> Altura = 1
#########################################################
#### Parâmetros iniciais ####
# Feixe:
E0=100000 #eV
E0=99000 #eV
mi=1.
Zi=1.
# Alvo:
dedx=192. # eV/nm
dw2dx=20740. # eV^2/nm
Theta_in=0.
Theta_out=70.
dedx=230. # eV/nm
dw2dx=19000. # eV^2/nm
Theta_in=7.
Theta_out=58.67
Theta_s = PI - (Theta_in + Theta_out)*PI/180
espessura=2
espessura=19
# Cálculo:
deltaE=20000
passoE=20
passox=0.02
passox=0.01
e = arange(E0-deltaE, E0, passoE)
totalbessel=e*0
totalbessel_ls=e*0
totalgaussiana=e*0
resolucao = 350 # Largura a meia algura da Gaussiana que representa a resolução experimental (FWHM) em eV
## TODO: - Otimizar os cálculos de cada componente: evitar cálculos (de zeros) em energias desnecessários
## - Definir o intervalo de plotagem do espectro independente do deltaE
......@@ -69,7 +91,9 @@ totalgaussiana=e*0
## Define a composição da amostra ##
elementos=[[72.,178.49,217.],[14.,28.,121.],[8.,16.,102]] # Hf, Si e O
elementos=[[57.,138.9,192.,"la.prof"],[38.,87.62,165.,"sr.prof"],[22.,47.867,142.,"sr.prof"]] # La, Sr,ti
#elementos=[[72.,178.49,217.,"hf.txt"],[14.,28.,121.,"si.txt"],[8.,16.,102,"o.txt"]] # Hf, Si e O
#elementos=[[72.,178.49,17.],[14.,28.,21.],[8.,16.,12]] # Hf, Si e O
#elementos=[[72.,178.49,217.]]#,[14.,28.,121.],[8.,16.,102]] # Hf, Si e O
# Esta lista deve incluir a distribuição
......@@ -82,6 +106,8 @@ for componente in elementos:
mt=componente[1];
Zt=componente[0];
alpha_lineshape=componente[2]
perfil = loadtxt(componente[3], unpack=True)
# plot(perfil)
k = ((sqrt (mt**2+mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
dedxef=dedx*(k/cos(Theta_in*PI/180.) + 1/cos(Theta_out*PI/180.))
......@@ -108,7 +134,7 @@ for componente in elementos:
#
# CSRf(Z1, Z2, Ein, M1, M2, Theta):
secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox
secchoque=CSRf(Zi, Zt,k*E0-dedxef*x,mi,mt, Theta_s)*passox*perfil[int(x/passox)]
# É multiplicada pelo passo para manter a escala do gráfico
# independente do mesmo.
......@@ -132,7 +158,7 @@ for componente in elementos:
# Gaussiana:
Em=k*E0-x*dedxef
sigma=x*dw2dxef
sigma=(sqrt(x*dw2dxef)+resolucao/2.35482)**2
gaussianacamada = exp((-(e-Em)**2)/(2.*sigma))/sqrt(2.*PI*sigma)
sinalgaussiana=gaussianacamada*secchoque+sinalgaussiana
#plt.plot(e,gaussianacamada)
......@@ -146,19 +172,56 @@ for componente in elementos:
totalbessel=sinalbessel+totalbessel
totalbessel_ls=sinalbessel_ls+totalbessel_ls
totalbessel_ls_gauss,gshift=convoluiGauss(totalbessel_ls,resolucao/2.35482,passoE)
totalgaussiana=sinalgaussiana+totalgaussiana
#####################################
# Mostra o Gráfico
bes=plt.plot(e,totalbessel,lw=2,color="red")
besls=plt.plot(e,totalbessel_ls,lw=2,color="blue")
#bes=plt.plot(e,totalbessel,lw=2,color="red")
#besls=plt.plot(e,totalbessel_ls,lw=2,color="blue")
beslsgaus=plt.plot(e-gshift,totalbessel_ls_gauss,lw=2,color="orange")
gaus=plt.plot(e,totalgaussiana,lw=2,color="green")
#plt.plot(e,totalbessel*10,lw=2, color="blue")
#plt.plot(e,totalgaussiana*10,lw=2, color="green")
plt.legend([bes,besls,gaus],["Bessel","Bessel+lineshape","Gaussiana"],loc= "upper left", shadow=True)
#plt.legend([bes,besls,beslsgaus,gaus],["Bessel","Bessel+lineshape","Bessel+ls+expres","Gaussiana"],loc= "upper left", shadow=True)
plt.legend([beslsgaus,gaus],["Bessel+ls+expres","Gaussiana"],loc= "upper left", shadow=True)
#plt.legend([bes,besls,gaus],["Bessel","Bessel+lineshape","Gaussiana"],loc= "upper left", shadow=True)
plt.xlabel("Energia dos Ions (eV)")
plt.ylabel("Rendimento (contagens)")
#print totalbessel
################################################
# Testando convolução com deltas de Dirac
#figure()
#xxy=e*0
#xxy[00]=1. # Correspondente a uma delta em 80000
#xxy[100]=1. # Correspondente a uma delta em 82000
#xxy[300]=1.
#xxy[400]=1.
#xxy[500]=1.
# Convoluindo as deltas com uma gaussiana de FWHM = 550 eV
#xxyz,g_shift=convoluiGauss(xxy,resolucao/2.35482,passoE)
# Convoluindo a convolução novamente -> FWHM Res = 777.81 eV = sqrt(550^2+550^2)
#xxyz2,g_shift2=convoluiGauss(xxyz,resolucao/2.35482,passoE)
#deltas = plot(e,xxy,lw=2)
#convo1 = plot(e-g_shift,xxyz+0.001,lw=2)
#convo2 = plot(e-g_shift-g_shift2,xxyz2+0.002,lw=2)
#convo2_analitica=plot(e,exp(-(e-82000)**2.0/(2.*(777.81/2.35482)**2))+0.002,lw=2)
#plt.legend([deltas,convo1,convo2,convo2_analitica],["Deltas","Conv Gauss 1x 550 eV","Conv Gauss 2x 550 eV","Conv Analitica (777.81)"],loc= "upper right", shadow=True)
#print e-g_shift, xxyz
#print str(g_shift) + " <- g_shift"
#print len(e)
#print len(xxy)
#print len(xxyz)
show()
##################
......@@ -40,10 +40,46 @@ def convoluiF0(espectro, alpha,passoE):
return temp[::-1]*areaantes
####################################################################################################################################################
def spectro(p, modelo, ion, Emin, Emax, EPasso, passox, LScontrol):
# Função de convolução com Gaussiana
def convoluiGauss(espectro, Sigma,passoE):
areaantes=espectro.sum()
temp=zeros(len(espectro),float)
wg4_step=4.0*Sigma/passoE;
wg2=Sigma**2
eixo=arange(-6*Sigma,6*Sigma,passoE)
eixo=arange(0,len(temp)*passoE,passoE)
for i in range(0,len(eixo)):
for j in range(0,i):
temp[i]+=espectro[j]*exp(-1.0*pow(passoE*(i-j-wg4_step),2.0)/(wg2*2.0));
return temp*areaantes/temp.sum(), wg4_step*passoE
####################################################################################################################################################
def spectro(p, modelo, ion, Emin, Emax, EPasso, passox, LScontrol, resol):
e = arange(Emin, Emax, EPasso)
Y = e*0.
if modelo == 'bessel':
Y = e*0.
try:
YY
except NameError:
YY, = plot(e,e*0., 'b')
else:
YY.set_xdata(e)
elif modelo == 'gaussiana':
X = e*0.
try:
XX
except NameError:
XX, = plot(e,e*0., 'r')
else:
XX.set_xdata(e)
dose = 2e34
Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180
......@@ -98,12 +134,20 @@ def spectro(p, modelo, ion, Emin, Emax, EPasso, passox, LScontrol):
if LScontrol == 1:
besself = convoluiF0(besself, alpha_lineshape, EPasso)
Y = besself + Y
Y, gshiftY = convoluiGauss(Y,resol/2.35482,EPasso)
elif modelo == 'gaussiana':
if LScontrol == 1:
sinalgaussiana = convoluiF0(sinalgaussiana, alpha_lineshape, EPasso)
Y = sinalgaussiana + Y
plt.plot(e,Y*dose)
X = sinalgaussiana + X
X, gshiftX = convoluiGauss(X,resol/2.35482,EPasso)
if modelo == 'bessel':
YY.set_ydata(Y*dose)
ylim([0,max(Y*dose)*1.1])
elif modelo == 'gaussiana':
XX.set_ydata(X*dose)
ylim([0,max(X*dose)*1.1])
#plt.plot(e,Y*dose)
suptitle('Signal received', fontsize=12)
xlabel("Enegy (eV)")
ylabel("Yield")
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment