a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Charm_Kpnunu.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012 HEPfit Collaboration
3 *
4 *
5 * For the licensing terms see doc/COPYING.
6 */
7
8#include "Charm_Kpnunu.h"
9#include "StandardModel.h"
10
12: model(model_i),evoCkpnn(1, NDR, NNLO, NLO_QED11), dcp(3, 0.) ,dcb(2, 0.),
13 CW0b(2,0.),CW1b(2,0.),CW2b(2,0.),CWeb(2,0.),CWesb(2,0.),
14 U0_4b(2,0.),U0_5b(2,0.),J1_4b(2,0.),J2_4b(2,0.),J1_5b(2,0.),J2_5b(2,0.),
15 R0_4b(2,0.),R0_5b(2,0.),R1_4b(2,0.),R1_5b(2,0.),
16 CW0p(3,0.),CW1p(3,0.),CW2p(3,0.),CWep(3,0.),CWesp(3,0.),
17 U0_4p(3,0.),U0_5p(3,0.),J1_4p(3,0.),J2_4p(3,0.),J1_5p(3,0.),J2_5p(3,0.),
18 R0_4p(3,0.),R0_5p(3,0.),R1_4p(3,0.),R1_5p(3,0.)
19{
25 kc= pow(etac, 24. / 25.);
26 xc_mc_qed=sqrt(2.) * model.sW2_ND() * model.getGF() / M_PI / model.alphaMz() * mc_mc * mc_mc ;
27 L = log(model.getMuc() * model.getMuc() / mc_mc / mc_mc);
28
29 xi1c=15212. / 1875. * (etac - 1.) / etac;
30 xi2c=966966391. / 10546875. - 231404944. / 3515625. * (1. / etac) - 272751559. / 10546875. *
31 (1. / etac)*(1. / etac) - 128. / 5. * (1. - (1. / etac)*(1. / etac)) * gsl_sf_zeta_int(3);;
32 xice= 8. / 3. / (11. - 2. / 3. * 3. ) * (etac - 1.);
33 xices=(32. / 9. / (11. - 2. / 3. * 3. ) - (-8. / 9. * (1. + 2./4. ) ) * 8. / (11. - 2. / 3. * 3. ) / (11. - 2. / 3. * 3. ) - (102. - 38. / 3. * 3. ) * 8. / 3. / (11. - 2. / 3. * 3. ) / (11. - 2. / 3. * 3. )) *
34 log(etac) + 8. / 3. / (11. - 2. / 3. * 3. ) * (8. * (102. - 38. / 3. * 3. ) / (11. - 2. / 3. * 3. ) / (11. - 2. / 3. * 3. ) - (404. / 3. - 40. / 9. * 3.) / (11. - 2. / 3. * 3. )) * (1. - 1. / etac) * (1. - etac);
35
36 //box
37 CW0b = CWin_muw(LO,NO_QED ,0);
38 CW1b = CWin_muw(NLO,NO_QED ,0);
40 CWeb = CWin_muw(LO,LO_QED ,0);
42
43 U0_4b = RGevol_J(LO,4,0);
44 U0_5b = RGevol_J(LO,5,0);
45 J1_4b = RGevol_J(NLO,4,0);
46 J2_4b = RGevol_J(NNLO,4,0);
47 J1_5b = RGevol_J(NLO,5,0);
48 J2_5b = RGevol_J(NNLO,5,0);
49
50 R0_4b = RGevol_R(LO_QED,4,0);
51 R0_5b = RGevol_R(LO_QED,5,0);
54
55 //penguin
56 CW0p = CWin_muw(LO,NO_QED ,1);
57 CW1p = CWin_muw(NLO,NO_QED ,1);
59 CWep = CWin_muw(LO,LO_QED ,1);
61
62 U0_4p = RGevol_J(LO,4,1);
63 U0_5p = RGevol_J(LO,5,1);
64 J1_4p = RGevol_J(NLO,4,1);
65 J2_4p = RGevol_J(NNLO,4,1);
66 J1_5p = RGevol_J(NLO,5,1);
67 J2_5p = RGevol_J(NNLO,5,1);
68
69 R0_4p = RGevol_R(LO_QED,4,1);
70 R0_5p = RGevol_R(LO_QED,5,1);
72 R1_5p = RGevol_R(NLO_QED11,5,1);
73}
74
76{}
77
78gslpp::matrix<double> Charm_Kpnunu::RGevolP(int nf, orders order)
79{
80
81 gslpp::matrix<double> evo(3, 3, 0.);
82
83 switch (order) {
84 case (LO):
85 switch (nf) {
86 case(5): //U5P(0)
87 evo(0, 0) = pow(etab, 6. / 23.);
88 evo(1, 1) = pow(etab, -12. / 23.);
89 evo(2, 0) = 12. / 5. * (pow(etab, 6. / 23.) - pow(etab, 1. / 23.));
90 evo(2, 1) = 6. / 13. * (pow(etab, -12. / 23.) - pow(etab, 1. / 23.));
91 evo(2, 2) = pow(etab, 1. / 23.);
92
93 return (evo);
94
95 case(4): //U4P(0)
96 evo(0, 0) = pow(etacb, 6. / 25.);
97 evo(1, 1) = pow(etacb, -12. / 25.);
98 evo(2, 0) = 12. / 7. * (pow(etacb, 6. / 25.) - pow(etacb, -1. / 25.));
99 evo(2, 1) = 6. / 11. * (pow(etacb, -12. / 25.) - pow(etacb, -1. / 25.));
100 evo(2, 2) = pow(etacb, -1. / 25.);
101
102 return (evo);
103
104 default:
105 std::stringstream out;
106 out << nf;
107 throw std::runtime_error("Charm_Kpnunu::RgevolP() " + out.str() + " wrong number of falvours");
108 }
109 case(NLO):
110 switch (nf) {
111 case(5): //J5P(1)
112 evo(0, 0) = 5165. / 3174.;
113 evo(1, 1) = -2267. / 1587.;
114 evo(2, 0) = -15857. / 1587.;
115 evo(2, 1) = 15305. / 3174.;
116 evo(2, 2) = -14924. / 1587.;
117
118 return (evo);
119
120 case(4): //J4P(1)
121 evo(0, 0) = 6719. / 3750.;
122 evo(1, 1) = -3569. / 1875.;
123 evo(2, 0) = -15931. / 1875.;
124 evo(2, 1) = 5427. / 1250.;
125 evo(2, 2) = -15212. / 1875.;
126
127 return (evo);
128
129 default:
130 std::stringstream out;
131 out << nf;
132 throw std::runtime_error("Charm_Kpnunu::RgevolP() " + out.str() + " wrong number of flavours");
133 }
134 case(NNLO):
135 switch (nf) {
136 case(5): //J5P(2)
137 evo(0, 0) = -7.35665;
138 evo(1, 1) = -54.9107;
139 evo(2, 0) = 17.7699;
140 evo(2, 1) = -1.7514;
141 evo(2, 2) = 18.3025;
142
143 return (evo);
144
145 case(4): //J4P(2)
146 evo(0, 0) = -10.2451;
147 evo(1, 1) = -50.3422;
148 evo(2, 0) = 8.0325;
149 evo(2, 1) = -0.3657;
150 evo(2, 2) = 4.91177;
151
152 return (evo);
153
154 default:
155 std::stringstream out;
156 out << nf;
157 throw std::runtime_error("Charm_Kpnunu::RgevolP() " + out.str() + " wrong number of flavours");
158 }
159 default:
160 std::stringstream out;
161 out << order;
162 throw std::runtime_error("Charm_Kpnunu::RgevolP() " + out.str() + " wrong order assignment ");
163 }
164}
165
166gslpp::vector<double> Charm_Kpnunu::ThresholdCp(orders order)
167{
168
169 double mub = model.getMub();
172 double Lb = log(mub * mub / Mb / Mb);
173
174 switch (order) {
175 case(LO):
176 dcp = 0.;
177
178 return (dcp);
179
180 case(NLO):
181 dcp = 0.;
182
183 return (dcp);
184
185 case(NNLO):
186 dcp(0) = -pow(etab, 6. / 23.)*(2. / 3. * Lb * ((631. + 9699.) / 6348. * (1 - etab) + etab * CW1p(0))
187 - 2. / 3. * (59. / 36. + 1. / 3. * Lb + Lb * Lb));
188 dcp(1) = -pow(etab, -12. / 23.)*(2. / 3. * Lb * ((631. - 9699.) / 6348. * (1 - etab) + etab * CW1p(1))
189 + 4. / 3. * (59. / 36. + 1. / 3. * Lb + Lb * Lb));
190 dcp(2) = (-2. / 3.) * Lb * ((284704. / 2645. + 694522. / 20631. * etab) * pow(etab, 1. / 23.)
191 -(1033492. / 7935. + 8264. / 529. * etab) * pow(etab, 6. / 23.) + (3058. / 1587. + 18136. / 6877. * etab)
192 * pow(etab, -12. / 23.) + etab * (pow(etab, 1. / 23.) * CW1p(2) + 48. / 5. * (pow(etab, 6. / 23.)
193 - pow(etab, 1. / 23.)) * CW1p(0) + 24. / 13. * (pow(etab, -12. / 23.) - pow(etab, 1. / 23.)) * CW1p(1)));
194
195 return (dcp);
196
197 default:
198 std::stringstream out;
199 out << order;
200 throw std::runtime_error("Charm_Kpnunu::MatchingCp() " + out.str() + " wrong order assignment ");
201 }
202}
203
205{
206 switch (order) {
207 case(LO):{
208 return (C(LO,NO_QED,1)(2));
209 }
210 case(NLO):{
211 double rhoP1p = 4 * (1. - L) + 4. * log(kc);
212 double rhoP1m = -2. * (1. - L) - 2 * log(kc);
213 gslpp::vector<double> C_LO = C(LO,NO_QED,1);
214
215 return (C(NLO,NO_QED,1)(2) + 4. * (C_LO(0) * rhoP1p +C_LO(1) * rhoP1m) + xi1c * C_LO(2));
216 }
217 case(NNLO):{
218 double rhoP1p = 4 * (1. - L) + 4. * log(kc);
219 double rhoP1m = -2. * (1. - L) - 2 * log(kc);
220 double rhoP2p = 11. + 20. * L - 12. * L * L - 20. * log(kc) - 12. * log(kc) * log(kc) +
221 24. * log(kc) * L + 4. * xi1c;
222 double rhoP2m = -7. + 12. * L + 12. * L * L - 12. * log(kc) + 12. * log(kc) * log(kc)
223 - 24. * log(kc) * L - 2. * xi1c;
224 gslpp::vector<double> C_NLO = C(NLO,NO_QED,1);
225 gslpp::vector<double> C_LO = C(LO,NO_QED,1);
226
227 return (C(NNLO,NO_QED,1)(2) + 4. * (C_NLO(0) * rhoP1p + C_NLO(1) * rhoP1m + C_LO(0) * rhoP2p
228 + C_LO(1) * rhoP2m) + xi1c * (C_NLO(2) + 4. * (C_LO(0) * rhoP1p + C_LO(1)
229 * rhoP1m)) + xi2c * C_LO(2));
230 }
231 default:{
232 std::stringstream out;
233 out << order;
234 throw std::runtime_error("Charm_Kpnunu::C_P() order" + out.str() + " not implemented");
235 }
236 }
237}
238
239gslpp::matrix<double> Charm_Kpnunu::RGevolB(int nf, orders order)
240{
241
242 gslpp::matrix<double> evo(2, 2, 0.);
243
244 switch (order) {
245 case (LO):
246 switch (nf) {
247 case(5):
248 evo(0, 0) = 1.;
249 evo(0, 1) = 0.;
250 evo(1, 0) = 12. * (1. - pow(etab, 1. / 23.));
251 evo(1, 1) = pow(etab, 1. / 23.);
252
253 return (evo);
254
255 case(4):
256 evo(0, 0) = 1.;
257 evo(0, 1) = 0.;
258 evo(1, 0) = -12. * (1. - pow(etacb, -1. / 25.));
259 evo(1, 1) = pow(etacb, -1. / 25.);
260
261 return (evo);
262
263 default:
264 std::stringstream out;
265 out << nf;
266 throw std::runtime_error("Charm_Kpnunu::RgevolB() " + out.str() + " wrong number of falvours");
267 }
268 case(NLO):
269 switch (nf) {
270 case(5):
271 evo(0, 0) = 0.;
272 evo(0, 1) = 0.;
273 evo(1, 0) = 2402. / 1587.;
274 evo(1, 1) = -14924. / 1587.;
275
276 return (evo);
277
278 case(4):
279 evo(0, 0) = 0.;
280 evo(0, 1) = 0.;
281 evo(1, 0) = 581. / 1875.;
282 evo(1, 1) = -15212. / 1875.;
283
284 return (evo);
285
286 default:
287 std::stringstream out;
288 out << nf;
289 throw std::runtime_error("Charm_Kpnunu::RgevolB() " + out.str() + " wrong number of flavours");
290 }
291 case(NNLO):
292 switch (nf) {
293 case(5):
294 evo(0, 0) = 0.;
295 evo(0, 1) = 0.;
296 evo(1, 0) = 1296371522. / 39457581. - 34624. / 1081. * gsl_sf_zeta_int(3);
297 evo(1, 1) = -177621017. / 7555707. + 800. / 23. * gsl_sf_zeta_int(3);
298
299 return (evo);
300
301 case(4):
302 evo(0, 0) = 0.;
303 evo(0, 1) = 0.;
304 evo(1, 0) = 684990354. / 19140625. - 6976. / 245. * gsl_sf_zeta_int(3);
305 evo(1, 1) = -272751559. / 10546875. + 128. / 5. * gsl_sf_zeta_int(3);
306
307 return (evo);
308
309 default:
310 std::stringstream out;
311 out << nf;
312 throw std::runtime_error("Charm_Kpnunu::RgevolB() " + out.str() + " wrong number of flavours");
313 }
314 default:
315 std::stringstream out;
316 out << order;
317 throw std::runtime_error("Charm_Kpnunu::RgevolB() " + out.str() + " wrong order assignment ");
318 }
319}
320
321gslpp::vector<double> Charm_Kpnunu::ThresholdCb(orders order)
322{
323
324 double mub = model.getMub();
327 double Lb = log(mub * mub / Mb / Mb);
328
329 switch (order) {
330 case(LO):
331 dcb = 0.;
332
333 return (dcb);
334
335 case(NLO):
336 dcb = 0.;
337
338 return (dcb);
339
340 case(NNLO):
341 dcb(0) = 0.;
342 dcb(1) = -2. / 3. * Lb * ((238784. / 529. - 9608. / 1587 * etab) * pow(etab, 1. / 23.) - 1336. / 3. +
343 pow(etab, 24. / 23.) * CW1b(1));
344
345 return (dcb);
346
347 default:
348 std::stringstream out;
349 out << order;
350 throw std::runtime_error("Charm_Kpnunu::MatchingCp() " + out.str() + " wrong order assignment ");
351 }
352}
353
355{
356 switch (order) {
357 case(LO):{
358 return (C(LO,NO_QED,0)(1));
359 }
360 case(NLO): {
361 double rhoB1_e = 5. + 4. * L - 4. * log(kc);
362
363 return (C(NLO,NO_QED,0)(1) + 4. * rhoB1_e + xi1c * C(LO,NO_QED,0)(1));
364 }
365 case(NNLO): {
366 double rhoB1_e = 5. + 4. * L - 4. * log(kc);
367 double rhoB2_e = -2. * model.getCF()*(9. - L - 6. * L * L) - 8. / 3. * log(kc) + 16. * log(kc) * log(kc)
368 - 32. * log(kc) * L - 4. * xi1c;
369
370 return (C(NNLO,NO_QED,0)(1) + 4. * rhoB2_e + xi1c * C(NLO,NO_QED,0)(1) + 4. * xi1c * rhoB1_e + xi2c * C(LO,NO_QED,0)(1));
371 }
372 default:{
373 std::stringstream out;
374 out << order;
375 throw std::runtime_error("Charm_Kpnunu::C_Be() order" + out.str() + " not implemented");
376 }
377 }
378}
379
381{
382
383 switch (order) {
384 case(LO):{
385 return ( C(LO,NO_QED,0)(1));
386 }
387 case(NLO):{
389 double rhoB1_t = 5. + 4. * L + 4. * x_t * log(x_t) / (1. - x_t) + 4. / (x_t - kc)*(kc * log(kc) -
390 x_t * (1. - kc) / (1. - x_t) * log(x_t));
391 return (C(NLO,NO_QED,0)(1) + 4. * rhoB1_t + xi1c * C(LO,NO_QED,0)(1));
392 }
393 case(NNLO):{
395 double rhoB1_t = 5. + 4. * L + 4. * x_t * log(x_t) / (1. - x_t) + 4. / (x_t - kc)*(kc * log(kc) -
396 x_t * (1. - kc) / (1. - x_t) * log(x_t));
397
398 double rhoB2_t = -2. * model.getCF()*((9. + 7. * x_t) / (1. - x_t) + x_t * (3. + 13. * x_t) /
399 (1. - x_t) / (1. - x_t) * log(x_t) - 12. * x_t / (1. - x_t) * gsl_sf_dilog(1. - x_t)
400 - ((1. - 13. * x_t) / (1. - x_t) - 12. * x_t * x_t / (1. - x_t) / (1. - x_t) * log(x_t))
401 * L - 6. * L * L) + 32. / (x_t - kc)*(4. * x_t * (1 - kc) / 3. / (1. - x_t) - x_t
402 * ((x_t * (13. - 29. * x_t)) + kc * (3. + 29. * x_t * x_t) - kc * kc * (3. + 13. * x_t))
403 / 12. / (x_t - kc) / (1. - x_t) / (1. - x_t) * log(x_t) + kc * (17. * x_t - kc) / 12.
404 / (x_t - kc) * log(kc) + x_t * x_t / (x_t - kc) * log(x_t) * log(kc) - (x_t * x_t +
405 2. * x_t * kc - kc * kc) / 2. / (x_t - kc) * log(kc) * log(kc) - x_t * (x_t - kc) /
406 (1. - x_t) * gsl_sf_dilog(1. - x_t) - x_t * gsl_sf_dilog(1. - x_t / kc))
407 + 32. / (x_t - kc) / (1. - x_t)*(x_t * (1. - kc) + kc * (1. - x_t)*(2. * x_t - kc)
408 / (x_t - kc) * log(kc) - x_t * x_t * (1. - kc)*(1. - 2. * x_t + kc) / (x_t - kc)
409 / (1. - x_t) * log(x_t)) * L + 4. * kc / (x_t - kc)*(1. - x_t / (x_t - kc) * log(x_t)
410 + x_t / (x_t - kc) * log(kc)) * xi1c;
411
412 return (C(NNLO,NO_QED,0)(1) + 4. * rhoB2_t + xi1c * C(NLO,NO_QED,0)(1) + 4. * xi1c * rhoB1_t + xi2c * C(LO,NO_QED,0)(1));
413 }
414 default:{
415 std::stringstream out;
416 out << order;
417 throw std::runtime_error("Charm_Kpnunu::C_Bt() order" + out.str() + " not implemented");
418 }
419 }
420}
421
422
423
424
425
427{
428
429 switch (order_qed)
430 {
431 case(NO_QED):{
432 return (0.);
433 }
434 case(LO_QED):{
435 double rhop_e = 0.;
436 double rhom_e = 0.;
437
438 double xi_ce= 8. / 3. / (11. - 2. / 3. * 3.) * (etac - 1.);
439
440 gslpp::vector<double> C0 = C(LO,NO_QED,1);
441 gslpp::vector<double> Ce = C(LO,LO_QED,1);
442
443 return (Ce(2) + C0(2) * xi_ce + 4. * C0(0) / 4. * rhop_e + 4. * C0(0) / 4. * rhom_e);
444 }
445 case(NLO_QED11):{
446 double rhop_e = 0.;
447 double rhom_e = 0.;
448 double rhop_es = 4. * xice;
449 double rhom_es = -2. * xice;
450 double rhop_1 = 4. * (1. - L) + 4. * log(kc);
451 double rhom_1 = -2. * (1. - L) - 2. * log(kc);
452
453 gslpp::vector<double> C0 = C(LO,NO_QED,1);
454 gslpp::vector<double> C1 = C(NLO,NO_QED,1);
455 gslpp::vector<double> Ce = C(LO,LO_QED,1);
456 gslpp::vector<double> Ces = C(LO,NLO_QED11,1);
457
458 return (Ces(2) + Ce(2) * xi1c + C1(2) * xice + C0(2) * xices +
459 4. * (rhop_es + rhop_e * xi1c + rhop_1 * xice) * C0(0) / 4. +
460 4. * (rhom_es + rhom_e * xi1c + rhom_1 * xice) * C0(1) / 4. +
461 4. * rhop_1 * Ce(0) / 4. +
462 4. * rhom_1 * Ce(1) / 4. +
463 4. * rhop_e * C1(0) / 4. + 4. * rhom_e * C1(1) / 4.);
464 }
465 default:{
466 std::stringstream out;
467 out << order_qed;
468 throw std::runtime_error("Charm_Kpnunu::C_P_qed() order" + out.str() + " not implemented");
469 }
470 }
471}
472
473
474
475
477{
478 switch (order_qed)
479 {
480 case(NO_QED):{
481 return (0.);
482 }
483 case(LO_QED):{
484 double rho_e = 0.;
485
486 gslpp::vector<double> C0 = C(LO,NO_QED,0);
487 gslpp::vector<double> Ce = C(LO,LO_QED,0);
488
489 return (Ce(1) + C0(1) * xice + 4. * C0(0) / 4. * rho_e);
490 }
491 case(NLO_QED11):{
492 double rho_e = 0.;
493 double rho_es = -4. * xice;
494 double rho_1 = 5. + 4. * L - 4. * log(kc);
495
496 gslpp::vector<double> C0 = C(LO,NO_QED,0);
497 gslpp::vector<double> Ce = C(LO,LO_QED,0);
498 gslpp::vector<double> C1 = C(NLO,NO_QED,0);
499 gslpp::vector<double> Ces = C(LO,NLO_QED11,0);
500
501 return (Ces(1) + Ce(1) * xi1c + C1(1) * xice + C0(1) * xices +
502 4. * C0(0) / 4. * rho_es + 4. * C0(0) / 4. * rho_e * xi1c +
503 8. * Ce(0) / 4. * rho_1 + 4. * C0(0) / 4. * rho_1 * xice);
504 }
505
506 default:{
507 std::stringstream out;
508 out << order_qed;
509 throw std::runtime_error("Charm_Kpnunu::C_Be_qed() order" + out.str() + " not implemented");
510 }
511 }
512
513}
514
516{
517
518 switch (order_qed)
519 {
520 case(NO_QED):{
521 return (0.);
522 }
523 case(LO_QED):{
524 double rho_e = 0.;
525
526 gslpp::vector<double> C0 = C(LO,NO_QED,0);
527 gslpp::vector<double> Ce = C(LO,LO_QED,0);
528
529
530 return (Ce(1) + C0(1) * xice + 4. * C0(0) / 4. * rho_e);
531 }
532 case(NLO_QED11):{
534
535 double rho_e = 0.;
536 double rho_es = -4. * kc * xice * (kc - xt * (1 - log(xt / kc) ) ) / (kc - xt) / (kc - xt) ;
537 double rho_1 = 5. + 4. * L + 4. * xt / (1. - xt) + 4. / (xt - kc)*(kc * log(kc)
538 - xt * (1. - kc) / (1. - xt) * log(xt));
539
540 gslpp::vector<double> C0 = C(LO,NO_QED,0);
541 gslpp::vector<double> Ce = C(LO,LO_QED,0);
542 gslpp::vector<double> C1 = C(NLO,NO_QED,0);
543 gslpp::vector<double> Ces = C(LO,NLO_QED11,0);
544
545 return (Ces(1) + Ce(1) * xi1c + C1(1) * xice + C0(1) * xices +
546 4. * C0(0) / 4. * rho_es + 4. * C0(0) / 4. * rho_e * xi1c +
547 8. * Ce(0) / 4. * rho_1 + 4. * C0(0) / 4. * rho_1 * xice);
548 }
549
550 default:{
551 std::stringstream out;
552 out << order_qed;
553 throw std::runtime_error("Charm_Kpnunu::C_Be_qed() order" + out.str() + " not implemented");
554 }
555 }
556}
557
558std::vector<WilsonCoefficient>& Charm_Kpnunu::EVOCkpnn()
559{
560 double co = 4. * model.getGF() / sqrt(2.) * model.alphaMz() / 2. / M_PI / model.sW2_ND() ; //SM prefactor as in eq. (1.1) of arXiv:0805.4119
561 double cpeng = kc * xc_mc_qed / 32. ;
562 double cbox = 2. * cpeng ;
563 gslpp::complex lam_c = model.getCKM().computelamc();
564
565 vevoCkpnn.clear();
566
568
569 switch (evoCkpnn.getOrder()) {
570 case NNLO:
571 evoCkpnn.setCoeff(0, co * lam_c * model.Als(model.getMuc(),FULLNLO) / 4. / M_PI * (cpeng * C_P(NNLO) + cbox * (2. / 3. * C_Be(NNLO) + 1. / 3. * C_Bt(NNLO))) , NNLO);
572 case NLO:
573 evoCkpnn.setCoeff(0, co * lam_c * (cpeng * C_P(NLO) + cbox * (2. / 3. * C_Be(NLO) + 1. / 3. * C_Bt(NLO))) , NLO);
574 case LO:
575 evoCkpnn.setCoeff(0, co * lam_c * 4. * M_PI / model.Als(model.getMuc(),FULLNLO) * (cpeng * C_P(LO) + cbox * (2. / 3. * C_Be(LO) + 1. / 3. * C_Bt(LO))) , LO);
576 break;
577 default:
578 std::stringstream out;
579 out << evoCkpnn.getOrder();
580 throw std::runtime_error("Charm_Kpnunu::EVOCkpnn(): order " + out.str() + " not implemented");
581 }
582
583 switch (evoCkpnn.getOrder_qed()) {
584 case NLO_QED11:
585 evoCkpnn.setCoeff(0, co * lam_c * model.Ale(model.getMuc(),FULLNLO) / model.Als(model.getMuc(),FULLNLO) * (cpeng * C_P_qed(NLO_QED11) + cbox * (2. / 3. * C_Be_qed(NLO_QED11) + 1. / 3. * C_Bt_qed(NLO_QED11))) , NLO_QED11);
586 case LO_QED:
587 evoCkpnn.setCoeff(0, co * lam_c * 4. * M_PI * model.Ale(model.getMuc(),FULLNLO) / model.Als(model.getMuc(),FULLNLO) / model.Als(model.getMuc(),FULLNLO) * (cpeng * C_P_qed(LO_QED) + cbox * (2. / 3. * C_Be_qed(LO_QED) + 1. / 3. * C_Bt_qed(LO_QED))) , LO_QED);
588 //case NO_QED:
589 // evoCkpnn.setCoeff(0, 0. , NO_QED);
590 break;
591 default:
592 std::stringstream out;
593 out << evoCkpnn.getOrder_qed();
594 throw std::runtime_error("Charm_Kpnunu::EVOCkpnn(): qed order " + out.str() + " not implemented");
595 }
596
597 vevoCkpnn.push_back(evoCkpnn);
598 return(vevoCkpnn);
599}
600
601
602gslpp::matrix<double> Charm_Kpnunu::ADM(orders order,orders_qed order_qed, double nf,int contribution){
603
604 double nfu,nfd;
605
606 if(nf==1){
607 nfu=1.;
608 nfd=0;
609 }else{
610 nfu = floor(nf/2.) ;
611 nfd = nf-nfu ;
612 }
613
614
615 // QED ADMs are taken from: 0805.4119v1
616 // QCD ADMs are taken from: hep-ph/0603079v3
617 // There is a "transpose difference" between the penguin ADM of the 2 articles. For now the QCD article is used for the definition.
618 switch (contribution)
619 {
620 case (0): // BOX
621 {
622 switch (order_qed)
623 {
624 case (NO_QED): // no qed => go to qcd
625 {
626 switch (order)
627 {
628 case (LO): // 0
629 {
630 double gamma_nuB = -8. ;
631 double gamma_M = 8. ;
632 double beta = 11. - 2. / 3. * nf ;
633 double gamma_nu = 2. * ( gamma_M - beta ) ;
634
635 gslpp::matrix<double> gamma_T(2, 0.);
636
637 gamma_T(0,0) = 0.;
638 gamma_T(1,0) = gamma_nuB;
639 gamma_T(1,1) = gamma_nu;
640
641 return (gamma_T);
642
643 break;
644
645 }
646 case (NLO): // 1
647 {
648 double gamma_nuB = 8. * model.getCF() ;
649 double gamma_M = 404. / 3. - 40. / 9. * nf ;
650 double beta = 102. - 38. / 3. * nf ;
651 double gamma_nu = 2. * ( gamma_M - beta ) ;
652
653 gslpp::matrix<double> gamma_T(2, 0.);
654
655 gamma_T(0,0) = 0.;
656 gamma_T(1,0) = gamma_nuB;
657 gamma_T(1,1) = gamma_nu;
658
659 return (gamma_T);
660
661 break;
662
663 }
664 case(NNLO): // 2
665 {
666 double gamma_nuB = 2. * model.getCF() * (69. / 3. - 458./3.*3. - (48. / 3. - 96. * 3.) * gsl_sf_zeta_int(3) + 38. / 3. * nf) ;
667 double gamma_M = 2498. - (4432./27. + 320.*gsl_sf_zeta_int(3)/3.) * nf - (280. * nf * nf )/81. ;
668 double beta = 2857./2. - 5033.*nf/18. + 325.*nf*nf/54. ;
669 double gamma_nu = 2. * ( gamma_M - beta ) ;
670
671 gslpp::matrix<double> gamma_T(2, 0.);
672
673 gamma_T(0,0) = 0.;
674 gamma_T(1,0) = gamma_nuB;
675 gamma_T(1,1) = gamma_nu;
676
677 return (gamma_T);
678
679 break;
680
681 }
682 default:
683 std::stringstream out;
684 out << order;
685 throw std::runtime_error("Charm_Kpnunu::ADM: order " + out.str() + " not implemented\n");
686
687 }
688 }
689 case (LO_QED): // e
690 {
691 double gamma_W = -4. ;
692 double gamma_nuB = 0. ;
693 double gamma_M = 8. / 3. ;
694 double beta = 0. ;
695 double gamma_nu = 2. * ( gamma_M - beta ) ;
696
697 gslpp::matrix<double> gamma_T(2, 0.);
698
699 gamma_T(0,0) = 2. * gamma_W;
700 gamma_T(1,0) = gamma_nuB;
701 gamma_T(1,1) = gamma_nu;
702
703 return (gamma_T);
704
705 break;
706 }
707 case (NLO_QED11): // es
708 {
709 double gamma_W = 4. ;
710 double gamma_nuB = -316. / 9. ;
711 double gamma_M = 32. / 9. ;
712 double beta = -8. / 9. * (nfu + nfd/4. ) ;
713 double gamma_nu = 2. * ( gamma_M - beta ) ;
714
715 gslpp::matrix<double> gamma_T(2, 0.);
716
717 gamma_T(0,0) = 2. * gamma_W;
718 gamma_T(1,0) = gamma_nuB;
719 gamma_T(1,1) = gamma_nu;
720
721 return (gamma_T);
722
723 break;
724 }
725 default:
726 std::stringstream out, out2;
727 out << order_qed;
728 out2 << contribution;
729 throw std::runtime_error("Charm_Kpnunu::ADM: order_qed " + out.str() + " of contribution " + out2.str() + " not implemented\n");
730 }
731 }
732 case (1): // PENGUIN
733 {
734 switch (order_qed)
735 {
736 case (NO_QED): // no qed => go to qcd
737 {
738 switch (order)
739 {
740 case (LO): // 0
741 {
742 gslpp::matrix<double> gamma_pm(2,0.); // from ref 29 of BG08
743 gamma_pm(0,0)= +6.*(1.-1./3.);
744 gamma_pm(0,1)= 0.;
745 gamma_pm(1,0)= 0.;
746 gamma_pm(1,1)= -6.*(1.+1./3.);
747
748 double gamma_pnu = -0.5 * ( -4. * ( 1. + 3. ) ) ;
749 double gamma_mnu = -0.5 * ( -4. * ( 1. - 3. ) ) ;
750 double gamma_M = 8. ;
751 double beta = 11. - 2. / 3. * nf ;
752 double gamma_nu = 2. * ( gamma_M - beta ) ;
753
754 gslpp::matrix<double> gamma_T(3, 0.);
755
756 gamma_T(0,0) = gamma_pm(0,0);
757 gamma_T(1,0) = gamma_pm(0,1); // inverted beacuse is transposed
758 gamma_T(0,1) = gamma_pm(1,0);
759 gamma_T(1,1) = gamma_pm(1,1);
760
761 gamma_T(2,0) = gamma_pnu;
762 gamma_T(2,1) = gamma_mnu;
763
764 gamma_T(2,2) = gamma_nu;
765
766 return (gamma_T);
767 break;
768
769 }
770 case (NLO): // 1
771 {
772 gslpp::matrix<double> gamma_pm(2,0.); // from ref 29 of BG08
773 gamma_pm(0,0)= (-21./2. + 2.*nf/3.) * (1. - 1./3.);
774 gamma_pm(0,1)= 0.;
775 gamma_pm(1,0)= 0.;
776 gamma_pm(1,1)= (-21./2. - 2.*nf/3.) * (1. + 1./3.);
777
778 double gamma_pnu = -0.5 * ( 16. * (2. - 11.) ) ;
779 double gamma_mnu = -0.5 * ( 16. * (2. + 11.) ) ;
780 double gamma_M = 404. / 3. - 40. / 9. * nf ;
781 double beta = 102. - 38. / 3. * nf ;
782 double gamma_nu = 2. * ( gamma_M - beta ) ;
783
784 gslpp::matrix<double> gamma_T(3, 0.);
785
786 gamma_T(0,0) = gamma_pm(0,0);
787 gamma_T(1,0) = gamma_pm(0,1); // inverted beacuse is transposed
788 gamma_T(0,1) = gamma_pm(1,0);
789 gamma_T(1,1) = gamma_pm(1,1);
790
791 gamma_T(2,0) = gamma_pnu;
792 gamma_T(2,1) = gamma_mnu;
793
794 gamma_T(2,2) = gamma_nu;
795
796 return (gamma_T);
797 break;
798
799 }
800 case(NNLO): // 2
801 {
802 gslpp::matrix<double> gamma_pm(2,0.); // from ref 29 of BG08
803 gamma_pm(0,0)= 1./300. * (349049. + 201485.) - 1./1350. * (115577. - 9795.)*nf - 130./27. * (1. - 1./3.)*nf*nf - (672. + 80. * (1. - 1./3.)*nf ) * gsl_sf_zeta_int(3);
804 gamma_pm(0,1)= 0.;
805 gamma_pm(1,0)= 0.;
806 gamma_pm(1,1)= 1./300. * (349049. - 201485.) - 1./1350. * (115577. + 9795.)*nf + 130./27. * (1. + 1./3.)*nf*nf + (672. + 80. * (1. + 1./3.)*nf ) * gsl_sf_zeta_int(3);
807
808 double gamma_pnu = -0.5 * ( -2./225. * (45124. + 484917.) + 32.*(13. + 15.)*gsl_sf_zeta_int(3) + 144.*nf ) ;
809 double gamma_mnu = -0.5 * ( -2./225. * (45124. - 484917.) + 32.*(13. - 15.)*gsl_sf_zeta_int(3) - 144.*nf );
810 double gamma_M = 2498. - (4432./27. + 320.*gsl_sf_zeta_int(3)/3.) * nf - (280. * nf * nf )/81. ;
811 double beta = 2857./2. - 5033.*nf/18. + 325.*nf*nf/54. ;
812 double gamma_nu = 2. * ( gamma_M - beta ) ; ;
813
814 gslpp::matrix<double> gamma_T(3, 0.);
815
816 gamma_T(0,0) = gamma_pm(0,0);
817 gamma_T(1,0) = gamma_pm(0,1); // inverted beacuse is transposed
818 gamma_T(0,1) = gamma_pm(1,0);
819 gamma_T(1,1) = gamma_pm(1,1);
820
821 gamma_T(2,0) = gamma_pnu;
822 gamma_T(2,1) = gamma_mnu;
823
824 gamma_T(2,2) = gamma_nu;
825
826 return (gamma_T);
827 break;
828
829 }
830 default:
831 std::stringstream out;
832 out << order;
833 throw std::runtime_error("Charm_Kpnunu::ADM: order " + out.str() + " not implemented\n");
834
835 }
836 }
837 case (LO_QED): // e
838 {
839 gslpp::matrix<double> gamma_pm(2,0.); // from ref 29 of BG08
840 gamma_pm(0,0)= -8./3.;
841 gamma_pm(0,1)= 0.;
842 gamma_pm(1,0)= 0.;
843 gamma_pm(1,1)= -8./3.;
844
845 double gamma_pnu = 0. ;
846 double gamma_mnu = 0. ;
847 double gamma_M = 8. / 3. ;
848 double beta = 0. ;
849 double gamma_nu = 2. * ( gamma_M - beta ) ;
850
851 gslpp::matrix<double> gamma_T(3, 0.);
852
853 gamma_T(0,0) = gamma_pm(0,0);
854 gamma_T(1,0) = gamma_pm(0,1);
855 gamma_T(0,1) = gamma_pm(1,0);
856 gamma_T(1,1) = gamma_pm(1,1);
857
858 gamma_T(2,0) = gamma_pnu;
859 gamma_T(2,1) = gamma_mnu;
860
861 gamma_T(2,2) = gamma_nu;
862
863 return (gamma_T);
864
865 break;
866 }
867 case (NLO_QED11): // es
868 {
869 gslpp::matrix<double> gamma_12(2,0.); // from ref 29 of BG08 (look out for the different basis)
870 gslpp::matrix<complex> trash(2,0.);
871 gslpp::vector<complex> ev_gamma12(2,0.);
872
873 gamma_12.assign(0,0,194./9.);
874 gamma_12.assign(0,1,-2./3.);
875 gamma_12.assign(1,0,25./3.);
876 gamma_12.assign(1,1,-49./9.) ;
877
878 gamma_12.eigensystem(trash,ev_gamma12);
879
880 double gamma_pnu = 0. ;
881 double gamma_mnu = 0. ;
882 double gamma_M = 32. / 9. ;
883 double beta = -8. / 9. * (nfu + nfd/4. ) ;
884 double gamma_nu = 2. * ( gamma_M - beta ) ;
885
886 gslpp::matrix<double> gamma_T(3, 0.);
887
888 gamma_T(0,0) = ev_gamma12(0).real();
889 gamma_T(1,1) = ev_gamma12(1).real();
890
891 gamma_T(2,0) = gamma_pnu;
892 gamma_T(2,1) = gamma_mnu;
893
894 gamma_T(2,2) = gamma_nu;
895
896 return (gamma_T);
897
898 break;
899 }
900 default:
901 std::stringstream out, out2;
902 out << order_qed;
903 out2 << contribution;
904 throw std::runtime_error("Charm_Kpnunu::ADM: order_qed " + out.str() + " of contribution " + out2.str() + " not implemented\n");
905 break;
906 }
907 break;
908 }
909 default:
910 std::stringstream out ;
911 out << contribution;
912 throw std::runtime_error("Charm_Kpnunu::ADM: contribution " + out.str() + " not implemented\n");
913 break;
914 }
915
916
917}
918
919gslpp::vector<double> Charm_Kpnunu::CWin_muw(orders order,orders_qed order_qed, int contribution)
920{
921 switch (contribution)
922 {
923 case(0): // box
924 {
925 switch (order_qed)
926 {
927 case(NO_QED): // no qed => go to qcd
928 {
929 switch (order)
930 {
931
932 case(LO): // 0
933 {
934 gslpp::vector<double> CWin(2, 0.);
935 double Cnu= 0. ;
936
937 CWin(0)= 4. ;
938 CWin(1)= Cnu ;
939
940 return (CWin);
941 break;
942
943 }
944 case (NLO):
945 {
946 double MW = model.Mw();
947 double LW = log(model.getMuw() * model.getMuw() / MW / MW);
948
949 gslpp::vector<double> CWin(2, 0.);
950 double Cnu= -4. * (9. + 4. * LW) ;
951
952 CWin(0)= 0. ;
953 CWin(1)= Cnu ;
954
955 return (CWin);
956 break;
957
958 }
959 case (NNLO):
960 {
961 double MW = model.Mw();
962 double LW = log(model.getMuw() * model.getMuw() / MW / MW);
963
964 gslpp::vector<double> CWin(2, 0.);
965 double Cnu= -8. * model.getCF()*(20. + 2. * M_PI * M_PI + 25. * LW + 6. * LW * LW) ;
966
967 CWin(0)= 0. ;
968 CWin(1)= Cnu ;
969
970 return (CWin);
971 break;
972
973 }
974
975 default:
976 std::stringstream out, out2;
977 out << order;
978 out2 << contribution;
979 throw std::runtime_error("Charm_Kpnunu::CWin_muw: order " + out.str() + " of contribution " + out2.str() +" not implemented\n");
980 break;
981
982 }
983 }
984 case(LO_QED):// e
985 {
986 gslpp::vector<double> CWin(2, 0.);
987 double CW= 0.;
988 double Cnu= 0. ;
989
990 CWin(0)= 4. * CW * CW ;
991 CWin(1)= Cnu ;
992
993 return (CWin);
994 break;
995 }
996 case(NLO_QED11):// es
997 {
998 double MZ = model.getMz();
999 double LZ = log(model.getMuw() * model.getMuw() / MZ / MZ);
1000
1001 gslpp::vector<double> CWin(2, 0.);
1002 double CW= -11. / 3. - 2. * LZ ;
1003 double Cnu= 0. ;
1004
1005 CWin(0)= 4. * CW * CW ;
1006 CWin(1)= Cnu ;
1007
1008 return (CWin);
1009 break;
1010 }
1011 default:
1012 std::stringstream out, out2;
1013 out << order_qed;
1014 out2 << contribution;
1015 throw std::runtime_error("Charm_Kpnunu::Cwin_qed: order_qed " + out.str() + " of contribution " + out2.str() +" not implemented\n");
1016 break;
1017 }
1018 break;
1019 }
1020
1021
1022 case(1): // penguin
1023 {
1024 switch (order_qed)
1025 {
1026 case(NO_QED): // no qed => go to qcd
1027 {
1028 switch (order)
1029 {
1030
1031 case(LO):
1032 {
1033 gslpp::vector<double> CWin(3, 0.);
1034 double Cp= 1. ;
1035 double Cm= 1. ;
1036 double Cnu= 0. ;
1037
1038 CWin(0)= 4. * Cp ;
1039 CWin(1)= 4. * Cm ;
1040 CWin(2)= Cnu ;
1041
1042 return (CWin);
1043
1044 break;
1045
1046 }
1047 case (NLO):
1048 {
1049 double MW = model.Mw();
1050 double LW = log(model.getMuw() * model.getMuw() / MW / MW);
1051 gslpp::vector<double> CWin(3, 0.);
1052
1053 double Cp= 0.5*(1. - 1./3.) * (11. + 6. * LW) ;
1054 double Cm= -0.5*(1. + 1./3.) * (11. + 6. * LW) ;
1055 double Cnu= 8. * (2. + LW) ;
1056
1057 CWin(0)= 4. * Cp ;
1058 CWin(1)= 4. * Cm ;
1059 CWin(2)= Cnu ;
1060
1061 return (CWin);
1062
1063 break;
1064
1065 }
1066 case (NNLO):
1067 {
1068 double MW = model.Mw();
1069 double LW = log(model.getMuw() * model.getMuw() / MW / MW);
1070 double x = model.getMatching().x_t(model.getMut());
1071
1072 gslpp::vector<double> CWin(3, 0.);
1073
1074 double Cp= (-(135677. - 124095.) / 3600. + 58. / 18. * M_PI * M_PI - 0.5 * (2. / 3.)*
1075 (112. / 9. + 32. * x + (20. / 3. + 16. * x) * log(x) - (8. + 16. * x) *
1076 sqrt(4. * x - 1.) * gsl_sf_clausen(2. * asin(1. / (2. * sqrt(x)))))
1077 +(5. / 36. * 238. * LW) + 58. / 6. * LW * LW) ;
1078 double Cm= (-(135677. + 124095.) / 3600. - 44. / 18. * M_PI * M_PI + 0.5 * (4. / 3.)*
1079 (112. / 9. + 32. * x + (20. / 3. + 16. * x) * log(x) - (8. + 16. * x) *
1080 sqrt(4. * x - 1.) * gsl_sf_clausen(2. * asin(1. / (2. * sqrt(x)))))
1081 - (5. / 36. * 260. * LW) - 44. / 6. * LW * LW);
1082 double Cnu= 4. * model.getCF() * (33. + 4. * M_PI * M_PI + 34. * LW + 12. * LW * LW) ;
1083
1084 CWin(0)= 4. * Cp ;
1085 CWin(1)= 4. * Cm ;
1086 CWin(2)= Cnu ;
1087
1088 return (CWin);
1089
1090 break;
1091
1092 }
1093
1094 default:
1095 std::stringstream out, out2;
1096 out << order;
1097 out2 << contribution;
1098 throw std::runtime_error("Charm_Kpnunu::CWin_muw: order " + out.str() + " of contribution " + out2.str() +" not implemented\n");
1099 break;
1100
1101 }
1102 }
1103 case(LO_QED): // e
1104 {
1105 gslpp::vector<double> CWin(3, 0.);
1106 double Cp= 0. ;
1107 double Cm= 0. ;
1108 double CA= 0. ;
1109 double Cnu= 0. ;
1110
1111 CWin(0)= 4. * Cp * CA ;
1112 CWin(1)= 4. * Cm * CA ;
1113 CWin(2)= Cnu ;
1114
1115 return (CWin);
1116
1117 break;
1118 }
1119 case(NLO_QED11): // es
1120 {
1121 double mt = model.getMatching().getMt_mut();
1122 double MH = model.getMHl();
1123 double MW = model.Mw();
1124 double MZ = model.getMz();
1125 double LZ = log(model.getMuw() * model.getMuw() / MZ / MZ);
1126
1127 double mt2 = mt * mt;
1128 double MW2 = MW * MW;
1129 double MZ2 = MZ * MZ;
1130 double MH2 = MH * MH;
1131 double MH4 = MH2 * MH2;
1132 double sw2 = model.sW2_MSbar_Approx();
1133 double sw4 = sw2 * sw2;
1134 double cw2 = 1. - sw2;
1135
1136 gslpp::vector<double> CWin(3, 0.);
1137 double Cp= -22. / 9. - 4. / 3. * LZ ;
1138 double Cm= -22. / 9. - 4. / 3. * LZ ;
1139 double CA= 3. * mt2 / 4. / sw2 / MW2 + (11. * sw2 - 6.) / 4. / sw2 / cw2 -
1140 3. / 4. * (MW2 - cw2 * MH2) / (MH2 - MW2) / sw4 * log(MW2 / MZ2) +
1141 3. * MH4 / 4. / (MH2 - MW2) / (MW2 - cw2 * MH2) * log(MH2 / MZ2);
1142 double Cnu= 0. ;
1143
1144 CWin(0)= 4. * Cp * CA ;
1145 CWin(1)= 4. * Cm * CA ;
1146 CWin(2)= Cnu ;
1147
1148 return (CWin);
1149 break;
1150 }
1151 default:
1152 std::stringstream out, out2;
1153 out << order_qed;
1154 out2 << contribution;
1155 throw std::runtime_error("Charm_Kpnunu::Cwin_qed: order_qed " + out.str() + " of contribution " + out2.str() +" not implemented\n");
1156 break;
1157 }
1158 break;
1159 }
1160 default:
1161 std::stringstream out ;
1162 out << contribution;
1163 throw std::runtime_error("Charm_Kpnunu::Cwin_qed: contribution " + out.str() +" not implemented\n");
1164 break;
1165
1166 }
1167
1168
1169}
1170
1171gslpp::matrix<double> Charm_Kpnunu::RGevol_J(orders order, int nf, int contribution)
1172{
1173 int dimension;
1174
1175 switch (contribution)
1176 {
1177 case(0): // box
1178 dimension=2;
1179 break;
1180
1181 case(1): // penguin
1182 dimension=3;
1183 break;
1184
1185 default:
1186 std::stringstream out ;
1187 out << contribution;
1188 throw std::runtime_error("Charm_Kpnunu::RGevol_J : contribution " + out.str() +" not implemented\n");
1189 break;
1190
1191 }
1192
1193 gslpp::matrix<gslpp::complex> V(dimension,0.) , V_inv(dimension,0.) ;
1194 gslpp::matrix<double> G1(dimension,0.) ;
1195 gslpp::matrix<double> S1(dimension,0.) ;
1196 gslpp::matrix<double> J1(dimension,0.) ;
1197 gslpp::matrix<double> ADM_0(dimension,0.) , ADM_1(dimension,0.) , ADM_2(dimension,0.);
1198
1199
1200 gslpp::vector<gslpp::complex> e(dimension,0.) ;
1201
1202 double beta_0 = 11. - 2. / 3. * nf ;
1203 double beta_1 = 102. - 38. / 3. * nf ;
1204 double beta_2 = 2857. / 2. - 5033. / 18. * nf + 325. / 54. * nf * nf ;
1205
1206
1207 ADM_0 = ADM(LO,NO_QED,nf,contribution); // remember that these matrices are all transposed as in the reference
1208 ADM_1 = ADM(NLO,NO_QED,nf,contribution);
1209 ADM_2 = ADM(NNLO,NO_QED,nf,contribution);
1210
1211 ADM_0.eigensystem(V,e);
1212
1213 V_inv=V.inverse();
1214
1215 G1 = V_inv.real() * ADM_1 * V.real() ;
1216
1217 //S1
1218 for(unsigned int i = 0; i < G1.size_i(); i++){
1219 for (unsigned int j = 0; j < G1.size_j(); j++) {
1220 S1.assign(i , j, beta_1 / 2. / beta_0 / beta_0 * e(i).real() * (i==j) - G1(i,j) / (2. * beta_0 + e(i).real() - e(j).real() ));
1221 }
1222 }
1223
1224 J1 = V.real() * S1 * V_inv.real() ;
1225
1226 if (order==LO){
1227 switch (contribution){
1228 case(0): //box
1229 return RGevolB(nf,LO);
1230 case(1): // penguin
1231 return RGevolP(nf,LO);
1232 default:
1233 std::stringstream out ;
1234 out << contribution;
1235 throw std::runtime_error("Charm_Kpnunu::RGevol_J : contribution " + out.str() +" not implemented\n");
1236 }
1237 } else if (order==NLO){
1238 return J1;
1239 } else if (order==NNLO){
1240 gslpp::matrix<double> S2(dimension,0.) , G2(dimension,0.) , J2(dimension,0.);
1241 G2 = V_inv.real() * ADM_2 * V.real() ;
1242 for(unsigned int i = 0; i < G1.size_i(); i++){
1243 for (unsigned int j = 0; j < G1.size_j(); j++) {
1244 double term=0. ;
1245 for (unsigned int k = 0; k < G1.size_j(); k++){
1246 term += ( 1. + e(i).real() / 2. / beta_0 - e(k).real() / 2. / beta_0 ) / ( 2. + e(i).real() / 2. / beta_0 - e(j).real() / 2. / beta_0 ) * ( S1(i,k) * S1(k,j) - beta_1 / beta_0 * S1(i,j) * (j==k) );
1247 }
1248 S2.assign(i , j , beta_2 / 4. / beta_0 / beta_0 * e(i).real() * (i==j) + term - G2(i,j) / (4. * beta_0 + e(i).real() - e(j).real() ));
1249 }
1250 }
1251 J2 = V.real() * S2 * V_inv.real() ;
1252 return J2;
1253
1254 } else{
1255 std::stringstream out ;
1256 out << order;
1257 throw std::runtime_error("Charm_Kpnunu::RGevol_J : order " + out.str() +" not implemented\n");
1258 }
1259}
1260
1261
1262gslpp::matrix<double> Charm_Kpnunu::RGevol_R(orders_qed order_qed, int nf, int contribution)
1263{
1264 int dimension;
1265
1266 switch (contribution)
1267 {
1268 case(0): // box
1269 dimension=2;
1270 break;
1271
1272 case(1): // penguin
1273 dimension=3;
1274 break;
1275
1276 default:
1277 std::stringstream out ;
1278 out << contribution;
1279 throw std::runtime_error("Charm_Kpnunu::RGevol_K : contribution " + out.str() +" not implemented\n");
1280 break;
1281 }
1282
1283
1284 gslpp::matrix<gslpp::complex> V(dimension,0.) , V_inv(dimension,0.) ;
1285 gslpp::matrix<double> M_0(dimension,0.) , K_0(dimension,0.) ;
1286 gslpp::matrix<double> ADM_0(dimension,0.) , ADM_1(dimension,0.) , ADM_2(dimension,0.), ADM_e(dimension,0.), ADM_es(dimension,0.);
1287 gslpp::matrix<double> M_1(dimension,0.), S_1(dimension,0.) , K1_1(dimension,0.) , K2_1(dimension,0.) , K3_1(dimension,0.), G_1(dimension,0.) ;
1288
1289
1290 gslpp::vector<gslpp::complex> e(dimension,0.) ;
1291
1292 double beta_0 = 11. - 2. / 3. * nf ;
1293 double beta_1 = 102. - 38. / 3. * nf ;
1294 //double beta_2 = 2857. / 2. - 5033. / 18. * nf + 325. / 54. * nf * nf ;
1295
1296
1297 ADM_0 = ADM(LO,NO_QED,nf,contribution); // remember that these matrices are all transposed as in the reference
1298 ADM_e = ADM(LO,LO_QED,nf,contribution); // remember that these matrices are all transposed as in the reference
1299 ADM_es = ADM(LO,NLO_QED11,nf,contribution); // remember that these matrices are all transposed as in the reference
1300 ADM_1 = ADM(NLO,NO_QED,nf,contribution);
1301 ADM_2 = ADM(NNLO,NO_QED,nf,contribution);
1302
1303 ADM_0.eigensystem(V,e);
1304 V_inv=V.inverse();
1305
1306
1307 M_0 = V_inv.real() * ADM_e * V.real();
1308 M_1 = V_inv.real() * ( (ADM_es - beta_1 / beta_0 * ADM_e) + (ADM_e * RGevol_J(NLO,nf,contribution) - RGevol_J(NLO,nf,contribution) * ADM_e) ) * V.real();
1309 G_1 = V_inv.real() * ADM_1 * V.real();
1310 //S1
1311 for(unsigned int i = 0; i < G_1.size_i(); i++){
1312 for (unsigned int j = 0; j < G_1.size_j(); j++) {
1313 S_1.assign(i , j, beta_1 / 2. / beta_0 / beta_0 * e(i).real() * (i==j) - G_1(i,j) / (2. * beta_0 + e(i).real() - e(j).real() ));
1314 }
1315 }
1316 switch(nf){
1317 case(4):
1318 for(unsigned int i = 0; i < M_0.size_i(); i++){
1319 for (unsigned int j = 0; j < M_0.size_j(); j++) {
1320 if (e(i).real() != e(j).real() + 2. * beta_0){
1321 K_0.assign(i , j, M_0(i,j) / (e(i).real() / 2. / beta_0 - e(j).real() / 2. / beta_0 - 1.) * ( pow(etacb,e(j).real() /2. / beta_0) - pow(etacb,e(i).real() /2. / beta_0) / etacb) );
1322 } else {
1323 K_0.assign(i , j, M_0(i,j) * pow(etacb, e(j).real() / 2. / beta_0) * log( 1. / etacb) );
1324 }
1325 }
1326 }
1327 break;
1328 case(5):
1329 for(unsigned int i = 0; i < M_0.size_i(); i++){
1330 for (unsigned int j = 0; j < M_0.size_j(); j++) {
1331 if (e(i).real() != e(j).real() + 2. * beta_0){
1332 K_0.assign(i , j, M_0(i,j) / (e(i).real() / 2. / beta_0 - e(j).real() / 2. / beta_0 - 1.) * ( pow(etab,e(j).real() /2. / beta_0) / etacb - pow(etab,e(i).real() /2. / beta_0) / etacb / etab) );
1333 } else {
1334 K_0.assign(i , j, M_0(i,j) * pow(etab, e(j).real() / 2. / beta_0) * log( 1. / etab) / etacb );
1335 }
1336 }
1337 }
1338 break;
1339 default:
1340 std::stringstream out ;
1341 out << nf;
1342 throw std::runtime_error("Charm_Kpnunu::RGevol_K : nf " + out.str() +" not implemented\n");
1343 break;
1344 }
1345
1346
1347 switch (nf)
1348 {
1349 case(4): //b -> c
1350 switch(order_qed){
1351 case(LO_QED):
1352 return ((-2. * M_PI / beta_0) * (V.real() * K_0 * V_inv.real()) );
1353 break;
1354 case(NLO_QED11):
1355 for(unsigned int i = 0; i < M_1.size_i(); i++){
1356 for (unsigned int j = 0; j < M_1.size_j(); j++) {
1357 if (i!=j){
1358 K1_1.assign(i , j, M_1(i,j) / (e(i).real() / 2. / beta_0 - e(j).real() / 2. / beta_0) * ( pow(etacb,e(j).real() /2. / beta_0) - pow(etacb,e(i).real() /2. / beta_0) ) );
1359 } else {
1360 K1_1.assign(i , j, M_1(i,j) * pow(etacb, e(i).real() / 2. / beta_0) * log( 1. / etacb) );
1361 }
1362 }
1363 }
1364 K2_1 = - etacb * K_0 * S_1;
1365 K3_1 = S_1 * K_0;
1366
1367 return ( ( -0.5 / beta_0) * V.real() * (K1_1 + K2_1 + K3_1) * V_inv.real());
1368 break;
1369 default:
1370 std::stringstream out ;
1371 out << order_qed;
1372 throw std::runtime_error("Charm_Kpnunu::RGevol_K : order " + out.str() +" not implemented\n");
1373 break;
1374 }
1375 case(5): //w -> b
1376 switch(order_qed){
1377 case(LO_QED):
1378 return ((-2. * M_PI / beta_0) * (V.real() * K_0 * V_inv.real()) );
1379 break;
1380 case(NLO_QED11):
1381 for(unsigned int i = 0; i < M_1.size_i(); i++){
1382 for (unsigned int j = 0; j < M_1.size_j(); j++) {
1383 if (i!=j){
1384 K1_1.assign(i , j, M_1(i,j) / (e(i).real() / 2. / beta_0 - e(j).real() / 2. / beta_0) * ( pow(etab,e(j).real() /2. / beta_0) - pow(etab,e(i).real() /2. / beta_0) ) );
1385 } else {
1386 K1_1.assign(i , j, M_1(i,j) * pow(etab, e(i).real() / 2. / beta_0) * log( 1. / etab) );
1387 }
1388 }
1389 }
1390 K2_1 = - etacb * etab * K_0 * S_1;
1391 K3_1 = etacb * S_1 * K_0;
1392
1393 return ( ( -0.5 / beta_0) * V.real() * (K1_1 + K2_1 + K3_1) * V_inv.real());
1394 break;
1395 default:
1396 std::stringstream out ;
1397 out << order_qed;
1398 throw std::runtime_error("Charm_Kpnunu::RGevol_K : order " + out.str() +" not implemented\n");
1399 break;
1400 }
1401
1402 default:
1403 std::stringstream out ;
1404 out << nf;
1405 throw std::runtime_error("Charm_Kpnunu::RGevol_K : nf " + out.str() +" not implemented\n");
1406 break;
1407 }
1408
1409
1410}
1411
1412
1413
1414
1415gslpp::vector<double> Charm_Kpnunu::C(orders order,orders_qed order_qed, int contribution){
1416
1417 switch(contribution){
1418 case(0):{ //box
1419 switch(order_qed){
1420 case (NO_QED):{
1421 switch (order)
1422 {
1423 case (LO):{
1424 return U0_4b * U0_5b * CW0b ;
1425 break;
1426 }
1427 case (NLO):{
1428 return J1_4b * U0_4b * U0_5b * CW0b + etacb * U0_4b * (J1_5b - J1_4b) * U0_5b * CW0b +
1429 etab * etacb * U0_4b * U0_5b * (CW1b - J1_5b * CW0b) ;
1430
1431 break;
1432 }
1433 case(NNLO):{
1434 return J2_4b * U0_4b * U0_5b * CW0b + etacb * J1_4b * U0_4b * (J1_5b - J1_4b) * U0_5b * CW0b +
1435 etab * etacb * J1_4b * U0_4b * U0_5b * (CW1b - J1_5b * CW0b) + etacb * etacb * U0_4b *
1436 (J2_5b - J2_4b - J1_4b * (J1_5b - J1_4b)) * U0_5b * CW0b + etab * etacb * etacb *
1437 U0_4b * (J1_5b - J1_4b) * U0_5b * (CW1b - J1_5b * CW0b) + etab * etab * etacb * etacb * U0_4b * U0_5b *
1438 (CW2b - J1_5b * CW1b - (J2_5b - J1_5b * J1_5b) * CW0b) - ThresholdCb(NNLO);
1439
1440 break;
1441 }
1442 default:{
1443 std::stringstream out ;
1444 out << order;
1445 throw std::runtime_error("Charm_Kpnunu::C : order " + out.str() +" not implemented\n");
1446 break;
1447 }
1448
1449 }
1450 }
1451 case(LO_QED):{
1452 return ( (U0_4b*R0_5b/(4.*M_PI) + R0_4b*U0_5b/(4.*M_PI))*CW0b + U0_4b*U0_5b*CWeb );
1453 break;
1454 }
1455 case(NLO_QED11):{
1456 return ( (U0_4b*R1_5b + (J1_4b*U0_4b*R0_5b)/(4*M_PI) - (U0_4b*J1_4b*R0_5b)/(4*M_PI) + etacb/(4*M_PI)*(R0_4b*J1_5b*U0_5b) -
1457 (etacb*etab)/(4*M_PI)*(R0_4b*U0_5b*J1_5b) + R1_4b*U0_5b)*CW0b + ( (etacb*etab)/(4*M_PI)*U0_4b*R0_5b + (etacb*etab)/(4*M_PI)*R0_4b*U0_5b)*CW1b +
1459 break;
1460 }
1461 default:{
1462 std::stringstream out ;
1463 out << order_qed;
1464 throw std::runtime_error("Charm_Kpnunu::C : order_qed " + out.str() +" not implemented\n");
1465 break;
1466 }
1467 }
1468 }
1469 case(1):{
1470 switch(order_qed){
1471 case (NO_QED):{
1472 switch (order)
1473 {
1474 case (LO):{
1475 return U0_4p * U0_5p * CW0p ;
1476 break;
1477 }
1478 case (NLO):{
1479 return J1_4p * U0_4p * U0_5p * CW0p + etacb * U0_4p * (J1_5p - J1_4p) * U0_5p * CW0p +
1480 etab * etacb * U0_4p * U0_5p * (CW1p - J1_5p * CW0p) ;
1481
1482 break;
1483 }
1484 case(NNLO):{
1485 return J2_4p * U0_4p * U0_5p * CW0p + etacb * J1_4p * U0_4p * (J1_5p - J1_4p) * U0_5p * CW0p +
1486 etab * etacb * J1_4p * U0_4p * U0_5p * (CW1p - J1_5p * CW0p) + etacb * etacb * U0_4p *
1487 (J2_5p - J2_4p - J1_4p * (J1_5p - J1_4p)) * U0_5p * CW0p + etab * etacb * etacb *
1488 U0_4p * (J1_5p - J1_4p) * U0_5p * (CW1p - J1_5p * CW0p) + etab * etab * etacb * etacb * U0_4p * U0_5p *
1489 (CW2p - J1_5p * CW1p - (J2_5p - J1_5p * J1_5p) * CW0p) - ThresholdCp(NNLO);
1490
1491 break;
1492 }
1493 default:{
1494 std::stringstream out ;
1495 out << order;
1496 throw std::runtime_error("Charm_Kpnunu::C : order " + out.str() +" not implemented\n");
1497 break;
1498 }
1499
1500 }
1501 }
1502 case(LO_QED):{
1503 return ( (U0_4p*R0_5p/(4.*M_PI) + R0_4p*U0_5p/(4.*M_PI))*CW0p + U0_4p*U0_5p*CWep );
1504 break;
1505 }
1506 case(NLO_QED11):{
1507 return ( (U0_4p*R1_5p + (J1_4p*U0_4p*R0_5p)/(4*M_PI) - (U0_4p*J1_4p*R0_5p)/(4*M_PI) + etacb/(4*M_PI)*(R0_4p*J1_5p*U0_5p) -
1508 (etacb*etab)/(4*M_PI)*(R0_4p*U0_5p*J1_5p) + R1_4p*U0_5p)*CW0p + ( (etacb*etab)/(4*M_PI)*U0_4p*R0_5p + (etacb*etab)/(4*M_PI)*R0_4p*U0_5p)*CW1p +
1510 break;
1511 }
1512 default:{
1513 std::stringstream out ;
1514 out << order_qed;
1515 throw std::runtime_error("Charm_Kpnunu::C : order_qed " + out.str() +" not implemented\n");
1516 break;
1517 }
1518 }
1519
1520 }
1521 default:{
1522 std::stringstream out ;
1523 out << contribution;
1524 throw std::runtime_error("Charm_Kpnunu::C : contribution " + out.str() +" not implemented\n");
1525 break;
1526 }
1527 }
1528
1529}
@ NNLO
Definition: OrderScheme.h:36
@ LO
Definition: OrderScheme.h:34
@ NLO
Definition: OrderScheme.h:35
@ FULLNNLO
Definition: OrderScheme.h:39
@ FULLNLO
Definition: OrderScheme.h:38
@ NDR
Definition: OrderScheme.h:21
@ NLO_QED11
Definition: OrderScheme.h:59
@ LO_QED
Definition: OrderScheme.h:58
@ NO_QED
Definition: OrderScheme.h:57
const gslpp::complex computelamc() const
The product of the CKM elements .
Definition: CKM.cpp:147
~Charm_Kpnunu()
destructor
gslpp::vector< double > CW1p
Definition: Charm_Kpnunu.h:200
gslpp::vector< double > CW1b
Definition: Charm_Kpnunu.h:197
gslpp::matrix< double > RGevol_R(orders_qed order_qed, int nf, int contribution)
This method returns the R matrices useful for the RGE of the charm penguin/box diagrams in qed.
gslpp::matrix< double > J2_5b
Definition: Charm_Kpnunu.h:198
double C_Be(orders order)
gslpp::vector< double > CWesp
Definition: Charm_Kpnunu.h:200
gslpp::matrix< double > R0_4b
Definition: Charm_Kpnunu.h:198
std::vector< WilsonCoefficient > vevoCkpnn
Definition: Charm_Kpnunu.h:204
gslpp::vector< double > CWesb
Definition: Charm_Kpnunu.h:197
gslpp::matrix< double > J2_5p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > R1_5p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > ADM(orders order, orders_qed order_qed, double nf, int contribution)
Reference: 0805.4119v1 and hep-ph/0603079v3. To retrieve QCD ADM just put NO_QED to order_qed and the...
gslpp::vector< double > CWep
Definition: Charm_Kpnunu.h:200
gslpp::matrix< double > J1_5p
Definition: Charm_Kpnunu.h:201
gslpp::vector< double > CW2b
Definition: Charm_Kpnunu.h:197
gslpp::matrix< double > R0_5p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > J2_4p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > R1_4p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > J2_4b
Definition: Charm_Kpnunu.h:198
gslpp::vector< double > CW0p
Definition: Charm_Kpnunu.h:200
double C_Be_qed(orders_qed order_qed)
gslpp::vector< double > dcp
Definition: Charm_Kpnunu.h:193
gslpp::matrix< double > R0_5b
Definition: Charm_Kpnunu.h:198
gslpp::matrix< double > RGevol_J(orders order, int nf, int contribution)
This method returns the J matrices useful for the RGE of the charm penguin/box diagrams.
gslpp::vector< double > ThresholdCb(orders order)
gslpp::vector< double > CW2p
Definition: Charm_Kpnunu.h:200
double C_Bt_qed(orders_qed order_qed)
const StandardModel & model
Definition: Charm_Kpnunu.h:191
gslpp::vector< double > ThresholdCp(orders order)
gslpp::matrix< double > RGevolB(int nf, orders order)
gslpp::vector< double > CW0b
Definition: Charm_Kpnunu.h:197
gslpp::matrix< double > RGevolP(int nf, orders order)
gslpp::vector< double > C(orders order, orders_qed order_qed, int contribution)
This method solves the RGE equation: d/dlog(mu) C = gamma_transp C from mu_w to mu_c relevant for the...
double C_Bt(orders order)
gslpp::vector< double > CWin_muw(orders order, orders_qed order_qed, int contribution)
0805.4119v1 and hep-ph/0603079v3. To retrieve QCD CW at MuW just put NO_QED to order_qed and the qcd_...
double xc_mc_qed
Definition: Charm_Kpnunu.h:194
gslpp::matrix< double > J1_5b
Definition: Charm_Kpnunu.h:198
gslpp::vector< double > dcb
Definition: Charm_Kpnunu.h:193
gslpp::matrix< double > J1_4b
Definition: Charm_Kpnunu.h:198
double C_P(orders order)
gslpp::matrix< double > U0_5b
Definition: Charm_Kpnunu.h:198
gslpp::matrix< double > R1_5b
Definition: Charm_Kpnunu.h:198
WilsonCoefficient evoCkpnn
Definition: Charm_Kpnunu.h:192
Charm_Kpnunu(const StandardModel &model_i)
constructor
std::vector< WilsonCoefficient > & EVOCkpnn()
gslpp::matrix< double > J1_4p
Definition: Charm_Kpnunu.h:201
double C_P_qed(orders_qed order_qed)
gslpp::matrix< double > U0_4b
Definition: Charm_Kpnunu.h:198
gslpp::matrix< double > R1_4b
Definition: Charm_Kpnunu.h:198
gslpp::vector< double > CWeb
Definition: Charm_Kpnunu.h:197
gslpp::matrix< double > R0_4p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > U0_5p
Definition: Charm_Kpnunu.h:201
gslpp::matrix< double > U0_4p
Definition: Charm_Kpnunu.h:201
double getMass_scale() const
A get method to access the scale at which the particle mass is defined.
Definition: Particle.h:133
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
const double getMuc() const
A get method to access the threshold between four- and three-flavour theory in GeV.
Definition: QCD.h:591
const double getCF() const
A get method to access the Casimir factor of QCD.
Definition: QCD.h:618
const double Mrun(const double mu, const double m, const quark q, const orders order=FULLNNLO) const
Computes a running quark mass from .
Definition: QCD.cpp:1353
@ BOTTOM
Definition: QCD.h:329
@ CHARM
Definition: QCD.h:326
const double getMut() const
A get method to access the threshold between six- and five-flavour theory in GeV.
Definition: QCD.h:573
@ TAU
Definition: QCD.h:316
const Particle & getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:536
const double getMub() const
A get method to access the threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:582
A model class for the Standard Model.
const double getMuw() const
A get method to retrieve the matching scale around the weak scale.
const Particle & getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
const double sW2_MSbar_Approx() const
The (approximated formula for the) square of the sine of the weak mixing angle in the MSbar scheme,...
const double getMz() const
A get method to access the mass of the boson .
const double Ale(double mu, orders order, bool Nf_thr=true) const
The running electromagnetic coupling in the scheme.
const CKM & getCKM() const
A get method to retrieve the member object of type CKM.
virtual const double getMHl() const
A get method to retrieve the Higgs mass .
virtual const double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
virtual StandardModelMatching & getMatching() const
A get method to access the member reference of type StandardModelMatching.
const double getGF() const
A get method to retrieve the Fermi constant .
const double Als(const double mu, const orders order, const bool Nf_thr, const bool qed_flag) const
The running QCD coupling in the scheme including QED corrections.
const double sW2_ND() const
The square of the sine of the weak mixing angle in the MSbar-ND scheme (w/o decoupling $\alpha\ln(m_t...
virtual const double alphaMz() const
The electromagnetic coupling at the -mass scale, .
void setCoeff(const gslpp::vector< gslpp::complex > &z, orders order_i)
virtual void setMu(double mu)
orders_qed getOrder_qed() const
orders getOrder() const
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:33
orders_qed
An enum type for orders in electroweak.
Definition: OrderScheme.h:56