script.py 2.17 KB
Newer Older
1 2 3 4 5 6 7 8 9
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.

10 11 12
def k(n, m, x):
    return (m*x)**n*exp(-m*x)/math.factorial(n)  

13 14 15 16 17 18 19
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.
               

20
x = arange(0,1.,0.01)
21 22 23 24 25 26 27
dedx = 250.
dW2dx = 26500.
mt=178.
mi=1.
Theta_in = 0.
Theta_out = 70.
Theta_s = PI - (Theta_in + Theta_out)*PI/180.
28 29 30
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.))
31 32
alpha=dedxef*(2./dw2dxef)
E0 = 100000.
33
m=alpha*dedxef
34 35

for j in x:
36 37 38 39 40 41 42
    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)
43
    espectro.clear()
44 45 46 47 48 49
    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()))
50
    espectro2.Convolve(espectro)
51
    espectro2.multiply(CS)
52
    espectro3.soma(espectro2.spectro)
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    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)
68 69 70 71 72 73 74 75 76 77 78


amostra = np.loadtxt('monoteste.sim')
xxx= amostra[:,0]*1000
ccc=amostra[:,1]

plot (xxx, ccc, 'r.')



show()