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