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

Arquivo de testes p metodo monocamadas

parent c017ed87
......@@ -200,7 +200,7 @@ OFlabentries[13].insert(0,param['FWHM0'])
OFlabentries[14].insert(0,70000)
OFlabentries[15].insert(0,100000)
OFlabentries[16].insert(0,20)
OFlabentries[17].insert(0,0.01)
OFlabentries[17].insert(0,0.05)
OFlabentries[21].insert(0,param['E0'])
OFlabentries[22].insert(0, ionb['mass'])
OFlabentries[23].insert(0, ionb['Z'])
......@@ -273,6 +273,7 @@ def init_calc():
ionb['Z'] = float(OFlabentries[23].get())
espectro = spectro(int(float(OFlabentries[21].get())/float(OFlabentries[16].get())), int(OFlabentries[16].get()))
aux = spectro(int(float(OFlabentries[21].get())/float(OFlabentries[16].get())), int(OFlabentries[16].get()))
aux2 = spectro(int(float(OFlabentries[21].get())/float(OFlabentries[16].get())), int(OFlabentries[16].get()))
doseold = float(OFlabentries[24].get())
LSi = int(LSvar.get())
......@@ -285,7 +286,7 @@ def init_calc():
calc(param, str(modelvar.get()), ionb, OFlabentries, LSi, CGi, 'layer', espectro, OFlabentries[14].get(), Plot, f2, FT, Plotbuttons, Pcontrol)
elif int(modevar.get()) == 1:
spectrolayer(param, str(modelvar.get()), ionb, OFlabentries, LSi, CGi, 'layer', espectro, OFlabentries[14].get(), Plot, f2, FT, Plotbuttons, Pcontrol,aux)
spectrolayer(param, str(modelvar.get()), ionb, OFlabentries, LSi, CGi, 'layer', espectro, OFlabentries[14].get(), Plot, f2, FT, Plotbuttons, Pcontrol,aux,aux2)
FT = False
......
......@@ -30,7 +30,7 @@ def kinematic(mt, mi, esp):
def spect_file(x, y, shift, dose):
data2 = open('temp.sim', 'w')
for i in range(len(y)):
print>>data2, x[i]-shift,y[i]*dose
print>>data2, x[i]/1000.,y[i]
data2.close()
def set_lim(x, x2, y, shift):
......@@ -143,17 +143,16 @@ def calc(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo, espectr
####################################################################################################################################################
####################################################################################################################################################
def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo, espectro, emin, plote, fig, firsttime, botoes, indice, aux):
def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo, espectro, emin, plote, fig, firsttime, botoes, indice, aux, aux2):
EPasso, passox, e, resol, dose, Emin, Emax = variables(OFlabentries)
gshift=0.
curdepht=passox #Avoid zero division
espectro.clear()
espectro.spectro[0] = 1.
Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180
Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180.
mi = ion['mass']
first = False
deltaE = 0
for camada in Ldict:
......@@ -168,12 +167,12 @@ def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo,
dedx = float(camada.LandE[3].get())
dW2dx = float(camada.LandE[4].get())
k = ((sqrt (mt**2-mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2
dedxef=dedx*(k/cos(p['Theta_in']*PI/180.) + 1/cos(p['Theta_out']*PI/180.))
k = ((sqrt (mt**2.-mi**2.*(sin(Theta_s)**2.)) + mi*cos(Theta_s) ) / (mi+mt) )**2.
dedxef=dedx*(k/cos(p['Theta_in']*PI/180.) + 1./cos(p['Theta_out']*PI/180.))
dw2dxef=dW2dx*(k**2./cos(p['Theta_in']*PI/180.)+1./cos(p['Theta_out']*PI/180.))
sinalbessel = e*0
sinalgaussiana = e*0
sinalbessel = e*0.
sinalgaussiana = e*0.
# Calcula a contribuição de cada passo
for x in arange(curdepht, xfinal, passox):
......@@ -187,16 +186,13 @@ def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo,
alpha=dedxef*(2./dw2dxef)
m=alpha*dedxef
lbd=m*x*alpha
aux.spectro[int((k*p['E0']-deltaE)/EPasso-1)] = CS*exp(-m*x)+x*m*alpha*exp(-m*x)*i1(0)/sqrt(0.00000000001)*passox*passox
besselcamada = lbd*exp(-m*x-alpha*(k*p['E0']+deltaE-e))*i1(2.*sqrt(lbd*(k*p['E0']+deltaE-e)))/(sqrt(lbd*(k*p['E0']+deltaE-e)+.0000001))
aux.soma( list(nan_to_num(besselcamada*CS)))
if first:
espectro.Convolve(aux)
else:
espectro.soma(aux.spectro)
first = True
a = ((p['E0'] - k*dedxef)/p['E0'])
aux.clear()
deltaE = deltaE + dedxef*passox
aux.spectro[0] = a*1./EPasso #[int((k*p['E0'])/EPasso-1)] = exp(-m*x)+x*m*alpha*exp(-m*x)*i1(0)/sqrt(0.00000000001)*passox*passox
#lbd*exp(-m*x-alpha*(k*p['E0']-e))*i1(2.*sqrt(lbd*(k*p['E0']-e)))/(sqrt(lbd*(k*p['E0']-e)+.0000001))
aux.soma( alpha*exp(-alpha*arange(0,10.,passox)*dedxef))
espectro.Convolve(aux)
#espectro.soma(aux2.spectro*CS)
# Gaussiana:
elif modelo == 'gaussiana':
......@@ -219,8 +215,7 @@ def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo,
espectro.ConvolveGauss(resol/2.35482)
espectro.multiply(dose)
espectro.eshift = espectro.eshift -(1-k)*p['E0'] - EPasso
print espectro.eshift
espectro.spectro[0]=0.
plote[indice].set_xdata(espectro.axisenergy())
plote[indice].set_ydata(espectro.spectro)
......
from pylab import *
from spec import *
PI=math.pi
espectro = spectro(2000, 50)
espectro2 = spectro(2000, 50)
espectro3 = spectro(2000, 50)
espectro2.spectro[0] = 1.
def CSRf(Z1, Z2, Ein, M1, M2, Theta):
if M1 < M2:
return ((Z1*Z2*4.8e-20*1.e8)/(4.*Ein))**2. * sin(Theta)**-4.
else:
return 0.
x = arange(0,1.,0.05)
k = 1.
dedx = 250.
dW2dx = 26500.
mt=178.
mi=1.
Theta_in = 0.
Theta_out = 70.
Theta_s = PI - (Theta_in + Theta_out)*PI/180.
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.))
dw2dxef=dW2dx*(k**2./cos(Theta_in*PI/180.)+1./cos(Theta_out*PI/180.))
alpha=dedxef*(2./dw2dxef)
E0 = 100000.
for j in x:
CS = k*CSRf(1., 72., k*E0-dedxef*j, mi, mt, Theta_s)*0.05
a = ((E0 - 2*k*dedxef)/E0)
espectro.clear()
espectro.soma(CS*alpha*exp(-alpha*x*dedxef))
espectro.spectro[0] = CS*a*1./50. #- espectro.norma() #espectro.spectro[0] + (1.-espectro.norma())/50.
espectro2.Convolve(espectro)
espectro3.soma(espectro2.spectro)
#plot(espectro3.axisenergy(),espectro3.reverse())
espectro2.multiply(99000.)
plot(espectro2.axisenergy(),espectro2.reverse())
amostra = np.loadtxt('monoteste.sim')
xxx= amostra[:,0]*1000
ccc=amostra[:,1]
plot (xxx, ccc, 'r.')
show()
......@@ -40,6 +40,12 @@ class spectro:
self.spectro = zeros(self.size)
self.eshift = 0.
def reverse(self):
temp = list()
for j in range(self.size):
temp.insert(j, self.spectro[int(self.size-j-1)])
return temp
def multiply(self, factor): ### Multiplica o espectro por um fator ###
for j in range(self.size):
self.spectro[j] = self.spectro[j]*factor
......
This diff is collapsed.
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