a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
GeneralTHDMRunner.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 "GeneralTHDMRunner.h"
9 #include <gsl/gsl_errno.h>
10 #include <gsl/gsl_matrix.h>
11 #include <gsl/gsl_odeiv2.h>
12 
13 GeneralTHDMRunner::GeneralTHDMRunner(const StandardModel& SM_i) : myGTHDM(static_cast<const GeneralTHDM*> (&SM_i))
14 {}
15 
16 GeneralTHDMRunner::~GeneralTHDMRunner()
17 {};
18 
19 int RGEsGTHDM(double t, const double y[], double beta[], void *flags)
20 {
21  (void)(t); /* avoid unused parameter warning */
22  int ord = *(int *)flags;
23  double g1=y[0];
24  double g2=y[1];
25  double g3=y[2];
26  double v1=y[3];
27  double v2=y[4];
28  double etaU1=y[5];
29  double etaU2=y[6];
30  double etaD1=y[7];
31  double etaD2=y[8];
32  double etaL1=y[9];
33  double etaL2=y[10];
34  double m11sq=y[11];
35  double m12sq=y[12];
36  double m22sq=y[13];
37  double la1=y[14];
38  double la2=y[15];
39  double la3=y[16];
40  double la4=y[17];
41  double la5=y[18];
42  double la6=y[19];
43  double la7=y[20];
44 
45 // double pi=M_PI;
46 
47  //beta_g1
48  beta[0] = 7.*g1*g1*g1;
49  //beta_g2
50  beta[1] = -3.*g2*g2*g2;
51  //beta_g3
52  beta[2] = -7.*g3*g3*g3;
53  //beta_v1
54  beta[3] = (-3.*etaD1*etaD1-etaL1*etaL1-3.*etaU1*etaU1+0.75*g1*g1+2.25*g2*g2)*v1
55  +(-3.*etaD1*etaD2-etaL1*etaL2-3.*etaU1*etaU2)*v2;
56  //beta_v2
57  beta[4] = (-3.*etaD1*etaD2-etaL1*etaL2-3.*etaU1*etaU2)*v1
58  +(-3.*etaD2*etaD2-etaL2*etaL2-3.*etaU2*etaU2+0.75*g1*g1+2.25*g2*g2)*v2;
59  //beta_etaU1
60  beta[5] = 1.5*etaD1*etaD1*etaU1+0.5*etaD2*etaD2*etaU1+etaL1*etaL1*etaU1+4.5*etaU1*etaU1*etaU1
61  +etaD1*etaD2*etaU2+etaL1*etaL2*etaU2+4.5*etaU1*etaU2*etaU2
62  -1.41667*etaU1*g1*g1-2.25*etaU1*g2*g2-8.*etaU1*g3*g3;
63  //beta_etaU2
64  beta[6] = etaD1*etaD2*etaU1+etaL1*etaL2*etaU1+0.5*etaD1*etaD1*etaU2
65  +etaU2*(1.5*etaD2*etaD2+etaL2*etaL2+4.5*etaU1*etaU1+4.5*etaU2*etaU2
66  -1.41667*g1*g1-2.25*g2*g2-8.*g3*g3);
67  //beta_etaD1
68  beta[7] = 4.5*etaD1*etaD1*etaD1+etaD2*(etaL1*etaL2+etaU1*etaU2)
69  +etaD1*(4.5*etaD2*etaD2+etaL1*etaL1+1.5*etaU1*etaU1+0.5*etaU2*etaU2
70  -0.416667*g1*g1-2.25*g2*g2-8.*g3*g3);
71  //beta_etaD2
72  beta[8] = 4.5*etaD1*etaD1*etaD2+etaD1*(etaL1*etaL2+etaU1*etaU2)
73  +etaD2*(4.5*etaD2*etaD2+etaL2*etaL2+0.5*etaU1*etaU1+1.5*etaU2*etaU2
74  -0.416667*g1*g1-2.25*g2*g2-8.*g3*g3);
75  //beta_etaL1
76  beta[9] = 3.*etaD1*etaD1*etaL1+2.5*etaL1*etaL1*etaL1+3.*etaD1*etaD2*etaL2
77  +3.*etaL2*etaU1*etaU2+etaL1*(2.5*etaL2*etaL2+3.*etaU1*etaU1-3.75*g1*g1-2.25*g2*g2);
78  //beta_etaL2
79  beta[10] = 3.*etaD1*etaD2*etaL1+3.*etaD2*etaD2*etaL2+2.5*etaL1*etaL1*etaL2
80  +2.5*etaL2*etaL2*etaL2+3.*etaL1*etaU1*etaU2+3.*etaL2*etaU2*etaU2
81  -3.75*etaL2*g1*g1-2.25*etaL2*g2*g2;
82  //beta_m11sq
83  beta[11] = 6.*etaD1*etaD1*m11sq+2.*etaL1*etaL1*m11sq+6.*etaU1*etaU1*m11sq
84  -1.5*g1*g1*m11sq-4.5*g2*g2*m11sq+6.*la1*m11sq
85  -6.*etaD1*etaD2*m12sq-2.*etaL1*etaL2*m12sq-6.*etaU1*etaU2*m12sq
86  -12.*la6*m12sq+4.*la3*m22sq+2.*la4*m22sq;
87  //beta_m12sq
88  beta[12] = (3.*etaD1*etaD1+3.*etaD2*etaD2+etaL1*etaL1+etaL2*etaL2+3.*etaU1*etaU1
89  +3.*etaU2*etaU2+2.*la3+4.*la4)*m12sq
90  -3.*etaD1*etaD2*(m11sq+m22sq)-etaL1*etaL2*(m11sq+m22sq)
91  -1.5*(4.*la6*m11sq+g1*g1*m12sq+3.*g2*g2*m12sq-4.*la5*m12sq
92  +4.*la7*m22sq+2.*etaU1*etaU2*(m11sq+m22sq));
93  //beta_m22sq
94  beta[13] = 4.*la3*m11sq+2.*la4*m11sq-6.*etaD1*etaD2*m12sq-2.*etaL1*etaL2*m12sq
95  -6.*etaU1*etaU2*m12sq-12.*la7*m12sq+6.*etaD2*etaD2*m22sq+2.*etaL2*etaL2*m22sq
96  +6.*etaU2*etaU2*m22sq-1.5*g1*g1*m22sq-4.5*g2*g2*m22sq+6.*la2*m22sq;
97  //beta_la1
98  beta[14] = -12.*etaD1*etaD1*etaD1*etaD1-4.*etaL1*etaL1*etaL1*etaL1-12.*etaU1*etaU1*etaU1*etaU1
99  +0.75*g1*g1*g1*g1+1.5*g1*g1*g2*g2+2.25*g2*g2*g2*g2
100  +12.*etaD1*etaD1*la1+4.*etaL1*etaL1*la1+12.*etaU1*etaU1*la1
101  -3.*g1*g1*la1-9.*g2*g2*la1
102  +12.*la1*la1+4.*la3*la3+4.*la3*la4+2.*la4*la4
103  +12.*etaD1*etaD2*la6+4.*etaL1*etaL2*la6+12.*etaU1*etaU2*la6
104  +24.*la6*la6+2.*la5*la5;
105  //beta_la2
106  beta[15] = -12.*etaD2*etaD2*etaD2*etaD2-4.*etaL2*etaL2*etaL2*etaL2-12.*etaU2*etaU2*etaU2*etaU2
107  +0.75*g1*g1*g1*g1+1.5*g1*g1*g2*g2+2.25*g2*g2*g2*g2
108  +12.*etaD2*etaD2*la2+4.*etaL2*etaL2*la2+12.*etaU2*etaU2*la2
109  -3.*g1*g1*la2-9.*g2*g2*la2
110  +12.*la2*la2+4.*la3*la3+4.*la3*la4+2.*la4*la4
111  +12.*etaD1*etaD2*la7+4.*etaL1*etaL2*la7+12.*etaU1*etaU2*la7
112  +24.*la7*la7+2.*la5*la5;
113  //beta_la3
114  beta[16] = -12.*etaD1*etaD1*etaD2*etaD2-4.*etaL1*etaL1*etaL2*etaL2-12.*etaD2*etaD2*etaU1*etaU1
115  +24.*etaD1*etaD2*etaU1*etaU2-12.*etaD1*etaD1*etaU2*etaU2-12.*etaU1*etaU1*etaU2*etaU2
116  +0.75*(g1*g1*g1*g1+3.*g2*g2*g2*g2+g1*g1*(-2.*g2*g2-4.*la3)-12.*g2*g2*la3)
117  +2.*(la3*(3.*etaD1*etaD1+3.*etaD2*etaD2+etaL1*etaL1+etaL2*etaL2
118  +3.*(etaU1*etaU1+etaU2*etaU2+la1+la2)+2.*la3)
119  +(la1+la2)*la4+la4*la4)
120  +(3.*etaD1*etaD2+etaL1*etaL2+3.*etaU1*etaU2)*(la6+la7)
121  +la7*(3.*etaD1*etaD2+etaL1*etaL2+3.*etaU1*etaU2+8.*la6+4.*la7)
122  +la6*(3.*etaD1*etaD2+etaL1*etaL2+3.*etaU1*etaU2+4.*la6+8.*la7)+2.*la5*la5;
123  //beta_la4
124  beta[17] = 12.*etaD2*etaD2*etaU1*etaU1-12.*etaU1*etaU1*etaU2*etaU2
125  +3.*g1*g1*g2*g2+6.*etaD2*etaD2*la4+2.*etaL2*etaL2*la4+6.*etaU1*etaU1*la4
126  +6.*etaU2*etaU2*la4-3.*g1*g1*la4-9.*g2*g2*la4
127  +2.*la1*la4+2.*la2*la4+8.*la3*la4+4.*la4*la4
128  +etaL1*etaL1*(-4.*etaL2*etaL2+2.*la4)
129  +etaD1*etaD1*(-12.*etaD2*etaD2+12.*etaU2*etaU2+6.*la4)
130  +8.*la5*la5+6.*etaU1*etaU2*la6+10.*la6*la6+6.*etaU1*etaU2*la7
131  +4.*la6*la7+10.*la7*la7+etaL1*etaL2*(2.*la6+2.*la7)
132  +etaD1*etaD2*(-24.*etaU1*etaU2+6.*la6+6.*la7);
133  //beta_la5
134  beta[18] = -12.*etaU1*etaU1*etaU2*etaU2+6.*etaD2*etaD2*la5+2.*etaL2*etaL2*la5
135  +6.*etaU1*etaU1*la5+6.*etaU2*etaU2*la5-3.*g1*g1*la5-9.*g2*g2*la5
136  +2.*la1*la5+2.*la2*la5+8.*la3*la5+12.*la4*la5
137  +etaL1*etaL1*(-4.*etaL2*etaL2+2.*la5)
138  +etaD1*etaD1*(-12.*etaD2*etaD2+6.*la5)
139  +6.*etaU1*etaU2*la6+10.*la6*la6+6.*etaU1*etaU2*la7
140  +4.*la6*la7+10.*la7*la7+etaL1*etaL2*(2.*la6+2.*la7)
141  +etaD1*etaD2*(6.*la6+6.*la7);
142  //beta_la6
143  beta[19] = -12.*etaD1*etaD1*etaD1*etaD2-4.*etaL1*etaL1*etaL1*etaL2-12.*etaU1*etaU1*etaU1*etaU2
144  +3.*etaU1*etaU2*la1+3.*etaU1*etaU2*la3+3.*etaU1*etaU2*la4
145  +3.*etaU1*etaU2*la5+etaL1*etaL2*(la1+la3+la4+la5)
146  +3.*etaD1*etaD2*(la1+la3+la4+la5)+9.*etaD1*etaD1*la6
147  +3.*etaD2*etaD2*la6+3.*etaL1*etaL1*la6+etaL2*etaL2*la6
148  +9.*etaU1*etaU1*la6+3.*etaU2*etaU2*la6-3.*g1*g1*la6-9.*g2*g2*la6
149  +12.*la1*la6+6.*la3*la6+8.*la4*la6+10.*la5*la6
150  +6.*la3*la7+4.*la4*la7+2.*la5*la7;
151  //beta_la7
152  beta[20] = -12.*etaU1*etaU2*etaU2*etaU2+3.*etaU1*etaU2*la2+3.*etaU1*etaU2*la3
153  +3.*etaU1*etaU2*la4+3.*etaU1*etaU2*la5
154  +etaL1*etaL2*(-4.*etaL2*etaL2+la2+la3+la4+la5)
155  +etaD1*etaD2*(-12.*etaD2*etaD2+3.*la2+3.*la3+3.*la4+3.*la5)
156  +6.*la3*la6+4.*la4*la6+2.*la5*la6+3.*etaD1*etaD1*la7
157  +9.*etaD2*etaD2*la7+etaL1*etaL1*la7+3.*etaL2*etaL2*la7
158  +3.*etaU1*etaU1*la7+9.*etaU2*etaU2*la7
159  -3.*g1*g1*la7-9.*g2*g2*la7+12.*la2*la7
160  +6.*la3*la7+8.*la4*la7+10.*la5*la7;
161 
162  if(ord==1){
163  //beta_g1
164  beta[0] += 0.;
165  //beta_g2
166  beta[1] += 0.;
167  //beta_g3
168  beta[2] += 0.;
169  //beta_v1
170  beta[3] += 0.;
171  //beta_v2
172  beta[4] += 0.;
173  //beta_etaU1
174  beta[5] += 0.;
175  //beta_etaU2
176  beta[6] += 0.;
177  //beta_etaD1
178  beta[7] += 0.;
179  //beta_etaD2
180  beta[8] += 0.;
181  //beta_etaL1
182  beta[9] += 0.;
183  //beta_etaL2
184  beta[10] += 0.;
185  //beta_m11sq
186  beta[11] += 0.;
187  //beta_m12sq
188  beta[12] += 0.;
189  //beta_m22sq
190  beta[13] += 0.;
191  //beta_la1
192  beta[14] += 0.;
193  //beta_la2
194  beta[15] += 0.;
195  //beta_la3
196  beta[16] += 0.;
197  //beta_la4
198  beta[17] += 0.;
199  //beta_la5
200  beta[18] += 0.;
201  //beta_la6
202  beta[19] += 0.;
203  //beta_la7
204  beta[20] += 0.;
205  }
206 
207  return 0;
208 }
209 
210 int JacobianGTHDM (double t, const double y[], double *dfdy, double dfdt[], void *order)
211 {
212  return 0;
213 }
214 
215 int RGEcheckGTHDM(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
216 {
217  int check=0;
218 
219  //perturbativity of the Yukawa couplings
220  for(int i=5;i<11;i++)
221  {
222  if(fabs(InitialValues[i])>sqrt(4.0*M_PI)) check=1;
223  }
224 
225  //perturbativity of the quartic Higgs couplings
226  for(int i=14;i<21;i++)
227  {
228  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
229  }
230 
231  double la1Q = InitialValues[14];
232  double la2Q = InitialValues[15];
233 // double la3Q = InitialValues[16];
234 // double la4Q = InitialValues[17];
235 // double la5Q = InitialValues[18];
236 // double la6Q = InitialValues[19];
237 // double la7Q = InitialValues[20];
238 
239  //positivity checks
240  if(la1Q<0.0) check=1;
241  if(la2Q<0.0) check=1;
242 // if(la3Q+sqrt(la1Q*la2Q)<0.0) check=1;
243 // if(la3Q+2.0*la4Q+sqrt(la1Q*la2Q)<0.0) check=1;
244 
246  //NLO unitarity//
248 
249 // double pi=M_PI;
250 // gslpp::matrix<gslpp::complex> Smatrix1(4,4,0.), Smatrix2(4,4,0.);
251 // gslpp::matrix<gslpp::complex> Seigenvectors1(4,4,0.), Seigenvectors2(4,4,0.);
252 // gslpp::vector<double> Seigenvalues1(4,0.), Seigenvalues2(4,0.);
253 // gslpp::vector<gslpp::complex> unitarityeigenvalues(11,0.);
254 
255  if(t1>tNLOuni)
256  {
257  } //end of the if(t1>tNLOuni)
258  return check;
259 }
260 
261 double GeneralTHDMRunner::RGEGeneralTHDMRunner(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double NLOuniscale)
262 {
263  //Define which stepping function should be used
264  const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk4;
265 
266  //Allocate space for the stepping function
267  gsl_odeiv2_step * s = gsl_odeiv2_step_alloc(T, NumberOfRGEs);
268 
269  //Define the absolute (A) and relative (R) error on y at each step.
270  //The real error will be compared to the following error estimate:
271  // A + R * |y_i|
272  gsl_odeiv2_control * c = gsl_odeiv2_control_y_new(1e-6, 0.0);
273 
274  //Allocate space for the evolutor
275  gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(NumberOfRGEs);
276 
277  //Definition of the RGE system (the Jacobian is not necessary for the RK4 method; it's an empty function here)
278  gsl_odeiv2_system RGEsystem = {RGEsGTHDM, JacobianGTHDM, NumberOfRGEs, &order};
279 
280  //Set starting and end point as natural logarithmic scales (conversion from decadic log scale)
281  double t1 = Q1*log(10.0);
282  double t2 = Q2*log(10.0);
283  double tNLOuni = NLOuniscale*log(10.0);
284 
285  //Set initial step size
286  double InitialStepSize = 1e-6;
287 
288  //Run!
289  while (t1 < t2)
290  {
291  int status = gsl_odeiv2_evolve_apply (e, c, s, &RGEsystem, &t1, t2, &InitialStepSize, InitialValues);
292  if(status != GSL_SUCCESS) break;
293 
294  //intermediate checks if appropriate
295  if(RGEcheckGTHDM(InitialValues,t1,Rpeps,tNLOuni) != 0) break;
296  }
297 
298  gsl_odeiv2_evolve_free (e);
299  gsl_odeiv2_control_free (c);
300  gsl_odeiv2_step_free (s);
301 
302  //Return the decadic log scale at which the evolution stopped
303  return t1/log(10.0);
304 }
GeneralTHDMRunner::GeneralTHDMRunner
GeneralTHDMRunner(const StandardModel &SM_i)
GeneralTHDMRunner constructor.
Definition: GeneralTHDMRunner.cpp:13
GeneralTHDMRunner::RGEGeneralTHDMRunner
virtual double RGEGeneralTHDMRunner(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double tNLOuni)
Definition: GeneralTHDMRunner.cpp:261
GeneralTHDMRunner.h
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
GeneralTHDMRunner::~GeneralTHDMRunner
virtual ~GeneralTHDMRunner()
Runner destructor.
Definition: GeneralTHDMRunner.cpp:16
RGEsGTHDM
int RGEsGTHDM(double t, const double y[], double beta[], void *flags)
Definition: GeneralTHDMRunner.cpp:19
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
GeneralTHDM
A class for general Two-Higgs-Doublet models.
Definition: GeneralTHDM.h:463
RGEcheckGTHDM
int RGEcheckGTHDM(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: GeneralTHDMRunner.cpp:215
JacobianGTHDM
int JacobianGTHDM(double t, const double y[], double *dfdy, double dfdt[], void *order)
Definition: GeneralTHDMRunner.cpp:210