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