Belle II Software  release-06-01-15
ME-pair.cc
1 #include <math.h>
2 #include <cstdlib>
3 #include <cstdio>
4 #include <iostream>
5 #include <complex>
6 #include "ME-pair.hh"
7 #include "TString.h"
8 #include "TLorentzVector.h"
9 using namespace std;
10 
11 //----------------------------------------------------------------------
12 //
13 //
14 // PHOB: PHotosBoost
15 //
16 // Purpose: Boosts VEC to (MODE=1) rest frame of PBOOS1;
17 // or back (MODE=1)
18 //
19 // Input Parameters: MODE,PBOOS1,VEC
20 //
21 // Output Parameters: VEC
22 //
23 // Author(s): Created at: 08/12/05
24 // Z. Was Last Update: 13/06/13
25 //
26 //----------------------------------------------------------------------
27 
28 void PHOB(int MODE, double PBOOS1[4], double vec[4])
29 {
30  double BET1[3], GAM1, PB;
31  static int j0 = 1;
32  int J;
33 
34 
35  PB = sqrt(PBOOS1[4 - j0] * PBOOS1[4 - j0] - PBOOS1[3 - j0] * PBOOS1[3 - j0] - PBOOS1[2 - j0] * PBOOS1[2 - j0] - PBOOS1[1 - j0] *
36  PBOOS1[1 - j0]);
37  for (J = 1; J < 4; J++) {
38  if (MODE == 1) BET1[J - j0] = -PBOOS1[J - j0] / PB;
39  else BET1[J - j0] = PBOOS1[J - j0] / PB;
40  }
41 
42  GAM1 = PBOOS1[4 - j0] / PB;
43 
44  //--
45  //-- Boost vector
46 
47  PB = BET1[1 - j0] * vec[1 - j0] + BET1[2 - j0] * vec[2 - j0] + BET1[3 - j0] * vec[3 - j0];
48 
49  for (J = 1; J < 4; J++) vec[J - j0] = vec[J - j0] + BET1[J - j0] * (vec[4 - j0] + PB / (GAM1 + 1.0));
50  vec[4 - j0] = GAM1 * vec[4 - j0] + PB;
51  //--
52 }
53 
54 void coutLorentzVector(const TLorentzVector& p4, bool inMeV = false)
55 {
56  double factor = (inMeV) ? 1000. : 1;
57  TString units = (inMeV) ? "MeV" : "GeV";
58  // cout << Form("M,E,Rho,Et,Pt,Eta,Theta,Phi,X,Y,Z (%s) = ",units.Data())
59  cout << Form("M,E,Rho,Theta,Phi,X,Y,Z (%s) = ", units.Data())
60  << Form("% 10.5f", p4.M() / factor) << " "
61  << Form("% 10.5f", p4.E() / factor) << " "
62  << Form("% 10.5f", p4.Rho() / factor) << " "
63  << Form("% 10.5f", p4.Theta() * (180. / acos(-1.))) << " "
64  << Form("% 10.5f", p4.Phi() * (180. / acos(-1.))) << " "
65  << Form("% 10.5f", p4.Px() / factor) << " "
66  << Form("% 10.5f", p4.Py() / factor) << " "
67  << Form("% 10.5f", p4.Pz() / factor) << " "
68  << endl;
69 }
70 double dilog(double X)
71 {
72 // CERN C304 VERSION 29/07/71 DILOG 59 C
73  double Z, T, S, Y, B, A, di;
74  Z = -1.64493406684822;
75  if (X < -1.0) goto case1;
76  if (X <= 0.5) goto case2;
77  if (X == 1.0) goto case3;
78  if (X <= 2.0) goto case4;
79  Z = 3.2898681336964;
80 case1:
81  T = 1.0 / X;
82  S = -0.5;
83  Z = Z - 0.5 * log(abs(X)) * log(abs(X));
84  goto case5;
85 case2:
86  T = X;
87  S = 0.5;
88  Z = 0.0;
89  goto case5;
90 case3:
91  di = 1.64493406684822;
92  return di;
93 case4:
94  T = 1.0 - X;
95  S = -0.5;
96  Z = 1.64493406684822 - log(X) * log(abs(T));
97 case5:
98  Y = 2.66666666666666 * T + 0.66666666666666;
99  B = 0.000000000000001;
100  A = Y * B + 0.000000000000004;
101  B = Y * A - B + 0.000000000000011;
102  A = Y * B - A + 0.000000000000037;
103  B = Y * A - B + 0.000000000000121;
104  A = Y * B - A + 0.000000000000398;
105  B = Y * A - B + 0.000000000001312;
106  A = Y * B - A + 0.000000000004342;
107  B = Y * A - B + 0.000000000014437;
108  A = Y * B - A + 0.000000000048274;
109  B = Y * A - B + 0.000000000162421;
110  A = Y * B - A + 0.000000000550291;
111  B = Y * A - B + 0.000000001879117;
112  A = Y * B - A + 0.000000006474338;
113  B = Y * A - B + 0.000000022536705;
114  A = Y * B - A + 0.000000079387055;
115  B = Y * A - B + 0.000000283575385;
116  A = Y * B - A + 0.000001029904264;
117  B = Y * A - B + 0.000003816329463;
118  A = Y * B - A + 0.000014496300557;
119  B = Y * A - B + 0.000056817822718;
120  A = Y * B - A + 0.000232002196094;
121  B = Y * A - B + 0.001001627496164;
122  A = Y * B - A + 0.004686361959447;
123  B = Y * A - B + 0.024879322924228;
124  A = Y * B - A + 0.166073032927855;
125  A = Y * A - B + 1.93506430086996;
126  di = S * T * (A - B) + Z;
127  return di;
128 }
129 
130 //***********************************************************************************************************
131 // demo pi- l- l+ ME, based on paper arxiv.org/abs/1306.1732 (updated to also generate dark-photon model)
132 //***********************************************************************************************************
133 void eepi_(const double& m_dark, const double& dark_width, const double& Mlep, const double& Mpion,
134  const double* pt, const double* pn, const double* p1, const double* p2, const double* p3,
135  double& amplit, double* hv, int& iflag)
136 {
137  if (iflag == 0) {
138  (Mlep < 0.1) ? e_all++ : m_all++;
139  }
140 
141  double p23[4]; // momentum of l-l+ pair
142  p23[3] = p3[3] + p2[3];
143  p23[2] = p3[2] + p2[2];
144  p23[1] = p3[1] + p2[1];
145  p23[0] = p3[0] + p2[0];
146  double m2ee = p23[3] * p23[3] - p23[2] * p23[2] - p23[1] * p23[1] - p23[0] * p23[0];
147 
148  double col[4];
149  col[3] = p1[3] * p23[3];
150  col[2] = p1[2] * p23[2];
151  col[1] = p1[1] * p23[1];
152  col[0] = p1[0] * p23[0];
153  double pk = col[3] - col[2] - col[1] - col[0];
154 
155  double ptk[4];
156  ptk[3] = pt[3] * p23[3];
157  ptk[2] = pt[2] * p23[2];
158  ptk[1] = pt[1] * p23[1];
159  ptk[0] = pt[0] * p23[0];
160  double ptp23 = ptk[3] - ptk[2] - ptk[1] - ptk[0];
161 
162  double kpn[4];
163  kpn[3] = pn[3] * p23[3];
164  kpn[2] = pn[2] * p23[2];
165  kpn[1] = pn[1] * p23[1];
166  kpn[0] = pn[0] * p23[0];
167  double pnk = kpn[3] - kpn[2] - kpn[1] - kpn[0];
168 
169  double ptq[4];
170  ptq[3] = pn[3] * pt[3];
171  ptq[2] = pn[2] * pt[2];
172  ptq[1] = pn[1] * pt[1];
173  ptq[0] = pn[0] * pt[0];
174  double ptpn = ptq[3] - ptq[2] - ptq[1] - ptq[0];
175 
176  double pp[4];
177  pp[3] = p2[3] * p3[3];
178  pp[2] = p2[2] * p3[2];
179  pp[1] = p2[1] * p3[1];
180  pp[0] = p2[0] * p3[0];
181  double p2p3 = pp[3] - pp[2] - pp[1] - pp[0];
182 
183  double tmp[4];
184  tmp[3] = p2[3] * pt[3];
185  tmp[2] = p2[2] * pt[2];
186  tmp[1] = p2[1] * pt[1];
187  tmp[0] = p2[0] * pt[0];
188  double p2pt = tmp[3] - tmp[2] - tmp[1] - tmp[0];
189  tmp[3] = p3[3] * pt[3];
190  tmp[2] = p3[2] * pt[2];
191  tmp[1] = p3[1] * pt[1];
192  tmp[0] = p3[0] * pt[0];
193  double p3pt = tmp[3] - tmp[2] - tmp[1] - tmp[0];
194  tmp[3] = p2[3] * pn[3];
195  tmp[2] = p2[2] * pn[2];
196  tmp[1] = p2[1] * pn[1];
197  tmp[0] = p2[0] * pn[0];
198  double p2pn = tmp[3] - tmp[2] - tmp[1] - tmp[0];
199  tmp[3] = p3[3] * pn[3];
200  tmp[2] = p3[2] * pn[2];
201  tmp[1] = p3[1] * pn[1];
202  tmp[0] = p3[0] * pn[0];
203  double p3pn = tmp[3] - tmp[2] - tmp[1] - tmp[0];
204  tmp[3] = p2[3] * p1[3];
205  tmp[2] = p2[2] * p1[2];
206  tmp[1] = p2[1] * p1[1];
207  tmp[0] = p2[0] * p1[0];
208  double p2p1 = tmp[3] - tmp[2] - tmp[1] - tmp[0];
209  tmp[3] = p3[3] * p1[3];
210  tmp[2] = p3[2] * p1[2];
211  tmp[1] = p3[1] * p1[1];
212  tmp[0] = p3[0] * p1[0];
213  double p3p1 = tmp[3] - tmp[2] - tmp[1] - tmp[0];
214  tmp[3] = p1[3] * pn[3];
215  tmp[2] = p1[2] * pn[2];
216  tmp[1] = p1[1] * pn[1];
217  tmp[0] = p1[0] * pn[0];
218  double p1pn = tmp[3] - tmp[2] - tmp[1] - tmp[0];
219  tmp[3] = p1[3] * pt[3];
220  tmp[2] = p1[2] * pt[2];
221  tmp[1] = p1[1] * pt[1];
222  tmp[0] = p1[0] * pt[0];
223  double p1pt = tmp[3] - tmp[2] - tmp[1] - tmp[0];
224 
225  double ptsq = pt[3] * pt[3] - pt[2] * pt[2] - pt[1] * pt[1] - pt[0] * pt[0];
226  double p1sq = p1[3] * p1[3] - p1[2] * p1[2] - p1[1] * p1[1] - p1[0] * p1[0];
227 
228  // eq. 31 of https://arxiv.org/pdf/1306.1732.pdf
229  double d1, d2; //denominator parts
230  d1 = m2ee - 2.*ptp23;
231  d2 = m2ee + 2.*pk;
232  double c1, c2, c3, c4, c5, c6, c7, c8, c9, a; // 9 contractions needed
233 
234  a = (p2p3 + Mlep * Mlep); // part of l_munu, used in all contractions
235 
236  c1 = p2pt * p3pn + p3pt * p2pn - p2p3 * ptpn - ptpn * a;
237 //c1 = -(c1 * 2. + a * ptpn) * m2ee / d1 / d1; factor of 4 was missing
238  c1 = -(c1 * 2. + 4.*a * ptpn) * m2ee / d1 / d1; // ZBW 16 Jan 2021
239 
240  c2 = p2p1 * p3pn + p2pn * p3p1 - p1pn * a;
241  c2 = c2 * 4.*ptp23 / d1 / d2;
242 
243  c3 = p2pt * p3pn + p2pn * p3pt - ptpn * a;
244  c3 = c3 * 4.*ptp23 / d1 / d1;
245 
246 //c4 = 2.*p2p3 - a; factor of 4 was missing
247  c4 = 2.*p2p3 - 4.* a; // ZBW 16 Jan 2021
248  c4 = -c4 * 2.*pnk * ptp23 / d1 / d1;
249 
250  c5 = p2p1 * p3pt + p2pt * p3p1 - p1pt * a;
251  c5 = -c5 * 4.*pnk / d1 / d2;
252 
253  c6 = 2.*p2pt * p3pt - ptsq * a;
254  c6 = -c6 * 4.*pnk / d1 / d1;
255 
256  c7 = p2p1 * p3pt + p2pt * p3p1 - p1pt * a;
257  c7 = c7 * 8.*ptpn / d1 / d2;
258 
259  c8 = 2.*p2p1 * p3p1 - p1sq * a;
260  c8 = c8 * 4.*ptpn / d2 / d2;
261 
262  c9 = 2.*p2pt * p3pt - ptsq * a;
263  c9 = c9 * 4.*ptpn / d1 / d1;
264 
265  double wigner = (m2ee - m_dark * m_dark) * (m2ee - m_dark * m_dark) + m_dark * m_dark * dark_width * dark_width;
266 
267  amplit = (c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9) / wigner;
268 
269  hv[0] = 0.0;
270  hv[1] = 0.0;
271  hv[2] = 0.0;
272  hv[3] = 1.0;
273 
274  if (amplit < 0.0 && iflag == 0) {
275  //
276  (Mlep < 0.1) ? e_neg++ : m_neg++;
277  /*
278  cout << "trouble eepi: "
279  << " m2ee = " << m2ee
280  << " c1 = " << c1
281  << " c2 = " << c2
282  << " c3 = " << c3
283  << " c4 = " << c3
284  << " c5 = " << c5
285  << " c6 = " << c6
286  << " c7 = " << c7
287  << " c8 = " << c8
288  << " c9 = " << c9
289  << " c19 = " << c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9
290  << " wigner = " << wigner
291  << " amplit = " << amplit;
292  //
293  if (Mlep<0.1) {
294  cout << " e_neg = " << e_neg << " e_all = " << e_all;
295  } else {
296  cout << " m_neg = " << m_neg << " m_all = " << m_all;
297  }
298  cout << endl;
299  //
300  TLorentzVector p4t(pt); cout <<"p4tau: ";coutLorentzVector(p4t);
301  TLorentzVector p4n(pn); cout <<"p4neu: ";coutLorentzVector(p4n);
302  TLorentzVector p41(p1); cout <<"p4pi: ";coutLorentzVector(p41);
303  TLorentzVector p42(p2); cout <<"p4l+: ";coutLorentzVector(p42);
304  TLorentzVector p43(p3); cout <<"p4l-: ";coutLorentzVector(p43);
305  */
306  // If something is wrong and we need to understand what it is
307  // let's correct kinematic locally
308  double pt_new[4], pn_new[4], p1_new[4], p2_new[4], p3_new[4];
309  for (int i = 0; i < 4 ; i++) {
310  pt_new[i] = pt[i];
311  pn_new[i] = pn[i];
312  p1_new[i] = p1[i];
313  p2_new[i] = p2[i];
314  p3_new[i] = p3[i];
315  }
316  pn_new[3] = sqrt(pn[2] * pn[2] + pn[1] * pn[1] + pn[0] * pn[0]);
317  p1_new[3] = sqrt(Mpion * Mpion + p1[2] * p1[2] + p1[1] * p1[1] + p1[0] * p1[0]);
318  p2_new[3] = sqrt(Mlep * Mlep + p2[2] * p2[2] + p2[1] * p2[1] + p2[0] * p2[0]);
319  p3_new[3] = sqrt(Mlep * Mlep + p3[2] * p3[2] + p3[1] * p3[1] + p3[0] * p3[0]); // so the leptons and neutrino are on mass shell
320  for (int i = 0; i < 4 ; i++) {
321  pt_new[i] = pn_new[i] + p1_new[i] + p2_new[i] + p3_new[i]; // re-inforce energy-momenentum conservation
322  }
323  //
324  PHOB(1, pt_new, pn_new); // try as second step : this should put all 4 vectors
325  PHOB(1, pt_new, p1_new); // into tau rest frame, but neutrino may get transverse components
326  PHOB(1, pt_new, p2_new); //
327  PHOB(1, pt_new, p3_new);
328  /*
329  cout << "p4 transformed ... " << endl;
330  TLorentzVector p4t_new(pt_new);cout <<"p4tau: ";coutLorentzVector(p4t_new);
331  TLorentzVector p4n_new(pn_new);cout <<"p4neu: ";coutLorentzVector(p4n_new);
332  TLorentzVector p41_new(p1_new);cout <<"p4pi: ";coutLorentzVector(p41_new);
333  TLorentzVector p42_new(p2_new);cout <<"p4l+: ";coutLorentzVector(p42_new);
334  TLorentzVector p43_new(p3_new);cout <<"p4l-: ";coutLorentzVector(p43_new);
335  */
336  double amplit_redone;
337  iflag = 1;
338  eepi_(m_dark, dark_width, Mlep, Mpion, pt_new, pn_new, p1_new, p2_new, p3_new, amplit_redone, hv, iflag);
339  /*
340  cout << "amplit = " << amplit << " amplit_redone = " << amplit_redone ;//<< endl;
341  if (Mlep < 0.1) {
342  cout << " e_neg = " << e_neg << " e_all = " << e_all;
343  } else {
344  cout << " m_neg = " << m_neg << " m_all = " << m_all;
345  }
346  cout << endl;
347  */
348  // The last resort : for now at least
349  //
350  amplit = 0.0;
351  }
352 }
353 
354 //**************************************************************************
355 // demo mu nu_mu ee nu_tau ME, based on paper by S.Antropov https://arxiv.org/pdf/1912.11376.pdf *
356 //**************************************************************************
357 void eemu_(const double& m_dark, const double& dark_width, const double& Mtau,
358  const double* pt, const double* pn, const double* p1, const double* p2, const double* p3, const double* p4,
359  double& amplit, double* hv)
360 {
361  // we assume order: neutrino, muon, electron, positron:
362 
363  double p34[4];
364  p34[3] = p3[3] + p4[3];
365  p34[2] = p3[2] + p4[2];
366  p34[1] = p3[1] + p4[1];
367  p34[0] = p3[0] + p4[0];
368 
369  double col[4];
370  col[3] = p2[3] * p34[3];
371  col[2] = p2[2] * p34[2];
372  col[1] = p2[1] * p34[1];
373  col[0] = p2[0] * p34[0];
374 
375  double pk = col[3] - col[2] - col[1] - col[0];
376  double m2ee = p34[3] * p34[3] - p34[2] * p34[2] - p34[1] * p34[1] - p34[0] * p34[0];
377 
378  double ptk[4];
379  ptk[3] = pt[3] * p34[3];
380  ptk[2] = pt[2] * p34[2];
381  ptk[1] = pt[1] * p34[1];
382  ptk[0] = pt[0] * p34[0];
383 
384  double ptp34 = ptk[3] - ptk[2] - ptk[1] - ptk[0];
385 
386  double PP[4];
387  PP[3] = p2[3] / (pk + m2ee / 2.) - pt[3] / (ptp34 - m2ee / 2.);
388  PP[2] = p2[2] / (pk + m2ee / 2.) - pt[2] / (ptp34 - m2ee / 2.);
389  PP[1] = p2[1] / (pk + m2ee / 2.) - pt[1] / (ptp34 - m2ee / 2.);
390  PP[0] = p2[0] / (pk + m2ee / 2.) - pt[0] / (ptp34 - m2ee / 2.);
391 
392  double PPsq;
393  PPsq = PP[3] * PP[3] - PP[2] * PP[2] - PP[1] * PP[1] - PP[0] * PP[0];
394 
395  double p3PP;
396  p3PP = p3[3] * PP[3] - p3[2] * PP[2] - p3[1] * PP[1] - p3[0] * PP[0];
397 
398  double p4PP;
399  p4PP = p4[3] * PP[3] - p4[2] * PP[2] - p4[1] * PP[1] - p4[0] * PP[0];
400 
401  double wigner = (m2ee - m_dark * m_dark) * (m2ee - m_dark * m_dark) + m_dark * m_dark * dark_width * dark_width;
402 
403  double ampser = (4.*p3PP * p4PP - m2ee * PPsq) / (2.*wigner) ;
404 
405  double ak0 = 0.001 * Mtau;
406  hv[0] = 0.0;
407  hv[1] = 0.0;
408  hv[2] = 0.0;
409  hv[3] = 1.0;
410  amplit = nunul_(Mtau, p2, pn, p1, ak0, hv) * ampser;
411 
412 }
413 
414 double nunul_(const double& Mtau, const double* pl, const double* pnu, const double* pnl, double ak0, double* hv)
415 {
416 
417  double Mnu = 0.0; // 0.01
418  double Msq = Mtau * Mtau;
419  double pr[4];
420  double rxnl[3], rxnu[3], rxl[3];
421  pr[0] = 0.0;
422  pr[1] = 0.0;
423  pr[2] = 0.0;
424  pr[3] = Mtau;
425  for (int i = 0; i < 3; i++) {
426  rxnl[i] = pr[3] * pnl[3] - pr[2] * pnl[2] - pr[1] * pnl[1] - pr[0] * pnl[0];
427  rxnu[i] = pr[3] * pnu[3] - pr[2] * pnu[2] - pr[1] * pnu[1] - pr[0] * pnu[0];
428  rxl[i] = pr[3] * pl[3] - pr[2] * pl[2] - pr[1] * pl[1] - pr[0] * pl[0];
429  }
430  // QUASI TWO-BODY VARIABLES
431  double u0, u3, w3, w0, up, um, wp, wm, yu, yw, eps2, eps, y, al;
432  u0 = pl[3] / Mtau;
433  u3 = sqrt(pl[0] * pl[0] + pl[1] * pl[1] + pl[2] * pl[2]) / Mtau;
434  w3 = u3;
435  w0 = (pnu[3] + pnl[3]) / Mtau; //does this require ee pair energy?
436  up = u0 + u3;
437  um = u0 - u3;
438  wp = w0 + w3;
439  wm = w0 - w3;
440  yu = log(up / um) / 2.;
441  if (isnan(yu) || isinf(yu)) yu = 0;
442  yw = log(wp / wm) / 2.;
443  if (isnan(yw) || isinf(yw)) yw = 0;
444  eps2 = u0 * u0 - u3 * u3;
445  eps = sqrt(eps2);
446  y = w0 * w0 - w3 * w3;
447  al = ak0 / Mtau;
448 
449  // FORMFACTORS
450  double f0, fp, fm, f3;
451  f0 = 2.*u0 / u3 * (dilog(1.0 - (um * wm / (up * wp))) - dilog(1.0 - wm / wp)
452  + dilog(1.0 - um / up) - 2.*yu + 2.*log(up) * (yw + yu))
453  + 1.0 / y * (2.*u3 * yu + (1.0 - eps2 - 2 * y) * log(eps)) + 2.0 - 4.*(u0 / u3 * yu - 1.0) * log(2.*al);
454  fp = yu / (2 * u3) * (1. + (1. - eps2) / y) + log(eps) / y;
455  fm = yu / (2 * u3) * (1. - (1. - eps2) / y) - log(eps) / y;
456  f3 = eps2 * (fp + fm) / 2.;
457 
458  // SCALAR PRODUCTS OF FOUR-MOMENTA
459  double lxnl, lxnu, nuxnl, mxnl, mxnu, mxl;
460  lxnu = pl[3] * pnu[3] - pl[2] * pnu[2] - pl[1] * pnu[1] - pl[0] * pnu[0];
461  lxnl = pl[3] * pnl[3] - pl[2] * pnl[2] - pl[1] * pnl[1] - pl[0] * pnl[0];
462  nuxnl = pnu[3] * pnl[3] - pnu[2] * pnl[2] - pnu[1] * pnl[1] - pnu[0] * pnl[0];
463 
464  mxnu = Mtau * pnu[3];
465  mxnl = Mtau * pnl[3];
466  mxl = Mtau * pl[3];
467  // DECAY DIFFERENTIAL WIDTH WITHOUT AND WITH POLARIZATION
468  double c3, am3, xm3;
469  c3 = 1. ;
470  xm3 = -(f0 * lxnu * mxnl + fp * eps2 * mxnu * mxnl + fm * lxnu * lxnl + f3 * Msq * nuxnl);
471  if (isnan(xm3) || isinf(xm3)) xm3 = 0.0;
472  am3 = Mnu * xm3 * c3;
473 
474  // V-A AND V+A COUPLINGS, BUT IN THE BORN PART ONLY
475  double brak, gv, ga, born, xm3pol[3], am3pol[3], bornp[3];
476  gv = 1.;
477  ga = -1.;
478  brak = (gv + ga) * (gv + ga) * mxl * nuxnl
479  + (gv - ga) * (gv - ga) * mxnl * lxnu
480  - (gv * gv - ga * ga) * Mtau * Mnu * lxnl;
481  born = brak;
482  for (int i = 0; i < 3; i++) {
483  xm3pol[i] = -(f0 * lxnu * rxnl[i] + fp * eps2 * mxnu * rxnl[i]
484  + fm * lxnu * (lxnl + (rxnl[i] * mxl - mxnl * rxl[i]) / Msq)
485  + f3 * (Msq * nuxnl + mxnu * rxnl[i] - rxnu[i] * mxnl));
486  if (isnan(xm3pol[i]) || isinf(xm3pol[i])) xm3pol[i] = 0.0;
487  am3pol[i] = Mnu * xm3pol[i] * c3;
488  // V-A AND V+A COUPLINGS, BUT IN THE BORN PART ONLY
489  bornp[i] = born + ((gv + ga) * (gv + ga) * Mtau * nuxnl * pl[i]
490  - (gv + ga) * (gv + ga) * Mtau * lxnu * pnu[i]
491  + (gv * gv - ga * ga) * Mnu * mxnl * pl[i]
492  - (gv * gv - ga * ga) * Mnu * mxl * pnu[i]);
493  hv[i] = (bornp[i] + am3pol[i]) / (born + am3) - 1.0;
494  }
495 
496  double tbh;
497  tbh = born + am3;
498  if (isnan(tbh) || isinf(tbh)) {
499  // cout << f0 << yw << wm << wp << " " << w0 << " " << w3 << " " << xm3pol[3] << endl;
500  tbh = 0.0;
501  return tbh;
502  }
503  if (tbh / born < 0.1) {
504  // cout << "ERROR in nunul (matrix element)" << tbh / born << endl;
505  tbh = 0.0;
506  return tbh;
507  }
508 // if(tbh!=0.0) cout<<"amplit= "<< tbh << endl;
509  return tbh;
510 }
511