gambit is hosted by Hepforge, IPPP Durham
 GAMBIT  v1.5.0-2191-ga4742ac a Global And Modular Bsm Inference Tool
flav_loop_functions.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
17
18 #ifndef __flav_loop_functions_hpp
19 #define __flav_loop_functions_hpp
20
21 namespace Gambit
22 {
23
24  namespace FlavBit
25  {
26  using namespace std;
27
28  // Loop functions for LFV diagrams
29  namespace LoopFunctions
30  {
31  double G1(const double x)
32  {
33  if(x == 0)
34  return -7./12.;
35  if(x == 1)
36  return -5./12.;
37  else
38  return (-7. + 33.*x - 57.*x*x + 31.*x*x*x + 6.*x*x*(1. - 3*x)*log(x))/(12.*pow(1.-x,4));
39  }
40
41  double G1(const double a, const double b, const double c)
42  {
43  if(b == c and b != 0)
44  return G1(a/b)/b;
45  else
46  return 0; // TODO: 2C12 + 2C22 - C1 or 2C12 + C11 - C2
47  }
48
49  double MFVV(const double a, const double b)
50  {
51  if(a == b)
52  return 1. / (3. * b);
53  else if(a == 0)
54  return 5. / (9. * b);
55  else
56  return (6.*a*a*(a-3.*b)*log(a/b) - (a-b)*(5.*a*a - 22.*a*b + 5.*b*b))/(9.*pow(a-b,4));
57  }
58
59  double B1(const double a, const double b, const double Q)
60  {
61  if(a == b)
62  return 0.5 * log(b / pow(Q,2));
63  else if(a == 0)
64  return -0.25 + 0.5*log(b / pow(Q,2));
65  else
66  return -0.5 + 0.5*log(b / pow(Q,2)) - (a*a - b*b + 2.*a*a*log(b/a)) / (4.*pow(a-b,2));
67  }
68
69  double B0(const double a, const double b, const double Q)
70  {
71  // TODO: behaviour when a = 0 and b = 0 undefined
72  if(a == 0 and b == 0)
73  return 0;
74  else if(a == b)
75  return -log(b / pow(Q,2));
76  else if(a == 0)
77  return 1. - log(b / pow(Q,2));
78  else if(b == 0)
79  return 1. - log(a) - log(1./pow(Q,2));
80  else
81  return 1. - log(b / pow(Q,2)) + 1./(b-a) * a * log(a/b);
82  }
83
84  double C0(const double a, const double b, const double c)
85  {
86  // TODO: behaviour when two paramers are 0 undefined, set it to zero
87  if(a == 0 and b == 0 and c == 0)
88  return 0;
89  if(a == 0 and b == 0)
90  return 0;
91  else if(a == 0 and c == 0)
92  return 0;
93  else if(b == 0 and c == 0)
94  return 0;
95  else if(c == 0)
96  return C0(a,c,b);
97  else if(a == b and b == c)
98  return - 1./(2*c);
99  else if(a == b)
100  return (-b + c- c*log(c/b)) / pow(b-c,2);
101  else if(a == c and b != 0)
102  return C0(a,c,b);
103  else if(a == c and b == 0)
104  return -1./c;
105  else if(b == c and a != 0)
106  return (a - c + a*log(c/a)) / pow(a-c,2);
107  else if(b == c and a == 0)
108  return -1./c;
109  else if(a == 0)
110  return (-log(b) + log(c)) / (b-c);
111  else if(b == 0)
112  return log(c/a)/(a-c);
113  else
114  return -1. / (a-b)*(a-c)*(b-c)*( b*(c-a)*log(b/a) + c*(a-b)*log(c/a));
115  }
116
117  double C00(const double a, const double b, const double c, const double Q)
118  {
119  // TODO: behaviour when all three parameters are zero is undefined, set it to zero
120  if(a == 0 and b == 0 and c == 0)
121  return 0;
122  else if(b == 0 and c == 0)
123  return 0.125*(3. - 2.*log(a/pow(Q,2)));
124  else if(a == 0 and b == 0)
125  return 0.125*(3. - 2.*log(c) - 2.*log(1./pow(Q,2)));
126  else if(c == 0)
127  return C00(a,c,b,Q);
128  else if(a == b and b == c)
129  return -0.25*log(c/pow(Q,2));
130  else if(a == b)
131  return - (2.*c*c*log(c/b) + (b-c)*(-b + 3.*c + 2.*(b-c)*log(b/pow(Q,2))))/(8.*pow(b-c,2));
132  else if(a == c and b != 0)
133  return C00(a,c,b,Q);
134  else if(a == c and b == 0)
135  return 0.125*(1. - 2.*log(c/pow(Q,2)));
136  else if(b == c and a != 0)
137  return (2.*(2.*a-c)*c*log(c/a)-(a-c)*(-3.*a+c+2.*(a-c)*log(a/pow(Q,2))))/(8.*pow(a-c,2));
138  else if(b == c and a == 0)
139  return 0.125*(1. - 2.*log(c) - 2.*log(1./pow(Q,2)));
140  else if(a == 0)
141  return -(2.*b*log(b) - 2.*c*log(c) + (b-c)*(-3.+2.*log(1./pow(Q,2))))/(8.*(b-c));
142  else if(b == 0)
143  return (2.*c*log(c/a) - (a-c)*(-3. + 2.*log(a/pow(Q,2))))/(8.*(a-c));
144  else
145  return 1. / (8.*(a-b)*(a-c)*(b-c)) * ( (c-a)*((a-b)*(2.*log(a/pow(Q,2))-3.)*(b-c) - 2.*b*b*log(b/a)) + 2.*c*c*(b-a)*log(c/a));
146  }
147
148  // Finite combination of loop functions that appears in VZw10
149  double B02C00C0(const double a, const double b, const double c, const double Q)
150  {
151  if(a == 0 and b == 0)
152  return 0.25*(1.0 - 2.0*log(c) - 2.0*log(1 / pow(Q,2)));
153  else
154  return B0(a,b,Q) - 2*C00(a,b,c,Q) + C0(a,b,c)*c;
155  }
156
157  double D0(const double a, const double b, const double c, const double d)
158  {
159  //TODO: behaviour when two or more parameters are zero is undefined, set it to zero
160  if((!a and !b) or (!b and !c) or (!b and !d) or (!c and !d))
161  return 0;
162  else if(c == 0)
163  return D0(a,c,b,d);
164  else if(d == 0)
165  return D0(a,d,c,b);
166  else if(a == b and b == c and c == d)
167  return 1. / (6.*d*d);
168  else if(a == b and b == c)
169  return D0(a,d,c,d);
170  else if(a == b and b == d)
171  return D0(a,c,b,d);
172  else if(a == c and c == d and b == 0)
173  return 1. / (2.*c*c);
174  else if(a == c and c == d and b != 0)
175  return (-b*b + c*c + 2.*b*c*log(b/c)) / (2.*c*pow(c-b,3));
176  else if(b == c and c == d and a == 0)
177  return 1. / (2.*d*d);
178  else if(b == c and c == d and a != 0)
179  return (a*a - d*d + 2.*a*d*log(d/a)) / (2.*d*pow(a-d,3));
180  else if(a == d and b == c)
181  return (-2*c + 2*d + (c+d)*log(c/d)) / pow(c-d,3);
182  else if(a == d and b == 0)
183  return (c - d -d*log(c/d)) / (pow(c-d,2)*d);
184  else if(a == d)
185  return 1./ ((b-d)*(d-c))-(b*log(b/d))/((b-c)*pow(b-d,2))+(c*log(c/d))/((b-c)*pow(c-d,2));
186  else if(a == c)
187  return D0(a,b,d,c);
188  else if(a == b)
189  return D0(a,d,c,b);
190  else if(b == c)
191  return D0(a,d,c,b);
192  else if(b == d)
193  return D0(a,c,b,d);
194  else if(c == d and b == 0)
195  return (a - d + d*log(d/a)) / (d*pow(a-d,2));
196  else if(c == d and a == 0)
197  return (b - d + d*log(d/b)) / (d*pow(b-d,2));
198  else if(c == d)
199  return (b*pow(a-d,2)*log(b/a) - (a-b)*( (a-d)*(b-d) + (a*b-d*d)*log(d/a) )) / ((a-b)*pow(a-d,2)*pow(b-d,2));
200  else if(b == 0)
201  return log(c/a)/((a-c)*(c-d)) + log(d/a)/((a-d)*(d-c));
202  else if(a == 0)
203  return ((d-c)*log(b) + (b-d)*log(c) + (c-b)*log(d))/((b-c)*(b-d)*(c-d));
204  else
205  return -(b*log(b/a)/((b-a)*(b-c)*(b-d)) + c*log(c/a)/((c-a)*(c-b)*(c-d)) + d*log(d/a)/((d-a)*(d-b)*(d-c)));
206  }
207
208  double D27(const double a, const double b, const double c, const double d)
209  {
210  //TODO: behaviour when three or more parameters are zero is undefined, set it to zero
211  if((!a and !b and !c) or (!a and !b and !d) or (!a and !c and !d) or (!b and !c and !d))
212  return 0;
213  if(a == b and b == c and c == d)
214  return -1./(12.*d);
215  if(a == d and c == d and b == 0)
216  return -1. / (8.*d);
217  if(a == d and c == d)
218  return (3.*b*b - 4.*b*d + d*d - 2.*b*b*log(b/d))/(8.*pow(b-d,3));
219  if(b == c and c == d)
220  return D27(b,a,c,d);
221  if(a == b and b == c)
222  return D27(a,d,c,b);
223  if(a == b and b == d)
224  return D27(a,c,b,d);
225  if(a == b and c == d and b == 0)
226  return -1./(4.*d);
227  if(a == b and c == d and d == 0)
228  return -1./(4.*b);
229  if(a == b and c == d)
230  return (-b*b + d*d -2.*b*d*log(d/b)) / (4.*pow(b-d,3));
231  if(a == c and b == d)
232  return D27(a,c,b,d);
233  if(a == d and b == c)
234  return D27(a,b,d,c);
235  if(a == b and b == 0)
236  return log(d/c)/(4.*(c-d));
237  if(a == b and c == 0)
238  return - (b -d +d*log(d/b))/(4.*pow(b-d,2));
239  if(a == b and d == 0)
240  return D27(a,b,d,c);
241  if(a == b)
242  return 0.25*(-c*c*log(c/b)/(pow(b-c,2)*(c-d)) + (b*(d-b)/(b-c) + d*d*log(d/b)/(c-d))/pow(b-d,2));
243  if(a == c)
244  return D27(a,c,b,d);
245  if(a == d)
246  return D27(a,d,c,b);
247  if(b == c and a == 0)
248  return (-c+d+d*log(c/d))/(4.*pow(c-d,2));
249  if(b == c and c == 0)
250  return log(d/a)/(4.*(a-d));
251  if(b == c and d == 0)
252  return (a-c+a*log(c/a))/(4.*pow(a-c,2));
253  if(b == c)
254  return (c*(a-d)*(a*(c-2.*d)+c*d)*log(c/a)+(a-c)*(c*(a-d)*(c-d)+(a-c)*d*d*log(d/a)))/(4.*pow(a-c,2)*(a-d)*pow(c-d,2));
255  if(b == d)
256  return D27(a,b,d,c);
257  if(c == d)
258  return D27(a,c,d,b);
259  if(a == 0)
260  return (b*(-c+d)*log(b)+c*(b-d)*log(c)+(-b+c)*d*log(d))/(4.*(b-c)*(b-d)*(c-d));
261  if(b == 0)
262  return ((c*log(c/a))/((a - c)*(c - d)) + (d*log(d/a))/((a - d)*(-c + d)))/4.;
263  if(c == 0)
264  return D27(a,c,b,d);
265  if(d == 0)
266  return D27(a,d,c,b);
267  else
268  return -0.25*(b*b*log(b/a)/((b-a)*(b-c)*(b-d)) + c*c*log(c/a)/((c-a)*(c-b)*(c-d)) + d*d*log(d/a)/((d-a)*(d-b)*(d-c)));
269  }
270
271  double IC0D0(const double a, const double b, const double c, const double d)
272  {
273  return C0(a,b,c) + d*D0(a,b,c,d);
274  }
275  }
276
277  // Loop function for RK
278  namespace LoopFunctions
279  {
280  double E(const double x, const double y)
281  {
282  if(x == 0 or y == 0)
283  return 0.0;
284  if(x == y)
285  return (x*(-4.0 + 15.0*x - 12.0*pow(x,2) + pow(x,3) + 6.0*pow(x,2)*log(x)))/ (4.*pow(-1.0 + x,3));
286  return x*y*(-3.0/(4.0*(1.0-x)*(1.0 - y)) + ((0.25 - 3.0/(4.0*pow(-1.0 + x,2)) - 3.0/(2.0*(-1.0 + x)))*log(x))/(x - y) + ((0.25 - 3.0/(4.0*pow(-1.0 + y,2)) - 3.0/(2.0*(-1 + y)))*log(y))/(-x + y));
287  }
288
289  }
290
291  // Vertices for LFV diagrams
292  namespace Vertices
293  {
294  // Fermion-vector vertices
295  complex<double> VpL(int i, int j, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U)
296  {
297  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
298  return -1. / sqrt(2) * g2 * U(i,j);
299
300  }
301
302  double EL(int i,int j, SMInputs sminputs)
303  {
304  if(i != j) return 0;
305
306  double e = sqrt(4. * pi / sminputs.alphainv);
307  double g1 = e * sminputs.mZ / sminputs.mW;
308  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
309  double cw = sminputs.mW / sminputs.mZ;
310  double sw = sqrt(1. - cw*cw);
311  return 0.5 * (-g1*sw + g2*cw);
312
313  }
314
315  double ER(int i, int j, SMInputs sminputs)
316  {
317  if(i != j) return 0;
318
319  double e = sqrt(4. * pi / sminputs.alphainv);
320  double g1 = e * sminputs.mZ / sminputs.mW;
321  double cw = sminputs.mW / sminputs.mZ;
322  double sw = sqrt(1. - cw*cw);
323  return - g1*sw;
324  }
325
326  complex<double> VL(int i, int j, SMInputs sminputs)
327  {
328  double e = sqrt(4. * pi / sminputs.alphainv);
329  double g1 = e * sminputs.mZ / sminputs.mW;
330  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
331  double cw = sminputs.mW / sminputs.mZ;
332  double sw = sqrt(1. - cw*cw);
333
334  if(i == j)
335  return -0.5*(g1*sw + g2*cw);
336  else
337  return 0.;
338  }
339
340  complex<double> VR(int i, int j, SMInputs sminputs)
341  {
342  double e = sqrt(4. * pi / sminputs.alphainv);
343  double g1 = e * sminputs.mZ / sminputs.mW;
344  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
345  double cw = sminputs.mW / sminputs.mZ;
346  double sw = sqrt(1. - cw*cw);
347
348  if(i == j)
349  return 0.5*(g1*sw + g2*cw);
350  else
351  return 0.;
352  }
353
354  complex<double> DL(int i, int j, SMInputs sminputs)
355  {
356  double e = sqrt(4. * pi / sminputs.alphainv);
357  double g1 = e * sminputs.mZ / sminputs.mW;
358  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
359  double cw = sminputs.mW / sminputs.mZ;
360  double sw = sqrt(1. - cw*cw);
361
362  if(i == j)
363  return 1./6. * (3.*g2*cw + g1*sw);
364  else
365  return 0;
366  }
367
368  complex<double> DR(int i, int j, SMInputs sminputs)
369  {
370  double e = sqrt(4. * pi / sminputs.alphainv);
371  double g1 = e * sminputs.mZ / sminputs.mW;
372  double cw = sminputs.mW / sminputs.mZ;
373  double sw = sqrt(1. - cw*cw);
374
375  if(i == j)
376  return -1./3.*g1*sw;
377  else
378  return 0.;
379  }
380
381  complex<double> UL(int i, int j, SMInputs sminputs)
382  {
383  double e = sqrt(4. * pi / sminputs.alphainv);
384  double g1 = e * sminputs.mZ / sminputs.mW;
385  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
386  double cw = sminputs.mW / sminputs.mZ;
387  double sw = sqrt(1. - cw*cw);
388
389  if(i == j)
390  return -1./6. * (3.*g2*cw - g1*sw);
391  else
392  return 0;
393  }
394
395  complex<double> UR(int i, int j, SMInputs sminputs)
396  {
397  double e = sqrt(4. * pi / sminputs.alphainv);
398  double g1 = e * sminputs.mZ / sminputs.mW;
399  double cw = sminputs.mW / sminputs.mZ;
400  double sw = sqrt(1. - cw*cw);
401
402  if(i == j)
403  return 2./3.*g1*sw;
404  else
405  return 0.;
406  }
407
408  complex<double> VuL(int i, int j, SMInputs sminputs)
409  {
410  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
411  Eigen::Matrix3cd VCKM;
412  double lambda = sminputs.CKM.lambda, A = sminputs.CKM.A;
413  double rhobar = sminputs.CKM.rhobar, etabar = sminputs.CKM.etabar;
414  complex<double> I = {0,1};
415
416  complex<double> Vub = real(rhobar + I*etabar)*sqrt(1.-A*A*pow(lambda,4))/(sqrt(1.-pow(lambda,2))*(1.- A*A*pow(lambda,4)*(rhobar+I*etabar)));
417  double rho = real(Vub);
418  double eta = imag(Vub);
419
420  VCKM << 1. - 0.5*pow(lambda,2), lambda, A*pow(lambda,3)*(rho - I*eta),
421  -lambda, 1. - 0.5*pow(lambda,2), A*pow(lambda,2),
422  A*pow(lambda,3)*(1. - eta - I*eta), -A*pow(lambda,2), 1;
423
424  return -1./sqrt(2) * g2 * VCKM(i,j);
425  }
426
427  // Vector vertices
428  double Fw(SMInputs sminputs)
429  {
430  return sqrt(4.* pi/ sminputs.alphainv);
431  }
432
433  double Zww(SMInputs sminputs)
434  {
435  double g2 = sminputs.mW * sqrt( 8. * sminputs.GF / sqrt(2));
436  return -g2 * sminputs.mW / sminputs.mZ;
437  }
438
439  // Scalar vertices
440  double HL(int i, int j, SMInputs sminputs)
441  {
442  double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
443
444  if(i == 0 and j == 0)
445  return -1. / vev * sminputs.mE;
446  if(i == 1 and j == 1)
447  return -1. / vev * sminputs.mMu;
448  if(i == 2 and j == 2)
449  return -1. / vev * sminputs.mTau;
450  else
451  return 0;
452  }
453
454  double HR(int i, int j, SMInputs sminputs)
455  {
456  return HL(i, j , sminputs);
457  }
458
459  double HdL(int i, int j, SMInputs sminputs)
460  {
461  double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
462
463  if(i == 0 and j == 0)
464  return -1. / vev * sminputs.mD;
465  else if(i == 1 and j == 1)
466  return -1. / vev * sminputs.mS;
467  else if(i == 2 and j == 2)
468  return -1. / vev * sminputs.mBmB;
469  else
470  return 0;
471  }
472
473  double HdR(int i, int j, SMInputs sminputs)
474  {
475  return HdL(i, j, sminputs);
476  }
477
478  double HuL(int i, int j, SMInputs sminputs)
479  {
480  double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
481
482  if(i == 0 and j == 0)
483  return -1. / vev * sminputs.mU;
484  else if(i == 1 and j == 1)
485  return -1. / vev * sminputs.mCmC;
486  else if(i == 2 and j == 2)
487  return -1. / vev * sminputs.mT;
488  else
489  return 0;
490  }
491
492  double HuR(int i, int j, SMInputs sminputs)
493  {
494  return HuL(i, j, sminputs);
495  }
496
497
498  }
499
500  // Penguin contributions
501  namespace Penguins
502  {
503  // Fotonic penguins
504
505  complex<double> A1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
506  {
507  complex<double> a1r = {0,0};
508
509  for(int a=0; a<6; a++)
510  {
511  a1r += Vertices::Fw(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * LoopFunctions::MFVV(pow(mnu[a],2), pow(sminputs.mW,2));
512  }
513
514  return a1r;
515  }
516
517  complex<double> A2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
518  {
519  complex<double> a2l = {0,0};
520  double mW = sminputs.mW;
521
522  for(int a=0; a<6; a++)
523  a2l += -2. * Vertices::Fw(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * LoopFunctions::G1(pow(mnu[a],2), mW*mW, mW*mW) * ml[beta];
524
525  return a2l;
526  }
527
528  complex<double> A2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
529  {
530  complex<double> a2r = {0,0};
531  double mW = sminputs.mW;
532  for(int a=0; a<6; a++)
533  a2r += -2. * Vertices::Fw(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * LoopFunctions::G1(pow(mnu[a],2), mW*mW, mW*mW) * ml[alpha];
534
535  return a2r;
536  }
537
538  // Z penguins
539  complex<double> VZw2w4LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
540  {
541  complex<double> vzll = {0,0};
542
543  for(int a=0; a<6; a++)
544  for(int c=0; c<3; c++)
545  {
546  // Use MZ for the renormalization scale Q
547  if(beta == c)
548  vzll += Vertices::EL(beta,c, sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* pow(ml[alpha],2) / (pow(ml[alpha],2) - pow(ml[c],2));
549  if(alpha == c)
550  vzll += Vertices::EL(alpha,c, sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* pow(ml[beta],2) / (pow(ml[beta],2) - pow(ml[c],2));
551  }
552
553  return vzll;
554  }
555
556  complex<double> VZw2w4LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
557  {
558  return VZw2w4LL(alpha,beta,sminputs,U,ml,mnu);
559  }
560
561  complex<double> VZw2w4RR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
562  {
563  complex<double> vzrr = {0,0};
564
565  for(int a=0; a<6; a++)
566  for(int c=0; c<3; c++)
567  {
568  if(beta == c)
569  vzrr += Vertices::ER(beta,c,sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* ml[c]*ml[alpha] / (pow(ml[alpha],2) - pow(ml[c],2));
570  if(alpha == c)
571  vzrr += Vertices::ER(alpha,c, sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),pow(sminputs.mW,2),sminputs.mZ))* ml[c]*ml[beta] / (pow(ml[beta],2) - pow(ml[c],2));
572  }
573
574  return vzrr;
575  }
576
577  complex<double> VZw2w4RL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
578  {
579  return VZw2w4RR(alpha, beta, sminputs, U, ml, mnu);
580  }
581
582  complex<double> VZw8LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
583  {
584  complex<double> vzll = {0,0};
585  double mW = sminputs.mW;
586
587  // Use MZ as the renormalization scale Q
588  for(int a=0; a<6; a++)
589  vzll += Vertices::Zww(sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * (1. - 2.*(LoopFunctions::B0(mW*mW,mW*mW,sminputs.mZ) + 2.*LoopFunctions::C00(pow(mnu[a],2),mW*mW,mW*mW,sminputs.mZ) + LoopFunctions::C0(pow(mnu[a],2),mW*mW,mW*mW)*pow(mnu[a],2)));
590
591  return vzll;
592  }
593
594  complex<double> VZw8LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
595  {
596  return VZw8LL(alpha, beta, sminputs, U, mnu);
597  }
598
599  complex<double> VZw10LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
600  {
601  complex<double> vzll = {0,0};
602  double mW = sminputs.mW;
603
604  // Use MZ as the renormalization scale Q
605  for(int b=0; b<6; b++)
606  {
607  // Use different loop function in case that mnu[b] 0
608  if(mnu[b])
609  vzll += - Vertices::VpL(alpha,b,sminputs,U) * conj(Vertices::VpL(beta,b,sminputs,U)) * (2.* Vertices::VR(b,b,sminputs) * LoopFunctions::C0(pow(mnu[b],2),pow(mnu[b],2),mW*mW) * mnu[b] * mnu[b] + Vertices::VL(b,b,sminputs) * (1. - 2.*(LoopFunctions::B0(pow(mnu[b],2),pow(mnu[b],2),sminputs.mZ) - 2.*LoopFunctions::C00(pow(mnu[b],2),pow(mnu[b],2),mW*mW,sminputs.mZ) + LoopFunctions::C0(pow(mnu[b],2),pow(mnu[b],2),mW*mW)*mW*mW)));
610  else
611  vzll += - Vertices::VpL(alpha,b,sminputs,U) * conj(Vertices::VpL(beta,b,sminputs,U)) * Vertices::VL(b,b,sminputs) * (1. - 2.*(LoopFunctions::B02C00C0(pow(mnu[b],2),pow(mnu[b],2),mW*mW,sminputs.mZ)));
612  }
613
614  return vzll;
615  }
616
617  complex<double> VZw10LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
618  {
619  return VZw10LL(alpha, beta, sminputs, U, mnu);
620  }
621
622  // Sum over Z penguins
623  complex<double> VZsumLL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
624  {
625  return 1. / (16.*pow(pi,2)) * (VZw2w4LL(alpha, beta, sminputs, U, ml, mnu) + VZw8LL(alpha, beta, sminputs, U, mnu) + VZw10LL(alpha, beta, sminputs, U, mnu));
626  }
627
628  complex<double> VZsumLR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
629  {
630  return 1. / (16.*pow(pi,2)) * (VZw2w4LR(alpha, beta, sminputs, U, ml, mnu) + VZw8LR(alpha, beta, sminputs, U, mnu) + VZw10LR(alpha, beta, sminputs, U, mnu));
631  }
632
633  complex<double> VZsumRL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
634  {
635  return 1. / (16.*pow(pi,2)) * (VZw2w4RL(alpha, beta, sminputs, U, ml, mnu));
636  }
637
638  complex<double> VZsumRR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
639  {
640  return 1. / (16.*pow(pi,2)) * (VZw2w4RR(alpha, beta, sminputs, U, ml, mnu));
641  }
642
643  // Scalar penguins
644  complex<double> Shw2w4LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
645  {
646  complex<double> shll = {0,0};
647  double mW = sminputs.mW;
648
649  // Use mZ for the renormalisation scale Q
650  for(int a=0; a<6; a++)
651  for(int c=0; c<3; c++)
652  {
653  if(beta == c)
654  shll += - (Vertices::HL(beta,c,sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * pow(ml[alpha],2))/(pow(ml[alpha],2) - pow(ml[c],2));
655  if(alpha == c)
656  shll += - (Vertices::HL(alpha,c,sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * pow(ml[beta],2))/(pow(ml[beta],2) - pow(ml[c],2));
657  }
658
659  return shll;
660  }
661
662  complex<double> Shw2w4LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
663  {
664  return Shw2w4LL(alpha, beta, sminputs, U, ml, mnu);
665  }
666
667  complex<double> Shw2w4RR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
668  {
669  complex<double> shrr = {0,0};
670  double mW = sminputs.mW;
671
672  // Use mZ for the renormalisation scale Q
673  for(int a=0; a<6; a++)
674  for(int c=0; c<3; c++)
675  {
676  if(beta == c)
677  shrr += - (Vertices::HR(beta,c,sminputs) * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * ml[c]*ml[alpha])/(pow(ml[alpha],2) - pow(ml[c],2));
678  if(alpha == c)
679  shrr += - (Vertices::HR(alpha,c,sminputs) * Vertices::VpL(beta,a,sminputs,U) * conj(Vertices::VpL(c,a,sminputs,U)) * (1. + 2.* LoopFunctions::B1(pow(mnu[a],2),mW*mW, sminputs.mZ)) * ml[c]*ml[beta])/(pow(ml[beta],2) - pow(ml[c],2));
680  }
681
682  return shrr;
683  }
684
685  complex<double> Shw2w4RL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
686  {
687  return Shw2w4RR(alpha, beta, sminputs, U, ml, mnu);
688  }
689
690  // Sum over scalar penguins
691  complex<double> ShsumLL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
692  {
693  return 1. / (16.*pow(pi,2)) * Shw2w4LL(alpha, beta, sminputs, U, ml, mnu);
694  }
695
696  complex<double> ShsumLR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
697  {
698  return 1. / (16.*pow(pi,2)) * Shw2w4LR(alpha, beta, sminputs, U, ml, mnu);
699  }
700
701  complex<double> ShsumRL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
702  {
703  return 1. / (16.*pow(pi,2)) * Shw2w4RL(alpha, beta, sminputs, U, ml, mnu);
704  }
705
706  complex<double> ShsumRR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
707  {
708  return 1. / (16.*pow(pi,2)) * Shw2w4RR(alpha, beta, sminputs, U, ml, mnu);
709  }
710
711  }
712
713  // Box contributions
714  namespace Boxes
715  {
716  complex<double> Vw4lLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
717  {
718  complex<double> vll = {0,0};
719  double mW = sminputs.mW;
720
721  for(int a=0; a<6; a++)
722  for(int c=0; c<6; c++)
723  vll += -4. * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(beta,a,sminputs,U)) * Vertices::VpL(gamma,c,sminputs,U) * conj(Vertices::VpL(delta,c,sminputs,U)) * (LoopFunctions::IC0D0(pow(mnu[c],2),mW*mW, mW*mW, pow(mnu[a],2)) - 3. * LoopFunctions::D27(pow(mnu[a],2),pow(mnu[c],2),mW*mW,mW*mW));
724
725  return vll;
726  }
727
728  complex<double> Vw8lLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
729  {
730  complex<double> vll = {0,0};
731  double mW = sminputs.mW;
732
733  for(int a=0; a<6; a++)
734  for(int c=0; c<6; c++)
735  vll += -2. * Vertices::VpL(alpha,a,sminputs,U) * conj(Vertices::VpL(delta,c,sminputs,U)) * Vertices::VpL(gamma,a,sminputs,U) * conj(Vertices::VpL(beta,c,sminputs,U)) * mnu[a] * mnu[c] * LoopFunctions::D0(pow(mnu[a],2),pow(mnu[c],2),mW*mW,mW*mW);
736
737  return vll;
738  }
739
740  complex<double> Vw4lpLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
741  {
742  return Vw4lLL(alpha, delta, gamma, beta, sminputs, U, mnu);
743  }
744
745  complex<double> Vw8lpLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
746  {
747  return Vw8lLL(alpha, delta, gamma, beta, sminputs, U, mnu);
748  }
749
750  complex<double> Vw4dLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
751  {
752  complex<double> vll = {0,0};
753  double mW = sminputs.mW;
754  vector<double> mu = {sminputs.mU, sminputs.mCmC, sminputs.mT};
755
756  for(int a=0; a<6; a++)
757  for(int c=0; c<3; c++)
758  vll += -4.*Vertices::VpL(alpha,a,sminputs,U)*conj(Vertices::VpL(beta,a,sminputs,U))*Vertices::VuL(gamma,c,sminputs)*conj(Vertices::VuL(delta,c,sminputs))*(LoopFunctions::IC0D0(pow(mu[c],2),mW*mW, mW*mW, pow(mnu[a],2)) - 3.*LoopFunctions::D27(pow(mnu[a],2),pow(mu[c],2),mW*mW,mW*mW));
759
760  return vll;
761  }
762
763  complex<double> Vw4uLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
764  {
765  complex<double> vll = {0,0};
766  double mW = sminputs.mW;
767  vector<double> md = {sminputs.mD, sminputs.mS, sminputs.mBmB};
768
769  for(int a=0; a<6; a++)
770  for(int c=0; c<3; c++)
771  vll += 16.*Vertices::VpL(alpha,a,sminputs,U)*conj(Vertices::VpL(beta,a,sminputs,U))*Vertices::VuL(delta,c,sminputs)*conj(Vertices::VuL(gamma,c,sminputs))*LoopFunctions::D27(pow(mnu[a],2),pow(md[c],2),mW*mW,mW*mW);
772
773  return vll;
774  }
775
776  // Sum over boxes
777  complex<double> VsumlLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
778  {
779  return 1. / (16.*pow(pi,2)) *( Vw4lLL(alpha, beta, gamma, delta, sminputs, U, mnu) + Vw8lLL(alpha, beta, gamma, delta, sminputs, U, mnu) + Vw4lpLL(alpha, beta, gamma, delta, sminputs, U, mnu) + Vw8lpLL(alpha, beta, gamma, delta, sminputs, U, mnu));
780  }
781
782  complex<double> VsumdLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
783  {
784  return 1./(16.*pow(pi,2)) *Vw4dLL(alpha, beta, gamma, delta, sminputs, U, mnu);
785  }
786
787  complex<double> VsumuLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
788  {
789  return 1./(16.*pow(pi,2)) *Vw4uLL(alpha, beta, gamma, delta, sminputs, U, mnu);
790  }
791
792  } // Diagrams
793
794
795  // Form factors for LFV diagrams
796  namespace FormFactors
797  {
798
799  complex<double> K1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> mnu)
800  {
801  double e = sqrt(4. * pi / sminputs.alphainv);
802
803  return 1. / (16*pow(pi,2)*e) * Penguins::A1R(alpha, beta, sminputs, U, mnu);
804  }
805
806  complex<double> K2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
807  {
808  double e = sqrt(4. * pi / sminputs.alphainv);
809
810  return 1. / (2. * 16.*pow(pi,2) * e * ml[alpha] ) * Penguins::A2L(alpha, beta, sminputs, U, ml, mnu);
811  }
812
813  complex<double> K2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
814  {
815  double e = sqrt(4. * pi / sminputs.alphainv);
816
817  return 1. / (2. * 16.*pow(pi,2)* e * ml[alpha] ) * Penguins::A2R(alpha, beta, sminputs, U, ml, mnu);
818  }
819
820  complex<double> AVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
821  {
822  return Penguins::VZsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::EL(gamma,delta,sminputs) / pow(sminputs.mZ,2) + Boxes::VsumlLL(alpha,beta,gamma,delta,sminputs,U,mnu);
823  }
824
825  complex<double> AVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
826  {
827  return Penguins::VZsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::ER(gamma,delta,sminputs) / pow(sminputs.mZ,2);
828  }
829
830  complex<double> AVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
831  {
832  return Penguins::VZsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::EL(gamma,delta,sminputs) / pow(sminputs.mZ,2);
833  }
834
835  complex<double> AVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
836  {
837  return Penguins::VZsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::ER(gamma,delta,sminputs) / pow(sminputs.mZ,2);
838  }
839
840  complex<double> ASLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
841  {
842  return Penguins::ShsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HL(gamma,delta,sminputs) / pow(mh,2);
843  }
844
845  complex<double> ASLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
846  {
847  return Penguins::ShsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HR(gamma,delta,sminputs) / pow(mh,2);
848  }
849
850  complex<double> ASRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
851  {
852  return Penguins::ShsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HL(gamma,delta,sminputs) / pow(mh,2);
853  }
854
855  complex<double> ASRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
856  {
857  return Penguins::ShsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HR(gamma,delta,sminputs) / pow(mh,2);
858  }
859
860  complex<double> BVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
861  {
862  return Penguins::VZsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::DL(gamma,delta,sminputs) / pow(sminputs.mZ,2) + Boxes::VsumdLL(alpha,beta,gamma,delta,sminputs,U,mnu);
863  }
864
865  complex<double> BVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
866  {
867  return Penguins::VZsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::DR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
868  }
869
870  complex<double> BVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
871  {
872  return Penguins::VZsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::DL(gamma,delta,sminputs) / pow(sminputs.mZ,2);
873  }
874
875  complex<double> BVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
876  {
877  return Penguins::VZsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::DR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
878  }
879
880  complex<double> BSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
881  {
882  return Penguins::ShsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdL(gamma,delta,sminputs) / pow(mh,2);
883  }
884
885  complex<double> BSLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
886  {
887  return Penguins::ShsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdR(gamma,delta,sminputs) / pow(mh,2);
888  }
889
890  complex<double> BSRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
891  {
892  return Penguins::ShsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdL(gamma,delta,sminputs) / pow(mh,2);
893  }
894
895  complex<double> BSRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
896  {
897  return Penguins::ShsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HdR(gamma,delta,sminputs) / pow(mh,2);
898  }
899
900  complex<double> CVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
901  {
902  return Penguins::VZsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::UL(gamma,delta,sminputs) / pow(sminputs.mZ,2) + Boxes::VsumuLL(alpha,beta,gamma,delta,sminputs,U,mnu);
903  }
904
905  complex<double> CVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
906  {
907  return Penguins::VZsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::UR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
908  }
909
910  complex<double> CVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
911  {
912  return Penguins::VZsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::UL(gamma,delta,sminputs) / pow(sminputs.mZ,2);
913  }
914
915  complex<double> CVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu)
916  {
917  return Penguins::VZsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::UR(gamma,delta,sminputs) / pow(sminputs.mZ,2);
918  }
919
920  complex<double> CSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
921  {
922  return Penguins::ShsumLL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuL(gamma,delta,sminputs) / pow(mh,2);
923  }
924
925  complex<double> CSLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
926  {
927  return Penguins::ShsumLR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuR(gamma,delta,sminputs) / pow(mh,2);
928  }
929
930  complex<double> CSRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
931  {
932  return Penguins::ShsumRL(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuL(gamma,delta,sminputs) / pow(mh,2);
933  }
934
935  complex<double> CSRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix<complex<double>,3,6> U, vector<double> ml, vector<double> mnu, double mh)
936  {
937  return Penguins::ShsumRR(alpha,beta,sminputs,U,ml,mnu)*Vertices::HuR(gamma,delta,sminputs) / pow(mh,2);
938  }
939
940  } // Form Factors
941  }
942 }
943
944 #endif //#defined __flav_loop_functions_hpp__
complex< double > VL(int i, int j, SMInputs sminputs)
complex< double > Shw2w4RR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > VZw10LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > A2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > BVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > Vw4dLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > A1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > CSRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
double G1(const double a, const double b, const double c)
double EL(int i, int j, SMInputs sminputs)
complex< double > Shw2w4LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
START_MODEL beta
Definition: Axions.hpp:36
double HdR(int i, int j, SMInputs sminputs)
complex< double > AVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > Shw2w4LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > VZw2w4RR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > Vw8lpLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > ShsumRL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
Funk delta(std::string arg, double pos, double width)
Definition: daFunk.hpp:902
double HL(int i, int j, SMInputs sminputs)
complex< double > Vw4lpLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > VZsumLR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
double HdL(int i, int j, SMInputs sminputs)
complex< double > VR(int i, int j, SMInputs sminputs)
STL namespace.
complex< double > ShsumLL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > AVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > UL(int i, int j, SMInputs sminputs)
complex< double > VZw10LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > VuL(int i, int j, SMInputs sminputs)
double Fw(SMInputs sminputs)
double D0(const double a, const double b, const double c, const double d)
double B02C00C0(const double a, const double b, const double c, const double Q)
complex< double > VsumlLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > A2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
START_MODEL b
Definition: demo.hpp:270
const double mW
Definition: topness.h:40
complex< double > VZsumRR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > BSLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > K1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > VZw8LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
double ER(int i, int j, SMInputs sminputs)
complex< double > ASRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > CSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
double C0(const double a, const double b, const double c)
START_MODEL alpha
complex< double > Vw4lLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > VpL(int i, int j, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U)
complex< double > DR(int i, int j, SMInputs sminputs)
complex< double > VZw2w4LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > CSLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > CVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
double HuL(int i, int j, SMInputs sminputs)
double HuR(int i, int j, SMInputs sminputs)
double Zww(SMInputs sminputs)
const double pi
complex< double > VZw8LL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > VsumdLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
double HR(int i, int j, SMInputs sminputs)
const double mu
Definition: SM_Z.hpp:42
complex< double > Vw4uLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > BSRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > ASLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
double IC0D0(const double a, const double b, const double c, const double d)
complex< double > ASRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > BSRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
double B0(const double a, const double b, const double Q)
complex< double > ShsumRR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > K2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > AVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
double D27(const double a, const double b, const double c, const double d)
complex< double > DL(int i, int j, SMInputs sminputs)
double C00(const double a, const double b, const double c, const double Q)
complex< double > AVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > CVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > Vw8lLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
double MFVV(const double a, const double b)
complex< double > VZw2w4RL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > UR(int i, int j, SMInputs sminputs)
complex< double > Shw2w4RL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > K2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
complex< double > VZw2w4LR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > BVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > BVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
double B1(const double a, const double b, const double Q)
complex< double > VZsumRL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > ASLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > VsumuLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
complex< double > BVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)
complex< double > CVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > VZsumLL(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > ShsumLR(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > CSRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
Container class for Standard Model input information (defined as in SLHA2)
Definition: sminputs.hpp:29
complex< double > CVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > BSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)