CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
EvtBtoXsgammaKagan Class Reference

#include <EvtBtoXsgammaKagan.hh>

+ Inheritance diagram for EvtBtoXsgammaKagan:

Public Member Functions

 EvtBtoXsgammaKagan ()
 
virtual ~EvtBtoXsgammaKagan ()
 
void init (int, double *)
 
void computeHadronicMass (int, double *)
 
void getDefaultHadronicMass ()
 
double GetMass (int code)
 
double CalcAlphaS (double)
 
void CalcWilsonCoeffs ()
 
void CalcDelta ()
 
double Fz (double)
 
- Public Member Functions inherited from EvtBtoXsgammaAbsModel
 EvtBtoXsgammaAbsModel ()
 
virtual ~EvtBtoXsgammaAbsModel ()
 
virtual void init (int, double *)
 
virtual double GetMass (int code)=0
 

Detailed Description

Definition at line 29 of file EvtBtoXsgammaKagan.hh.

Constructor & Destructor Documentation

◆ EvtBtoXsgammaKagan()

EvtBtoXsgammaKagan::EvtBtoXsgammaKagan ( )
inline

Definition at line 33 of file EvtBtoXsgammaKagan.hh.

33{}

◆ ~EvtBtoXsgammaKagan()

EvtBtoXsgammaKagan::~EvtBtoXsgammaKagan ( )
virtual

Definition at line 59 of file EvtBtoXsgammaKagan.cc.

59 {
60 delete [] massHad;
61 delete [] brHad;
62}

Member Function Documentation

◆ CalcAlphaS()

double EvtBtoXsgammaKagan::CalcAlphaS ( double  scale)

Definition at line 438 of file EvtBtoXsgammaKagan.cc.

438 {
439
440 double v = 1. -_beta0*(_alphasmZ/(2.*EvtConst::pi))*(log(_mZ/scale));
441 return (_alphasmZ/v)*(1. - ((_beta1/_beta0)*(_alphasmZ/(4.*EvtConst::pi))*(log(v)/v)));
442
443}
v
**********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
EvtConst::pi
static const double pi
Definition: EvtConst.hh:28

Referenced by computeHadronicMass().

◆ CalcDelta()

void EvtBtoXsgammaKagan::CalcDelta ( )

Definition at line 572 of file EvtBtoXsgammaKagan.cc.

572 {
573
574 double cDelta77 = (1. + (_alphasmu/(2.*EvtConst::pi)) *(_r7 - (16./3.) + _gam77*log(_mb/_mu)) + ( (pow((1. - _z),4.)/_fz) - 1.)*(6.*_lam2/pow(_mb,2.)) + (_alphasmubar/(2.*EvtConst::pi))*_kappabar )*pow(_c70mu,2.);
575
576 double cDelta27 = ((_alphasmu/(2.*EvtConst::pi))*(_rer2 + _gam27*log(_mb/_mu)) - (_lam2/(9.*_z*pow(_mb,2.))))*_c2mu*_c70mu;
577
578 double cDelta78 = (_alphasmu/(2.*EvtConst::pi))*(_rer8 + _gam87*log(_mb/_mu))*_c70mu*_c80mu;
579
580 _cDeltatot = cDelta77 + cDelta27 + cDelta78 + (_alphasmu/(2.*EvtConst::pi))*_c71mu*_c70mu + (_alpha/_alphasmu)*(2.*_c7emmu*_c70mu - _kSLemmu*pow(_c70mu,2.));
581
582}

Referenced by computeHadronicMass().

◆ CalcWilsonCoeffs()

void EvtBtoXsgammaKagan::CalcWilsonCoeffs ( )

Definition at line 445 of file EvtBtoXsgammaKagan.cc.

445 {
446
447 double mtatmw=_mt*pow((_alphasmW/_alphasmt),(12./23.))*(1 + (12./23.)*((253./18.) - (116./23.))*((_alphasmW - _alphasmt)/(4.0*EvtConst::pi)) - (4./3.)*(_alphasmt/EvtConst::pi));
448 double xt=pow(mtatmw,2.)/pow(_mW,2.);
449
450
451
452 /////LO
453 _c2mu = .5*pow(_etamu,(-12./23.)) + .5*pow(_etamu,(6./23.));
454
455 double c7mWsm = ((3.*pow(xt,3.) - 2.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
456 + ((-8.*pow(xt,3.) - 5.*pow(xt,2.) + 7.*xt)/(24.*pow((xt - 1.),3.) )) ;
457
458 double c8mWsm = ((-3.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
459 + ((- pow(xt,3.) + 5.*pow(xt,2.) + 2.*xt)/(8.*pow((xt - 1.),3.)));
460
461 double c7constmu = (626126./272277.)*pow(_etamu,(14./23.))
462 - (56281./51730.)*pow(_etamu,(16./23.)) - (3./7.)*pow(_etamu,(6./23.))
463 - (1./14.)*pow(_etamu,(-12./23.)) - .6494*pow(_etamu,.4086) - .038*pow(_etamu,-.423)
464 - .0186*pow(_etamu,-.8994) - .0057*pow(_etamu,.1456);
465
466 _c70mu = c7mWsm*pow(_etamu,(16./23.)) + (8./3.)*(pow(_etamu,(14./23.))
467 -pow(_etamu,(16./23.)))*c8mWsm + c7constmu;
468
469 double c8constmu = (313063./363036.)*pow(_etamu,(14./23.))
470 -.9135*pow(_etamu,.4086) + .0873*pow(_etamu,-.423) - .0571*pow(_etamu,-.8994)
471 + .0209*pow(_etamu,.1456);
472
473 _c80mu = c8mWsm*pow(_etamu,(14./23.)) + c8constmu;
474
475 //Compute the dilogarithm (PolyLog(2,x)) with the Simpson integrator
476 //The dilogarithm is defined as: Li_2(x)=Int_0^x(-log(1.-z)/z)
477 //however, Mathematica implements it as Sum[z^k/k^2,{k,1,Infinity}], so, althought the two
478 //results are similar and both implemented in the program, we prefer to use the
479 //one closer to the Mathematica implementation as identical to what used by the theorists.
480
481 // EvtItgFunction *myDiLogFunc = new EvtItgFunction(&diLogFunc, 0., 1.-1./xt);
482 //EvtItgAbsIntegrator *myDiLogSimp = new EvtItgSimpsonIntegrator(*myDiLogFunc, 1.0e-4, 50);
483 //double li2 = myDiLogSimp->evaluate(1.0e-20,1.-1./xt);
484
485 double li2=diLogMathematica(1.-1./xt);
486
487double c7mWsm1 = ( (-16. *pow(xt,4.) -122. *pow(xt,3.) + 80. *pow(xt,2.) -8. *xt)/
488(9. *pow((xt -1.),4.)) * li2 +
489(6. *pow(xt,4.) + 46. *pow(xt,3.) -28. *pow(xt,2.))/(3. *pow((xt-1.),5.)) *pow(log(xt),2.)
490+ (-102. *pow(xt,5.) -588. *pow(xt,4.) -2262. *pow(xt,3.) + 3244. *pow(xt,2.) -1364. *xt
491+ 208.)/(81. *pow((xt-1),5.)) *log(xt)
492+ (1646. *pow(xt,4.) + 12205. *pow(xt,3.) -10740. *pow(xt,2.) + 2509. *xt -436.)/
493(486. *pow((xt-1),4.)) );
494
495double c8mWsm1 = ((-4. *pow(xt,4.) + 40. *pow(xt,3.) + 41. *pow(xt,2.) + xt)/
496(6. *pow((xt-1.),4.)) * li2
497+ (-17. *pow(xt,3.) -31. *pow(xt,2.))/(2. *pow((xt-1.),5.) ) *pow(log(xt),2.)
498+ (-210. *pow(xt,5.) + 1086. *pow(xt,4.) + 4893. *pow(xt,3.) + 2857. *pow(xt,2.)
499-1994. *xt + 280.)/(216. *pow((xt-1),5.)) *log(xt)
500+ (737. *pow(xt,4.) -14102. *pow(xt,3.) -28209. *pow(xt,2.) + 610. *xt -508.)/
501(1296. *pow((xt-1),4.)) );
502
503double E1 = (xt *(18. -11. *xt -pow(xt,2.))/(12.*pow( (1. -xt),3.))
504+ pow(xt,2.)* (15. -16. *xt + 4. *pow(xt,2.))/(6. *pow((1. -xt),4.)) *log(xt)
505-2./3. *log(xt) );
506
507double e1 = 4661194./816831.;
508double e2 = -8516./2217. ;
509double e3 = 0.;
510double e4 = 0.;
511double e5 = -1.9043;
512double e6 = -.1008;
513double e7 = .1216;
514double e8 = .0183;
515
516double f1 = -17.3023;
517double f2 = 8.5027;
518double f3 = 4.5508;
519double f4 = .7519;
520double f5 = 2.004;
521double f6 = .7476;
522double f7 = -.5385;
523double f8 = .0914;
524
525double g1 = 14.8088;
526double g2 = -10.809;
527double g3 = -.874;
528double g4 = .4218;
529double g5 = -2.9347;
530double g6 = .3971;
531double g7 = .1600;
532double g8 = .0225;
533
534
535double c71constmu = ((e1 *_etamu *E1 + f1 + g1 *_etamu) *pow(_etamu,(14./23.))
536+ (e2 *_etamu *E1 + f2 + g2 *_etamu) *pow(_etamu,(16./23.))
537+ (e3 *_etamu *E1 + f3 + g3 *_etamu) *pow(_etamu,(6./23.))
538+ (e4 *_etamu *E1 + f4 + g4 *_etamu) *pow(_etamu,(-12./23.))
539+ (e5 *_etamu *E1 + f5 + g5 *_etamu) *pow(_etamu,.4086)
540+ (e6 *_etamu *E1 + f6 + g6 *_etamu) *pow(_etamu,(-.423))
541+ (e7 *_etamu *E1 + f7 + g7 *_etamu) *pow(_etamu,(-.8994))
542+ (e8 *_etamu *E1 + f8 + g8 *_etamu) *pow(_etamu,.1456 ));
543
544double c71pmu = ( ((297664./14283. *pow(_etamu,(16./23.))
545-7164416./357075. *pow(_etamu,(14./23.))
546+ 256868./14283. *pow(_etamu,(37./23.)) - 6698884./357075. *pow(_etamu,(39./23.)))
547*(c8mWsm))
548+ 37208./4761. *(pow(_etamu,(39./23.)) - pow(_etamu,(16./23.))) *(c7mWsm)
549+ c71constmu );
550
551_c71mu = (_alphasmW/_alphasmu *(pow(_etamu,(16./23.))* c7mWsm1 + 8./3. *(pow(_etamu,(14./23.))
552- pow(_etamu,(16./23.)) ) *c8mWsm1 ) + c71pmu);
553
554_c7emmu = ((32./75. *pow(_etamu,(-9./23.)) - 40./69. *pow(_etamu,(-7./23.)) +
555 88./575. *pow(_etamu,(16./23.))) *c7mWsm + (-32./575. *pow(_etamu,(-9./23.)) +
556 32./1449. *pow(_etamu,(-7./23.)) + 640./1449.*pow(_etamu,(14./23.)) -
557 704./1725.*pow(_etamu,(16./23.)) ) *c8mWsm
558 - 190./8073.*pow(_etamu,(-35./23.)) - 359./3105. *pow(_etamu,(-17./23.)) +
559 4276./121095. *pow(_etamu,(-12./23.)) + 350531./1009125.*pow(_etamu,(-9./23.))
560 + 2./4347. *pow(_etamu,(-7./23.)) - 5956./15525. *pow(_etamu,(6./23.)) +
561 38380./169533. *pow(_etamu,(14./23.)) - 748./8625. *pow(_etamu,(16./23.)));
562
563// Wilson coefficients values as according to Kagan's program
564// _c2mu=1.10566;
565//_c70mu=-0.314292;
566// _c80mu=-0.148954;
567// _c71mu=0.480964;
568// _c7emmu=0.0323219;
569
570}
f1
TFile * f1
Definition: DataBase/tau_mode.c:4
g1
TF1 * g1
Definition: DataBase/tau_mode.c:66
e1
Double_t e1
Definition: DataBase/tau_mode.c:7
e2
Double_t e2
Definition: DataBase/tau_mode.c:7

Referenced by computeHadronicMass().

◆ computeHadronicMass()

void EvtBtoXsgammaKagan::computeHadronicMass ( int  nArg,
double *  args 
)

Definition at line 133 of file EvtBtoXsgammaKagan.cc.

133 {
134
135 //Input parameters
136 int fermiFunction = (int)args[1];
137 _mB = args[2];
138 _mb = args[3];
139 _mu = args[4];
140 _lam1 = args[5];
141 _delta = args[6];
142 _z = args[7];
143 _nIntervalS = args[8];
144 _nIntervalmH = args[9];
145 std::vector<double> mHVect(int(_nIntervalmH+1.0));
146 massHad = new double[int(_nIntervalmH+1.0)];
147 brHad = new double[int(_nIntervalmH+1.0)];
148 intervalMH=_nIntervalmH;
149
150 //Going to have to add a new entry into the data file - takes ages...
151 report(WARNING,"EvtGen") << "EvtBtoXsgammaKagan: calculating new hadronic mass spectra. This takes a while..." << endl;
152
153 //Now need to compute the mHVect vector for
154 //the current parameters
155
156 //A few more parameters
157 double _mubar = _mu;
158 _mW = 80.33;
159 _mt = 175.0;
160 _alpha = 1./137.036;
161 _lambdabar = _mB - _mb;
162 _kappabar = 3.382 - 4.14*(sqrt(_z) - 0.29);
163 _fz=Fz(_z);
164 _rer8 = (44./9.) - (8./27.)*pow(EvtConst::pi,2.);
165 _r7 = (-10./3.) - (8./9.)*pow(EvtConst::pi,2.);
166 _rer2 = -4.092 + 12.78*(sqrt(_z) -.29);
167 _gam77 = 32./3.;
168 _gam27 = 416./81.;
169 _gam87 = -32./9.;
170 _lam2 = .12;
171 _beta0 = 23./3.;
172 _beta1 = 116./3.;
173 _alphasmZ = .118;
174 _mZ = 91.187;
175 _ms = _mb/50.;
176
177 double eGammaMin = 0.5*_mB*(1. - _delta);
178 double eGammaMax = 0.5*_mB;
179 double yMin = 2.*eGammaMin/_mB;
180 double yMax = 2.*eGammaMax/_mB;
181 double _CKMrat= 0.976;
182 double Nsl = 1.0;
183
184 //Calculate alpha the various scales
185 _alphasmW = CalcAlphaS(_mW);
186 _alphasmt = CalcAlphaS(_mt);
187 _alphasmu = CalcAlphaS(_mu);
188 _alphasmubar = CalcAlphaS(_mubar);
189
190 //Calculate the Wilson Coefficients and Delta
191 _etamu = _alphasmW/_alphasmu;
192 _kSLemmu = (12./23.)*((1./_etamu) -1.);
193 CalcWilsonCoeffs();
194 CalcDelta();
195
196 //Build s22 and s27 vector - saves time because double
197 //integration is required otherwise
198 std::vector<double> s22Coeffs(int(_nIntervalS+1.0));
199 std::vector<double> s27Coeffs(int(_nIntervalS+1.0));
200 std::vector<double> s28Coeffs(int(_nIntervalS+1.0));
201
202 double dy = (yMax - yMin)/_nIntervalS;
203 double yp = yMin;
204
205 std::vector<double> sCoeffs(1);
206 sCoeffs[0] = _z;
207
208 //Define s22 and s27 functions
209 EvtItgPtrFunction *mys22Func = new EvtItgPtrFunction(&s22Func, 0., yMax+0.1, sCoeffs);
210 EvtItgPtrFunction *mys27Func = new EvtItgPtrFunction(&s27Func, 0., yMax+0.1, sCoeffs);
211
212 //Use a simpson integrator
213 EvtItgAbsIntegrator *mys22Simp = new EvtItgSimpsonIntegrator(*mys22Func, 1.0e-4, 20);
214 EvtItgAbsIntegrator *mys27Simp = new EvtItgSimpsonIntegrator(*mys27Func, 1.0e-4, 50);
215
216 int i;
217
218 for (i=0;i<int(_nIntervalS+1.0);i++) {
219
220 s22Coeffs[i] = (16./27.)*mys22Simp->evaluate(1.0e-20,yp);
221 s27Coeffs[i] = (-8./9.)*_z*mys27Simp->evaluate(1.0e-20,yp);
222 s28Coeffs[i] = -s27Coeffs[i]/3.;
223 yp = yp + dy;
224
225 }
226
227 delete mys22Func;
228 delete mys27Func;
229 delete mys22Simp;
230 delete mys27Simp;
231
232 //Define functions and vectors used to calculate mHVect. Each function takes a set
233 //of vectors which are used as the function coefficients
234 std::vector<double> FermiCoeffs(6);
235 std::vector<double> varCoeffs(3);
236 std::vector<double> DeltaCoeffs(1);
237 std::vector<double> s88Coeffs(2);
238 std::vector<double> sInitCoeffs(3);
239
240 varCoeffs[0] = _mB;
241 varCoeffs[1] = _mb;
242 varCoeffs[2] = 0.;
243
244 DeltaCoeffs[0] = _alphasmu;
245
246 s88Coeffs[0] = _mb;
247 s88Coeffs[1] = _ms;
248
249 sInitCoeffs[0] = _nIntervalS;
250 sInitCoeffs[1] = yMin;
251 sInitCoeffs[2] = yMax;
252
253 FermiCoeffs[0]=fermiFunction;
254 FermiCoeffs[1]=0.0;
255 FermiCoeffs[2]=0.0;
256 FermiCoeffs[3]=0.0;
257 FermiCoeffs[4]=0.0;
258 FermiCoeffs[5]=0.0;
259
260 //Coefficients for gamma function
261 std::vector<double> gammaCoeffs(6);
262 gammaCoeffs[0]=76.18009172947146;
263 gammaCoeffs[1]=-86.50532032941677;
264 gammaCoeffs[2]=24.01409824083091;
265 gammaCoeffs[3]=-1.231739572450155;
266 gammaCoeffs[4]=0.1208650973866179e-2;
267 gammaCoeffs[5]=-0.5395239384953e-5;
268
269 //Calculate quantities for the fermi function to be used
270 //Distinguish among the different shape functions
271 if (fermiFunction == 1) {
272
273 FermiCoeffs[1]=_lambdabar;
274 FermiCoeffs[2]=(-3.*pow(_lambdabar,2.)/_lam1) - 1.;
275 FermiCoeffs[3]=_lam1;
276 FermiCoeffs[4]=1.0;
277
278 EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiExpFunc, -_mb, _mB-_mb, FermiCoeffs);
279 EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
280 FermiCoeffs[4]=myNormSimp->normalisation();
281 delete myNormFunc; myNormFunc=0;
282 delete myNormSimp; myNormSimp=0;
283
284 } else if (fermiFunction == 2) {
285
286 double a = EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(_lambdabar, _lam1, _mb, gammaCoeffs);
287 FermiCoeffs[1]=_lambdabar;
288 FermiCoeffs[2]=a;
289 FermiCoeffs[3]= EvtBtoXsgammaFermiUtil::Gamma((2.0 + a)/2., gammaCoeffs)/
290 EvtBtoXsgammaFermiUtil::Gamma((1.0 + a)/2., gammaCoeffs);
291 FermiCoeffs[4]=1.0;
292
293 EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiGaussFunc, -_mb, _mB-_mb, FermiCoeffs);
294 EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
295 FermiCoeffs[4]=myNormSimp->normalisation();
296 delete myNormFunc; myNormFunc=0;
297 delete myNormSimp; myNormSimp=0;
298
299 }
300 else if (fermiFunction == 3) {
301
302 double rho = EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(_lambdabar, _lam1);
303 FermiCoeffs[1]=_mB;
304 FermiCoeffs[2]=_mb;
305 FermiCoeffs[3]= rho;
306 FermiCoeffs[4]=_lambdabar;
307 FermiCoeffs[5]=1.0;
308
309 EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiRomanFunc, -_mb, _mB-_mb, FermiCoeffs);
310 EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
311 FermiCoeffs[5]=myNormSimp->normalisation();
312 delete myNormFunc; myNormFunc=0;
313 delete myNormSimp; myNormSimp=0;
314
315 }
316
317 //Define functions
318 EvtItgThreeCoeffFcn* myDeltaFermiFunc = new EvtItgThreeCoeffFcn(&DeltaFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, DeltaCoeffs);
319 EvtItgThreeCoeffFcn* mys88FermiFunc = new EvtItgThreeCoeffFcn(&s88FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, s88Coeffs);
320 EvtItgTwoCoeffFcn* mys77FermiFunc = new EvtItgTwoCoeffFcn(&s77FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
321 EvtItgTwoCoeffFcn* mys78FermiFunc = new EvtItgTwoCoeffFcn(&s78FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
322 EvtItgFourCoeffFcn* mys22FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s22Coeffs);
323 EvtItgFourCoeffFcn* mys27FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s27Coeffs);
324 EvtItgFourCoeffFcn* mys28FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s28Coeffs);
325
326 //Define integrators
327 EvtItgSimpsonIntegrator* myDeltaFermiSimp =
328 new EvtItgSimpsonIntegrator(*myDeltaFermiFunc, 1.0e-4, 40);
329 EvtItgSimpsonIntegrator* mys77FermiSimp =
330 new EvtItgSimpsonIntegrator(*mys77FermiFunc, 1.0e-4, 40);
331 EvtItgSimpsonIntegrator* mys88FermiSimp =
332 new EvtItgSimpsonIntegrator(*mys88FermiFunc, 1.0e-4, 40);
333 EvtItgSimpsonIntegrator* mys78FermiSimp =
334 new EvtItgSimpsonIntegrator(*mys78FermiFunc, 1.0e-4, 40);
335 EvtItgSimpsonIntegrator* mys22FermiSimp =
336 new EvtItgSimpsonIntegrator(*mys22FermiFunc, 1.0e-4, 40);
337 EvtItgSimpsonIntegrator* mys27FermiSimp =
338 new EvtItgSimpsonIntegrator(*mys27FermiFunc, 1.0e-4, 40);
339 EvtItgSimpsonIntegrator* mys28FermiSimp =
340 new EvtItgSimpsonIntegrator(*mys28FermiFunc, 1.0e-4, 40);
341
342 //Finally calculate mHVect for the range of hadronic masses
343 double mHmin = sqrt(_mB*_mB - 2.*_mB*eGammaMax);
344 double mHmax = sqrt(_mB*_mB - 2.*_mB*eGammaMin);
345 double dmH = (mHmax - mHmin)/_nIntervalmH;
346
347 double mH=mHmin;
348
349 //Calculating the Branching Fractions
350 for (i=0;i<int(_nIntervalmH+1.0);i++) {
351
352 double ymH = 1. - ((mH*mH)/(_mB*_mB));
353
354 //Need to set ymH as one of the input parameters
355 myDeltaFermiFunc->setCoeff(2, 2, ymH);
356 mys77FermiFunc->setCoeff(2, 2, ymH);
357 mys88FermiFunc->setCoeff(2, 2, ymH);
358 mys78FermiFunc->setCoeff(2, 2, ymH);
359 mys22FermiFunc->setCoeff(2, 2, ymH);
360 mys27FermiFunc->setCoeff(2, 2, ymH);
361 mys28FermiFunc->setCoeff(2, 2, ymH);
362
363 //Integrate
364
365 double deltaResult = myDeltaFermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
366 double s77Result = mys77FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
367 double s88Result = mys88FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
368 double s78Result = mys78FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
369 double s22Result = mys22FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
370 double s27Result = mys27FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
371
372 double py = (pow(_CKMrat,2.)*(6./_fz)*(_alpha/EvtConst::pi)*(deltaResult*_cDeltatot + (_alphasmu/EvtConst::pi)*(s77Result*pow(_c70mu,2.) + s27Result*_c2mu*(_c70mu - _c80mu/3.) + s78Result*_c70mu*_c80mu + s22Result*_c2mu*_c2mu + s88Result*_c80mu*_c80mu ) ) );
373
374 mHVect[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
375
376 massHad[i] = mH;
377 brHad[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
378
379 mH = mH+dmH;
380
381 }
382
383 //Clean up
384 delete myDeltaFermiFunc; myDeltaFermiFunc=0;
385 delete mys88FermiFunc; mys88FermiFunc=0;
386 delete mys77FermiFunc; mys77FermiFunc=0;
387 delete mys78FermiFunc; mys78FermiFunc=0;
388 delete mys22FermiFunc; mys22FermiFunc=0;
389 delete mys27FermiFunc; mys27FermiFunc=0;
390 delete mys28FermiFunc; mys28FermiFunc=0;
391
392 delete myDeltaFermiSimp; myDeltaFermiSimp=0;
393 delete mys77FermiSimp; mys77FermiSimp=0;
394 delete mys88FermiSimp; mys88FermiSimp=0;
395 delete mys78FermiSimp; mys78FermiSimp=0;
396 delete mys22FermiSimp; mys22FermiSimp=0;
397 delete mys27FermiSimp; mys27FermiSimp=0;
398 delete mys28FermiSimp; mys28FermiSimp=0;
399
400}
report
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
WARNING
@ WARNING
Definition: EvtReport.hh:50
EvtBtoXsgammaFermiUtil::FermiGaussFunc
static double FermiGaussFunc(double, std::vector< double > const &coeffs)
Definition: EvtBtoXsgammaFermiUtil.cc:48
EvtBtoXsgammaFermiUtil::FermiRomanFunc
static double FermiRomanFunc(double, std::vector< double > const &coeffs)
Definition: EvtBtoXsgammaFermiUtil.cc:187
EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot
static double FermiRomanFuncRoot(double, double)
Definition: EvtBtoXsgammaFermiUtil.cc:163
EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot
static double FermiGaussFuncRoot(double, double, double, std::vector< double > &coeffs)
Definition: EvtBtoXsgammaFermiUtil.cc:55
EvtBtoXsgammaFermiUtil::FermiExpFunc
static double FermiExpFunc(double var, const std::vector< double > &coeffs)
Definition: EvtBtoXsgammaFermiUtil.cc:40
EvtBtoXsgammaFermiUtil::Gamma
static double Gamma(double, const std::vector< double > &coeffs)
Definition: EvtBtoXsgammaFermiUtil.cc:100
EvtBtoXsgammaKagan::CalcAlphaS
double CalcAlphaS(double)
Definition: EvtBtoXsgammaKagan.cc:438
EvtItgAbsIntegrator::evaluate
double evaluate(double lower, double upper) const
Definition: EvtItgAbsIntegrator.cc:52
EvtItgAbsIntegrator::normalisation
double normalisation() const
Definition: EvtItgAbsIntegrator.cc:47
EvtItgFourCoeffFcn::setCoeff
virtual void setCoeff(int, int, double)
Definition: EvtItgFourCoeffFcn.cc:48
EvtItgThreeCoeffFcn::setCoeff
virtual void setCoeff(int, int, double)
Definition: EvtItgThreeCoeffFcn.cc:50
EvtItgTwoCoeffFcn::setCoeff
virtual void setCoeff(int, int, double)
Definition: EvtItgTwoCoeffFcn.cc:47

Referenced by init().

◆ Fz()

double EvtBtoXsgammaKagan::Fz ( double  z)

Definition at line 703 of file EvtBtoXsgammaKagan.cc.

703 {
704
705 return (1. -8.*z + 8.*pow(z,3.) - pow(z,4.) - 12.*pow(z,2.)*log(z));
706}

Referenced by computeHadronicMass().

◆ getDefaultHadronicMass()

void EvtBtoXsgammaKagan::getDefaultHadronicMass ( )

Definition at line 117 of file EvtBtoXsgammaKagan.cc.

117 {
118
119 massHad = new double[81];
120 brHad = new double[81];
121
122 double mass[81] = { 0, 0.0625995, 0.125199, 0.187798, 0.250398, 0.312997, 0.375597, 0.438196, 0.500796, 0.563395, 0.625995, 0.688594, 0.751194, 0.813793, 0.876392, 0.938992, 1.00159, 1.06419, 1.12679, 1.18939, 1.25199, 1.31459, 1.37719, 1.43979, 1.50239, 1.56499, 1.62759, 1.69019, 1.75278, 1.81538, 1.87798, 1.94058, 2.00318, 2.06578, 2.12838, 2.19098, 2.25358, 2.31618, 2.37878, 2.44138, 2.50398, 2.56658, 2.62918, 2.69178, 2.75438, 2.81698, 2.87958, 2.94217, 3.00477, 3.06737, 3.12997, 3.19257, 3.25517, 3.31777, 3.38037, 3.44297, 3.50557, 3.56817, 3.63077, 3.69337, 3.75597, 3.81857, 3.88117, 3.94377, 4.00637, 4.06896, 4.13156, 4.19416, 4.25676, 4.31936, 4.38196, 4.44456, 4.50716, 4.56976, 4.63236, 4.69496, 4.75756, 4.82016, 4.88276, 4.94536, 5.00796};
123
124 double br[81] = { 0, 1.03244e-09, 3.0239e-08, 1.99815e-07, 7.29392e-07, 1.93129e-06, 4.17806e-06, 7.86021e-06, 1.33421e-05, 2.09196e-05, 3.07815e-05, 4.29854e-05, 5.74406e-05, 7.3906e-05, 9.2003e-05, 0.000111223, 0.000130977, 0.000150618, 0.000169483, 0.000186934, 0.000202392, 0.000215366, 0.000225491, 0.000232496, 0.000236274, 0.000236835, 0.000234313, 0.000228942, 0.000221042, 0.000210994, 0.000199215, 0.000186137, 0.000172194, 0.000157775, 0.000143255, 0.000128952, 0.000115133, 0.000102012, 8.97451e-05, 7.84384e-05, 6.81519e-05, 5.89048e-05, 5.06851e-05, 4.34515e-05, 3.71506e-05, 3.1702e-05, 2.70124e-05, 2.30588e-05, 1.96951e-05, 1.68596e-05, 1.44909e-05, 1.25102e-05, 1.08596e-05, 9.48476e-06, 8.34013e-06, 7.38477e-06, 6.58627e-06, 5.91541e-06, 5.35022e-06, 4.87047e-06, 4.46249e-06, 4.11032e-06, 3.80543e-06, 3.54051e-06, 3.30967e-06, 3.10848e-06, 2.93254e-06, 2.78369e-06, 2.65823e-06, 2.55747e-06, 2.51068e-06, 2.57179e-06, 2.74684e-06, 3.02719e-06, 3.41182e-06, 3.91387e-06, 4.56248e-06, 5.40862e-06, 6.53915e-06, 8.10867e-06, 1.04167e-05 };
125
126 for(int i=0; i<81; i++){
127 massHad[i] = mass[i];
128 brHad[i] = br[i];
129 }
130 intervalMH=80;
131}
mass
double mass
Definition: CosmicGenerator.cxx:138

Referenced by init().

◆ GetMass()

double EvtBtoXsgammaKagan::GetMass ( int  code)
virtual

Implements EvtBtoXsgammaAbsModel.

Definition at line 402 of file EvtBtoXsgammaKagan.cc.

402 {
403
404// Get hadronic mass for the event according to the hadronic mass spectra computed in computeHadronicMass
405 double mass=0.0;
406 // double min=0.6373; // usually just above K pi threshold for Xsd/u
407 double min=_mHmin;
408 if(bbprod)min=1.1;
409 // double max=4.5;
410 double max=_mHmax;
411 double xbox(0), ybox(0);
412 double boxheight(0);
413 double trueHeight(0);
414 double boxwidth=max-min;
415
416 for (int i=0;i<int(intervalMH+1.0);i++) {
417 if(brHad[i]>boxheight)boxheight=brHad[i];
418 }
419 while ((mass > max) || (mass < min)){
420 xbox = EvtRandom::Flat(boxwidth)+min;
421 ybox=EvtRandom::Flat(boxheight);
422 trueHeight=0.0;
423 for (int i=0;i<int(intervalMH+1.0);i++) {
424 if(massHad[i]>=xbox&& trueHeight==0.0){
425 trueHeight=(brHad[i]+brHad[i+1])/2.;
426 }
427 }
428 if (ybox>trueHeight) {
429 mass=0.0;
430 } else {
431 mass=xbox;
432 }
433 }
434
435 return mass;
436}
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:73

◆ init()

void EvtBtoXsgammaKagan::init ( int  nArg,
double *  args 
)
virtual

Reimplemented from EvtBtoXsgammaAbsModel.

Definition at line 64 of file EvtBtoXsgammaKagan.cc.

64 {
65
66 if ((nArg) > 12 || (nArg > 1 && nArg <10) || nArg == 11){
67
68 report(ERROR,"EvtGen") << "EvtBtoXsgamma generator model "
69 << "EvtBtoXsgammaKagan expected "
70 << "either 1(default config) or "
71 << "10 (default mass range) or "
72 << "12 (user range) arguments but found: "
73 <<nArg<<endl;
74 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
75 ::abort();
76 }
77
78 if(nArg == 1){
79 bbprod = true;
80 getDefaultHadronicMass();
81 }else{
82 bbprod = false;
83 computeHadronicMass(nArg, args);
84 }
85
86 double mHminLimit=0.6373;
87 double mHmaxLimit=4.5;
88
89 if (nArg>10){
90 _mHmin = args[10];
91 _mHmax = args[11];
92 if (_mHmin > _mHmax){
93 report(ERROR,"EvtGen") << "Minimum hadronic mass exceeds maximum "
94 << endl;
95 report(ERROR,"EvtGen") << "Will terminate execution!" << endl;
96 ::abort();
97 }
98 if (_mHmin < mHminLimit){
99 report(ERROR,"EvtGen") << "Minimum hadronic mass below K pi threshold"
100 << endl;
101 report(ERROR,"EvtGen") << "Resetting to K pi threshold" << endl;
102 _mHmin = mHminLimit;
103 }
104 if (_mHmax > mHmaxLimit){
105 report(ERROR,"EvtGen") << "Maximum hadronic mass above 4.5 GeV/c^2"
106 << endl;
107 report(ERROR,"EvtGen") << "Resetting to 4.5 GeV/c^2" << endl;
108 _mHmax = mHmaxLimit;
109 }
110 }else{
111 _mHmin=mHminLimit; // usually just above K pi threshold for Xsd/u
112 _mHmax=mHmaxLimit;
113 }
114
115}
ERROR
@ ERROR
Definition: EvtReport.hh:49
EvtBtoXsgammaKagan::computeHadronicMass
void computeHadronicMass(int, double *)
Definition: EvtBtoXsgammaKagan.cc:133
The documentation for this class was generated from the following files: