a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EWSM_Output.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2013 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include "EWSM_Output.h"
9 #include "EWSMcache.h"
10 #include "EWSMOneLoopEW.h"
11 #include "EWSMTwoLoopQCD.h"
12 #include "EWSMThreeLoopQCD.h"
13 #include "EWSMTwoLoopEW.h"
14 #include "EWSMThreeLoopEW2QCD.h"
15 #include "EWSMThreeLoopEW.h"
16 #include "EWSMApproximateFormulae.h"
17 #include "EWSMOneLoopEW_HV.h"
18 #include "EWSMTwoFermionsLEP2.h"
19 
20 EWSM_Output::EWSM_Output(const StandardModel& SM_in)
21 : StandardModel(SM_in)
22 {
23 }
24 
25 void EWSM_Output::outputEachDeltaR(const double Mw_i) const
26 {
27  std::cout << "Mw_SM = " << Mw() << std::endl;
28  std::cout << "DeltaR_SM() = " << DeltaR() << std::endl;
29  std::cout << "DeltaRbar_SM() = " << DeltaRbar() << std::endl;
30  std::cout << "Mw(input) = " << Mw_i << std::endl;
31 
32  double cW2_TMP = Mw_i * Mw_i / getMz() / getMz();
33  double sW2_TMP = 1.0 - cW2_TMP;
34 
35  double DeltaRho[StandardModel::orders_EW_size];
36  DeltaRho[StandardModel::EW1] = getMyOneLoopEW()->DeltaRho(Mw_i);
37  DeltaRho[StandardModel::EW1QCD1] = getMyTwoLoopQCD()->DeltaRho(Mw_i);
38  DeltaRho[StandardModel::EW1QCD2] = getMyThreeLoopQCD()->DeltaRho(Mw_i);
39  DeltaRho[StandardModel::EW2] = getMyTwoLoopEW()->DeltaRho(Mw_i);
40  DeltaRho[StandardModel::EW2QCD1] = getMyThreeLoopEW2QCD()->DeltaRho(Mw_i);
41  DeltaRho[StandardModel::EW3] = getMyThreeLoopEW()->DeltaRho(Mw_i);
42 
43  double DeltaR_rem[StandardModel::orders_EW_size];
44  DeltaR_rem[StandardModel::EW1] = getMyOneLoopEW()->DeltaR_rem(Mw_i);
45  DeltaR_rem[StandardModel::EW1QCD1] = getMyTwoLoopQCD()->DeltaR_rem(Mw_i);
46  DeltaR_rem[StandardModel::EW1QCD2] = getMyThreeLoopQCD()->DeltaR_rem(Mw_i);
47  DeltaR_rem[StandardModel::EW2] = getMyTwoLoopEW()->DeltaR_rem(Mw_i);
48  DeltaR_rem[StandardModel::EW2QCD1] = getMyThreeLoopEW2QCD()->DeltaR_rem(Mw_i);
49  DeltaR_rem[StandardModel::EW3] = getMyThreeLoopEW()->DeltaR_rem(Mw_i);
50 
51  double f_AlphaToGF = sqrt(2.0) * getGF() * pow(getMz(), 2.0) * sW2_TMP * cW2_TMP / M_PI / getAle();
52  //f_AlphaToGF = 1.0; /* for test */
53  double DeltaRho_sum = f_AlphaToGF * DeltaRho[StandardModel::EW1]
54  + f_AlphaToGF * DeltaRho[StandardModel::EW1QCD1]
55  + f_AlphaToGF * DeltaRho[StandardModel::EW1QCD2]
56  + pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2]
57  + pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2QCD1]
58  + pow(f_AlphaToGF, 3.0) * DeltaRho[StandardModel::EW3];
59  double DeltaRho_G = f_AlphaToGF * DeltaRho[StandardModel::EW1];
60 
61  if (getFlagMw().compare("NORESUM") == 0) {
62 
63  f_AlphaToGF = 1.0;
64  DeltaRho[StandardModel::EW1QCD2] *= f_AlphaToGF;
65  DeltaRho[StandardModel::EW2QCD1] *= pow(f_AlphaToGF, 2.0);
66  DeltaRho[StandardModel::EW3] *= pow(f_AlphaToGF, 3.0);
67 
68  // Full EW one-loop contribution (without the full DeltaAlphaL5q)
69  double DeltaR_EW1 = -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW1] + DeltaR_rem[StandardModel::EW1];
70 
71  // Full EW two-loop contribution with reducible corrections
72  double DeltaR_EW2_rem = getMyApproximateFormulae()->DeltaR_TwoLoopEW_rem(Mw_i);
73 
74  // EW two-loop irreducible contributions with large-mt expansion
75  double DeltaR_EW2_old_red = DeltaAlphaL5q() * DeltaAlphaL5q()
76  - 2.0 * cW2_TMP / sW2_TMP * DeltaAlphaL5q() * DeltaRho[StandardModel::EW1]
77  + pow(cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW1], 2.0);
78  double DeltaR_EW2_old_irred = -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW2] + DeltaR_rem[StandardModel::EW2];
79 
80  // Delta r, including the full EW two-loop contribution
81  double deltaR = DeltaAlphaL5q();
82  for (int j = 0; j < StandardModel::orders_EW_size; ++j) {
83  deltaR += -cW2_TMP / sW2_TMP * DeltaRho[(StandardModel::orders_EW)j];
84  deltaR += DeltaR_rem[(StandardModel::orders_EW)j];
85  }
86  deltaR -= -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW2];
87  deltaR -= DeltaR_rem[StandardModel::EW2];
88  deltaR += DeltaAlphaL5q() * DeltaAlphaL5q() + 2.0 * DeltaAlphaL5q() * DeltaR_EW1 + DeltaR_EW2_rem;
89 
90  std::cout << "(1+dr) - 1 = " << deltaR << std::endl;
91  std::cout << " EW1 = " << DeltaAlphaL5q() + DeltaR_EW1 << std::endl;
92  std::cout << " DeltaAlphaL5q = " << DeltaAlphaL5q() << std::endl;
93  std::cout << " dR = " << DeltaR_EW1 << std::endl;
94  std::cout << " EW1QCD1 = " << -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW1QCD1] + DeltaR_rem[StandardModel::EW1QCD1] << std::endl;
95  std::cout << " EW2(full) = " << DeltaR_EW2_rem + DeltaAlphaL5q() * DeltaAlphaL5q() + 2.0 * DeltaAlphaL5q() * DeltaR_EW1 << std::endl;
96  std::cout << " dAle*dAle = " << DeltaAlphaL5q() * DeltaAlphaL5q() << std::endl;
97  std::cout << " 2*dAle*dR = " << 2.0 * DeltaAlphaL5q() * DeltaR_EW1 << std::endl;
98  std::cout << " others = " << DeltaR_EW2_rem << std::endl;
99  std::cout << " EW1QCD2 = " << -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW1QCD2] + DeltaR_rem[StandardModel::EW1QCD2] << std::endl;
100  std::cout << " EW2QCD1 = " << -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW2QCD1] + DeltaR_rem[StandardModel::EW2QCD1] << std::endl;
101  std::cout << " EW3 = " << -cW2_TMP / sW2_TMP * DeltaRho[StandardModel::EW3] + DeltaR_rem[StandardModel::EW3] << std::endl;
102  std::cout << " EW2(old,irreducible) = " << DeltaR_EW2_old_irred << std::endl;
103  std::cout << " EW2(old,red+irred) = " << DeltaR_EW2_old_red + DeltaR_EW2_old_irred << std::endl;
104  std::cout << " EW2(old,red+irred-dAle*dAle-2*dAle*dR) = "
105  << DeltaR_EW2_old_red + DeltaR_EW2_old_irred
106  - DeltaAlphaL5q() * DeltaAlphaL5q()
107  - 2.0 * DeltaAlphaL5q() * DeltaR_EW1 << std::endl;
108 
109  } else if (getFlagMw().compare("OMSI") == 0) {
110 
111  // R = 1/(1 - Delta r)
112  double R = 1.0 / (1.0 + cW2_TMP / sW2_TMP * DeltaRho_sum)
113  / (1.0 - DeltaAlphaL5q() - DeltaR_rem[StandardModel::EW1] - DeltaR_rem[StandardModel::EW1QCD1] - DeltaR_rem[StandardModel::EW2]);
114 
115  std::cout << "1/(1-dr) - 1 (exact) = " << R - 1.0 << std::endl;
116  std::cout << " --> dr = " << 1.0 - 1.0 / R << std::endl;
117 
118  // each contribution
119  double DeltaR_EW1 = DeltaAlphaL5q() - cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1] + DeltaR_rem[StandardModel::EW1];
120  double DeltaR_EW1QCD1 = -cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1QCD1] + DeltaR_rem[StandardModel::EW1QCD1];
121  double DeltaR_EW2 = -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2]
122  + DeltaR_rem[StandardModel::EW2]
123  + cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1]
124  *(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1])
125  + DeltaR_EW1*DeltaR_EW1;
126  double DeltaR_EW1QCD2 = -cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1QCD2];
127  double DeltaR_EW2QCD1 = -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2QCD1]
128  + cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1QCD1]
129  *(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1])
130  + 2.0 * DeltaR_EW1*DeltaR_EW1QCD1;
131  double DeltaR_EW3 = -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 3.0) * DeltaRho[StandardModel::EW3]
132  + cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2]
133  *(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1])
134  + pow(DeltaR_EW1, 3.0)
135  + 2.0 * DeltaR_EW1 * (DeltaR_EW2 - DeltaR_EW1 * DeltaR_EW1);
136 
137  std::cout << " EW1 = " << DeltaR_EW1 << std::endl;
138  std::cout << " DeltaAlphaL5q = " << DeltaAlphaL5q() << std::endl;
139  std::cout << " -cW2/sW2*dRho1= " << -cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1] << std::endl;
140  std::cout << " DeltaR1_rem = " << DeltaR_rem[StandardModel::EW1] << std::endl;
141  std::cout << " EW1QCD1 = " << DeltaR_EW1QCD1 << std::endl;
142  std::cout << " EW2(full) = " << DeltaR_EW2 << std::endl;
143  std::cout << " EW1*EW1 = " << DeltaR_EW1 * DeltaR_EW1 << std::endl;
144  std::cout << " dAle*dAle = " << DeltaAlphaL5q() * DeltaAlphaL5q() << std::endl;
145  std::cout << " others = " << DeltaR_EW1 * DeltaR_EW1 - DeltaAlphaL5q() * DeltaAlphaL5q() << std::endl;
146  std::cout << " -cW2/sW2*dRho2= " << -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2] << std::endl;
147  std::cout << " DeltaR2_rem = " << DeltaR_rem[StandardModel::EW2] << std::endl;
148  std::cout << " others = " << cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1]*(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1]) << std::endl;
149  std::cout << " EW1QCD2 = " << DeltaR_EW1QCD2 << std::endl;
150  std::cout << " EW2QCD1 = " << DeltaR_EW2QCD1 << std::endl;
151  std::cout << " EW3 = " << DeltaR_EW3 << std::endl;
152  std::cout << " -cW2/sW2*dRho3= " << -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 3.0) * DeltaRho[StandardModel::EW3] << std::endl;
153  std::cout << " EW1^3 = " << pow(DeltaR_EW1, 3.0) << std::endl;
154  std::cout << " 2*EW1*(EW2-EW1^2)=" << 2.0 * DeltaR_EW1 * (DeltaR_EW2 - DeltaR_EW1 * DeltaR_EW1) << std::endl;
155  std::cout << " others = " << cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2]*(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1]) << std::endl;
156 
157  } else if (getFlagMw().compare("OMSII") == 0) {
158 
159  // R = 1/(1 - Delta r)
160  double R = 1.0 / ((1.0 + cW2_TMP / sW2_TMP * DeltaRho_sum)*(1.0 - DeltaAlphaL5q())
161  - (1.0 + cW2_TMP / sW2_TMP * DeltaRho_G) * DeltaR_rem[StandardModel::EW1]
162  - DeltaR_rem[StandardModel::EW1QCD1] - DeltaR_rem[StandardModel::EW2]);
163 
164  // each contribution
165  double DeltaR_EW1 = DeltaAlphaL5q() - cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1] + DeltaR_rem[StandardModel::EW1];
166  double DeltaR_EW1QCD1 = -cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1QCD1] + DeltaR_rem[StandardModel::EW1QCD1];
167  double DeltaR_EW2 = -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2]
168  + DeltaR_rem[StandardModel::EW2]
169  + cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1]
170  *(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1])
171  + DeltaR_EW1*DeltaR_EW1;
172  double DeltaR_EW1QCD2 = -cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1QCD2];
173  double DeltaR_EW2QCD1 = -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2QCD1]
174  + cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1QCD1] * DeltaAlphaL5q()
175  + 2.0 * DeltaR_EW1*DeltaR_EW1QCD1;
176  double DeltaR_EW3 = -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 3.0) * DeltaRho[StandardModel::EW3]
177  + cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2] * DeltaAlphaL5q()
178  + pow(DeltaR_EW1, 3.0)
179  + 2.0 * DeltaR_EW1 * (DeltaR_EW2 - DeltaR_EW1 * DeltaR_EW1);
180 
181  std::cout << "1/(1-dr) - 1 (exact) = " << R - 1.0 << std::endl;
182  std::cout << " --> dr = " << 1.0 - 1.0 / R << std::endl;
183  std::cout << "1/(1-dr) - 1 (sum of expanded terms) = "
184  << DeltaR_EW1 + DeltaR_EW1QCD1 + DeltaR_EW2 + DeltaR_EW1QCD2
185  + DeltaR_EW2QCD1 + DeltaR_EW3 << std::endl;
186  std::cout << " EW1 = " << DeltaR_EW1 << std::endl;
187  std::cout << " DeltaAlphaL5q = " << DeltaAlphaL5q() << std::endl;
188  std::cout << " -cW2/sW2*dRho1= " << -cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1] << std::endl;
189  std::cout << " DeltaR1_rem = " << DeltaR_rem[StandardModel::EW1] << std::endl;
190  std::cout << " EW1QCD1 = " << DeltaR_EW1QCD1 << std::endl;
191  std::cout << " EW2(full) = " << DeltaR_EW2 << std::endl;
192  std::cout << " EW1*EW1 = " << DeltaR_EW1 * DeltaR_EW1 << std::endl;
193  std::cout << " dAle*dAle = " << DeltaAlphaL5q() * DeltaAlphaL5q() << std::endl;
194  std::cout << " others = " << DeltaR_EW1 * DeltaR_EW1 - DeltaAlphaL5q() * DeltaAlphaL5q() << std::endl;
195  std::cout << " -cW2/sW2*dRho2= " << -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2] << std::endl;
196  std::cout << " DeltaR2_rem = " << DeltaR_rem[StandardModel::EW2] << std::endl;
197  std::cout << " others = " << cW2_TMP / sW2_TMP * f_AlphaToGF * DeltaRho[StandardModel::EW1]*(DeltaAlphaL5q() + DeltaR_rem[StandardModel::EW1]) << std::endl;
198  std::cout << " EW1QCD2 = " << DeltaR_EW1QCD2 << std::endl;
199  std::cout << " EW2QCD1 = " << DeltaR_EW2QCD1 << std::endl;
200  std::cout << " EW3 = " << DeltaR_EW3 << std::endl;
201  std::cout << " -cW2/sW2*dRho3= " << -cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 3.0) * DeltaRho[StandardModel::EW3] << std::endl;
202  std::cout << " EW1^3 = " << pow(DeltaR_EW1, 3.0) << std::endl;
203  std::cout << " 2*EW1*(EW2-EW1^2)=" << 2.0 * DeltaR_EW1 * (DeltaR_EW2 - DeltaR_EW1 * DeltaR_EW1) << std::endl;
204  std::cout << " others = " << cW2_TMP / sW2_TMP * pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2] * DeltaAlphaL5q() << std::endl;
205 
206  } else
207  std::cout << "EWSM_Output::outputEachDeltaR(): Not implemented for schemeMw="
208  << getFlagMw() << std::endl;
209 }
210 
211 void EWSM_Output::outputEachDeltaRhoZ_l(const QCD::lepton l, const double Mw_i) const
212 {
213  std::cout << "================================================" << std::endl;
214  std::cout << "rhoZ_l[(QCD::lepton)" << l << "]" << std::endl;
215  std::cout << "Mw(input) = " << Mw_i << std::endl;
216 
217  double cW2_TMP = Mw_i * Mw_i / getMz() / getMz();
218  double sW2_TMP = 1.0 - cW2_TMP;
219 
220  double DeltaRho[StandardModel::orders_EW_size];
221  DeltaRho[StandardModel::EW1] = getMyOneLoopEW()->DeltaRho(Mw_i);
222  DeltaRho[StandardModel::EW1QCD1] = getMyTwoLoopQCD()->DeltaRho(Mw_i);
223  DeltaRho[StandardModel::EW1QCD2] = getMyThreeLoopQCD()->DeltaRho(Mw_i);
224  DeltaRho[StandardModel::EW2] = getMyTwoLoopEW()->DeltaRho(Mw_i);
225  DeltaRho[StandardModel::EW2QCD1] = getMyThreeLoopEW2QCD()->DeltaRho(Mw_i);
226  DeltaRho[StandardModel::EW3] = getMyThreeLoopEW()->DeltaRho(Mw_i);
227 
228  /* compute delta rho_rem^f */
229  gslpp::complex deltaRho_rem_f[StandardModel::orders_EW_size];
230  deltaRho_rem_f[StandardModel::EW1] = getMyOneLoopEW()->deltaRho_rem_f(getLeptons(l), Mw_i);
231 #ifdef WITHIMTWOLOOPQCD
232  deltaRho_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaRho_rem_f(getLeptons(l), Mw_i).real(),
233  getMyTwoLoopQCD()->deltaRho_rem_f(getLeptons(l), Mw_i).imag(), false);
234 #else
235  deltaRho_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaRho_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
236 #endif
237  deltaRho_rem_f[StandardModel::EW1QCD2] = gslpp::complex(getMyThreeLoopQCD()->deltaRho_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
238  deltaRho_rem_f[StandardModel::EW2] = gslpp::complex(getMyTwoLoopEW()->deltaRho_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
239  deltaRho_rem_f[StandardModel::EW2QCD1] = gslpp::complex(getMyThreeLoopEW2QCD()->deltaRho_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
240  deltaRho_rem_f[StandardModel::EW3] = gslpp::complex(getMyThreeLoopEW()->deltaRho_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
241 
242  /* compute Delta rbar_rem */
243  double DeltaRbar_rem = getMyOneLoopEW()->DeltaRbar_rem(Mw_i);
244 
245  double f_AlphaToGF = sqrt(2.0) * getGF() * pow(getMz(), 2.0)
246  * sW2_TMP * cW2_TMP / M_PI / getAle();
247 
248  /* Re[rho_Z^f] with or without resummation */
249  double deltaRho_rem_f_real[StandardModel::orders_EW_size];
250  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
251  deltaRho_rem_f_real[j] = deltaRho_rem_f[j].real();
252 
253  double dummy[StandardModel::orders_EW_size];
254  outputEachDeltaRhoZ(f_AlphaToGF, DeltaRho, deltaRho_rem_f_real,
255  DeltaRbar_rem, false, dummy, 0.0);
256 
257  /* Im[rho_Z^f] without resummation */
258  double ImRhoZf = 0.0;
259  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
260  ImRhoZf += deltaRho_rem_f[j].imag();
261  std::cout << "ImRhoZf(with alpha)=" << ImRhoZf << std::endl;
262 
263  std::cout << "================================================" << std::endl;
264 }
265 
266 void EWSM_Output::outputEachDeltaRhoZ_q(const QCD::quark q, const double Mw_i) const
267 {
268  std::cout << "================================================" << std::endl;
269  std::cout << "rhoZ_q[(QCD::quark)" << q << "]" << std::endl;
270  std::cout << "Mw(input) = " << Mw_i << std::endl;
271 
272  double cW2_TMP = Mw_i * Mw_i / getMz() / getMz();
273  double sW2_TMP = 1.0 - cW2_TMP;
274 
275  double DeltaRho[StandardModel::orders_EW_size];
276  DeltaRho[StandardModel::EW1] = getMyOneLoopEW()->DeltaRho(Mw_i);
277  DeltaRho[StandardModel::EW1QCD1] = getMyTwoLoopQCD()->DeltaRho(Mw_i);
278  DeltaRho[StandardModel::EW1QCD2] = getMyThreeLoopQCD()->DeltaRho(Mw_i);
279  DeltaRho[StandardModel::EW2] = getMyTwoLoopEW()->DeltaRho(Mw_i);
280  DeltaRho[StandardModel::EW2QCD1] = getMyThreeLoopEW2QCD()->DeltaRho(Mw_i);
281  DeltaRho[StandardModel::EW3] = getMyThreeLoopEW()->DeltaRho(Mw_i);
282 
283  /* compute delta rho_rem^f */
284  gslpp::complex deltaRho_rem_f[StandardModel::orders_EW_size];
285  deltaRho_rem_f[StandardModel::EW1] = getMyOneLoopEW()->deltaRho_rem_f(getQuarks(q), Mw_i);
286 #ifdef WITHIMTWOLOOPQCD
287  deltaRho_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaRho_rem_f(getQuarks(q), Mw_i).real(),
288  getMyTwoLoopQCD()->deltaRho_rem_f(getQuarks(q), Mw_i).imag(), false);
289 #else
290  deltaRho_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaRho_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
291 #endif
292  deltaRho_rem_f[StandardModel::EW1QCD2] = gslpp::complex(getMyThreeLoopQCD()->deltaRho_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
293  deltaRho_rem_f[StandardModel::EW2] = gslpp::complex(getMyTwoLoopEW()->deltaRho_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
294  deltaRho_rem_f[StandardModel::EW2QCD1] = gslpp::complex(getMyThreeLoopEW2QCD()->deltaRho_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
295  deltaRho_rem_f[StandardModel::EW3] = gslpp::complex(getMyThreeLoopEW()->deltaRho_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
296 
297  /* compute Delta rbar_rem */
298  double DeltaRbar_rem = getMyOneLoopEW()->DeltaRbar_rem(Mw_i);
299 
300  /* conversion factor */
301  double f_AlphaToGF = sqrt(2.0) * getGF() * pow(getMz(), 2.0)
302  * sW2_TMP * cW2_TMP / M_PI / getAle();
303 
304  /* Zbb */
305  bool bool_Zbb = false;
306  if (q == QCD::BOTTOM) bool_Zbb = true;
307  double ZbbSubtract = 0.0;
308  if (bool_Zbb)
309  ZbbSubtract = -getAle() / 4.0 / M_PI / sW2_TMP
310  * pow(getMtpole() / Mw_i, 2.0);
311  double taub[StandardModel::orders_EW_size];
312  double Xt = getMyEWSMcache()->Xt_alpha(Mw_i);
313  if (bool_Zbb) {
314  taub[StandardModel::EW1] = -2.0 * Xt;
315  taub[StandardModel::EW1QCD1] = 2.0 / 3.0 * M_PI * Xt * getMyEWSMcache()->alsMt();
316  taub[StandardModel::EW2] = -2.0 * Xt * Xt * getMyTwoLoopEW()->tau_2();
317  }
318 
319  /* Re[rho_Z^f] with or without resummation */
320  double deltaRho_rem_f_real[StandardModel::orders_EW_size];
321  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
322  deltaRho_rem_f_real[j] = deltaRho_rem_f[j].real();
323 
324  outputEachDeltaRhoZ(f_AlphaToGF, DeltaRho, deltaRho_rem_f_real,
325  DeltaRbar_rem, bool_Zbb, taub, ZbbSubtract);
326 
327  /* Im[rho_Z^f] without resummation */
328  double ImRhoZf = 0.0;
329  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
330  ImRhoZf += deltaRho_rem_f[j].imag();
331  std::cout << "ImRhoZf(with alpha)=" << ImRhoZf << std::endl;
332 
333  std::cout << "================================================" << std::endl;
334 }
335 
336 void EWSM_Output::outputEachDeltaRhoZ(const double f_AlphaToGF,
337  const double DeltaRho[StandardModel::orders_EW_size],
338  const double deltaRho_rem[StandardModel::orders_EW_size],
339  const double DeltaRbar_rem,
340  const bool bool_Zbb,
341  const double taub[StandardModel::orders_EW_size],
342  const double ZbbSubtract) const
343 {
344  if (getFlagRhoZ().compare("APPROXIMATEFORMULA") == 0) {
345 
346  } else if (getFlagRhoZ().compare("NORESUM") == 0) {
347  std::cout << "Leading contributions: alpha or Gmu" << std::endl;
348  std::cout << " DeltaRho[EW1]=" << DeltaRho[StandardModel::EW1] << " "
349  << f_AlphaToGF * DeltaRho[StandardModel::EW1] << std::endl;
350  std::cout << " DeltaRho[EW1QCD1]=" << DeltaRho[StandardModel::EW1QCD1] << " "
351  << f_AlphaToGF * DeltaRho[StandardModel::EW1QCD1] << std::endl;
352  std::cout << " DeltaRho[EW1QCD2]=" << DeltaRho[StandardModel::EW1QCD2] << " "
353  << f_AlphaToGF * DeltaRho[StandardModel::EW1QCD2] << std::endl;
354  std::cout << " DeltaRho[EW2]=" << DeltaRho[StandardModel::EW2] << " "
355  << pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2] << std::endl;
356  std::cout << " DeltaRho[EW2QCD1]=" << DeltaRho[StandardModel::EW2QCD1] << " "
357  << pow(f_AlphaToGF, 2.0) * DeltaRho[StandardModel::EW2QCD1] << std::endl;
358  std::cout << " DeltaRho[EW3]=" << DeltaRho[StandardModel::EW3] << " "
359  << pow(f_AlphaToGF, 3.0) * DeltaRho[StandardModel::EW3] << std::endl;
360  std::cout << "Remainder contributions: alpha or Gmu" << std::endl;
361  std::cout << " DeltaRbar_rem[EW1]=" << DeltaRbar_rem << " "
362  << f_AlphaToGF * DeltaRbar_rem << std::endl;
363  std::cout << " deltaRho_rem[EW1]=" << deltaRho_rem[StandardModel::EW1] - ZbbSubtract << " "
364  << f_AlphaToGF * (deltaRho_rem[StandardModel::EW1] - ZbbSubtract) << std::endl;
365  std::cout << " deltaRho_rem[EW1QCD1]=" << deltaRho_rem[StandardModel::EW1QCD1] << " "
366  << f_AlphaToGF * deltaRho_rem[StandardModel::EW1QCD1] << std::endl;
367  std::cout << " deltaRho_rem[EW2]=" << deltaRho_rem[StandardModel::EW2] << " "
368  << pow(f_AlphaToGF, 2.0) * deltaRho_rem[StandardModel::EW2] << std::endl;
369  if (bool_Zbb) {
370  std::cout << "Taub: alpha or Gmu" << std::endl;
371  std::cout << " taub[EW1]=" << taub[StandardModel::EW1] << " "
372  << f_AlphaToGF * taub[StandardModel::EW1] << std::endl;
373  std::cout << " taub[EW1QCD1]=" << taub[StandardModel::EW1QCD1] << " "
374  << f_AlphaToGF * taub[StandardModel::EW1QCD1] << std::endl;
375  std::cout << " taub[EW2]=" << taub[StandardModel::EW2] << " "
376  << pow(f_AlphaToGF, 2.0) * taub[StandardModel::EW2] << std::endl;
377  }
378  std::cout << "Each order: alpha or Gmu" << std::endl;
379  double dRho_EW1 = DeltaRho[StandardModel::EW1] + deltaRho_rem[StandardModel::EW1] - ZbbSubtract;
380  double dRho_EW1QCD1 = DeltaRho[StandardModel::EW1QCD1] + deltaRho_rem[StandardModel::EW1QCD1];
381  double dRho_EW2 = DeltaRho[StandardModel::EW1] * DeltaRho[StandardModel::EW1]
382  + (deltaRho_rem[StandardModel::EW1] - ZbbSubtract) * DeltaRho[StandardModel::EW1]
383  - DeltaRho[StandardModel::EW1] * DeltaRbar_rem
384  + DeltaRho[StandardModel::EW2]
385  + deltaRho_rem[StandardModel::EW2];
386  double dRho_EW1QCD2 = DeltaRho[StandardModel::EW1QCD2];
387  double dRho_EW2QCD1 = DeltaRho[StandardModel::EW2QCD1] + deltaRho_rem[StandardModel::EW1QCD1] * DeltaRho[StandardModel::EW1];
388  double dRho_EW3 = DeltaRho[StandardModel::EW3];
389  //
390  double dRho_EW1_TMP = dRho_EW1;
391  double dRho_EW1QCD1_TMP = dRho_EW1QCD1;
392  double dRho_EW2_TMP = dRho_EW2;
393  double dRho_EW1QCD2_TMP = dRho_EW1QCD2;
394  double dRho_EW2QCD1_TMP = dRho_EW2QCD1;
395  double dRho_EW3_TMP = dRho_EW3;
396  //
397  if (bool_Zbb) {
398  dRho_EW1 = dRho_EW1_TMP + 2.0 * taub[StandardModel::EW1];
399  dRho_EW1QCD1 = dRho_EW1QCD1_TMP + 2.0 * taub[StandardModel::EW1QCD1];
400  dRho_EW2 = dRho_EW2_TMP + 2.0 * taub[StandardModel::EW2] + taub[StandardModel::EW1] * taub[StandardModel::EW1]
401  + dRho_EW1_TMP * 2.0 * taub[StandardModel::EW1];
402  dRho_EW1QCD2 = dRho_EW1QCD2_TMP;
403  dRho_EW2QCD1 = dRho_EW2QCD1_TMP + dRho_EW1_TMP * 2.0 * taub[StandardModel::EW1QCD1]
404  + dRho_EW1QCD1_TMP * 2.0 * taub[StandardModel::EW1] + 2.0 * taub[StandardModel::EW1] * taub[StandardModel::EW1QCD1];
405  dRho_EW3 = dRho_EW3_TMP + dRho_EW1_TMP * 2.0 * taub[StandardModel::EW2]
406  + dRho_EW2_TMP * 2.0 * taub[StandardModel::EW1] + 2.0 * taub[StandardModel::EW2] * taub[StandardModel::EW1];
407  }
408  std::cout << " EW1: " << dRho_EW1 << " " << f_AlphaToGF * dRho_EW1 << std::endl;
409  std::cout << " EW1QCD1: " << dRho_EW1QCD1 << " " << f_AlphaToGF * dRho_EW1QCD1 << std::endl;
410  std::cout << " EW2: " << dRho_EW2 << " " << pow(f_AlphaToGF, 2.0) * dRho_EW2 << std::endl;
411  std::cout << " EW1QCD2: " << dRho_EW1QCD2 << " " << f_AlphaToGF * dRho_EW1QCD2 << std::endl;
412  std::cout << " EW2QCD1: " << dRho_EW2QCD1 << " " << pow(f_AlphaToGF, 2.0) * dRho_EW2QCD1 << std::endl;
413  std::cout << " EW3: " << dRho_EW3 << " " << pow(f_AlphaToGF, 3.0) * dRho_EW3 << std::endl;
414  std::cout << "Total contribution: alpha or Gmu" << std::endl;
415  std::cout << " rhoZ="
416  << 1.0 + dRho_EW1 + dRho_EW1QCD1 + dRho_EW2
417  + dRho_EW1QCD2 + dRho_EW2QCD1 + dRho_EW3
418  << " "
419  << 1.0 + f_AlphaToGF * dRho_EW1 + f_AlphaToGF * dRho_EW1QCD1
420  + pow(f_AlphaToGF, 2.0) * dRho_EW2
421  + f_AlphaToGF * dRho_EW1QCD2 + pow(f_AlphaToGF, 2.0) * dRho_EW2QCD1
422  + pow(f_AlphaToGF, 3.0) * dRho_EW3
423  << std::endl;
424  if (bool_Zbb) {
425  std::cout << " rhoZ(taub resummed)="
426  << (1.0 + dRho_EW1_TMP + dRho_EW1QCD1_TMP + dRho_EW2_TMP
427  + dRho_EW1QCD2_TMP + dRho_EW2QCD1_TMP + dRho_EW3_TMP)
428  * pow(1.0 + taub[StandardModel::EW1] + taub[StandardModel::EW1QCD1] + taub[StandardModel::EW2], 2.0)
429  << " "
430  << (1.0 + f_AlphaToGF * dRho_EW1_TMP + f_AlphaToGF * dRho_EW1QCD1_TMP
431  + pow(f_AlphaToGF, 2.0) * dRho_EW2_TMP
432  + f_AlphaToGF * dRho_EW1QCD2_TMP + pow(f_AlphaToGF, 2.0) * dRho_EW2QCD1_TMP
433  + pow(f_AlphaToGF, 3.0) * dRho_EW3_TMP)
434  * pow(1.0 + f_AlphaToGF * taub[StandardModel::EW1] + f_AlphaToGF * taub[StandardModel::EW1QCD1]
435  + pow(f_AlphaToGF, 2.0) * taub[StandardModel::EW2], 2.0)
436  << std::endl;
437  }
438  } else
439  std::cout << "EWSM_Output::outputEachDeltaRhoZ(): Not implemented for schemeRhoZ="
440  << getFlagRhoZ() << std::endl;
441 }
442 
443 void EWSM_Output::outputEachDeltaKappaZ_l(const QCD::lepton l, const double Mw_i) const
444 {
445  std::cout << "================================================" << std::endl;
446  std::cout << "kappaZ_l[(QCD::lepton)" << l << "]" << std::endl;
447  std::cout << "Mw(input) = " << Mw_i << std::endl;
448 
449  double cW2_TMP = Mw_i * Mw_i / getMz() / getMz();
450  double sW2_TMP = 1.0 - cW2_TMP;
451 
452  double DeltaRho[StandardModel::orders_EW_size];
453  DeltaRho[StandardModel::EW1] = getMyOneLoopEW()->DeltaRho(Mw_i);
454  DeltaRho[StandardModel::EW1QCD1] = getMyTwoLoopQCD()->DeltaRho(Mw_i);
455  DeltaRho[StandardModel::EW1QCD2] = getMyThreeLoopQCD()->DeltaRho(Mw_i);
456  DeltaRho[StandardModel::EW2] = getMyTwoLoopEW()->DeltaRho(Mw_i);
457  DeltaRho[StandardModel::EW2QCD1] = getMyThreeLoopEW2QCD()->DeltaRho(Mw_i);
458  DeltaRho[StandardModel::EW3] = getMyThreeLoopEW()->DeltaRho(Mw_i);
459 
460  /* compute delta kappa_rem^f */
461  gslpp::complex deltaKappa_rem_f[StandardModel::orders_EW_size];
462  deltaKappa_rem_f[StandardModel::EW1] = getMyOneLoopEW()->deltaKappa_rem_f(getLeptons(l), Mw_i);
463 #ifdef WITHIMTWOLOOPQCD
464  deltaKappa_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaKappa_rem_f(getLeptons(l), Mw_i).real(),
465  getMyTwoLoopQCD()->deltaKappa_rem_f(getLeptons(l), Mw_i).imag(), false);
466 #else
467  deltaKappa_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaKappa_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
468 #endif
469  deltaKappa_rem_f[StandardModel::EW1QCD2] = gslpp::complex(getMyThreeLoopQCD()->deltaKappa_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
470  deltaKappa_rem_f[StandardModel::EW2] = gslpp::complex(getMyTwoLoopEW()->deltaKappa_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
471  deltaKappa_rem_f[StandardModel::EW2QCD1] = gslpp::complex(getMyThreeLoopEW2QCD()->deltaKappa_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
472  deltaKappa_rem_f[StandardModel::EW3] = gslpp::complex(getMyThreeLoopEW()->deltaKappa_rem_f(getLeptons(l), Mw_i).real(), 0.0, false);
473 
474  /* compute Delta rbar_rem */
475  double DeltaRbar_rem = getMyOneLoopEW()->DeltaRbar_rem(Mw_i);
476 
477  /* conversion factor */
478  double f_AlphaToGF = sqrt(2.0) * getGF() * pow(getMz(), 2.0)
479  * sW2_TMP * cW2_TMP / M_PI / getAle();
480 
481  /* Re[Kappa_Z^f] with or without resummation */
482  double deltaKappa_rem_f_real[StandardModel::orders_EW_size];
483  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
484  deltaKappa_rem_f_real[j] = deltaKappa_rem_f[j].real();
485 
486  /* O(alpha^2) correction to Re[kappa_Z^f] from the Z-gamma mixing */
487  double ReKappaZf = resumKappaZ(DeltaRho, deltaKappa_rem_f_real,
488  DeltaRbar_rem, false);
489  double Zgamma_EW2 = 35.0 * alphaMz() * alphaMz() / 18.0 / sW2_TMP
490  * (1.0 - 8.0 / 3.0 * ReKappaZf * sW2_TMP);
491 
492  double dummy[StandardModel::orders_EW_size];
493  outputEachDeltaKappaZ(f_AlphaToGF, cW2_TMP / sW2_TMP,
494  DeltaRho, deltaKappa_rem_f_real,
495  DeltaRbar_rem, false, dummy, 0.0, Zgamma_EW2);
496 
497  /* Im[kappa_Z^f] without resummation */
498  double ImKappaZf = 0.0;
499  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
500  ImKappaZf += deltaKappa_rem_f[j].imag();
501  std::cout << "ImKappaZf(with alpha)=" << ImKappaZf << std::endl;
502 
503  std::cout << "================================================" << std::endl;
504 }
505 
506 void EWSM_Output::outputEachDeltaKappaZ_q(const QCD::quark q, const double Mw_i) const
507 {
508  std::cout << "================================================" << std::endl;
509  std::cout << "kappaZ_q[(QCD::quark)" << q << "]" << std::endl;
510  std::cout << "Mw(input) = " << Mw_i << std::endl;
511 
512  double cW2_TMP = Mw_i * Mw_i / getMz() / getMz();
513  double sW2_TMP = 1.0 - cW2_TMP;
514 
515  double DeltaRho[StandardModel::orders_EW_size];
516  DeltaRho[StandardModel::EW1] = getMyOneLoopEW()->DeltaRho(Mw_i);
517  DeltaRho[StandardModel::EW1QCD1] = getMyTwoLoopQCD()->DeltaRho(Mw_i);
518  DeltaRho[StandardModel::EW1QCD2] = getMyThreeLoopQCD()->DeltaRho(Mw_i);
519  DeltaRho[StandardModel::EW2] = getMyTwoLoopEW()->DeltaRho(Mw_i);
520  DeltaRho[StandardModel::EW2QCD1] = getMyThreeLoopEW2QCD()->DeltaRho(Mw_i);
521  DeltaRho[StandardModel::EW3] = getMyThreeLoopEW()->DeltaRho(Mw_i);
522 
523  /* compute delta kappa_rem^f */
524  gslpp::complex deltaKappa_rem_f[StandardModel::orders_EW_size];
525  deltaKappa_rem_f[StandardModel::EW1] = getMyOneLoopEW()->deltaKappa_rem_f(getQuarks(q), Mw_i);
526 #ifdef WITHIMTWOLOOPQCD
527  deltaKappa_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaKappa_rem_f(getQuarks(q), Mw_i).real(),
528  getMyTwoLoopQCD()->deltaKappa_rem_f(getQuarks(q), Mw_i).imag(), false);
529 #else
530  deltaKappa_rem_f[StandardModel::EW1QCD1] = gslpp::complex(getMyTwoLoopQCD()->deltaKappa_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
531 #endif
532  deltaKappa_rem_f[StandardModel::EW1QCD2] = gslpp::complex(getMyThreeLoopQCD()->deltaKappa_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
533  deltaKappa_rem_f[StandardModel::EW2] = gslpp::complex(getMyTwoLoopEW()->deltaKappa_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
534  deltaKappa_rem_f[StandardModel::EW2QCD1] = gslpp::complex(getMyThreeLoopEW2QCD()->deltaKappa_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
535  deltaKappa_rem_f[StandardModel::EW3] = gslpp::complex(getMyThreeLoopEW()->deltaKappa_rem_f(getQuarks(q), Mw_i).real(), 0.0, false);
536 
537  /* compute Delta rbar_rem */
538  double DeltaRbar_rem = getMyOneLoopEW()->DeltaRbar_rem(Mw_i);
539 
540  /* conversion factor */
541  double f_AlphaToGF = sqrt(2.0) * getGF() * pow(getMz(), 2.0)
542  * sW2_TMP * cW2_TMP / M_PI / getAle();
543 
544  /* Zbb */
545  bool bool_Zbb = false;
546  if (q == QCD::BOTTOM) bool_Zbb = true;
547  double ZbbSubtract = 0.0;
548  if (bool_Zbb)
549  ZbbSubtract = getAle() / 8.0 / M_PI / sW2_TMP
550  * pow(getMtpole() / Mw_i, 2.0);
551  double taub[StandardModel::orders_EW_size];
552  double Xt = getMyEWSMcache()->Xt_alpha(Mw_i);
553  if (bool_Zbb) {
554  taub[StandardModel::EW1] = -2.0 * Xt;
555  taub[StandardModel::EW1QCD1] = 2.0 / 3.0 * M_PI * Xt * getMyEWSMcache()->alsMt();
556  taub[StandardModel::EW2] = -2.0 * Xt * Xt * getMyTwoLoopEW()->tau_2();
557  }
558 
559  /* Re[Kappa_Z^f] with or without resummation */
560  double deltaKappa_rem_f_real[StandardModel::orders_EW_size];
561  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
562  deltaKappa_rem_f_real[j] = deltaKappa_rem_f[j].real();
563 
564  /* O(alpha^2) correction to Re[kappa_Z^f] from the Z-gamma mixing */
565  double ReKappaZf = resumKappaZ(DeltaRho, deltaKappa_rem_f_real,
566  DeltaRbar_rem, bool_Zbb);
567  double Zgamma_EW2 = 35.0 * alphaMz() * alphaMz() / 18.0 / sW2_TMP
568  * (1.0 - 8.0 / 3.0 * ReKappaZf * sW2_TMP);
569 
570  outputEachDeltaKappaZ(f_AlphaToGF, cW2_TMP / sW2_TMP,
571  DeltaRho, deltaKappa_rem_f_real,
572  DeltaRbar_rem, bool_Zbb, taub, ZbbSubtract, Zgamma_EW2);
573 
574  /* Im[kappa_Z^f] without resummation */
575  double ImKappaZf = 0.0;
576  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
577  ImKappaZf += deltaKappa_rem_f[j].imag();
578  std::cout << "ImKappaZf(with alpha)=" << ImKappaZf << std::endl;
579 
580  std::cout << "================================================" << std::endl;
581 }
582 
583 void EWSM_Output::outputEachDeltaKappaZ(const double f_AlphaToGF,
584  const double cW2overSW2,
585  const double DeltaRho[StandardModel::orders_EW_size],
586  const double deltaKappa_rem[StandardModel::orders_EW_size],
587  const double DeltaRbar_rem,
588  const bool bool_Zbb,
589  const double taub[StandardModel::orders_EW_size],
590  const double ZbbSubtract,
591  const double Zgamma_EW2) const
592 {
593  /* rescale */
594  double DeltaRho_new[StandardModel::orders_EW_size];
595  for (int j = 0; j < StandardModel::orders_EW_size; ++j)
596  DeltaRho_new[j] = cW2overSW2 * DeltaRho[j];
597 
598  if (getFlagKappaZ().compare("APPROXIMATEFORMULA") == 0) {
599  std::cout << "Delta kappaZb (from the approximate formula of sin2thb) = "
600  << kappaZ_f(getQuarks(QCD::BOTTOM)) - 1.0 << std::endl;
601  } else if (getFlagKappaZ().compare("NORESUM") == 0) {
602  std::cout << "Leading contributions: alpha or Gmu" << std::endl;
603  std::cout << " DeltaRho[EW1]=" << DeltaRho_new[StandardModel::EW1] << " "
604  << f_AlphaToGF * DeltaRho_new[StandardModel::EW1] << std::endl;
605  std::cout << " DeltaRho[EW1QCD1]=" << DeltaRho_new[StandardModel::EW1QCD1] << " "
606  << f_AlphaToGF * DeltaRho_new[StandardModel::EW1QCD1] << std::endl;
607  std::cout << " DeltaRho[EW1QCD2]=" << DeltaRho_new[StandardModel::EW1QCD2] << " "
608  << f_AlphaToGF * DeltaRho_new[StandardModel::EW1QCD2] << std::endl;
609  std::cout << " DeltaRho[EW2]=" << DeltaRho_new[StandardModel::EW2] << " "
610  << pow(f_AlphaToGF, 2.0) * DeltaRho_new[StandardModel::EW2] << std::endl;
611  std::cout << " DeltaRho[EW2QCD1]=" << DeltaRho_new[StandardModel::EW2QCD1] << " "
612  << pow(f_AlphaToGF, 2.0) * DeltaRho_new[StandardModel::EW2QCD1] << std::endl;
613  std::cout << " DeltaRho[EW3]=" << DeltaRho_new[StandardModel::EW3] << " "
614  << pow(f_AlphaToGF, 3.0) * DeltaRho_new[StandardModel::EW3] << std::endl;
615  std::cout << "EW2 from Z-gamma = " << Zgamma_EW2 << std::endl;
616  std::cout << "Remainder contributions: alpha or Gmu" << std::endl;
617  std::cout << " DeltaRbar_rem[EW1]=" << DeltaRbar_rem << " "
618  << f_AlphaToGF * DeltaRbar_rem << std::endl;
619  std::cout << " deltaKappa_rem[EW1]=" << deltaKappa_rem[StandardModel::EW1] - ZbbSubtract << " "
620  << f_AlphaToGF * (deltaKappa_rem[StandardModel::EW1] - ZbbSubtract) << std::endl;
621  std::cout << " deltaKappa_rem[EW1QCD1]=" << deltaKappa_rem[StandardModel::EW1QCD1] << " "
622  << f_AlphaToGF * deltaKappa_rem[StandardModel::EW1QCD1] << std::endl;
623  std::cout << " deltaKappa_rem[EW1QCD2]=" << deltaKappa_rem[StandardModel::EW1QCD2] << " "
624  << f_AlphaToGF * deltaKappa_rem[StandardModel::EW1QCD2] << std::endl;
625  std::cout << " deltaKappa_rem[EW2]=" << deltaKappa_rem[StandardModel::EW2] << " "
626  << pow(f_AlphaToGF, 2.0) * deltaKappa_rem[StandardModel::EW2] << std::endl;
627  if (bool_Zbb) {
628  std::cout << "Taub: alpha or Gmu" << std::endl;
629  std::cout << " taub[EW1]=" << taub[StandardModel::EW1] << " "
630  << f_AlphaToGF * taub[StandardModel::EW1] << std::endl;
631  std::cout << " taub[EW1QCD1]=" << taub[StandardModel::EW1QCD1] << " "
632  << f_AlphaToGF * taub[StandardModel::EW1QCD1] << std::endl;
633  std::cout << " taub[EW2]=" << taub[StandardModel::EW2] << " "
634  << pow(f_AlphaToGF, 2.0) * taub[StandardModel::EW2] << std::endl;
635  }
636  std::cout << "Each order: alpha or Gmu" << std::endl;
637  double dKappa_EW1 = DeltaRho_new[StandardModel::EW1] + deltaKappa_rem[StandardModel::EW1] - ZbbSubtract;
638  double dKappa_EW1QCD1 = DeltaRho_new[StandardModel::EW1QCD1] + deltaKappa_rem[StandardModel::EW1QCD1];
639  double dKappa_EW2 = (deltaKappa_rem[StandardModel::EW1] - ZbbSubtract) * DeltaRho_new[StandardModel::EW1]
640  - DeltaRho_new[StandardModel::EW1] * DeltaRbar_rem
641  + DeltaRho_new[StandardModel::EW2]
642  + deltaKappa_rem[StandardModel::EW2];
643  double dKappa_EW1QCD2 = DeltaRho_new[StandardModel::EW1QCD2] + deltaKappa_rem[StandardModel::EW1QCD2];
644  double dKappa_EW2QCD1 = DeltaRho_new[StandardModel::EW2QCD1] + deltaKappa_rem[StandardModel::EW1QCD1] * DeltaRho_new[StandardModel::EW1];
645  double dKappa_EW3 = DeltaRho_new[StandardModel::EW3];
646  double dKappa_EW2QCD2 = deltaKappa_rem[StandardModel::EW1QCD2] * DeltaRho_new[StandardModel::EW1];
647  //
648  double dKappa_EW1_TMP = dKappa_EW1;
649  double dKappa_EW1QCD1_TMP = dKappa_EW1QCD1;
650  double dKappa_EW2_TMP = dKappa_EW2;
651  double dKappa_EW1QCD2_TMP = dKappa_EW1QCD2;
652  double dKappa_EW2QCD1_TMP = dKappa_EW2QCD1;
653  double dKappa_EW3_TMP = dKappa_EW3;
654  double dKappa_EW2QCD2_TMP = dKappa_EW2QCD2;
655  //
656  if (bool_Zbb) {
657  dKappa_EW1 = dKappa_EW1_TMP - taub[StandardModel::EW1];
658  dKappa_EW1QCD1 = dKappa_EW1QCD1_TMP - taub[StandardModel::EW1QCD1];
659  dKappa_EW2 = dKappa_EW2_TMP - taub[StandardModel::EW2] + taub[StandardModel::EW1] * taub[StandardModel::EW1]
660  - dKappa_EW1_TMP * taub[StandardModel::EW1];
661  dKappa_EW1QCD2 = dKappa_EW1QCD2_TMP;
662  dKappa_EW2QCD1 = dKappa_EW2QCD1_TMP - dKappa_EW1_TMP * taub[StandardModel::EW1QCD1]
663  - dKappa_EW1QCD1_TMP * taub[StandardModel::EW1] + 2.0 * taub[StandardModel::EW1] * taub[StandardModel::EW1QCD1];
664  dKappa_EW3 = dKappa_EW3_TMP - dKappa_EW1_TMP * taub[StandardModel::EW2]
665  - dKappa_EW2_TMP * taub[StandardModel::EW1] + 2.0 * taub[StandardModel::EW2] * taub[StandardModel::EW1];
666  dKappa_EW2QCD2 = dKappa_EW2QCD2_TMP - dKappa_EW1QCD2_TMP * taub[StandardModel::EW1]
667  - dKappa_EW1QCD1_TMP * taub[StandardModel::EW1QCD1]
668  + taub[StandardModel::EW1QCD1] * taub[StandardModel::EW1QCD1];
669  }
670  std::cout << " EW1: " << dKappa_EW1 << " " << f_AlphaToGF * dKappa_EW1 << std::endl;
671  std::cout << " EW1QCD1: " << dKappa_EW1QCD1 << " " << f_AlphaToGF * dKappa_EW1QCD1 << std::endl;
672  std::cout << " EW2: " << dKappa_EW2 + Zgamma_EW2
673  << " " << pow(f_AlphaToGF, 2.0) * dKappa_EW2 + Zgamma_EW2 << std::endl;
674  std::cout << " EW1QCD2: " << dKappa_EW1QCD2 << " " << f_AlphaToGF * dKappa_EW1QCD2 << std::endl;
675  std::cout << " EW2QCD1: " << dKappa_EW2QCD1 << " " << pow(f_AlphaToGF, 2.0) * dKappa_EW2QCD1 << std::endl;
676  std::cout << " EW3: " << dKappa_EW3 << " " << pow(f_AlphaToGF, 3.0) * dKappa_EW3 << std::endl;
677  std::cout << " EW2QCD2: " << dKappa_EW2QCD2 << " " << pow(f_AlphaToGF, 2.0) * dKappa_EW2QCD2 << std::endl;
678  std::cout << "Total contribution: alpha or Gmu" << std::endl;
679  std::cout << " kappaZ="
680  << 1.0 + dKappa_EW1 + dKappa_EW1QCD1 + dKappa_EW2
681  + dKappa_EW1QCD2 + dKappa_EW2QCD1 + dKappa_EW3
682  + dKappa_EW2QCD2 + Zgamma_EW2
683  << " "
684  << 1.0 + f_AlphaToGF * dKappa_EW1 + f_AlphaToGF * dKappa_EW1QCD1
685  + pow(f_AlphaToGF, 2.0) * dKappa_EW2 + f_AlphaToGF * dKappa_EW1QCD2
686  + pow(f_AlphaToGF, 2.0) * dKappa_EW2QCD1
687  + pow(f_AlphaToGF, 3.0) * dKappa_EW3
688  + pow(f_AlphaToGF, 2.0) * dKappa_EW2QCD2
689  + Zgamma_EW2 << std::endl;
690  if (bool_Zbb) {
691  std::cout << " kappaZ(taub resummed)="
692  << (1.0 + dKappa_EW1_TMP + dKappa_EW1QCD1_TMP
693  + dKappa_EW2_TMP + dKappa_EW1QCD2_TMP
694  + dKappa_EW2QCD1_TMP + dKappa_EW3_TMP
695  + dKappa_EW2QCD2_TMP)
696  / (1.0 + taub[StandardModel::EW1] + taub[StandardModel::EW1QCD1] + taub[StandardModel::EW2])
697  + Zgamma_EW2
698  << " "
699  << (1.0 + f_AlphaToGF * dKappa_EW1_TMP
700  + f_AlphaToGF * dKappa_EW1QCD1_TMP
701  + pow(f_AlphaToGF, 2.0) * dKappa_EW2_TMP
702  + f_AlphaToGF * dKappa_EW1QCD2_TMP
703  + pow(f_AlphaToGF, 2.0) * dKappa_EW2QCD1_TMP
704  + pow(f_AlphaToGF, 3.0) * dKappa_EW3_TMP
705  + pow(f_AlphaToGF, 2.0) * dKappa_EW2QCD2_TMP)
706  / (1.0 + f_AlphaToGF * taub[StandardModel::EW1]
707  + f_AlphaToGF * taub[StandardModel::EW1QCD1]
708  + pow(f_AlphaToGF, 2.0) * taub[StandardModel::EW2])
709  + Zgamma_EW2
710  << std::endl;
711  }
712  } else
713  std::cout << "EWSM_Output::outputEachDeltaKappaZ(): Not implemented for schemeKappaZ="
714  << getFlagKappaZ() << std::endl;
715 }
716 
717 
718 
EWSMThreeLoopEW::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMThreeLoopEW.cpp:109
StandardModel::getMyOneLoopEW
EWSMOneLoopEW * getMyOneLoopEW() const
A get method to retrieve the member pointer of type EWSMOneLoopEW,.
Definition: StandardModel.h:970
StandardModel::EW2
Two-loop of .
Definition: StandardModel.h:499
StandardModel::taub
double taub() const
Top-mass corrections to the vertex, denoted by .
Definition: StandardModel.cpp:2104
EWSMThreeLoopQCD.h
QCD::BOTTOM
Definition: QCD.h:329
StandardModel::DeltaRbar
virtual double DeltaRbar() const
The SM prediction for derived from that for the -boson mass.
Definition: StandardModel.cpp:1146
EWSMThreeLoopEW::DeltaRho
double DeltaRho(const double Mw_i) const
Leading three-loop contribution of to , denoted as .
Definition: EWSMThreeLoopEW.cpp:70
EWSMOneLoopEW::DeltaRbar_rem
double DeltaRbar_rem(const double Mw_i) const
.
Definition: EWSMOneLoopEW.cpp:73
EWSM_Output::outputEachDeltaRhoZ
void outputEachDeltaRhoZ(const double f_AlphaToGF, const double DeltaRho[StandardModel::orders_EW_size], const double deltaRho_rem[StandardModel::orders_EW_size], const double DeltaRbar_rem, const bool bool_Zbb, const double taub[StandardModel::orders_EW_size], const double ZbbSubtract) const
Definition: EWSM_Output.cpp:336
EWSMOneLoopEW::DeltaRho
double DeltaRho(const double Mw_i) const
Leading one-loop contribution of to , denoted as .
Definition: EWSMOneLoopEW.cpp:43
StandardModel::getFlagRhoZ
std::string getFlagRhoZ() const
A method to retrieve the model flag RhoZ.
Definition: StandardModel.h:672
StandardModel::EW1
One-loop of .
Definition: StandardModel.h:496
StandardModel::getMyEWSMcache
EWSMcache * getMyEWSMcache() const
A get method to retrieve the member pointer of type EWSMcache.
Definition: StandardModel.h:961
StandardModel::getMyTwoLoopQCD
EWSMTwoLoopQCD * getMyTwoLoopQCD() const
Definition: StandardModel.h:1015
StandardModel::getMyThreeLoopQCD
EWSMThreeLoopQCD * getMyThreeLoopQCD() const
Definition: StandardModel.h:1005
StandardModel::alphaMz
double alphaMz() const
The electromagnetic coupling at the -mass scale, .
Definition: StandardModel.cpp:893
StandardModel::DeltaAlphaL5q
double DeltaAlphaL5q() const
The sum of the leptonic and the five-flavour hadronic corrections to the electromagnetic coupling at...
Definition: StandardModel.cpp:856
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
EWSM_Output.h
EWSMOneLoopEW::deltaKappa_rem_f
gslpp::complex deltaKappa_rem_f(const Particle f, const double Mw_i) const
Remainder contribution of to the effective couplings , denoted as .
Definition: EWSMOneLoopEW.cpp:143
StandardModel::getMyThreeLoopEW
EWSMThreeLoopEW * getMyThreeLoopEW() const
Definition: StandardModel.h:995
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:477
EWSMOneLoopEW::deltaRho_rem_f
gslpp::complex deltaRho_rem_f(const Particle f, const double Mw_i) const
Remainder contribution of to the effective couplings , denoted as .
Definition: EWSMOneLoopEW.cpp:113
EWSMThreeLoopEW2QCD::DeltaRho
double DeltaRho(const double Mw_i) const
Leading three-loop contribution of to , denoted as .
Definition: EWSMThreeLoopEW2QCD.cpp:29
StandardModel::DeltaR
virtual double DeltaR() const
The SM prediction for derived from that for the boson mass.
Definition: StandardModel.cpp:1036
EWSMTwoLoopQCD::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:44
EWSMThreeLoopEW.h
EWSM_Output::outputEachDeltaR
void outputEachDeltaR(const double Mw_i) const
Definition: EWSM_Output.cpp:25
StandardModel::EW1QCD1
Two-loop of .
Definition: StandardModel.h:497
EWSMOneLoopEW.h
EWSM_Output::outputEachDeltaKappaZ_q
void outputEachDeltaKappaZ_q(const QCD::quark q, const double Mw_i) const
Definition: EWSM_Output.cpp:506
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
StandardModel::getGF
double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:739
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
EWSMThreeLoopQCD::DeltaRho
double DeltaRho(const double Mw_i) const
Leading three-loop QCD contribution of to , denoted as .
Definition: EWSMThreeLoopQCD.cpp:42
EWSM_Output::outputEachDeltaKappaZ
void outputEachDeltaKappaZ(const double f_AlphaToGF, const double cW2overSW2, const double DeltaRho[StandardModel::orders_EW_size], const double deltaKappa_rem[StandardModel::orders_EW_size], const double DeltaRbar_rem, const bool bool_Zbb, const double taub[StandardModel::orders_EW_size], const double ZbbSubtract, const double Zgamma_EW2) const
Definition: EWSM_Output.cpp:583
EWSMApproximateFormulae.h
StandardModel::EW3
Three-loop of .
Definition: StandardModel.h:501
EWSM_Output::outputEachDeltaRhoZ_q
void outputEachDeltaRhoZ_q(const QCD::quark q, const double Mw_i) const
Definition: EWSM_Output.cpp:266
QCD::getMtpole
double getMtpole() const
A get method to access the pole mass of the top quark.
Definition: QCD.h:588
EWSMTwoLoopQCD::DeltaRho
double DeltaRho(const double Mw_i) const
Leading two-loop QCD contribution of to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:38
EWSMcache.h
EWSMOneLoopEW_HV.h
EWSMThreeLoopEW2QCD::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMThreeLoopEW2QCD.cpp:63
StandardModel::getMyTwoLoopEW
EWSMTwoLoopEW * getMyTwoLoopEW() const
Definition: StandardModel.h:1010
EWSMTwoLoopEW::tau_2
double tau_2() const
The function .
Definition: EWSMTwoLoopEW.cpp:151
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
StandardModel::getMyApproximateFormulae
EWSMApproximateFormulae * getMyApproximateFormulae() const
A get method to retrieve the member pointer of type EWSMApproximateFormulae.
Definition: StandardModel.h:979
EWSMApproximateFormulae::DeltaR_TwoLoopEW_rem
double DeltaR_TwoLoopEW_rem(const double Mw_i) const
.
Definition: EWSMApproximateFormulae.cpp:252
QCD::quark
quark
An enum type for quarks.
Definition: QCD.h:323
EWSM_Output::outputEachDeltaRhoZ_l
void outputEachDeltaRhoZ_l(const QCD::lepton l, const double Mw_i) const
Definition: EWSM_Output.cpp:211
EWSMThreeLoopEW2QCD.h
EWSMTwoLoopEW::DeltaRho
double DeltaRho(const double Mw_i) const
Leading two-loop contribution of to , denoted as .
Definition: EWSMTwoLoopEW.cpp:54
StandardModel::resumKappaZ
double resumKappaZ(const double DeltaRho[orders_EW_size], const double deltaKappa_rem[orders_EW_size], const double DeltaRbar_rem, const bool bool_Zbb) const
A method to compute the real part of the effetvive coupling from , and .
Definition: StandardModel.cpp:2027
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:721
EWSMOneLoopEW::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMOneLoopEW.cpp:49
StandardModel::getFlagKappaZ
std::string getFlagKappaZ() const
A method to retrieve the model flag KappaZ.
Definition: StandardModel.h:682
EWSMThreeLoopQCD::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMThreeLoopQCD.cpp:48
EWSMTwoFermionsLEP2.h
StandardModel::getFlagMw
std::string getFlagMw() const
A method to retrieve the model flag Mw.
Definition: StandardModel.h:662
EWSM_Output::outputEachDeltaKappaZ_l
void outputEachDeltaKappaZ_l(const QCD::lepton l, const double Mw_i) const
Definition: EWSM_Output.cpp:443
EWSMTwoLoopEW.h
StandardModel::EW2QCD1
Three-loop of .
Definition: StandardModel.h:500
StandardModel::Mw
virtual double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel.cpp:970
EWSMcache::alsMt
double alsMt() const
The strong coupling at NNLO.
Definition: EWSMcache.h:378
StandardModel::kappaZ_f
virtual gslpp::complex kappaZ_f(const Particle f) const
The effective leptonic neutral-current coupling in the SM.
Definition: StandardModel.cpp:1644
StandardModel::EW1QCD2
Three-loop of .
Definition: StandardModel.h:498
EWSMcache::Xt_alpha
double Xt_alpha(const double Mw_i) const
The quantity with the coupling .
Definition: EWSMcache.h:355
EWSMTwoLoopQCD.h
StandardModel::getAle
double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:748
EWSM_Output::EWSM_Output
EWSM_Output(const StandardModel &SM_in)
Constructor.
Definition: EWSM_Output.cpp:20
StandardModel::orders_EW_size
The size of this enum.
Definition: StandardModel.h:502
EWSMTwoLoopEW::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMTwoLoopEW.cpp:81
StandardModel::getMyThreeLoopEW2QCD
EWSMThreeLoopEW2QCD * getMyThreeLoopEW2QCD() const
Definition: StandardModel.h:1000
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
StandardModel::getLeptons
Particle getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:712
StandardModel::orders_EW
orders_EW
An enumerated type representing perturbative orders of radiative corrections to EW precision observab...
Definition: StandardModel.h:495