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 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. else: return 0. x = arange(0,1.,0.01) dedx = 250. dW2dx = 26500. mt=178. mi=1. Theta_in = 0. Theta_out = 70. Theta_s = PI - (Theta_in + 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 = 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.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) amostra = np.loadtxt('monoteste.sim') xxx= amostra[:,0]*1000 ccc=amostra[:,1] plot (xxx, ccc, 'r.') show()