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