Implementado modelo monocamadas

Lento - sugerido multithread para acelerar

calc.py

@@ -24,6 +24,11 @@ def kinematic(mt, mi, esp):
     MM = mi/mt
     return ((sqrt (1. - MM**2.*sin(esp)**2.) + MM*cos(esp) ) / (1. + MM) )**2.
 
+####################################################################################################################################################
+# Coeficientes de Poisson
+def kpoisson(n, m, x):
+    return (m*x)**n*exp(-m*x)/math.factorial(n)  
+
 ####################################################################################################################################################
 # Funcoes em comum
 
@@ -146,10 +151,13 @@ 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, aux2):
 
     EPasso, passox, e, resol, dose, Emin, Emax = variables(OFlabentries)
+    passox = 0.005
     gshift=0.
     curdepht=passox    #Avoid zero division
     espectro.clear()
-    espectro.spectro[0] = 1.
+    aux.clear()
+    aux2.clear()
+    aux2.spectro[0] = 1./EPasso
 
     Theta_s = PI - (p['Theta_in'] + p['Theta_out'])*PI/180.
     mi = ion['mass']
@@ -170,9 +178,6 @@ def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo,
             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.
             
             # Calcula a contribuição de cada passo
             for x in arange(curdepht, xfinal, passox): 
@@ -181,41 +186,31 @@ def spectrolayer(p, modelo, ion, OFlabentries, LScontrol, CGausscontrol, metodo,
                 CS = c*k*CSRf(ion['Z'], camada.Elements[componente]['Z'], k*p['E0']-dedxef*x, mi, mt, Theta_s)*passox 
                 # É multiplicada pelo passo para manter a escala do gráfico independente do mesmo.
 
-                # Bessel:  
-                if modelo == 'bessel':
-                    alpha=dedxef*(2./dw2dxef)
-                    m=alpha*dedxef
-                    lbd=m*x*alpha
-                    a = ((p['E0'] - k*dedxef)/p['E0'])
-                    aux.clear()
-                    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)
+                alpha=dedxef*(2./dw2dxef)
+                m=alpha*dedxef
 
-                # Gaussiana:
-                elif modelo == 'gaussiana':
-                    Em=k*p['E0']-x*dedxef
-                    sigma=x*dw2dxef
-                    gaussianacamada = exp((-(e-Em)**2.)/(2.*sigma))/sqrt(2.*PI*sigma) 
-                    espectro.soma( list(nan_to_num(gaussianacamada*CS)) )
-          
-            if modelo == 'bessel':
-                linecolor='r'
+                aux.clear()
+                for y in range(20):
+                    aux.soma(kpoisson((y+1.),m,passox)*alpha*exp(-alpha*aux.axisenergy()))
+                aux.spectro[0] = kpoisson(0.,m,passox)*1./EPasso
 
-            elif modelo == 'gaussiana':         
-                linecolor='b'
+                aux2.Convolve(aux)
+                aux2.multiply(CS)
+                espectro.soma(aux2.spectro)
+                aux2.multiply(1./CS)
+         
+            linecolor='b'
 
         curdepht = xfinal
 
+    espectro.reverse()
     if LScontrol == 1:
         espectro.ConvolveF0(1./alpha_lineshape)
     if CGausscontrol == 1:
         espectro.ConvolveGauss(resol/2.35482) 
 
-    espectro.multiply(dose)
-    espectro.spectro[0]=0.
+    espectro.eshift = espectro.eshift - (1-k)*p['E0'] - EPasso
+    espectro.multiply(dose*0.93)
     plote[indice].set_xdata(espectro.axisenergy())           
     plote[indice].set_ydata(espectro.spectro)
 

script.py

@@ -7,6 +7,9 @@ espectro2 = spectro(2000, 50)
 espectro3 = spectro(2000, 50)
 espectro2.spectro[0] = 1.
 
+def k(n, m, x):
+    return (m*x)**n*exp(-m*x)/math.factorial(n)  
+
 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.
@@ -14,8 +17,7 @@ def CSRf(Z1, Z2, Ein, M1, M2, Theta):
         return 0.
                
 
-x = arange(0,1.,0.05)
-k = 1.
+x = arange(0,1.,0.01)
 dedx = 250.
 dW2dx = 26500.
 mt=178.
@@ -23,26 +25,48 @@ 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.))
+kc = ((sqrt (mt**2-mi**2*(sin(Theta_s)**2)) + mi*cos(Theta_s) ) / (mi+mt) )**2	
+dedxef=dedx*(kc/cos(Theta_in*PI/180.) + 1./cos(Theta_out*PI/180.))
+dw2dxef=dW2dx*(kc**2./cos(Theta_in*PI/180.)+1./cos(Theta_out*PI/180.))
 alpha=dedxef*(2./dw2dxef)
 E0 = 100000.
-
+m=alpha*dedxef
 
 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)
+    CS = kc*CSRf(1., 72., kc*E0-dedxef*j, mi, mt, Theta_s)*0.05 
+    a = k(0.,m,0.01)
+    b = k(1.,m,0.01)
+    c = k(2.,m,0.01)
+    d = k(3.,m,0.01)
+    e = k(4.,m,0.01)
+    f = k(5.,m,0.01)
     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.
+    espectro.spectro[0] = a*1./50.
+    espectro.soma(b*alpha*exp(-alpha*espectro.axisenergy()))
+    espectro.soma(c*alpha*exp(-alpha*espectro.axisenergy()))
+    espectro.soma(d*alpha*exp(-alpha*espectro.axisenergy()))
+    espectro.soma(e*alpha*exp(-alpha*espectro.axisenergy()))
+    espectro.soma(f*alpha*exp(-alpha*espectro.axisenergy()))
     espectro2.Convolve(espectro)
+    espectro2.multiply(CS)
     espectro3.soma(espectro2.spectro)
+    espectro2.multiply(1./CS)
+
+    #print '#',int(j/0.05), ' Area     ',' Media    ','Variancia     ','\nkmda', espectro2.norma(), espectro2.firstmoment(), espectro2.secondmoment()
+    #print 'conv', espectro2.norma(), espectro2.firstmoment(), espectro2.secondmoment()
+    #print 'soma', espectro3.norma(), espectro3.firstmoment(), espectro3.secondmoment(), '\n'
+espectro3.reverse()
+espectro3.ConvolveF0(1./200.)
+espectro3.ConvolveGauss(250./2.35482) 
+
+#plot(espectro2.axisenergy(),espectro2.spectro)
+espectro3.multiply(22.e33)
+
+#espectro3.eshift = espectro3.eshift - 1600.
+
+plot(espectro3.axisenergy(),espectro3.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]

spec.py

@@ -19,7 +19,7 @@ class spectro:
         self.size = N                      ### Tamanho do Vetor                    ###
 
     def norma(self): ### Retorna a area total do espectro ###
-        return sum(self.spectro)/self.stepsize
+        return sum(self.spectro)*self.stepsize
 
     def axisenergy(self): ### Retorna o vetor com as energias analogas a self.spectro ###
         ae = arange(self.eshift,((self.size*self.stepsize)+self.eshift-1),self.stepsize)
@@ -44,7 +44,7 @@ class spectro:
         temp = list()
         for j in range(self.size):
             temp.insert(j, self.spectro[int(self.size-j-1)])
-        return temp
+        self.spectro = temp
     
     def multiply(self, factor): ### Multiplica o espectro por um fator ###
         for j in range(self.size):
@@ -93,16 +93,25 @@ class spectro:
         self.normalize()
 
     def firstmoment(self): ### Retorna o primeiro momento da distribuicao ###
-        moment=0;
-        for i in arange(0,self.size,1):
-            moment = moment + (1.*i*self.stepsize)*self.spectro[i]
-        return  moment/self.norma()
+        moment=0.;
+        i=0;
+        #for i in arange(0,self.size,1):
+        #    moment = moment + (1.*i*self.stepsize)*self.spectro[i]
+        while moment < self.norma()/2.:
+            moment = moment +self.spectro[i]*self.stepsize
+            i= i +1
+        return   i*self.stepsize
 
     def secondmoment(self): ### Retorna o segundo momento da distribuicao ###
         moment=0;
         for i in arange(0,self.size,1):
              moment = moment + (1.*i*self.stepsize)**2*self.spectro[i]
-        return  moment/self.norma()
+        #a = list()
+        #for k in range(self.size):
+        #    a.insert(k, [self.axisenergy()[k],self.spectro[k]])
+        #return  np.var(a)#
+        return moment/self.norma()
+        #print self.axisenergy()
 
     def thirdmoment(self): ### Retorna o terceiro momento da distribuicao ###
         moment=0;