evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtbsToLLLLAmp.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
21 #include "EvtGenModels/EvtbsToLLLLAmp.hh"
22 
23 #include "EvtGenBase/EvtAmp.hh"
24 #include "EvtGenBase/EvtComplex.hh"
25 #include "EvtGenBase/EvtDiracSpinor.hh"
26 #include "EvtGenBase/EvtGenKine.hh"
27 #include "EvtGenBase/EvtId.hh"
28 #include "EvtGenBase/EvtIdSet.hh"
29 #include "EvtGenBase/EvtPDL.hh"
30 #include "EvtGenBase/EvtParticle.hh"
31 #include "EvtGenBase/EvtPatches.hh"
32 #include "EvtGenBase/EvtReport.hh"
33 #include "EvtGenBase/EvtScalarParticle.hh"
34 #include "EvtGenBase/EvtTensor4C.hh"
35 #include "EvtGenBase/EvtVector4C.hh"
36 
37 #include "EvtGenModels/EvtbTosllWilsCoeffNLO.hh"
38 #include "EvtGenModels/Evtbs2llGammaFFMNT.hh"
39 
40 #include <cstdlib>
41 
42 // input: *parent - the pointer to the parent particle (B-meson, the
43 // object of the EvtParticle class);
44 // *formFactors - the pointer to instance of EvtbTosllGammaFF class object;
45 // *WilsCoeff - the pointer to the Standart Model Wilson Coefficients class;
46 // mu - the scale parameter, GeV;
47 // Nf - number of "effective" flavors (for b-quark Nf=5);
48 // res_swch - resonant switching parameter:
49 // = 0 the resonant contribution switched OFF,
50 // = 1 the resonant contribution switched ON;
51 // ias - switching parameter for \alpha_s(M_Z) value:
52 // = 0 PDG 1sigma minimal alpha_s(M_Z),
53 // = 1 PDG average value alpha_s(M_Z),
54 // = 2 PDG 1sigma maximal alpha_s(M_Z).
55 // Wolfenstein parameterization for CKM matrix
56 // CKM_A, CKM_lambda, CKM_barrho, CKM_bareta
57 
58 void EvtbsToLLLLAmp::CalcAmp( EvtParticle* parent, EvtAmp& amp,
59  Evtbs2llGammaFF* formFactors,
60  EvtbTosllWilsCoeffNLO* WilsCoeff, double mu,
61  int Nf, int res_swch, int ias, double CKM_A,
62  double CKM_lambda, double CKM_barrho,
63  double CKM_bareta )
64 {
65  // FILE *mytest;
66 
67  int il1 = 0, il2 = 1, il3 = 2,
68  il4 = 3; // leptons are the first, second, thirds
69  // and fourth daughter particles
70 
71  EvtComplex unit1( 1.0, 0.0 ); // real unit
72  EvtComplex uniti( 0.0, 1.0 ); // imaginary unit
73 
74  double M1 = parent->mass(); // B - meson mass, GeV
75  double ml = parent->getDaug( il1 )->mass(); // leptonic mass, GeV
76  double mq = 0.0; // light quark mass from the dispersion QM, GeV
77  double mc = formFactors->getQuarkMass(
78  4 ); // m_c mass from the dispersion QM, GeV
79  double mb = formFactors->getQuarkMass(
80  5 ); // m_b mass from the dispersion QM, GeV
81  double Mw = 80.403; // GeV W-boson mass, GeV
82  double mt = 174.2; // GeV t-quark mass, GeV
83  double fb = 0.0; // leptonic decay constant of B-meson, Gev
84 
85  EvtComplex Vtb, Vtq, Vub, Vuq, Vcb,
86  Vcq; // V_{tb}, V_{tq}, V_{ub}, V_{uq}, V_{cb}, V_{cq}
87  EvtComplex CKM_factor; // V^*_{tq}*V_{tb}, where q={d,s}
88  EvtComplex lambda_qu; // V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}, where q={d,s}
89  EvtComplex lambda_qc; // V^*_{cq}*V_{cb}/V^*_{tq}*V_{tb}, where q={d,s}
90  double Relambda_qu, Imlambda_qu;
91 
92  //
93  // Setting of the mq and CKM matrix elements for different Bq-mesons tipes
94  //
95  EvtId idparent = parent->getId(); // Bq-meson Id
96  EvtId IdMu1, IdMu2, IdMu3, IdMu4;
97 
98  if ( idparent == EvtPDL::getId( std::string( "B_s0" ) ) ||
99  idparent == EvtPDL::getId( std::string( "anti-B_s0" ) ) ) {
100  mq = formFactors->getQuarkMass( 3 ); // m_s mass from the dispersion QM
101  fb = 0.24; // leptonic decay constant
102  // V_{ts}
103  Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
104  pow( CKM_lambda, 2.0 ) *
105  ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
106  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
107  Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
108  // V_{us}
109  Vuq = CKM_lambda * unit1;
110  // V_{cs}
111  Vcq = unit1 - 0.5 * pow( CKM_lambda, 2.0 ) -
112  0.125 * pow( CKM_lambda, 4.0 ) * ( 1.0 + 4.0 * pow( CKM_A, 2.0 ) );
113  }
114 
115  if ( idparent == EvtPDL::getId( std::string( "B0" ) ) ||
116  idparent == EvtPDL::getId( std::string( "anti-B0" ) ) ) {
117  mq = formFactors->getQuarkMass( 2 ); // m_d mass from the dispersion QM
118  fb = 0.20; // leptonic decay constant
119  // V_{td}
120  Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
121  ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
122  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
123  Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
124  // V_{ud}
125  Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
126  0.125 * pow( CKM_lambda, 4.0 ) );
127  // V_{cd}
128  Vcq = unit1 *
129  ( -CKM_lambda +
130  0.5 * pow( CKM_A, 2.0 ) * pow( CKM_lambda, 5.0 ) *
131  ( 1.0 - 2.0 * ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
132  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ) ) );
133  }
134 
135  if ( mq < 0.001 ) {
136  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
137  << "\n\n The function EvtbsToLLLLAmp::CalcAmp(...)"
138  << "\n Error in the mq setting!"
139  << "\n mq = " << mq << "< 0.001"
140  << "\n idparent = " << idparent << std::endl;
141  ::abort();
142  }
143 
144  Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
145  2.0 ) ); // V_{tb}
146  Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
147  ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
148  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); // V_{ub}
149  Vcb = unit1 * CKM_A * pow( CKM_lambda, 2.0 ); // V_{cb}
150 
151  CKM_factor = conj( Vtq ) * Vtb; // V^*_{tq}*V_{tb}
152 
153  lambda_qu = conj( Vuq ) * Vub /
154  CKM_factor; // V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}
155  Relambda_qu = real( lambda_qu );
156  Imlambda_qu = imag( lambda_qu );
157 
158  lambda_qc = conj( Vcq ) * Vcb /
159  CKM_factor; // V^*_{cq}*V_{cb}/V^*_{tq}*V_{tb}
160 
161  //
162  // Setting the leptonic kinematical properties
163  //
164 
165  // to find charges of ell^+ and ell^- in the B-meson daughters
166  int charge1 = ( EvtPDL::chg3( parent->getDaug( il1 )->getId() ) ) / 3;
167  int charge2 = ( EvtPDL::chg3( parent->getDaug( il2 )->getId() ) ) / 3;
168  int charge3 = ( EvtPDL::chg3( parent->getDaug( il3 )->getId() ) ) / 3;
169  int charge4 = ( EvtPDL::chg3( parent->getDaug( il4 )->getId() ) ) / 3;
170  if ( ( abs( charge1 ) != 1 ) || ( abs( charge2 ) != 1 ) ||
171  ( abs( charge3 ) != 1 ) || ( abs( charge4 ) != 1 ) ||
172  ( charge1 + charge2 + charge3 + charge4 != 0 ) ) {
173  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
174  << "\n\n The function EvtbsToLLLLAmp::CalcAmp(...)"
175  << "\n Error in the leptonic charge definition!"
176  << "\n charge1 =" << charge1
177  << "\n charge2 =" << charge2
178  << "\n charge3 =" << charge3
179  << "\n charge4 =" << charge4
180  << "\n number of daughters =" << parent->getNDaug() << std::endl;
181  ::abort();
182  }
183 
184  EvtParticle* lep1Plus = 0;
185  EvtParticle* lep1Minus = 0;
186  EvtParticle* lep2Plus = 0;
187  EvtParticle* lep2Minus = 0;
188 
189  EvtVector4R p; // B-meson momentum in the B-rest frame
190  EvtVector4R q; // first transition 4-momentum in the B-rest frame
191  EvtVector4R k; // second transition 4-momentum in the B-rest frame
192 
193  double q2; // Mandelstam variable s=q^2
194 
195  // Nondimentional 4-momentums
196  EvtVector4R hatp;
197  EvtVector4R hatq;
198  EvtVector4R hatk;
199 
200  EvtVector4R qsecond; // first transition 4-momentum in the B-rest frame
201  EvtVector4R ksecond; // second transition 4-momentum in the B-rest frame
202 
203  double q2second; // Mandelstam variable s=q^2
204 
205  // Nondimentional 4-momentums
206  EvtVector4R hatpsecond;
207  EvtVector4R hatqsecond;
208  EvtVector4R hatksecond;
209 
210  EvtVector4R k_1; // 4-momentum of ell^+ in the B-rest frame
211  EvtVector4R k_2; // 4-momentum of ell^- in the B-rest frame
212  EvtVector4R k_3; // 4-momentum of ell^+ in the B-rest frame
213  EvtVector4R k_4; // 4-momentum of ell^- in the B-rest frame
214 
215  k_1.set( 0.0, 0.0, 0.0, 0.0 );
216  k_2.set( 0.0, 0.0, 0.0, 0.0 );
217  k_3.set( 0.0, 0.0, 0.0, 0.0 );
218  k_4.set( 0.0, 0.0, 0.0, 0.0 );
219 
220  if ( ( charge1 + charge2 == 0 ) && ( charge3 + charge4 == 0 ) ) {
221  // positive charged lepton 1
222  lep1Plus = ( charge1 > charge2 ) ? parent->getDaug( il1 )
223  : parent->getDaug( il2 );
224  // negative charged lepton 1
225  lep1Minus = ( charge1 < charge2 ) ? parent->getDaug( il1 )
226  : parent->getDaug( il2 );
227  if ( charge1 > charge2 ) {
228  k_1 = parent->getDaug( il1 )->getP4();
229  k_2 = parent->getDaug( il2 )->getP4();
230  IdMu1 = parent->getDaug( il1 )->getId();
231  IdMu2 = parent->getDaug( il2 )->getId();
232  } else {
233  k_1 = parent->getDaug( il2 )->getP4();
234  k_2 = parent->getDaug( il1 )->getP4();
235  IdMu1 = parent->getDaug( il2 )->getId();
236  IdMu2 = parent->getDaug( il1 )->getId();
237  }
238  // positive charged lepton 2
239  lep2Plus = ( charge3 > charge4 ) ? parent->getDaug( il3 )
240  : parent->getDaug( il4 );
241  // negative charged lepton 2
242  lep2Minus = ( charge3 < charge4 ) ? parent->getDaug( il3 )
243  : parent->getDaug( il4 );
244  if ( charge3 > charge4 ) {
245  k_3 = parent->getDaug( il3 )->getP4();
246  k_4 = parent->getDaug( il4 )->getP4();
247  IdMu3 = parent->getDaug( il3 )->getId();
248  IdMu4 = parent->getDaug( il4 )->getId();
249  } else {
250  k_3 = parent->getDaug( il4 )->getP4();
251  k_4 = parent->getDaug( il3 )->getP4();
252  IdMu3 = parent->getDaug( il4 )->getId();
253  IdMu4 = parent->getDaug( il3 )->getId();
254  }
255  }
256  if ( ( charge1 + charge3 == 0 ) && ( charge2 + charge4 == 0 ) ) {
257  // positive charged lepton 1
258  lep1Plus = ( charge1 > charge3 ) ? parent->getDaug( il1 )
259  : parent->getDaug( il3 );
260  // negative charged lepton 1
261  lep1Minus = ( charge1 < charge3 ) ? parent->getDaug( il1 )
262  : parent->getDaug( il3 );
263  if ( charge1 > charge3 ) {
264  k_1 = parent->getDaug( il1 )->getP4();
265  k_2 = parent->getDaug( il3 )->getP4();
266  IdMu1 = parent->getDaug( il1 )->getId();
267  IdMu2 = parent->getDaug( il3 )->getId();
268  } else {
269  k_1 = parent->getDaug( il3 )->getP4();
270  k_2 = parent->getDaug( il1 )->getP4();
271  IdMu1 = parent->getDaug( il3 )->getId();
272  IdMu2 = parent->getDaug( il1 )->getId();
273  }
274  // positive charged lepton 2
275  lep2Plus = ( charge2 > charge4 ) ? parent->getDaug( il2 )
276  : parent->getDaug( il4 );
277  // negative charged lepton 2
278  lep2Minus = ( charge2 < charge4 ) ? parent->getDaug( il2 )
279  : parent->getDaug( il4 );
280  if ( charge2 > charge4 ) {
281  k_3 = parent->getDaug( il2 )->getP4();
282  k_4 = parent->getDaug( il4 )->getP4();
283  IdMu3 = parent->getDaug( il2 )->getId();
284  IdMu4 = parent->getDaug( il4 )->getId();
285  } else {
286  k_3 = parent->getDaug( il4 )->getP4();
287  k_4 = parent->getDaug( il2 )->getP4();
288  IdMu3 = parent->getDaug( il4 )->getId();
289  IdMu4 = parent->getDaug( il2 )->getId();
290  }
291  }
292 
293  p = parent->getP4Restframe();
294  hatp = p / M1;
295 
296  //
297  // The calculation of the FIRST part of the amplitude
298  //
299 
300  q = k_1 + k_2;
301  k = k_3 + k_4;
302  q2 = q.mass2(); // Mandelstam variable s=q^2
303  hatq = q / M1;
304  hatk = k / M1;
305 
306  // The Wilson Coefficients preparation according to the paper
307  // A.J.Buras, M.Munz, Phys.Rev.D52, p.189 (1995)
308  double c1, c2;
309  EvtComplex a1, c7gam, c9eff_b2q, c9eff_barb2barq, c10a;
310 
311  // Excluded of the J/psi and psi' resonant area
312  if ( ( res_swch == 1 ) && ( q2 >= 9.199 ) && ( q2 <= 15.333 ) ) {
313  c1 = 0.0;
314  c2 = 0.0;
315  a1 = unit1 * 0.0;
316  c7gam = unit1 * 0.0;
317  c9eff_b2q = unit1 * 0.0;
318  c9eff_barb2barq = unit1 * 0.0;
319  c10a = unit1 * 0.0;
320  } else {
321  c1 = WilsCoeff->C1( mu, Mw, Nf, ias );
322  c2 = WilsCoeff->C2( mu, Mw, Nf, ias );
323  a1 = unit1 * ( c1 + c2 / 3.0 );
324  c7gam = WilsCoeff->GetC7Eff( mu, Mw, mt, Nf, ias );
325  c9eff_b2q = WilsCoeff->GetC9Eff( 0, res_swch, ias, Nf, q2, mb, mq, mc, mu,
326  mt, Mw, ml, Relambda_qu, Imlambda_qu );
327  c9eff_barb2barq = WilsCoeff->GetC9Eff( 1, res_swch, ias, Nf, q2, mb, mq,
328  mc, mu, mt, Mw, ml, Relambda_qu,
329  Imlambda_qu );
330  c10a = WilsCoeff->GetC10Eff( mt, Mw );
331  }
332 
333  EvtComplex Fv,
334  Fa; // The change of the sign is included in the amplitudes definition!
335  EvtComplex Ftv_b2q, Ftv_barb2barq;
336  EvtComplex Fta_b2q, Fta_barb2barq;
337 
338  // Excluded of the J/psi and psi' resonant area
339  if ( ( res_swch == 1 ) && ( q2 >= 9.199 ) && ( q2 <= 15.333 ) ) {
340  fb = 0.0;
341  Fa = unit1 * 0.0;
342  Fv = unit1 * 0.0;
343  Fta_b2q = unit1 * 0.0;
344  Fta_barb2barq = unit1 * 0.0;
345  Ftv_b2q = unit1 * 0.0;
346  Ftv_barb2barq = unit1 * 0.0;
347  } else {
348  if ( fb < 0.01 ) {
349  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
350  << "\n\n The function EvtbsToLLLLAmp::CalcAmp(...)"
351  << "\n Leptonic decay constant fb is not uninitialized in this function!"
352  << " fb = " << fb << std::endl;
353  ::abort();
354  }
355 
356  // For \bar B^0_q -> l^+ l^- gamma
357  formFactors->getPhotonFF( 0, fb, parent->getId(), q2, M1, mb, mq, c7gam,
358  a1, lambda_qu, lambda_qc, Fv, Fa, Ftv_b2q,
359  Fta_b2q );
360 
361  // For B^0_q -> l^+ l^- gamma
362  formFactors->getPhotonFF( 1, fb, parent->getId(), q2, M1, mb, mq, c7gam,
363  a1, lambda_qu, lambda_qc, Fv, Fa,
364  Ftv_barb2barq, Fta_barb2barq );
365  }
366 
367  // The functions for the hadronic matrix element calculation
368  EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c_b2q, c_barb2barq;
369  EvtComplex e_b2q, e_barb2barq, f_b2q, f_barb2barq, g_b2q, g_barb2barq;
370 
371  a_b2q = c9eff_b2q * Fv + 2.0 * c7gam * Ftv_b2q * mb * M1 / q2;
372  a_barb2barq = c9eff_barb2barq * Fv +
373  2.0 * c7gam * Ftv_barb2barq * mb * M1 / q2;
374 
375  b_b2q = ( c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2 ) *
376  ( hatp * hatk );
377  b_barb2barq = ( c9eff_barb2barq * Fa +
378  2.0 * c7gam * Fta_barb2barq * mb * M1 / q2 ) *
379  ( hatp * hatk );
380 
381  c_b2q = c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2;
382  c_barb2barq = c9eff_barb2barq * Fa +
383  2.0 * c7gam * Fta_barb2barq * mb * M1 / q2;
384 
385  e_b2q = c10a * Fv;
386  e_barb2barq = e_b2q;
387 
388  f_b2q = c10a * Fa * ( hatp * hatk );
389  f_barb2barq = f_b2q;
390 
391  g_b2q = c10a * Fa;
392  g_barb2barq = g_b2q;
393 
394  //
395  // The calculation of the SECOND part of the amplitude
396  //
397 
398  qsecond = k_1 + k_4;
399  ksecond = k_3 + k_2;
400  q2second = qsecond.mass2(); // Mandelstam variable s=q^2
401  hatqsecond = qsecond / M1;
402  hatksecond = ksecond / M1;
403 
404  // Excluded of the J/psi and psi' resonant area
405  if ( ( res_swch == 1 ) && ( q2second >= 9.199 ) && ( q2second <= 15.333 ) ) {
406  c1 = 0.0;
407  c2 = 0.0;
408  a1 = unit1 * 0.0;
409  c7gam = unit1 * 0.0;
410  c9eff_b2q = unit1 * 0.0;
411  c9eff_barb2barq = unit1 * 0.0;
412  c10a = unit1 * 0.0;
413  } else {
414  c1 = WilsCoeff->C1( mu, Mw, Nf, ias );
415  c2 = WilsCoeff->C2( mu, Mw, Nf, ias );
416  a1 = unit1 * ( c1 + c2 / 3.0 );
417  c7gam = WilsCoeff->GetC7Eff( mu, Mw, mt, Nf, ias );
418  c9eff_b2q = WilsCoeff->GetC9Eff( 0, res_swch, ias, Nf, q2second, mb, mq,
419  mc, mu, mt, Mw, ml, Relambda_qu,
420  Imlambda_qu );
421  c9eff_barb2barq = WilsCoeff->GetC9Eff( 1, res_swch, ias, Nf, q2second,
422  mb, mq, mc, mu, mt, Mw, ml,
423  Relambda_qu, Imlambda_qu );
424  c10a = WilsCoeff->GetC10Eff( mt, Mw );
425  }
426 
427  // Excluded of the J/psi and psi' resonant area
428  if ( ( res_swch == 1 ) && ( q2second >= 9.199 ) && ( q2second <= 15.333 ) ) {
429  fb = 0.0;
430  Fa = unit1 * 0.0;
431  Fv = unit1 * 0.0;
432  Fta_b2q = unit1 * 0.0;
433  Fta_barb2barq = unit1 * 0.0;
434  Ftv_b2q = unit1 * 0.0;
435  Ftv_barb2barq = unit1 * 0.0;
436  } else {
437  if ( fb < 0.01 ) {
438  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
439  << "\n\n The function EvtbsToLLLLAmp::CalcAmp(...)"
440  << "\n Leptonic decay constant fb is not uninitialized in this function!"
441  << " fb = " << fb << std::endl;
442  ::abort();
443  }
444 
445  // For \bar B^0_q -> l^+ l^- gamma
446  formFactors->getPhotonFF( 0, fb, parent->getId(), q2second, M1, mb, mq,
447  c7gam, a1, lambda_qu, lambda_qc, Fv, Fa,
448  Ftv_b2q, Fta_b2q );
449 
450  // For B^0_q -> l^+ l^- gamma
451  formFactors->getPhotonFF( 1, fb, parent->getId(), q2second, M1, mb, mq,
452  c7gam, a1, lambda_qu, lambda_qc, Fv, Fa,
453  Ftv_barb2barq, Fta_barb2barq );
454  }
455 
456  // The functions for the hadronic matrix element calculation
457  EvtComplex a_b2qsecond, a_barb2barqsecond, b_b2qsecond, b_barb2barqsecond,
458  c_b2qsecond, c_barb2barqsecond;
459  EvtComplex e_b2qsecond, e_barb2barqsecond, f_b2qsecond, f_barb2barqsecond,
460  g_b2qsecond, g_barb2barqsecond;
461 
462  a_b2qsecond = c9eff_b2q * Fv + 2.0 * c7gam * Ftv_b2q * mb * M1 / q2second;
463  a_barb2barqsecond = c9eff_barb2barq * Fv +
464  2.0 * c7gam * Ftv_barb2barq * mb * M1 / q2second;
465 
466  b_b2qsecond = ( c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2second ) *
467  ( hatpsecond * hatksecond );
468  b_barb2barqsecond = ( c9eff_barb2barq * Fa +
469  2.0 * c7gam * Fta_barb2barq * mb * M1 / q2second ) *
470  ( hatpsecond * hatksecond );
471 
472  c_b2qsecond = c9eff_b2q * Fa + 2.0 * c7gam * Fta_b2q * mb * M1 / q2second;
473  c_barb2barqsecond = c9eff_barb2barq * Fa +
474  2.0 * c7gam * Fta_barb2barq * mb * M1 / q2second;
475 
476  e_b2qsecond = c10a * Fv;
477  e_barb2barqsecond = e_b2qsecond;
478 
479  f_b2qsecond = c10a * Fa * ( hatpsecond * hatksecond );
480  f_barb2barqsecond = f_b2qsecond;
481 
482  g_b2qsecond = c10a * Fa;
483  g_barb2barqsecond = g_b2qsecond;
484 
485  EvtTensor4C T1, T2; // Tensor structures for
486  EvtTensor4C T1second,
487  T2second; // the hadronic matrix element calculation
488 
489  // B - and barB - mesons descriptors
490  static EvtIdSet bmesons( "anti-B0", "anti-B_s0" );
491  static EvtIdSet bbarmesons( "B0", "B_s0" );
492 
493  EvtId parentID = parent->getId();
494 
495  if ( bmesons.contains( parentID ) ) {
496  // The amplitude for the decay barB -> gamma ell^+ ell^- or
497  // b \bar q -> gamma ell^+ ell^-
498 
499  T1 = a_b2q * unit1 * dual( EvtGenFunctions::directProd( hatq, hatk ) ) -
500  b_b2q * uniti * EvtTensor4C::g() +
501  c_b2q * uniti * EvtGenFunctions::directProd( hatk, hatq );
502 
503  T2 = e_b2q * unit1 * dual( EvtGenFunctions::directProd( hatp, hatk ) ) -
504  f_b2q * uniti * EvtTensor4C::g() +
505  g_b2q * uniti * EvtGenFunctions::directProd( hatk, hatq );
506 
507  T1second = a_b2qsecond * unit1 *
508  dual( EvtGenFunctions::directProd( hatqsecond,
509  hatksecond ) ) -
510  b_b2qsecond * uniti * EvtTensor4C::g() +
511  c_b2qsecond * uniti *
512  EvtGenFunctions::directProd( hatksecond, hatqsecond );
513 
514  T2second = e_b2qsecond * unit1 *
515  dual( EvtGenFunctions::directProd( hatpsecond,
516  hatksecond ) ) -
517  f_b2qsecond * uniti * EvtTensor4C::g() +
518  g_b2qsecond * uniti *
519  EvtGenFunctions::directProd( hatksecond, hatqsecond );
520 
521  int i1, i2, i3, i4; // leptonic spin structures counters
522  int leptonicspin[4]; // array for the saving of the leptonic spin configuration
523 
524  // Tables for correspondings
525  // l^+(k_1) && lep1Plus && k_1 && i1
526  // l^-(k_2) && lep1Minus && k_2 && i2
527  // l^+(k_3) && lep2Plus && k_3 && i3
528  // l^-(k_4) && lep2Minus && k_4 && i4
529 
530  for ( i2 = 0; i2 < 2; i2++ ) {
531  leptonicspin[0] = i2;
532  for ( i1 = 0; i1 < 2; i1++ ) {
533  leptonicspin[1] = i1;
534  for ( i4 = 0; i4 < 2; i4++ ) {
535  leptonicspin[2] = i4;
536  for ( i3 = 0; i3 < 2; i3++ ) {
537  leptonicspin[3] = i3;
538  EvtVector4C VL2L1, AL2L1, VL4L3;
539  EvtVector4C E1, E2;
540  EvtVector4C VL2L1second, AL2L1second, VL4L3second;
541  EvtVector4C E1second, E2second;
542 
543  VL2L1 = EvtLeptonVCurrent( lep1Minus->spParent( i2 ),
544  lep1Plus->spParent( i1 ) );
545  AL2L1 = EvtLeptonACurrent( lep1Minus->spParent( i2 ),
546  lep1Plus->spParent( i1 ) );
547  VL4L3 = EvtLeptonVCurrent( lep2Minus->spParent( i4 ),
548  lep2Plus->spParent( i3 ) );
549  E1 = T1.cont2( VL4L3 );
550  E2 = T2.cont2( VL4L3 );
551 
552  VL2L1second = EvtLeptonVCurrent(
553  lep2Minus->spParent( i2 ), lep1Plus->spParent( i1 ) );
554  AL2L1second = EvtLeptonACurrent(
555  lep2Minus->spParent( i2 ), lep1Plus->spParent( i1 ) );
556  VL4L3second = EvtLeptonVCurrent(
557  lep1Minus->spParent( i4 ), lep2Plus->spParent( i3 ) );
558  E1second = T1second.cont2( VL4L3second );
559  E2second = T2second.cont2( VL4L3second );
560 
561  amp.vertex( leptonicspin,
562  CKM_factor * ( VL2L1 * E1 + AL2L1 * E2 +
563  VL2L1second * E1second +
564  AL2L1second * E2second ) );
565 
566  // EvtGenReport(EVTGEN_ERROR,"EvtGen")
567  // << "\n\n ============================================================================"
568  // << "\n The matrix element (first + second) = "
569  // << CKM_factor*(VL2L1*E1+AL2L1*E2+VL2L1second*E1second+AL2L1second*E2second)
570  // << "\n The matrix element (only first) = "
571  // << CKM_factor*(VL2L1*E1+AL2L1*E2)
572  // << "============================================================================\n\n"
573  // << std::endl;
574  }
575  }
576  }
577  }
578 
579  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n The function EvtbsToLLLLAmp::CalcAmp(...) passed with arguments:"
580  // << "\n ============================================================================"
581  // << "\n Input parameters:"
582  // << "\n mu = " << mu
583  // << "\n Nf =" << Nf
584  // << "\n res_swch = " << res_swch
585  // << "\n ias = " << ias
586  // << "\n CKM_A = " << CKM_A
587  // << "\n CKM_lambda = " << CKM_lambda
588  // << "\n CKM_barrho = " << CKM_barrho
589  // << "\n CKM_bareta = " << CKM_bareta
590  // << "\n CKM_factor = " << CKM_factor
591  // << "\n ============================================================================"
592  // << "\n Kinematics:"
593  // << "\n k_1 = " << k_1
594  // << "\n m_ell_1 =" << parent->getDaug(il1)->mass()
595  // << "\n k_2 = " << k_2
596  // << "\n m_ell_2 =" << parent->getDaug(il2)->mass()
597  // << "\n k_3 = " << k_3
598  // << "\n m_ell_3 =" << parent->getDaug(il3)->mass()
599  // << "\n k_4 = " << k_4
600  // << "\n m_ell_4 =" << parent->getDaug(il4)->mass()
601  // << "\n p = " << p
602  // << "\n q = " << q
603  // << "\n k = " << k
604  // << "\n ============================================================================"
605  // << "\n Form-factors"
606  // << "\n Fv = " << Fv
607  // << "\n Fa = " << Fa
608  // << "\n Ftv_b2q = " << Ftv_b2q
609  // << "\n Fta_b2q = " << Fta_b2q
610  // << "\n Ftv_barb2barq = " << Ftv_barb2barq
611  // << "\n Fta_barb2barq = " << Fta_barb2barq
612  // << "\n fb = " << fb
613  // << "\n ============================================================================"
614  // << "\n Wilson Coefficients:"
615  // << "\n Re(c7gam) = " << real(c7gam)
616  // << " Im(c7gam) = " << imag(c7gam)
617  // << "\n Re(c9eff_b2q) = " << real(c9eff_b2q)
618  // << " Im(c9eff_b2q) = " << imag(c9eff_b2q)
619  // << "\n Re(c9eff_barb2barq) = " << real(c9eff_barb2barq)
620  // << " Im(c9eff_barb2barq) = " << imag(c9eff_barb2barq)
621  // << "\n Re(c10a) = " << real(c10a)
622  // << " Im(c10a) = " << imag(c10a)
623  // << "\n ============================================================================"
624  // << "\n Functions in the matrix element:"
625  // << "\n a_b2q = " << a_b2q
626  // << "\n b_b2q = " << b_b2q
627  // << "\n c_b2q = " << c_b2q
628  // << "\n e_b2q = " << e_b2q
629  // << "\n f_b2q = " << f_b2q
630  // << "\n g_b2q = " << g_b2q
631  // << "\n ============================================================================"
632  // << "\n Partical Properties:"
633  // << "\n IdB = " << idparent << " == " << EvtPDL::getId(std::string("anti-B_s0"))
634  // << "\n IdMu1 = " << IdMu1 << " == " << EvtPDL::getId(std::string("mu+"))
635  // << "\n IdMu2 = " << IdMu2 << " == " << EvtPDL::getId(std::string("mu-"))
636  // << "\n IdMu3 = " << IdMu3 << " == " << EvtPDL::getId(std::string("mu+"))
637  // << "\n IdMu4 = " << IdMu4 << " == " << EvtPDL::getId(std::string("mu-"))
638  // << "\n\n\n\n"
639  // << std::endl;
640 
641  } else {
642  if ( bbarmesons.contains( parentID ) ) {
643  // The amplitude for the decay B -> gamma ell^+ ell^- or
644  // q bar b -> gamma ell^+ ell^-
645 
646  T1 = -a_barb2barq * unit1 *
647  dual( EvtGenFunctions::directProd( hatq, hatk ) ) -
648  b_barb2barq * uniti * EvtTensor4C::g() +
649  c_barb2barq * uniti * EvtGenFunctions::directProd( hatk, hatq );
650 
651  T2 = -e_barb2barq * unit1 *
652  dual( EvtGenFunctions::directProd( hatq, hatk ) ) -
653  f_barb2barq * uniti * EvtTensor4C::g() +
654  g_barb2barq * uniti * EvtGenFunctions::directProd( hatk, hatq );
655 
656  T1second = -a_barb2barqsecond * unit1 *
657  dual( EvtGenFunctions::directProd( hatqsecond,
658  hatksecond ) ) -
659  b_barb2barqsecond * uniti * EvtTensor4C::g() +
660  c_barb2barqsecond * uniti *
661  EvtGenFunctions::directProd( hatksecond, hatqsecond );
662 
663  T2second = -e_barb2barqsecond * unit1 *
664  dual( EvtGenFunctions::directProd( hatpsecond,
665  hatksecond ) ) -
666  f_barb2barqsecond * uniti * EvtTensor4C::g() +
667  g_barb2barqsecond * uniti *
668  EvtGenFunctions::directProd( hatksecond, hatqsecond );
669 
670  int i1, i2, i3, i4; // leptonic spin structures counters
671  int leptonicspin[4]; // array for the saving of the leptonic spin configuration
672 
673  // Tables for correspondings
674  // l^+(k_1) && lep1Plus && k_1 && i1
675  // l^-(k_2) && lep1Minus && k_2 && i2
676  // l^+(k_3) && lep2Plus && k_3 && i3
677  // l^-(k_4) && lep2Minus && k_4 && i4
678 
679  for ( i2 = 1; i2 < 0; i2-- ) {
680  leptonicspin[0] = i2;
681  for ( i1 = 1; i1 < 0; i1-- ) {
682  leptonicspin[1] = i1;
683  for ( i4 = 1; i4 < 0; i4-- ) {
684  leptonicspin[2] = i4;
685  for ( i3 = 1; i3 < 0; i3-- ) {
686  leptonicspin[3] = i3;
687  EvtVector4C VL2L1, AL2L1, VL4L3;
688  EvtVector4C E1, E2;
689  EvtVector4C VL2L1second, AL2L1second, VL4L3second;
690  EvtVector4C E1second, E2second;
691 
692  VL2L1 = EvtLeptonVCurrent( lep1Minus->spParent( i2 ),
693  lep1Plus->spParent( i1 ) );
694  AL2L1 = EvtLeptonACurrent( lep1Minus->spParent( i2 ),
695  lep1Plus->spParent( i1 ) );
696  VL4L3 = EvtLeptonVCurrent( lep2Minus->spParent( i4 ),
697  lep2Plus->spParent( i3 ) );
698  E1 = T1.cont2( VL4L3 );
699  E2 = T2.cont2( VL4L3 );
700 
701  VL2L1second =
702  EvtLeptonVCurrent( lep2Minus->spParent( i2 ),
703  lep1Plus->spParent( i1 ) );
704  AL2L1second =
705  EvtLeptonACurrent( lep2Minus->spParent( i2 ),
706  lep1Plus->spParent( i1 ) );
707  VL4L3second =
708  EvtLeptonVCurrent( lep1Minus->spParent( i4 ),
709  lep2Plus->spParent( i3 ) );
710  E1second = T1second.cont2( VL4L3second );
711  E2second = T2second.cont2( VL4L3second );
712 
713  amp.vertex( leptonicspin,
714  conj( CKM_factor ) *
715  ( VL2L1 * E1 + AL2L1 * E2 +
716  VL2L1second * E1second +
717  AL2L1second * E2second ) );
718  }
719  }
720  }
721  }
722 
723  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n The function EvtbsToLLLLAmp::CalcAmp(...) passed with arguments:"
724  // << "\n ============================================================================"
725  // << "\n Input parameters:"
726  // << "\n mu = " << mu
727  // << "\n Nf =" << Nf
728  // << "\n res_swch = " << res_swch
729  // << "\n ias = " << ias
730  // << "\n CKM_A = " << CKM_A
731  // << "\n CKM_lambda = " << CKM_lambda
732  // << "\n CKM_barrho = " << CKM_barrho
733  // << "\n CKM_bareta = " << CKM_bareta
734  // << "\n CKM_factor = " << CKM_factor
735  // << "\n ============================================================================"
736  // << "\n Kinematics:"
737  // << "\n k_1 = " << k_1
738  // << "\n m_ell_1 =" << parent->getDaug(il1)->mass()
739  // << "\n k_2 = " << k_2
740  // << "\n m_ell_2 =" << parent->getDaug(il2)->mass()
741  // << "\n k_3 = " << k_3
742  // << "\n m_ell_3 =" << parent->getDaug(il3)->mass()
743  // << "\n k_4 = " << k_4
744  // << "\n m_ell_4 =" << parent->getDaug(il4)->mass()
745  // << "\n p = " << p
746  // << "\n q = " << q
747  // << "\n k = " << k
748  // << "\n ============================================================================"
749  // << "\n Form-factors"
750  // << "\n Fv = " << Fv
751  // << "\n Fa = " << Fa
752  // << "\n Ftv_b2q = " << Ftv_b2q
753  // << "\n Fta_b2q = " << Fta_b2q
754  // << "\n Ftv_barb2barq = " << Ftv_barb2barq
755  // << "\n Fta_barb2barq = " << Fta_barb2barq
756  // << "\n fb = " << fb
757  // << "\n ============================================================================"
758  // << "\n Wilson Coefficients:"
759  // << "\n Re(c7gam) = " << real(c7gam)
760  // << " Im(c7gam) = " << imag(c7gam)
761  // << "\n Re(c9eff_b2q) = " << real(c9eff_b2q)
762  // << " Im(c9eff_b2q) = " << imag(c9eff_b2q)
763  // << "\n Re(c9eff_barb2barq) = " << real(c9eff_barb2barq)
764  // << " Im(c9eff_barb2barq) = " << imag(c9eff_barb2barq)
765  // << "\n Re(c10a) = " << real(c10a)
766  // << " Im(c10a) = " << imag(c10a)
767  // << "\n ============================================================================"
768  // << "\n Functions in the matrix element:"
769  // << "\n a_barb2barq = " << a_barb2barq
770  // << "\n b_barb2barq = " << b_barb2barq
771  // << "\n c_barb2barq = " << c_barb2barq
772  // << "\n e_barb2barq = " << e_barb2barq
773  // << "\n f_barb2barq = " << f_barb2barq
774  // << "\n g_barb2barq = " << g_barb2barq
775  // << "\n ============================================================================"
776  // << "\n Partical Properties:"
777  // << "\n IdB = " << idparent << " == " << EvtPDL::getId(std::string("B_s0"))
778  // << "\n IdMu1 = " << IdMu1 << " == " << EvtPDL::getId(std::string("mu+"))
779  // << "\n IdMu2 = " << IdMu2 << " == " << EvtPDL::getId(std::string("mu-"))
780  // << "\n IdMu3 = " << IdMu3 << " == " << EvtPDL::getId(std::string("mu+"))
781  // << "\n IdMu4 = " << IdMu4 << " == " << EvtPDL::getId(std::string("mu-"))
782  // << "\n\n\n\n"
783  // << std::endl;
784  } else {
785  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
786  << "\n\n The function EvtbsToLLLLAmp::CalcAmp(...)"
787  << "\n Wrong Bq-meson number" << std::endl;
788  ::abort();
789  }
790  }
791 }
792 
793 //
794 // The decays Bq -> ell^+ ell^- ell^+ ell^- maximum probability calculation
795 //
796 double EvtbsToLLLLAmp::CalcMaxProb(
797  // EvtId parnum,
798  // EvtId l1num, EvtId l2num,
799  // EvtId l3num, EvtId l4num,
800  // Evtbs2llGammaFF *formFactors,
801  // EvtbTosllWilsCoeffNLO *WilsCoeff,
802  // double mu, int Nf,
803  // int res_swch, int ias,
804  // double CKM_A, double CKM_lambda,
805  // double CKM_barrho, double CKM_bareta
806 )
807 {
808  double maxfoundprob = 5.0; // maximum of the probability
809 
810  return maxfoundprob;
811 }
812 
813 // Triangular function
814 double EvtbsToLLLLAmp::lambda( double a, double b, double c )
815 {
816  double l;
817 
818  l = pow( a, 2.0 ) + pow( b, 2.0 ) + pow( c, 2.0 ) - 2.0 * a * b -
819  2.0 * a * c - 2.0 * b * c;
820 
821  return l;
822 }
EvtGenFunctions::directProd
EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
Definition: EvtTensor3C.cpp:154
EvtTensor4C::g
static const EvtTensor4C & g()
Definition: EvtTensor4C.cpp:44
EvtbTosllWilsCoeffNLO::GetC7Eff
EvtComplex GetC7Eff(double mu, double Mw, double mt, int Nf, int ias)
Definition: EvtbTosllWilsCoeffNLO.cpp:1230
dual
EvtTensor4C dual(const EvtTensor4C &t2)
Definition: EvtTensor4C.cpp:364
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cpp:172
Evtbs2llGammaFF::getPhotonFF
virtual void getPhotonFF(int, double, EvtId, double, double, double, double, EvtComplex, EvtComplex, EvtComplex, EvtComplex, EvtComplex &, EvtComplex &, EvtComplex &, EvtComplex &)
Definition: Evtbs2llGammaFF.hh:33
EvtTensor4C::cont2
EvtVector4C cont2(const EvtVector4C &v4) const
Definition: EvtTensor4C.cpp:461
EvtbsToLLLLAmp::lambda
double lambda(double a, double b, double c)
Definition: EvtbsToLLLLAmp.cpp:814
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cpp:130
a1
const double a1
Definition: EvtLi2Spence.cpp:23
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:167
EvtPDL::chg3
static int chg3(EvtId i)
Definition: EvtPDL.cpp:372
EvtParticle::spParent
virtual EvtDiracSpinor spParent(int) const
Definition: EvtParticle.cpp:665
EvtId
Definition: EvtId.hh:27
EvtLeptonACurrent
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:224
EvtParticle::getNDaug
size_t getNDaug() const
Definition: EvtParticle.cpp:160
EvtAmp::vertex
void vertex(const EvtComplex &amp)
Definition: EvtAmp.cpp:454
EvtbTosllWilsCoeffNLO::C2
double C2(double mu, double Mw, int Nf, int ias)
Definition: EvtbTosllWilsCoeffNLO.cpp:136
EvtVector4R::mass2
double mass2() const
Definition: EvtVector4R.hh:100
EvtbTosllWilsCoeffNLO::GetC9Eff
EvtComplex GetC9Eff(int decay_id, int res_swch, int ias, int Nf, double q2, double m2, double md, double mc, double mu, double mt, double Mw, double ml, double Relambda_qu, double Imlambda_qu)
Definition: EvtbTosllWilsCoeffNLO.cpp:1195
EvtAmp
Definition: EvtAmp.hh:30
EvtbTosllWilsCoeffNLO::GetC10Eff
EvtComplex GetC10Eff(double mt, double Mw)
Definition: EvtbTosllWilsCoeffNLO.cpp:1250
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:201
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cpp:150
imag
double imag(const EvtComplex &c)
Definition: EvtComplex.hh:235
EvtParticle::mass
double mass() const
Definition: EvtParticle.cpp:165
EvtbTosllWilsCoeffNLO::C1
double C1(double mu, double Mw, int Nf, int ias)
Definition: EvtbTosllWilsCoeffNLO.cpp:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtLeptonVCurrent
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
Definition: EvtDiracSpinor.cpp:209
EvtIdSet::contains
int contains(const EvtId id)
Definition: EvtIdSet.cpp:422
EvtbsToLLLLAmp::CalcAmp
void CalcAmp(EvtParticle *parent, EvtAmp &amp, Evtbs2llGammaFF *formFactors, EvtbTosllWilsCoeffNLO *WilsCoeff, double mu, int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta)
Definition: EvtbsToLLLLAmp.cpp:58
EvtParticle::getP4Restframe
EvtVector4R getP4Restframe() const
Definition: EvtParticle.cpp:793
real
double real(const EvtComplex &c)
Definition: EvtComplex.hh:230
Evtbs2llGammaFF::getQuarkMass
virtual double getQuarkMass(int)
Definition: Evtbs2llGammaFF.hh:44
Generated by   doxygen 1.8.15