a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
THDMMatching.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2015 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include "THDMMatching.h"
9 #include "THDM.h"
10 #include <gsl/gsl_integration.h>
11 #include <gsl/gsl_sf_dilog.h>
12 #include <math.h>
13 #include <stdexcept>
14 
15 THDMMatching::THDMMatching(const THDM & THDM_i) :
16 
17  StandardModelMatching(THDM_i),
18  myTHDM(THDM_i),
19  myCKM(3, 3, 0.),
20  mcdbs2(5, NDR, NLO),
21  mcbtaunu(3, NDR, LO),
22  mcbsg(8, NDR, NNLO),
23  mcprimebsg(8, NDR, NNLO),
24  mcgminus2mu(2,NDR,NLO)
25 {}
26 
27 std::vector<WilsonCoefficient>& THDMMatching::CMdbs2() {
28 
29  double Mut = myTHDM.getMut();
30  double xt = x_t(Mut);
31  double GF=myTHDM.getGF();
32  double MW=myTHDM.Mw();
33  gslpp::complex co = GF / 4. / M_PI * MW * myTHDM.getCKM().computelamt_s();
34  double tanb = myTHDM.gettanb();
35  double mHp2=myTHDM.getmHp2();
36  double xHW=mHp2/(MW*MW);
37  double xtH=xt/xHW;
38  double SWH=xtH*((2.0*xHW-8.0)*log(xtH)/((1.0-xHW)*(1.0-xtH)*(1.0-xtH))+6.0*xHW*log(xt)/((1.0-xHW)*(1.0-xt)*(1.0-xt))-(8.0-2.0*xt)/((1.0-xt)*(1.0-xtH)))/(tanb*tanb);
39  double SHH=xtH*((1.0+xtH)/((1.0-xtH)*(1.0-xtH))+2.0*xtH*log(xtH)/((1.0-xtH)*(1.0-xtH)*(1.0-xtH)))/(tanb*tanb*tanb*tanb);
40 
41  vmcds = StandardModelMatching::CMdbs2();
42  mcdbs2.setMu(Mut);
43 
44  switch (mcdbs2.getOrder()) {
45  case NNLO:
46  case NLO:
47  case LO:
48  mcdbs2.setCoeff(0, co * co * xt * (SWH+SHH), LO);
49  break;
50  default:
51  std::stringstream out;
52  out << mcdbs2.getOrder();
53  throw std::runtime_error("THDMMatching::CMdbs2(): order " + out.str() + "not implemented");
54  }
55 
56  vmcds.push_back(mcdbs2);
57  return(vmcds);
58 }
59 
60 std::vector<WilsonCoefficient>& THDMMatching::CMbtaunu(QCD::meson meson_i) {
61 
62  double Muw = myTHDM.getMuw();
63  double GF = myTHDM.getGF();
64  myCKM = myTHDM.getVCKM();
65  double mB = myTHDM.getMesons(meson_i).getMass();
66  double tanb = myTHDM.gettanb();
67  double mHp2=myTHDM.getmHp2();
68 
69  vmcbtaunu = StandardModelMatching::CMbtaunu(meson_i);
70  mcbtaunu.setMu(Muw);
71 
72  switch (mcbtaunu.getOrder()) {
73  case NNLO:
74  case NLO:
75  case LO:
76  if (meson_i == QCD::B_P) mcbtaunu.setCoeff(0, -4.*GF * myCKM(0,2) / sqrt(2.) * mB*mB*tanb*tanb/mHp2, LO);
77  else if (meson_i == QCD::B_C) mcbtaunu.setCoeff(0, -4.*GF * myCKM(1,2) / sqrt(2.) * mB*mB*tanb*tanb/mHp2, LO);
78  throw std::runtime_error("THDMMatching::CMbtaunu(): meson not implemented");
79  break;
80  default:
81  std::stringstream out;
82  out << mcbtaunu.getOrder();
83  throw std::runtime_error("THDMMatching::CMbtaunu(): order " + out.str() + "not implemented");
84  }
85 
86  vmcbtaunu.push_back(mcbtaunu);
87  return(vmcbtaunu);
88 
89 }
90 
91 std::vector<WilsonCoefficient>& THDMMatching::CMbsg()
92 {
93  double Muw = myTHDM.getMuw();
94 
95  double Mut = myTHDM.getMut();
96  double mt = myTHDM.Mrun(Mut, myTHDM.getQuarks(QCD::TOP).getMass_scale(),
97  myTHDM.getQuarks(QCD::TOP).getMass(), FULLNNLO);
98  double mHp=myTHDM.getmHp();
99 
100  double tanb = myTHDM.gettanb();
101 
102  gslpp::complex co = 1.; // (- 4. * GF / sqrt(2)) * SM.computelamt_s(); THIS SHOULD ALREADY BE IMPLEMENTED IN THE OBSERVABLE
103 
104  vmcbsg = StandardModelMatching::CMbsg();
105  mcbsg.setMu(Muw);
106 
107  switch (mcbsg.getOrder()) {
108  case NNLO:
109  for (int j=0; j<8; j++){
110  mcbsg.setCoeff(j, co * myTHDM.Alstilde5(Muw) * myTHDM.Alstilde5(Muw) * setWCbsg(j, tanb, mt, mHp, Muw, NNLO), NNLO);
111  }
112  case NLO:
113  for (int j=0; j<8; j++){
114  mcbsg.setCoeff(j, co * myTHDM.Alstilde5(Muw) * setWCbsg(j, tanb, mt, mHp, Muw, NLO), NLO);
115  }
116  case LO:
117  for (int j=0; j<8; j++){
118  mcbsg.setCoeff(j, co * setWCbsg(j, tanb, mt, mHp, Muw, LO), LO);
119  }
120  break;
121  default:
122  std::stringstream out;
123  out << mcbsg.getOrder();
124  throw std::runtime_error("THDMMatching::CMbsg(): order " + out.str() + "not implemented");
125  }
126 
127  vmcbsg.push_back(mcbsg);
128  return(vmcbsg);
129 }
130 
131 
132  std::vector<WilsonCoefficient>& THDMMatching::CMprimebsg()
133 {
134  double Muw = myTHDM.getMuw();
135 
136  vmcprimebsg = StandardModelMatching::CMprimebsg();
137  mcprimebsg.setMu(Muw);
138 
139  switch (mcprimebsg.getOrder()) {
140  case NNLO:
141  for (int j=0; j<8; j++){
142  mcprimebsg.setCoeff(j, 0., NNLO);//* CHECK ORDER *//
143  }
144  case NLO:
145  for (int j=0; j<8; j++){
146  mcprimebsg.setCoeff(j, 0., NLO);//* CHECK ORDER *//
147  }
148  case LO:
149  for (int j=0; j<8; j++){
150  mcprimebsg.setCoeff(j, 0., LO);
151  }
152  break;
153  default:
154  std::stringstream out;
155  out << mcprimebsg.getOrder();
156  throw std::runtime_error("THDMMatching::CMprimebsg(): order " + out.str() + "not implemented");
157  }
158 
159  vmcprimebsg.push_back(mcprimebsg);
160  return(vmcprimebsg);
161 }
162 
163 
164 
165 /*******************************************************************************
166  * Wilson coefficients calculus, Misiak base for b -> s gamma *
167  * ****************************************************************************/
168 
169 double THDMMatching::setWCbsg(int i, double tan, double mt, double mhp, double mu, orders order)
170 {
171  if ( tanbsg == tan && mtbsg == mt && mhpbsg == mhp && mubsg == mu){
172  switch (order){
173  case NNLO:
174  return ( CWbsgArrayNNLO[i] );
175  case NLO:
176  return ( CWbsgArrayNLO[i] );
177  break;
178  case LO:
179  return ( CWbsgArrayLO[i] );
180  break;
181  default:
182  std::stringstream out;
183  out << order;
184  throw std::runtime_error("order" + out.str() + "not implemeted");
185  }
186  }
187 
188  tanbsg = tan; mtbsg = mt; mhpbsg = mhp; mubsg = mu;
189 
190  double x = mt*mt/mhp/mhp;
191 
192  double x2 = x*x;
193  double x3 = x2*x;
194  double x4 = x3*x;
195  double x5 = x4*x;
196  double xm = 1. - x;
197  double xm2 = xm*xm;
198  double xm3 = xm2*xm;
199  double xm4 = xm3*xm;
200  double xm6 = xm4*xm2;
201  double xm8 = xm4*xm4;
202  double xo = 1. - 1./x;
203  double xo2 = xo*xo;
204  double xo4 = xo*xo2*xo;
205  double xo6 = xo4*xo2;
206  double xo8 = xo4*xo4;
207  double Lx = log(x);
208  double Lx2 = Lx*Lx;
209  double Lx3 = Lx2*Lx;
210  double tan2 = tan*tan;
211  double Li2 = gsl_sf_dilog(1.-1./x);
212 
213  double lstmu = 2. * log(mu/mt);
214 
215  double n70ct = - ( (7. - 5.*x - 8.*x2)/(36.*xm3) + (2.*x - 3.*x2)*Lx/(6.*xm4) ) * x/2.;
216  double n70fr = ( (3.*x - 5.*x2)/(6.*xm2) + (2.*x - 3.*x2)*Lx/(3.*xm3) ) / 2.;
217 
218  double n80ct = - ( (2. + 5.*x - x2)/(12.*xm3) + (x*Lx)/(2.*xm4) ) * x/2.;
219  double n80fr = ( (3.*x - x2)/(2.*xm2) + (x*Lx)/xm3 ) / 2.;
220 
221  double n41ct = (-16.*x + 29.*x2 - 7.*x3)/(36.*xm3) + (-2.*x + 3.*x2)*Lx/(6.*xm4);
222 
223 
224  double n71ct = (797.*x - 5436.*x2 + 7569.*x3 - 1202.*x4)/(486.*xm4) +
225  (36.*x2 - 74.*x3 + 16.*x4)*Li2/(9.*xm4) +
226  ((7.*x - 463.*x2 + 807.*x3 - 63.*x4)*Lx)/(81.*xm3*xm2);
227  double cd71ct = (-31.*x - 18.*x2 + 135.*x3 - 14.*x4)/(27.*xm4) + (-28.*x2 + 46.*x3 + 6.*x4)*Lx/(9.*xm3*xm2);
228  double n71fr = (28.*x - 52.*x2 + 8.*x3)/(3.*xm3) + (-48.*x + 112.*x2 - 32.*x3)*Li2/(9.*xm3) +
229  (66.*x - 128.*x2 + 14.*x3)*Lx/(9.*xm4);
230  double cd71fr = (42.*x - 94.*x2 + 16.*x3)/(9.*xm3) + (32.*x - 56.*x2 - 12.*x3)*Lx/(9.*xm4);
231 
232  double n81ct = (1130.*x - 18153.*x2 + 7650.*x3 - 4451.*x4)/(1296.*xm4) +
233  (30.*x2 - 17.*x3 + 13.*x4)*Li2/(6.*xm4) +
234  (-2.*x - 2155.*x2 + 321.*x3 - 468.*x4)*Lx/(216.*xm3*xm2);
235  double cd81ct = (-38.*x - 261.*x2 + 18.*x3 - 7.*x4)/(36.*xm4) + (-31.*x2 - 17.*x3)*Lx/(6.*xm3*xm2);
236  double n81fr = (143.*x - 44.*x2 + 29.*x3)/(8.*xm3) + (-36.*x + 25.*x2 - 17.*x3)*Li2/(6.*xm3) +
237  (165.*x - 7.*x2 + 34.*x3)*Lx/(12.*xm4);
238  double cd81fr = (81.*x - 16.*x2 + 7.*x3)/(6.*xm3) + (19.*x + 17.*x2)*Lx/(3.*xm4);
239 
240  double n32ct = (10.*x4 + 30.*x2 - 20.*x)/(27.*xm4) * Li2 +
241  (30.*x3 - 66.*x2 - 56.*x)/(81.*xm4) * Lx + (6.*x3 - 187.*x2 + 213.*x)/(-81.*xm3);
242  double cd32ct = (-30.*x2 + 20.*x)/(27.*xm4)*Lx + (-35.*x3 + 145.*x2 - 80.*x)/(-81.*xm3);
243 
244  double n42ct = (515.*x4 - 906.*x3 + 99.*x2 + 182.*x)/(54.*xm4) * Li2 +
245  (1030.*x4 - 2763.*x3 - 15.*x2 + 980.*x)/(-108.*xm3*xm2)*Lx +
246  (-29467.*x4 + 68142.*x3 - 6717.*x2 - 18134.*x)/(1944.*xm4);
247  double cd42ct = (-375.*x3 - 95.*x2 + 182.*x)/(-54.*xm3*xm2)*Lx +
248  (133.*x4 - 108.*x3 + 4023.*x2 - 2320.*x)/(324.*xm4);
249 
250  double cd72ct = -(x * (67930.*x4 - 470095.*x3 + 1358478.*x2 - 700243.*x + 54970.))/(-2187.*xm3*xm2) +
251  (x * (10422.*x4 - 84390.*x3 + 322801.*x2 - 146588.*x + 1435.))/(729.*xm6)*Lx +
252  (2.*x2 * (260.*x3 - 1515.*x2 + 3757.*x - 1446.))/(-27. * xm3*xm2) * Li2;
253  double ce72ct = (x * (-518.*x4 + 3665.*x3 - 17397.*x2 + 3767.*x + 1843.))/(-162.*xm3*xm2) +
254  (x2 * (-63.*x3 + 532.*x2 + 2089.*x - 1118.))/(27.*xm6)*Lx;
255  double cd72fr = -(x * (3790.*x3 - 22511.*x2 + 53614.*x - 21069.))/(81.*xm4) -
256  (2.*x * (-1266.*x3 + 7642.*x2 - 21467.*x + 8179.))/(-81.*xm3*xm2)*Lx +
257  (8.*x * (139.*x3 - 612.*x2 + 1103.*x - 342.))/(27.*xm4) * Li2;
258  double ce72fr = -(x * (284.*x3 - 1435.*x2 + 4304.*x - 1425.))/(27.*xm4) -
259  (2.*x * (63.*x3 - 397.*x2 - 970.*x + 440.))/(-27.*xm3*xm2)*Lx;
260 
261  double cd82ct = -(x * (522347.*x4 - 2423255.*x3 + 2706021.*x2 - 5930609.*x + 148856))/(-11664.*xm3*xm2) +
262  (x * (51948.*x4 - 233781.*x3 + 48634.*x2 - 698693.*x + 2452.))/(1944.*xm6)*Lx +
263  (x2 * (481.*x3 - 1950.*x2 + 1523.*x - 2550.))/(-18.*xm3*xm2) * Li2;
264  double ce82ct = (x * (-259.*x4 + 1117.*x3 + 2925.*x2 + 28411.*x + 2366.))/(-216.*xm3*xm2) -
265  (x2 * (139.*x2 + 2938.*x + 2683.))/(36.*xm6)*Lx;
266  double cd82fr = -(x * (1463.*x3 - 5794.*x2 + 5543.*x - 15036.))/(27.*xm4) -
267  (x * (-1887.*x3 + 7115.*x2 + 2519.*x + 19901.))/(-54.*xm3*xm2)*Lx -
268  (x * (-629.*x3 + 2178.*x2 - 1729.*x + 2196.))/(18.*xm4) * Li2;
269  double ce82fr = -(x * (259.*x3 - 947.*x2 - 251.*x - 5973.))/(36.*xm4)-
270  (x * (139.*x2 + 2134.*x + 1183.))/(-18.*xm3*xm2)*Lx;
271 
272  double n72ct = 0.;
273  if (mhp < 50.)
274  n72ct = -274.2/x5 - 72.13*Lx/x5 + 24.41/x4 - 168.3*Lx/x4 + 79.15/x3 -
275  103.8*Lx/x3 + 47.09/x2 - 38.12*Lx/x2 + 15.35/x - 8.753*Lx/x + 3.970;
276  else if (mhp < mt)
277  n72ct = 1.283 + 0.7158 * xo + 0.4119 * xo2 + 0.2629 * xo*xo2 + 0.1825 * xo4 +
278  0.1347 * xo*xo4 + 0.1040 * xo6 + 0.0831 * xo*xo6 + 0.06804 * xo8 +
279  0.05688 * xo*xo8 + 0.04833 * xo2*xo8 + 0.04163 * xo*xo2*xo8 + 0.03625 * xo4*xo8 +
280  0.03188 * xo*xo4*xo8 + 0.02827 * xo6*xo8 + 0.02525 * xo*xo6*xo8 + 0.02269 * xo8*xo8;
281  else if (mhp < 520.)
282  n72ct = 1.283 - 0.7158 * xm - 0.3039 * xm2 - 0.1549 * xm3 - 0.08625 * xm4 -
283  0.05020 * xm3*xm2 - 0.02970 * xm6 - 0.01740 * xm3*xm4 - 0.009752 * xm8 -
284  0.004877 * xm3*xm6 - 0.001721 * xm2*xm8 + 0.0003378 * xm3*xm8 + 0.001679 * xm4*xm8 +
285  0.002542 * xm*xm4*xm8 + 0.003083 * xm6*xm8 + 0.003404 * xm*xm6*xm8 + 0.003574 * xm8*xm8;
286  else n72ct = -823.0*x5 + 42.30*x5*Lx3 - 412.4*x5*Lx2 - 3362*x5*Lx -
287  1492*x4 - 23.26*x4*Lx3 - 541.4*x4*Lx2 - 2540*x4*Lx -
288  1158*x3 - 34.50*x3*Lx3 - 348.2*x3*Lx2 - 1292*x3*Lx -
289  480.9*x2 - 20.73*x2*Lx3 - 112.4*x2*Lx2 - 396.1*x2*Lx -
290  8.278*x + 0.9225*x*Lx2 + 4.317*x*Lx;
291 
292  double n72fr = 0.;
293  if (mhp < 50.)
294  n72fr = -( 194.3/x5 + 101.1*Lx/x5 - 24.97/x4 + 168.4*Lx/x4 - 78.90/x3 +
295  106.2*Lx/x3 - 49.32/x2 + 38.43*Lx/x2 - 12.91/x + 9.757*Lx/x + 8.088 );
296  else if (mhp < mt)
297  n72fr = -( 12.82 - 1.663 * xo - 0.8852 * xo2 - 0.4827 * xo*xo2 - 0.2976 * xo4 -
298  0.2021 * xo*xo4 - 0.1470 * xo6 - 0.1125 * xo*xo6 - 0.08931 * xo8 -
299  0.07291 * xo*xo8 - 0.06083 * xo2*xo8 - 0.05164 * xo*xo2*xo8 - 0.04446 * xo4*xo8 -
300  0.03873 * xo*xo4*xo8 - 0.03407 * xo6*xo8 - 0.03023 * xo*xo6*xo8 - 0.02702 * xo8*xo8 );
301  else if (mhp < 400.)
302  n72fr = -( 12.82 + 1.663 * xm + 0.7780 * xm2 + 0.3755 * xm3 + 0.1581 * xm4 +
303  0.03021 * xm3*xm2 - 0.04868 * xm6 - 0.09864 * xm3*xm4 - 0.1306 * xm8 -
304  0.1510 * xm3*xm6 - 0.1637 * xm2*xm8 - 0.1712 * xm3*xm8 - 0.1751 * xm4*xm8 -
305  0.1766 * xm*xm4*xm8 - 0.1763 * xm6*xm8 - 0.1748 * xm*xm6*xm8 - 0.1724 * xm8*xm8 );
306  else n72fr = -( 2828 * x5 - 66.63 * x5*Lx3 + 469.4 * x5*Lx2 + 1986 * x5*Lx +
307  1480 * x4 + 36.08 * x4*Lx3 + 323.2 * x4*Lx2 + 169.9 * x4*Lx +
308  166.7 * x3 + 19.73 * x3*Lx3 - 46.61 * x3*Lx2 - 826.2 * x3*Lx -
309  524.1 * x2 - 8.889 * x2*Lx3 - 195.7 * x2*Lx2 - 870.3 * x2*Lx -
310  572.2 * x - 20.94 * x*Lx3 - 123.5 * x*Lx2 - 453.5 * x*Lx );
311 
312  double n82ct = 0.;
313  if (mhp < 50.)
314  n82ct = 826.2/x5 - 300.7*Lx/x5 + 96.35/x4 + 91.89*Lx/x4 - 66.39/x3 +
315  78.58*Lx/x3 - 39.76/x2 + 20.02*Lx/x2 - 5.214/x + 2.278;
316  else if (mhp < mt)
317  n82ct = 1.188 + 0.4078 * xo + 0.2002 * xo2 + 0.1190 * xo*xo2 + 0.07861 * xo4 +
318  0.05531 * xo*xo4 + 0.04061 * xo6 + 0.03075 * xo*xo6 + 0.02386 * xo8 +
319  0.01888 * xo*xo8 + 0.01520 * xo2*xo8 + 0.01241 * xo*xo2*xo8 + 0.01026 * xo4*xo8 +
320  0.008575 * xo*xo4*xo8 + 0.007238 * xo6*xo8 + 0.006164 * xo*xo6*xo8 + 0.005290 * xo8*xo8;
321  else if (mhp < 600.)
322  n82ct = 1.188 - 0.4078 * xm - 0.2076 * xm2 - 0.1265 * xm3 - 0.08570 * xm4 -
323  0.06204 * xm3*xm2 - 0.04689 * xm6 - 0.03652 * xm3*xm4 - 0.02907 * xm8 -
324  0.02354 * xm3*xm6 - 0.01933 * xm2*xm8 - 0.01605 * xm3*xm8 - 0.01345 * xm4*xm8 -
325  0.01137 * xm*xm4*xm8 - 0.009678 * xm6*xm8 - 0.008293 * xm*xm6*xm8 - 0.007148 * xm8*xm8;
326  else n82ct = -19606 * x5 - 226.7 * x5*Lx3 - 5251 * x5*Lx2 - 26090 * x5*Lx -
327  9016 * x4 - 143.4 * x4*Lx3 - 2244 * x4*Lx2 - 10102 * x4*Lx -
328  3357 * x3 - 66.32 * x3*Lx3 - 779.6 * x3*Lx2 - 3077 * x3*Lx -
329  805.5 * x2 - 22.98 * x2*Lx3 - 169.1 * x2*Lx2 - 602.7 * x2*Lx +
330  0.7437 * x + 0.6908 * x*Lx2 + 3.238 * x*Lx;
331 
332  double n82fr = 0.;
333  if (mhp < 50.)
334  n82fr = -( -1003/x5 + 476.9*Lx/x5 - 205.7/x4 - 71.62*Lx/x4 + 62.26/x3 -
335  110.7*Lx/x3 + 63.74/x2 - 35.42*Lx/x2 + 10.89/x - 3.174 );
336  else if (mhp < mt)
337  n82fr = -( -0.6110 - 1.095 * xo - 0.4463 * xo2 - 0.2568 * xo*xo2 - 0.1698 * xo4 -
338  0.1197 * xo*xo4 - 0.08761 * xo6 - 0.06595 * xo*xo6 - 0.05079 * xo8 -
339  0.03987 * xo*xo8 - 0.03182 * xo2*xo8 - 0.02577 * xo*xo2*xo8 - 0.02114 * xo4*xo8 -
340  0.01754 * xo*xo4*xo8 - 0.01471 * xo6*xo8 - 0.01244 * xo*xo6*xo8 - 0.01062 * xo8*xo8 );
341  else if (mhp < 520.)
342  n82fr = -( -0.6110 + 1.095 * xm + 0.6492 * xm2 + 0.4596 * xm3 + 0.3569 * xm4 +
343  0.2910 * xm3*xm2 + 0.2438 * xm6 + 0.2075 * xm3*xm4 + 0.1785 * xm8 +
344  0.1546 * xm3*xm6 + 0.1347 * xm2*xm8 + 0.1177 * xm3*xm8 + 0.1032 * xm4*xm8 +
345  0.09073 * xm*xm4*xm8 + 0.07987 * xm6*xm8 + 0.07040 * xm*xm6*xm8 + 0.06210 * xm8*xm8 );
346  else n82fr = -( -15961 * x5 + 1003 * x5*Lx3 - 2627 * x5*Lx2 - 29962 * x5*Lx -
347  11683 * x4 + 54.66 * x4*Lx3 - 2777 * x4*Lx2 - 17770 * x4*Lx -
348  6481 * x3 - 40.68 * x3*Lx3 - 1439 * x3*Lx2 - 7906 * x3*Lx -
349  2943 * x2 - 31.83 * x2*Lx3 - 612.6 * x2*Lx2 - 2770 * x2*Lx -
350  929.8 * x - 19.80 * x*Lx3 - 174.7 * x*Lx2 - 658.4 * x*Lx );
351 
352 
353  switch (order){
354  case NNLO:
355  CWbsgArrayNNLO[2] = ( n32ct + cd32ct * lstmu ) / tan2;
356  CWbsgArrayNNLO[3] = ( n42ct + cd42ct * lstmu ) / tan2;
357  CWbsgArrayNNLO[4] = 2./15. * n41ct / tan2 - CWbsgArrayNNLO[2] / 10.;
358  CWbsgArrayNNLO[5] = 1./4. * n41ct / tan2 - CWbsgArrayNNLO[2] * 3./16.;
359  CWbsgArrayNNLO[6] = ( n72ct + cd72ct * lstmu + ce72ct * lstmu * lstmu)/tan2 +
360  n72fr + cd72fr * lstmu + ce72fr * lstmu * lstmu
361  - 1./3.*CWbsgArrayNNLO[2] - 4./9.*CWbsgArrayNNLO[3] - 20./3.*CWbsgArrayNNLO[4] - 80./9.*CWbsgArrayNNLO[5];
362  CWbsgArrayNNLO[7] = ( n82ct + cd82ct * lstmu + ce82ct * lstmu * lstmu)/tan2 +
363  n82fr + cd82fr * lstmu + ce82fr * lstmu * lstmu
364  + CWbsgArrayNNLO[2] - 1./6.*CWbsgArrayNNLO[3] - 20.*CWbsgArrayNNLO[4] - 10./3.*CWbsgArrayNNLO[5];
365  case NLO:
366  CWbsgArrayNLO[3] = n41ct / tan2;
367  CWbsgArrayNLO[6] = ( n71ct + cd71ct * lstmu ) / tan2 + n71fr + cd71fr * lstmu - 4./9.*CWbsgArrayNLO[3];
368  CWbsgArrayNLO[7] = ( n81ct + cd81ct * lstmu ) / tan2 + n81fr + cd81fr * lstmu - 1./6.*CWbsgArrayNLO[3];
369  case LO:
370  CWbsgArrayLO[6] = n70ct /tan2 + n70fr ;
371  CWbsgArrayLO[7] = n80ct /tan2 + n80fr ;
372  break;
373  default:
374  std::stringstream out;
375  out << order;
376  throw std::runtime_error("order" + out.str() + "not implemeted");
377  }
378 
379  /*std::cout << "CWbsgArrayLO[6] = " << CWbsgArrayLO[6] << std::endl;
380  std::cout << "CWbsgArrayLO[7] = " << CWbsgArrayLO[7] << std::endl << std::endl;
381  std::cout << "CWbsgArrayNLO[3] = " << CWbsgArrayNLO[3] << std::endl;
382  std::cout << "CWbsgArrayNLO[6] = " << CWbsgArrayNLO[6] << std::endl;
383  std::cout << "CWbsgArrayNLO[7] = " << CWbsgArrayNLO[7] << std::endl << std::endl;
384  std::cout << "CWbsgArrayNNLO[2] = " << CWbsgArrayNNLO[2] << std::endl;
385  std::cout << "CWbsgArrayNNLO[3] = " << CWbsgArrayNNLO[3] << std::endl;
386  std::cout << "CWbsgArrayNNLO[4] = " << CWbsgArrayNNLO[4] << std::endl;
387  std::cout << "CWbsgArrayNNLO[5] = " << CWbsgArrayNNLO[5] << std::endl;
388  std::cout << "CWbsgArrayNNLO[6] = " << CWbsgArrayNNLO[6] << std::endl;
389  std::cout << "CWbsgArrayNNLO[7] = " << CWbsgArrayNNLO[7] << std::endl << std::endl;*/
390 
391 
392  switch (order){
393  case NNLO:
394  return ( CWbsgArrayNNLO[i] );
395  case NLO:
396  return ( CWbsgArrayNLO[i] );
397  break;
398  case LO:
399  return ( CWbsgArrayLO[i] );
400  break;
401  default:
402  std::stringstream out;
403  out << order;
404  throw std::runtime_error("order" + out.str() + "not implemeted");
405  }
406 }
407 
408 
409 /*******************************************************************************
410  * Muon g-2 (according to arxiv:1409.3199 *
411  * ****************************************************************************/
412 
413 double THDMMatching::gminus2muLO() {
414 
415  double gminus2muLO;
416 
417  double pi=M_PI;
418  double GF=myTHDM.getGF();
419  double mMU=myTHDM.getLeptons(StandardModel::MU).getMass();
420  std::string modelflag=myTHDM.getModelTypeflag();
421  double mHl2=myTHDM.getmHl2();
422  double mHh2=myTHDM.getmHh2();
423  double mA2=myTHDM.getmA2();
424  double mHp2=myTHDM.getmHp2();
425  double tanb=myTHDM.gettanb();
426  double sinb=tanb/sqrt(1.0+tanb*tanb);
427  double cosb=1.0/sqrt(1.0+tanb*tanb);
428  double bma=myTHDM.getbma();
429  double sina=sinb*cos(bma)-cosb*sin(bma);
430  double cosa=cosb*cos(bma)+sinb*sin(bma);
431 
432  double ymu_h, ymu_H, ymu_A;
433  double rmu_hSM, rmu_h, rmu_H, rmu_A, rmu_Hp;
434  double part_hSM, part_h, part_H, part_A, part_Hp;
435 
436  if( modelflag == "type1" ) {
437  ymu_h=cosa/sinb;
438  ymu_H=sina/sinb;
439  ymu_A=-1.0/tanb;
440  }
441  else if( modelflag == "type2" ) {
442  ymu_h=-sina/cosb;
443  ymu_H=cosa/cosb;
444  ymu_A=tanb;
445  }
446  else if( modelflag == "typeX" ) {
447  ymu_h=-sina/cosb;
448  ymu_H=cosa/cosb;
449  ymu_A=tanb;
450  }
451  else if( modelflag == "typeY" ) {
452  ymu_h=cosa/sinb;
453  ymu_H=sina/sinb;
454  ymu_A=-1.0/tanb;
455  }
456  else {
457  throw std::runtime_error("modelflag can be only any of \"type1\", \"type2\", \"typeX\" or \"typeY\"");
458  }
459 
460  if( mHl2<1.0 || mHh2<1.0 || mA2<1.0 || mHp2<1.0)
461  {
462  throw std::runtime_error("The implemented approximation for g-2_\\mu only works for Higgs masses above 1 GeV.");
463  }
464  rmu_hSM=mMU*mMU/15647.5081;
465  rmu_h=mMU*mMU/mHl2;
466  rmu_H=mMU*mMU/mHh2;
467  rmu_A=mMU*mMU/mA2;
468  rmu_Hp=mMU*mMU/mHp2;
469 
470  //https://arxiv.org/pdf/1409.3199.pdf
471 
472  part_hSM=rmu_hSM*(-7.0/6.0-log(rmu_hSM));
473  part_h=ymu_h*ymu_h*rmu_h*(-7.0/6.0-log(rmu_h));
474  part_H=ymu_H*ymu_H*rmu_H*(-7.0/6.0-log(rmu_H));
475  part_A=ymu_A*ymu_A*rmu_A*(11.0/6.0+log(rmu_A));
476  part_Hp=-ymu_A*ymu_A*rmu_Hp*1.0/6.0;
477 
478  gminus2muLO=GF*mMU*mMU/(4.0*pi*pi*sqrt(2.0)) * (-part_hSM+part_h+part_H+part_A+part_Hp);
479 
480  return(gminus2muLO);
481 }
482 
483 struct __f_params{
484  double a;
485 };
486 
487 double __fS_integrand(double x, void* p){
488  __f_params &params= *reinterpret_cast<__f_params *>(p);
489  double integ = (2.0*x*(1.0-x)-1.0) * log(x*(1.0-x)/params.a) / (x*(1.0-x)-params.a);
490  return integ;
491 }
492 
493 double __fPS_integrand(double x, void* p){
494  __f_params &params= *reinterpret_cast<__f_params *>(p);
495  double integ = log(x*(1.0-x)/params.a) / (x*(1.0-x)-params.a);
496  return integ;
497 }
498 
499 double gscalar(double r)
500 {
501  __f_params params;
502  params.a=r;
503 
504  double result;
505  gsl_integration_glfixed_table * w
506  = gsl_integration_glfixed_table_alloc(100);
507  gsl_function F;
508 
509  F.function = &__fS_integrand;
510  F.params = reinterpret_cast<void *>(&params);
511 
512  result = gsl_integration_glfixed (&F, 0, 1, w);
513 
514  gsl_integration_glfixed_table_free (w);
515 
516  return result;
517 }
518 
519 
520 double gpseudoscalar(double r)
521 {
522  __f_params params;
523  params.a=r;
524 
525  double result;
526  gsl_integration_glfixed_table * w
527  = gsl_integration_glfixed_table_alloc(100);
528  gsl_function F;
529 
530  F.function = &__fPS_integrand;
531  F.params = reinterpret_cast<void *>(&params);
532 
533  result = gsl_integration_glfixed (&F, 0, 1, w);
534 
535  gsl_integration_glfixed_table_free (w);
536 
537  return result;
538 }
539 
540 double THDMMatching::gminus2muNLO() {
541 
542  double gminus2muNLO;
543 
544  double pi=M_PI;
545  double GF=myTHDM.getGF();
546  double aem=myTHDM.getAle();
547  double mMU=myTHDM.getLeptons(StandardModel::MU).getMass();
548  double mTAU=myTHDM.getLeptons(StandardModel::TAU).getMass();
549  double mt=myTHDM.getQuarks(QCD::TOP).getMass();
550  double mb=myTHDM.getQuarks(QCD::BOTTOM).getMass();
551  std::string modelflag=myTHDM.getModelTypeflag();
552  double mHl2=myTHDM.getmHl2();
553  double mHh2=myTHDM.getmHh2();
554  double mA2=myTHDM.getmA2();
555  double mHp2=myTHDM.getmHp2();
556  double tanb=myTHDM.gettanb();
557  double sinb=tanb/sqrt(1.0+tanb*tanb);
558  double cosb=1.0/sqrt(1.0+tanb*tanb);
559  double bma=myTHDM.getbma();
560  double sina=sinb*cos(bma)-cosb*sin(bma);
561  double cosa=cosb*cos(bma)+sinb*sin(bma);
562 
563  double ymu_h, ymu_H, ymu_A, yb_h, yb_H, yb_A, yt_h, yt_H, yt_A;
564  double rtau_hSM, rtau_h, rtau_H, rtau_A, rt_hSM, rt_h, rt_H, rt_A, rb_hSM, rb_h, rb_H, rb_A;
565  double part_hSM, part_h, part_H, part_A;
566 
567  yt_h=cosa/sinb;
568  yt_H=sina/sinb;
569  yt_A=1.0/tanb;
570 
571  if( modelflag == "type1" ) {
572  ymu_h=cosa/sinb;
573  ymu_H=sina/sinb;
574  ymu_A=-1.0/tanb;
575  yb_h=ymu_h;
576  yb_H=ymu_H;
577  yb_A=ymu_A;
578  }
579  else if( modelflag == "type2" ) {
580  ymu_h=-sina/cosb;
581  ymu_H=cosa/cosb;
582  ymu_A=tanb;
583  yb_h=ymu_h;
584  yb_H=ymu_H;
585  yb_A=ymu_A;
586  }
587  else if( modelflag == "typeX" ) {
588  ymu_h=-sina/cosb;
589  ymu_H=cosa/cosb;
590  ymu_A=tanb;
591  yb_h=yt_h;
592  yb_H=yt_H;
593  yb_A=-yt_A;
594  }
595  else if( modelflag == "typeY" ) {
596  ymu_h=cosa/sinb;
597  ymu_H=sina/sinb;
598  ymu_A=-1.0/tanb;
599  yb_h=-sina/cosb;
600  yb_H=cosa/cosb;
601  yb_A=tanb;
602  }
603  else {
604  throw std::runtime_error("modelflag can be only any of \"type1\", \"type2\", \"typeX\" or \"typeY\"");
605  }
606 
607  if( mHl2<mMU*mMU || mHh2<mMU*mMU || mA2<mMU*mMU || mHp2<mMU*mMU)
608  {
609  throw std::runtime_error("The implemented two-loop approximation for g-2_\\mu only works for Higgs masses above m_mu^2.");
610  }
611  rtau_hSM=mTAU*mTAU/15647.5081;
612  rtau_h=mTAU*mTAU/mHl2;
613  rtau_H=mTAU*mTAU/mHh2;
614  rtau_A=mTAU*mTAU/mA2;
615  rt_hSM=mt*mt/15647.5081;
616  rt_h=mt*mt/mHl2;
617  rt_H=mt*mt/mHh2;
618  rt_A=mt*mt/mA2;
619  rb_hSM=mb*mb/15647.5081;
620  rb_h=mb*mb/mHl2;
621  rb_H=mb*mb/mHh2;
622  rb_A=mb*mb/mA2;
623 
624  part_hSM=rtau_hSM*gscalar(rtau_hSM)+(4.0/3.0)*rt_hSM*gscalar(rt_hSM)+(1.0/3.0)*rb_hSM*gscalar(rb_hSM);
625  part_h=ymu_h*(ymu_h*rtau_h*gscalar(rtau_h)+(4.0/3.0)*yt_h*rt_h*gscalar(rt_h)+(1.0/3.0)*yb_h*rb_h*gscalar(rb_h));
626  part_H=ymu_H*(ymu_H*rtau_H*gscalar(rtau_H)+(4.0/3.0)*yt_H*rt_H*gscalar(rt_H)+(1.0/3.0)*yb_H*rb_H*gscalar(rb_H));
627  part_A=ymu_A*(ymu_A*rtau_A*gpseudoscalar(rtau_A)+(4.0/3.0)*yt_A*rt_A*gpseudoscalar(rt_A)+(1.0/3.0)*yb_A*rb_A*gpseudoscalar(rb_A));
628 
629  gminus2muNLO=GF*mMU*mMU/(4.0*pi*pi*pi*sqrt(2.0)) * aem * (-part_hSM+part_h+part_H+part_A);
630 
631  return(gminus2muNLO);
632 }
633 
634 std::vector<WilsonCoefficient>& THDMMatching::CMgminus2mu() {
635 
636  vmcgminus2mu = StandardModelMatching::CMgminus2mu();
637 
638  double gminus2muLOvalue=gminus2muLO();
639  double gminus2muNLOvalue=gminus2muNLO();
640 
641  switch (mcgminus2mu.getOrder()) {
642  case LO:
643  mcgminus2mu.setCoeff(0, gminus2muLOvalue, LO); //g-2_muR
644  mcgminus2mu.setCoeff(1, 0., LO); //g-2_muL
645  break;
646  case NLO:
647  mcgminus2mu.setCoeff(0, gminus2muLOvalue+gminus2muNLOvalue, NLO); //g-2_muR
648  mcgminus2mu.setCoeff(1, 0., NLO); //g-2_muL
649  break;
650  case NNLO:
651  default:
652  std::stringstream out;
653  out << mcgminus2mu.getOrder();
654  throw std::runtime_error("THDMMatching::CMgminus2mu(): order " + out.str() + " not implemented.\nOnly leading order (LO) or next-to-leading order (NLO) are allowed.");
655  }
656 
657  vmcgminus2mu.push_back(mcgminus2mu);
658  return(vmcgminus2mu);
659 
660 }
QCD::TAU
Definition: QCD.h:316
__fS_integrand
double __fS_integrand(double x, void *p)
Definition: THDMMatching.cpp:485
gslpp::cos
complex cos(const complex &z)
Definition: gslpp_complex.cpp:429
__f_params
Definition: THDMMatching.cpp:481
THDMMatching::mcgminus2mu
WilsonCoefficient mcgminus2mu
Definition: THDMMatching.h:89
QCD::BOTTOM
Definition: QCD.h:329
Li2
cd Li2(cd x)
Definition: hpl.h:1011
THDMMatching::mcprimebsg
WilsonCoefficient mcprimebsg
Definition: THDMMatching.h:89
THDMMatching::gminus2muLO
virtual double gminus2muLO()
Calculates amplitudes for at one loop from .
Definition: THDMMatching.cpp:411
__f_params::a
double a
Definition: THDMMatching.cpp:482
THDMMatching::setWCbsg
double setWCbsg(int i, double tan, double mt, double mhp, double mu, orders order)
Definition: THDMMatching.cpp:168
StandardModelMatching::CMdbs2
virtual std::vector< WilsonCoefficient > & CMdbs2()
,
Definition: StandardModelMatching.cpp:1049
THDM
A base class for symmetric Two-Higgs-Doublet models.
Definition: THDM.h:120
CKM::computelamt_s
gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:136
gslpp::sin
complex sin(const complex &z)
Definition: gslpp_complex.cpp:420
THDMMatching::CMgminus2mu
virtual std::vector< WilsonCoefficient > & CMgminus2mu()
Wilson coefficient for .
Definition: THDMMatching.cpp:632
THDM::getmA2
double getmA2() const
A method get the squared mass of the pseudoscalar Higgs A.
Definition: THDM.h:423
LO
Definition: OrderScheme.h:33
THDMMatching::mcbtaunu
WilsonCoefficient mcbtaunu
Definition: THDMMatching.h:89
THDM::gettanb
double gettanb() const
A method get .
Definition: THDM.h:283
THDMMatching::CMbtaunu
virtual std::vector< WilsonCoefficient > & CMbtaunu(QCD::meson meson_i)
Definition: THDMMatching.cpp:60
NDR
Definition: OrderScheme.h:21
ModelMatching::CMbsg
virtual std::vector< WilsonCoefficient > & CMbsg()=0
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
THDM::getmHp2
double getmHp2() const
A method get the squared charged Higgs mass.
Definition: THDM.h:457
Particle::getMass_scale
double getMass_scale() const
A get method to access the scale at which the particle mass is defined.
Definition: Particle.h:133
THDMMatching.h
THDMMatching::tanbsg
double tanbsg
Definition: THDMMatching.h:92
THDMMatching::CMprimebsg
virtual std::vector< WilsonCoefficient > & CMprimebsg()
operator basis: current current; qcd penguins; magnetic and chromomagnetic penguins; semileptonic
Definition: THDMMatching.cpp:132
THDM::getmHh2
double getmHh2() const
A method get the squared mass of the "non-125 GeV" neutral scalar Higgs.
Definition: THDM.h:365
THDMMatching::CWbsgArrayNNLO
double CWbsgArrayNNLO[8]
Definition: THDMMatching.h:91
StandardModel::Alstilde5
double Alstilde5(const double mu) const
The value of at any scale with the number of flavours and full EW corrections.
Definition: StandardModel.cpp:899
THDMMatching::mcdbs2
WilsonCoefficient mcdbs2
Definition: THDMMatching.h:89
THDMMatching::myTHDM
const THDM & myTHDM
Definition: THDMMatching.h:87
__fPS_integrand
double __fPS_integrand(double x, void *p)
Definition: THDMMatching.cpp:491
WilsonTemplate::getOrder
orders getOrder() const
Definition: WilsonTemplate.h:65
StandardModelMatching
A class for the matching in the Standard Model.
Definition: StandardModelMatching.h:26
WilsonCoefficient::setCoeff
void setCoeff(const gslpp::vector< gslpp::complex > &z, orders order_i)
Definition: WilsonCoefficient.h:34
THDMMatching::mtbsg
double mtbsg
Definition: THDMMatching.h:92
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
gslpp::tan
complex tan(const complex &z)
Definition: gslpp_complex.cpp:438
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
QCD::B_P
Definition: QCD.h:344
QCD::TOP
Definition: QCD.h:328
StandardModel::getGF
double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:739
THDMMatching::THDMMatching
THDMMatching(const THDM &THDM_i)
Definition: THDMMatching.cpp:15
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
QCD::getMesons
Meson getMesons(const QCD::meson m) const
A get method to access a meson as an object of the type Meson.
Definition: QCD.h:524
QCD::getMut
double getMut() const
A get method to access the threshold between six- and five-flavour theory in GeV.
Definition: QCD.h:561
THDM::getbma
double getbma() const
A method get .
Definition: THDM.h:307
THDMMatching::CWbsgArrayNLO
double CWbsgArrayNLO[8]
Definition: THDMMatching.h:91
THDMMatching::mcbsg
WilsonCoefficient mcbsg
Definition: THDMMatching.h:89
NNLO
Definition: OrderScheme.h:35
QCD::B_C
Definition: QCD.h:346
QCD::getQuarks
Particle getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:534
WilsonTemplate::setMu
virtual void setMu(double mu)
Definition: WilsonTemplate.h:92
StandardModel::getVCKM
gslpp::matrix< gslpp::complex > getVCKM() const
A get method to retrieve the CKM matrix.
Definition: StandardModel.h:870
THDM::getModelTypeflag
std::string getModelTypeflag() const
A method get the THDM model type.
Definition: THDM.h:242
THDMMatching::mubsg
double mubsg
Definition: THDMMatching.h:92
ModelMatching::CMprimebsg
virtual std::vector< WilsonCoefficient > & CMprimebsg()=0
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
gscalar
double gscalar(double r)
Definition: THDMMatching.cpp:497
THDMMatching::CWbsgArrayLO
double CWbsgArrayLO[8]
Definition: THDMMatching.h:91
THDM::getmHl2
double getmHl2() const
A method get the squared mass of the lighter neutral scalar Higgs.
Definition: THDM.h:339
QCD::Mrun
double Mrun(const double mu, const double m, const orders order=FULLNNLO) const
Computes a running quark mass from .
Definition: QCD.cpp:1064
gpseudoscalar
double gpseudoscalar(double r)
Definition: THDMMatching.cpp:518
NLO
Definition: OrderScheme.h:34
THDM.h
StandardModel::Mw
virtual double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel.cpp:970
THDMMatching::mhpbsg
double mhpbsg
Definition: THDMMatching.h:92
THDMMatching::gminus2muNLO
virtual double gminus2muNLO()
Calculates amplitudes for at approximate two-loop from .
Definition: THDMMatching.cpp:538
THDMMatching::CMdbs2
virtual std::vector< WilsonCoefficient > & CMdbs2()
Definition: THDMMatching.cpp:27
FULLNNLO
Definition: OrderScheme.h:38
THDM::getmHp
double getmHp() const
A method get the charged Higgs mass.
Definition: THDM.h:471
StandardModel::getAle
double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:748
THDMMatching::myCKM
gslpp::matrix< gslpp::complex > myCKM
Definition: THDMMatching.h:88
StandardModel::getCKM
CKM getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:879
StandardModel::getMuw
double getMuw() const
A get method to retrieve the matching scale around the weak scale.
Definition: StandardModel.h:938
QCD::MU
Definition: QCD.h:314
THDMMatching::CMbsg
virtual std::vector< WilsonCoefficient > & CMbsg()
operator basis: current current; qcd penguins; magnetic and chromomagnetic penguins; semileptonic
Definition: THDMMatching.cpp:91
StandardModel::getLeptons
Particle getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:712