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.
EvtbTosllVectorAmpNewExt.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 
22 
23 #include "EvtGenBase/EvtAmp.hh"
24 #include "EvtGenBase/EvtComplex.hh"
26 #include "EvtGenBase/EvtGenKine.hh"
27 #include "EvtGenBase/EvtId.hh"
28 #include "EvtGenBase/EvtIdSet.hh"
29 #include "EvtGenBase/EvtPDL.hh"
31 #include "EvtGenBase/EvtPatches.hh"
32 #include "EvtGenBase/EvtReport.hh"
37 
41 
42 #include <cstdlib>
43 
44 //
45 // The main functiom for the amplitude calculation
46 //
47 // input: *parent - the pointer to the parent particle (B-meson, the
48 // object of the EvtParticle class);
49 // *formFactors - the pointer to instance of EvtbTosllFFNew class object;
50 // *WilsCoeff - the pointer to the Standart Model Wilson Coefficients class;
51 // mu - the scale parameter, GeV;
52 // Nf - number of "effective" flavors (for b-quark Nf=5);
53 // res_swch - resonant switching parameter:
54 // = 0 the resonant contribution switched OFF,
55 // = 1 the resonant contribution switched ON;
56 // ias - switching parameter for \alpha_s(M_Z) value:
57 // = 0 PDG 1sigma minimal alpha_s(M_Z),
58 // = 1 PDG average value alpha_s(M_Z),
59 // = 2 PDG 1sigma maximal alpha_s(M_Z).
60 // Wolfenstein parameterization for CKM matrix
61 // CKM_A, CKM_lambda, CKM_barrho, CKM_bareta
62 // Multiplication coefficients
63 // A7 = ReA7 + i*ImA7
64 // A10 = ReA10 + i*ImA10
65 //
66 // return: amp - amplitude for the decay B -> V ell^+ ell^-
67 //
68 // Note: in our calculations we assume, that V-meson is the first
69 // daughter particle (iV=0) and leptons are the second and thirds
70 // daughter particles (il1=1 and il2=2).
71 //
73  EvtParticle* parent, EvtAmp& amp, EvtbTosllFFNew* formFactors,
74  EvtbTosllWilsCoeffNLO* WilsCoeff, double mu, int Nf, int res_swch, int ias,
75  double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta,
76  double ReA7, double ImA7, double ReA10, double ImA10 )
77 {
78  // FILE *mytest;
79 
80  EvtComplex unit1( 1.0, 0.0 ); // real unit
81  EvtComplex uniti( 0.0, 1.0 ); // imaginary unit
82 
83  EvtComplex A7 = ReA7 * unit1 + ImA7 * uniti;
84  EvtComplex A10 = ReA10 * unit1 + ImA10 * uniti;
85 
86  int iV = 0; // V-meson is the first daughter particle
87  int il1 = 1,
88  il2 = 2; // leptons are the second and thirds daughter particles
89 
90  // transition momentum of the leptonic pair q=k1+k2 or q=p1-p2
91  EvtVector4R q = parent->getDaug( il1 )->getP4() +
92  parent->getDaug( il2 )->getP4();
93 
94  // Mandelstam variable t=q^2
95  double q2 = q.mass2();
96 
97  double M1 = parent->mass(); // B - meson mass
98  double M2 = parent->getDaug( iV )->mass(); // V - meson mass
99  double ml = parent->getDaug( il1 )->mass(); // leptonic mass
100  double ms = 0.0; // light quark mass from the dispersion QM
101  double mc = formFactors->getQuarkMass( 4 ); // m_c mass from the dispersion QM
102  double mb = formFactors->getQuarkMass( 5 ); // m_b mass from the dispersion QM
103  // double Mw = EvtPDL::getNominalMass("W+"); // W-boson mass
104  // double mt = EvtPDL::getNominalMass("t"); // t-quark mass
105  double Mw = 80.403; // GeV W-boson mass
106  double mt = 174.2; // GeV t-quark mass
107 
108  EvtComplex Vtb, Vtq, Vub, Vuq; // V_{tb}, V_{tq}, V_{ub} and V_{uq}
109  EvtComplex CKM_factor; // V^*_{tq}*V_{tb}, where q={d,s}
110  EvtComplex lambda_qu; // V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}, where q={d,s}
111  double Relambda_qu, Imlambda_qu;
112 
113  EvtId idparent = parent->getId(); // B-meson Id
114  EvtId iddaught = parent->getDaug( iV )->getId(); // The vector meson Id
115 
116  // set of the light quark mass value
117  if ( ( idparent == EvtPDL::getId( std::string( "B+" ) ) &&
118  iddaught == EvtPDL::getId( std::string( "K*+" ) ) ) ||
119  ( idparent == EvtPDL::getId( std::string( "B-" ) ) &&
120  iddaught == EvtPDL::getId( std::string( "K*-" ) ) ) ||
121  ( idparent == EvtPDL::getId( std::string( "B0" ) ) &&
122  iddaught == EvtPDL::getId( std::string( "K*0" ) ) ) ||
123  ( idparent == EvtPDL::getId( std::string( "anti-B0" ) ) &&
124  iddaught == EvtPDL::getId( std::string( "anti-K*0" ) ) ) ||
125  ( idparent == EvtPDL::getId( std::string( "B_s0" ) ) &&
126  iddaught == EvtPDL::getId( std::string( "phi" ) ) ) ||
127  ( idparent == EvtPDL::getId( std::string( "anti-B_s0" ) ) &&
128  iddaught == EvtPDL::getId( std::string( "phi" ) ) ) ||
129  ( idparent == EvtPDL::getId( std::string( "B+" ) ) &&
130  iddaught == EvtPDL::getId( std::string( "K_1+" ) ) ) ||
131  ( idparent == EvtPDL::getId( std::string( "B-" ) ) &&
132  iddaught == EvtPDL::getId( std::string( "K_1-" ) ) ) ||
133  ( idparent == EvtPDL::getId( std::string( "B0" ) ) &&
134  iddaught == EvtPDL::getId( std::string( "K_10" ) ) ) ||
135  ( idparent == EvtPDL::getId( std::string( "anti-B0" ) ) &&
136  iddaught == EvtPDL::getId( std::string( "anti-K_10" ) ) ) ||
137  ( idparent == EvtPDL::getId( std::string( "B+" ) ) &&
138  iddaught == EvtPDL::getId( std::string( "K'_1+" ) ) ) ||
139  ( idparent == EvtPDL::getId( std::string( "B-" ) ) &&
140  iddaught == EvtPDL::getId( std::string( "K'_1-" ) ) ) ||
141  ( idparent == EvtPDL::getId( std::string( "B0" ) ) &&
142  iddaught == EvtPDL::getId( std::string( "K'_10" ) ) ) ||
143  ( idparent == EvtPDL::getId( std::string( "anti-B0" ) ) &&
144  iddaught == EvtPDL::getId( std::string( "anti-K'_10" ) ) ) ) {
145  ms = formFactors->getQuarkMass( 3 ); // m_s mass from the dispersion QM
146  // V_{ts}
147  Vtq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) +
148  pow( CKM_lambda, 2.0 ) *
149  ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
150  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
151  Vtq = -CKM_A * pow( CKM_lambda, 2.0 ) * Vtq;
152  // V_{us}
153  Vuq = CKM_lambda * unit1;
154  }
155 
156  if ( ( idparent == EvtPDL::getId( std::string( "B+" ) ) &&
157  iddaught == EvtPDL::getId( std::string( "rho+" ) ) ) ||
158  ( idparent == EvtPDL::getId( std::string( "B-" ) ) &&
159  iddaught == EvtPDL::getId( std::string( "rho-" ) ) ) ||
160  ( idparent == EvtPDL::getId( std::string( "B0" ) ) &&
161  iddaught == EvtPDL::getId( std::string( "rho0" ) ) ) ||
162  ( idparent == EvtPDL::getId( std::string( "anti-B0" ) ) &&
163  iddaught == EvtPDL::getId( std::string( "rho0" ) ) ) ||
164  ( idparent == EvtPDL::getId( std::string( "B0" ) ) &&
165  iddaught == EvtPDL::getId( std::string( "omega" ) ) ) ||
166  ( idparent == EvtPDL::getId( std::string( "anti-B0" ) ) &&
167  iddaught == EvtPDL::getId( std::string( "omega" ) ) ) ||
168  ( idparent == EvtPDL::getId( std::string( "B_s0" ) ) &&
169  iddaught == EvtPDL::getId( std::string( "anti-K*0" ) ) ) ||
170  ( idparent == EvtPDL::getId( std::string( "anti-B_s0" ) ) &&
171  iddaught == EvtPDL::getId( std::string( "K*0" ) ) ) ) {
172  ms = formFactors->getQuarkMass( 2 ); // m_d mass from the dispersion QM
173  // V_{td}
174  Vtq = unit1 - ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) ) *
175  ( CKM_barrho * unit1 + CKM_bareta * uniti ) /
176  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) );
177  Vtq = CKM_A * pow( CKM_lambda, 3.0 ) * Vtq;
178  // V_{ud}
179  Vuq = unit1 * ( 1.0 - 0.5 * pow( CKM_lambda, 2.0 ) -
180  0.125 * pow( CKM_lambda, 4.0 ) );
181  }
182 
183  if ( ms < 0.001 ) {
184  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
185  << "\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...)"
186  << "\n Error in the model set!"
187  << " ms = " << ms << std::endl;
188  ::abort();
189  }
190 
191  Vtb = unit1 * ( 1.0 - 0.5 * pow( CKM_A * CKM_lambda * CKM_lambda,
192  2.0 ) ); // V_{tb}
193  Vub = CKM_A * pow( CKM_lambda, 3.0 ) *
194  ( CKM_barrho * unit1 - CKM_bareta * uniti ) /
195  sqrt( 1.0 - pow( CKM_lambda, 2.0 ) ); // V_{ub}
196 
197  CKM_factor = conj( Vtq ) * Vtb; // V^*_{tq}*V_{tb}
198 
199  lambda_qu = conj( Vuq ) * Vub /
200  CKM_factor; // V^*_{uq}*V_{ub}/V^*_{tq}*V_{tb}
201  Relambda_qu = real( lambda_qu );
202  Imlambda_qu = imag( lambda_qu );
203 
204  double a1, a2, a0, v, t1, t2, t3; // B -> V transition form-factors
205 
206  // To get the B -> V transition form-factors
207  formFactors->getVectorFF( parent->getId(), parent->getDaug( iV )->getId(),
208  q2, a1, a2, a0, v, t1, t2, t3 );
209 
210  // The Wilson Coefficients preparation according to the paper
211  // A.J.Buras, M.Munz, Phys.Rev.D52, p.189 (1995)
212  EvtComplex c7gam = WilsCoeff->GetC7Eff( mu, Mw, mt, Nf, ias );
213  c7gam = c7gam * A7;
214  EvtComplex c9eff_b2q = WilsCoeff->GetC9Eff( 0, res_swch, ias, Nf, q2, mb,
215  ms, mc, mu, mt, Mw, ml,
216  Relambda_qu, Imlambda_qu );
217  EvtComplex c9eff_barb2barq = WilsCoeff->GetC9Eff( 1, res_swch, ias, Nf, q2,
218  mb, ms, mc, mu, mt, Mw, ml,
219  Relambda_qu, Imlambda_qu );
220  EvtComplex c10a = WilsCoeff->GetC10Eff( mt, Mw );
221  c10a = c10a * A10;
222 
223  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << "\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...) passed."
224  // << "\n Particle masses:"
225  // << "\n B - meson mass M1 = " << M1
226  // << "\n V - meson mass M2 = " << M2
227  // << "\n leptonic mass ml = " << ml
228  // << "\n light quark mass = " << ms
229  // << "\n c - quark mass mc = " << mc
230  // << "\n b - quark mass mb = " << mb
231  // << "\n t - quark mass mt = " << mt
232  // << "\n W - boson mass Mw = " << Mw
233  // << "\n ============================================================================"
234  // << "\n Input parameters:"
235  // << "\n scale parameter mu = " << mu
236  // << "\n number of flavors Nf = " << Nf
237  // << "\n resonant switching = " << res_swch
238  // << "\n parameter for alpha_s(M_Z) = " << ias
239  // << "\n ============================================================================"
240  // << "\n Vector form-factors at q^2 = " << q2
241  // << " for B -> V transition:"
242  // << "\n v = " << v
243  // << "\n a0 = " << a0
244  // << "\n a1 = " << a1
245  // << "\n a2 = " << a2
246  // << "\n t1 = " << t1
247  // << "\n t2 = " << t2
248  // << "\n t3 = " << t3
249  // << "\n ============================================================================"
250  // << "\n Wilson Coefficients:"
251  // << "\n Re(c7gam) = " << real(c7gam) << " Im(c7gam) = " << imag(c7gam)
252  // << "\n Re(c9eff_b2q) = " << real(c9eff_b2q)
253  // << " Im(c9eff_b2q) = " << imag(c9eff_b2q)
254  // << "\n Re(c9eff_barb2barq) = " << real(c9eff_barb2barq)
255  // << " Im(c9eff_barb2barq) = " << imag(c9eff_barb2barq)
256  // << "\n Re(c10a) = " << real(c10a) << " Im(c10a) = " << imag(c10a)
257  // << std::endl;
258 
259  // mytest = fopen("output.txt","a");
260  //
261  // if(mytest != NULL){
262  // fprintf(mytest,"%lf\n",q2);
263  // fclose(mytest);
264  // }
265  // else{
266  // EvtGenReport(EVTGEN_ERROR,"EvtGen") << "\n Error in writing to file.\n"
267  // << std::endl;
268  // return;
269  // }
270 
271  // 4- momentum of the B-meson in the the B-meson rest frame
272  EvtVector4R p1 = parent->getP4Restframe();
273  EvtVector4R hatp1 = p1 / M1;
274  // 4-momentum of the V-meson in the B-meson rest frame
275  EvtVector4R p2 = parent->getDaug( 0 )->getP4();
276  EvtVector4R hatp2 = p2 / M1;
277 
278  double hats = q2 / pow( M1, 2 );
279  double hatM2 = M2 / M1;
280  double hatmb = mb / M1;
281  double hatms = ms / M1;
282 
283  // Hadronic matrix element coefficients according to the paper
284  // A. Ali, A. Salim Safir, Eur.Phys.J.C25, pp.583-601 (2002)
285  // with m_s.NE.0
286  EvtComplex a_b2q, a_barb2barq, b_b2q, b_barb2barq, c_b2q, c_barb2barq, e, f,
287  g, h;
288 
289  a_b2q = 2.0 * c9eff_b2q * v / ( 1.0 + hatM2 ) +
290  4.0 * ( hatmb + hatms ) * c7gam * t1 / hats;
291  a_barb2barq = 2.0 * c9eff_barb2barq * v / ( 1.0 + hatM2 ) +
292  4.0 * ( hatmb + hatms ) * c7gam * t1 / hats;
293 
294  b_b2q = ( c9eff_b2q * a1 +
295  2.0 * ( hatmb - hatms ) * ( 1.0 - hatM2 ) * c7gam * t2 / hats ) *
296  ( 1.0 + hatM2 );
297  b_barb2barq = ( c9eff_barb2barq * a1 + 2.0 * ( hatmb - hatms ) *
298  ( 1.0 - hatM2 ) * c7gam * t2 /
299  hats ) *
300  ( 1.0 + hatM2 );
301 
302  c_b2q = ( c9eff_b2q * ( 1.0 - hatM2 ) * a2 +
303  2.0 * ( hatmb - hatms ) * ( 1.0 - pow( hatM2, 2 ) ) * c7gam * t2 /
304  hats +
305  2.0 * ( hatmb - hatms ) * c7gam * t3 ) /
306  ( 1 - pow( hatM2, 2 ) );
307 
308  c_barb2barq = ( c9eff_barb2barq * ( 1.0 - hatM2 ) * a2 +
309  2.0 * ( hatmb - hatms ) * ( 1.0 - pow( hatM2, 2 ) ) *
310  c7gam * t2 / hats +
311  2.0 * ( hatmb - hatms ) * c7gam * t3 ) /
312  ( 1 - pow( hatM2, 2 ) );
313 
314  e = 2.0 * c10a * v / ( 1 + hatM2 );
315 
316  f = ( 1.0 + hatM2 ) * c10a * a1;
317 
318  g = c10a * a2 / ( 1 + hatM2 );
319 
320  h = ( ( 1.0 + hatM2 ) * a1 - ( 1.0 - hatM2 ) * a2 - 2.0 * hatM2 * a0 ) *
321  c10a / hats;
322 
323  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << " a_b2q = " << a_b2q
324  // << " a_barb2barq = " << a_barb2barq
325  // << " b_b2q = " << b_b2q
326  // << " b_barb2barq = " << b_barb2barq
327  // << " c_b2q = " << c_b2q
328  // << " c_barb2barq = " << c_barb2barq
329  // << " e = " << e
330  // << " f = " << f
331  // << " g = " << g
332  // << " h = " << h
333  // << std::endl;
334 
335  // to find ell^+ and ell^- in the B-meson daughters
336  int charge1 = EvtPDL::chg3( parent->getDaug( 1 )->getId() );
337  int charge2 = EvtPDL::chg3( parent->getDaug( 2 )->getId() );
338 
339  EvtParticle* lepPlus = 0;
340  EvtParticle* lepMinus = 0;
341 
342  lepPlus = ( charge1 > charge2 ) ? parent->getDaug( 1 ) : parent->getDaug( 2 );
343  lepMinus = ( charge1 < charge2 ) ? parent->getDaug( 1 )
344  : parent->getDaug( 2 );
345 
346  EvtTensor4C T1, T2; // hadronic matrix element tensor structures
347 
348  EvtVector4C epsV; // vector meson polarisation vector
349 
350  int i; // vector meson polarisations counter
351 
352  EvtVector4C lvc11, lvc12; // spin structures for
353  EvtVector4C lvc21, lvc22; // the leptonic vector current
354 
355  EvtVector4C lac11, lac12; // spin structures for
356  EvtVector4C lac21, lac22; // the leptonic axial current
357 
358  // B - and barB - mesons descriptors
359  EvtIdSet bmesons( "B-", "anti-B0", "anti-B_s0", "B_c-" );
360  EvtIdSet bbarmesons( "B+", "B0", "B_s0", "B_c+" );
361 
362  EvtId parentID = parent->getId();
363 
364  if ( bmesons.contains( parentID ) ) {
365  // The amplitude for the decay barB -> barV ell^+ ell^-
366  // (b -> q ell^+ ell^- transition)
367 
368  T1 = -a_b2q * unit1 * dual( EvtGenFunctions::directProd( hatp1, hatp2 ) ) -
369  b_b2q * uniti * EvtTensor4C::g() +
370  c_b2q * uniti *
371  EvtGenFunctions::directProd( ( hatp1 + hatp2 ), hatp1 );
372 
373  T2 = -e * unit1 * dual( EvtGenFunctions::directProd( hatp1, hatp2 ) ) -
374  f * uniti * EvtTensor4C::g() +
375  g * uniti * EvtGenFunctions::directProd( ( hatp1 + hatp2 ), hatp1 ) +
376  h * uniti * EvtGenFunctions::directProd( ( hatp1 - hatp2 ), hatp1 );
377 
378  lvc11 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
379  lepMinus->spParent( 0 ) );
380  lvc21 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
381  lepMinus->spParent( 0 ) );
382  lvc12 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
383  lepMinus->spParent( 1 ) );
384  lvc22 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
385  lepMinus->spParent( 1 ) );
386 
387  lac11 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
388  lepMinus->spParent( 0 ) );
389  lac21 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
390  lepMinus->spParent( 0 ) );
391  lac12 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
392  lepMinus->spParent( 1 ) );
393  lac22 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
394  lepMinus->spParent( 1 ) );
395 
396  // summing up vector meson polarisations \epsilon^*_{\nu}(i)
397  for ( i = 0; i < 3; i++ ) {
398  EvtVector4C epsV = parent->getDaug( 0 )->epsParent( i ).conj();
399 
400  EvtVector4C E1 = M1 * T1.cont2( epsV );
401  EvtVector4C E2 = M1 * T2.cont2( epsV );
402 
403  amp.vertex( i, 0, 0, CKM_factor * ( lvc11 * E1 + lac11 * E2 ) );
404  amp.vertex( i, 0, 1, CKM_factor * ( lvc12 * E1 + lac12 * E2 ) );
405  amp.vertex( i, 1, 0, CKM_factor * ( lvc21 * E1 + lac21 * E2 ) );
406  amp.vertex( i, 1, 1, CKM_factor * ( lvc22 * E1 + lac22 * E2 ) );
407  }
408 
409  } else {
410  if ( bbarmesons.contains( parentID ) ) {
411  // The amplitude for the decay B -> V ell^+ ell^-
412  // (barb -> barq ell^+ ell^- transition)
413 
414  T1 = a_barb2barq * unit1 *
415  dual( EvtGenFunctions::directProd( hatp1, hatp2 ) ) -
416  b_barb2barq * uniti * EvtTensor4C::g() +
417  c_barb2barq * uniti *
418  EvtGenFunctions::directProd( ( hatp1 + hatp2 ), hatp1 );
419 
420  T2 = e * unit1 * dual( EvtGenFunctions::directProd( hatp1, hatp2 ) ) -
421  f * uniti * EvtTensor4C::g() +
422  g * uniti *
423  EvtGenFunctions::directProd( ( hatp1 + hatp2 ), hatp1 ) +
424  h * uniti *
425  EvtGenFunctions::directProd( ( hatp1 - hatp2 ), hatp1 );
426 
427  lvc11 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
428  lepMinus->spParent( 1 ) );
429  lvc21 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
430  lepMinus->spParent( 1 ) );
431  lvc12 = EvtLeptonVCurrent( lepPlus->spParent( 1 ),
432  lepMinus->spParent( 0 ) );
433  lvc22 = EvtLeptonVCurrent( lepPlus->spParent( 0 ),
434  lepMinus->spParent( 0 ) );
435 
436  lac11 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
437  lepMinus->spParent( 1 ) );
438  lac21 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
439  lepMinus->spParent( 1 ) );
440  lac12 = EvtLeptonACurrent( lepPlus->spParent( 1 ),
441  lepMinus->spParent( 0 ) );
442  lac22 = EvtLeptonACurrent( lepPlus->spParent( 0 ),
443  lepMinus->spParent( 0 ) );
444 
445  // summing up vector meson polarisations \epsilon^*_{\nu}(i)
446  for ( i = 0; i < 3; i++ ) {
447  EvtVector4C barepsV = parent->getDaug( 0 )->epsParent( i ).conj();
448 
449  EvtVector4C E3 = M1 * T1.cont2( barepsV );
450  EvtVector4C E4 = M1 * T2.cont2( barepsV );
451 
452  amp.vertex( i, 0, 0,
453  conj( CKM_factor ) * ( lvc11 * E3 + lac11 * E4 ) );
454  amp.vertex( i, 0, 1,
455  conj( CKM_factor ) * ( lvc12 * E3 + lac12 * E4 ) );
456  amp.vertex( i, 1, 0,
457  conj( CKM_factor ) * ( lvc21 * E3 + lac21 * E4 ) );
458  amp.vertex( i, 1, 1,
459  conj( CKM_factor ) * ( lvc22 * E3 + lac22 * E4 ) );
460  }
461 
462  } else {
463  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
464  << "\n\n The function EvtbTosllVectorAmpNew::CalcAmp(...)"
465  << "\n Wrong B-meson number" << std::endl;
466  ::abort();
467  }
468  }
469 
470  // Test of the signature for Levi-Civita tensor
471  // EvtVector4C Vec0, Vec1, Vec2, Vec3;
472  // EvtTensor4C Ttest;
473  // Vec0.set(1.0,0.0,0.0,0.0);
474  // Vec1.set(0.0,1.0,0.0,0.0);
475  // Vec2.set(0.0,0.0,1.0,0.0);
476  // Vec3.set(0.0,0.0,0.0,1.0);
477  // Ttest=dual(directProd(Vec2,Vec3));
478  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << "\n\n\n e^{0123} =" << Ttest.get(0,1) << std::endl;
479  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << " e^{1023} =" << Ttest.get(1,0) << std::endl;
480  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << " e^{1123} =" << Ttest.get(1,1) << "\n" << std::endl;
481  // EvtVector4C Vtest=Ttest.cont2(Vec1);
482  // for(i=0;i<=3;i++){
483  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << " Vtest =" << Vtest.get(i) << std::endl;
484  // }
485  // EvtComplex Atest;
486  // Atest=Vec0*Vtest;
487  // EvtGenReport(EVTGEN_NOTICE,"EvtGen") << "\n Atest =" << Atest << "\n\n\n" << std::endl;
488 }
489 
490 //
491 // The decays B -> V ell^+ ell^- maximum probability calculation for the
492 // d^2\Gamma/dq^2 d\cos\theta distribution.
493 //
494 // \theta - the angle between the vector meson and ell^- directions in the
495 // B-meson rest frame.
496 //
497 // If ias=0 (nonresonant case), the maximum is achieved at q2=q2_min=4*ml^2.
498 // If ias=1 (resonant case), the maximum is achieved at q2=M^2_{J/\psi}.
499 //
501  EvtId parnum, EvtId mesnum, EvtId l1num, EvtId l2num,
502  EvtbTosllFFNew* formFactors, EvtbTosllWilsCoeffNLO* WilsCoeff, double mu,
503  int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda,
504  double CKM_barrho, double CKM_bareta, double ReA7, double ImA7,
505  double ReA10, double ImA10 )
506 {
507  double maxfoundprob = -100.0; // maximum of the probability
508  int katmax = 0;
509 
510  double M1 = EvtPDL::getMeanMass( parnum ); // B - meson mass
511  double M2 = EvtPDL::getMeanMass( mesnum ); // V - meson mass
512  double ml = EvtPDL::getMeanMass( l1num ); // leptonic mass
513 
514  if ( res_swch == 0 ) {
515  // B-meson rest frame particles and they kinematics inicialization
516  double s_min, t_for_s;
517  s_min = 4.0 * pow( ml, 2.0 );
518  t_for_s = 0.5 *
519  ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - 2.0 * pow( ml, 2.0 ) );
520 
521  double EV, El2;
522  EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s_min ) /
523  ( 2.0 * M1 ); // V-meson energy
524  El2 = ( s_min + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
525  ( 2.0 * M1 ); // ell^- energy
526 
527  double modV, modl2;
528  modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
529  modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
530 
531  double cosVellminus; // angle between the vector meson and ell^- directions
532  cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 -
533  t_for_s ) /
534  ( 2.0 * modV * modl2 );
535  if ( ( fabs( cosVellminus ) > 1.0 ) &&
536  ( fabs( cosVellminus ) <= 1.0001 ) ) {
537  // EvtGenReport(EVTGEN_DEBUG,"EvtGen")
538  // << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
539  // << "\n cos(theta) = " << cosVellminus
540  // << std::endl;
541  cosVellminus = cosVellminus / fabs( cosVellminus );
542  }
543  if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
544  cosVellminus = cosVellminus / fabs( cosVellminus );
545  EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
546  << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
547  << "\n modV = " << modV << "\n modl2 = " << modl2
548  << "\n cos(theta) = " << cosVellminus
549  << "\n t_for_s = " << t_for_s << "\n s_min = " << s_min
550  << "\n EV = " << EV << "\n El2 = " << El2
551  << "\n M2 = " << M2 << "\n ml = " << ml
552  << std::endl;
553  }
554  if ( fabs( cosVellminus ) > 1.0001 ) {
555  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
556  << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
557  << "\n |cos(theta)| = " << fabs( cosVellminus ) << " > 1"
558  << "\n s_min = " << s_min << "\n t_for_s = " << t_for_s
559  << "\n EV = " << EV << "\n El2 = " << El2
560  << "\n modV = " << modV << "\n modl2 = " << modl2
561  << "\n M2 = " << M2 << "\n ml = " << ml << std::endl;
562  ::abort();
563  }
564 
565  EvtVector4R p1, p2, k1, k2;
566  p1.set( M1, 0.0, 0.0, 0.0 );
567  p2.set( EV, modV, 0.0, 0.0 );
568  k2.set( El2, modl2 * cosVellminus,
569  -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 );
570  k1 = p1 - p2 - k2;
571 
572  // EvtGenReport(EVTGEN_DEBUG,"EvtGen")
573  // << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
574  // << "\n mu =" << mu << " Nf =" << Nf
575  // << " res_swch =" << res_swch
576  // << " ias =" << ias
577  // << "\n CKM_A = " << CKM_A
578  // << " CKM_lambda = " << CKM_lambda
579  // << "\n CKM_barrho = " << CKM_barrho
580  // << " CKM_bareta = " << CKM_bareta
581  // << "\n M1 = " << M1
582  // << "\n M2 = " << M2
583  // << "\n ml = " << ml
584  // << "\n s_min = " << s_min
585  // << "\n t_for_s = " << t_for_s
586  // << "\n EV = " << EV
587  // << "\n El1 = " << El1
588  // << "\n El2 = " << El2
589  // << "\n modV = " << modV
590  // << "\n modl1 = " << modl1
591  // << "\n modl2 = " << modl2
592  // << "\n cos(theta) = " << cosVellminus
593  // << "\n p1 =" << p1
594  // << "\n p2 =" << p2
595  // << "\n k1 =" << k1
596  // << "\n k2 =" << k2
597  // << std::endl;
598 
599  // B-meson state preparation at the rest frame of B-meson
600  EvtScalarParticle* scalar_part;
601  EvtParticle* root_part;
602  scalar_part = new EvtScalarParticle;
603 
604  scalar_part->noLifeTime();
605  scalar_part->init( parnum, p1 );
606  root_part = (EvtParticle*)scalar_part;
607  root_part->setDiagonalSpinDensity();
608 
609  // Amplitude initialization
610  EvtId listdaug[3];
611  listdaug[0] = mesnum;
612  listdaug[1] = l1num;
613  listdaug[2] = l2num;
614 
615  EvtAmp amp;
616  amp.init( parnum, 3, listdaug );
617 
618  // Daughters states preparation at the rest frame of B-meson
619  root_part->makeDaughters( 3, listdaug );
620 
621  EvtParticle *vect, *lep1, *lep2;
622  vect = root_part->getDaug( 0 );
623  lep1 = root_part->getDaug( 1 );
624  lep2 = root_part->getDaug( 2 );
625 
626  vect->noLifeTime();
627  lep1->noLifeTime();
628  lep2->noLifeTime();
629 
630  vect->init( mesnum, p2 );
631  lep1->init( l1num, k1 );
632  lep2->init( l2num, k2 );
633 
634  EvtSpinDensity rho;
635  rho.setDiag( root_part->getSpinStates() );
636 
637  // The amplitude calculation at the
638  // "maximum amplitude" kinematical configuration
639  CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch, ias,
640  CKM_A, CKM_lambda, CKM_barrho, CKM_bareta, ReA7, ImA7, ReA10,
641  ImA10 );
642 
643  // Now find the probability at this q2 and cos theta lepton point
644  maxfoundprob = rho.normalizedProb( amp.getSpinDensity() );
645 
646  delete scalar_part;
647  // delete root_part;
648  delete vect;
649  delete lep1;
650  delete lep2;
651 
652  } // if(res_swch==0)
653 
654  if ( res_swch == 1 ) {
655  double s, t_for_s; // Mandelstam variables
656  double t_plus, t_minus; // t-variable boundaries for current s-variable
657  double dt;
658  int k;
659 
660  s = pow( 3.09688, 2.0 ); // s = (M_{J/\psi})^2
661 
662  t_plus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
663  t_plus = t_plus + sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
664  sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
665  t_plus *= 0.5;
666 
667  t_minus = pow( M1, 2.0 ) + pow( M2, 2.0 ) + 2.0 * pow( ml, 2.0 ) - s;
668  t_minus = t_minus -
669  sqrt( 1.0 - 4.0 * pow( ml, 2.0 ) / s ) *
670  sqrt( lambda( s, pow( M1, 2.0 ), pow( M2, 2.0 ) ) );
671  t_minus *= 0.5;
672 
673  dt = ( t_plus - t_minus ) / 1000.0;
674 
675  // The maximum probability calculation
676  for ( k = 0; k <= 1000; k++ ) {
677  t_for_s = t_plus - dt * ( (double)k );
678 
679  if ( ( t_for_s < t_minus ) && ( t_for_s >= ( 0.9999 * t_minus ) ) ) {
680  t_for_s = t_minus;
681  }
682  if ( t_for_s < ( 0.9999 * t_minus ) ) {
683  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
684  << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
685  << "\n t_for_s = " << t_for_s << " < t_minus = " << t_minus
686  << " ! "
687  << "\n t_plus = " << t_plus << "\n dt = " << dt
688  << "\n k = " << k << "\n s = " << s
689  << "\n M1 = " << M1 << "\n M2 = " << M2
690  << "\n ml = " << ml << std::endl;
691  ::abort();
692  }
693 
694  // B-meson rest frame particles and they kinematics inicialization
695  double EV, El2;
696  EV = ( pow( M1, 2.0 ) + pow( M2, 2.0 ) - s ) /
697  ( 2.0 * M1 ); // V-meson energy
698  El2 = ( s + t_for_s - pow( M2, 2.0 ) - pow( ml, 2.0 ) ) /
699  ( 2.0 * M1 ); // ell^- energy
700 
701  double modV, modl2;
702  modV = sqrt( pow( EV, 2.0 ) - pow( M2, 2.0 ) );
703  modl2 = sqrt( pow( El2, 2.0 ) - pow( ml, 2.0 ) );
704 
705  double cosVellminus; // angle between the vector meson and ell^- directions
706  cosVellminus = ( pow( M2, 2.0 ) + pow( ml, 2.0 ) + 2.0 * EV * El2 -
707  t_for_s ) /
708  ( 2.0 * modV * modl2 );
709  if ( ( fabs( cosVellminus ) > 1.0 ) &&
710  ( fabs( cosVellminus ) <= 1.0001 ) ) {
711  // EvtGenReport(EVTGEN_DEBUG,"EvtGen")
712  // << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
713  // << "\n cos(theta) = " << cosVellminus
714  // << std::endl;
715  cosVellminus = cosVellminus / fabs( cosVellminus );
716  }
717  if ( ( modV <= 0.000001 ) || ( modl2 <= 0.000001 ) ) {
718  cosVellminus = cosVellminus / fabs( cosVellminus );
719  EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
720  << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
721  << "\n modV = " << modV << "\n modl2 = " << modl2
722  << "\n cos(theta) = " << cosVellminus
723  << "\n s = " << s << "\n t_for_s = " << t_for_s
724  << "\n t_plus = " << t_plus
725  << "\n t_minus = " << t_minus << "\n dt = " << dt
726  << "\n EV = " << EV << "\n El2 = " << El2
727  << "\n M2 = " << M2 << "\n ml = " << ml
728  << std::endl;
729  }
730  if ( fabs( cosVellminus ) > 1.0001 ) {
731  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
732  << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
733  << "\n |cos(theta)| = " << fabs( cosVellminus ) << " > 1"
734  << "\n s = " << s << "\n t_for_s = " << t_for_s
735  << "\n EV = " << EV << "\n El2 = " << El2
736  << "\n modV = " << modV << "\n modl2 = " << modl2
737  << "\n M2 = " << M2 << "\n ml = " << ml
738  << std::endl;
739  ::abort();
740  }
741 
742  EvtVector4R p1, p2, k1, k2;
743  p1.set( M1, 0.0, 0.0, 0.0 );
744  p2.set( EV, modV, 0.0, 0.0 );
745  k2.set( El2, modl2 * cosVellminus,
746  -modl2 * sqrt( 1.0 - pow( cosVellminus, 2.0 ) ), 0.0 );
747  k1 = p1 - p2 - k2;
748 
749  // EvtGenReport(EVTGEN_DEBUG,"EvtGen")
750  // << "\n Debug in the function EvtbTosllVectorAmpNew::CalcMaxProb(...):"
751  // << "\n mu =" << mu << " Nf =" << Nf
752  // << " res_swch =" << res_swch
753  // << " ias =" << ias
754  // << "\n CKM_A = " << CKM_A
755  // << " CKM_lambda = " << CKM_lambda
756  // << "\n CKM_barrho = " << CKM_barrho
757  // << " CKM_bareta = " << CKM_bareta
758  // << "\n M1 = " << M1
759  // << "\n M2 = " << M2
760  // << "\n ml = " << ml
761  // << "\n s = " << s
762  // << "\n t_for_s = " << t_for_s
763  // << "\n EV = " << EV
764  // << "\n El1 = " << El1
765  // << "\n El2 = " << El2
766  // << "\n modV = " << modV
767  // << "\n modl1 = " << modl1
768  // << "\n modl2 = " << modl2
769  // << "\n cos(theta) = " << cosVellminus
770  // << "\n p1 =" << p1
771  // << "\n p2 =" << p2
772  // << "\n k1 =" << k1
773  // << "\n k2 =" << k2
774  // << std::endl;
775 
776  // B-meson state preparation at the rest frame of B-meson
777  EvtScalarParticle* scalar_part;
778  EvtParticle* root_part;
779  scalar_part = new EvtScalarParticle;
780 
781  scalar_part->noLifeTime();
782  scalar_part->init( parnum, p1 );
783  root_part = (EvtParticle*)scalar_part;
784  root_part->setDiagonalSpinDensity();
785 
786  // Amplitude initialization
787  EvtId listdaug[3];
788  listdaug[0] = mesnum;
789  listdaug[1] = l1num;
790  listdaug[2] = l2num;
791 
792  EvtAmp amp;
793  amp.init( parnum, 3, listdaug );
794 
795  // Daughters states preparation at the rest frame of B-meson
796  root_part->makeDaughters( 3, listdaug );
797 
798  EvtParticle *vect, *lep1, *lep2;
799  vect = root_part->getDaug( 0 );
800  lep1 = root_part->getDaug( 1 );
801  lep2 = root_part->getDaug( 2 );
802 
803  vect->noLifeTime();
804  lep1->noLifeTime();
805  lep2->noLifeTime();
806 
807  vect->init( mesnum, p2 );
808  lep1->init( l1num, k1 );
809  lep2->init( l2num, k2 );
810 
811  EvtSpinDensity rho;
812  rho.setDiag( root_part->getSpinStates() );
813 
814  // The amplitude calculation at the
815  // "maximum amplitude" kinematical configuration
816  CalcAmp( root_part, amp, formFactors, WilsCoeff, mu, Nf, res_swch,
817  ias, CKM_A, CKM_lambda, CKM_barrho, CKM_bareta, ReA7, ImA7,
818  ReA10, ImA10 );
819 
820  // Now find the probability at this q2 and cos theta lepton point
821  double nikmax = rho.normalizedProb( amp.getSpinDensity() );
822 
823  if ( nikmax > maxfoundprob ) {
824  maxfoundprob = nikmax;
825  katmax = k;
826  EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
827  << "\n maxfoundprob ( s =" << s << ", t = " << t_for_s
828  << " ) = " << maxfoundprob << "\n k =" << katmax
829  << std::endl;
830  }
831 
832  delete scalar_part;
833  // delete root_part;
834  delete vect;
835  delete lep1;
836  delete lep2;
837 
838  } // for(k=0; k<=1000; k++)
839 
840  } // if(res_swch==1)
841 
842  if ( maxfoundprob <= 0.0 ) {
843  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
844  << "\n\n In the function EvtbTosllVectorAmpNew::CalcMaxProb(...)"
845  << "\n maxfoundprob = " << maxfoundprob << " <0 or =0!"
846  << "\n res_swch = " << res_swch << std::endl;
847  ::abort();
848  }
849 
850  EvtGenReport( EVTGEN_NOTICE, "EvtGen" )
851  << "\n maxfoundprob (...) = " << maxfoundprob << std::endl;
852 
853  maxfoundprob *= 1.01;
854 
855  // EvtGenReport(EVTGEN_NOTICE,"EvtGen")
856  // << "\n ***************************************************************************"
857  // << "\n The function EvtbTosllVectorAmpNew::CalcMaxProb(...) passed with arguments:"
858  // << "\n mu =" << mu << " Nf =" << Nf
859  // << " res_swch =" << res_swch
860  // << " ias =" << ias
861  // << "\n CKM_A = " << CKM_A
862  // << " CKM_lambda = " << CKM_lambda
863  // << "\n CKM_barrho = " << CKM_barrho
864  // << " CKM_bareta = " << CKM_bareta
865  // << "\n The distribution maximum maxfoundprob =" << maxfoundprob
866  // << "\n k = " << katmax
867  // << "\n ***************************************************************************"
868  // << std::endl;
869 
870  return maxfoundprob;
871 }
872 
873 // Triangular function
874 double EvtbTosllVectorAmpNewExt::lambda( double a, double b, double c )
875 {
876  double l;
877 
878  l = pow( a, 2.0 ) + pow( b, 2.0 ) + pow( c, 2.0 ) - 2.0 * a * b -
879  2.0 * a * c - 2.0 * b * c;
880 
881  return l;
882 }
EvtTensor3C directProd(const EvtVector3C &c1, const EvtVector3C &c2)
double CalcMaxProb(EvtId parnum, EvtId mesnum, EvtId l1num, EvtId l2num, EvtbTosllFFNew *formFactors, EvtbTosllWilsCoeffNLO *WilsCoeff, double mu, int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta, double ReA7, double ImA7, double ReA10, double ImA10) override
void setDiag(int n)
static const EvtTensor4C & g()
Definition: EvtTensor4C.cpp:44
EvtComplex GetC7Eff(double mu, double Mw, double mt, int Nf, int ias)
EvtTensor4C dual(const EvtTensor4C &t2)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cpp:172
const double a2
EvtVector4C cont2(const EvtVector4C &v4) const
static double getMeanMass(EvtId i)
Definition: EvtPDL.cpp:314
int getSpinStates() const
EvtId getId() const
void makeDaughters(unsigned int ndaug, EvtId *id)
const double a1
void set(int i, double d)
Definition: EvtVector4R.hh:167
EvtSpinDensity getSpinDensity()
Definition: EvtAmp.cpp:142
double lambda(double a, double b, double c) override
static int chg3(EvtId i)
Definition: EvtPDL.cpp:372
virtual EvtDiracSpinor spParent(int) const
Definition: EvtId.hh:27
EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void vertex(const EvtComplex &amp)
Definition: EvtAmp.cpp:454
virtual double getQuarkMass(int)
double mass2() const
Definition: EvtVector4R.hh:100
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)
virtual EvtVector4C epsParent(int i) const
Definition: EvtAmp.hh:30
EvtComplex GetC10Eff(double mt, double Mw)
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
const EvtVector4R & getP4() const
double imag(const EvtComplex &c)
Definition: EvtComplex.hh:235
void init(EvtId part_n, double e, double px, double py, double pz)
double mass() const
virtual void getVectorFF(EvtId, EvtId, double, double &, double &, double &, double &, double &, double &, double &)
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor &d, const EvtDiracSpinor &dp)
void noLifeTime()
Definition: EvtParticle.hh:382
EvtVector4C conj() const
Definition: EvtVector4C.hh:202
int contains(const EvtId id)
Definition: EvtIdSet.cpp:422
EvtVector4R getP4Restframe() const
void init(EvtId p, int ndaug, EvtId *daug)
Definition: EvtAmp.cpp:68
double real(const EvtComplex &c)
Definition: EvtComplex.hh:230
double normalizedProb(const EvtSpinDensity &d)
void CalcAmp(EvtParticle *parent, EvtAmp &amp, EvtbTosllFFNew *formFactors, EvtbTosllWilsCoeffNLO *WilsCoeff, double mu, int Nf, int res_swch, int ias, double CKM_A, double CKM_lambda, double CKM_barrho, double CKM_bareta, double ReA7, double ImA7, double ReA10, double ImA10) override