flavour/eigen.cpp at master · leonardopedro/flavour · GitHub

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
#include <Eigen/Eigenvalues>
#include <iostream>
#include <vector>
using namespace std;
using namespace Eigen;

class points{
public:
	points(double p, double m, double e):prediction(p),measure(m),error(e){}
double prediction, measure, error;
};

double chi2(double * v){
	double buu = v[0], auu = v[1];
	double add = v[2], bdd = v[3];
	double eu =v[4], gu = v[5];
	double ed = v[6], gd = v[7];
complex<double> bu = buu*exp(complex<double>(0,gu*M_PI_2)), au = auu*exp(complex<double>(0,eu*M_PI_2));
complex<double> ad = add*exp(complex<double>(0,ed*M_PI_2)), bd = bdd*exp(complex<double>(0,gd*M_PI_2));
complex<double> cd = bd - ad;
complex<double> cu =bu + au;

//Matrix3cd X = Matrix3cd::Random(3,3);
Matrix3cd mu,md;
mu<<1,1,1.0+au+bu,1,1,1.0+bu,1.0+bu+au,1.0+bu,1.0+cu;
md<<1,1,1.0+ad+bd,1,1,1.0+bd,1.0+bd+ad,1.0+bd,1.0+cd;
	mu=mu*v[8];
	md=md*v[9];

Matrix3cd Hu = mu*mu.adjoint(), Hd = md*md.adjoint();
SelfAdjointEigenSolver<Matrix3cd> usolver(Hu), dsolver(Hd);
Vector3d Du=usolver.eigenvalues();//.asDiagonal();
Vector3d Dd=dsolver.eigenvalues();//.asDiagonal();
Matrix3cd VLd=dsolver.eigenvectors().adjoint();
Matrix3cd VLu=usolver.eigenvectors().adjoint();
Matrix3cd Vckm=VLu*VLd.adjoint();
double lambda=sqrt(norm(Vckm(0,1))/(norm(Vckm(0,0))+norm(Vckm(0,1))));
double A=abs(Vckm(1,2))/abs(Vckm(0,1))/lambda;
complex<double> rhoeta=-Vckm(0,0)*conj(Vckm(0,2))/Vckm(1,0)/conj(Vckm(1,2));
double rho=real(rhoeta), eta=imag(rhoeta);

vector<points> comp;
comp.push_back(points(lambda, 0.22535, 0.00065));
comp.push_back(points(A, 0.811, 0.017));
comp.push_back(points(rho, 0.131, 0.02));
comp.push_back(points(eta, 0.345, 0.014));
comp.push_back(points(eta, 0.345, 0.014));
comp.push_back(points(Du[0], 1.27e-3, 0.46e-3));
comp.push_back(points(Dd[0], 2.9e-3, 1.22e-3));
comp.push_back(points(Dd[1], 55e-3, 16e-3));
comp.push_back(points(Du[1], 0.619, 0.084));
comp.push_back(points(Dd[2], 2.89, 0.09));
comp.push_back(points(Du[2], 171.7, 3.0));

double chi=0;

for(uint i=0;i<comp.size();i++){
	chi+=pow((comp[i].prediction-comp[i].measure)/comp[i].error,2);
	}
return chi/comp.size();
}

int main(){
double buu = 0.0162, auu = 0.0009;
double add = 0.018, bdd = 0.09;
//double x = 1.0/3;
double eu = -1.8, gu = -1.0/2;
double ed = -1.0/2, gd = -2;
complex<double> bu = buu*exp(complex<double>(0,gu*M_PI_2)), au = auu*exp(complex<double>(0,eu*M_PI_2));
complex<double> ad = add*exp(complex<double>(0,ed*M_PI_2)), bd = bdd*exp(complex<double>(0,gd*M_PI_2));
complex<double> cd = exp(complex<double>(0,0.09/1.8));

//Matrix3cd X = Matrix3cd::Random(3,3);
Matrix3cd mu,md;
mu<<1,1,1.0+au+bu,1,1,1.0+bu,1.0+bu+au,1.0+bu,1.0+bu;
//md<<1,1,1.0+ad+bd,1,1,1.0+bd,1.0+bd+ad,1.0+bd,1.0+bd;
md<<1,1,1.0+ad+bd,1,1,1.0+bd,cd*(1.0+bd+ad),cd*(1.0+bd),cd*(1.0+bd);
mu=mu*57.4822;
md=md*1.0147;

Matrix3cd Hu = mu*mu.adjoint(), Hd = md*md.adjoint();
SelfAdjointEigenSolver<Matrix3cd> usolver(Hu), dsolver(Hd);
Vector3d Du=usolver.eigenvalues();//.asDiagonal();
Vector3d Dd=dsolver.eigenvalues();//.asDiagonal();
Matrix3cd VLd=dsolver.eigenvectors().adjoint();
Matrix3cd VLu=usolver.eigenvectors().adjoint();
Matrix3cd Vckm=VLu*VLd.adjoint();
double lambda=sqrt(norm(Vckm(0,1))/(norm(Vckm(0,0))+norm(Vckm(0,1))));
double A=abs(Vckm(1,2))/abs(Vckm(0,1))/lambda;
complex<double> rhoeta=-Vckm(0,0)*conj(Vckm(0,2))/Vckm(1,0)/conj(Vckm(1,2));
double rho=real(rhoeta), eta=imag(rhoeta);
double fu=173.5/sqrt(Du[2]), fd=2.89/sqrt(Dd[2]);

vector<points> comp;
comp.push_back(points(lambda, 0.22535, 0.00065));
comp.push_back(points(A, 0.811, 0.017));
comp.push_back(points(rho, 0.131, 0.02));
comp.push_back(points(eta, 0.345, 0.014));
comp.push_back(points(sqrt(Du[0]), 1.27e-3, 0.46e-3));
comp.push_back(points(sqrt(Dd[0]), 2.9e-3, 1.22e-3));
comp.push_back(points(sqrt(Dd[1]), 55e-3, 16e-3));
comp.push_back(points(sqrt(Du[1]), 0.619, 0.084));
comp.push_back(points(sqrt(Dd[2]), 2.89, 0.09));
comp.push_back(points(sqrt(Du[2]), 171.7, 3.0));

double chi=0;

for(uint i=0;i<comp.size();i++){
	double aa=(comp[i].prediction-comp[i].measure)/comp[i].error;
	cout<<i<<" "<<aa<<endl;
	chi+=pow((comp[i].prediction-comp[i].measure)/comp[i].error,2);
	}
chi/=comp.size();

cout<<lambda<<" "<<A<<" "<<rho<<" "<<eta<<endl;
cout<<sqrt(Du[2])<<" "<<sqrt(Du[1])<<" "<<sqrt(Du[0])<<endl;
cout<<sqrt(Dd[2])<<" "<<sqrt(Dd[1])<<" "<<sqrt(Dd[0])<<endl;
cout<<chi<<endl;
/*
cout << "The eigenvalues of mu are:" << endl << es.eigenvalues() << endl;
cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
double lambda = es.eigenvalues()[0];
cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
Vector3cd v = es.eigenvectors().col(0);
cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
cout << "... and A * v = " << endl << A * v << endl << endl;
Matrix3d D = es.eigenvalues().asDiagonal();
Matrix3cd V = es.eigenvectors();
cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
*/
return 0;
}