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
EvtBtoXsgammaKagan.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
25#include "EvtGenBase/EvtPDL.hh"
29
39
40#include <fstream>
41#include <stdlib.h>
42#include <string>
43using std::endl;
44using std::fstream;
45
48
49void EvtBtoXsgammaKagan::init( int nArg, double* args )
50{
51 if ( ( nArg ) > 12 || ( nArg > 1 && nArg < 10 ) || nArg == 11 ) {
52 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
53 << "EvtBtoXsgamma generator model "
54 << "EvtBtoXsgammaKagan expected "
55 << "either 1(default config) or "
56 << "10 (default mass range) or "
57 << "12 (user range) arguments but found: " << nArg << endl;
58 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
59 << "Will terminate execution!" << endl;
60 ::abort();
61 }
62
63 if ( nArg == 1 ) {
64 m_bbprod = true;
66 } else {
67 m_bbprod = false;
68 computeHadronicMass( nArg, args );
69 }
70
71 double mHminLimit = 0.6373;
72 double mHmaxLimit = 4.5;
73
74 if ( nArg > 10 ) {
75 m_mHmin = args[10];
76 m_mHmax = args[11];
77 if ( m_mHmin > m_mHmax ) {
78 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
79 << "Minimum hadronic mass exceeds maximum " << endl;
80 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
81 << "Will terminate execution!" << endl;
82 ::abort();
83 }
84 if ( m_mHmin < mHminLimit ) {
85 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
86 << "Minimum hadronic mass below K pi threshold" << endl;
87 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
88 << "Resetting to K pi threshold" << endl;
89 m_mHmin = mHminLimit;
90 }
91 if ( m_mHmax > mHmaxLimit ) {
92 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
93 << "Maximum hadronic mass above 4.5 GeV/c^2" << endl;
94 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
95 << "Resetting to 4.5 GeV/c^2" << endl;
96 m_mHmax = mHmaxLimit;
97 }
98 } else {
99 m_mHmin = mHminLimit; // usually just above K pi threshold for Xsd/u
100 m_mHmax = mHmaxLimit;
101 }
102}
103
105{
106 m_massHad = { 0, 0.0625995, 0.125199, 0.187798, 0.250398, 0.312997,
107 0.375597, 0.438196, 0.500796, 0.563395, 0.625995, 0.688594,
108 0.751194, 0.813793, 0.876392, 0.938992, 1.00159, 1.06419,
109 1.12679, 1.18939, 1.25199, 1.31459, 1.37719, 1.43979,
110 1.50239, 1.56499, 1.62759, 1.69019, 1.75278, 1.81538,
111 1.87798, 1.94058, 2.00318, 2.06578, 2.12838, 2.19098,
112 2.25358, 2.31618, 2.37878, 2.44138, 2.50398, 2.56658,
113 2.62918, 2.69178, 2.75438, 2.81698, 2.87958, 2.94217,
114 3.00477, 3.06737, 3.12997, 3.19257, 3.25517, 3.31777,
115 3.38037, 3.44297, 3.50557, 3.56817, 3.63077, 3.69337,
116 3.75597, 3.81857, 3.88117, 3.94377, 4.00637, 4.06896,
117 4.13156, 4.19416, 4.25676, 4.31936, 4.38196, 4.44456,
118 4.50716, 4.56976, 4.63236, 4.69496, 4.75756, 4.82016,
119 4.88276, 4.94536, 5.00796 };
120 m_brHad = { 0, 1.03244e-09, 3.0239e-08, 1.99815e-07, 7.29392e-07,
121 1.93129e-06, 4.17806e-06, 7.86021e-06, 1.33421e-05, 2.09196e-05,
122 3.07815e-05, 4.29854e-05, 5.74406e-05, 7.3906e-05, 9.2003e-05,
123 0.000111223, 0.000130977, 0.000150618, 0.000169483, 0.000186934,
124 0.000202392, 0.000215366, 0.000225491, 0.000232496, 0.000236274,
125 0.000236835, 0.000234313, 0.000228942, 0.000221042, 0.000210994,
126 0.000199215, 0.000186137, 0.000172194, 0.000157775, 0.000143255,
127 0.000128952, 0.000115133, 0.000102012, 8.97451e-05, 7.84384e-05,
128 6.81519e-05, 5.89048e-05, 5.06851e-05, 4.34515e-05, 3.71506e-05,
129 3.1702e-05, 2.70124e-05, 2.30588e-05, 1.96951e-05, 1.68596e-05,
130 1.44909e-05, 1.25102e-05, 1.08596e-05, 9.48476e-06, 8.34013e-06,
131 7.38477e-06, 6.58627e-06, 5.91541e-06, 5.35022e-06, 4.87047e-06,
132 4.46249e-06, 4.11032e-06, 3.80543e-06, 3.54051e-06, 3.30967e-06,
133 3.10848e-06, 2.93254e-06, 2.78369e-06, 2.65823e-06, 2.55747e-06,
134 2.51068e-06, 2.57179e-06, 2.74684e-06, 3.02719e-06, 3.41182e-06,
135 3.91387e-06, 4.56248e-06, 5.40862e-06, 6.53915e-06, 8.10867e-06,
136 1.04167e-05 };
137 m_massHad.resize( 81 );
138 m_brHad.resize( 81 );
139
140 m_intervalMH = 80;
141}
142
143void EvtBtoXsgammaKagan::computeHadronicMass( int /*nArg*/, double* args )
144{
145 //Input parameters
146 int fermiFunction = (int)args[1];
147 m_mB = args[2];
148 m_mb = args[3];
149 m_mu = args[4];
150 m_lam1 = args[5];
151 m_delta = args[6];
152 m_z = args[7];
153 m_nIntervalS = args[8];
154 m_nIntervalmH = args[9];
155 std::vector<double> mHVect( int( m_nIntervalmH + 1.0 ) );
156 m_massHad.clear();
157 m_massHad.resize( int( m_nIntervalmH + 1.0 ) );
158 m_brHad.clear();
159 m_brHad.resize( int( m_nIntervalmH + 1.0 ) );
161
162 //Going to have to add a new entry into the data file - takes ages...
163 EvtGenReport( EVTGEN_WARNING, "EvtGen" )
164 << "EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..."
165 << endl;
166
167 //Now need to compute the mHVect vector for
168 //the current parameters
169
170 //A few more parameters
171 double _mubar = m_mu;
172 m_mW = 80.33;
173 m_mt = 175.0;
174 m_alpha = 1. / 137.036;
176 m_kappabar = 3.382 - 4.14 * ( sqrt( m_z ) - 0.29 );
177 m_fz = Fz( m_z );
178 m_rer8 = ( 44. / 9. ) - ( 8. / 27. ) * pow( EvtConst::pi, 2. );
179 m_r7 = ( -10. / 3. ) - ( 8. / 9. ) * pow( EvtConst::pi, 2. );
180 m_rer2 = -4.092 + 12.78 * ( sqrt( m_z ) - .29 );
181 m_gam77 = 32. / 3.;
182 m_gam27 = 416. / 81.;
183 m_gam87 = -32. / 9.;
184 m_lam2 = .12;
185 m_beta0 = 23. / 3.;
186 m_beta1 = 116. / 3.;
187 m_alphasmZ = .118;
188 m_mZ = 91.187;
189 m_ms = m_mb / 50.;
190
191 double eGammaMin = 0.5 * m_mB * ( 1. - m_delta );
192 double eGammaMax = 0.5 * m_mB;
193 double yMin = 2. * eGammaMin / m_mB;
194 double yMax = 2. * eGammaMax / m_mB;
195 double _CKMrat = 0.976;
196 double Nsl = 1.0;
197
198 //Calculate alpha the various scales
202 m_alphasmubar = CalcAlphaS( _mubar );
203
204 //Calculate the Wilson Coefficients and Delta
206 m_kSLemmu = ( 12. / 23. ) * ( ( 1. / m_etamu ) - 1. );
208 CalcDelta();
209
210 //Build s22 and s27 vector - saves time because double
211 //integration is required otherwise
212 std::vector<double> s22Coeffs( int( m_nIntervalS + 1.0 ) );
213 std::vector<double> s27Coeffs( int( m_nIntervalS + 1.0 ) );
214 std::vector<double> s28Coeffs( int( m_nIntervalS + 1.0 ) );
215
216 double dy = ( yMax - yMin ) / m_nIntervalS;
217 double yp = yMin;
218
219 std::vector<double> sCoeffs( 1 );
220 sCoeffs[0] = m_z;
221
222 //Define s22 and s27 functions
223 auto mys22Func = EvtItgPtrFunction{ &s22Func, 0., yMax + 0.1, sCoeffs };
224 auto mys27Func = EvtItgPtrFunction{ &s27Func, 0., yMax + 0.1, sCoeffs };
225
226 //Use a simpson integrator
227 auto mys22Simp = EvtItgSimpsonIntegrator{ mys22Func, 1.0e-4, 20 };
228 auto mys27Simp = EvtItgSimpsonIntegrator{ mys27Func, 1.0e-4, 50 };
229
230 int i;
231
232 for ( i = 0; i < int( m_nIntervalS + 1.0 ); i++ ) {
233 s22Coeffs[i] = ( 16. / 27. ) * mys22Simp.evaluate( 1.0e-20, yp );
234 s27Coeffs[i] = ( -8. / 9. ) * m_z * mys27Simp.evaluate( 1.0e-20, yp );
235 s28Coeffs[i] = -s27Coeffs[i] / 3.;
236 yp = yp + dy;
237 }
238
239 //Define functions and vectors used to calculate mHVect. Each function takes a set
240 //of vectors which are used as the function coefficients
241 std::vector<double> FermiCoeffs( 6 );
242 std::vector<double> varCoeffs( 3 );
243 std::vector<double> DeltaCoeffs( 1 );
244 std::vector<double> s88Coeffs( 2 );
245 std::vector<double> sInitCoeffs( 3 );
246
247 varCoeffs[0] = m_mB;
248 varCoeffs[1] = m_mb;
249 varCoeffs[2] = 0.;
250
251 DeltaCoeffs[0] = m_alphasmu;
252
253 s88Coeffs[0] = m_mb;
254 s88Coeffs[1] = m_ms;
255
256 sInitCoeffs[0] = m_nIntervalS;
257 sInitCoeffs[1] = yMin;
258 sInitCoeffs[2] = yMax;
259
260 FermiCoeffs[0] = fermiFunction;
261 FermiCoeffs[1] = 0.0;
262 FermiCoeffs[2] = 0.0;
263 FermiCoeffs[3] = 0.0;
264 FermiCoeffs[4] = 0.0;
265 FermiCoeffs[5] = 0.0;
266
267 //Coefficients for gamma function
268 std::vector<double> gammaCoeffs( 6 );
269 gammaCoeffs[0] = 76.18009172947146;
270 gammaCoeffs[1] = -86.50532032941677;
271 gammaCoeffs[2] = 24.01409824083091;
272 gammaCoeffs[3] = -1.231739572450155;
273 gammaCoeffs[4] = 0.1208650973866179e-2;
274 gammaCoeffs[5] = -0.5395239384953e-5;
275
276 //Calculate quantities for the fermi function to be used
277 //Distinguish among the different shape functions
278 if ( fermiFunction == 1 ) {
279 FermiCoeffs[1] = m_lambdabar;
280 FermiCoeffs[2] = ( -3. * pow( m_lambdabar, 2. ) / m_lam1 ) - 1.;
281 FermiCoeffs[3] = m_lam1;
282 FermiCoeffs[4] = 1.0;
283
284 auto myNormFunc = std::make_unique<EvtItgPtrFunction>(
286 FermiCoeffs );
287 auto myNormSimp =
288 std::make_unique<EvtItgSimpsonIntegrator>( *myNormFunc, 1.0e-4, 40 );
289 FermiCoeffs[4] = myNormSimp->normalisation();
290
291 } else if ( fermiFunction == 2 ) {
293 m_lam1, m_mb,
294 gammaCoeffs );
295 FermiCoeffs[1] = m_lambdabar;
296 FermiCoeffs[2] = a;
297 FermiCoeffs[3] =
298 EvtBtoXsgammaFermiUtil::Gamma( ( 2.0 + a ) / 2., gammaCoeffs ) /
299 EvtBtoXsgammaFermiUtil::Gamma( ( 1.0 + a ) / 2., gammaCoeffs );
300 FermiCoeffs[4] = 1.0;
301
302 auto myNormFunc = std::make_unique<EvtItgPtrFunction>(
304 FermiCoeffs );
305 auto myNormSimp =
306 std::make_unique<EvtItgSimpsonIntegrator>( *myNormFunc, 1.0e-4, 40 );
307 FermiCoeffs[4] = myNormSimp->normalisation();
308
309 } else if ( fermiFunction == 3 ) {
311 m_lam1 );
312 FermiCoeffs[1] = m_mB;
313 FermiCoeffs[2] = m_mb;
314 FermiCoeffs[3] = rho;
315 FermiCoeffs[4] = m_lambdabar;
316 FermiCoeffs[5] = 1.0;
317
318 auto myNormFunc = std::make_unique<EvtItgPtrFunction>(
320 FermiCoeffs );
321 auto myNormSimp =
322 std::make_unique<EvtItgSimpsonIntegrator>( *myNormFunc, 1.0e-4, 40 );
323 FermiCoeffs[5] = myNormSimp->normalisation();
324 }
325
326 //Define functions
327 auto myDeltaFermiFunc = EvtItgThreeCoeffFcn{ &DeltaFermiFunc, -m_mb,
328 m_mB - m_mb, FermiCoeffs,
329 varCoeffs, DeltaCoeffs };
330 auto mys88FermiFunc = EvtItgThreeCoeffFcn{ &s88FermiFunc, -m_mb,
331 m_mB - m_mb, FermiCoeffs,
332 varCoeffs, s88Coeffs };
333 auto mys77FermiFunc = EvtItgTwoCoeffFcn{ &s77FermiFunc, -m_mb, m_mB - m_mb,
334 FermiCoeffs, varCoeffs };
335 auto mys78FermiFunc = EvtItgTwoCoeffFcn{ &s78FermiFunc, -m_mb, m_mB - m_mb,
336 FermiCoeffs, varCoeffs };
337 auto mys22FermiFunc = EvtItgFourCoeffFcn{ &sFermiFunc, -m_mb,
338 m_mB - m_mb, FermiCoeffs,
339 varCoeffs, sInitCoeffs,
340 s22Coeffs };
341 auto mys27FermiFunc = EvtItgFourCoeffFcn{ &sFermiFunc, -m_mb,
342 m_mB - m_mb, FermiCoeffs,
343 varCoeffs, sInitCoeffs,
344 s27Coeffs };
345 auto mys28FermiFunc = EvtItgFourCoeffFcn{ &sFermiFunc, -m_mb,
346 m_mB - m_mb, FermiCoeffs,
347 varCoeffs, sInitCoeffs,
348 s28Coeffs };
349
350 //Define integrators
351 auto myDeltaFermiSimp = EvtItgSimpsonIntegrator{ myDeltaFermiFunc, 1.0e-4,
352 40 };
353 auto mys77FermiSimp = EvtItgSimpsonIntegrator{ mys77FermiFunc, 1.0e-4, 40 };
354 auto mys88FermiSimp = EvtItgSimpsonIntegrator{ mys88FermiFunc, 1.0e-4, 40 };
355 auto mys78FermiSimp = EvtItgSimpsonIntegrator{ mys78FermiFunc, 1.0e-4, 40 };
356 auto mys22FermiSimp = EvtItgSimpsonIntegrator{ mys22FermiFunc, 1.0e-4, 40 };
357 auto mys27FermiSimp = EvtItgSimpsonIntegrator{ mys27FermiFunc, 1.0e-4, 40 };
358 auto mys28FermiSimp = EvtItgSimpsonIntegrator{ mys28FermiFunc, 1.0e-4, 40 };
359
360 //Finally calculate mHVect for the range of hadronic masses
361 double mHmin = sqrt( m_mB * m_mB - 2. * m_mB * eGammaMax );
362 double mHmax = sqrt( m_mB * m_mB - 2. * m_mB * eGammaMin );
363 double dmH = ( mHmax - mHmin ) / m_nIntervalmH;
364
365 double mH = mHmin;
366
367 //Calculating the Branching Fractions
368 for ( i = 0; i < int( m_nIntervalmH + 1.0 ); i++ ) {
369 double ymH = 1. - ( ( mH * mH ) / ( m_mB * m_mB ) );
370
371 //Need to set ymH as one of the input parameters
372 myDeltaFermiFunc.setCoeff( 2, 2, ymH );
373 mys77FermiFunc.setCoeff( 2, 2, ymH );
374 mys88FermiFunc.setCoeff( 2, 2, ymH );
375 mys78FermiFunc.setCoeff( 2, 2, ymH );
376 mys22FermiFunc.setCoeff( 2, 2, ymH );
377 mys27FermiFunc.setCoeff( 2, 2, ymH );
378 mys28FermiFunc.setCoeff( 2, 2, ymH );
379
380 //Integrate
381
382 double deltaResult = myDeltaFermiSimp.evaluate( ( m_mB * ymH - m_mb ),
383 m_mB - m_mb );
384 double s77Result = mys77FermiSimp.evaluate( ( m_mB * ymH - m_mb ),
385 m_mB - m_mb );
386 double s88Result = mys88FermiSimp.evaluate( ( m_mB * ymH - m_mb ),
387 m_mB - m_mb );
388 double s78Result = mys78FermiSimp.evaluate( ( m_mB * ymH - m_mb ),
389 m_mB - m_mb );
390 double s22Result = mys22FermiSimp.evaluate( ( m_mB * ymH - m_mb ),
391 m_mB - m_mb );
392 double s27Result = mys27FermiSimp.evaluate( ( m_mB * ymH - m_mb ),
393 m_mB - m_mb );
394 mys28FermiSimp.evaluate( ( m_mB * ymH - m_mb ), m_mB - m_mb );
395
396 double py =
397 ( pow( _CKMrat, 2. ) * ( 6. / m_fz ) * ( m_alpha / EvtConst::pi ) *
398 ( deltaResult * m_cDeltatot +
400 ( s77Result * pow( m_c70mu, 2. ) +
401 s27Result * m_c2mu * ( m_c70mu - m_c80mu / 3. ) +
402 s78Result * m_c70mu * m_c80mu + s22Result * m_c2mu * m_c2mu +
403 s88Result * m_c80mu * m_c80mu ) ) );
404
405 mHVect[i] = 2. * ( mH / ( m_mB * m_mB ) ) * 0.105 * Nsl * py;
406
407 m_massHad[i] = mH;
408 m_brHad[i] = 2. * ( mH / ( m_mB * m_mB ) ) * 0.105 * Nsl * py;
409
410 mH = mH + dmH;
411 }
412}
413
414double EvtBtoXsgammaKagan::GetMass( int /*Xscode*/ )
415{
416 // Get hadronic mass for the event according to the hadronic mass spectra computed in computeHadronicMass
417 double mass = 0.0;
418 double min = m_mHmin;
419 if ( m_bbprod )
420 min = 1.1;
421 // double max=4.5;
422 double max = m_mHmax;
423 double xbox( 0 ), ybox( 0 );
424 double boxheight( 0 );
425 double trueHeight( 0 );
426 double boxwidth = max - min;
427 double wgt( 0. );
428
429 for ( int i = 0; i < int( m_intervalMH + 1.0 ); i++ ) {
430 if ( m_brHad[i] > boxheight )
431 boxheight = m_brHad[i];
432 }
433 while ( ( mass > max ) || ( mass < min ) ) {
434 xbox = EvtRandom::Flat( boxwidth ) + min;
435 ybox = EvtRandom::Flat( boxheight );
436 trueHeight = 0.0;
437 // Correction by Peter Richardson
438 for ( int i = 1; i < int( m_intervalMH + 1.0 ); ++i ) {
439 if ( ( m_massHad[i] >= xbox ) && ( 0.0 == trueHeight ) ) {
440 wgt = ( xbox - m_massHad[i - 1] ) /
441 ( m_massHad[i] - m_massHad[i - 1] );
442 trueHeight = m_brHad[i - 1] +
443 wgt * ( m_brHad[i] - m_brHad[i - 1] );
444 }
445 }
446
447 if ( ybox > trueHeight ) {
448 mass = 0.0;
449 } else {
450 mass = xbox;
451 }
452 }
453
454 return mass;
455}
456
457double EvtBtoXsgammaKagan::CalcAlphaS( double scale )
458{
459 double v = 1. - m_beta0 * ( m_alphasmZ / ( 2. * EvtConst::pi ) ) *
460 ( log( m_mZ / scale ) );
461 return ( m_alphasmZ / v ) *
462 ( 1. - ( ( m_beta1 / m_beta0 ) *
463 ( m_alphasmZ / ( 4. * EvtConst::pi ) ) * ( log( v ) / v ) ) );
464}
465
467{
468 double mtatmw = m_mt * pow( ( m_alphasmW / m_alphasmt ), ( 12. / 23. ) ) *
469 ( 1 +
470 ( 12. / 23. ) * ( ( 253. / 18. ) - ( 116. / 23. ) ) *
471 ( ( m_alphasmW - m_alphasmt ) / ( 4.0 * EvtConst::pi ) ) -
472 ( 4. / 3. ) * ( m_alphasmt / EvtConst::pi ) );
473 double xt = pow( mtatmw, 2. ) / pow( m_mW, 2. );
474
476 m_c2mu = .5 * pow( m_etamu, ( -12. / 23. ) ) +
477 .5 * pow( m_etamu, ( 6. / 23. ) );
478
479 double c7mWsm = ( ( 3. * pow( xt, 3. ) - 2. * pow( xt, 2. ) ) /
480 ( 4. * pow( ( xt - 1. ), 4. ) ) ) *
481 log( xt ) +
482 ( ( -8. * pow( xt, 3. ) - 5. * pow( xt, 2. ) + 7. * xt ) /
483 ( 24. * pow( ( xt - 1. ), 3. ) ) );
484
485 double c8mWsm = ( ( -3. * pow( xt, 2. ) ) / ( 4. * pow( ( xt - 1. ), 4. ) ) ) *
486 log( xt ) +
487 ( ( -pow( xt, 3. ) + 5. * pow( xt, 2. ) + 2. * xt ) /
488 ( 8. * pow( ( xt - 1. ), 3. ) ) );
489
490 double c7constmu = ( 626126. / 272277. ) * pow( m_etamu, ( 14. / 23. ) ) -
491 ( 56281. / 51730. ) * pow( m_etamu, ( 16. / 23. ) ) -
492 ( 3. / 7. ) * pow( m_etamu, ( 6. / 23. ) ) -
493 ( 1. / 14. ) * pow( m_etamu, ( -12. / 23. ) ) -
494 .6494 * pow( m_etamu, .4086 ) -
495 .038 * pow( m_etamu, -.423 ) -
496 .0186 * pow( m_etamu, -.8994 ) -
497 .0057 * pow( m_etamu, .1456 );
498
499 m_c70mu = c7mWsm * pow( m_etamu, ( 16. / 23. ) ) +
500 ( 8. / 3. ) *
501 ( pow( m_etamu, ( 14. / 23. ) ) -
502 pow( m_etamu, ( 16. / 23. ) ) ) *
503 c8mWsm +
504 c7constmu;
505
506 double c8constmu = ( 313063. / 363036. ) * pow( m_etamu, ( 14. / 23. ) ) -
507 .9135 * pow( m_etamu, .4086 ) +
508 .0873 * pow( m_etamu, -.423 ) -
509 .0571 * pow( m_etamu, -.8994 ) +
510 .0209 * pow( m_etamu, .1456 );
511
512 m_c80mu = c8mWsm * pow( m_etamu, ( 14. / 23. ) ) + c8constmu;
513
514 //Compute the dilogarithm (PolyLog(2,x)) with the Simpson integrator
515 //The dilogarithm is defined as: Li_2(x)=Int_0^x(-log(1.-z)/z)
516 //however, Mathematica implements it as Sum[z^k/k^2,{k,1,Infinity}], so, althought the two
517 //results are similar and both implemented in the program, we prefer to use the
518 //one closer to the Mathematica implementation as identical to what used by the theorists.
519
520 // EvtItgFunction *myDiLogFunc = new EvtItgFunction(&diLogFunc, 0., 1.-1./xt);
521 //EvtItgAbsIntegrator *myDiLogSimp = new EvtItgSimpsonIntegrator(*myDiLogFunc, 1.0e-4, 50);
522 //double li2 = myDiLogSimp->evaluate(1.0e-20,1.-1./xt);
523
524 double li2 = diLogMathematica( 1. - 1. / xt );
525
526 double c7mWsm1 =
527 ( ( -16. * pow( xt, 4. ) - 122. * pow( xt, 3. ) + 80. * pow( xt, 2. ) -
528 8. * xt ) /
529 ( 9. * pow( ( xt - 1. ), 4. ) ) * li2 +
530 ( 6. * pow( xt, 4. ) + 46. * pow( xt, 3. ) - 28. * pow( xt, 2. ) ) /
531 ( 3. * pow( ( xt - 1. ), 5. ) ) * pow( log( xt ), 2. ) +
532 ( -102. * pow( xt, 5. ) - 588. * pow( xt, 4. ) -
533 2262. * pow( xt, 3. ) + 3244. * pow( xt, 2. ) - 1364. * xt + 208. ) /
534 ( 81. * pow( ( xt - 1 ), 5. ) ) * log( xt ) +
535 ( 1646. * pow( xt, 4. ) + 12205. * pow( xt, 3. ) -
536 10740. * pow( xt, 2. ) + 2509. * xt - 436. ) /
537 ( 486. * pow( ( xt - 1 ), 4. ) ) );
538
539 double c8mWsm1 =
540 ( ( -4. * pow( xt, 4. ) + 40. * pow( xt, 3. ) + 41. * pow( xt, 2. ) + xt ) /
541 ( 6. * pow( ( xt - 1. ), 4. ) ) * li2 +
542 ( -17. * pow( xt, 3. ) - 31. * pow( xt, 2. ) ) /
543 ( 2. * pow( ( xt - 1. ), 5. ) ) * pow( log( xt ), 2. ) +
544 ( -210. * pow( xt, 5. ) + 1086. * pow( xt, 4. ) +
545 4893. * pow( xt, 3. ) + 2857. * pow( xt, 2. ) - 1994. * xt + 280. ) /
546 ( 216. * pow( ( xt - 1 ), 5. ) ) * log( xt ) +
547 ( 737. * pow( xt, 4. ) - 14102. * pow( xt, 3. ) -
548 28209. * pow( xt, 2. ) + 610. * xt - 508. ) /
549 ( 1296. * pow( ( xt - 1 ), 4. ) ) );
550
551 double E1 = ( xt * ( 18. - 11. * xt - pow( xt, 2. ) ) /
552 ( 12. * pow( ( 1. - xt ), 3. ) ) +
553 pow( xt, 2. ) * ( 15. - 16. * xt + 4. * pow( xt, 2. ) ) /
554 ( 6. * pow( ( 1. - xt ), 4. ) ) * log( xt ) -
555 2. / 3. * log( xt ) );
556
557 double e1 = 4661194. / 816831.;
558 double e2 = -8516. / 2217.;
559 double e3 = 0.;
560 double e4 = 0.;
561 double e5 = -1.9043;
562 double e6 = -.1008;
563 double e7 = .1216;
564 double e8 = .0183;
565
566 double f1 = -17.3023;
567 double f2 = 8.5027;
568 double f3 = 4.5508;
569 double f4 = .7519;
570 double f5 = 2.004;
571 double f6 = .7476;
572 double f7 = -.5385;
573 double f8 = .0914;
574
575 double g1 = 14.8088;
576 double g2 = -10.809;
577 double g3 = -.874;
578 double g4 = .4218;
579 double g5 = -2.9347;
580 double g6 = .3971;
581 double g7 = .1600;
582 double g8 = .0225;
583
584 double c71constmu =
585 ( ( e1 * m_etamu * E1 + f1 + g1 * m_etamu ) *
586 pow( m_etamu, ( 14. / 23. ) ) +
587 ( e2 * m_etamu * E1 + f2 + g2 * m_etamu ) *
588 pow( m_etamu, ( 16. / 23. ) ) +
589 ( e3 * m_etamu * E1 + f3 + g3 * m_etamu ) * pow( m_etamu, ( 6. / 23. ) ) +
590 ( e4 * m_etamu * E1 + f4 + g4 * m_etamu ) *
591 pow( m_etamu, ( -12. / 23. ) ) +
592 ( e5 * m_etamu * E1 + f5 + g5 * m_etamu ) * pow( m_etamu, .4086 ) +
593 ( e6 * m_etamu * E1 + f6 + g6 * m_etamu ) * pow( m_etamu, ( -.423 ) ) +
594 ( e7 * m_etamu * E1 + f7 + g7 * m_etamu ) * pow( m_etamu, ( -.8994 ) ) +
595 ( e8 * m_etamu * E1 + f8 + g8 * m_etamu ) * pow( m_etamu, .1456 ) );
596
597 double c71pmu = ( ( ( 297664. / 14283. * pow( m_etamu, ( 16. / 23. ) ) -
598 7164416. / 357075. * pow( m_etamu, ( 14. / 23. ) ) +
599 256868. / 14283. * pow( m_etamu, ( 37. / 23. ) ) -
600 6698884. / 357075. * pow( m_etamu, ( 39. / 23. ) ) ) *
601 ( c8mWsm ) ) +
602 37208. / 4761. *
603 ( pow( m_etamu, ( 39. / 23. ) ) -
604 pow( m_etamu, ( 16. / 23. ) ) ) *
605 ( c7mWsm ) +
606 c71constmu );
607
609 ( pow( m_etamu, ( 16. / 23. ) ) * c7mWsm1 +
610 8. / 3. *
611 ( pow( m_etamu, ( 14. / 23. ) ) -
612 pow( m_etamu, ( 16. / 23. ) ) ) *
613 c8mWsm1 ) +
614 c71pmu );
615
616 m_c7emmu = ( ( 32. / 75. * pow( m_etamu, ( -9. / 23. ) ) -
617 40. / 69. * pow( m_etamu, ( -7. / 23. ) ) +
618 88. / 575. * pow( m_etamu, ( 16. / 23. ) ) ) *
619 c7mWsm +
620 ( -32. / 575. * pow( m_etamu, ( -9. / 23. ) ) +
621 32. / 1449. * pow( m_etamu, ( -7. / 23. ) ) +
622 640. / 1449. * pow( m_etamu, ( 14. / 23. ) ) -
623 704. / 1725. * pow( m_etamu, ( 16. / 23. ) ) ) *
624 c8mWsm -
625 190. / 8073. * pow( m_etamu, ( -35. / 23. ) ) -
626 359. / 3105. * pow( m_etamu, ( -17. / 23. ) ) +
627 4276. / 121095. * pow( m_etamu, ( -12. / 23. ) ) +
628 350531. / 1009125. * pow( m_etamu, ( -9. / 23. ) ) +
629 2. / 4347. * pow( m_etamu, ( -7. / 23. ) ) -
630 5956. / 15525. * pow( m_etamu, ( 6. / 23. ) ) +
631 38380. / 169533. * pow( m_etamu, ( 14. / 23. ) ) -
632 748. / 8625. * pow( m_etamu, ( 16. / 23. ) ) );
633
634 // Wilson coefficients values as according to Kagan's program
635 // m_c2mu=1.10566;
636 //m_c70mu=-0.314292;
637 // m_c80mu=-0.148954;
638 // m_c71mu=0.480964;
639 // m_c7emmu=0.0323219;
640}
641
643{
644 double cDelta77 = ( 1. +
645 ( m_alphasmu / ( 2. * EvtConst::pi ) ) *
646 ( m_r7 - ( 16. / 3. ) + m_gam77 * log( m_mb / m_mu ) ) +
647 ( ( pow( ( 1. - m_z ), 4. ) / m_fz ) - 1. ) *
648 ( 6. * m_lam2 / pow( m_mb, 2. ) ) +
649 ( m_alphasmubar / ( 2. * EvtConst::pi ) ) * m_kappabar ) *
650 pow( m_c70mu, 2. );
651
652 double cDelta27 = ( ( m_alphasmu / ( 2. * EvtConst::pi ) ) *
653 ( m_rer2 + m_gam27 * log( m_mb / m_mu ) ) -
654 ( m_lam2 / ( 9. * m_z * pow( m_mb, 2. ) ) ) ) *
655 m_c2mu * m_c70mu;
656
657 double cDelta78 = ( m_alphasmu / ( 2. * EvtConst::pi ) ) *
658 ( m_rer8 + m_gam87 * log( m_mb / m_mu ) ) * m_c70mu *
659 m_c80mu;
660
661 m_cDeltatot = cDelta77 + cDelta27 + cDelta78 +
662 ( m_alphasmu / ( 2. * EvtConst::pi ) ) * m_c71mu * m_c70mu +
663 ( m_alpha / m_alphasmu ) * ( 2. * m_c7emmu * m_c70mu -
664 m_kSLemmu * pow( m_c70mu, 2. ) );
665}
666
667double EvtBtoXsgammaKagan::Delta( double y, double alphasMu )
668{
669 //Fix for singularity at endpoint
670 if ( y >= 1.0 )
671 y = 0.9999999999;
672
673 return ( -4. * ( alphasMu / ( 3. * EvtConst::pi * ( 1. - y ) ) ) *
674 ( log( 1. - y ) + 7. / 4. ) *
675 exp( -2. * ( alphasMu / ( 3. * EvtConst::pi ) ) *
676 ( pow( log( 1. - y ), 2 ) + ( 7. / 2. ) * log( 1. - y ) ) ) );
677}
678
679double EvtBtoXsgammaKagan::s77( double y )
680{
681 //Fix for singularity at endpoint
682 if ( y >= 1.0 )
683 y = 0.9999999999;
684
685 return ( ( 1. / 3. ) *
686 ( 7. + y - 2. * pow( y, 2 ) - 2. * ( 1. + y ) * log( 1. - y ) ) );
687}
688
689double EvtBtoXsgammaKagan::s88( double y, double mb, double ms )
690{
691 //Fix for singularity at endpoint
692 if ( y >= 1.0 )
693 y = 0.9999999999;
694
695 return ( ( 1. / 27. ) * ( ( 2. * ( 2. - 2. * y + pow( y, 2 ) ) / y ) *
696 ( log( 1. - y ) + 2. * log( mb / ms ) ) -
697 2. * pow( y, 2 ) - y - 8. * ( ( 1. - y ) / y ) ) );
698}
699
700double EvtBtoXsgammaKagan::s78( double y )
701{
702 //Fix for singularity at endpoint
703 if ( y >= 1.0 )
704 y = 0.9999999999;
705
706 return ( ( 8. / 9. ) * ( ( ( 1. - y ) / y ) * log( 1. - y ) + 1. +
707 ( pow( y, 2 ) / 4. ) ) );
708}
709
710double EvtBtoXsgammaKagan::ReG( double y )
711{
712 if ( y < 4. )
713 return -2. * pow( atan( sqrt( y / ( 4. - y ) ) ), 2. );
714 else {
715 return 2. * ( pow( log( ( sqrt( y ) + sqrt( y - 4. ) ) / 2. ), 2. ) ) -
716 ( 1. / 2. ) * pow( EvtConst::pi, 2. );
717 }
718}
719
720double EvtBtoXsgammaKagan::ImG( double y )
721{
722 if ( y < 4. )
723 return 0.0;
724 else {
725 return ( -2. * EvtConst::pi * log( ( sqrt( y ) + sqrt( y - 4. ) ) / 2. ) );
726 }
727}
728
729double EvtBtoXsgammaKagan::s22Func( double y, const std::vector<double>& coeffs )
730{
731 //coeffs[0]=z
732 return ( 1. - y ) *
733 ( ( pow( coeffs[0], 2. ) / pow( y, 2. ) ) *
734 ( pow( ReG( y / coeffs[0] ), 2. ) +
735 pow( ImG( y / coeffs[0] ), 2. ) ) +
736 ( coeffs[0] / y ) * ReG( y / coeffs[0] ) + ( 1. / 4. ) );
737}
738
739double EvtBtoXsgammaKagan::s27Func( double y, const std::vector<double>& coeffs )
740{
741 //coeffs[0] = z
742 return ( ReG( y / coeffs[0] ) + y / ( 2. * coeffs[0] ) );
743}
744
746 const std::vector<double>& coeffs1,
747 const std::vector<double>& coeffs2,
748 const std::vector<double>& coeffs3 )
749{
750 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
751 //coeffs2[2]=ymH, coeffs3[0]=DeltaCoeff (alphasmu)
752
753 return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
754 Delta( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ), coeffs3[0] );
755}
756
758 const std::vector<double>& coeffs1,
759 const std::vector<double>& coeffs2 )
760{
761 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
762 //coeffs2[2]=ymH
763 return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
764 s77( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ) );
765}
766
768 const std::vector<double>& coeffs1,
769 const std::vector<double>& coeffs2,
770 const std::vector<double>& coeffs3 )
771{
772 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
773 //coeffs2[2]=ymH, coeffs3=s88 coeffs
774 return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
775 s88( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ), coeffs3[0],
776 coeffs3[1] );
777}
778
780 const std::vector<double>& coeffs1,
781 const std::vector<double>& coeffs2 )
782{
783 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
784 //coeffs2[2]=ymH
785 return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
786 s78( ( coeffs2[0] * coeffs2[2] ) / ( coeffs2[1] + y ) );
787}
788
790 const std::vector<double>& coeffs1,
791 const std::vector<double>& coeffs2,
792 const std::vector<double>& coeffs3,
793 const std::vector<double>& coeffs4 )
794{
795 //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
796 //coeffs2[2]=ymH, coeffs3[0]=nIntervals in s22 or s27 array, coeffs3[1]=yMin,
797 //coeffs3[2]=yMax, coeffs4=s22 or s27 array
798 return FermiFunc( y, coeffs1 ) * ( coeffs2[0] / ( coeffs2[1] + y ) ) *
799 GetArrayVal( coeffs2[0] * coeffs2[2] / ( coeffs2[1] + y ),
800 coeffs3[0], coeffs3[1], coeffs3[2], coeffs4 );
801}
802
803double EvtBtoXsgammaKagan::Fz( double z )
804{
805 return ( 1. - 8. * z + 8. * pow( z, 3. ) - pow( z, 4. ) -
806 12. * pow( z, 2. ) * log( z ) );
807}
808
809double EvtBtoXsgammaKagan::GetArrayVal( double xp, double nInterval, double xMin,
810 double xMax, std::vector<double> array )
811{
812 double dx = ( xMax - xMin ) / nInterval;
813 int bin1 = int( ( ( xp - xMin ) / ( xMax - xMin ) ) * nInterval );
814
815 double x1 = double( bin1 ) * dx + xMin;
816
817 if ( xp == x1 )
818 return array[bin1];
819
820 int bin2( 0 );
821 if ( xp > x1 ) {
822 bin2 = bin1 + 1;
823 } else if ( xp < x1 ) {
824 bin2 = bin1 - 1;
825 }
826
827 if ( bin1 <= 0 ) {
828 bin1 = 0;
829 bin2 = 1;
830 }
831
832 //If xp is in the last bin, always interpolate between the last two bins
833 if ( bin1 == (int)nInterval ) {
834 bin2 = (int)nInterval;
835 bin1 = (int)nInterval - 1;
836 x1 = double( bin1 ) * dx + xMin;
837 }
838
839 double x2 = double( bin2 ) * dx + xMin;
840 double y1 = array[bin1];
841
842 double y2 = array[bin2];
843 double m = ( y2 - y1 ) / ( x2 - x1 );
844 double c = y1 - m * x1;
845 double result = m * xp + c;
846
847 return result;
848}
849
850double EvtBtoXsgammaKagan::FermiFunc( double y, const std::vector<double>& coeffs )
851{
852 //Fermi shape functions :1=exponential, 2=gaussian, 3=roman
853 if ( int( coeffs[0] ) == 1 )
854 return EvtBtoXsgammaFermiUtil::FermiExpFunc( y, coeffs );
855 if ( int( coeffs[0] ) == 2 )
856 return EvtBtoXsgammaFermiUtil::FermiGaussFunc( y, coeffs );
857 if ( int( coeffs[0] ) == 3 )
858 return EvtBtoXsgammaFermiUtil::FermiRomanFunc( y, coeffs );
859 return 1.;
860}
861
863{
864 return -log( fabs( 1. - y ) ) / y;
865}
866
868{
869 double li2( 0 );
870 for ( int i = 1; i < 1000;
871 i++ ) { //the value 1000 should actually be Infinite...
872 li2 += pow( y, i ) / ( i * i );
873 }
874 return li2;
875}
EvtComplex exp(const EvtComplex &c)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
Definition EvtReport.cpp:32
@ EVTGEN_WARNING
Definition EvtReport.hh:50
@ EVTGEN_ERROR
Definition EvtReport.hh:49
static double FermiGaussFunc(double, std::vector< double > const &coeffs)
static double FermiRomanFunc(double, std::vector< double > const &coeffs)
static double FermiRomanFuncRoot(double, double)
static double FermiGaussFuncRoot(double, double, double, std::vector< double > &coeffs)
static double FermiExpFunc(double var, const std::vector< double > &coeffs)
static double Gamma(double, const std::vector< double > &coeffs)
static double s77(double)
static double ImG(double)
static double ReG(double)
static double DeltaFermiFunc(double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2, const std::vector< double > &coeffs3)
static double s77FermiFunc(double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2)
static double s78FermiFunc(double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2)
static double s27Func(double var, const std::vector< double > &coeffs)
static double s88(double, double, double)
static double s78(double)
double GetMass(int code) override
static double GetArrayVal(double, double, double, double, std::vector< double >)
std::vector< double > m_brHad
void init(int, double *) override
static double FermiFunc(double, const std::vector< double > &coeffs)
static double s22Func(double var, const std::vector< double > &coeffs)
static double sFermiFunc(double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2, const std::vector< double > &coeffs3, const std::vector< double > &coeffs4)
static double diLogMathematica(double)
void computeHadronicMass(int, double *)
static double diLogFunc(double)
std::vector< double > m_massHad
static double Delta(double, double)
static double s88FermiFunc(double, const std::vector< double > &coeffs1, const std::vector< double > &coeffs2, const std::vector< double > &coeffs3)
static const double pi
Definition EvtConst.hh:26
static double Flat()
Definition EvtRandom.cpp:95