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