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 
15 Runner::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 
26 Runner::~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 
37 double Runner::computeThValue()
38 {
39  return 0.;
40 }
41 
42 int 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 
412 int Jacobian (double t, const double y[], double *dfdy, double dfdt[], void *order)
413 {
414  return 0;
415 }
416 
417 int 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 
545 double 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 }
Runner::mylambda4
lambda4 * mylambda4
Definition: Runner.h:56
m11_2
An observable class for the quadratic Higgs potential coupling .
Definition: THDMquantities.h:336
Jacobian
int Jacobian(double t, const double y[], double *dfdy, double dfdt[], void *order)
Definition: Runner.cpp:412
THDM
A base class for symmetric Two-Higgs-Doublet models.
Definition: THDM.h:120
lambda3
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:428
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:477
lambda1
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:382
Runner::mylambda2
lambda2 * mylambda2
Definition: Runner.h:54
lambda4
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:451
Runner::computeThValue
double computeThValue()
Empty function.
Definition: Runner.cpp:37
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
Runner::mylambda5
lambda5 * mylambda5
Definition: Runner.h:57
RGEcheck
int RGEcheck(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: Runner.cpp:417
THDMquantities.h
lambda2
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:405
Runner::~Runner
~Runner()
Runner destructor.
Definition: Runner.cpp:26
RGEs
int RGEs(double t, const double y[], double beta[], void *flags)
Definition: Runner.cpp:42
m22_2
An observable class for the quadratic Higgs potential coupling .
Definition: THDMquantities.h:359
Runner::Runner
Runner(const StandardModel &SM_i)
Runner constructor.
Definition: Runner.cpp:15
Runner::mylambda1
lambda1 * mylambda1
Definition: Runner.h:53
ThObservable
A class for a model prediction of an observable.
Definition: ThObservable.h:25
lambda5
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:474
Runner::mylambda3
lambda3 * mylambda3
Definition: Runner.h:55
THDM.h
Runner.h
Runner::mym11_2
m11_2 * mym11_2
Definition: Runner.h:51
Runner::mym22_2
m22_2 * mym22_2
Definition: Runner.h:52
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