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)