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