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