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
13RunnerTHDMW::RunnerTHDMW(const StandardModel& SM_i) : myTHDMW(static_cast<const THDMW*> (&SM_i))
14{}
15
16RunnerTHDMW::~RunnerTHDMW()
17{};
18
19int 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
139int 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
260int 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
329int JacobianTHDMW (double t, const double y[], double *dfdy, double dfdt[], void *order)
330{
331 return 0;
332}
333
334int 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
550int 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
715int 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
837double 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
882double 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
927double 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}
RGEcheckTHDMW
int RGEcheckTHDMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: RunnerTHDMW.cpp:334
RGEsTHDMW
int RGEsTHDMW(double t, const double y[], double beta[], void *flags)
Definition: RunnerTHDMW.cpp:19
RGEcheckcustodialMW
int RGEcheckcustodialMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: RunnerTHDMW.cpp:715
RGEscustodialMW
int RGEscustodialMW(double t, const double y[], double beta[], void *flags)
Definition: RunnerTHDMW.cpp:260
RGEcheckMW
int RGEcheckMW(const double InitialValues[], const double t1, const double Rpeps, const double tNLOuni)
Definition: RunnerTHDMW.cpp:550
RGEsMW
int RGEsMW(double t, const double y[], double beta[], void *flags)
Definition: RunnerTHDMW.cpp:139
JacobianTHDMW
int JacobianTHDMW(double t, const double y[], double *dfdy, double dfdt[], void *order)
Definition: RunnerTHDMW.cpp:329
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
RunnerTHDMW::RunnerTHDMW
RunnerTHDMW(const StandardModel &SM_i)
RunnerTHDMW constructor.
Definition: RunnerTHDMW.cpp:13
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
virtual ~RunnerTHDMW()
Runner destructor.
Definition: RunnerTHDMW.cpp:16
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:509
THDMW
A base class for symmetric Two-Higgs-Doublet-Manohar-Wise models.
Definition: THDMW.h:233
s
Test Observable.
t
Test Observable.