1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
|
import numpy as np
from PhysConst import PhysicsConstants
def eigenvectors(M):
""" Calculates the eigenvectors and eigenvalues ordered by eigenvalue size
@type M : matrix
@param M : matrix M
@rtype : list
@return : [eigenvalues list, eigenvector list]
"""
D,V = np.linalg.eig(M)
DV = []
VT = V.T
for i,eigenvalue in enumerate(D):
DV.append([eigenvalue,VT[i]])
DV = sorted(DV,key = lambda x : x[0].real)#np.abs(x[0].real))
V2 = []
D2 = []
for e in DV:
V2.append(e[1])
D2.append(e[0])
return D2,V2
# General Rotation Matrix
def R(i,j,cp,param):
""" Rotation Matrix
Calculates the R_ij rotations. Also incorporates CP-phases when necesary.
@type i : int
@param i : i-column.
@type j : int
@param j : j-row.
@type cp : int
@param cp : if cp = 0 : no CP-phase. else CP-phase = CP_array[cp]
@rtype : numpy.array
@return : returns the R_ij rotation matrix.
"""
# if cp = 0 -> no complex phase
# R_ij, i<j
if(j<i):
# no funny business
l = i
i = j
j = l
# rotation matrix generator
R = np.zeros([param.numneu,param.numneu],complex)
# diagonal terms
for k in np.arange(0,param.numneu,1):
if(k != i-1 and k != j-1):
R[k,k] = 1.0
else :
R[k,k] = param.c[i,j]
# non-diagonal terms
if(cp != 0):
sd = np.sin(param.dcp[cp])
cd = np.cos(param.dcp[cp])
faseCP = complex(cd,sd)
else:
faseCP = complex(1.0,0.0)
R[i-1,j-1] = param.s[i,j]*faseCP.conjugate()
R[j-1,i-1] = -param.s[i,j]*faseCP
return R
def calcU(param):
""" Defining the mixing matrix parametrization.
@type param : PhysicsConstants
@param param : set of physical parameters to be used.
@rtype : None
@return : Sets mixing matrix.
"""
if(param.numneu == 2):
return self.R(1,2,0,param)
elif(param.numneu == 3):
return np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param)))
elif(param.numneu == 4):
return np.dot(R(3,4,0,param),np.dot(R(2,4,2,param),np.dot(R(1,4,0,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))))))
elif(param.numneu == 5):
return np.dot(R(4,5,0,param),np.dot(R(3,5,0,param),np.dot(R(2,5,0,param),np.dot(R(1,5,3,param),np.dot(R(3,4,0,param),np.dot(R(2,4,0,param),np.dot(R(1,4,2,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))))))))))
elif(param.numneu == 6):
# 3+3 twin-sterile-neutrino model
return np.dot(R(3,6,0,param),np.dot(R(2,5,0,param),np.dot(R(1,4,0,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))))))
else:
print "Sorry, too many neutrinos. Not yet implemented! =(."
quit()
# antineutrino case
if param.neutype == "antineutrino" :
return self.U.conjugate()
def massM2(param):
""" Mass term in the neutrino mass basis.
@type param : PhysicsConstants
@param param : set of physical parameters to be used.
@rtype : numpy array
@return : mass matrix in mass basis.
"""
M2 = np.zeros([param.numneu,param.numneu],complex)
for k in np.arange(1,param.numneu,1):
M2[k,k] = param.dmsq[k+1]
return M2
def flavorM2(param):
""" Mass term in the neutrino flavor basis.
@type param : PhysicsConstants
@param param : set of physical parameters to be used.
@rtype : numpy array
@return : mass matrix in flavor basis.
"""
U = calcU(param)
return np.dot(U,np.dot(massM2(param),U.conjugate().T))
class LVP(PhysicsConstants):
def __init__(self):
super(LVP,self).__init__()
self.th12 = 0.0
self.th13 = 0.0
self.th23 = 0.0
self.delta1 = 0.0
self.deltaCP = 0.0
super(LVP,self).Refresh()
# LVS
self.LVS21 = 0.0 #
self.LVS31 = 0.0 #
self.LVS41 = 0.0 #
self.LVS51 = 0.0 #
self.LVS61 = 0.0 #
# SQUARED MASS DIFFERENCE MATRIX
self.LVS = np.zeros([self.numneumax+2],float)
self.LVS[2] = self.LVS21
self.LVS[3] = self.LVS31
self.LVS[4] = self.LVS41
self.LVS[5] = self.LVS51
self.LVS[6] = self.LVS61
def Refresh(self):
super(LVP,self).Refresh()
LVS = self.LVS
LVS[2] = self.LVS21
LVS[3] = self.LVS31
LVS[4] = self.LVS41
LVS[5] = self.LVS51
LVS[6] = self.LVS61
def DiagonalMatrixLV(param):
DD = np.zeros([param.numneu,param.numneu],complex)
for k in np.arange(1,param.numneu,1):
DD[k,k] = param.LVS[k+1]
return DD
def LVTerm(LVparam):
U = calcU(LVparam)
return np.dot(U,np.dot(DiagonalMatrixLV(LVparam),U.conjugate().T))
|
