a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EWSUSY.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 <iostream>
9#include <iomanip>
10#include <stdexcept>
11#include <cmath>
12#include "EWSUSY.h"
13#include "SUSY.h"
14#include <EWSMcache.h>
15#include "EWSMOneLoopEW.h"
16#include "EWSMTwoLoopQCD.h"
17#include "EWSMThreeLoopQCD.h"
18#include "EWSMTwoLoopEW.h"
19#include "EWSMThreeLoopEW2QCD.h"
20#include "EWSMThreeLoopEW.h"
21#include "EWSMApproximateFormulae.h"
22#include "EWSMOneLoopEW_HV.h"
23/* BEGIN: REMOVE FROM THE PACKAGE */
24#include "FeynHiggsWrapper.h"
25#include "EWSMTwoFermionsLEP2.h"
26/* END: REMOVE FROM THE PACKAGE */
27
28
29const double EWSUSY::Mw_unphysical = 2.0;
30const double EWSUSY::RenormalizationScaleFactor = 1.0;
31//const double EWSUSY::RenormalizationScaleFactor = 2.0; // for debug
32
33EWSUSY::EWSUSY(const SUSY& SUSY_in)
34: PV(true), mySUSY(SUSY_in),
35 Yu(3,3,0.0), Yd(3,3,0.0), Yl(3,3,0.0),
36 Au(3,3,0.0), Ad(3,3,0.0), Al(3,3,0.0),
37 Zm(2,2,0.), Zp(2,2,0.), ZN(4,4,0.),
38 ZU(6,6,0.), ZD(6,6,0.), ZL(6,6,0.), Zne(6,6,0.),
39 ZR(2,2,0.), ZH(2,2,0.)
40{
41}
42
43void EWSUSY::SetRosiekParameters()
44{
45 Yu = mySUSY.getYu();
46 Yd = - mySUSY.getYd();
47 Yl = - mySUSY.getYe();
48
49 Au = - mySUSY.getTUhat().transpose();
50 Ad = mySUSY.getTDhat().transpose();
51 Al = mySUSY.getTEhat().transpose();
52
53 Zm = mySUSY.getU().hconjugate();
54 Zp = mySUSY.getV().hconjugate();
55 ZN = mySUSY.getN().hconjugate();
56
57 ZU = mySUSY.getRu().hconjugate();
58 ZD = mySUSY.getRd().transpose();
59 ZL = mySUSY.getRl().transpose();
60 Zne = mySUSY.getRn().hconjugate();
61
62 double sinAlpha = mySUSY.getSaeff().real(); /* Correct? */
63 double cosAlpha = sqrt(1.0 - sinAlpha*sinAlpha); /* -Pi/2 < alpha < 0 */
64 ZR.assign(0,0, cosAlpha);
65 ZR.assign(0,1, - sinAlpha);
66 ZR.assign(1,0, sinAlpha);
67 ZR.assign(1,1, cosAlpha);
68
69 ZH.assign(0,0, mySUSY.getSinb());
70 ZH.assign(0,1, - mySUSY.getCosb());
71 ZH.assign(1,0, mySUSY.getCosb());
72 ZH.assign(1,1, mySUSY.getSinb());
73
74 /* particle masses */
75 for (int I=0; I<3; ++I) {
76 /* up-type quarks */
77 m_u[I] = mySUSY.getMyEWSMcache()->mf(mySUSY.getQuarks((QCD::quark)(2*I)), mySUSY.getMz(), FULLNNLO);
78 /* down-type quarks */
79 m_d[I] = mySUSY.getMyEWSMcache()->mf(mySUSY.getQuarks((QCD::quark)(2*I + 1)), mySUSY.getMz(), FULLNNLO);
80 /* charged leptons */
81 m_l[I] = mySUSY.getLeptons((QCD::lepton)(2*I + 1)).getMass();
82 }
83 /* H^0_i = (H^0, h^0, A^0, G^0) */
84 mH02[0] = mySUSY.getMHh()*mySUSY.getMHh();
85 mH02[1] = mySUSY.getMHl()*mySUSY.getMHl();
86 mH02[2] = mySUSY.getMHa()*mySUSY.getMHa();
87 mH02[3] = mySUSY.getMz()*mySUSY.getMz(); /* mass squared of the neutral Goldstone boson */
88 for (int k=0; k<6; ++k) {
89 Msu2[k] = mySUSY.getMsu2()(k);
90 Msd2[k] = mySUSY.getMsd2()(k);
91 Mse2[k] = mySUSY.getMse2()(k);
92 }
93 for (int k=0; k<3; ++k)
94 Msn2[k] = mySUSY.getMsn2()(k);
95 for (int i=0; i<2; ++i)
96 mC[i] = mySUSY.getMch()(i);
97 for (int j=0; j<4; ++j)
98 mN[j] = mySUSY.getMneu()(j);
99}
100
101gslpp::complex EWSUSY::FA(const double mu, const double p2,
102 const double mi, const double mj,
103 const gslpp::complex cV_aij, const gslpp::complex cV_bji,
104 const gslpp::complex cA_aij, const gslpp::complex cA_bji) const
105{
106 double mu2 = mu*mu, mi2 = mi*mi, mj2 = mj*mj;
107
108 /* PV functions */
109 double A0i = PV.A0(mu2, mi2);
110 double A0j = PV.A0(mu2, mj2);
111 gslpp::complex B0 = PV.B0(mu2, p2, mi2, mj2);
112 gslpp::complex B00 = PV.B00(mu2, p2, mi2, mj2);
113
114 return ( -2.0*(cV_aij*cV_bji + cA_aij*cA_bji)
115 *(4.0*B00 + A0i + A0j + (p2 - mi*mi - mj*mj)*B0)
116 -4.0*(cV_aij*cV_bji - cA_aij*cA_bji)*mi*mj*B0 );
117}
118
119gslpp::complex EWSUSY::dFA(const double mu, const double p2,
120 const double mi, const double mj,
121 const gslpp::complex cV_aij, const gslpp::complex cV_bji,
122 const gslpp::complex cA_aij, const gslpp::complex cA_bji) const
123{
124 double mu2 = mu*mu, mi2 = mi*mi, mj2 = mj*mj;
125
126 /* PV functions */
127 gslpp::complex B0 = PV.B0(mu2, p2, mi2, mj2);
128 gslpp::complex B0p = PV.B0p(mu2, p2, mi2, mj2);
129 gslpp::complex B00p = PV.B00p(mu2, p2, mi2, mj2);
130
131 if (mi == mj && cA_aij == 0.0 && cA_bji == 0.0)
132 return ( -2.0*cV_aij*cV_bji*(4.0*B00p + p2*B0p + B0) );
133 else
134 return ( -2.0*(cV_aij*cV_bji + cA_aij*cA_bji)
135 *(4.0*B00p + (p2 - mi*mi - mj*mj)*B0p + B0)
136 -4.0*(cV_aij*cV_bji - cA_aij*cA_bji)*mi*mj*B0p );
137}
138
139gslpp::complex EWSUSY::PiT_Z(const double mu, const double p2, const double Mw_i) const
140{
141 double mu2 = mu*mu;
142 double e2 = 4.0*M_PI*mySUSY.getAle();
143 double e = sqrt(e2);
144 double Mz = mySUSY.getMz();
145 double Nc = mySUSY.getNc();
146
147 /* variables depending on Mw_i */
148 double Mw2 = Mw_i*Mw_i;
149 double mHp2[2] = {mySUSY.getMHp()*mySUSY.getMHp(), Mw2}; /* H^+_i = (H^+, G^+) */
150 double cW = Mw_i/Mz;
151 double cW2 = cW*cW;
152 double sW2 = 1.0 - cW2;
153 double sW = sqrt(sW2);
154 double g2sq = e2/sW2; /* g2 squared */
155 double e_4sc = e/4.0/sW/cW;
156 double e_2sc = 2.0*e_4sc;
157
158 gslpp::complex PiT_f = gslpp::complex(0.0, 0.0, false);
159 gslpp::complex PiT_sf = gslpp::complex(0.0, 0.0, false);
160 gslpp::complex PiT_ch = gslpp::complex(0.0, 0.0, false);
161 gslpp::complex PiT_WZH = gslpp::complex(0.0, 0.0, false);
162 double a0;
163 gslpp::complex b0, b00;
164 gslpp::complex cV_Zij, cV_Zji, cA_Zij, cA_Zji;
165 gslpp::matrix<double> Id6 = gslpp::matrix<double>::Id(6);
166 gslpp::matrix<double> Id2 = gslpp::matrix<double>::Id(2);
167
168 /* neutrino loops */
169 b0 = PV.B0(mu2, p2, 0.0, 0.0);
170 b00 = PV.B00(mu2, p2, 0.0, 0.0);
171 PiT_f += - 3.0/4.0*g2sq/cW2*(4.0*b00 + p2*b0);
172
173 /* other SM fermion loops */
174 gslpp::complex cV_Zee = - e_4sc*(1.0 - 4.0*sW2);
175 gslpp::complex cA_Zee = - e_4sc;
176 gslpp::complex cV_Zdd = - e_4sc*(1.0 - 4.0/3.0*sW2);
177 gslpp::complex cA_Zdd = - e_4sc;
178 gslpp::complex cV_Zuu = e_4sc*(1.0 - 8.0/3.0*sW2);
179 gslpp::complex cA_Zuu = e_4sc;
180 for (int I=0; I<3; ++I) {
181 /* charged leptons */
182 PiT_f += FA(mu, p2, m_l[I], m_l[I], cV_Zee, cV_Zee, cA_Zee, cA_Zee);
183
184 /* down-type quarks */
185 PiT_f += Nc*FA(mu, p2, m_d[I], m_d[I], cV_Zdd, cV_Zdd, cA_Zdd, cA_Zdd);
186
187 /* up-type quarks */
188 PiT_f += Nc*FA(mu, p2, m_u[I], m_u[I], cV_Zuu, cV_Zuu, cA_Zuu, cA_Zuu);
189 }
190
191 /* sneutrino loops */
192 gslpp::complex VZsnsn_II = e_2sc;
193 gslpp::complex VZZsnsn_II = e2/2.0/sW2/cW2;
194 for (int I=0; I<3; ++I) { /* I=0-2 for left-handed sneutrinos */
195 b00 = PV.B00(mu2, p2, Msn2[I], Msn2[I]);
196 PiT_sf += 4.0*VZsnsn_II.abs2()*b00;
197 a0 = PV.A0(mu2, Msn2[I]);
198 PiT_sf += VZZsnsn_II*a0;
199 }
200
201 /* charged-slepton loops */
202 gslpp::complex VZLL_mn, VZZLL_nn;
203 for (int n=0; n<6; ++n) {
204 for (int m=0; m<6; ++m) {
205 VZLL_mn = gslpp::complex(0.0, 0.0, false);
206 for (int I=0; I<3; ++I) /* sum over left-handed sleptons */
207 VZLL_mn += - e_2sc*ZL(I,n)*ZL(I,m).conjugate();
208 VZLL_mn += - e_2sc*(- 2.0*sW2*Id6(m,n));
209 b00 = PV.B00(mu2, p2, Mse2[m], Mse2[n]);
210 PiT_sf += 4.0*VZLL_mn.abs2()*b00;
211 }
212 VZZLL_nn = gslpp::complex(0.0, 0.0, false);
213 VZZLL_nn += 2.0*e2/cW2*sW2;
214 for (int I=0; I<3; ++I) /* sum over left-handed sleptons */
215 VZZLL_nn += 2.0*e2/cW2*(1.0 - 4.0*sW2)/4.0/sW2*ZL(I,n)*ZL(I,n).conjugate();
216 a0 = PV.A0(mu2, Mse2[n]);
217 PiT_sf += VZZLL_nn*a0;
218 }
219
220 /* down-type squark loops */
221 gslpp::complex VZDD_mn, VZZDD_nn;
222 for (int n=0; n<6; ++n) {
223 for (int m=0; m<6; ++m) {
224 VZDD_mn = gslpp::complex(0.0, 0.0, false);
225 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
226 VZDD_mn += - e_2sc*ZD(I,n)*ZD(I,m).conjugate();
227 VZDD_mn += - e_2sc*(- 2.0/3.0*sW2*Id6(m,n));
228 b00 = PV.B00(mu2, p2, Msd2[m], Msd2[n]);
229 PiT_sf += 4.0*Nc*VZDD_mn.abs2()*b00;
230 }
231 VZZDD_nn = gslpp::complex(0.0, 0.0, false);
232 VZZDD_nn += 2.0*e2/3.0/cW2*sW2/3.0;
233 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
234 VZZDD_nn += 2.0*e2/3.0/cW2*(3.0 - 4.0*sW2)/4.0/sW2*ZD(I,n)*ZD(I,n).conjugate();
235 a0 = PV.A0(mu2, Msd2[n]);
236 PiT_sf += Nc*VZZDD_nn*a0;
237 }
238
239 /* up-type squark loops */
240 gslpp::complex VZUU_mn, VZZUU_nn;
241 for (int n=0; n<6; ++n) {
242 for (int m=0; m<6; ++m) {
243 VZUU_mn = gslpp::complex(0.0, 0.0, false);
244 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
245 VZUU_mn += e_2sc*ZU(I,m).conjugate()*ZU(I,n);
246 VZUU_mn += e_2sc*(- 4.0/3.0*sW2*Id6(m,n));
247 b00 = PV.B00(mu2, p2, Msu2[m], Msu2[n]);
248 PiT_sf += 4.0*Nc*VZUU_mn.abs2()*b00;
249 }
250 VZZUU_nn = gslpp::complex(0.0, 0.0, false);
251 VZZUU_nn += 2.0*e2/3.0/cW2*4.0*sW2/3.0;
252 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
253 VZZUU_nn += 2.0*e2/3.0/cW2*(3.0 - 8.0*sW2)/4.0/sW2*ZU(I,n).conjugate()*ZU(I,n);
254 a0 = PV.A0(mu2, Msu2[n]);
255 PiT_sf += Nc*VZZUU_nn*a0;
256 }
257
258 /* chargino loops */
259 for (int i=0; i<2; ++i)
260 for (int j=0; j<2; ++j) {
261 cV_Zij = e_4sc*( Zp(0,j).conjugate()*Zp(0,i)
262 + Zm(0,j)*Zm(0,i).conjugate()
263 + 2.0*(cW2 - sW2)*Id2(j,i) );
264 cV_Zji = cV_Zij.conjugate();
265 cA_Zij = e_4sc*( Zp(0,j).conjugate()*Zp(0,i)
266 - Zm(0,j)*Zm(0,i).conjugate() );
267 cA_Zji = cA_Zij.conjugate();
268 PiT_ch += FA(mu, p2, mC[i], mC[j], cV_Zij, cV_Zji, cA_Zij, cA_Zji);
269 }
270
271 /* neutralino loops */
272 for (int i=0; i<4; ++i)
273 for (int j=0; j<4; ++j) {
274 cV_Zij = - e_4sc*( ZN(3,j).conjugate()*ZN(3,i)
275 - ZN(2,j).conjugate()*ZN(2,i)
276 - ZN(3,j)*ZN(3,i).conjugate()
277 + ZN(2,j)*ZN(2,i).conjugate() );
278 cV_Zji = cV_Zij.conjugate();
279 cA_Zij = - e_4sc*( ZN(3,j).conjugate()*ZN(3,i)
280 - ZN(2,j).conjugate()*ZN(2,i)
281 + ZN(3,j)*ZN(3,i).conjugate()
282 - ZN(2,j)*ZN(2,i).conjugate() );
283 cA_Zji = cA_Zij.conjugate();
284 PiT_ch += 0.5*FA(mu, p2, mN[i], mN[j], cV_Zij, cV_Zji, cA_Zij, cA_Zji);
285 }
286
287 /* charged-Higgs loops */
288 double cot_2thW = (cW2 - sW2)/(2.0*sW*cW);
289 for (int i=0; i<2; ++i) {
290 b00 = PV.B00(mu2, p2, mHp2[i], mHp2[i]);
291 a0 = PV.A0(mu2, mHp2[i]);
292 PiT_WZH += 2.0*e2*cot_2thW*cot_2thW*(2.0*b00 + a0);
293 }
294
295 /* neutral-Higgs loops */
296 double AM_ij;
297 for (int i=0; i<2; ++i)
298 for (int j=0; j<2; ++j) {
299 AM_ij = ZR(0,i)*ZH(0,j) - ZR(1,i)*ZH(1,j);
300 b00 = PV.B00(mu2, p2, mH02[i], mH02[j+2]);
301 PiT_WZH += g2sq/cW2*AM_ij*AM_ij*b00;
302 }
303 for (int j=0; j<4; ++j) {
304 a0 = PV.A0(mu2, mH02[j]);
305 PiT_WZH += g2sq/4.0/cW2*a0;
306 }
307
308 /* W-boson - charged-Goldstone-boson loop*/
309 b0 = PV.B0(mu2, p2, Mw2, Mw2);
310 PiT_WZH += - 2.0*g2sq*sW2*sW2*Mz*Mz*b0;
311
312 /* Z-boson - Higgs loops */
313 double CR_i;
314 for (int i=0; i<2; ++i) {
315 CR_i = mySUSY.v1()*ZR(0,i) + mySUSY.v2()*ZR(1,i);
316 b0 = PV.B0(mu2, p2, Mz*Mz, mH02[i]);
317 /* Mw^2/v^2 is substituted for g2^2/4 compared to the expression in the
318 * paper, in order to ensure the cancellation of the UV divergences in
319 * the case where Mw is not the tree-level value. */
320 PiT_WZH += - g2sq*Mw2/mySUSY.v()/mySUSY.v()/cW2/cW2*CR_i*CR_i*b0;
321 }
322
323 /* W-boson loops */
324 a0 = PV.A0(mu2, Mw2);
325 b0 = PV.B0(mu2, p2, Mw2, Mw2);
326 b00 = PV.B00(mu2, p2, Mw2, Mw2);
327 /* typo in the paper: a0^2 --> a0 in the first term */
328 PiT_WZH += 2.0*g2sq*cW2*(2.0*a0 + (2.0*p2 + Mw2)*b0 + 4.0*b00);
329
330 /* Sum of all contributions */
331 gslpp::complex PiT = PiT_f + PiT_sf + PiT_ch + PiT_WZH;
332
333 return ( PiT/16.0/M_PI/M_PI );
334}
335
336
337gslpp::complex EWSUSY::PiT_W(const double mu, const double p2, const double Mw_i) const
338{
339 double mu2 = mu*mu;
340 double e2 = 4.0*M_PI*mySUSY.getAle();
341 double e = sqrt(e2);
342 double Mz = mySUSY.getMz();
343 double Nc = mySUSY.getNc();
344
345 /* variables depending on Mw_i */
346 double Mw2 = Mw_i*Mw_i;
347 double mHp2[2] = {mySUSY.getMHp()*mySUSY.getMHp(), Mw2}; /* H^+_i = (H^+, G^+) */
348 double cW = Mw_i/Mz;
349 double cW2 = cW*cW;
350 double sW2 = 1.0 - cW2;
351 double sW = sqrt(sW2);
352 double g2sq = e2/sW2; /* g2 squared */
353 double e_2s = e/2.0/sW;
354 double e_sq2s = e/sqrt(2.0)/sW;
355 double e_2sq2s = e_sq2s/2.0;
356 double e2_2s2 = e2/2.0/sW2;
357
358 gslpp::complex PiT_f = gslpp::complex(0.0, 0.0, false);
359 gslpp::complex PiT_sf = gslpp::complex(0.0, 0.0, false);
360 gslpp::complex PiT_ch = gslpp::complex(0.0, 0.0, false);
361 gslpp::complex PiT_WZH = gslpp::complex(0.0, 0.0, false);
362 double a0;
363 gslpp::complex b0, b00;
364
365 /* SM fermion loops */
366 gslpp::complex cV_Wen = e_2sq2s;
367 gslpp::complex cA_Wen = e_2sq2s;
368 gslpp::complex cV_Wne = cV_Wen.conjugate();
369 gslpp::complex cA_Wne = cA_Wen.conjugate();
370 gslpp::complex cV_Wdu = e_2sq2s; /* no CKM */
371 gslpp::complex cA_Wdu = e_2sq2s; /* no CKM */
372 gslpp::complex cV_Wud = cV_Wdu.conjugate();
373 gslpp::complex cA_Wud = cA_Wdu.conjugate();
374 for (int I=0; I<3; ++I) {
375 /* leptons */
376 PiT_f += FA(mu, p2, m_l[I], 0.0, cV_Wen, cV_Wne, cA_Wen, cA_Wne);
377
378 /* quarks (no CKM) */
379 PiT_f += Nc*FA(mu, p2, m_d[I], m_u[I], cV_Wdu, cV_Wud, cA_Wdu, cA_Wud);
380 }
381
382 /* slepton loops */
383 gslpp::complex VWsnL_In, VWWsnsn_II, VWWLL_nn;
384 for (int I=0; I<3; ++I) { /* I=0-2 for left-handed sneutrinos */
385 for (int n=0; n<6; ++n) {
386 VWsnL_In = gslpp::complex(0.0, 0.0, false);
387 for (int J=0; J<3; ++J) /* sum over left-handed sleptons */
388 VWsnL_In += e_sq2s*Zne(J,I)*ZL(J,n);
389 b00 = PV.B00(mu2, p2, Msn2[I], Mse2[n]);
390 PiT_sf += 4.0*VWsnL_In.abs2()*b00;
391 }
392 VWWsnsn_II = e2_2s2;
393 a0 = PV.A0(mu2, Msn2[I]);
394 PiT_sf += VWWsnsn_II*a0;
395 }
396 for (int n=0; n<6; ++n) {
397 VWWLL_nn = gslpp::complex(0.0, 0.0, false);
398 for (int I=0; I<3; ++I) /* sum over left-handed sleptons */
399 VWWLL_nn += e2_2s2*ZL(I,n)*ZL(I,n).conjugate();
400 a0 = PV.A0(mu2, Mse2[n]);
401 PiT_sf += VWWLL_nn*a0;
402 }
403
404 /* squark loops (no CKM) */
405 gslpp::complex VWDU_nm, VWWDD_nn, VWWUU_nn;
406 for (int n=0; n<6; ++n) {
407 for (int m=0; m<6; ++m) {
408 VWDU_nm = gslpp::complex(0.0, 0.0, false);
409 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
410 VWDU_nm += e_sq2s*ZD(I,n).conjugate()*ZU(I,m).conjugate();
411 b00 = PV.B00(mu2, p2, Msd2[n], Msu2[m]);
412 PiT_sf += 4.0*Nc*VWDU_nm.abs2()*b00;
413 }
414 VWWDD_nn = gslpp::complex(0.0, 0.0, false);
415 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
416 VWWDD_nn += e2_2s2*ZD(I,n)*ZD(I,n).conjugate();
417 a0 = PV.A0(mu2, Msd2[n]);
418 PiT_sf += Nc*VWWDD_nn*a0;
419 VWWUU_nn = gslpp::complex(0.0, 0.0, false);
420 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
421 VWWUU_nn += e2_2s2*ZU(I,n).conjugate()*ZU(I,n);
422 a0 = PV.A0(mu2, Msu2[n]);
423 PiT_sf += Nc*VWWUU_nn*a0;
424 }
425
426 /* chargino - neutralino loops */
427 gslpp::complex cV_Wij, cV_Wji, cA_Wij, cA_Wji;
428 for (int i=0; i<2; ++i)
429 for (int j=0; j<4; ++j) {
430 /* W^+ + neutralino(j) -> chi^+(i) */
431 /* W^+ + chi^-(i) -> neutralino(j) */
432 cV_Wji = - e_2s*( ZN(1,j)*Zp(0,i).conjugate()
433 - ZN(3,j)*Zp(1,i).conjugate()/sqrt(2.0)
434 + ZN(1,j).conjugate()*Zm(0,i)
435 + ZN(2,j).conjugate()*Zm(1,i)/sqrt(2.0) );
436 cV_Wij = cV_Wji.conjugate();
437 cA_Wji = - e_2s*( ZN(1,j)*Zp(0,i).conjugate()
438 - ZN(3,j)*Zp(1,i).conjugate()/sqrt(2.0)
439 - ZN(1,j).conjugate()*Zm(0,i)
440 - ZN(2,j).conjugate()*Zm(1,i)/sqrt(2.0) );
441 cA_Wij = cA_Wji.conjugate();
442 PiT_ch += FA(mu, p2, mC[i], mN[j], cV_Wij, cV_Wji, cA_Wij, cA_Wji);
443 }
444
445 /* Higgs loops */
446 double AM_ij;
447 for (int i=0; i<2; ++i) {
448 for (int j=0; j<2; ++j) {
449 AM_ij = ZR(0,i)*ZH(0,j) - ZR(1,i)*ZH(1,j);
450 b00 = PV.B00(mu2, p2, mH02[i], mHp2[j]);
451 PiT_WZH += g2sq*AM_ij*AM_ij*b00;
452 }
453 b00 = PV.B00(mu2, p2, mH02[i+2], mHp2[i]);
454 PiT_WZH += g2sq*b00;
455 }
456 for (int i=0; i<4; ++i) {
457 a0 = PV.A0(mu2, mH02[i]);
458 PiT_WZH += g2sq/4.0*a0;
459 }
460 for (int i=0; i<2; ++i) {
461 a0 = PV.A0(mu2, mHp2[i]);
462 PiT_WZH += g2sq/2.0*a0;
463 }
464
465 /* photon - charged-Goldstone-boson loops */
466 b0 = PV.B0(mu2, p2, Mw2, 0.0);
467 PiT_WZH += - e2*Mw2*b0;
468
469 /* W-boson - Higgs loops */
470 double CR_i;
471 for (int i=0; i<2; ++i) {
472 CR_i = mySUSY.v1()*ZR(0,i) + mySUSY.v2()*ZR(1,i);
473 b0 = PV.B0(mu2, p2, Mw2, mH02[i]);
474 /* Mw^2/v^2 is substituted for g2^2/4 compared to the expression in the
475 * paper, in order to ensure the cancellation of the UV divergences in
476 * the case where Mw is not the tree-level value. */
477 PiT_WZH += - g2sq*Mw2/mySUSY.v()/mySUSY.v()*CR_i*CR_i*b0;
478 }
479
480 /* Z-boson - charged-Goldstone-boson loops */
481 b0 = PV.B0(mu2, p2, Mw2, Mz*Mz);
482 PiT_WZH += - e2*sW2*Mz*Mz*b0;
483
484 /* gauge-boson loops */
485 a0 = PV.A0(mu2, Mw2);
486 b0 = PV.B0(mu2, p2, Mw2, Mz*Mz);
487 b00 = PV.B00(mu2, p2, Mz*Mz, Mw2);
488 PiT_WZH += g2sq*cW2*(2.0*PV.A0(mu2, Mz*Mz) - a0
489 + (4.0*p2 + Mz*Mz + Mw2)*b0 + 8.0*b00);
490 //
491 PiT_WZH += 3.0*g2sq*a0;
492 //
493 b0 = PV.B0(mu2, p2, Mw2, 0.0);
494 b00 = PV.B00(mu2, p2, Mw2, 0.0);
495 PiT_WZH += e2*((4.0*p2 + Mw2)*b0 - a0 + 8.0*b00);
496
497 /* Sum of all contributions */
498 gslpp::complex PiT = PiT_f + PiT_sf + PiT_ch + PiT_WZH;
499
500 return ( PiT/16.0/M_PI/M_PI );
501}
502
503gslpp::complex EWSUSY::PiT_AZ(const double mu, const double p2, const double Mw_i) const
504{
505 double mu2 = mu*mu;
506 double e2 = 4.0*M_PI*mySUSY.getAle();
507 double e = sqrt(e2);
508 double Mz = mySUSY.getMz();
509 double Nc = mySUSY.getNc();
510
511 /* variables depending on Mw_i */
512 double Mw2 = Mw_i*Mw_i;
513 double mHp2[2] = {mySUSY.getMHp()*mySUSY.getMHp(), Mw2}; /* H^+_i = (H^+, G^+) */
514 double cW = Mw_i/Mz;
515 double cW2 = cW*cW;
516 double sW2 = 1.0 - cW2;
517 double sW = sqrt(sW2);
518 double g2 = e/sW;
519 double e_4sc = e/4.0/sW/cW;
520 double e_2sc = 2.0*e_4sc;
521 double e2_sc = e2/sW/cW;
522
523 gslpp::complex PiT_f = gslpp::complex(0.0, 0.0, false);
524 gslpp::complex PiT_sf = gslpp::complex(0.0, 0.0, false);
525 gslpp::complex PiT_ch = gslpp::complex(0.0, 0.0, false);
526 gslpp::complex PiT_WZH = gslpp::complex(0.0, 0.0, false);
527 double a0;
528 gslpp::complex b0, b00;
529
530 /* SM fermion loops */
531 gslpp::complex cV_Aee = - e;
532 gslpp::complex cA_Aee = 0.0;
533 gslpp::complex cV_Add = - e/3.0;
534 gslpp::complex cA_Add = 0.0;
535 gslpp::complex cV_Auu = 2.0/3.0*e;
536 gslpp::complex cA_Auu = 0.0;
537 gslpp::complex cV_Zee = - e_4sc*(1.0 - 4.0*sW2);
538 gslpp::complex cA_Zee = - e_4sc;
539 gslpp::complex cV_Zdd = - e_4sc*(1.0 - 4.0/3.0*sW2);
540 gslpp::complex cA_Zdd = - e_4sc;
541 gslpp::complex cV_Zuu = e_4sc*(1.0 - 8.0/3.0*sW2);
542 gslpp::complex cA_Zuu = e_4sc;
543 for (int I=0; I<3; ++I) {
544 /* charged leptons */
545 PiT_f += FA(mu, p2, m_l[I], m_l[I], cV_Aee, cV_Zee, cA_Aee, cA_Zee);
546
547 /* down-type quarks */
548 PiT_f += Nc*FA(mu, p2, m_d[I], m_d[I], cV_Add, cV_Zdd, cA_Add, cA_Zdd);
549
550 /* up-type quarks */
551 PiT_f += Nc*FA(mu, p2, m_u[I], m_u[I], cV_Auu, cV_Zuu, cA_Auu, cA_Zuu);
552 }
553
554 /* charged-slepton loops */
555 gslpp::complex VZLL_nn, VAZLL_nn;
556 for (int n=0; n<6; ++n) {
557 VZLL_nn = gslpp::complex(0.0, 0.0, false);
558 for (int I=0; I<3; ++I) /* sum over left-handed sleptons */
559 VZLL_nn += - e_2sc*ZL(I,n)*ZL(I,n).conjugate();
560 VZLL_nn += - e_2sc*(- 2.0*sW2);
561 b00 = PV.B00(mu2, p2, Mse2[n], Mse2[n]);
562 /* typo in the paper: e^2 --> e */
563 PiT_sf += - 4.0*e*VZLL_nn*b00;
564
565 VAZLL_nn = gslpp::complex(0.0, 0.0, false);
566 for (int I=0; I<3; ++I) /* sum over left-handed sleptons */
567 VAZLL_nn += e2_sc*ZL(I,n)*ZL(I,n).conjugate();
568 VAZLL_nn += e2_sc*(- 2.0*sW2);
569 a0 = PV.A0(mu2, Mse2[n]);
570 PiT_sf += VAZLL_nn*a0;
571 }
572
573 /* down-type squark loops */
574 gslpp::complex VZDD_nn, VAZDD_nn;
575 for (int n=0; n<6; ++n) {
576 VZDD_nn = gslpp::complex(0.0, 0.0, false);
577 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
578 VZDD_nn += - e_2sc*ZD(I,n)*ZD(I,n).conjugate();
579 VZDD_nn += - e_2sc*(- 2.0/3.0*sW2);
580 b00 = PV.B00(mu2, p2, Msd2[n], Msd2[n]);
581 /* typo in the paper: e^2 --> e */
582 PiT_sf += - 4.0*e*Nc/3.0*VZDD_nn*b00;
583
584 VAZDD_nn = gslpp::complex(0.0, 0.0, false);
585 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
586 VAZDD_nn += e2_sc/3.0*ZD(I,n)*ZD(I,n).conjugate();
587 VAZDD_nn += e2_sc/3.0*(- 2.0/3.0*sW2);
588 a0 = PV.A0(mu2, Msd2[n]);
589 PiT_sf += Nc*VAZDD_nn*a0;
590 }
591
592 /* up-type squark loops */
593 gslpp::complex VZUU_nn, VAZUU_nn;
594 for (int n=0; n<6; ++n) {
595 VZUU_nn = gslpp::complex(0.0, 0.0, false);
596 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
597 VZUU_nn += e_2sc*ZU(I,n).conjugate()*ZU(I,n);
598 VZUU_nn += e_2sc*(- 4.0/3.0*sW2);
599 b00 = PV.B00(mu2, p2, Msu2[n], Msu2[n]);
600 /* typo in the paper: e^2 --> e */
601 PiT_sf += 4.0*e*Nc*2.0/3.0*VZUU_nn*b00;
602
603 VAZUU_nn = gslpp::complex(0.0, 0.0, false);
604 for (int I=0; I<3; ++I) /* sum over left-handed squarks */
605 VAZUU_nn += 2.0*e2_sc/3.0*ZU(I,n).conjugate()*ZU(I,n);
606 VAZUU_nn += 2.0*e2_sc/3.0*(- 4.0/3.0*sW2);
607 a0 = PV.A0(mu2, Msu2[n]);
608 PiT_sf += Nc*VAZUU_nn*a0;
609 }
610
611 /* chargino loops */
612 gslpp::complex cV_Aii = e;
613 gslpp::complex cA_Aii = 0.0;
614 gslpp::complex cV_Zii, cA_Zii;
615 for (int i=0; i<2; ++i) {
616 cV_Zii = e_4sc*( Zp(0,i).conjugate()*Zp(0,i)
617 + Zm(0,i)*Zm(0,i).conjugate() + 2.0*(cW2 - sW2) );
618 cA_Zii = e_4sc*( Zp(0,i).conjugate()*Zp(0,i)
619 - Zm(0,i)*Zm(0,i).conjugate() );
620 PiT_ch += FA(mu, p2, mC[i], mC[i], cV_Aii, cV_Zii, cA_Aii, cA_Zii);
621 }
622
623 /* W-boson - charged-Goldstone-boson loops */
624 b0 = PV.B0(mu2, p2, Mw2, Mw2);
625 PiT_WZH += 2.0*e2*cW*sW*Mz*Mz*b0;
626
627 /* W-boson loops */
628 a0 = PV.A0(mu2, Mw2);
629 //b0 = PV.B0(mu2, p2, Mw2, Mw2); /* same as the above */
630 b00 = PV.B00(mu2, p2, Mw2, Mw2);
631 PiT_WZH += 2.0*e*g2*cW*(2.0*a0 + (2.0*p2 + Mw2)*b0 + 4.0*b00);
632
633 /* charged-Higgs loops */
634 double cot_2thW = (cW2 - sW2)/(2.0*sW*cW);
635 for (int i=0; i<2; ++i) {
636 a0 = PV.A0(mu2, mHp2[i]);
637 b00 = PV.B00(mu2, p2, mHp2[i], mHp2[i]);
638 PiT_WZH += 2.0*e2*cot_2thW*(2.0*b00 + a0);
639 }
640
641 /* Sum of all contributions */
642 gslpp::complex PiT = PiT_f + PiT_sf + PiT_ch + PiT_WZH;
643
644 return ( PiT/16.0/M_PI/M_PI );
645}
646
647gslpp::complex EWSUSY::PiTp_A(const double mu, const double p2, const double Mw_i) const
648{
649 double mu2 = mu*mu;
650 double e2 = 4.0*M_PI*mySUSY.getAle();
651 double e = sqrt(e2);
652 double Nc = mySUSY.getNc();
653
654 /* variables depending on Mw_i */
655 double Mw2 = Mw_i*Mw_i;
656 double mHp2[2] = {mySUSY.getMHp()*mySUSY.getMHp(), Mw2}; /* H^+_i = (H^+, G^+) */
657
658 gslpp::complex PiTp_f = gslpp::complex(0.0, 0.0, false);
659 gslpp::complex PiTp_sf = gslpp::complex(0.0, 0.0, false);
660 gslpp::complex PiTp_ch = gslpp::complex(0.0, 0.0, false);
661 gslpp::complex PiTp_WZH = gslpp::complex(0.0, 0.0, false);
662 gslpp::complex b0, b0p, b00p;
663
664 /* SM fermion loops */
665 gslpp::complex cV_Aee = - e;
666 gslpp::complex cA_Aee = 0.0;
667 gslpp::complex cV_Add = - e/3.0;
668 gslpp::complex cA_Add = 0.0;
669 gslpp::complex cV_Auu = 2.0/3.0*e;
670 gslpp::complex cA_Auu = 0.0;
671 for (int I=0; I<3; ++I) {
672 /* charged leptons */
673 PiTp_f += dFA(mu, p2, m_l[I], m_l[I], cV_Aee, cV_Aee, cA_Aee, cA_Aee);
674
675 /* down-type quarks */
676 PiTp_f += Nc*dFA(mu, p2, m_d[I], m_d[I], cV_Add, cV_Add, cA_Add, cA_Add);
677
678 /* up-type quarks */
679 PiTp_f += Nc*dFA(mu, p2, m_u[I], m_u[I], cV_Auu, cV_Auu, cA_Auu, cA_Auu);
680 }
681
682 /* charged-slepton loops */
683 for (int n=0; n<6; ++n) {
684 b00p = PV.B00p(mu2, p2, Mse2[n], Mse2[n]);
685 PiTp_sf += 4.0*e2*b00p;
686 }
687
688 /* down-type squark loops */
689 for (int n=0; n<6; ++n) {
690 b00p = PV.B00p(mu2, p2, Msd2[n], Msd2[n]);
691 PiTp_sf += 4.0*e2*Nc/3.0/3.0*b00p;
692 }
693
694 /* up-type squark loops */
695 for (int n=0; n<6; ++n) {
696 b00p = PV.B00p(mu2, p2, Msu2[n], Msu2[n]);
697 PiTp_sf += 4.0*e2*Nc*2.0/3.0*2.0/3.0*b00p;
698 }
699
700 /* chargino loops */
701 gslpp::complex cV_Aii = e;
702 gslpp::complex cA_Aii = 0.0;
703 for (int i=0; i<2; ++i)
704 PiTp_ch += dFA(mu, p2, mC[i], mC[i], cV_Aii, cV_Aii, cA_Aii, cA_Aii);
705
706 /* charged-Higgs loops */
707 for (int i=0; i<2; ++i) {
708 b00p = PV.B00p(mu2, p2, mHp2[i], mHp2[i]);
709 PiTp_WZH += 4.0*e2*b00p;
710 }
711
712 /* W-boson loops */
713 /* The Mw_i*Mw_i*b0p term, adding the corresponding contribution from
714 * the W-G loop below, differs from the one in the paper. */
715 b0 = PV.B0(mu2, p2, Mw2, Mw2);
716 b0p = PV.B0p(mu2, p2, Mw2, Mw2);
717 b00p = PV.B00p(mu2, p2, Mw2, Mw2);
718 PiTp_WZH += 2.0*e2*( (2.0*p2 + Mw2)*b0p + 2.0*b0 + 4.0*b00p);
719
720 /* W-boson - charged-Goldstone-boson loop */
721 //b0p = PV.B0p(mu2, p2, Mw2, Mw2); /* Same as the above */
722 PiTp_WZH += - 2.0*e2*Mw2*b0p;
723
724 /* Sum of all contributions */
725 gslpp::complex PiTp = PiTp_f + PiTp_sf + PiTp_ch + PiTp_WZH;
726
727 return ( PiTp/16.0/M_PI/M_PI );
728}
729
730double EWSUSY::PiThat_W_0(const double Mw_i) const
731{
732 /* Renormalization scale (varied for checking the cancellation of UV divergences */
733 double mu = Mw_i * RenormalizationScaleFactor;
734
735 double Mz = mySUSY.getMz();
736 double cW = Mw_i/Mz;
737 double cW2 = cW*cW;
738 double sW2 = 1.0 - cW2;
739 double sW = sqrt(sW2);
740
741 double PiThat = 0.0;
742
743 /* W self-energy */
744 PiThat += PiT_W(mu, 0.0, Mw_i).real();
745
746 /* W-mass counter term */
747 double delMw2 = PiT_W(mu, Mw_i*Mw_i, Mw_i).real();
748 PiThat -= delMw2;
749
750 /* counter term for e: (del e)/e */
751 double dele_over_e = PiTp_A(mu, 0.0, Mw_i).real()/2.0
752 + sW/cW*PiT_AZ(mu, 0.0, Mw_i).real()/Mz/Mz;
753 PiThat += 2.0*Mw_i*Mw_i*dele_over_e;
754
755 /* counter term for sW: (del sW)/sW */
756 double delSw_overSw = - cW2/2.0/sW2
757 *( PiT_W(mu, Mw_i*Mw_i, Mw_i).real()/Mw_i/Mw_i
758 - PiT_Z(mu, Mz*Mz, Mw_i).real()/Mz/Mz );
759 PiThat -= 2.0*Mw_i*Mw_i*delSw_overSw;
760
761 /* remaining counter terms,
762 * usually denoted by 2/(sW*cW)*PiT_AZ(0)/Mz/Mz. */
763 PiThat += - 2.0*Mw_i*Mw_i/(sW*cW)*PiT_AZ(mu, 0.0, Mw_i).real()/Mz/Mz;
764
765 return PiThat;
766}
767
768double EWSUSY::DeltaR_rem_SM(const double Mw_i) const
769{
770 double cW2 = Mw_i*Mw_i/mySUSY.getMz()/mySUSY.getMz();
771 double sW2 = 1.0 - cW2;
772
773 /* renormalized vertex corrections + box */
774 return ( mySUSY.getAle()/4.0/M_PI/sW2
775 *(6.0 + (7.0 - 4.0*sW2)/2.0/sW2*log(cW2)) );
776}
777
778double EWSUSY::DeltaR_boxLL_SUSY(const double Mw_i) const
779{
780 int M = 1; // MU
781 int N = 0; // ELECTRON
782 int J = 1; // NEUTRINO_2
783 int I = 0; // NEUTRINO_1
784
785 gslpp::complex a11 = gslpp::complex(0.0, 0.0, false);
786 gslpp::complex a12 = gslpp::complex(0.0, 0.0, false);
787
788 /* charged-lepton - sneutrino - chargino - neutralino loop */
789 for (int k=0; k<6; ++k)
790 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
791 for (int i=0; i<2; ++i)
792 for (int j=0; j<4; ++j) {
793 gslpp::complex FF = F(sqrt(Mse2[k]), sqrt(Msn2[K]), mC[i], mN[j]);
794 a11 += 0.5
795 *L_esnC(M, K, i, Mw_i)
796 *L_nLC(I, k, i, Mw_i)
797 *L_nsnN(J, K, j, Mw_i).conjugate()
798 *L_eLN(N, k, j, Mw_i).conjugate()
799 *mC[i]*mN[j]*FF;
800 a11 += 0.5
801 *L_eLN(M, k, j, Mw_i)
802 *L_nsnN(I, K, j, Mw_i)
803 *L_nLC(J, k, i, Mw_i).conjugate()
804 *L_esnC(N, K, i, Mw_i).conjugate()
805 *mC[i]*mN[j]*FF;
806 }
807
808 /* charged-lepton - charged-lepton - chargino - neutralino loop */
809 for (int k=0; k<6; ++k)
810 for (int l=0; l<6; ++l)
811 for (int i=0; i<2; ++i)
812 for (int j=0; j<4; ++j) {
813 a11 += L_eLN(M, k, j, Mw_i)
814 *L_nLC(J, k, i, Mw_i).conjugate()
815 *L_nLC(I, l, i, Mw_i)
816 *L_eLN(N, l, j, Mw_i).conjugate()
817 *H(sqrt(Mse2[k]), sqrt(Mse2[l]), mC[i], mN[j]);
818 }
819
820 /* sneutrino - sneutrino - chargino - neutralino loop */
821 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
822 for (int L=0; L<3; ++L) /* L=0-2 for left-handed sneutrinos */
823 for (int i=0; i<2; ++i)
824 for (int j=0; j<4; ++j) {
825 a11 += L_esnC(M, K, i, Mw_i)
826 *L_nsnN(J, K, j, Mw_i).conjugate()
827 *L_nsnN(I, L, j, Mw_i)
828 *L_esnC(N, L, i, Mw_i).conjugate()
829 *H(sqrt(Msn2[K]), sqrt(Msn2[L]), mC[i], mN[j]);
830 }
831
832 /* charged-lepton - sneutrino - chargino - chargino loop */
833 for (int k=0; k<6; ++k)
834 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
835 for (int i=0; i<2; ++i)
836 for (int j=0; j<2; ++j) {
837 a12 += 0.5
838 *L_esnC(M, K, i, Mw_i)
839 *L_nLC(I, k, i, Mw_i)
840 *L_nLC(J, k, j, Mw_i).conjugate()
841 *L_esnC(N, K, j, Mw_i).conjugate()
842 *mC[i]*mC[j]*F(sqrt(Mse2[k]), sqrt(Msn2[K]), mC[i], mC[j]);
843 }
844
845 /* charged-lepton - sneutrino - neutralino - neutralino loop */
846 for (int k=0; k<6; ++k)
847 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
848 for (int i=0; i<4; ++i)
849 for (int j=0; j<4; ++j) {
850 a12 += 0.5
851 *L_eLN(M, k, i, Mw_i)
852 *L_nsnN(I, K, i, Mw_i)
853 *L_nsnN(J, K, j, Mw_i).conjugate()
854 *L_eLN(N, k, j, Mw_i).conjugate()
855 *mN[i]*mN[j]*F(sqrt(Mse2[k]), sqrt(Msn2[K]), mN[i], mN[j]);
856 a12 += L_eLN(M, k, i, Mw_i)
857 *L_nsnN(J, K, i, Mw_i).conjugate()
858 *L_nsnN(I, K, j, Mw_i)
859 *L_eLN(N, k, j, Mw_i).conjugate()
860 *H(sqrt(Mse2[k]), sqrt(Msn2[K]), mN[i], mN[j]);
861 }
862
863 gslpp::complex a1 = (a11 + a12)/16.0/M_PI/M_PI;
864
865 double sW2 = 1.0 - Mw_i*Mw_i/mySUSY.getMz()/mySUSY.getMz();
866 return ( - sW2*Mw_i*Mw_i/2.0/M_PI/mySUSY.getAle()*a1.real() );
867}
868
869double EWSUSY::DeltaR_boxLR_SUSY(const double Mw_i) const
870{
871 int M = 1; // MU
872 int N = 0; // ELECTRON
873 int J = 1; // NEUTRINO_2
874 int I = 0; // NEUTRINO_1
875
876 gslpp::complex a21 = gslpp::complex(0.0, 0.0, false);
877 gslpp::complex a22 = gslpp::complex(0.0, 0.0, false);
878
879 /* charged-lepton - sneutrino - chargino - neutralino loop */
880 for (int k=0; k<6; ++k)
881 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
882 for (int i=0; i<2; ++i)
883 for (int j=0; j<4; ++j) {
884 gslpp::complex HH = H(sqrt(Mse2[k]), sqrt(Msn2[K]), mC[i], mN[j]);
885 a21 += - 2.0
886 *R_esnC(M, K, i)
887 *L_nLC(I, k, i, Mw_i)
888 *L_nsnN(J, K, j, Mw_i).conjugate()
889 *R_eLN(N, k, j, Mw_i).conjugate()
890 *HH;
891 a21 += - 2.0
892 *R_eLN(M, k, j, Mw_i)
893 *L_nsnN(I, K, j, Mw_i)
894 *L_nLC(J, k, i, Mw_i).conjugate()
895 *R_esnC(N, K, i).conjugate()
896 *HH;
897 }
898
899 /* charged-lepton - charged-lepton - chargino - neutralino loop */
900 for (int k=0; k<6; ++k)
901 for (int l=0; l<6; ++l)
902 for (int i=0; i<2; ++i)
903 for (int j=0; j<4; ++j) {
904 a21 += - 2.0
905 *R_eLN(M, k, j, Mw_i)
906 *L_nLC(J, k, i, Mw_i).conjugate()
907 *L_nLC(I, l, i, Mw_i)
908 *R_eLN(N, l, j, Mw_i).conjugate()
909 *H(sqrt(Mse2[k]), sqrt(Mse2[l]), mC[i], mN[j]);
910 }
911
912 /* sneutrino - sneutrino - chargino - neutralino loop */
913 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
914 for (int L=0; L<3; ++L) /* L=0-2 for left-handed sneutrinos */
915 for (int i=0; i<2; ++i)
916 for (int j=0; j<4; ++j) {
917 a21 += - 2.0
918 *R_esnC(M, K, i)
919 *L_nsnN(J, K, j, Mw_i).conjugate()
920 *L_nsnN(I, L, j, Mw_i)
921 *R_esnC(N, L, i).conjugate()
922 *H(sqrt(Msn2[K]), sqrt(Msn2[L]), mC[i], mN[j]);
923 }
924
925 /* charged-lepton - sneutrino - neutralino - neutralino loop */
926 for (int k=0; k<6; ++k)
927 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
928 for (int i=0; i<4; ++i)
929 for (int j=0; j<4; ++j) {
930 a22 += R_eLN(M, k, i, Mw_i)
931 *L_nsnN(J, K, i, Mw_i).conjugate()
932 *L_nsnN(I, K, j, Mw_i)
933 *R_eLN(N, k, j, Mw_i).conjugate()
934 *mN[i]*mN[j]*F(sqrt(Mse2[k]), sqrt(Msn2[K]), mN[i], mN[j]);
935 a22 += 2.0
936 *R_eLN(M, k, i, Mw_i)
937 *L_nsnN(I, K, i, Mw_i)
938 *L_nsnN(J, K, j, Mw_i).conjugate()
939 *R_eLN(N, k, j, Mw_i).conjugate()
940 *H(sqrt(Mse2[k]), sqrt(Msn2[K]), mN[i], mN[j]);
941 }
942
943 /* charged-lepton - sneutrino - chargino - chargino loop */
944 for (int k=0; k<6; ++k)
945 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
946 for (int i=0; i<2; ++i)
947 for (int j=0; j<2; ++j) {
948 a22 += 2.0
949 *R_esnC(M, K, i)
950 *L_nLC(I, k, i, Mw_i)
951 *L_nLC(J, k, j, Mw_i).conjugate()
952 *R_esnC(N, K, j).conjugate()
953 *H(sqrt(Mse2[k]), sqrt(Msn2[K]), mC[i], mC[j]);
954 }
955
956 gslpp::complex a2 = (a21 + a22)/16.0/M_PI/M_PI;
957
958 double sW2 = 1.0 - Mw_i*Mw_i/mySUSY.getMz()/mySUSY.getMz();
959 return ( sW2*Mw_i*Mw_i/4.0/M_PI/mySUSY.getAle()*a2.real() );
960}
961
962gslpp::complex EWSUSY::v(const double mu, const QCD::lepton M,
963 const QCD::lepton J, const double Mw_i) const
964{
965 int intM, intJ;
966 switch (M) {
967 case StandardModel::ELECTRON:
968 case StandardModel::MU:
969 case StandardModel::TAU:
970 intM = ((int)M - StandardModel::ELECTRON)/2;
971 break;
972 default:
973 throw std::runtime_error("EWSUSY::v(): Wrong argument!");
974 }
975 switch (J) {
976 case StandardModel::NEUTRINO_1:
977 case StandardModel::NEUTRINO_2:
978 case StandardModel::NEUTRINO_3:
979 intJ = ((int)J - StandardModel::NEUTRINO_1)/2;
980 break;
981 default:
982 throw std::runtime_error("EWSUSY::v(): Wrong argument!");
983 }
984
985 gslpp::complex v = gslpp::complex(0.0, 0.0, false);
986 gslpp::complex b0, ff;
987 gslpp::complex CL_ji, CR_ji; /* chargino-neutralino-W couplings */
988
989 /* charged-slepton - chargino - neutralino loops */
990 for (int k=0; k<6; ++k)
991 for (int j=0; j<4; ++j)
992 for (int i=0; i<2; ++i) {
993 CL_ji = ZN(1,j)*Zp(0,i).conjugate()
994 - ZN(3,j)*Zp(1,i).conjugate()/sqrt(2.0);
995 CR_ji = ZN(1,j).conjugate()*Zm(0,i)
996 + ZN(2,j).conjugate()*Zm(1,i)/sqrt(2.0);
997 b0 = PV.B0(mu*mu, 0.0, mC[i]*mC[i], mN[j]*mN[j]);
998 ff = f(sqrt(Mse2[k]), mC[i], mN[j]);
999 v += L_nLC(intJ, k, i, Mw_i).conjugate()*L_eLN(intM, k, j, Mw_i)
1000 *( sqrt(2.0)*CL_ji*mC[i]*mN[j]*ff
1001 - CR_ji/sqrt(2.0)*(b0 - 0.5 + Mse2[k]*ff) );
1002 }
1003
1004 /* sneutrino - neutralino - chargino loops */
1005 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
1006 for (int j=0; j<4; ++j)
1007 for (int i=0; i<2; ++i) {
1008 CL_ji = ZN(1,j)*Zp(0,i).conjugate()
1009 - ZN(3,j)*Zp(1,i).conjugate()/sqrt(2.0);
1010 CR_ji = ZN(1,j).conjugate()*Zm(0,i)
1011 + ZN(2,j).conjugate()*Zm(1,i)/sqrt(2.0);
1012 b0 = PV.B0(mu*mu, 0.0, mC[i]*mC[i], mN[j]*mN[j]);
1013 ff = f(sqrt(Msn2[K]), mC[i], mN[j]);
1014 v += L_nsnN(intJ, K, j, Mw_i).conjugate()*L_esnC(intM, K, i, Mw_i)
1015 *( - sqrt(2.0)*CR_ji*mC[i]*mN[j]*ff
1016 + CL_ji/sqrt(2.0)*(b0 - 0.5 + Msn2[K]*ff) );
1017 }
1018
1019 /* sneutrino - charged-slepton - neutralino loops */
1020 gslpp::matrix<gslpp::complex> ZneT_ZL = Zne.transpose()*ZL;
1021 for (int i=0; i<6; ++i)
1022 for (int j=0; j<4; ++j)
1023 for (int K=0; K<3; ++K) { /* K=0-2 for left-handed sneutrinos */
1024 b0 = PV.B0(mu*mu, 0.0, Mse2[i], Msn2[K]);
1025 ff = f(mN[j], sqrt(Mse2[i]), sqrt(Msn2[K]));
1026 v += 0.5*L_nsnN(intJ, K, j, Mw_i).conjugate()*L_eLN(intM, i, j, Mw_i)
1027 *ZneT_ZL(K, i).conjugate()*(b0 + 0.5 + mN[j]*mN[j]*ff);
1028 }
1029
1030 return ( v/16.0/M_PI/M_PI );
1031}
1032
1033gslpp::complex EWSUSY::delta_v(const double mu, const QCD::lepton M,
1034 const QCD::lepton J, const double Mw_i) const
1035{
1036 int intM, intJ;
1037 switch (M) {
1038 case StandardModel::ELECTRON:
1039 case StandardModel::MU:
1040 case StandardModel::TAU:
1041 intM = ((int)M - StandardModel::ELECTRON)/2;
1042 break;
1043 default:
1044 throw std::runtime_error("EWSUSY::delta_v(): Wrong argument!");
1045 }
1046 switch (J) {
1047 case StandardModel::NEUTRINO_1:
1048 case StandardModel::NEUTRINO_2:
1049 case StandardModel::NEUTRINO_3:
1050 intJ = ((int)J - StandardModel::NEUTRINO_1)/2;
1051 break;
1052 default:
1053 throw std::runtime_error("EWSUSY::delta_v(): Wrong argument!");
1054 }
1055
1056 gslpp::complex delv = gslpp::complex(0.0, 0.0, false);
1057 double muIR = mu; /* fictional scale, since B0p(0,m1^2,m2^2) is IR finite */
1058 gslpp::complex b0p, b0;
1059
1060 /* charged-slepton - neutralino loops */
1061 for (int k=0; k<6; ++k)
1062 for (int j=0; j<4; ++j) {
1063 b0p = PV.B0p(muIR*muIR, 0.0, Mse2[k], mN[j]*mN[j]);
1064 b0 = PV.B0(mu*mu, 0.0, Mse2[k], mN[j]*mN[j]);
1065 delv += 0.5*L_eLN(intM, k, j, Mw_i)*L_eLN(intJ, k, j, Mw_i).conjugate()
1066 *( (Mse2[k] - mN[j]*mN[j])*b0p - b0 );
1067 }
1068
1069 /* sneutrino - chargino loops */
1070 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
1071 for (int i=0; i<2; ++i) {
1072 b0p = PV.B0p(muIR*muIR, 0.0, Msn2[K], mC[i]*mC[i]);
1073 b0 = PV.B0(mu*mu, 0.0, Msn2[K], mC[i]*mC[i]);
1074 delv += 0.5*L_esnC(intM, K, i, Mw_i)*L_esnC(intJ, K, i, Mw_i).conjugate()
1075 *( (Msn2[K] - mC[i]*mC[i])*b0p - b0 );
1076 }
1077
1078 return ( delv/16.0/M_PI/M_PI );
1079}
1080
1081double EWSUSY::DeltaR_vertex_SUSY(const double Mw_i) const
1082{
1083 /* Renormalization scale (varied for checking the cancellation of UV divergences */
1084 double mu = Mw_i * RenormalizationScaleFactor;
1085
1086 return ( v(mu, mySUSY.ELECTRON, mySUSY.NEUTRINO_1, Mw_i).real()
1087 + delta_v(mu, mySUSY.ELECTRON, mySUSY.NEUTRINO_1, Mw_i).real()
1088 + v(mu, mySUSY.MU, mySUSY.NEUTRINO_2, Mw_i).real()
1089 + delta_v(mu, mySUSY.MU, mySUSY.NEUTRINO_2, Mw_i).real() );
1090}
1091
1092gslpp::complex EWSUSY::Sigma_nu_0(const double mu, const QCD::lepton I,
1093 const QCD::lepton J, const double Mw_i) const
1094{
1095 int intI, intJ;
1096 switch (I) {
1097 case StandardModel::NEUTRINO_1:
1098 case StandardModel::NEUTRINO_2:
1099 case StandardModel::NEUTRINO_3:
1100 intI = ((int)I - StandardModel::NEUTRINO_1)/2;
1101 break;
1102 default:
1103 throw std::runtime_error("EWSUSY::Sigma_nu(): Wrong argument!");
1104 }
1105 switch (J) {
1106 case StandardModel::NEUTRINO_1:
1107 case StandardModel::NEUTRINO_2:
1108 case StandardModel::NEUTRINO_3:
1109 intJ = ((int)J - StandardModel::NEUTRINO_1)/2;
1110 break;
1111 default:
1112 throw std::runtime_error("EWSUSY::Sigma_nu(): Wrong argument!");
1113 }
1114
1115 gslpp::complex Sigma = gslpp::complex(0.0, 0.0, false);
1116 double muIR = mu; /* fictional scale, since B0p(0,m1,m2) is IR finite */
1117 gslpp::complex b0p, b0;
1118
1119 /* charged-slepton - chargino loops */
1120 for (int k=0; k<6; ++k)
1121 for (int i=0; i<2; ++i) {
1122 b0p = PV.B0p(muIR*muIR, 0.0, Mse2[k], mC[i]*mC[i]);
1123 b0 = PV.B0(mu*mu, 0.0, Mse2[k], mC[i]*mC[i]);
1124 Sigma += 0.5*L_nLC(intI, k, i, Mw_i)*L_nLC(intJ, k, i, Mw_i).conjugate()
1125 *( (Mse2[k] - mC[i]*mC[i])*b0p - b0 );
1126 }
1127
1128 /* sneutrino - neutralino loops */
1129 for (int K=0; K<3; ++K) /* K=0-2 for left-handed sneutrinos */
1130 for (int j=0; j<4; ++j) {
1131 b0p = PV.B0p(muIR*muIR, 0.0, Msn2[K], mN[j]*mN[j]);
1132 b0 = PV.B0(mu*mu, 0.0, Msn2[K], mN[j]*mN[j]);
1133 Sigma += 0.5*L_nsnN(intI, K, j, Mw_i)*L_nsnN(intJ, K, j, Mw_i).conjugate()
1134 *( (Msn2[K] - mN[j]*mN[j])*b0p - b0 );
1135 }
1136
1137 return ( Sigma/16.0/M_PI/M_PI );
1138}
1139
1140double EWSUSY::DeltaR_neutrino_SUSY(const double Mw_i) const
1141{
1142 /* Renormalization scale (varied for checking the cancellation of UV divergences */
1143 double mu = Mw_i * RenormalizationScaleFactor;
1144
1145 return ( ( Sigma_nu_0(mu, mySUSY.NEUTRINO_1, mySUSY.NEUTRINO_1, Mw_i).real()
1146 - delta_v(mu, mySUSY.ELECTRON, mySUSY.NEUTRINO_1, Mw_i).real()
1147 + Sigma_nu_0(mu, mySUSY.NEUTRINO_2, mySUSY.NEUTRINO_2, Mw_i).real()
1148 - delta_v(mu, mySUSY.MU, mySUSY.NEUTRINO_2, Mw_i).real() )/2.0 );
1149}
1150
1151double EWSUSY::DeltaR_TOTAL_EW1(const double Mw_i) const
1152{
1153 double DeltaR = 0.0;
1154
1155 /* SM+SUSY renormalized W self energy */
1156 DeltaR += - PiThat_W_0(Mw_i)/Mw_i/Mw_i;
1157
1158 /* SM renormalized vertex + box */
1159 DeltaR += DeltaR_rem_SM(Mw_i);
1160
1161 /* SUSY box corrections */
1162 DeltaR += DeltaR_boxLL_SUSY(Mw_i);
1163 DeltaR += DeltaR_boxLR_SUSY(Mw_i);
1164
1165 /* SUSY renormalized vertex corrections */
1166 DeltaR += DeltaR_vertex_SUSY(Mw_i);
1167
1168 /* SUSY renormalized neutrino wave function */
1169 DeltaR += DeltaR_neutrino_SUSY(Mw_i);
1170
1171 /* Debug */
1172 //std::cout << "MSSM WSE = " << - PiThat_W_0(Mw_i)/Mw_i/Mw_i << std::endl;
1173 //std::cout << "SM VC+Box = " << DeltaR_rem_SM(Mw_i) << std::endl;
1174 //std::cout << "SUSY BoxLL = " << DeltaR_boxLL_SUSY(Mw_i) << std::endl;
1175 //std::cout << "SUSY BoxLR = " << DeltaR_boxLR_SUSY(Mw_i) << std::endl;
1176 //std::cout << "SUSY VC = " << DeltaR_vertex_SUSY(Mw_i) << std::endl;
1177 //std::cout << "SUSY nuSE = " << DeltaR_neutrino_SUSY(Mw_i) << std::endl;
1178
1179 return DeltaR;
1180}
1181
1182double EWSUSY::DeltaAlphaL5q_SM_EW1() const
1183{
1184 /* Renormalization scale (varied for checking the cancellation of UV divergences */
1185 double mu = mySUSY.getMz() * RenormalizationScaleFactor;
1186
1187 double Mz2 = mySUSY.getMz()*mySUSY.getMz();
1188 double e = sqrt(4.0*M_PI*mySUSY.getAle());
1189 double Nc = mySUSY.getNc();
1190
1191 double DelA_l = 0.0, DelA_d = 0.0, DelA_u = 0.0;
1192
1193 /* SM fermion loops */
1194 gslpp::complex cV_Aee = - e;
1195 gslpp::complex cA_Aee = 0.0;
1196 gslpp::complex cV_Add = - e/3.0;
1197 gslpp::complex cA_Add = 0.0;
1198 gslpp::complex cV_Auu = 2.0/3.0*e;
1199 gslpp::complex cA_Auu = 0.0;
1200 for (int I=0; I<3; ++I) {
1201 /* charged leptons */
1202 DelA_l += FA(mu, Mz2, m_l[I], m_l[I], cV_Aee, cV_Aee, cA_Aee, cA_Aee).real()/Mz2;
1203 DelA_l -= dFA(mu, 0.0, m_l[I], m_l[I], cV_Aee, cV_Aee, cA_Aee, cA_Aee).real();
1204
1205 /* down-type quarks */
1206 DelA_d += Nc*FA(mu, Mz2, m_d[I], m_d[I], cV_Add, cV_Add, cA_Add, cA_Add).real()/Mz2;
1207 DelA_d -= Nc*dFA(mu, 0.0, m_d[I], m_d[I], cV_Add, cV_Add, cA_Add, cA_Add).real();
1208
1209 /* up-type quarks, not including top quark */
1210 if (I!=3) {
1211 DelA_u += Nc*FA(mu, Mz2, m_u[I], m_u[I], cV_Auu, cV_Auu, cA_Auu, cA_Auu).real()/Mz2;
1212 DelA_u -= Nc*dFA(mu, 0.0, m_u[I], m_u[I], cV_Auu, cV_Auu, cA_Auu, cA_Auu).real();
1213 }
1214 }
1215
1216 /* Debug */
1217 //std::cout << "EWSUSY(l) " << DelA_l/16.0/M_PI/M_PI << std::endl;
1218 //std::cout << "EWSM(l) " << getMyOneLoopEW()->DeltaAlpha_l(Mz2) << std::endl;
1219 //std::cout << "EWSUSY(q) " << (DelA_d + DelA_u)/16.0/M_PI/M_PI << std::endl;
1220 //std::cout << "EWSM(had) " << getMyOneLoopEW()->DeltaAlpha_5q(Mz2) << std::endl;
1221
1222 return ( (DelA_l + DelA_d + DelA_u)/16.0/M_PI/M_PI );
1223}
1224
1225double EWSUSY::DeltaR_SUSY_EW1(const double Mw_i) const
1226{
1227 double cW2 = Mw_i*Mw_i/mySUSY.getMz()/mySUSY.getMz();
1228 double sW2 = 1.0 - cW2;
1229
1230 /* SM one-loop contributions */
1231 double DeltaAlphaL5q_EW1 = DeltaAlphaL5q_SM_EW1();
1232 double DeltaRho_EW1 = mySUSY.getMyOneLoopEW()->DeltaRho(Mw_i);
1233 double DeltaR_rem_EW1 = mySUSY.getMyOneLoopEW()->DeltaR_rem(Mw_i);
1234 double DeltaR_SM_EW1 = DeltaAlphaL5q_EW1 - cW2/sW2*DeltaRho_EW1 + DeltaR_rem_EW1;
1235
1236 /* Debug */
1237 //std::cout << std::endl;
1238 //std::cout << "DeltaAlphaL5q_EW1 = " << DeltaAlphaL5q_EW1 << std::endl;
1239 //std::cout << "-cW2/sW2*DeltaRho_EW1 = " << - cW2/sW2*DeltaRho_EW1 << std::endl;
1240 //std::cout << "DeltaR_rem_EW1 = " << DeltaR_rem_EW1 << std::endl;
1241 //std::cout << "DeltaR_SM_EW1 = " << DeltaR_SM_EW1 << std::endl;
1242
1243 return ( DeltaR_TOTAL_EW1(Mw_i) - DeltaR_SM_EW1 );
1244}
1245
1246double EWSUSY::Mw_MSSM_TMP(const double Mw_i) const
1247{
1248 if (mySUSY.getFlagMw().compare("NORESUM") != 0)
1249 throw std::runtime_error("EWSUSY::Mw_SUSY(): Scheme for Mw is not applicable");
1250
1251 double cW2 = Mw_i*Mw_i/mySUSY.getMz()/mySUSY.getMz();
1252 double sW2 = 1.0 - cW2;
1253 if (sW2 < 0.0)
1254 throw std::runtime_error("EWSUSY::Mw_SUSY(): negative sW2");
1255
1256 /* SM contributions to Delta r */
1257 double dAleL5q = mySUSY.DeltaAlphaL5q();
1258 double DeltaRho[StandardModel::orders_EW_size], DeltaRho_sum = 0.0;
1259 double DeltaR_rem[StandardModel::orders_EW_size], DeltaR_rem_sum = 0.0;
1260 mySUSY.ComputeDeltaRho(Mw_i, DeltaRho);
1261 mySUSY.ComputeDeltaR_rem(Mw_i, DeltaR_rem);
1262 for (int j=0; j<StandardModel::orders_EW_size; ++j) {
1263 /* excluding EW two-loop contributions, which will be added below */
1264 if (j!=(int)StandardModel::EW2) {
1265 DeltaRho_sum += DeltaRho[(StandardModel::orders_EW)j];
1266 DeltaR_rem_sum += DeltaR_rem[(StandardModel::orders_EW)j];
1267 }
1268 }
1269
1270 /* Full EW one-loop contribution (without the full DeltaAlphaL5q) */
1271 double DeltaR_EW1 = - cW2/sW2*DeltaRho[StandardModel::EW1] + DeltaR_rem[StandardModel::EW1];
1272
1273 /* Full EW two-loop contribution with reducible corrections */
1274 double DeltaR_EW2 = dAleL5q*dAleL5q + 2.0*dAleL5q*DeltaR_EW1
1275 + mySUSY.getMyApproximateFormulae()->DeltaR_TwoLoopEW_rem(Mw_i);
1276
1277 /* R = 1 + Delta r */
1278 double R = 1.0 + dAleL5q - cW2/sW2*DeltaRho_sum + DeltaR_rem_sum + DeltaR_EW2;
1279
1280 /* SUSY contribution */
1281 R += DeltaR_SUSY_EW1(Mw_i);
1282
1283 /* the W-boson mass in the complex pole scheme */
1284 double tmp = 4.0*M_PI*mySUSY.getAle()/sqrt(2.0)/mySUSY.getGF()/mySUSY.Mzbar()/mySUSY.Mzbar();
1285 if (tmp*R > 1.0) throw std::runtime_error("EWSUSY::Mw(): Negative (1-tmp*R)");
1286 double Mwbar = mySUSY.Mzbar()/sqrt(2.0) * sqrt(1.0 + sqrt(1.0 - tmp*R));
1287
1288 /* complex-pole/fixed-width scheme --> experimental/running-width scheme */
1289 double Mw_exp = mySUSY.MwFromMwbar(Mwbar);
1290
1291 if (Mw_exp >= mySUSY.getMz()) {
1292 std::cout << "WARNING: Mw > Mz in EWSUSY::Mw_MSSM_TMP" << std::endl;
1293 //std::cout << Mw_i << std::endl;
1294 //std::cout << Mw_exp << std::endl;
1295 //std::cout << - PiThat_W_0(Mw_i)/Mw_i/Mw_i << std::endl;
1296 //double mu = Mw_i, Mz = mySUSY.getMz();
1297 //std::cout << " " << PiT_W(mu, 0.0, Mw_i).real() << std::endl;
1298 //std::cout << " " << PiT_W(mu, Mw_i*Mw_i, Mw_i).real() << std::endl;
1299 //std::cout << " " << PiT_W(mu, Mw_i*Mw_i, Mw_i).real() << std::endl;
1300 //std::cout << " " << PiT_Z(mu, Mz*Mz, Mw_i).real() << std::endl;
1301 //std::cout << " " << PiT_AZ(mu, 0.0, Mw_i).real() << std::endl;
1302 //std::cout << " " << PiT_AZ(mu, 0.0, Mw_i).real() << std::endl;
1303 //std::cout << " " << PiTp_A(mu, 0.0, Mw_i).real() << std::endl;
1304 //std::cout << DeltaR_rem_SM(Mw_i) << std::endl;
1305 //std::cout << DeltaR_boxLL_SUSY(Mw_i) << std::endl;
1306 //std::cout << DeltaR_boxLR_SUSY(Mw_i) << std::endl;
1307 //std::cout << DeltaR_vertex_SUSY(Mw_i) << std::endl;
1308 //std::cout << DeltaR_neutrino_SUSY(Mw_i) << std::endl;
1309
1310 return Mw_unphysical;
1311 } else
1312 return Mw_exp;
1313}
1314
1315double EWSUSY::Mw_MSSM() const
1316{
1317 /* initial value for Mw */
1318 double Mw_org;
1319 Mw_org = mySUSY.Mw_tree();
1320 /* BEGIN: REMOVE FROM THE PACKAGE */
1321#if FEYNHIGGS
1322 if(mySUSY.IsFlag_FH()) Mw_org = mySUSY.getMyFH()->getMw_FHinput();
1323#endif
1324 /* END: REMOVE FROM THE PACKAGE */
1325
1326 double Mw = Mw_MSSM_TMP(Mw_org);
1327 //std::cout << std::endl << std::setprecision(12)
1328 // << "EWSUSY::Mw_MSSM(): Mw_org = " << Mw_org
1329 // << " Mw_new = " << Mw << std::endl;
1330
1331 if (Mw == Mw_unphysical) return Mw_unphysical;
1332
1333 /* iterations */
1334 while (fabs(Mw - Mw_org) > StandardModel::Mw_error) {
1335 Mw_org = Mw;
1336 Mw = Mw_MSSM_TMP(Mw);
1337 //std::cout << std::setprecision(12)
1338 // << "EWSUSY::Mw_MSSM(): Mw_org = " << Mw_org
1339 // << " Mw_new = " << Mw << std::endl;
1340
1341 if (Mw == Mw_unphysical) return Mw_unphysical;
1342 }
1343
1344 return Mw;
1345}
1346
1347gslpp::complex EWSUSY::L_esnC(const int N, const int K, const int i, const double Mw_i) const
1348{
1349 double e = sqrt(4.0*M_PI*mySUSY.getAle());
1350 double cW = Mw_i/mySUSY.getMz();
1351 double sW = sqrt(1.0 - cW*cW);
1352
1353 return ( e/sW*Zp(0,i)*Zne(N,K).conjugate() );
1354}
1355
1356gslpp::complex EWSUSY::R_esnC(const int N, const int K, const int i) const
1357{
1358 return ( Yl(N,N)*Zne(N,K).conjugate()*Zm(1,i).conjugate() );
1359}
1360
1361gslpp::complex EWSUSY::L_nLC(const int I, const int k, const int i, const double Mw_i) const
1362{
1363 double e = sqrt(4.0*M_PI*mySUSY.getAle());
1364 double cW = Mw_i/mySUSY.getMz();
1365 double sW = sqrt(1.0 - cW*cW);
1366
1367 return ( e/sW*ZL(I,k)*Zm(0,i) + Yl(I,I)*ZL(I+3,k)*Zm(1,i) );
1368}
1369
1370gslpp::complex EWSUSY::L_nsnN(const int J, const int K, const int j, const double Mw_i) const
1371{
1372 double e = sqrt(4.0*M_PI*mySUSY.getAle());
1373 double cW = Mw_i/mySUSY.getMz();
1374 double sW = sqrt(1.0 - cW*cW);
1375
1376 return ( - e/sqrt(2.0)/sW/cW*Zne(J,K).conjugate()*(ZN(0,j)*sW - ZN(1,j)*cW) );
1377}
1378
1379gslpp::complex EWSUSY::L_eLN(const int N, const int k, const int j, const double Mw_i) const
1380{
1381 double e = sqrt(4.0*M_PI*mySUSY.getAle());
1382 double cW = Mw_i/mySUSY.getMz();
1383 double sW = sqrt(1.0 - cW*cW);
1384
1385 return ( - e/sqrt(2.0)/sW/cW*ZL(N,k)*(ZN(0,j)*sW + ZN(1,j)*cW)
1386 - Yl(N,N)*ZL(N+3,k)*ZN(2,j) );
1387}
1388
1389gslpp::complex EWSUSY::R_eLN(const int N, const int k, const int j, const double Mw_i) const
1390{
1391 double e = sqrt(4.0*M_PI*mySUSY.getAle());
1392 double cW = Mw_i/mySUSY.getMz();
1393
1394 return ( e*sqrt(2.0)/cW*ZL(N+3,k)*ZN(0,j).conjugate()
1395 - Yl(N,N)*ZL(N,k)*ZN(2,j).conjugate() );
1396}
1397
1398gslpp::complex EWSUSY::F(const double m1, const double m2, const double m3,
1399 const double m4) const
1400{
1401 return PV.D0(0.0, 0.0, m1*m1, m2*m2, m3*m3, m4*m4);
1402}
1403
1404gslpp::complex EWSUSY::H(const double m1, const double m2, const double m3,
1405 const double m4) const
1406{
1407 return PV.D00(0.0, 0.0, m1*m1, m2*m2, m3*m3, m4*m4);
1408}
1409
1410gslpp::complex EWSUSY::f(const double m1, const double m2, const double m3) const
1411{
1412 return ( - PV.C0(0.0, m1*m1, m2*m2, m3*m3) );
1413}
1414
1415
EWSMOneLoopEW_HV.h
FULLNNLO
@ FULLNNLO
Definition: OrderScheme.h:39
EWSUSY::PiTp_A
gslpp::complex PiTp_A(const double mu, const double p2, const double Mw_i) const
The derivative of the transverse part of the photon self-energy with respect to , ,...
Definition: EWSUSY.cpp:647
EWSUSY::v
gslpp::complex v(const double mu, const QCD::lepton M, const QCD::lepton J, const double Mw_i) const
Definition: EWSUSY.cpp:962
EWSUSY::DeltaR_TOTAL_EW1
double DeltaR_TOTAL_EW1(const double Mw_i) const
The total one-loop contribution to in the MSSM.
Definition: EWSUSY.cpp:1151
EWSUSY::EWSUSY
EWSUSY(const SUSY &SUSY_in)
Constructor.
Definition: EWSUSY.cpp:33
EWSUSY::PiT_Z
gslpp::complex PiT_Z(const double mu, const double p2, const double Mw_i) const
The transverse part of the Z-boson self-energy, , in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:139
EWSUSY::DeltaR_rem_SM
double DeltaR_rem_SM(const double Mw_i) const
The SM one-loop renormalized vertex and box corrections to in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:768
EWSUSY::DeltaR_neutrino_SUSY
double DeltaR_neutrino_SUSY(const double Mw_i) const
The renormalized SUSY neutrino wave-function contribution to in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:1140
EWSUSY::DeltaAlphaL5q_SM_EW1
double DeltaAlphaL5q_SM_EW1() const
The SM one-loop leptonic and five-flavour-hadronic corrections to at Z-mass scale.
Definition: EWSUSY.cpp:1182
EWSUSY::Sigma_nu_0
gslpp::complex Sigma_nu_0(const double mu, const QCD::lepton I, const QCD::lepton J, const double Mw_i) const
The SUSY neutrino self-energy at zero momentum transfer in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:1092
EWSUSY::PiThat_W_0
double PiThat_W_0(const double Mw_i) const
The renormalized transverse W-boson self-energy at zero momentum transefer in the 't Hooft-Feynman ga...
Definition: EWSUSY.cpp:730
EWSUSY::dFA
gslpp::complex dFA(const double mu, const double p2, const double mi, const double mj, const gslpp::complex cV_aij, const gslpp::complex cV_bji, const gslpp::complex cA_aij, const gslpp::complex cA_bji) const
The derivative of with respect to .
Definition: EWSUSY.cpp:119
EWSUSY::DeltaR_vertex_SUSY
double DeltaR_vertex_SUSY(const double Mw_i) const
The renormalized SUSY vertex corrections to in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:1081
EWSUSY::DeltaR_boxLL_SUSY
double DeltaR_boxLL_SUSY(const double Mw_i) const
The LL SUSY box corrections to in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:778
EWSUSY::FA
gslpp::complex FA(const double mu, const double p2, const double mi, const double mj, const gslpp::complex cV_aij, const gslpp::complex cV_bji, const gslpp::complex cA_aij, const gslpp::complex cA_bji) const
Fermionic contribuiton to the transverse part of a gauge-boson self-energy, .
Definition: EWSUSY.cpp:101
EWSUSY::PiT_W
gslpp::complex PiT_W(const double mu, const double p2, const double Mw_i) const
The transverse part of the W-boson self-energy, , in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:337
EWSUSY::delta_v
gslpp::complex delta_v(const double mu, const QCD::lepton M, const QCD::lepton J, const double Mw_i) const
Definition: EWSUSY.cpp:1033
EWSUSY::SetRosiekParameters
void SetRosiekParameters()
Sets parameters in Rosiek's notation.
Definition: EWSUSY.cpp:43
EWSUSY::PiT_AZ
gslpp::complex PiT_AZ(const double mu, const double p2, const double Mw_i) const
The transverse part of the self-energy, , for the mixing between photon and Z boson in the 't Hooft-F...
Definition: EWSUSY.cpp:503
EWSUSY::DeltaR_SUSY_EW1
double DeltaR_SUSY_EW1(const double Mw_i) const
The one-loop SUSY contribution to .
Definition: EWSUSY.cpp:1225
EWSUSY::DeltaR_boxLR_SUSY
double DeltaR_boxLR_SUSY(const double Mw_i) const
The LR SUSY box corrections to in the 't Hooft-Feynman gauge.
Definition: EWSUSY.cpp:869
Mw
An observable class for the -boson mass.
Definition: Mw.h:22
QCD::quark
quark
An enum type for quarks.
Definition: QCD.h:323
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
QCD::NEUTRINO_2
@ NEUTRINO_2
Definition: QCD.h:313
QCD::NEUTRINO_1
@ NEUTRINO_1
Definition: QCD.h:311
QCD::MU
@ MU
Definition: QCD.h:314
QCD::ELECTRON
@ ELECTRON
Definition: QCD.h:312
QCD::NEUTRINO_3
@ NEUTRINO_3
Definition: QCD.h:315
QCD::TAU
@ TAU
Definition: QCD.h:316
SUSY
A base class for SUSY models.
Definition: SUSY.h:33
StandardModel::Mw_error
static const double Mw_error
The target accuracy of the iterative calculation of the -boson mass in units of GeV.
Definition: StandardModel.h:552
StandardModel::orders_EW
orders_EW
An enumerated type representing perturbative orders of radiative corrections to EW precision observab...
Definition: StandardModel.h:527
StandardModel::EW1
@ EW1
One-loop of .
Definition: StandardModel.h:528
StandardModel::EW2
@ EW2
Two-loop of .
Definition: StandardModel.h:531
StandardModel::orders_EW_size
@ orders_EW_size
The size of this enum.
Definition: StandardModel.h:534