Lesson: 標準微中子震盪, 本程式可用以計算不同初始比例與能量的微中子, 在經過不同L/E後的比例

result










source Code Pure_NeutrinoOsci.py

#copyright: TsungChe Liu, tcliu@ntu.edu.tw
import numpy as np
import matplotlib.pyplot as plt

def transitionProb(nu_alpha, nu_beta, L, E):# L:Km, Eng:GeV
	a = nu_alpha
	b = nu_beta
	if a > 2 or b > 2 or a < 0 or b < 0:
		return 0
	prob = 0
	if(a==b):
		prob = 1;
	for j in range(3):
		for i in range(j):
			oscPartReal = pow(np.sin(1.27*dmMatric[i][j]*L/E), 2)
			oscPartIm = np.sin(2.54*dmMatric[i][j]*L/E)
			prob -= 4*np.real(np.conjugate(U[a][i])*U[b][i]*U[a][j]*np.conjugate(U[b][j]))*oscPartReal
			prob += 2*np.imag(np.conjugate(U[a][i])*U[b][i]*U[a][j]*np.conjugate(U[b][j]))*oscPartIm
	return prob


def main(source, upperL, energy):
	E = energy#GeV
	lengthTravel= range(0, upperL, 1)#km
	flavorsR = np.zeros((upperL, 3))
	for L in lengthTravel:
		proM=np.array([[transitionProb (0,0,L,E),  -transitionProb (0,1,L,E), -transitionProb (0,2,L,E)],
						[-transitionProb (1,0,L,E), transitionProb (1,1,L,E), -transitionProb (1,2,L,E)],
						[-transitionProb (2,0,L,E),-transitionProb (2,1,L,E),  transitionProb (2,2,L,E)]])
		flavorsR[L] = np.matmul(source, np.absolute(proM))
	plt.figure(1)
	plt.plot(flavorsR[:, 0])
	plt.plot(flavorsR[:, 1])
	plt.plot(flavorsR[:, 2])
	plt.plot(flavorsR[:, 0]+flavorsR[:, 1]+flavorsR[:, 2])
	plt.title('Initial flavor ratio is [0.333333, 0.666667,0.] (pion source) ')
	plt.legend(['e flavor','u flavor','tau flabor',' total'])
	plt.xlabel('L/E (km/GeV)')
	plt.ylabel('Flavor Ratio')
	plt.show()


if __name__ == '__main__':
	upperL = 35000 #km
	energy = 1 #GeV 

	dm12 = 7.59E-5
	dm23 = 2.32E-3
	dm13 = dm12+dm23
	dmMatric=np.array(
				[[0,    dm12, dm13],
		         [dm12, 0,    dm23],
		         [dm13, dm23, 0]])
	sqs12 = pow(np.sin(np.arcsin(np.sqrt(0.861))/2),2)#	sqs12 = 0.304
	sqs13 = pow(np.sin(np.arcsin(np.sqrt(0.1))/2),2)#	sqs13 = 0.02189
	sqs23 = pow(np.sin(np.arcsin(np.sqrt(0.97))/2),2)#	sqs23 = 0.579
	c12= np.sqrt(1-sqs12);
	c13= np.sqrt(1-sqs13);
	c23= np.sqrt(1-sqs23);

	delta = 0 #1.08*np.pi #0.6*np.pi
	s12 = np.sqrt(sqs12);
	s13 = np.sqrt(sqs13);
	s23 = np.sqrt(sqs23);
	cp = np.exp(delta *1j)
	cpm = np.exp(-delta *1j)

	R23 = np.array([[1,   0,   0], 
					[0, c23, s23], 
					[0,-s23, c23]])
	R13 = np.array([[c13,     0, s13*cpm], 
					[0,       1,       0], 
					[-s13*cp, 0,    c13]])
	R12 = np.array([[c12,  s12,  0], 
					[-s12, c12,  0], 
					[0,      0,  1]])
	U = np.matmul(np.matmul(R23, R13),R12)
	source = [0.333333, 0.666667,0.]
	main(source, upperL, energy) 
# type 0 : eSource = np.array([1., 0., 0. ]); #flavor muon( );
# type 0 : muonD = np.array([0., 1., 0. ]); #flavor muon( );
# type 0 : pion = np.array([0.333333, 0.666666,0.]);#flavor pion( );