a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Runner.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 "Runner.h"
9#include <gsl/gsl_errno.h>
10#include <gsl/gsl_matrix.h>
11#include <gsl/gsl_odeiv2.h>
12#include "THDM.h"
13#include "THDMquantities.h"
14
15Runner::Runner(const StandardModel& SM_i) : ThObservable(SM_i), myTHDM(static_cast<const THDM*> (&SM_i))
16{
17 mym11_2=new m11_2(SM_i);
18 mym22_2=new m22_2(SM_i);
19 mylambda1=new lambda1(SM_i);
20 mylambda2=new lambda2(SM_i);
21 mylambda3=new lambda3(SM_i);
22 mylambda4=new lambda4(SM_i);
23 mylambda5=new lambda5(SM_i);
24}
25
26Runner::~Runner()
27{
28 delete mym11_2;
29 delete mym22_2;
30 delete mylambda1;
31 delete mylambda2;
32 delete mylambda3;
33 delete mylambda4;
34 delete mylambda5;
35};
36
37double Runner::computeThValue()
38{
39 return 0.;
40}
41
42int RGEs(double t, const double y[], double beta[], void *flags)
43{
44 (void)(t); /* avoid unused parameter warning */
45 int flag = *(int *)flags;
46 int ord = flag%3;
47 int type = flag-ord;
48 double Yb1=0;
49 double Yb2=0;
50 double Ytau1=0;
51 double Ytau2=0;
52 double g1=y[0];
53 double g2=y[1];
54 double g3=y[2];
55 double Yt=y[3];
56 if( type == 0 ) {//type I
57 Yb2=y[4];
58 Ytau2=y[5];
59 }
60 else if( type == 3 ) {//type II
61 Yb1=y[4];
62 Ytau1=y[5];
63 }
64 else if( type == 6 ) {//type X
65 Yb2=y[4];
66 Ytau1=y[5];
67 }
68 else if( type == 9 ) {//type Y
69 Yb1=y[4];
70 Ytau2=y[5];
71 }
72 else {
73 throw std::runtime_error("type should only be any of 0, 3, 6, 9");
74 }
75 double m11_2=y[6];
76 double m22_2=y[7];
77 double m12_2=y[8];
78 double la1=y[9];
79 double la2=y[10];
80 double la3=y[11];
81 double la4=y[12];
82 double la5=y[13];
83
84 double pi=M_PI;
85
86 //beta_g1
87 beta[0] = 7.0*g1*g1*g1/(16.0*pi*pi);
88 //beta_g2
89 beta[1] = -3.0*g2*g2*g2/(16.0*pi*pi);
90 //beta_g3
91 beta[2] = -7.0*g3*g3*g3/(16.0*pi*pi);
92 //beta_Yt
93 beta[3] = Yt*(-17.0*g1*g1+3.0*(-9.0*g2*g2-32.0*g3*g3+2.0*Yb1*Yb1+6.0*Yb2*Yb2+18.0*Yt*Yt+4.0*Ytau2*Ytau2))/(192.0*pi*pi);
94 //beta_Yb
95 beta[4] = ((-5.0*g1*g1-27.0*g2*g2-96.0*g3*g3)*(Yb1+Yb2)+(6.0*Yb1+18.0*Yb2)*Yt*Yt
96 +54.0*(Yb1*Yb1*Yb1+Yb2*Yb2*Yb2+Yb1*Yb2*Yb2)+12.0*(Yb1*Ytau1*Ytau1+Yb2*Ytau2*Ytau2))/(192.0*pi*pi);
97 //beta_Ytau
98 beta[5] = ((-15.0*g1*g1-9.0*g2*g2)*(Ytau1+Ytau2)+Ytau2*(12.0*Yb2*Yb2+12.0*Yt*Yt+10.0*Ytau2*Ytau2)
99 +Ytau1*(12.0*Yb1*Yb1+10.0*Ytau1*Ytau1))/(64.0*pi*pi);
100 //beta_m11_2
101 beta[6] = (-3.0*g1*g1*m11_2-9.0*g2*g2*m11_2
102 +4.0*(m11_2*(3.0*Yb1*Yb1+Ytau1*Ytau1+3.0*la1)+m22_2*(2.0*la3+la4)))/(32.0*pi*pi);
103 //beta_m22_2
104 beta[7] = (-3.0*g1*g1*m22_2-9.0*g2*g2*m22_2
105 +4.0*(m22_2*(3.0*Yb2*Yb2+3.0*Yt*Yt+Ytau2*Ytau2+3.0*la2)+m11_2*(2.0*la3+la4)))/(32.0*pi*pi);
106 //beta_m12_2
107 beta[8] = m12_2*(-3.0*g1*g1-9.0*g2*g2
108 +2.0*(3.0*Yb1*Yb1+3.0*Yb2*Yb2+3.0*Yt*Yt+Ytau1*Ytau1+Ytau2*Ytau2+2.0*la3+4.0*la4+6.0*la5))/(32.0*pi*pi);
109 //beta_lambda_1
110 beta[9] = (3.0*g1*g1*g1*g1+9.0*g2*g2*g2*g2+6.0*g1*g1*(g2*g2-2.0*la1)-36.0*g2*g2*la1
111 +8.0*(6.0*la1*la1+2.0*la3*la3+2.0*la3*la4+la4*la4+la5*la5)
112 -16.0*(3.0*Yb1*Yb1*Yb1*Yb1+Ytau1*Ytau1*Ytau1*Ytau1)+16.0*la1*(3.0*Yb1*Yb1+Ytau1*Ytau1))/(64.0*pi*pi);
113 //beta_lambda_2
114 beta[10] = (3.0*g1*g1*g1*g1+9.0*g2*g2*g2*g2+6.0*g1*g1*(g2*g2-2.0*la2)-36.0*g2*g2*la2
115 -8.0*(6.0*Yb2*Yb2*Yb2*Yb2+6.0*Yt*Yt*Yt*Yt+2.0*Ytau2*Ytau2*Ytau2*Ytau2-6.0*Yb2*Yb2*la2
116 -6.0*Yt*Yt*la2-2.0*Ytau2*Ytau2*la2-6.0*la2*la2-2.0*la3*la3-2.0*la3*la4-la4*la4-la5*la5))/(64.0*pi*pi);
117 //beta_lambda_3
118 beta[11] = (3.0*g1*g1*g1*g1+9.0*g2*g2*g2*g2-36.0*g2*g2*la3-6.0*g1*g1*(g2*g2+2.0*la3)
119 +8.0*(3.0*Yb2*Yb2*la3+3.0*Yt*Yt*la3+Ytau1*Ytau1*la3+Ytau2*Ytau2*la3+3.0*la1*la3+3.0*la2*la3+2.0*la3*la3
120 +Yb1*Yb1*(-6.0*Yt*Yt+3.0*la3)+la1*la4+la2*la4+la4*la4+la5*la5))/(64.0*pi*pi);
121 //beta_lambda_4
122 beta[12] = (3.0*g1*g1*(g2*g2-la4)-9.0*g2*g2*la4+6.0*Yb2*Yb2*la4+6.0*Yt*Yt*la4+2.0*Ytau1*Ytau1*la4
123 +2.0*Ytau2*Ytau2*la4+2.0*la1*la4+2.0*la2*la4+8.0*la3*la4+4.0*la4*la4+6.0*Yb1*Yb1*(2.0*Yt*Yt+la4)
124 +8.0*la5*la5)/(16.0*pi*pi);
125 //beta_lambda_5
126 beta[13] = (-3.0*g1*g1-9.0*g2*g2
127 +2.0*(3.0*Yb1*Yb1+3.0*Yb2*Yb2+3.0*Yt*Yt+Ytau1*Ytau1+Ytau2*Ytau2+la1+la2+4.0*la3+6.0*la4))*la5/(16.0*pi*pi);
128
129 if(ord==1||ord==2){
130 //beta_g1
131 beta[0] += (g1*g1*g1*(88.0*g3*g3-17.0*Yt*Yt))/(1536.0*pi*pi*pi*pi);
132 //beta_g2
133 beta[1] += (3.0*g2*g2*g2*(8.0*g3*g3-Yt*Yt))/(512.0*pi*pi*pi*pi);
134 //beta_g3
135 beta[2] += -((g3*g3*g3*(13.0*g3*g3+Yt*Yt))/(128.0*pi*pi*pi*pi));
136 //beta_Yt
137 beta[3] += (Yt*(-216.0*g3*g3*g3*g3+72.0*g3*g3*Yt*Yt-24.0*Yt*Yt*Yt*Yt-12.0*Yt*Yt*la2
138 +3.0*la2*la2+2.0*la3*la3+2.0*la3*la4+2.0*la4*la4+3.0*la5*la5))/(512.0*pi*pi*pi*pi);
139 //beta_Yb
140 beta[4] += (-1296.0*g3*g3*g3*g3*(Yb1+Yb2)+16.0*g3*g3*(4.0*Yb1+3.0*Yb2)*Yt*Yt
141 -3.0*(2.0*Yb1*(5.0*Yt*Yt*Yt*Yt-3.0*la1*la1-2.0*la3*la3
142 +4.0*Yt*Yt*(la3-la4)-2.0*la3*la4
143 -2.0*la4*la4-3.0*la5*la5)
144 +Yb2*(Yt*Yt*Yt*Yt-2.0*(3.0*la2*la2+2.0*la3*la3+2.0*la3*la4
145 +2.0*la4*la4+3.0*la5*la5))))
146 /(3072.0*pi*pi*pi*pi);
147 //beta_Ytau
148 beta[5] += (80.0*g3*g3*Yt*Yt*Ytau2-27.0*Yt*Yt*Yt*Yt*Ytau2+6.0*Ytau1*la1*la1+6.0*Ytau2*la2*la2
149 +(Ytau1+Ytau2)*(4.0*la3*la3+4.0*la3*la4+4.0*la4*la4+6.0*la5*la5))
150 /(1024.0*pi*pi*pi*pi);
151 //beta_m11_2
152 beta[6] += -(8.0*m22_2*(2.0*la3*la3+2.0*la3*la4+2.0*la4*la4
153 +3.0*Yt*Yt*(2.0*la3+la4)+3.0*la5*la5)
154 +m11_2*(30.0*la1*la1+4.0*la3*la3+4.0*la3*la4
155 +4.0*la4*la4+6.0*la5*la5))/(512.0*pi*pi*pi*pi);
156 //beta_m22_2
157 beta[7] += -((-80.0*g3*g3*m22_2*Yt*Yt
158 +8.0*m11_2*(2.0*la3*la3+2.0*la3*la4+2.0*la4*la4+3.0*la5*la5)
159 +m22_2*(27.0*Yt*Yt*Yt*Yt+72.0*Yt*Yt*la2+30.0*la2*la2+4.0*la3*la3
160 +4.0*la3*la4+4.0*la4*la4+6.0*la5*la5))/(512.0*pi*pi*pi*pi));
161 //beta_m12_2
162 beta[8] += (m12_2*(72.0*(la1*la1+la2*la2-4.0*la3*la4-8.0*la3*la5-8.0*la4*la5
163 +2.0*la5*la5-4.0*(la1+la2)*(la3+la4+la5))
164 +Yt*Yt*(960.0*g3*g3
165 -36.0*(9.0*Yt*Yt+8.0*(la3+2.0*la4+3.0*la5)))))
166 /(12288.0*pi*pi*pi*pi);
167 //beta_lambda_1
168 beta[9] += -(39.0*la1*la1*la1
169 +la1*(10.0*la3*la3+10.0*la3*la4+6.0*la4*la4+7.0*la5*la5)
170 +2.0*(4.0*la3*la3*la3+6.0*la3*la3*la4+8.0*la3*la4*la4
171 +3.0*la4*la4*la4+10.0*la3*la5*la5+11.0*la4*la5*la5
172 +3.0*Yt*Yt*(2.0*la3*la3+2.0*la3*la4+la4*la4+la5*la5)))
173 /(128.0*pi*pi*pi*pi);
174 //beta_lambda_2
175 beta[10] += -(-60.0*pow(Yt,6)+3.0*Yt*Yt*Yt*Yt*la2+72.0*Yt*Yt*la2*la2
176 +16.0*g3*g3*(4.0*Yt*Yt*Yt*Yt-5.0*Yt*Yt*la2)
177 +2.0*(39.0*la2*la2*la2+10.0*la2*la3*la3+8.0*la3*la3*la3
178 +10.0*la2*la3*la4+12.0*la3*la3*la4+6.0*la2*la4*la4
179 +16.0*la3*la4*la4+6.0*la4*la4*la4+7.0*la2*la5*la5
180 +20.0*la3*la5*la5+22.0*la4*la5*la5))/(256.0*pi*pi*pi*pi);
181 //beta_lambda_3
182 beta[11] += -(-80.0*g3*g3*Yt*Yt*la3+27.0*Yt*Yt*Yt*Yt*la3
183 +12.0*Yt*Yt*(2.0*la3*la3+la4*la4+2.0*la2*(3.0*la3+la4)+la5*la5)
184 +2.0*(la1*la1*(15.0*la3+4.0*la4)+la2*la2*(15.0*la3+4.0*la4)
185 +2.0*(la1+la2)*(18.0*la3*la3+8.0*la3*la4+7.0*la4*la4+9.0*la5*la5)
186 +2.0*(6.0*la3*la3*la3+2.0*la3*la3*la4+8.0*la3*la4*la4
187 +6.0*la4*la4*la4+9.0*la3*la5*la5+22.0*la4*la5*la5)))
188 /(512.0*pi*pi*pi*pi);
189 //beta_lambda_4
190 beta[12] += -(-80.0*g3*g3*Yt*Yt*la4+27.0*Yt*Yt*Yt*Yt*la4
191 +24.0*Yt*Yt*(la2*la4+2.0*la3*la4+la4*la4+2.0*la5*la5)
192 +2.0*(7.0*la1*la1*la4+7.0*la2*la2*la4+28.0*la3*la3*la4
193 +28.0*la3*la4*la4+48.0*la3*la5*la5+26.0*la4*la5*la5
194 +4.0*(la1+la2)*(10.0*la3*la4+5.0*la4*la4+6.0*la5*la5))
195 )/(512.0*pi*pi*pi*pi);
196 //beta_lambda_5
197 beta[13] += (la5*(80.0*g3*g3*Yt*Yt-3.0*Yt*Yt*Yt*Yt-24.0*Yt*Yt*(la2+2.0*la3+3.0*la4)
198 -2.0*(7.0*la1*la1+7.0*la2*la2+40.0*la1*la3+40.0*la2*la3+28.0*la3*la3
199 +44.0*la1*la4+44.0*la2*la4+76.0*la3*la4+32.0*la4*la4-6.0*la5*la5))
200 )/(512 *pi*pi*pi*pi);
201
202 if(ord==2){
203 //beta_g1
204 beta[0] += (g1*g1*g1*(208.0*g1*g1+108.0*g2*g2-15.0*Yb1*Yb1-15.0*Yb2*Yb2
205 -45.0*Ytau1*Ytau1-45.0*Ytau2*Ytau2))/(4608.0*pi*pi*pi*pi);
206 //beta_g2
207 beta[1] += (g2*g2*g2*(4.0*g1*g1+16.0*g2*g2-3.0*Yb1*Yb1-3.0*Yb2*Yb2
208 -Ytau1*Ytau1-Ytau2*Ytau2))/(512.0*pi*pi*pi*pi);
209 //beta_g3
210 beta[2] += (g3*g3*g3*(11.0*g1*g1+3.0*(9.0*g2*g2-4.0*(Yb1*Yb1+Yb2*Yb2))))/(1536.0*pi*pi*pi*pi);
211 //beta_Yt
212 beta[3] += (2534.0*g1*g1*g1*g1
213 +g1*g1*(-324.0*g2*g2+912.0*g3*g3-123.0*Yb1*Yb1+63.0*Yb2*Yb2
214 +3537.0*Yt*Yt+1350.0*Ytau2*Ytau2)
215 -9.0*(252.0*g2*g2*g2*g2
216 -9.0*g2*g2*(48.0*g3*g3+11.0*Yb1*Yb1+33.0*Yb2*Yb2+75.0*Yt*Yt+10.0*Ytau2*Ytau2)
217 -4.0*(16.0*g3*g3*(4.0*Yb1*Yb1+3.0*Yb2*Yb2)
218 -3.0*(10.0*Yb1*Yb1*Yb1*Yb1+Yb2*Yb2*Yb2*Yb2
219 +Yb2*Yb2*(11.0*Yt*Yt-5.0*Ytau2*Ytau2)
220 +9.0*Ytau2*Ytau2*(Yt*Yt+Ytau2*Ytau2)
221 +Yb1*Yb1*(10.0*Yt*Yt+3.0*Ytau1*Ytau1+8.0*la3-8.0*la4)))))
222 *Yt/(110592.0*pi*pi*pi*pi);
223 //beta_Yb
224 beta[4] += -(226.0*g1*g1*g1*g1*(Yb1+Yb2)
225 +3.0*g1*g1*(-711.0*Yb1*Yb1*Yb1-711.0*Yb2*Yb2*Yb2+324.0*g2*g2*(Yb1+Yb2)
226 -496.0*g3*g3*(Yb1+Yb2)+53.0*Yb1*Yt*Yt-273.0*Yb2*Yt*Yt
227 -450.0*Yb1*Ytau1*Ytau1-450.0*Yb2*Ytau2*Ytau2)
228 +27.0*(84.0*g2*g2*g2*g2*(Yb1+Yb2)
229 -3.0*g2*g2*(75.0*Yb1*Yb1*Yb1+75.0*Yb2*Yb2*Yb2+48.0*g3*g3*(Yb1+Yb2)
230 +11.0*Yb1*Yt*Yt+33.0*Yb2*Yt*Yt+10.0*Yb1*Ytau1*Ytau1
231 +10.0*Yb2*Ytau2*Ytau2)
232 +4.0*(48.0*pow(Yb1,5)-144.0*g3*g3*(Yb1*Yb1*Yb1+Yb2*Yb2*Yb2)
233 +3.0*Yb1*(3.0*Ytau1*Ytau1*Ytau1*Ytau1+Yt*Yt*Ytau2*Ytau2)
234 +Yb1*Yb1*Yb1*(10.0*Yt*Yt+9.0*Ytau1*Ytau1+24.0*la1)
235 +Yb2*(48.0*Yb2*Yb2*Yb2*Yb2-5.0*Yt*Yt*Ytau2*Ytau2
236 +9.0*Ytau2*Ytau2*Ytau2*Ytau2
237 +Yb2*Yb2*(11.0*Yt*Yt+9.0*Ytau2*Ytau2+24.0*la2)))))
238 /(110592.0*pi*pi*pi*pi);
239 //beta_Ytau
240 beta[5] += (966.0*g1*g1*g1*g1*(Ytau1+Ytau2)
241 +g1*g1*(50.0*Yb1*Yb1*Ytau1+537.0*Ytau1*Ytau1*Ytau1+50.0*Yb2*Yb2*Ytau2
242 +170.0*Yt*Yt*Ytau2+537.0*Ytau2*Ytau2*Ytau2+108.0*g2*g2*(Ytau1+Ytau2))
243 -3.0*(84.0*g2*g2*g2*g2*(Ytau1+Ytau2)
244 -15.0*g2*g2*(6.0*Yb1*Yb1*Ytau1+11.0*Ytau1*Ytau1*Ytau1+6.0*Yb2*Yb2*Ytau2
245 +6.0*Yt*Yt*Ytau2+11.0*Ytau2*Ytau2*Ytau2)
246 +4.0*(-80.0*g3*g3*(Yb1*Yb1*Ytau1+Yb2*Yb2*Ytau2)
247 +3.0*(9.0*Yb1*Yb1*Yb1*Yb1*Ytau1+4.0*pow(Ytau1,5)+9.0*Yb2*Yb2*Yb2*Yb2*Ytau2
248 -2.0*Yb2*Yb2*Yt*Yt*Ytau2+9.0*Yb2*Yb2*Ytau2*Ytau2*Ytau2
249 +9.0*Yt*Yt*Ytau2*Ytau2*Ytau2+4.0*pow(Ytau2,5)
250 +3.0*Yb1*Yb1*(3.0*Ytau1*Ytau1*Ytau1+Yt*Yt*(Ytau1+Ytau2))
251 +8.0*Ytau1*Ytau1*Ytau1*la1+8.0*Ytau2*Ytau2*Ytau2*la2))))
252 /(12288.0*pi*pi*pi*pi);
253 //beta_m11_2
254 beta[6] += (3.0*g1*g1*g1*g1*(193.0*m11_2+40.0*m22_2)
255 +2.0*g1*g1*(45.0*g2*g2*m11_2
256 +2.0*(m11_2*(25.0*Yb1*Yb1+75.0*Ytau1*Ytau1+144.0*la1)
257 +48.0*m22_2*(2.0*la3+la4)))
258 -3.0*(3.0*g2*g2*g2*g2*(41.0*m11_2-40.0*m22_2)
259 -12.0*g2*g2*(m11_2*(15.0*Yb1*Yb1+5.0*Ytau1*Ytau1+48.0*la1)
260 +16.0*m22_2*(2.0*la3+la4))
261 +8.0*(-80.0*g3*g3*m11_2*Yb1*Yb1
262 +3.0*m11_2*(9.0*Yb1*Yb1*Yb1*Yb1+3.0*Ytau1*Ytau1*Ytau1*Ytau1
263 +8.0*Ytau1*Ytau1*la1+3.0*Yb1*Yb1*(Yt*Yt+8.0*la1))
264 +8.0*m22_2*(3.0*Yb2*Yb2+Ytau2*Ytau2)*(2.0*la3+la4))))
265 /(12288.0*pi*pi*pi*pi);
266 //beta_m22_2
267 beta[7] += (3.0*g1*g1*g1*g1*(40.0*m11_2+193.0*m22_2)
268 +2.0*g1*g1*(45.0*g2*g2*m22_2
269 +2.0*(m22_2*(25.0*Yb2*Yb2+85.0*Yt*Yt+75.0*Ytau2*Ytau2+144.0*la2)
270 +48.0*m11_2*(2.0*la3+la4)))
271 +3.0*(3.0*g2*g2*g2*g2*(40.0*m11_2-41.0*m22_2)
272 +12.0*g2*g2*(m22_2*(15.0*Yb2*Yb2+15.0*Yt*Yt+5.0*Ytau2*Ytau2+48.0*la2)
273 +16.0*m11_2*(2.0*la3+la4))
274 +8.0*(80.0*g3*g3*m22_2*Yb2*Yb2
275 -3.0*m22_2*(9.0*Yb2*Yb2*Yb2*Yb2+3.0*Yb1*Yb1*Yt*Yt
276 +3.0*Ytau2*Ytau2*Ytau2*Ytau2+8.0*Ytau2*Ytau2*la2
277 +2.0*Yb2*Yb2*(7.0*Yt*Yt+12.0*la2))
278 -8.0*m11_2*(3.0*Yb1*Yb1+Ytau1*Ytau1)*(2.0*la3+la4))))
279 /(12288.0*pi*pi*pi*pi);
280 //beta_m12_2
281 beta[8] += (9.0*(51.0*g1*g1*g1*g1+10.0*g1*g1*g2*g2-81.0*g2*g2*g2*g2)
282 -324.0*(Yb1*Yb1*Yb1*Yb1+Yb2*Yb2*Yb2*Yb2)
283 +2.0*(85.0*g1*g1+135.0*g2*g2-432.0*Yb1*Yb1)*Yt*Yt
284 -108.0*(Ytau1*Ytau1*Ytau1*Ytau1+Ytau2*Ytau2*Ytau2*Ytau2)
285 +2.0*(Yb1*Yb1+Yb2*Yb2)*(25.0*g1*g1
286 +3.0*(45.0*g2*g2+4.0*(40.0*g3*g3+3.0*Yt*Yt-12.0*la3
287 -24.0*la4-36.0*la5)))
288 +192.0*(g1*g1+3.0*g2*g2)*(la3+2.0*la4+3.0*la5)
289 +6.0*(Ytau1*Ytau1+Ytau2*Ytau2)*(25.0*g1*g1+15.0*g2*g2
290 -16.0*(la3+2.0*la4+3.0*la5)))
291 *m12_2/(12288.0*pi*pi*pi*pi);
292 //beta_lambda_1
293 beta[9] += (-393.0*pow(g1,6)
294 +g1*g1*g1*g1*(-573.0*g2*g2+60.0*Yb1*Yb1-300.0*Ytau1*Ytau1
295 +651.0*la1+120.0*la3+60.0*la4)
296 +3.0*(291.0*pow(g2,6)
297 -3.0*g2*g2*g2*g2*(12.0*Yb1*Yb1+4.0*Ytau1*Ytau1+17.0*la1-40.0*la3-20.0*la4)
298 +12.0*g2*g2*(15.0*Yb1*Yb1*la1+5.0*Ytau1*Ytau1*la1
299 +4.0*(9.0*la1*la1+(2.0*la3+la4)*(2.0*la3+la4)))
300 -8.0*(-60.0*pow(Yb1,6)-20.0*pow(Ytau1,6)+Ytau1*Ytau1*Ytau1*Ytau1*la1
301 +24.0*Ytau1*Ytau1*la1*la1+3.0*Yb1*Yb1*Yb1*Yb1*(-4.0*Yt*Yt+la1)
302 +9.0*Yb1*Yb1*la1*(Yt*Yt+8.0*la1)
303 +16.0*g3*g3*(4.0*Yb1*Yb1*Yb1*Yb1-5.0*Yb1*Yb1*la1)
304 +24.0*Yb2*Yb2*la3*la3+8.0*Ytau2*Ytau2*la3*la3+24.0*Yb2*Yb2*la3*la4
305 +8.0*Ytau2*Ytau2*la3*la4+12.0*Yb2*Yb2*la4*la4+4.0*Ytau2*Ytau2*la4*la4
306 +12.0*Yb2*Yb2*la5*la5+4.0*Ytau2*Ytau2*la5*la5))
307 +g1*g1*(-303.0*g2*g2*g2*g2
308 +6.0*g2*g2*(36.0*Yb1*Yb1+44.0*Ytau1*Ytau1+39.0*la1+20.0*la4)
309 +4.0*(16.0*Yb1*Yb1*Yb1*Yb1+25.0*Yb1*Yb1*la1
310 +3.0*(-16.0*Ytau1*Ytau1*Ytau1*Ytau1+25.0*Ytau1*Ytau1*la1+36.0*la1*la1
311 +16.0*la3*la3+16.0*la3*la4+8.0*la4*la4-4.0*la5*la5))))
312 /(6144.0*pi*pi*pi*pi);
313 //beta_lambda_2
314 beta[10] += (-393.0*pow(g1,6)
315 +g1*g1*g1*g1*(-573.0*g2*g2+60.0*Yb2*Yb2-228.0*Yt*Yt-300.0*Ytau2*Ytau2
316 +651.0*la2+120.0*la3+60.0*la4)
317 +g1*g1*(-303.0*g2*g2*g2*g2
318 +6.0*g2*g2*(36.0*Yb2*Yb2+84.0*Yt*Yt+44.0*Ytau2*Ytau2+39.0*la2+20.0*la4)
319 +4.0*(16.0*Yb2*Yb2*Yb2*Yb2-32.0*Yt*Yt*Yt*Yt-48.0*Ytau2*Ytau2*Ytau2*Ytau2
320 +25.0*Yb2*Yb2*la2+85.0*Yt*Yt*la2+75.0*Ytau2*Ytau2*la2
321 +108.0*la2*la2+48.0*la3*la3+48.0*la3*la4+24.0*la4*la4-12.0*la5*la5))
322 +3.0*(291.0*pow(g2,6)
323 -3.0*g2*g2*g2*g2*(12.0*Yb2*Yb2+12.0*Yt*Yt+4.0*Ytau2*Ytau2
324 +17.0*la2-40.0*la3-20.0*la4)
325 +12.0*g2*g2*(15.0*Yb2*Yb2*la2+15.0*Yt*Yt*la2+5.0*Ytau2*Ytau2*la2
326 +36.0*la2*la2+16.0*la3*la3+16.0*la3*la4+4.0*la4*la4)
327 -8.0*(-60.0*pow(Yb2,6)-12.0*Yb1*Yb1*Yt*Yt*Yt*Yt-20.0*pow(Ytau2,6)
328 +9.0*Yb1*Yb1*Yt*Yt*la2+Ytau2*Ytau2*Ytau2*Ytau2*la2
329 +24.0*Ytau2*Ytau2*la2*la2+3.0*Yb2*Yb2*Yb2*Yb2*(4.0*Yt*Yt+la2)
330 +16.0*g3*g3*(4.0*Yb2*Yb2*Yb2*Yb2-5.0*Yb2*Yb2*la2)
331 +6.0*Yb2*Yb2*(2.0*Yt*Yt*Yt*Yt+7.0*Yt*Yt*la2+12.0*la2*la2)
332 +24.0*Yb1*Yb1*la3*la3+8.0*Ytau1*Ytau1*la3*la3+24.0*Yb1*Yb1*la3*la4
333 +8.0*Ytau1*Ytau1*la3*la4+12.0*Yb1*Yb1*la4*la4+4.0*Ytau1*Ytau1*la4*la4
334 +12.0*Yb1*Yb1*la5*la5+4.0*Ytau1*Ytau1*la5*la5)))/(6144.0*pi*pi*pi*pi);
335 //beta_lambda_3
336 beta[11] += (-393.0*pow(g1,6)
337 +3.0*g1*g1*g1*g1*(101.0*g2*g2+10.0*Yb1*Yb1+10.0*Yb2*Yb2-38.0*Yt*Yt-50.0*Ytau1*Ytau1
338 -50.0*Ytau2*Ytau2+30.0*la1+30.0*la2+197.0*la3+20.0*la4)
339 +g1*g1*(33.0*g2*g2*g2*g2-6.0*g2*g2*(18.0*Yb1*Yb1+18.0*Yb2*Yb2+42.0*Yt*Yt
340 +22.0*Ytau1*Ytau1+22.0*Ytau2*Ytau2+10.0*la1
341 +10.0*la2-11.0*la3+12.0*la4)
342 +2.0*(25.0*Yb2*Yb2*la3+85.0*Yt*Yt*la3+75.0*Ytau1*Ytau1*la3
343 +75.0*Ytau2*Ytau2*la3+144.0*la1*la3+144.0*la2*la3+24*la3*la3
344 +Yb1*Yb1*(-16.0*Yt*Yt+25.0*la3)
345 +48.0*la1*la4+48.0*la2*la4-24.0*la4*la4+48.0*la5*la5))
346 +3.0*(291.0*pow(g2,6)-3.0*g2*g2*g2*g2*(6.0*Yb1*Yb1+6.0*Yb2*Yb2+6.0*Yt*Yt
347 +2.0*Ytau1*Ytau1+2.0*Ytau2*Ytau2
348 -30.0*la1-30.0*la2+37.0*la3-20.0*la4)
349 +6.0*g2*g2*(15.0*Yb1*Yb1*la3+15.0*Yb2*Yb2*la3+15.0*Yt*Yt*la3
350 +5.0*Ytau1*Ytau1*la3+5.0*Ytau2*Ytau2*la3+48.0*la1*la3
351 +48.0*la2*la3+8.0*la3*la3+24.0*la1*la4+24.0*la2*la4
352 -16.0*la3*la4+8.0*la4*la4)
353 -4.0*(-9.0*Yb1*Yb1*Yb1*Yb1*(8.0*Yt*Yt-3.0*la3)+27.0*Yb2*Yb2*Yb2*Yb2*la3
354 +42.0*Yb2*Yb2*Yt*Yt*la3+9.0*Ytau1*Ytau1*Ytau1*Ytau1*la3
355 +9.0*Ytau2*Ytau2*Ytau2*Ytau2*la3+24.0*Ytau1*Ytau1*la1*la3
356 +72.0*Yb2*Yb2*la2*la3+24.0*Ytau2*Ytau2*la2*la3+24.0*Yb2*Yb2*la3*la3
357 +8.0*Ytau1*Ytau1*la3*la3+8.0*Ytau2*Ytau2*la3*la3
358 +16.0*g3*g3*(Yb1*Yb1*(8.0*Yt*Yt-5.0*la3)-5.0*Yb2*Yb2*la3)
359 +48.0*Yb2*Yb2*Yt*Yt*la4+8.0*Ytau1*Ytau1*la1*la4+24.0*Yb2*Yb2*la2*la4
360 +8.0*Ytau2*Ytau2*la2*la4+12.0*Yb2*Yb2*la4*la4+4.0*Ytau1*Ytau1*la4*la4
361 +4.0*Ytau2*Ytau2*la4*la4+12.0*Yb2*Yb2*la5*la5+4.0*Ytau1*Ytau1*la5*la5
362 +4.0*Ytau2*Ytau2*la5*la5
363 -6.0*Yb1*Yb1*(12.0*Yt*Yt*Yt*Yt+5.0*Yt*Yt*la3
364 -2.0*(6.0*la1*la3+2.0*la3*la3+2.0*la1*la4+la4*la4+la5*la5)))))
365 /(6144.0*pi*pi*pi*pi);
366 //beta_lambda_4
367 beta[12] += (g1*g1*g1*g1*(-876.0*g2*g2+471.0*la4)
368 +2.0*g1*g1*(-168.0*g2*g2*g2*g2+25.0*Yb2*Yb2*la4+85.0*Yt*Yt*la4
369 +75.0*Ytau1*Ytau1*la4+75.0*Ytau2*Ytau2*la4
370 +48.0*la1*la4+48.0*la2*la4+48.0*la3*la4+96.0*la4*la4
371 +Yb1*Yb1*(16.0*Yt*Yt+25.0*la4)
372 +3.0*g2*g2*(36.0*Yb1*Yb1+36.0*Yb2*Yb2+84.0*Yt*Yt
373 +44.0*Ytau1*Ytau1+44.0*Ytau2*Ytau2
374 +20.0*la1+20.0*la2+8.0*la3+51.0*la4)
375 +192.0*la5*la5)
376 -3.0*(231.0*g2*g2*g2*g2*la4-90.0*g2*g2*Yb2*Yb2*la4
377 +108.0*Yb2*Yb2*Yb2*Yb2*la4-90.0*g2*g2*Yt*Yt*la4
378 -216.0*Yb2*Yb2*Yt*Yt*la4-30.0*g2*g2*Ytau1*Ytau1*la4
379 +36.0*Ytau1*Ytau1*Ytau1*Ytau1*la4-30.0*g2*g2*Ytau2*Ytau2*la4
380 +36.0*Ytau2*Ytau2*Ytau2*Ytau2*la4+32.0*Ytau1*Ytau1*la1*la4
381 +96.0*Yb2*Yb2*la2*la4+32.0*Ytau2*Ytau2*la2*la4-288.0*g2*g2*la3*la4
382 +192.0*Yb2*Yb2*la3*la4+64.0*Ytau1*Ytau1*la3*la4+64.0*Ytau2*Ytau2*la3*la4
383 -144.0*g2*g2*la4*la4+96.0*Yb2*Yb2*la4*la4+32.0*Ytau1*Ytau1*la4*la4
384 +32.0*Ytau2*Ytau2*la4*la4+12.0*Yb1*Yb1*Yb1*Yb1*(16.0*Yt*Yt+9.0*la4)
385 -64.0*g3*g3*(5.0*Yb2*Yb2*la4+Yb1*Yb1*(8.0*Yt*Yt+5.0*la4))
386 -432.0*g2*g2*la5*la5+192.0*Yb2*Yb2*la5*la5+64.0*Ytau1*Ytau1*la5*la5
387 +64.0*Ytau2*Ytau2*la5*la5
388 +6.0*Yb1*Yb1*(32.0*Yt*Yt*Yt*Yt-15.0*g2*g2*la4
389 +4.0*Yt*Yt*(8.0*la3+11.0*la4)
390 +16.0*(la1*la4+2.0*la3*la4+la4*la4+2.0*la5*la5))))
391 /(6144.0*pi*pi*pi*pi);
392 //beta_lambda_5
393 beta[13] += (471.0*g1*g1*g1*g1
394 +2.0*g1*g1*(57.0*g2*g2+25.0*Yb1*Yb1+25.0*Yb2*Yb2+85.0*Yt*Yt+75.0*Ytau1*Ytau1
395 +75.0*Ytau2*Ytau2-24.0*la1-24.0*la2+192.0*la3+288.0*la4)
396 -3.0*(231.0*g2*g2*g2*g2
397 -6.0*g2*g2*(15.0*Yb1*Yb1+15.0*Yb2*Yb2+15.0*Yt*Yt+5.0*Ytau1*Ytau1
398 +5.0*Ytau2*Ytau2+48.0*la3+96.0*la4)
399 +4.0*(3.0*Yb1*Yb1*Yb1*Yb1+3.0*Yb2*Yb2*Yb2*Yb2-80.0*g3*g3*(Yb1*Yb1+Yb2*Yb2)
400 -6.0*Yb2*Yb2*Yt*Yt+Ytau1*Ytau1*Ytau1*Ytau1+Ytau2*Ytau2*Ytau2*Ytau2
401 +8.0*Ytau1*Ytau1*la1+24.0*Yb2*Yb2*la2+8.0*Ytau2*Ytau2*la2+48.0*Yb2*Yb2*la3
402 +16.0*Ytau1*Ytau1*la3+16.0*Ytau2*Ytau2*la3+72.0*Yb2*Yb2*la4
403 +24.0*Ytau1*Ytau1*la4+24.0*Ytau2*Ytau2*la4
404 +6.0*Yb1*Yb1*(11.0*Yt*Yt+4.0*la1+8.0*la3+12.0*la4))))
405 *la5/(6144.0*pi*pi*pi*pi);
406 }
407 }
408
409 return 0;
410}
411
412int Jacobian (double t, const double y[], double *dfdy, double dfdt[], void *order)
413{
414 return 0;
415}
416
417int RGEcheck(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
418{
419 int check=0;
420
421 //perturbativity of the Yukawa couplings
422 for(int i=3;i<6;i++)
423 {
424 if(fabs(InitialValues[i])>sqrt(4.0*M_PI)) check=1;
425 }
426 //perturbativity of the quartic Higgs couplings
427 for(int i=9;i<14;i++)
428 {
429 if(fabs(InitialValues[i])>4.0*M_PI) check=1;
430 }
431 //positivity
432 if(InitialValues[9]<0.0) check=1;
433 if(InitialValues[10]<0.0) check=1;
434 if(InitialValues[11]+sqrt(fabs(InitialValues[9]*InitialValues[10]))<0.0) check=1;
435 if(InitialValues[11]+InitialValues[12]-fabs(InitialValues[13])+sqrt(fabs(InitialValues[9]*InitialValues[10]))<0.0) check=1;
436 //NLO unitarity
437 double YtQ = InitialValues[3];
438 if(t1>tNLOuni)
439 {
440 double la1Q = InitialValues[9];
441 double la2Q = InitialValues[10];
442 double la3Q = InitialValues[11];
443 double la4Q = InitialValues[12];
444 double la5Q = InitialValues[13];
445
446 double betalambda1 = 6.0*la1Q*la1Q + 2.0*la3Q*la3Q + 2.0*la3Q*la4Q + la4Q*la4Q + la5Q*la5Q;
447// + 6.0*la1Q*Yb1Q*Yb1Q + 2.0*la1Q*Ytau1Q*Ytau1Q
448// - 6.0*Yb1Q*Yb1Q*Yb1Q*Yb1Q - 2.0*Ytau1Q*Ytau1Q*Ytau1Q*Ytau1Q;
449 double betalambda2 = 6.0*la2Q*la2Q + 2.0*la3Q*la3Q + 2.0*la3Q*la4Q + la4Q*la4Q + la5Q*la5Q + 6.0*la2Q*YtQ*YtQ - 6.0*YtQ*YtQ*YtQ*YtQ;
450// + 6.0*la2Q*Yb2Q*Yb2Q + 2.0*la2Q*Ytau2Q*Ytau2Q
451// - 6.0*Yb2Q*Yb2Q*Yb2Q*Yb2Q - 2.0*Ytau2Q*Ytau2Q*Ytau2Q*Ytau2Q;
452 double betalambda3 = 3.0*la1Q*la3Q + 3.0*la2Q*la3Q + 2.0*la3Q*la3Q + la1Q*la4Q + la2Q*la4Q + la4Q*la4Q + la5Q*la5Q + 3.0*la3Q*YtQ*YtQ;
453// + 3.0*la3Q*Yb1Q*Yb1Q + 3.0*la3Q*Yb2Q*Yb2Q + la3Q*Ytau1Q*Ytau1Q + la3Q*Ytau2Q*Ytau2Q
454// - 6.0*Yb1Q*Yb1Q*(Yb2Q*Yb2Q + YtQ*YtQ) - 2.0*Ytau1Q*Ytau1Q*Ytau2Q*Ytau2Q;
455 double betalambda4 = la1Q*la4Q + la2Q*la4Q + 4.0*la3Q*la4Q + 2.0*la4Q*la4Q + 4.0*la5Q*la5Q + 3.0*la4Q*YtQ*YtQ;
456// + 3.0*la4Q*Yb1Q*Yb1Q + 3.0*la4Q*Yb2Q*Yb2Q + la4Q*Ytau1Q*Ytau1Q + la4Q*Ytau2Q*Ytau2Q
457// - 2.0*Ytau1Q*Ytau1Q*Ytau2Q*Ytau2Q - 6.0*Yb1Q*Yb1Q*(Yb2Q*Yb2Q - YtQ*YtQ);
458 double betalambda5 = la5Q*(la1Q + la2Q + 4.0*la3Q + 6.0*la4Q + 3.0*YtQ*YtQ);
459// + 3.0*la5Q*Yb1Q*Yb1Q + la5Q*Ytau1Q*Ytau1Q + la5Q*(3.0*Yb2Q*Yb2Q + Ytau2Q*Ytau2Q)
460// - 6.0*Yb1Q*Yb1Q*Yb2Q*Yb2Q - 2.0*Ytau1Q*Ytau1Q*Ytau2Q*Ytau2Q;
461
462 double uniLO1 = -3.0*la1Q/(16.0*M_PI);
463 gslpp::complex uniNLO1 = -3.0*la1Q/(16.0*M_PI) +9.0*betalambda1/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(9.0*la1Q*la1Q+(2.0*la3Q+la4Q)*(2.0*la3Q+la4Q))/(256.0*M_PI*M_PI*M_PI);
464 double uniLO2 = -3.0*la2Q/(16.0*M_PI);
465 gslpp::complex uniNLO2 = -3.0*la2Q/(16.0*M_PI) +9.0*betalambda2/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(9.0*la2Q*la2Q+(2.0*la3Q+la4Q)*(2.0*la3Q+la4Q))/(256.0*M_PI*M_PI*M_PI);
466 double uniLO3 = -(2.0*la3Q+la4Q)/(16.0*M_PI);
467 gslpp::complex uniNLO3 = -(2.0*la3Q+la4Q)/(16.0*M_PI) +3.0*(2.0*betalambda3+betalambda4)/(256.0*M_PI*M_PI*M_PI) +3.0*(gslpp::complex::i()*M_PI-1.0)*(la1Q+la2Q)*(2.0*la3Q+la4Q)/(256.0*M_PI*M_PI*M_PI);
468 double uniLO4 = -(la3Q+2.0*la4Q)/(16.0*M_PI);
469 gslpp::complex uniNLO4 = -(la3Q+2.0*la4Q)/(16.0*M_PI) +3.0*(betalambda3+2.0*betalambda4)/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la3Q*la3Q+4.0*la3Q*la4Q+4.0*la4Q*la4Q+9.0*la5Q*la5Q)/(256.0*M_PI*M_PI*M_PI);
470 double uniLO5 = -3.0*la5Q/(16.0*M_PI);
471 gslpp::complex uniNLO5 = -3.0*la5Q/(16.0*M_PI) +9.0*betalambda5/(256.0*M_PI*M_PI*M_PI) +3.0*(gslpp::complex::i()*M_PI-1.0)*(la3Q+2.0*la4Q)*la5Q/(128.0*M_PI*M_PI*M_PI);
472 double uniLO6 = -la1Q/(16.0*M_PI);
473 gslpp::complex uniNLO6 = -la1Q/(16.0*M_PI) +3.0*betalambda1/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la1Q*la1Q+la4Q*la4Q)/(256.0*M_PI*M_PI*M_PI);
474 double uniLO7 = -la2Q/(16.0*M_PI);
475 gslpp::complex uniNLO7 = -la2Q/(16.0*M_PI) +3.0*betalambda2/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la2Q*la2Q+la4Q*la4Q)/(256.0*M_PI*M_PI*M_PI);
476 double uniLO8 = -la4Q/(16.0*M_PI);
477 gslpp::complex uniNLO8 = -la4Q/(16.0*M_PI) +3.0*betalambda4/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la1Q+la2Q)*la4Q/(256.0*M_PI*M_PI*M_PI);
478 double uniLO10 = -la3Q/(16.0*M_PI);
479 gslpp::complex uniNLO10 = -la3Q/(16.0*M_PI) +3.0*betalambda3/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la3Q*la3Q+la5Q*la5Q)/(256.0*M_PI*M_PI*M_PI);
480 double uniLO11 = -la5Q/(16.0*M_PI);
481 gslpp::complex uniNLO11 = -la5Q/(16.0*M_PI) +3.0*betalambda5/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*la3Q*la5Q/(128.0*M_PI*M_PI*M_PI);
482 double uniLO14 = -(la3Q-la4Q)/(16.0*M_PI);
483 gslpp::complex uniNLO14 = -(la3Q-la4Q)/(16.0*M_PI) +3.0*(betalambda3-betalambda4)/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la3Q-la4Q)*(la3Q-la4Q)/(256.0*M_PI*M_PI*M_PI);
484 double uniLO15 = -la1Q/(16.0*M_PI);
485 gslpp::complex uniNLO15 = -la1Q/(16.0*M_PI) +3.0*betalambda1/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la1Q*la1Q+la5Q*la5Q)/(256.0*M_PI*M_PI*M_PI);
486 double uniLO16 = -la2Q/(16.0*M_PI);
487 gslpp::complex uniNLO16 = -la2Q/(16.0*M_PI) +3.0*betalambda2/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la2Q*la2Q+la5Q*la5Q)/(256.0*M_PI*M_PI*M_PI);
488 double uniLO17 = -la5Q/(16.0*M_PI);
489 gslpp::complex uniNLO17 = -la5Q/(16.0*M_PI) +3.0*betalambda5/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la1Q+la2Q)*la5Q/(256.0*M_PI*M_PI*M_PI);
490 double uniLO24 = -(la3Q+la4Q)/(16.0*M_PI);
491 gslpp::complex uniNLO24 = -(la3Q+la4Q)/(16.0*M_PI) +3.0*(betalambda3+betalambda4)/(256.0*M_PI*M_PI*M_PI) +(gslpp::complex::i()*M_PI-1.0)*(la3Q+la4Q)*(la3Q+la4Q)/(256.0*M_PI*M_PI*M_PI);
492
493 double uniLOev1 = 0.5*(uniLO1+uniLO2+sqrt((uniLO1-uniLO2)*(uniLO1-uniLO2)+4.0*uniLO3*uniLO3));
494 double uniLOev2 = 0.5*(uniLO1+uniLO2-sqrt((uniLO1-uniLO2)*(uniLO1-uniLO2)+4.0*uniLO3*uniLO3));
495 double uniLOev3 = uniLO4+uniLO5;
496 double uniLOev4 = uniLO4-uniLO5;
497 double uniLOev5 = 0.5*(uniLO6+uniLO7+sqrt((uniLO6-uniLO7)*(uniLO6-uniLO7)+4.0*uniLO8*uniLO8));
498 double uniLOev6 = 0.5*(uniLO6+uniLO7-sqrt((uniLO6-uniLO7)*(uniLO6-uniLO7)+4.0*uniLO8*uniLO8));
499 double uniLOev9 = uniLO10+uniLO11;
500 double uniLOev10 = uniLO10-uniLO11;
501 double uniLOev13 = 0.5*(uniLO15+uniLO16+sqrt((uniLO15-uniLO16)*(uniLO15-uniLO16)+4.0*uniLO17*uniLO17));
502 double uniLOev14 = 0.5*(uniLO15+uniLO16-sqrt((uniLO15-uniLO16)*(uniLO15-uniLO16)+4.0*uniLO17*uniLO17));
503 gslpp::complex uniNLOev1wowfr = 0.5*(uniNLO1+uniNLO2+sqrt((uniNLO1-uniNLO2)*(uniNLO1-uniNLO2)+4.0*uniNLO3*uniNLO3));
504 gslpp::complex uniNLOev2wowfr = 0.5*(uniNLO1+uniNLO2-sqrt((uniNLO1-uniNLO2)*(uniNLO1-uniNLO2)+4.0*uniNLO3*uniNLO3));
505 gslpp::complex uniNLOev3wowfr = uniNLO4+uniNLO5;
506 gslpp::complex uniNLOev4wowfr = uniNLO4-uniNLO5;
507 gslpp::complex uniNLOev5wowfr = 0.5*(uniNLO6+uniNLO7+sqrt((uniNLO6-uniNLO7)*(uniNLO6-uniNLO7)+4.0*uniNLO8*uniNLO8));
508 gslpp::complex uniNLOev6wowfr = 0.5*(uniNLO6+uniNLO7-sqrt((uniNLO6-uniNLO7)*(uniNLO6-uniNLO7)+4.0*uniNLO8*uniNLO8));
509 gslpp::complex uniNLOev9wowfr = uniNLO10+uniNLO11;
510 gslpp::complex uniNLOev10wowfr = uniNLO10-uniNLO11;
511 gslpp::complex uniNLOev13wowfr = 0.5*(uniNLO15+uniNLO16+sqrt((uniNLO15-uniNLO16)*(uniNLO15-uniNLO16)+4.0*uniNLO17*uniNLO17));
512 gslpp::complex uniNLOev14wowfr = 0.5*(uniNLO15+uniNLO16-sqrt((uniNLO15-uniNLO16)*(uniNLO15-uniNLO16)+4.0*uniNLO17*uniNLO17));
513
514 if( (uniNLOev1wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
515 if( (uniNLOev2wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
516 if( (uniNLOev3wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
517 if( (uniNLOev4wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
518 if( (uniNLOev5wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
519 if( (uniNLOev6wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
520 if( (uniNLOev9wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
521 if( (uniNLOev10wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
522 if( (uniNLOev13wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
523 if( (uniNLOev14wowfr - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
524 if( (uniNLO14 - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
525 if( (uniNLO24 - 0.5*gslpp::complex::i()).abs() > 0.5) check=1;
526
527 if( fabs(uniLOev1) > Rpeps && (uniNLOev1wowfr/uniLOev1-1.0).abs() > 1.0) check=1;
528 if( fabs(uniLOev2) > Rpeps && (uniNLOev2wowfr/uniLOev2-1.0).abs() > 1.0) check=1;
529 if( fabs(uniLOev3) > Rpeps && (uniNLOev3wowfr/uniLOev3-1.0).abs() > 1.0) check=1;
530 if( fabs(uniLOev4) > Rpeps && (uniNLOev4wowfr/uniLOev4-1.0).abs() > 1.0) check=1;
531 if( fabs(uniLOev5) > Rpeps && (uniNLOev5wowfr/uniLOev5-1.0).abs() > 1.0) check=1;
532 if( fabs(uniLOev6) > Rpeps && (uniNLOev6wowfr/uniLOev6-1.0).abs() > 1.0) check=1;
533 if( fabs(uniLOev9) > Rpeps && (uniNLOev9wowfr/uniLOev9-1.0).abs() > 1.0) check=1;
534 if( fabs(uniLOev10) > Rpeps && (uniNLOev10wowfr/uniLOev10-1.0).abs() > 1.0) check=1;
535 if( fabs(uniLOev13) > Rpeps && (uniNLOev13wowfr/uniLOev13-1.0).abs() > 1.0) check=1;
536 if( fabs(uniLOev14) > Rpeps && (uniNLOev14wowfr/uniLOev14-1.0).abs() > 1.0) check=1;
537 if( fabs(uniLO14) > Rpeps && (uniNLO14/uniLO14-1.0).abs() > 1.0) check=1;
538 if( fabs(uniLO24) > Rpeps && (uniNLO24/uniLO24-1.0).abs() > 1.0) check=1;
539
540 }
541
542 return check;
543}
544
545double Runner::RGERunner(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double NLOuniscale)
546{
547 //Define which stepping function should be used
548 const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk4;
549
550 //Allocate space for the stepping function
551 gsl_odeiv2_step * s = gsl_odeiv2_step_alloc(T, NumberOfRGEs);
552
553 //Define the absolute (A) and relative (R) error on y at each step.
554 //The real error will be compared to the following error estimate:
555 // A + R * |y_i|
556 gsl_odeiv2_control * c = gsl_odeiv2_control_y_new(1e-6, 0.0);
557
558 //Allocate space for the evolutor
559 gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(NumberOfRGEs);
560
561 //Definition of the RGE system (the Jacobian is not necessary for the RK4 method; it's an empty function here)
562 gsl_odeiv2_system RGEsystem = {RGEs, Jacobian, NumberOfRGEs, &order};
563
564 //Set starting and end point as natural logarithmic scales (conversion from decadic log scale)
565 double t1 = Q1*log(10.0);
566 double t2 = Q2*log(10.0);
567 double tNLOuni = NLOuniscale*log(10.0);
568
569 //Set initial step size
570 double InitialStepSize = 1e-6;
571
572 //Run!
573 while (t1 < t2)
574 {
575 int status = gsl_odeiv2_evolve_apply (e, c, s, &RGEsystem, &t1, t2, &InitialStepSize, InitialValues);
576 if(status != GSL_SUCCESS) break;
577
578 //intermediate checks if appropriate
579 if(RGEcheck(InitialValues,t1,Rpeps,tNLOuni) != 0) break;
580 }
581
582 gsl_odeiv2_evolve_free (e);
583 gsl_odeiv2_control_free (c);
584 gsl_odeiv2_step_free (s);
585
586 //Return the decadic log scale at which the evolution stopped
587 return t1/log(10.0);
588}
RGEs
int RGEs(double t, const double y[], double beta[], void *flags)
Definition: Runner.cpp:42
RGEcheck
int RGEcheck(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: Runner.cpp:417
Jacobian
int Jacobian(double t, const double y[], double *dfdy, double dfdt[], void *order)
Definition: Runner.cpp:412
Runner::RGERunner
virtual double RGERunner(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double tNLOuni)
Definition: Runner.cpp:545
Runner::mylambda3
lambda3 * mylambda3
Definition: Runner.h:55
Runner::mylambda4
lambda4 * mylambda4
Definition: Runner.h:56
Runner::mylambda5
lambda5 * mylambda5
Definition: Runner.h:57
Runner::computeThValue
double computeThValue()
Empty function.
Definition: Runner.cpp:37
Runner::mym22_2
m22_2 * mym22_2
Definition: Runner.h:52
Runner::mylambda1
lambda1 * mylambda1
Definition: Runner.h:53
Runner::mylambda2
lambda2 * mylambda2
Definition: Runner.h:54
Runner::Runner
Runner(const StandardModel &SM_i)
Runner constructor.
Definition: Runner.cpp:15
Runner::mym11_2
m11_2 * mym11_2
Definition: Runner.h:51
Runner::~Runner
~Runner()
Runner destructor.
Definition: Runner.cpp:26
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:509
THDM
A base class for symmetric Two-Higgs-Doublet models.
Definition: THDM.h:120
ThObservable
A class for a model prediction of an observable.
Definition: ThObservable.h:25
lambda1
parameter of the Higgs potential
Definition: THDMquantities.h:382
lambda2
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:405
lambda3
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:428
lambda4
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:451
lambda5
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:474
m11_2
An observable class for the quadratic Higgs potential coupling .
Definition: THDMquantities.h:336
m22_2
An observable class for the quadratic Higgs potential coupling .
Definition: THDMquantities.h:359
s
Test Observable.
t
Test Observable.