a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
RunnerTHDMW.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 "RunnerTHDMW.h"
9 #include <gsl/gsl_errno.h>
10 #include <gsl/gsl_matrix.h>
11 #include <gsl/gsl_odeiv2.h>
12 
13 RunnerTHDMW::RunnerTHDMW(const StandardModel& SM_i) : myTHDMW(static_cast<const THDMW*> (&SM_i))
14 {}
15 
16 RunnerTHDMW::~RunnerTHDMW()
17 {};
18 
19 int RGEsTHDMW(double t, const double y[], double beta[], void *flags)
20 {
21  (void)(t); /* avoid unused parameter warning */
22  int ord = *(int *)flags;
23 // int type = flag-ord;
24 // double Yb1=0;
25 // double Yb2=0;
26 // double Ytau1=0;
27 // double Ytau2=0;
28  double la1=y[0];
29  double la2=y[1];
30  double la3=y[2];
31  double la4=y[3];
32  double mu1=y[4];
33  double mu3=y[5];
34  double mu4=y[6];
35  double nu1=y[7];
36  double om1=y[8];
37  double ka1=y[9];
38  double nu2=y[10];
39  double om2=y[11];
40  double ka2=y[12];
41  double nu4=y[13];
42  double om4=y[14];
43 
44  double pi=M_PI;
45 
46  //beta_la1
47  beta[0] = (12.0*la1*la1 + 4.0*la3*la3 + 4.0*la3*la4 + 4.0*la4*la4
48  + 8.0*nu1*nu1 + 8.0*nu1*nu2 + 8.0*nu2*nu2)/(16.0*pi*pi);
49  //beta_la2
50  beta[1] = (12.0*la2*la2 + 4.0*la3*la3 + 4.0*la3*la4 + 4.0*la4*la4
51  + 8.0*om1*om1 + 8.0*om1*om2 + 8.0*om2*om2)/(16.0*pi*pi);
52  //beta_la3
53  beta[2] = (4.0*la3*la3 + 4.0*la4*la4 + (la1+la2)*(6.0*la3+2.0*la4)
54  + 8.0*ka2*ka2 + 8.0*nu1*om1 + 4.0*nu2*om1 + 4.0*nu1*om2)/(16.0*pi*pi);
55  //beta_la4
56  beta[3] = (la1*la4 + la2*la4 + 4.0*la3*la4 + 6.0*la4*la4
57  + 4.0*ka1*ka1 + 4.0*ka1*ka2 + 2.0*ka2*ka2 + 2.0*nu2*om2)/(8.0*pi*pi);
58  //beta_mu1
59  beta[4] = (11.0*mu1*mu1 + 3.0*mu1*mu4 + mu1*(2.0*mu1+6.0*mu3+3.0*mu4)
60  + 3.0*nu4*nu4 + 3.0*om4*om4)/(16.0*pi*pi);
61  //beta_mu3
62  beta[5] = (18.0*ka1*ka1 + 18.0*ka1*ka2 + 134.0*mu1*mu1 + 6.0*mu1*(39.0*mu3 + 22.0*mu4)
63  + 3.0*(30.0*mu3*mu3 + 39.0*mu3*mu4 + 9.0*mu4*mu4
64  + 3.0*nu1*nu1 + 3.0*nu1*nu2 - 5.0*nu4*nu4
65  + 3.0*om1*om1 + 3.0*om1*om2 - 5.0*om4*om4))/(72.0*pi*pi);
66  //beta_mu4
67  beta[6] = (18.0*ka2*ka2 + 4.0*mu1*mu1 + 156.0*mu1*mu4 + 54.0*mu3*mu4 + 144.0*mu4*mu4
68  + 9.0*nu2*nu2 + 6.0*nu4*nu4 + 9.0*om2*om2 + 6.0*om4*om4)/(144.0*pi*pi);
69  //beta_nu1
70  beta[7] = (6.0*ka1*ka1 + 6.0*ka2*ka2 + 18.0*la1*nu1
71  + 78.0*mu1*nu1 + 51.0*mu3*nu1 + 39.0*mu4*nu1 + 6.0*nu1*nu1
72  + 6.0*la1*nu2 + 32.0*mu1*nu2 + 24.0*mu3*nu2 + 6.0*mu4*nu2
73  + 6.0*nu2*nu2 + 10.0*nu4*nu4
74  + 12.0*la3*om1 + 6.0*la4*om1 + 6.0*la3*om2)/(48.0*pi*pi);
75  //beta_om1
76  beta[8] = (6.0*ka1*ka1 + 6.0*ka2*ka2
77  + 12.0*la3*nu1 + 6.0*la4*nu1 + 6.0*la3*nu2
78  + 18.0*la2*om1 + 78.0*mu1*om1 + 51.0*mu3*om1 + 39.0*mu4*om1 + 6.0*om1*om1
79  + 6.0*la2*om2 + 32.0*mu1*om2 + 24.0*mu3*om2 + 6.0*mu4*om2 + 6.0*om2*om2
80  + 10.0*om4*om4)/(48.0*pi*pi);
81  //beta_ka1
82  beta[9] = (6.0*ka1*(2.0*la3 + 10.0*la4 + 18.0*mu1 + 17.0*mu3 + 13.0*mu4 + 2.0*nu1 + 2.0*om1)
83  + ka2*(24.0*la4 + 64.0*mu1 + 48.0*mu3 + 24.0*mu4 + 9.0*nu2 + 9.0*om2)
84  + 20.0*nu4*om4)/(96.0*pi*pi);
85  //beta_nu2
86  beta[10] = (4.0*ka1*ka2 + 6.0*ka2*ka2 + 2.0*la1*nu2 + ((14.0*mu1)/3.0 + mu3 + 9.0*mu4)*nu2
87  + 4.0*nu1*nu2 + 6.0*nu2*nu2 + (25.0*nu4*nu4)/3.0 + 2.0*la4*om2)/(16.0*pi*pi);
88  //beta_om2
89  beta[11] = (4.0*ka1*ka2 + 6.0*ka2*ka2 + 2.0*la4*nu2 + 2.0*la2*om2
90  + ((14.0*mu1)/3.0 + mu3 + 9.0*mu4)*om2 + 4.0*om1*om2 + 6.0*om2*om2
91  + (25.0*om4*om4)/3.0)/(16.0*pi*pi);
92  //beta_ka2
93  beta[12] = (ka2*(6.0*la3 + 6.0*la4 + 14.0*mu1 + 3.0*mu3 + 27.0*mu4
94  + 6.0*nu1 + 12.0*nu2 + 6.0*om1 + 12.0*om2)
95  + 6.0*ka1*(nu2 + om2) + 42.0*nu4*om4)/(48.0*pi*pi);
96  //beta_nu4
97  beta[13] = (11.0*mu1*nu4 + 3.0*mu3*nu4 + 9.0*mu4*nu4 + 3.0*nu1*nu4 + 9.0*nu2*nu4
98  + 3.0*ka1*om4 + 9.0*ka2*om4)/(16.0*pi*pi);
99  //beta_om4
100  beta[14] = (3.0*ka1*nu4 + 9.0*ka2*nu4
101  + (11.0*mu1 + 3.0*(mu3 + 3.0*mu4 + om1 + 3.0*om2))*om4)/(16.0*pi*pi);
102 
103  if(ord==1){
104  //beta_la1
105  beta[0] += 0.0;
106  //beta_la2
107  beta[1] += 0.0;
108  //beta_la3
109  beta[2] += 0.0;
110  //beta_la4
111  beta[3] += 0.0;
112  //beta_mu1
113  beta[4] += 0.0;
114  //beta_mu3
115  beta[5] += 0.0;
116  //beta_mu4
117  beta[6] += 0.0;
118  //beta_nu1
119  beta[7] += 0.0;
120  //beta_om1
121  beta[8] += 0.0;
122  //beta_ka1
123  beta[9] += 0.0;
124  //beta_nu2
125  beta[10] += 0.0;
126  //beta_om2
127  beta[11] += 0.0;
128  //beta_ka2
129  beta[12] += 0.0;
130  //beta_nu4
131  beta[13] += 0.0;
132  //beta_om4
133  beta[14] += 0.0;
134  }
135 
136  return 0;
137 }
138 
139 int RGEsMW(double t, const double y[], double beta[], void *flags)
140 {
141  (void)(t); /* avoid unused parameter warning */
142  int ord = *(int *)flags;
143 // int type = flag-ord;
144 // double Yb1=0;
145 // double Yb2=0;
146 // double Ytau1=0;
147 // double Ytau2=0;
148  double la1=y[0];
149  double nu1=y[1];
150  double nu2=y[2];
151  double nu3=y[3];
152  double nu4=y[4];
153  double nu5=y[5];
154  double mu1=y[6];
155  double mu2=y[7];
156  double mu3=y[8];
157  double mu4=y[9];
158  double mu5=y[10];
159  double mu6=y[11];
160 
161  double pi=M_PI;
162 
163  //RGE taken from 1303.4848
164  //The lambda1 running is taken from (4.1) in 1303.4848. The rest stems from the appendix.
165 
166  //beta_la1
167  beta[0] = (12.0*la1*la1 + 8.0*nu1*nu1 + 8.0*nu1*nu2 + 4.0*nu2*nu2 + 16.0*nu3*nu3)/(16.0*pi*pi);
168  //beta_nu1
169  beta[1] = (2.0*nu1*nu1 + nu2*nu2 + 4.0*nu3*nu3 + 2.0*la1*(3.0*nu1+nu2)
170  + (7.0*nu4*nu4 - 4.0*nu4*nu5 + 7.0*nu5*nu5)/3.0
171  + nu1*(8.0*mu1 + 8.0*mu2 + 17.0*mu3 + 10.0*mu4 + 3.0*mu5 + 5.0*mu6)
172  + nu2*(8.0*mu1 + 8.0*mu2 + 24.0*mu3 + 3.0*mu4
173  + 3.0*mu5 + 8.0*mu6)/3.0)/(16.0*pi*pi);
174  //beta_nu2
175  beta[2] = (2.0*nu2*nu2 + 4.0*nu1*nu2 + 16.0*nu3*nu3 + 2.0*la1*nu2
176  + (4.0*nu4*nu4 + 17.0*nu4*nu5 + 4.0*nu5*nu5)/3.0
177  + nu2*(8.0*mu1 + 8.0*mu2 + 3.0*mu3 + 24.0*mu4
178  + 3.0*mu5 - mu6)/3.0)/(16.0*pi*pi);
179  //beta_nu3
180  beta[3] = (2.0*nu3*(la1 + 2.0*nu1 + 3.0*nu2)
181  + (17.0*nu4*nu4 + 16.0*nu4*nu5 + 17.0*nu5*nu5)/12.0
182  + nu3*(-mu1 - mu2 + 3.0*mu3 + 3.0*mu4
183  + 24.0*mu5 + 8.0*mu6)/3.0)/(16.0*pi*pi);
184  //beta_nu4
185  beta[4] = (8.0*nu3*nu4 + 2.0*nu3*nu5
186  + nu5*(2.0*nu2 - mu2 + 2.0*mu4 + 4.0*mu5 + mu6)
187  + nu4*(3.0*nu1 + 2.0*nu2 + 6.0*mu1 + 2.0*mu2 + 3.0*mu3
188  + 2.0*mu4 + mu5 + mu6))/(16.0*pi*pi);
189  //beta_nu5
190  beta[5] = (2.0*nu3*nu4 + 8.0*nu3*nu5
191  + nu4*(2.0*nu2 - mu1 + 2.0*mu4 + 4.0*mu5 + mu6)
192  + nu5*(3.0*nu1 + 2.0*nu2 + 6.0*mu1 + 2.0*mu2 + 3.0*mu3
193  + 2.0*mu4 + mu5 + mu6))/(16.0*pi*pi);
194  //beta_mu1
195  beta[6] = (3.0*nu4*nu4 + 7.0*mu1*mu1
196  + mu1*(6.0*mu2 + 6.0*mu3 + 4.0*mu4 - mu5 - 2.0*mu6)
197  + mu2*(4.0*mu4 - mu5)
198  - 2.0*mu4*mu6 + 2.0*mu5*mu6 + mu6*mu6)/(16.0*pi*pi);
199  //beta_mu2
200  beta[7] = (3.0*nu5*nu5 + 7.0*mu2*mu2
201  + mu2*(6.0*mu1 + 6.0*mu3 + 4.0*mu4 - mu5 - 2.0*mu6)
202  + mu1*(4.0*mu4 - mu5)
203  - 2.0*mu4*mu6 + 2.0*mu5*mu6 + mu6*mu6)/(16.0*pi*pi);
204  //beta_mu3
205  beta[8] = (20.0*mu3*mu3
206  + mu3*(288.0*mu1 + 288.0*mu2 + 360.0*mu4 + 108.0*mu5 + 180.0*mu6)/18.0
207  + (36.0*nu1*nu1 + 36.0*nu1*nu2 - 24.0*nu4*nu4 - 12.0*nu4*nu5
208  - 24.0*nu5*nu5 + 62.0*mu1*mu1 + 64.0*mu1*mu2 + 62.0*mu2*mu2
209  + (96.0*mu4 + 18.0*mu5 + 58.0*mu6)*(mu1 + mu2)
210  + 54.0*mu4*mu4 + 36.0*mu4*mu5 + 132.0*mu4*mu6 + 18.0*mu5*mu5
211  + 18.0*mu5*mu6 + 29.0*mu6*mu6)/18.0)/(16.0*pi*pi);
212  //beta_mu4
213  beta[9] = (nu2*nu2 - (nu4*nu4 - 4.0*nu4*nu5 + nu5*nu5)/3.0 + 10.0*mu4*mu4 /*mu4??*/
214  + mu5*(mu1 + mu2 + mu6)
215  + mu4*(4.0*(4.0*mu1 + 4.0*mu2 + mu6)/3.0 + 2.0*mu5 + 6.0*mu4)
216  + 4.0*mu5*mu5
217  + (mu1*mu1 + mu2*mu2 - 4.0*mu6*(mu1+mu2) - 2.0*mu6*mu6)/9.0
218  + 26.0/9.0*mu1*mu2)/(16.0*pi*pi);
219  //beta_mu5
220  beta[10] = (4.0*nu3*nu3 - (nu4*nu4 - 4.0*nu4*nu5 + nu5*nu5)/3.0
221  + mu5*((mu1 + mu2 + 19.0*mu6)/3.0 + 8.0*mu4 + 6.0*mu3)
222  + 2.0*mu4*mu6 + 8.0*mu5*mu5
223  + (mu1*mu1 + mu2*mu2 - 4.0*mu6*(mu1+mu2) + 7.0*mu6*mu6)/9.0
224  - 10.0/9.0*mu1*mu2)/(16.0*pi*pi);
225  //beta_mu6
226  beta[11] = (0.5*mu6*mu6 + 3.0*nu4*nu4 + 3.0*nu5*nu5
227  - 2.0*(mu1*mu1 + mu2*mu2) + 6.0*mu5*(mu1 + mu2)
228  + 7.0*mu6*(mu1 + mu2 + mu3))/(16.0*pi*pi);
229 
230  if(ord==1){
231  //beta_la1
232  beta[0] += 0.0;
233  //beta_nu1
234  beta[1] += 0.0;
235  //beta_nu2
236  beta[2] += 0.0;
237  //beta_nu3
238  beta[3] += 0.0;
239  //beta_nu4
240  beta[4] += 0.0;
241  //beta_nu5
242  beta[5] += 0.0;
243  //beta_mu1
244  beta[6] += 0.0;
245  //beta_mu2
246  beta[7] += 0.0;
247  //beta_mu3
248  beta[8] += 0.0;
249  //beta_mu4
250  beta[9] += 0.0;
251  //beta_mu5
252  beta[10] += 0.0;
253  //beta_mu6
254  beta[11] += 0.0;
255  }
256 
257  return 0;
258 }
259 
260 int RGEscustodialMW(double t, const double y[], double beta[], void *flags)
261 {
262  (void)(t); /* avoid unused parameter warning */
263  int ord = *(int *)flags;
264 // int type = flag-ord;
265 // double Yb1=0;
266 // double Yb2=0;
267 // double Ytau1=0;
268 // double Ytau2=0;
269  double la1=y[0];
270  double nu1=y[1];
271  double nu2=y[2];
272  double nu4=y[3];
273  double mu1=y[4];
274  double mu3=y[5];
275  double mu4=y[6];
276 
277  double pi=M_PI;
278 
279  //RGE taken from 1303.4848
280 
281  //beta_la1
282  //This was taken from the custodial THDMW!
283  beta[0] = (12.0*la1*la1 + 8.0*nu1*nu1 + 8.0*nu1*nu2 + 8.0*nu2*nu2)/(16.0*pi*pi);
284  //beta_nu1
285  beta[1] = (2.0*nu1*nu1 + 2.0*nu2*nu2 + 2.0*la1*(3.0*nu1+nu2)
286  + 10.0*nu4*nu4/3.0
287  + nu1*(26.0*mu1 + 17.0*mu3 + 13.0*mu4)
288  + nu2*(32.0*mu1 + 24.0*mu3 + 6.0*mu4)/3.0)/(16.0*pi*pi);
289  //beta_nu2
290  beta[2] = (4.0*nu1*nu2 + 6.0*nu2*nu2 + 2.0*la1*nu2 + 25.0*nu4*nu4/3.0
291  + nu2*(14.0*mu1 + 3.0*mu3 + 27.0*mu4)/3.0)/(16.0*pi*pi);
292  //beta_nu4
293  beta[3] = (3.0*nu1*nu4 + 9.0*nu2*nu4
294  + nu4*(11.0*mu1 + 3.0*mu3 + 9.0*mu4))/(16.0*pi*pi);
295  //beta_mu1
296  beta[4] = (3.0*nu4*nu4 + 13.0*mu1*mu1
297  + 6.0*mu1*(mu3+mu4))/(16.0*pi*pi);
298  //beta_mu3
299  beta[5] = (20.0*mu3*mu3
300  + mu3*(52.0*mu1 + 26.0*mu4)
301  + 2.0*nu1*nu1 + 2.0*nu1*nu2 - 10.0*nu4*nu4/3.0
302  + 268.0*mu1*mu1/9.0 + 88.0*mu1*mu4/3.0 + 6.0*mu4*mu4)/(16.0*pi*pi);
303  //beta_mu4
304  //This was taken from the custodial THDMW, because I think the formula from 1303.4848 is incorrect
305  beta[6] = (nu2*nu2 + 2.0*nu4*nu4/3.0 + 16.0*mu4*mu4
306  + 52.0*mu1*mu4/3.0 + 6.0*mu3*mu4
307  + 4.0/9.0*mu1*mu1)/(16.0*pi*pi);
308 
309  if(ord==1){
310  //beta_la1
311  beta[0] += 0.0;
312  //beta_nu1
313  beta[1] += 0.0;
314  //beta_nu2
315  beta[2] += 0.0;
316  //beta_nu4
317  beta[3] += 0.0;
318  //beta_mu1
319  beta[4] += 0.0;
320  //beta_mu3
321  beta[5] += 0.0;
322  //beta_mu4
323  beta[6] += 0.0;
324  }
325 
326  return 0;
327 }
328 
329 int JacobianTHDMW (double t, const double y[], double *dfdy, double dfdt[], void *order)
330 {
331  return 0;
332 }
333 
334 int RGEcheckTHDMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
335 {
336  int check=0;
337 
338  //perturbativity of the Yukawa couplings
339 // for(int i=3;i<6;i++)
340 // {
341 // if(fabs(InitialValues[i])>sqrt(4.0*M_PI)) check=1;
342 // }
343 
344  //perturbativity of the quartic Higgs couplings
345  for(int i=0;i<15;i++)
346  {
347  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
348  }
349 
350  double la1Q = InitialValues[0];
351  double la2Q = InitialValues[1];
352  double la3Q = InitialValues[2];
353  double la4Q = InitialValues[3];
354  double mu1Q = InitialValues[4];
355  double mu3Q = InitialValues[5];
356  double mu4Q = InitialValues[6];
357  double nu1Q = InitialValues[7];
358  double om1Q = InitialValues[8];
359  double ka1Q = InitialValues[9];
360  double nu2Q = InitialValues[10];
361  double om2Q = InitialValues[11];
362  double ka2Q = InitialValues[12];
363  double nu4Q = InitialValues[13];
364  double om4Q = InitialValues[14];
365 
366  //positivity checks
367  double muAtimes2=4.0*mu1Q+2.0*mu3Q+4.0*mu4Q;
368  if(la1Q<0.0) check=1;
369  if(la2Q<0.0) check=1;
370  if(muAtimes2<0.0) check=1;
371  if(5.0*mu1Q+3.0*mu3Q+3.0*mu4Q-fabs(mu1Q)<0.0) check=1;
372  if(sqrt(la1Q*muAtimes2)+nu1Q<0.0) check=1;
373  if(sqrt(la1Q*muAtimes2)+nu1Q+2.0*nu2Q<0.0) check=1;
374  if(la1Q+0.25*muAtimes2+nu1Q+2.0*nu2Q-2.0*fabs(nu4Q)/sqrt(3.0)<0.0) check=1;
375  if(la3Q+sqrt(la1Q*la2Q)<0.0) check=1;
376  if(la3Q+2.0*la4Q+sqrt(la1Q*la2Q)<0.0) check=1;
377  if(sqrt(la2Q*muAtimes2)+om1Q<0.0) check=1;
378  if(sqrt(la2Q*muAtimes2)+om1Q+2.0*om2Q<0.0) check=1;
379  if(la2Q+0.25*muAtimes2+om1Q+2.0*om2Q-2.0*fabs(om4Q)/sqrt(3.0)<0.0) check=1;
380 
382  //NLO unitarity//
384 
385  double pi=M_PI;
386  gslpp::matrix<gslpp::complex> Smatrix1(4,4,0.), Smatrix2(4,4,0.);
387  gslpp::matrix<gslpp::complex> Seigenvectors1(4,4,0.), Seigenvectors2(4,4,0.);
388  gslpp::vector<double> Seigenvalues1(4,0.), Seigenvalues2(4,0.);
389  gslpp::vector<gslpp::complex> unitarityeigenvalues(11,0.);
390 
391  if(t1>tNLOuni)
392  {
393 
394  //LO part
395  Smatrix1.assign(0,0, 3.0*la1Q/(16.0*pi));
396  Smatrix1.assign(0,1, (2.0*la3Q+la4Q)/(16.0*pi));
397  Smatrix1.assign(1,0, Smatrix1(0,1));
398  Smatrix1.assign(0,3, (2.0*nu1Q+nu2Q)/(8.0*sqrt(2.0)*pi));
399  Smatrix1.assign(3,0, Smatrix1(0,3));
400  Smatrix1.assign(1,1, 3.0*la2Q/(16.0*pi));
401  Smatrix1.assign(1,3, (2.0*om1Q+om2Q)/(8.0*sqrt(2.0)*pi));
402  Smatrix1.assign(3,1, Smatrix1(1,3));
403  Smatrix1.assign(2,2, (la3Q+5.0*la4Q)/(16.0*pi));
404  Smatrix1.assign(2,3, (4.0*ka1Q+2.0*ka2Q)/(16.0*pi));
405  Smatrix1.assign(3,2, Smatrix1(2,3));
406  Smatrix1.assign(3,3, (26.0*mu1Q+17.0*mu3Q+13.0*mu4Q)/(32.0*pi));
407 
408  Smatrix2.assign(0,0, la1Q/(16.0*pi));
409  Smatrix2.assign(0,1, la4Q/(16.0*pi));
410  Smatrix2.assign(1,0, Smatrix2(0,1));
411  Smatrix2.assign(0,3, nu2Q/(8.0*sqrt(2.0)*pi));
412  Smatrix2.assign(3,0, Smatrix2(0,3));
413  Smatrix2.assign(1,1, la2Q/(16.0*pi));
414  Smatrix2.assign(1,3, om2Q/(8.0*sqrt(2.0)*pi));
415  Smatrix2.assign(3,1, Smatrix2(1,3));
416  Smatrix2.assign(2,2, (la3Q+la4Q)/(16.0*pi));
417  Smatrix2.assign(2,3, ka2Q/(8.0*pi));
418  Smatrix2.assign(3,2, Smatrix2(2,3));
419  Smatrix2.assign(3,3, (14.0*mu1Q+3.0*mu3Q+27.0*mu4Q)/(96.0*pi));
420 
421  Smatrix1.eigensystem(Seigenvectors1, Seigenvalues1);
422  Smatrix2.eigensystem(Seigenvectors2, Seigenvalues2);
423 
424  for (int i=0; i < 4; i++) {
425  unitarityeigenvalues.assign(i, Seigenvalues1(i));
426  unitarityeigenvalues.assign(4+i, Seigenvalues2(i));
427  }
428  unitarityeigenvalues.assign(8, (la3Q-la4Q)/(16.0*pi));
429  unitarityeigenvalues.assign(9, sqrt(15.0)*nu4Q/(16.0*pi));
430  unitarityeigenvalues.assign(10, sqrt(15.0)*om4Q/(16.0*pi));
431 
432  //beta_la1*16pi^2
433  double betala1 = 12.0*la1Q*la1Q + 4.0*la3Q*la3Q + 4.0*la3Q*la4Q + 4.0*la4Q*la4Q
434  + 8.0*nu1Q*nu1Q + 8.0*nu1Q*nu2Q + 8.0*nu2Q*nu2Q;
435  //beta_la2*16pi^2
436  double betala2 = 12.0*la2Q*la2Q + 4.0*la3Q*la3Q + 4.0*la3Q*la4Q + 4.0*la4Q*la4Q
437  + 8.0*om1Q*om1Q + 8.0*om1Q*om2Q + 8.0*om2Q*om2Q;
438  //beta_la3*16pi^2
439  double betala3 = 4.0*la3Q*la3Q + 4.0*la4Q*la4Q + (la1Q+la2Q)*(6.0*la3Q+2.0*la4Q)
440  + 8.0*ka2Q*ka2Q + 8.0*nu1Q*om1Q + 4.0*nu2Q*om1Q + 4.0*nu1Q*om2Q;
441  //beta_la4*16pi^2
442  double betala4 = (la1Q*la4Q + la2Q*la4Q + 4.0*la3Q*la4Q + 6.0*la4Q*la4Q
443  + 4.0*ka1Q*ka1Q + 4.0*ka1Q*ka2Q + 2.0*ka2Q*ka2Q + 2.0*nu2Q*om2Q)*2.0;
444  //beta_mu1*16pi^2
445  double betamu1 = 11.0*mu1Q*mu1Q + 3.0*mu1Q*mu4Q + mu1Q*(2.0*mu1Q+6.0*mu3Q+3.0*mu4Q)
446  + 3.0*nu4Q*nu4Q + 3.0*om4Q*om4Q;
447  //beta_mu3*16pi^2
448  double betamu3 = (18.0*ka1Q*ka1Q + 18.0*ka1Q*ka2Q + 134.0*mu1Q*mu1Q + 6.0*mu1Q*(39.0*mu3Q + 22.0*mu4Q)
449  + 3.0*(30.0*mu3Q*mu3Q + 39.0*mu3Q*mu4Q + 9.0*mu4Q*mu4Q
450  + 3.0*nu1Q*nu1Q + 3.0*nu1Q*nu2Q - 5.0*nu4Q*nu4Q
451  + 3.0*om1Q*om1Q + 3.0*om1Q*om2Q - 5.0*om4Q*om4Q))/4.5;
452  //beta_mu4*16pi^2
453  double betamu4 = (18.0*ka2Q*ka2Q + 4.0*mu1Q*mu1Q + 156.0*mu1Q*mu4Q + 54.0*mu3Q*mu4Q + 144.0*mu4Q*mu4Q
454  + 9.0*nu2Q*nu2Q + 6.0*nu4Q*nu4Q + 9.0*om2Q*om2Q + 6.0*om4Q*om4Q)/9.0;
455  //beta_nu1*16pi^2
456  double betanu1 = (6.0*ka1Q*ka1Q + 6.0*ka2Q*ka2Q + 18.0*la1Q*nu1Q
457  + 78.0*mu1Q*nu1Q + 51.0*mu3Q*nu1Q + 39.0*mu4Q*nu1Q + 6.0*nu1Q*nu1Q
458  + 6.0*la1Q*nu2Q + 32.0*mu1Q*nu2Q + 24.0*mu3Q*nu2Q + 6.0*mu4Q*nu2Q
459  + 6.0*nu2Q*nu2Q + 10.0*nu4Q*nu4Q
460  + 12.0*la3Q*om1Q + 6.0*la4Q*om1Q + 6.0*la3Q*om2Q)/3.0;
461  //beta_om1*16pi^2
462  double betaom1 = (6.0*ka1Q*ka1Q + 6.0*ka2Q*ka2Q
463  + 12.0*la3Q*nu1Q + 6.0*la4Q*nu1Q + 6.0*la3Q*nu2Q
464  + 18.0*la2Q*om1Q + 78.0*mu1Q*om1Q + 51.0*mu3Q*om1Q + 39.0*mu4Q*om1Q + 6.0*om1Q*om1Q
465  + 6.0*la2Q*om2Q + 32.0*mu1Q*om2Q + 24.0*mu3Q*om2Q + 6.0*mu4Q*om2Q + 6.0*om2Q*om2Q
466  + 10.0*om4Q*om4Q)/3.0;
467  //beta_ka1*16pi^2
468  double betaka1 = (6.0*ka1Q*(2.0*la3Q + 10.0*la4Q + 18.0*mu1Q + 17.0*mu3Q + 13.0*mu4Q + 2.0*nu1Q + 2.0*om1Q)
469  + ka2Q*(24.0*la4Q + 64.0*mu1Q + 48.0*mu3Q + 24.0*mu4Q + 9.0*nu2Q + 9.0*om2Q)
470  + 20.0*nu4Q*om4Q)/6.0;
471  //beta_nu2*16pi^2
472  double betanu2 = 4.0*ka1Q*ka2Q + 6.0*ka2Q*ka2Q + 2.0*la1Q*nu2Q + ((14.0*mu1Q)/3.0 + mu3Q + 9.0*mu4Q)*nu2Q
473  + 4.0*nu1Q*nu2Q + 6.0*nu2Q*nu2Q + (25.0*nu4Q*nu4Q)/3.0 + 2.0*la4Q*om2Q;
474  //beta_om2*16pi^2
475  double betaom2 = 4.0*ka1Q*ka2Q + 6.0*ka2Q*ka2Q + 2.0*la4Q*nu2Q + 2.0*la2Q*om2Q
476  + ((14.0*mu1Q)/3.0 + mu3Q + 9.0*mu4Q)*om2Q + 4.0*om1Q*om2Q + 6.0*om2Q*om2Q
477  + (25.0*om4Q*om4Q)/3.0;
478  //beta_ka2*16pi^2
479  double betaka2 = (ka2Q*(6.0*la3Q + 6.0*la4Q + 14.0*mu1Q + 3.0*mu3Q + 27.0*mu4Q
480  + 6.0*nu1Q + 12.0*nu2Q + 6.0*om1Q + 12.0*om2Q)
481  + 6.0*ka1Q*(nu2Q + om2Q) + 42.0*nu4Q*om4Q)/3.0;
482  //beta_nu4*16pi^2
483  double betanu4 = 11.0*mu1Q*nu4Q + 3.0*mu3Q*nu4Q + 9.0*mu4Q*nu4Q + 3.0*nu1Q*nu4Q + 9.0*nu2Q*nu4Q
484  + 3.0*ka1Q*om4Q + 9.0*ka2Q*om4Q;
485  //beta_om4*16pi^2
486  double betaom4 = 3.0*ka1Q*nu4Q + 9.0*ka2Q*nu4Q
487  + (11.0*mu1Q + 3.0*(mu3Q + 3.0*mu4Q + om1Q + 3.0*om2Q))*om4Q;
488 
489 // diagonalization
490  gslpp::matrix<gslpp::complex> Sbmatrix1(4,4,0.), Sbmatrix2(4,4,0.);
491  gslpp::matrix<gslpp::complex> Seigenvectors1T(4,4,0.), Seigenvectors2T(4,4,0.);
492  gslpp::vector<gslpp::complex> Sbeigenvalues1(4,0.), Sbeigenvalues2(4,0.);
493  gslpp::vector<gslpp::complex> betaeigenvalues(11,0.);
494  gslpp::vector<gslpp::complex> NLOunitarityeigenvalues(11,0.);
495 
496  Sbmatrix1.assign(0,0, 3.0*betala1/(16.0*pi));
497  Sbmatrix1.assign(0,1, (2.0*betala3+betala4)/(16.0*pi));
498  Sbmatrix1.assign(1,0, Sbmatrix1(0,1));
499  Sbmatrix1.assign(0,3, (2.0*betanu1+betanu2)/(8.0*sqrt(2.0)*pi));
500  Sbmatrix1.assign(3,0, Sbmatrix1(0,3));
501  Sbmatrix1.assign(1,1, 3.0*betala2/(16.0*pi));
502  Sbmatrix1.assign(1,3, (2.0*betaom1+betaom2)/(8.0*sqrt(2.0)*pi));
503  Sbmatrix1.assign(3,1, Sbmatrix1(1,3));
504  Sbmatrix1.assign(2,2, (betala3+5.0*betala4)/(16.0*pi));
505  Sbmatrix1.assign(2,3, (4.0*betaka1+2.0*betaka2)/(16.0*pi));
506  Sbmatrix1.assign(3,2, Sbmatrix1(2,3));
507  Sbmatrix1.assign(3,3, (26.0*betamu1+17.0*betamu3+13.0*betamu4)/(32.0*pi));
508 
509  Sbmatrix2.assign(0,0, betala1/(16.0*pi));
510  Sbmatrix2.assign(0,1, betala4/(16.0*pi));
511  Sbmatrix2.assign(1,0, Sbmatrix2(0,1));
512  Sbmatrix2.assign(0,3, betanu2/(8.0*sqrt(2.0)*pi));
513  Sbmatrix2.assign(3,0, Sbmatrix2(0,3));
514  Sbmatrix2.assign(1,1, betala2/(16.0*pi));
515  Sbmatrix2.assign(1,3, betaom2/(8.0*sqrt(2.0)*pi));
516  Sbmatrix2.assign(3,1, Sbmatrix2(1,3));
517  Sbmatrix2.assign(2,2, (betala3+betala4)/(16.0*pi));
518  Sbmatrix2.assign(2,3, betaka2/(8.0*pi));
519  Sbmatrix2.assign(3,2, Sbmatrix2(2,3));
520  Sbmatrix2.assign(3,3, (14.0*betamu1+3.0*betamu3+27.0*betamu4)/(96.0*pi));
521 
522  Seigenvectors1T=Seigenvectors1.hconjugate();
523  Seigenvectors2T=Seigenvectors2.hconjugate();
524 
525  for (int i=0; i < 4; i++) {
526  for (int k=0; k < 4; k++) {
527  for (int l=0; l < 4; l++) {
528  Sbeigenvalues1.assign(i, Sbeigenvalues1(i) + Seigenvectors1T(i,k) * Sbmatrix1(k,l) * Seigenvectors1(l,i) );
529  Sbeigenvalues2.assign(i, Sbeigenvalues2(i) + Seigenvectors2T(i,k) * Sbmatrix2(k,l) * Seigenvectors2(l,i) );
530  }
531  }
532  betaeigenvalues.assign(i, -1.5 * Sbeigenvalues1(i));
533  betaeigenvalues.assign(i+4, -1.5 * Sbeigenvalues2(i));
534  }
535 
536  betaeigenvalues.assign(8, -1.5 * (betala3-betala4)/(16.0*pi));
537  betaeigenvalues.assign(9, -1.5 * sqrt(15.0)*betanu4/(16.0*pi));
538  betaeigenvalues.assign(10, -1.5 * sqrt(15.0)*betaom4/(16.0*pi));
539 
540  for (int i=0; i < 11; i++) {
541  NLOunitarityeigenvalues.assign(i, -(gslpp::complex::i()-1.0/pi)*unitarityeigenvalues(i)*unitarityeigenvalues(i) + betaeigenvalues(i) );
542  if( ( unitarityeigenvalues(i) + NLOunitarityeigenvalues(i).real() ).abs() > 0.5) check=1;
543  if( (unitarityeigenvalues(i)).abs() > Rpeps && (NLOunitarityeigenvalues(i)/unitarityeigenvalues(i)).abs() > 1.0) check=1;
544  }
545 
546  } //end of the if(t1>tNLOuni)
547  return check;
548 }
549 
550 int RGEcheckMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
551 {
552  int check=0;
553 
554  //perturbativity of the quartic Higgs couplings
555  for(int i=0;i<12;i++)
556  {
557  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
558  }
559 
560  double la1Q = InitialValues[0];
561  double nu1Q = InitialValues[1];
562  double nu2Q = InitialValues[2];
563  double nu3Q = InitialValues[3];
564  double nu4Q = InitialValues[4];
565  double nu5Q = InitialValues[5];
566  double mu1Q = InitialValues[6];
567  double mu2Q = InitialValues[7];
568  double mu3Q = InitialValues[8];
569  double mu4Q = InitialValues[9];
570  double mu5Q = InitialValues[10];
571  double mu6Q = InitialValues[11];
572 
573  //positivity checks
574  double muAtimes2=mu1Q+mu2Q+mu6Q+2.0*(mu3Q+mu4Q+mu5Q);
575  if(la1Q<0.0) check=1;
576  if(muAtimes2<0.0) check=1;
577  if(mu1Q+mu2Q+mu3Q+mu4Q<0.0) check=1;
578  if(14.0*(mu1Q+mu2Q)+5.0*mu6Q+24.0*(mu3Q+mu4Q)-3.0*fabs(2.0*(mu1Q+mu2Q)-mu6Q)<0.0) check=1;
579  if(5.0*(mu1Q+mu2Q+mu6Q)+6.0*(2.0*mu3Q+mu4Q+mu5Q)-fabs(mu1Q+mu2Q+mu6Q)<0.0) check=1;
580  if(sqrt(la1Q*muAtimes2)+nu1Q<0.0) check=1;
581  if(sqrt(la1Q*muAtimes2)+nu1Q+nu2Q-2.0*fabs(nu3Q)<0.0) check=1;
582  if(la1Q+0.25*muAtimes2+nu1Q+nu2Q+2.0*nu3Q-fabs(nu4Q+nu5Q)/sqrt(3.0)<0.0) check=1;
583 
585  //NLO unitarity//
587 
588  double pi=M_PI;
589  gslpp::vector<gslpp::complex> unitarityeigenvalues(8,0.);
590 
591  if(t1>tNLOuni)
592  {
593  //LO part
594  double muA = 4.0*mu1Q+4.0*mu2Q+8.5*mu3Q+5.0*mu4Q+1.5*mu5Q+2.5*mu6Q;
595  double muB = (4.0*mu1Q+4.0*mu2Q+1.5*mu3Q+12.0*mu4Q+1.5*mu5Q-0.5*mu6Q)/3.0;
596  double muC = (-0.5*mu1Q-0.5*mu2Q+1.5*mu3Q+1.5*mu4Q+12.0*mu5Q+4.0*mu6Q)/3.0;
597  double MA1 = 3.0*la1Q + muA - sqrt(9.0*la1Q*la1Q-6.0*la1Q*muA+muA*muA+32.0*nu1Q*nu1Q+32.0*nu1Q*nu2Q+8.0*nu2Q*nu2Q);
598  double MA2 = 3.0*la1Q + muA + sqrt(9.0*la1Q*la1Q-6.0*la1Q*muA+muA*muA+32.0*nu1Q*nu1Q+32.0*nu1Q*nu2Q+8.0*nu2Q*nu2Q);
599  double MB1 = la1Q + muB - sqrt(la1Q*la1Q-2.0*la1Q*muB+muB*muB+8.0*nu2Q*nu2Q);
600  double MB2 = la1Q + muB + sqrt(la1Q*la1Q-2.0*la1Q*muB+muB*muB+8.0*nu2Q*nu2Q);
601  double MC1 = la1Q + muC - sqrt(la1Q*la1Q-2.0*la1Q*muC+muC*muC+32.0*nu3Q*nu3Q);
602  double MC2 = la1Q + muC + sqrt(la1Q*la1Q-2.0*la1Q*muC+muC*muC+32.0*nu3Q*nu3Q);
603  unitarityeigenvalues.assign(0, MA1/(32.0*pi));
604  unitarityeigenvalues.assign(1, MA2/(32.0*pi));
605  unitarityeigenvalues.assign(2, MB1/(32.0*pi));
606  unitarityeigenvalues.assign(3, MB2/(32.0*pi));
607  unitarityeigenvalues.assign(4, MC1/(32.0*pi));
608  unitarityeigenvalues.assign(5, MC2/(32.0*pi));
609  unitarityeigenvalues.assign(6, la1Q/(16.0*pi));
610  unitarityeigenvalues.assign(7, sqrt(15.0)*(nu4Q+nu5Q)/(64.0*pi));
611 
612  //NLO part
613  //beta_la1*16pi^2
614  double betala1 = 12.0*la1Q*la1Q + 8.0*nu1Q*nu1Q + 8.0*nu1Q*nu2Q + 4.0*nu2Q*nu2Q + 16.0*nu3Q*nu3Q;
615  //beta_nu1*16pi^2
616  double betanu1 = 2.0*nu1Q*nu1Q + nu2Q*nu2Q + 4.0*nu3Q*nu3Q + 2.0*la1Q*(3.0*nu1Q+nu2Q)
617  + (7.0*nu4Q*nu4Q - 4.0*nu4Q*nu5Q + 7.0*nu5Q*nu5Q)/3.0
618  + nu1Q*(8.0*mu1Q + 8.0*mu2Q + 17.0*mu3Q + 10.0*mu4Q + 3.0*mu5Q + 5.0*mu6Q)
619  + nu2Q*(8.0*mu1Q + 8.0*mu2Q + 24.0*mu3Q + 3.0*mu4Q
620  + 3.0*mu5Q + 8.0*mu6Q)/3.0;
621  //beta_nu2*16pi^2
622  double betanu2 = 2.0*nu2Q*nu2Q + 4.0*nu1Q*nu2Q + 16.0*nu3Q*nu3Q + 2.0*la1Q*nu2Q
623  + (4.0*nu4Q*nu4Q + 17.0*nu4Q*nu5Q + 4.0*nu5Q*nu5Q)/3.0
624  + nu2Q*(8.0*mu1Q + 8.0*mu2Q + 3.0*mu3Q + 24.0*mu4Q
625  + 3.0*mu5Q - mu6Q)/3.0;
626  //beta_nu3*16pi^2
627  double betanu3 = 2.0*nu3Q*(la1Q + 2.0*nu1Q + 3.0*nu2Q)
628  + (17.0*nu4Q*nu4Q + 16.0*nu4Q*nu5Q + 17.0*nu5Q*nu5Q)/12.0
629  + nu3Q*(-mu1Q - mu2Q + 3.0*mu3Q + 3.0*mu4Q
630  + 24.0*mu5Q + 8.0*mu6Q)/3.0;
631  //beta_nu4*16pi^2
632  double betanu4 = 8.0*nu3Q*nu4Q + 2.0*nu3Q*nu5Q
633  + nu5Q*(2.0*nu2Q - mu2Q + 2.0*mu4Q + 4.0*mu5Q + mu6Q)
634  + nu4Q*(3.0*nu1Q + 2.0*nu2Q + 6.0*mu1Q + 2.0*mu2Q + 3.0*mu3Q
635  + 2.0*mu4Q + mu5Q + mu6Q);
636  //beta_nu5*16pi^2
637  double betanu5 = 2.0*nu3Q*nu4Q + 8.0*nu3Q*nu5Q
638  + nu4Q*(2.0*nu2Q - mu1Q + 2.0*mu4Q + 4.0*mu5Q + mu6Q)
639  + nu5Q*(3.0*nu1Q + 2.0*nu2Q + 6.0*mu1Q + 2.0*mu2Q + 3.0*mu3Q
640  + 2.0*mu4Q + mu5Q + mu6Q);
641  //beta_mu1*16pi^2
642  double betamu1 = 3.0*nu4Q*nu4Q + 7.0*mu1Q*mu1Q
643  + mu1Q*(6.0*mu2Q + 6.0*mu3Q + 4.0*mu4Q - mu5Q - 2.0*mu6Q)
644  + mu2Q*(4.0*mu4Q - mu5Q)
645  - 2.0*mu4Q*mu6Q + 2.0*mu5Q*mu6Q + mu6Q*mu6Q;
646  //beta_mu2*16pi^2
647  double betamu2 = 3.0*nu5Q*nu5Q + 7.0*mu2Q*mu2Q
648  + mu2Q*(6.0*mu1Q + 6.0*mu3Q + 4.0*mu4Q - mu5Q - 2.0*mu6Q)
649  + mu1Q*(4.0*mu4Q - mu5Q)
650  - 2.0*mu4Q*mu6Q + 2.0*mu5Q*mu6Q + mu6Q*mu6Q;
651  //beta_mu3*16pi^2
652  double betamu3 = 20.0*mu3Q*mu3Q
653  + mu3Q*(288.0*mu1Q + 288.0*mu2Q + 360.0*mu4Q + 108.0*mu5Q + 180.0*mu6Q)/18.0
654  + (36.0*nu1Q*nu1Q + 36.0*nu1Q*nu2Q - 24.0*nu4Q*nu4Q - 12.0*nu4Q*nu5Q
655  - 24.0*nu5Q*nu5Q + 62.0*mu1Q*mu1Q + 64.0*mu1Q*mu2Q + 62.0*mu2Q*mu2Q
656  + (96.0*mu4Q + 18.0*mu5Q + 58.0*mu6Q)*(mu1Q + mu2Q)
657  + 54.0*mu4Q*mu4Q + 36.0*mu4Q*mu5Q + 132.0*mu4Q*mu6Q + 18.0*mu5Q*mu5Q
658  + 18.0*mu5Q*mu6Q + 29.0*mu6Q*mu6Q)/18.0;
659  //beta_mu4*16pi^2
660  double betamu4 = nu2Q*nu2Q - (nu4Q*nu4Q - 4.0*nu4Q*nu5Q + nu5Q*nu5Q)/3.0 + 10.0*mu4Q*mu4Q /*mu4Q??*/
661  + mu5Q*(mu1Q + mu2Q + mu6Q)
662  + mu4Q*(4.0*(4.0*mu1Q + 4.0*mu2Q + mu6Q)/3.0 + 2.0*mu5Q + 6.0*mu4Q)
663  + 4.0*mu5Q*mu5Q
664  + (mu1Q*mu1Q + mu2Q*mu2Q - 4.0*mu6Q*(mu1Q+mu2Q) - 2.0*mu6Q*mu6Q)/9.0
665  + 26.0/9.0*mu1Q*mu2Q;
666  //beta_mu5*16pi^2
667  double betamu5 = 4.0*nu3Q*nu3Q - (nu4Q*nu4Q - 4.0*nu4Q*nu5Q + nu5Q*nu5Q)/3.0
668  + mu5Q*((mu1Q + mu2Q + 19.0*mu6Q)/3.0 + 8.0*mu4Q + 6.0*mu3Q)
669  + 2.0*mu4Q*mu6Q + 8.0*mu5Q*mu5Q
670  + (mu1Q*mu1Q + mu2Q*mu2Q - 4.0*mu6Q*(mu1Q+mu2Q) + 7.0*mu6Q*mu6Q)/9.0
671  - 10.0/9.0*mu1Q*mu2Q;
672  //beta_mu6*16pi^2
673  double betamu6 = 0.5*mu6Q*mu6Q + 3.0*nu4Q*nu4Q + 3.0*nu5Q*nu5Q
674  - 2.0*(mu1Q*mu1Q + mu2Q*mu2Q) + 6.0*mu5Q*(mu1Q + mu2Q)
675  + 7.0*mu6Q*(mu1Q + mu2Q + mu3Q);
676 
677  double betamuA = 4.0*betamu1+4.0*betamu2+8.5*betamu3+5.0*betamu4+1.5*betamu5+2.5*betamu6;
678  double betamuB = (4.0*betamu1+4.0*betamu2+1.5*betamu3+12.0*betamu4+1.5*betamu5-0.5*betamu6)/3.0;
679  double betamuC = (-0.5*betamu1-0.5*betamu2+1.5*betamu3+1.5*betamu4+12.0*betamu5+4.0*betamu6)/3.0;
680  double betaMA1 = 3.0*betala1 + betamuA
681  - sqrt(9.0*betala1*betala1-6.0*betala1*betamuA+betamuA*betamuA
682  +32.0*betanu1*betanu1+32.0*betanu1*betanu2+8.0*betanu2*betanu2);
683  double betaMA2 = 3.0*betala1 + betamuA
684  + sqrt(9.0*betala1*betala1-6.0*betala1*betamuA+betamuA*betamuA
685  +32.0*betanu1*betanu1+32.0*betanu1*betanu2+8.0*betanu2*betanu2);
686  double betaMB1 = betala1 + betamuB - sqrt(betala1*betala1-2.0*betala1*betamuB+betamuB*betamuB+8.0*betanu2*betanu2);
687  double betaMB2 = betala1 + betamuB + sqrt(betala1*betala1-2.0*betala1*betamuB+betamuB*betamuB+8.0*betanu2*betanu2);
688  double betaMC1 = betala1 + betamuC - sqrt(betala1*betala1-2.0*betala1*betamuC+betamuC*betamuC+32.0*betanu3*betanu3);
689  double betaMC2 = betala1 + betamuC + sqrt(betala1*betala1-2.0*betala1*betamuC+betamuC*betamuC+32.0*betanu3*betanu3);
690 
691  gslpp::vector<gslpp::complex> betaeigenvalues(8,0.);
692  gslpp::vector<gslpp::complex> NLOunitarityeigenvalues(8,0.);
693 
694  betaeigenvalues.assign(0, -1.5 * betaMA1/(32.0*pi));
695  betaeigenvalues.assign(1, -1.5 * betaMA2/(32.0*pi));
696  betaeigenvalues.assign(2, -1.5 * betaMB1/(32.0*pi));
697  betaeigenvalues.assign(3, -1.5 * betaMB2/(32.0*pi));
698  betaeigenvalues.assign(4, -1.5 * betaMC1/(32.0*pi));
699  betaeigenvalues.assign(5, -1.5 * betaMC2/(32.0*pi));
700  betaeigenvalues.assign(6, -1.5 * betala1/(16.0*pi));
701  betaeigenvalues.assign(7, -1.5 * sqrt(15.0)*(betanu4+betanu5)/(64.0*pi));
702 
703  for (int i=0; i < 8; i++) {
704  NLOunitarityeigenvalues.assign(i, -(gslpp::complex::i()-1.0/pi)*unitarityeigenvalues(i)*unitarityeigenvalues(i) + betaeigenvalues(i) );
705  if( ( unitarityeigenvalues(i) + NLOunitarityeigenvalues(i).real() ).abs() > 0.5) check=1;
706  if( (unitarityeigenvalues(i)).abs() > Rpeps && (NLOunitarityeigenvalues(i)/unitarityeigenvalues(i)).abs() > 1.0) check=1;
707  }
708 
709  } //end of the if(t1>tNLOuni)
710 
711 
712  return check;
713 }
714 
715 int RGEcheckcustodialMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
716 {
717  int check=0;
718 
719  //perturbativity of the quartic Higgs couplings
720  for(int i=0;i<7;i++)
721  {
722  if(fabs(InitialValues[i])>4.0*M_PI) check=1;
723  }
724 
725  double la1Q = InitialValues[0];
726  double nu1Q = InitialValues[1];
727  double nu2Q = InitialValues[2];
728  double nu4Q = InitialValues[3];
729  double mu1Q = InitialValues[4];
730  double mu3Q = InitialValues[5];
731  double mu4Q = InitialValues[6];
732 
733  //positivity checks
734  double muAtimes2=4.0*mu1Q+2.0*mu3Q+4.0*mu4Q;
735  if(la1Q<0.0) check=1;
736  if(muAtimes2<0.0) check=1;
737  if(5.0*mu1Q+3.0*mu3Q+3.0*mu4Q-fabs(mu1Q)<0.0) check=1;
738  if(sqrt(la1Q*muAtimes2)+nu1Q+nu2Q-fabs(nu2Q)<0.0) check=1;
739  if(la1Q+0.25*muAtimes2+nu1Q+2.0*nu2Q-2.0*fabs(nu4Q)/sqrt(3.0)<0.0) check=1;
740 
741  //unitarity checks
742 
743  double pi=M_PI;
744  gslpp::matrix<gslpp::complex> Smatrix1(2,2,0.), Smatrix2(2,2,0.);
745  gslpp::matrix<gslpp::complex> Seigenvectors1(2,2,0.), Seigenvectors2(2,2,0.);
746  gslpp::vector<double> Seigenvalues1(2,0.), Seigenvalues2(2,0.);
747  gslpp::vector<gslpp::complex> unitarityeigenvalues(5,0.);
748 
749  if(t1>tNLOuni)
750  {
751 
752  //LO part
753  Smatrix1.assign(0,0, 3.0*la1Q/(16.0*pi));
754  Smatrix1.assign(0,1, (2.0*nu1Q+nu2Q)/(8.0*sqrt(2.0)*pi));
755  Smatrix1.assign(1,0, Smatrix1(0,1));
756  Smatrix1.assign(1,1, (26.0*mu1Q+17.0*mu3Q+13.0*mu4Q)/(32.0*pi));
757 
758  Smatrix2.assign(0,0, la1Q/(16.0*pi));
759  Smatrix2.assign(0,1, nu2Q/(8.0*sqrt(2.0)*pi));
760  Smatrix2.assign(1,0, Smatrix2(0,1));
761  Smatrix2.assign(1,1, (14.0*mu1Q+3.0*mu3Q+27.0*mu4Q)/(96.0*pi));
762 
763  Smatrix1.eigensystem(Seigenvectors1, Seigenvalues1);
764  Smatrix2.eigensystem(Seigenvectors2, Seigenvalues2);
765 
766  for (int i=0; i < 2; i++) {
767  unitarityeigenvalues.assign(i, Seigenvalues1(i));
768  unitarityeigenvalues.assign(2+i, Seigenvalues2(i));
769  }
770  unitarityeigenvalues.assign(4, sqrt(15.0)*nu4Q/(16.0*pi));
771 
772  //beta_la1*16pi^2
773  double betala1 = 12.0*la1Q*la1Q + 8.0*nu1Q*nu1Q + 8.0*nu1Q*nu2Q + 8.0*nu2Q*nu2Q;
774  //beta_mu1*16pi^2
775  double betamu1 = 13.0*mu1Q*mu1Q + 6.0*mu1Q*mu3Q + 6.0*mu1Q*mu4Q + 3.0*nu4Q*nu4Q;
776  //beta_mu3*16pi^2
777  double betamu3 = (134.0*mu1Q*mu1Q + 6.0*mu1Q*(39.0*mu3Q + 22.0*mu4Q)
778  + 3.0*(30.0*mu3Q*mu3Q + 39.0*mu3Q*mu4Q + 9.0*mu4Q*mu4Q
779  + 3.0*nu1Q*nu1Q + 3.0*nu1Q*nu2Q - 5.0*nu4Q*nu4Q))/4.5;
780  //beta_mu4*16pi^2
781  double betamu4 = (4.0*mu1Q*mu1Q + 156.0*mu1Q*mu4Q + 54.0*mu3Q*mu4Q + 144.0*mu4Q*mu4Q
782  + 9.0*nu2Q*nu2Q + 6.0*nu4Q*nu4Q)/9.0;
783  //beta_nu1*16pi^2
784  double betanu1 = (18.0*la1Q*nu1Q
785  + 78.0*mu1Q*nu1Q + 51.0*mu3Q*nu1Q + 39.0*mu4Q*nu1Q + 6.0*nu1Q*nu1Q
786  + 6.0*la1Q*nu2Q + 32.0*mu1Q*nu2Q + 24.0*mu3Q*nu2Q + 6.0*mu4Q*nu2Q
787  + 6.0*nu2Q*nu2Q + 10.0*nu4Q*nu4Q)/3.0;
788  //beta_nu2*16pi^2
789  double betanu2 = 2.0*la1Q*nu2Q + ((14.0*mu1Q)/3.0 + mu3Q + 9.0*mu4Q)*nu2Q
790  + 4.0*nu1Q*nu2Q + 6.0*nu2Q*nu2Q + (25.0*nu4Q*nu4Q)/3.0;
791  //beta_nu4*16pi^2
792  double betanu4 = 11.0*mu1Q*nu4Q + 3.0*mu3Q*nu4Q + 9.0*mu4Q*nu4Q + 3.0*nu1Q*nu4Q + 9.0*nu2Q*nu4Q;
793 
794 // diagonalization
795  gslpp::matrix<gslpp::complex> Sbmatrix1(2,2,0.), Sbmatrix2(2,2,0.);
796  gslpp::matrix<gslpp::complex> Seigenvectors1T(2,2,0.), Seigenvectors2T(2,2,0.);
797  gslpp::vector<gslpp::complex> Sbeigenvalues1(2,0.), Sbeigenvalues2(2,0.);
798  gslpp::vector<gslpp::complex> betaeigenvalues(5,0.);
799  gslpp::vector<gslpp::complex> NLOunitarityeigenvalues(5,0.);
800 
801  Sbmatrix1.assign(0,0, 3.0*betala1/(16.0*pi));
802  Sbmatrix1.assign(0,1, (2.0*betanu1+betanu2)/(8.0*sqrt(2.0)*pi));
803  Sbmatrix1.assign(1,0, Sbmatrix1(0,1));
804  Sbmatrix1.assign(1,1, (26.0*betamu1+17.0*betamu3+13.0*betamu4)/(32.0*pi));
805 
806  Sbmatrix2.assign(0,0, betala1/(16.0*pi));
807  Sbmatrix2.assign(0,1, betanu2/(8.0*sqrt(2.0)*pi));
808  Sbmatrix2.assign(1,0, Sbmatrix2(0,1));
809  Sbmatrix2.assign(1,1, (14.0*betamu1+3.0*betamu3+27.0*betamu4)/(96.0*pi));
810 
811  Seigenvectors1T=Seigenvectors1.hconjugate();
812  Seigenvectors2T=Seigenvectors2.hconjugate();
813 
814  for (int i=0; i < 2; i++) {
815  for (int k=0; k < 2; k++) {
816  for (int l=0; l < 2; l++) {
817  Sbeigenvalues1.assign(i, Sbeigenvalues1(i) + Seigenvectors1T(i,k) * Sbmatrix1(k,l) * Seigenvectors1(l,i) );
818  Sbeigenvalues2.assign(i, Sbeigenvalues2(i) + Seigenvectors2T(i,k) * Sbmatrix2(k,l) * Seigenvectors2(l,i) );
819  }
820  }
821  betaeigenvalues.assign(i, -1.5 * Sbeigenvalues1(i));
822  betaeigenvalues.assign(i+2, -1.5 * Sbeigenvalues2(i));
823  }
824 
825  betaeigenvalues.assign(4, -1.5 * sqrt(15.0)*betanu4/(16.0*pi));
826 
827  for (int i=0; i < 5; i++) {
828  NLOunitarityeigenvalues.assign(i, -(gslpp::complex::i()-1.0/pi)*unitarityeigenvalues(i)*unitarityeigenvalues(i) + betaeigenvalues(i) );
829  if( ( unitarityeigenvalues(i) + NLOunitarityeigenvalues(i).real() ).abs() > 0.5) check=1;
830  if( (unitarityeigenvalues(i)).abs() > Rpeps && (NLOunitarityeigenvalues(i)/unitarityeigenvalues(i)).abs() > 1.0) check=1;
831  }
832 
833  } //end of the if(t1>tNLOuni)
834  return check;
835 }
836 
837 double RunnerTHDMW::RGERunnerTHDMW(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double NLOuniscale)
838 {
839  //Define which stepping function should be used
840  const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk4;
841 
842  //Allocate space for the stepping function
843  gsl_odeiv2_step * s = gsl_odeiv2_step_alloc(T, NumberOfRGEs);
844 
845  //Define the absolute (A) and relative (R) error on y at each step.
846  //The real error will be compared to the following error estimate:
847  // A + R * |y_i|
848  gsl_odeiv2_control * c = gsl_odeiv2_control_y_new(1e-6, 0.0);
849 
850  //Allocate space for the evolutor
851  gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(NumberOfRGEs);
852 
853  //Definition of the RGE system (the Jacobian is not necessary for the RK4 method; it's an empty function here)
854  gsl_odeiv2_system RGEsystem = {RGEsTHDMW, JacobianTHDMW, NumberOfRGEs, &order};
855 
856  //Set starting and end point as natural logarithmic scales (conversion from decadic log scale)
857  double t1 = Q1*log(10.0);
858  double t2 = Q2*log(10.0);
859  double tNLOuni = NLOuniscale*log(10.0);
860 
861  //Set initial step size
862  double InitialStepSize = 1e-6;
863 
864  //Run!
865  while (t1 < t2)
866  {
867  int status = gsl_odeiv2_evolve_apply (e, c, s, &RGEsystem, &t1, t2, &InitialStepSize, InitialValues);
868  if(status != GSL_SUCCESS) break;
869 
870  //intermediate checks if appropriate
871  if(RGEcheckTHDMW(InitialValues,t1,Rpeps,tNLOuni) != 0) break;
872  }
873 
874  gsl_odeiv2_evolve_free (e);
875  gsl_odeiv2_control_free (c);
876  gsl_odeiv2_step_free (s);
877 
878  //Return the decadic log scale at which the evolution stopped
879  return t1/log(10.0);
880 }
881 
882 double RunnerTHDMW::RGERunnerMW(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double NLOuniscale)
883 {
884  //Define which stepping function should be used
885  const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk4;
886 
887  //Allocate space for the stepping function
888  gsl_odeiv2_step * s = gsl_odeiv2_step_alloc(T, NumberOfRGEs);
889 
890  //Define the absolute (A) and relative (R) error on y at each step.
891  //The real error will be compared to the following error estimate:
892  // A + R * |y_i|
893  gsl_odeiv2_control * c = gsl_odeiv2_control_y_new(1e-6, 0.0);
894 
895  //Allocate space for the evolutor
896  gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(NumberOfRGEs);
897 
898  //Definition of the RGE system (the Jacobian is not necessary for the RK4 method; it's an empty function here)
899  gsl_odeiv2_system RGEsystem = {RGEsMW, JacobianTHDMW, NumberOfRGEs, &order};
900 
901  //Set starting and end point as natural logarithmic scales (conversion from decadic log scale)
902  double t1 = Q1*log(10.0);
903  double t2 = Q2*log(10.0);
904  double tNLOuni = NLOuniscale*log(10.0);
905 
906  //Set initial step size
907  double InitialStepSize = 1e-6;
908 
909  //Run!
910  while (t1 < t2)
911  {
912  int status = gsl_odeiv2_evolve_apply (e, c, s, &RGEsystem, &t1, t2, &InitialStepSize, InitialValues);
913  if(status != GSL_SUCCESS) break;
914 
915  //intermediate checks if appropriate
916  if(RGEcheckMW(InitialValues,t1,Rpeps,tNLOuni) != 0) break;
917  }
918 
919  gsl_odeiv2_evolve_free (e);
920  gsl_odeiv2_control_free (c);
921  gsl_odeiv2_step_free (s);
922 
923  //Return the decadic log scale at which the evolution stopped
924  return t1/log(10.0);
925 }
926 
927 double RunnerTHDMW::RGERunnercustodialMW(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double NLOuniscale)
928 {
929  //Define which stepping function should be used
930  const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk4;
931 
932  //Allocate space for the stepping function
933  gsl_odeiv2_step * s = gsl_odeiv2_step_alloc(T, NumberOfRGEs);
934 
935  //Define the absolute (A) and relative (R) error on y at each step.
936  //The real error will be compared to the following error estimate:
937  // A + R * |y_i|
938  gsl_odeiv2_control * c = gsl_odeiv2_control_y_new(1e-6, 0.0);
939 
940  //Allocate space for the evolutor
941  gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(NumberOfRGEs);
942 
943  //Definition of the RGE system (the Jacobian is not necessary for the RK4 method; it's an empty function here)
944  gsl_odeiv2_system RGEsystem = {RGEscustodialMW, JacobianTHDMW, NumberOfRGEs, &order};
945 
946  //Set starting and end point as natural logarithmic scales (conversion from decadic log scale)
947  double t1 = Q1*log(10.0);
948  double t2 = Q2*log(10.0);
949  double tNLOuni = NLOuniscale*log(10.0);
950 
951  //Set initial step size
952  double InitialStepSize = 1e-6;
953 
954  //Run!
955  while (t1 < t2)
956  {
957  int status = gsl_odeiv2_evolve_apply (e, c, s, &RGEsystem, &t1, t2, &InitialStepSize, InitialValues);
958  if(status != GSL_SUCCESS) break;
959 
960  //intermediate checks if appropriate
961  if(RGEcheckcustodialMW(InitialValues,t1,Rpeps,tNLOuni) != 0) break;
962  }
963 
964  gsl_odeiv2_evolve_free (e);
965  gsl_odeiv2_control_free (c);
966  gsl_odeiv2_step_free (s);
967 
968  //Return the decadic log scale at which the evolution stopped
969  return t1/log(10.0);
970 }
RunnerTHDMW.h
RGEcheckMW
int RGEcheckMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: RunnerTHDMW.cpp:550
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
gslpp::matrix< gslpp::complex >
RunnerTHDMW::RGERunnerMW
virtual double RGERunnerMW(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double tNLOuni)
Definition: RunnerTHDMW.cpp:882
RunnerTHDMW::RunnerTHDMW
RunnerTHDMW(const StandardModel &SM_i)
RunnerTHDMW constructor.
Definition: RunnerTHDMW.cpp:13
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
RGEcheckcustodialMW
int RGEcheckcustodialMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: RunnerTHDMW.cpp:715
THDMW
A base class for symmetric Two-Higgs-Doublet-Manohar-Wise models.
Definition: THDMW.h:233
RunnerTHDMW::~RunnerTHDMW
virtual ~RunnerTHDMW()
Runner destructor.
Definition: RunnerTHDMW.cpp:16
RunnerTHDMW::RGERunnercustodialMW
virtual double RGERunnercustodialMW(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double tNLOuni)
Definition: RunnerTHDMW.cpp:927
RunnerTHDMW::RGERunnerTHDMW
virtual double RGERunnerTHDMW(double InitialValues[], unsigned long int NumberOfRGEs, double Q1, double Q2, int order, double Rpeps, double tNLOuni)
Definition: RunnerTHDMW.cpp:837
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
RGEsMW
int RGEsMW(double t, const double y[], double beta[], void *flags)
Definition: RunnerTHDMW.cpp:139
RGEcheckTHDMW
int RGEcheckTHDMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: RunnerTHDMW.cpp:334
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
RGEsTHDMW
int RGEsTHDMW(double t, const double y[], double beta[], void *flags)
Definition: RunnerTHDMW.cpp:19
JacobianTHDMW
int JacobianTHDMW(double t, const double y[], double *dfdy, double dfdt[], void *order)
Definition: RunnerTHDMW.cpp:329
RGEscustodialMW
int RGEscustodialMW(double t, const double y[], double beta[], void *flags)
Definition: RunnerTHDMW.cpp:260
gslpp::vector< gslpp::complex >