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