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