a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
ScalarPotential.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2016 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include "ScalarPotential.h"
9 
10 SUSYScalarPotential::SUSYScalarPotential(const StandardModel& SM_i) :
11  mySUSY(static_cast<const SUSY&> (SM_i))
12 
13 {}
14 
15 gslpp::vector<double> SUSYScalarPotential::coefficients()
16 {
17  gslpp::matrix<gslpp::complex> MsLhat2(3,3,0);
18  gslpp::matrix<gslpp::complex> MsEhat2(3,3,0);
19  MsLhat2 = mySUSY.getMsLhat2();
20  MsEhat2 = mySUSY.getMsEhat2();
21  double mstauLsq = MsLhat2(2,2).real();
22  double msmuRsq = MsEhat2(1,1).real();
23 // mhusq = ((mAsq/(1 + tb**2) + (self.mzsq/2)*((1 - tb**2)/(1 + tb**2)) - mu**2))
24 // mhdsq = ((mAsq * tb**2/(1 + tb**2) - (self.mzsq/2)*((1 - tb**2)/(1 + tb**2)) - mu**2))
25 // bmu = (mAsq) * tb/(1 + tb**2)
26 // double mA = mySUSY.getMHa();
27  double mA = 1000.;
28  double tanb = mySUSY.getTanb();
29  double MZ = mySUSY.getMz();
30  gslpp::complex muH = mySUSY.getMuH();
31  double musq = muH.abs2();
32  double mHdsq = mA*mA*tanb*tanb/(1.0+tanb*tanb)-0.5*MZ*MZ*(1.0-tanb*tanb)/(1.0+tanb*tanb)-musq;
33  double vd = mySUSY.v1();
34  double Mw = mySUSY.Mw_tree();
35  double sw2 = mySUSY.StandardModel::sW2(Mw);
36  double ttw2 = sw2/(1.0 - sw2);
37  double g2sq = 8. * mySUSY.getGF() / sqrt(2.0) * Mw * Mw;
38  double g1sq = g2sq*ttw2;
39  double mMU = mySUSY.getLeptons(StandardModel::MU).getMass();
40  double mTAU = mySUSY.getLeptons(StandardModel::TAU).getMass();
41  gslpp::matrix<gslpp::complex> TEhat = mySUSY.getTEhat();
42  double Al23 = TEhat(1,2).abs();
43 
44  std::cout << "tanb = " << tanb << std::endl;
45  std::cout << "MZ = " << MZ << std::endl;
46  std::cout << "mHdsq = " << mHdsq << std::endl;
47  std::cout << "g1 = " << sqrt(g1sq) << std::endl;
48  std::cout << "g2 = " << sqrt(g2sq) << std::endl;
49  std::cout << "Al23 = " << Al23 << std::endl;
50  std::cout << "vd = " << vd << std::endl;
51  std::cout << "ytau = " << mMU/vd*sqrt(2.) << std::endl;
52  std::cout << "ymu = " << mTAU/vd*sqrt(2.) << std::endl;
53 
54 
55  gslpp::vector<double> coefficients(35, 0.);
56 
57 // Let's call the scalar fields x1, x2 and x3.
58 // The coefficients will we the ones corresponding to:
59 // 1,
60 // x1^1, x2^1, x3^1,
61 // x1^2, x1*x2, x1*x3, x2^2, x2*x3, x3^2,
62 // x1^3, x1^2*x2, x1^2*x3, x1*x2^2, x1*x2*x3, x1*x3^2, x2^3, x2^2*x3, x2*x3^2, x3^3,
63 // x1^4, x1^3*x2, x1^3*x3, x1^2*x2^2, x1^2*x2*x3, x1^2*x3^2, x1*x2^3, x1*x2^2*x3, x1*x2*x3^2, x1*x3^3, x2^4, x2^3*x3, x2^2*x3^2, x2*x3^3, x3^4
64 
65 // Let's take x1=Hd, x2=stauL, x3=smuR
66  coefficients(0)=(mHdsq+musq+0.125*(g1sq+g2sq)*vd*vd)*vd*vd; //x1^0=x2^0=x3^0
67  coefficients(1)=-vd*(2.0*(mHdsq+musq)+0.5*vd*vd*(g1sq+g2sq)); //x1^1
68 // double x21 = 0.0;
69 // double x31 = 0.0;
70  coefficients(4)=mHdsq+musq+0.75*vd*vd*(g1sq+g2sq); //x1^2
71 // double x11x21 = 0.0;
72 // double x11x31 = 0.0;
73  double x22 = mstauLsq+mTAU*mTAU+0.25*vd*vd*(g1sq-g2sq);
74  coefficients(7)=x22;
75  double x21x31 = 2.0*Al23*vd;
76  coefficients(8)=x21x31;
77  double x32 = msmuRsq+mMU*mMU-0.5*vd*vd*g1sq;
78  coefficients(9)=x32;
79  double x13 = -0.5*vd*(g1sq+g2sq);
80  coefficients(10)=x13;
81 // double x12x21 = 0.0;
82 // double x12x31 = 0.0;
83  double x11x22 = -2.0*mTAU*mTAU/vd-0.5*vd*(g1sq-g2sq);
84  coefficients(13)=x11x22;
85  double x11x21x31 = -2.0*Al23;
86  coefficients(14)=x11x21x31;
87  double x11x32 = -2.0*mMU*mMU/vd+vd*g1sq;
88  coefficients(15)=x11x32;
89 // double x23 = 0.0;
90 // double x22x31 = 0.0;
91 // double x21x32 = 0.0;
92 // double x33 = 0.0;
93  double x14 = 0.125*(g1sq+g2sq);
94  coefficients(20)=x14;
95 // double x13x21 = 0.0;
96 // double x13x31 = 0.0;
97  double x12x22 = mTAU*mTAU/(vd*vd)+0.25*(g1sq-g2sq);
98  coefficients(23)=x12x22;
99 // double x12x21x31 = 0.0;
100  double x12x32 = mMU*mMU/(vd*vd)-0.5*g1sq;
101  coefficients(25)=x12x32;
102 // double x11x23 = 0.0;
103 // double x11x22x31 = 0.0;
104 // double x11x21x32 = 0.0;
105 // double x11x33 = 0.0;
106  double x24 = 0.125*(g1sq+g2sq);
107  coefficients(30)=x24;
108 // double x23x31 = 0.0;
109  double x22x32 = -0.5*g1sq;
110  coefficients(32)=x22x32;
111 // double x21x33 = 0.0;
112  double x34 = 0.5*g1sq;
113  coefficients(34)=x34;
114 
115  return coefficients;
116 }
117 
118 double SUSYScalarPotential::potential(gslpp::vector<double> coefficients, double field1, double field2, double field3)
119 {
120  gslpp::vector<double> a=coefficients;
121  double x1=field1;
122  double x2=field2;
123  double x3=field3;
124  double pot=a(0)
125  +a(1)*x1 +a(2)*x2 +a(3)*x3
126  +a(4)*x1*x1 +a(5)*x1*x2 +a(6)*x1*x3 +a(7)*x2*x2 +a(8)*x2*x3 +a(9)*x3*x3
127  +a(10)*x1*x1*x1 +a(11)*x1*x1*x2 +a(12)*x1*x1*x3 +a(13)*x1*x2*x2 +a(14)*x1*x2*x3 +a(15)*x1*x3*x3
128  +a(16)*x2*x2*x2 +a(17)*x2*x2*x3 +a(18)*x2*x3*x3 +a(19)*x3*x3*x3
129  +a(20)*x1*x1*x1*x1 +a(21)*x1*x1*x1*x2 +a(22)*x1*x1*x1*x3 +a(23)*x1*x1*x2*x2 +a(24)*x1*x1*x2*x3
130  +a(25)*x1*x1*x3*x3 +a(26)*x1*x2*x2*x2 +a(27)*x1*x2*x2*x3 +a(28)*x1*x2*x3*x3
131  +a(29)*x1*x3*x3*x3 +a(30)*x2*x2*x2*x2 +a(31)*x2*x2*x2*x3 +a(32)*x2*x2*x3*x3
132  +a(33)*x2*x3*x3*x3 +a(34)*x3*x3*x3*x3;
133 // double pot=0.0017924;
134  return pot;
135 }
136 
137 gslpp::vector<double> SUSYScalarPotential::potentialderivative(gslpp::vector<double> coefficients, double field1, double field2, double field3)
138 {
139  gslpp::vector<double> dV(3, 0.);
140  gslpp::vector<double> a=coefficients;
141  double x1=field1;
142  double x2=field2;
143  double x3=field3;
144  std::cout << "coefficients = " << coefficients << std::endl;
145  std::cout << "x1 = " << x1 << std::endl;
146  std::cout << "x2 = " << x2 << std::endl;
147  std::cout << "x3 = " << x3 << std::endl;
148 
149  //derivative wrt x1
150  dV(0)=a(1) + 2.0*a(4)*x1 + 3.0*a(10)*x1*x1 + 4.0*a(20)*x1*x1*x1
151  + a(5)*x2 + 2.0*a(11)*x1*x2 + 3.0*a(21)*x1*x1*x2 + a(13)*x2*x2 + 2.0*a(23)*x1*x2*x2 + a(26)*x2*x2*x2
152  + a(6)*x3 + 2.0*a(12)*x1*x3 + 3.0*a(22)*x1*x1*x3 + a(14)*x2*x3 + 2.0*a(24)*x1*x2*x3
153  + a(27)*x2*x2*x3 + a(15)*x3*x3 + 2.0*a(25)*x1*x3*x3 + a(28)*x2*x3*x3 + a(29)*x3*x3*x3;
154  //derivative wrt x2
155  dV(1)=a(2) + a(5)*x1 + a(11)*x1*x1 + a(21)*x1*x1*x1
156  + 2.0*a(7)*x2 + 2.0*a(13)*x1*x2 + 2.0*a(23)*x1*x1*x2 + 3.0*a(16)*x2*x2 + 3.0*a(26)*x1*x2*x2 + 4.0*a(30)*x2*x2*x2
157  + a(8)*x3 + a(14)*x1*x3 + a(24)*x1*x1*x3 + 2.0*a(17)*x2*x3 + 2.0*a(27)*x1*x2*x3
158  + 3.0*a(31)*x2*x2*x3 + a(18)*x3*x3 + a(28)*x1*x3*x3 + 2.0*a(32)*x2*x3*x3 + a(33)*x3*x3*x3;
159  //derivative wrt x3
160  dV(2)=a(3) + a(6)*x1 + a(12)*x1*x1 + a(22)*x1*x1*x1
161  + a(8)*x2 + a(14)*x1*x2 + a(24)*x1*x1*x2 + a(17)*x2*x2 + a(27)*x1*x2*x2 + a(31)*x2*x2*x2
162  + 2.0*a(9)*x3 + 2.0*a(15)*x1*x3 + 2.0*a(25)*x1*x1*x3 + 2.0*a(18)*x2*x3 + 2.0*a(28)*x1*x2*x3
163  + 2.0*a(32)*x2*x2*x3 + 3.0*a(19)*x3*x3 + 3.0*a(29)*x1*x3*x3 + 3.0*a(33)*x2*x3*x3 + 4.0*a(34)*x3*x3*x3;
164  return dV;
165 }
166 
167 gslpp::vector<double> SUSYScalarPotential::secondpotentialderivative(gslpp::vector<double> coefficients, double field1, double field2, double field3)
168 {
169  gslpp::vector<double> d2V(6, 0.);
170  gslpp::vector<double> a=coefficients;
171  double x1=field1;
172  double x2=field2;
173  double x3=field3;
174 
175  //derivative wrt x1*x1
176  d2V(0)=2.0*a(4) + 6.0*a(10)*x1 + 12.0*a(20)*x1*x1 + 2.0*a(11)*x2 + 6.0*a(21)*x1*x2
177  + 2.0*a(23)*x2*x2 + 2.0*a(12)*x3 + 6.0*a(22)*x1*x3 + 2.0*a(24)*x2*x3 + 2.0*a(25)*x3*x3;
178  //derivative wrt x1*x2
179  d2V(1)=a(5) + 2.0*a(11)*x1 + 3.0*a(21)*x1*x1 + 2.0*a(13)*x2 + 4.0*a(23)*x1*x2
180  + 3.0*a(26)*x2*x2 + a(14)*x3 + 2.0*a(24)*x1*x3 + 2.0*a(27)*x2*x3 + a(28)*x3*x3;
181  //derivative wrt x1*x3
182  d2V(2)=a(6) + 2.0*a(12)*x1 + 3.0*a(22)*x1*x1 + a(14)*x2 + 2.0*a(24)*x1*x2 + a(27)*x2*x2
183  + 2.0*a(15)*x3 + 4.0*a(25)*x1*x3 + 2.0*a(28)*x2*x3 + 3.0*a(29)*x3*x3;
184  //derivative wrt x2*x2
185  d2V(3)=2.0*a(7) + 2.0*a(13)*x1 + 2.0*a(23)*x1*x1 + 6.0*a(16)*x2 + 6.0*a(26)*x1*x2
186  + 12.0*a(30)*x2*x2 + 2.0*a(17)*x3 + 2.0*a(27)*x1*x3 + 6.0*a(31)*x2*x3 + 2.0*a(32)*x3*x3;
187  //derivative wrt x2*x3
188  d2V(4)=a(8) + a(14)*x1 + a(24)*x1*x1 + 2.0*a(17)*x2 + 2.0*a(27)*x1*x2 + 3.0*a(31)*x2*x2
189  + 2.0*a(18)*x3 + 2.0*a(28)*x1*x3 + 4.0*a(32)*x2*x3 + 3.0*a(33)*x3*x3;
190  //derivative wrt x3*x3
191  d2V(5)=2.0*a(9) + 2.0*a(15)*x1 + 2.0*a(25)*x1*x1 + 2.0*a(18)*x2 + 2.0*a(28)*x1*x2
192  + 2.0*a(32)*x2*x2 + 6.0*a(19)*x3 + 6.0*a(29)*x1*x3 + 6.0*a(33)*x2*x3 + 12.0*a(34)*x3*x3;
193  return d2V;
194 }
QCD::TAU
Definition: QCD.h:316
SUSYScalarPotential::secondpotentialderivative
gslpp::vector< double > secondpotentialderivative(gslpp::vector< double > coefficients, double field1, double field2, double field3)
Definition: ScalarPotential.cpp:167
ScalarPotential.h
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
gslpp::matrix< gslpp::complex >
gslpp::complex::abs2
double abs2() const
Definition: gslpp_complex.cpp:86
SUSYScalarPotential::SUSYScalarPotential
SUSYScalarPotential(const StandardModel &SM_i)
SUSYScalarPotential constructor.
Definition: ScalarPotential.cpp:10
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
SUSYScalarPotential::potentialderivative
gslpp::vector< double > potentialderivative(gslpp::vector< double > coefficients, double field1, double field2, double field3)
Definition: ScalarPotential.cpp:137
SUSYScalarPotential::coefficients
gslpp::vector< double > coefficients()
Definition: ScalarPotential.cpp:15
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
SUSYScalarPotential::potential
double potential(gslpp::vector< double > coefficients, double field1, double field2, double field3)
Definition: ScalarPotential.cpp:118
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
SUSYScalarPotential::mySUSY
const SUSY & mySUSY
Definition: ScalarPotential.h:48
Mw
An observable class for the -boson mass.
Definition: Mw.h:22
QCD::MU
Definition: QCD.h:314