Line data Source code
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 : //
31 : #include "EvtGenBase/EvtPatches.hh"
32 :
33 : #include <stdlib.h>
34 : #include "EvtGenModels/EvtBtoXsgamma.hh"
35 : #include "EvtGenBase/EvtRandom.hh"
36 : #include "EvtGenModels/EvtBtoXsgammaKagan.hh"
37 : #include <string>
38 : #include "EvtGenBase/EvtConst.hh"
39 : #include "EvtGenBase/EvtParticle.hh"
40 : #include "EvtGenBase/EvtGenKine.hh"
41 : #include "EvtGenBase/EvtPDL.hh"
42 : #include "EvtGenBase/EvtReport.hh"
43 : #include "EvtGenModels/EvtItgSimpsonIntegrator.hh"
44 : #include "EvtGenModels/EvtItgFunction.hh"
45 : #include "EvtGenModels/EvtItgPtrFunction.hh"
46 : #include "EvtGenModels/EvtItgTwoCoeffFcn.hh"
47 : #include "EvtGenModels/EvtItgThreeCoeffFcn.hh"
48 : #include "EvtGenModels/EvtItgFourCoeffFcn.hh"
49 : #include "EvtGenModels/EvtItgAbsIntegrator.hh"
50 : #include "EvtGenModels/EvtBtoXsgammaFermiUtil.hh"
51 :
52 : #include <fstream>
53 : using std::endl;
54 : using std::fstream;
55 :
56 : bool EvtBtoXsgammaKagan::bbprod = false;
57 : double EvtBtoXsgammaKagan::intervalMH = 0;
58 :
59 0 : EvtBtoXsgammaKagan::~EvtBtoXsgammaKagan(){
60 0 : delete [] massHad;
61 0 : delete [] brHad;
62 0 : }
63 :
64 : void EvtBtoXsgammaKagan::init(int nArg, double* args){
65 :
66 0 : if ((nArg) > 12 || (nArg > 1 && nArg <10) || nArg == 11){
67 :
68 0 : report(Severity::Error,"EvtGen") << "EvtBtoXsgamma generator model "
69 0 : << "EvtBtoXsgammaKagan expected "
70 0 : << "either 1(default config) or "
71 0 : << "10 (default mass range) or "
72 0 : << "12 (user range) arguments but found: "
73 0 : <<nArg<<endl;
74 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
75 0 : ::abort();
76 : }
77 :
78 0 : if(nArg == 1){
79 0 : bbprod = true;
80 0 : getDefaultHadronicMass();
81 0 : }else{
82 0 : bbprod = false;
83 0 : computeHadronicMass(nArg, args);
84 : }
85 :
86 : double mHminLimit=0.6373;
87 : double mHmaxLimit=4.5;
88 :
89 0 : if (nArg>10){
90 0 : _mHmin = args[10];
91 0 : _mHmax = args[11];
92 0 : if (_mHmin > _mHmax){
93 0 : report(Severity::Error,"EvtGen") << "Minimum hadronic mass exceeds maximum "
94 0 : << endl;
95 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!" << endl;
96 0 : ::abort();
97 : }
98 0 : if (_mHmin < mHminLimit){
99 0 : report(Severity::Error,"EvtGen") << "Minimum hadronic mass below K pi threshold"
100 0 : << endl;
101 0 : report(Severity::Error,"EvtGen") << "Resetting to K pi threshold" << endl;
102 0 : _mHmin = mHminLimit;
103 0 : }
104 0 : if (_mHmax > mHmaxLimit){
105 0 : report(Severity::Error,"EvtGen") << "Maximum hadronic mass above 4.5 GeV/c^2"
106 0 : << endl;
107 0 : report(Severity::Error,"EvtGen") << "Resetting to 4.5 GeV/c^2" << endl;
108 0 : _mHmax = mHmaxLimit;
109 0 : }
110 : }else{
111 0 : _mHmin=mHminLimit; // usually just above K pi threshold for Xsd/u
112 0 : _mHmax=mHmaxLimit;
113 : }
114 :
115 0 : }
116 :
117 : void EvtBtoXsgammaKagan::getDefaultHadronicMass(){
118 :
119 0 : massHad = new double[81];
120 0 : brHad = new double[81];
121 :
122 0 : 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 0 : 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 0 : for(int i=0; i<81; i++){
127 0 : massHad[i] = mass[i];
128 0 : brHad[i] = br[i];
129 : }
130 0 : intervalMH=80;
131 0 : }
132 :
133 : void EvtBtoXsgammaKagan::computeHadronicMass(int /*nArg*/, double* args){
134 :
135 : //Input parameters
136 0 : int fermiFunction = (int)args[1];
137 0 : _mB = args[2];
138 0 : _mb = args[3];
139 0 : _mu = args[4];
140 0 : _lam1 = args[5];
141 0 : _delta = args[6];
142 0 : _z = args[7];
143 0 : _nIntervalS = args[8];
144 0 : _nIntervalmH = args[9];
145 0 : std::vector<double> mHVect(int(_nIntervalmH+1.0));
146 0 : massHad = new double[int(_nIntervalmH+1.0)];
147 0 : brHad = new double[int(_nIntervalmH+1.0)];
148 0 : intervalMH=_nIntervalmH;
149 :
150 : //Going to have to add a new entry into the data file - takes ages...
151 0 : report(Severity::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 0 : double _mubar = _mu;
158 0 : _mW = 80.33;
159 0 : _mt = 175.0;
160 0 : _alpha = 1./137.036;
161 0 : _lambdabar = _mB - _mb;
162 0 : _kappabar = 3.382 - 4.14*(sqrt(_z) - 0.29);
163 0 : _fz=Fz(_z);
164 0 : _rer8 = (44./9.) - (8./27.)*pow(EvtConst::pi,2.);
165 0 : _r7 = (-10./3.) - (8./9.)*pow(EvtConst::pi,2.);
166 0 : _rer2 = -4.092 + 12.78*(sqrt(_z) -.29);
167 0 : _gam77 = 32./3.;
168 0 : _gam27 = 416./81.;
169 0 : _gam87 = -32./9.;
170 0 : _lam2 = .12;
171 0 : _beta0 = 23./3.;
172 0 : _beta1 = 116./3.;
173 0 : _alphasmZ = .118;
174 0 : _mZ = 91.187;
175 0 : _ms = _mb/50.;
176 :
177 0 : double eGammaMin = 0.5*_mB*(1. - _delta);
178 : double eGammaMax = 0.5*_mB;
179 0 : double yMin = 2.*eGammaMin/_mB;
180 0 : double yMax = 2.*eGammaMax/_mB;
181 : double _CKMrat= 0.976;
182 : double Nsl = 1.0;
183 :
184 : //Calculate alpha the various scales
185 0 : _alphasmW = CalcAlphaS(_mW);
186 0 : _alphasmt = CalcAlphaS(_mt);
187 0 : _alphasmu = CalcAlphaS(_mu);
188 0 : _alphasmubar = CalcAlphaS(_mubar);
189 :
190 : //Calculate the Wilson Coefficients and Delta
191 0 : _etamu = _alphasmW/_alphasmu;
192 0 : _kSLemmu = (12./23.)*((1./_etamu) -1.);
193 0 : CalcWilsonCoeffs();
194 0 : CalcDelta();
195 :
196 : //Build s22 and s27 vector - saves time because double
197 : //integration is required otherwise
198 0 : std::vector<double> s22Coeffs(int(_nIntervalS+1.0));
199 0 : std::vector<double> s27Coeffs(int(_nIntervalS+1.0));
200 0 : std::vector<double> s28Coeffs(int(_nIntervalS+1.0));
201 :
202 0 : double dy = (yMax - yMin)/_nIntervalS;
203 : double yp = yMin;
204 :
205 0 : std::vector<double> sCoeffs(1);
206 0 : sCoeffs[0] = _z;
207 :
208 : //Define s22 and s27 functions
209 0 : EvtItgPtrFunction *mys22Func = new EvtItgPtrFunction(&s22Func, 0., yMax+0.1, sCoeffs);
210 0 : EvtItgPtrFunction *mys27Func = new EvtItgPtrFunction(&s27Func, 0., yMax+0.1, sCoeffs);
211 :
212 : //Use a simpson integrator
213 0 : EvtItgAbsIntegrator *mys22Simp = new EvtItgSimpsonIntegrator(*mys22Func, 1.0e-4, 20);
214 0 : EvtItgAbsIntegrator *mys27Simp = new EvtItgSimpsonIntegrator(*mys27Func, 1.0e-4, 50);
215 :
216 : int i;
217 :
218 0 : for (i=0;i<int(_nIntervalS+1.0);i++) {
219 :
220 0 : s22Coeffs[i] = (16./27.)*mys22Simp->evaluate(1.0e-20,yp);
221 0 : s27Coeffs[i] = (-8./9.)*_z*mys27Simp->evaluate(1.0e-20,yp);
222 0 : s28Coeffs[i] = -s27Coeffs[i]/3.;
223 0 : yp = yp + dy;
224 :
225 : }
226 :
227 0 : delete mys22Func;
228 0 : delete mys27Func;
229 0 : delete mys22Simp;
230 0 : 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 0 : std::vector<double> FermiCoeffs(6);
235 0 : std::vector<double> varCoeffs(3);
236 0 : std::vector<double> DeltaCoeffs(1);
237 0 : std::vector<double> s88Coeffs(2);
238 0 : std::vector<double> sInitCoeffs(3);
239 :
240 0 : varCoeffs[0] = _mB;
241 0 : varCoeffs[1] = _mb;
242 0 : varCoeffs[2] = 0.;
243 :
244 0 : DeltaCoeffs[0] = _alphasmu;
245 :
246 0 : s88Coeffs[0] = _mb;
247 0 : s88Coeffs[1] = _ms;
248 :
249 0 : sInitCoeffs[0] = _nIntervalS;
250 0 : sInitCoeffs[1] = yMin;
251 0 : sInitCoeffs[2] = yMax;
252 :
253 0 : FermiCoeffs[0]=fermiFunction;
254 0 : FermiCoeffs[1]=0.0;
255 0 : FermiCoeffs[2]=0.0;
256 0 : FermiCoeffs[3]=0.0;
257 0 : FermiCoeffs[4]=0.0;
258 0 : FermiCoeffs[5]=0.0;
259 :
260 : //Coefficients for gamma function
261 0 : std::vector<double> gammaCoeffs(6);
262 0 : gammaCoeffs[0]=76.18009172947146;
263 0 : gammaCoeffs[1]=-86.50532032941677;
264 0 : gammaCoeffs[2]=24.01409824083091;
265 0 : gammaCoeffs[3]=-1.231739572450155;
266 0 : gammaCoeffs[4]=0.1208650973866179e-2;
267 0 : gammaCoeffs[5]=-0.5395239384953e-5;
268 :
269 : //Calculate quantities for the fermi function to be used
270 : //Distinguish among the different shape functions
271 0 : if (fermiFunction == 1) {
272 :
273 0 : FermiCoeffs[1]=_lambdabar;
274 0 : FermiCoeffs[2]=(-3.*pow(_lambdabar,2.)/_lam1) - 1.;
275 0 : FermiCoeffs[3]=_lam1;
276 0 : FermiCoeffs[4]=1.0;
277 :
278 0 : EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiExpFunc, -_mb, _mB-_mb, FermiCoeffs);
279 0 : EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
280 0 : FermiCoeffs[4]=myNormSimp->normalisation();
281 0 : delete myNormFunc; myNormFunc=0;
282 0 : delete myNormSimp; myNormSimp=0;
283 :
284 0 : } else if (fermiFunction == 2) {
285 :
286 0 : double a = EvtBtoXsgammaFermiUtil::FermiGaussFuncRoot(_lambdabar, _lam1, _mb, gammaCoeffs);
287 0 : FermiCoeffs[1]=_lambdabar;
288 0 : FermiCoeffs[2]=a;
289 0 : FermiCoeffs[3]= EvtBtoXsgammaFermiUtil::Gamma((2.0 + a)/2., gammaCoeffs)/
290 0 : EvtBtoXsgammaFermiUtil::Gamma((1.0 + a)/2., gammaCoeffs);
291 0 : FermiCoeffs[4]=1.0;
292 :
293 0 : EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiGaussFunc, -_mb, _mB-_mb, FermiCoeffs);
294 0 : EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
295 0 : FermiCoeffs[4]=myNormSimp->normalisation();
296 0 : delete myNormFunc; myNormFunc=0;
297 0 : delete myNormSimp; myNormSimp=0;
298 :
299 0 : }
300 0 : else if (fermiFunction == 3) {
301 :
302 0 : double rho = EvtBtoXsgammaFermiUtil::FermiRomanFuncRoot(_lambdabar, _lam1);
303 0 : FermiCoeffs[1]=_mB;
304 0 : FermiCoeffs[2]=_mb;
305 0 : FermiCoeffs[3]= rho;
306 0 : FermiCoeffs[4]=_lambdabar;
307 0 : FermiCoeffs[5]=1.0;
308 :
309 0 : EvtItgPtrFunction *myNormFunc = new EvtItgPtrFunction(&EvtBtoXsgammaFermiUtil::FermiRomanFunc, -_mb, _mB-_mb, FermiCoeffs);
310 0 : EvtItgAbsIntegrator *myNormSimp = new EvtItgSimpsonIntegrator(*myNormFunc, 1.0e-4, 40);
311 0 : FermiCoeffs[5]=myNormSimp->normalisation();
312 0 : delete myNormFunc; myNormFunc=0;
313 0 : delete myNormSimp; myNormSimp=0;
314 :
315 0 : }
316 :
317 : //Define functions
318 0 : EvtItgThreeCoeffFcn* myDeltaFermiFunc = new EvtItgThreeCoeffFcn(&DeltaFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, DeltaCoeffs);
319 0 : EvtItgThreeCoeffFcn* mys88FermiFunc = new EvtItgThreeCoeffFcn(&s88FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, s88Coeffs);
320 0 : EvtItgTwoCoeffFcn* mys77FermiFunc = new EvtItgTwoCoeffFcn(&s77FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
321 0 : EvtItgTwoCoeffFcn* mys78FermiFunc = new EvtItgTwoCoeffFcn(&s78FermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs);
322 0 : EvtItgFourCoeffFcn* mys22FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s22Coeffs);
323 0 : EvtItgFourCoeffFcn* mys27FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s27Coeffs);
324 0 : EvtItgFourCoeffFcn* mys28FermiFunc = new EvtItgFourCoeffFcn(&sFermiFunc, -_mb, _mB-_mb, FermiCoeffs, varCoeffs, sInitCoeffs, s28Coeffs);
325 :
326 : //Define integrators
327 : EvtItgSimpsonIntegrator* myDeltaFermiSimp =
328 0 : new EvtItgSimpsonIntegrator(*myDeltaFermiFunc, 1.0e-4, 40);
329 : EvtItgSimpsonIntegrator* mys77FermiSimp =
330 0 : new EvtItgSimpsonIntegrator(*mys77FermiFunc, 1.0e-4, 40);
331 : EvtItgSimpsonIntegrator* mys88FermiSimp =
332 0 : new EvtItgSimpsonIntegrator(*mys88FermiFunc, 1.0e-4, 40);
333 : EvtItgSimpsonIntegrator* mys78FermiSimp =
334 0 : new EvtItgSimpsonIntegrator(*mys78FermiFunc, 1.0e-4, 40);
335 : EvtItgSimpsonIntegrator* mys22FermiSimp =
336 0 : new EvtItgSimpsonIntegrator(*mys22FermiFunc, 1.0e-4, 40);
337 : EvtItgSimpsonIntegrator* mys27FermiSimp =
338 0 : new EvtItgSimpsonIntegrator(*mys27FermiFunc, 1.0e-4, 40);
339 : EvtItgSimpsonIntegrator* mys28FermiSimp =
340 0 : new EvtItgSimpsonIntegrator(*mys28FermiFunc, 1.0e-4, 40);
341 :
342 : //Finally calculate mHVect for the range of hadronic masses
343 0 : double mHmin = sqrt(_mB*_mB - 2.*_mB*eGammaMax);
344 0 : double mHmax = sqrt(_mB*_mB - 2.*_mB*eGammaMin);
345 0 : double dmH = (mHmax - mHmin)/_nIntervalmH;
346 :
347 : double mH=mHmin;
348 :
349 : //Calculating the Branching Fractions
350 0 : for (i=0;i<int(_nIntervalmH+1.0);i++) {
351 :
352 0 : double ymH = 1. - ((mH*mH)/(_mB*_mB));
353 :
354 : //Need to set ymH as one of the input parameters
355 0 : myDeltaFermiFunc->setCoeff(2, 2, ymH);
356 0 : mys77FermiFunc->setCoeff(2, 2, ymH);
357 0 : mys88FermiFunc->setCoeff(2, 2, ymH);
358 0 : mys78FermiFunc->setCoeff(2, 2, ymH);
359 0 : mys22FermiFunc->setCoeff(2, 2, ymH);
360 0 : mys27FermiFunc->setCoeff(2, 2, ymH);
361 0 : mys28FermiFunc->setCoeff(2, 2, ymH);
362 :
363 : //Integrate
364 :
365 0 : double deltaResult = myDeltaFermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
366 0 : double s77Result = mys77FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
367 0 : double s88Result = mys88FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
368 0 : double s78Result = mys78FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
369 0 : double s22Result = mys22FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
370 0 : double s27Result = mys27FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
371 0 : mys28FermiSimp->evaluate((_mB*ymH-_mb),_mB-_mb);
372 :
373 0 : 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 ) ) );
374 :
375 0 : mHVect[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
376 :
377 0 : massHad[i] = mH;
378 0 : brHad[i] = 2.*(mH/(_mB*_mB))*0.105*Nsl*py;
379 :
380 0 : mH = mH+dmH;
381 :
382 : }
383 :
384 : //Clean up
385 0 : delete myDeltaFermiFunc; myDeltaFermiFunc=0;
386 0 : delete mys88FermiFunc; mys88FermiFunc=0;
387 0 : delete mys77FermiFunc; mys77FermiFunc=0;
388 0 : delete mys78FermiFunc; mys78FermiFunc=0;
389 0 : delete mys22FermiFunc; mys22FermiFunc=0;
390 0 : delete mys27FermiFunc; mys27FermiFunc=0;
391 0 : delete mys28FermiFunc; mys28FermiFunc=0;
392 :
393 0 : delete myDeltaFermiSimp; myDeltaFermiSimp=0;
394 0 : delete mys77FermiSimp; mys77FermiSimp=0;
395 0 : delete mys88FermiSimp; mys88FermiSimp=0;
396 0 : delete mys78FermiSimp; mys78FermiSimp=0;
397 0 : delete mys22FermiSimp; mys22FermiSimp=0;
398 0 : delete mys27FermiSimp; mys27FermiSimp=0;
399 0 : delete mys28FermiSimp; mys28FermiSimp=0;
400 :
401 0 : }
402 :
403 : double EvtBtoXsgammaKagan::GetMass( int /*Xscode*/ ){
404 :
405 : // Get hadronic mass for the event according to the hadronic mass spectra computed in computeHadronicMass
406 : double mass=0.0;
407 0 : double min=_mHmin;
408 0 : if(bbprod)min=1.1;
409 : // double max=4.5;
410 0 : double max=_mHmax;
411 : double xbox(0), ybox(0);
412 : double boxheight(0);
413 : double trueHeight(0);
414 0 : double boxwidth=max-min;
415 : double wgt(0.);
416 :
417 0 : for (int i=0;i<int(intervalMH+1.0);i++) {
418 0 : if(brHad[i]>boxheight)boxheight=brHad[i];
419 : }
420 0 : while ((mass > max) || (mass < min)){
421 0 : xbox = EvtRandom::Flat(boxwidth)+min;
422 0 : ybox=EvtRandom::Flat(boxheight);
423 : trueHeight=0.0;
424 : // Correction by Peter Richardson
425 0 : for( int i = 1 ; i < int( intervalMH + 1.0 ) ; ++i ) {
426 0 : if ( ( massHad[i] >= xbox ) && ( 0.0 == trueHeight ) ) {
427 0 : wgt=(xbox-massHad[i-1])/(massHad[i]-massHad[i-1]);
428 0 : trueHeight=brHad[i-1]+wgt*(brHad[i]-brHad[i-1]);
429 0 : }
430 : }
431 :
432 0 : if (ybox>trueHeight) {
433 : mass=0.0;
434 0 : } else {
435 : mass=xbox;
436 : }
437 : }
438 :
439 0 : return mass;
440 : }
441 :
442 : double EvtBtoXsgammaKagan::CalcAlphaS(double scale) {
443 :
444 0 : double v = 1. -_beta0*(_alphasmZ/(2.*EvtConst::pi))*(log(_mZ/scale));
445 0 : return (_alphasmZ/v)*(1. - ((_beta1/_beta0)*(_alphasmZ/(4.*EvtConst::pi))*(log(v)/v)));
446 :
447 : }
448 :
449 : void EvtBtoXsgammaKagan::CalcWilsonCoeffs( ){
450 :
451 0 : 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));
452 0 : double xt=pow(mtatmw,2.)/pow(_mW,2.);
453 :
454 :
455 :
456 : /////LO
457 0 : _c2mu = .5*pow(_etamu,(-12./23.)) + .5*pow(_etamu,(6./23.));
458 :
459 0 : double c7mWsm = ((3.*pow(xt,3.) - 2.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
460 0 : + ((-8.*pow(xt,3.) - 5.*pow(xt,2.) + 7.*xt)/(24.*pow((xt - 1.),3.) )) ;
461 :
462 0 : double c8mWsm = ((-3.*pow(xt,2.))/(4.*pow((xt - 1.),4.)))*log(xt)
463 0 : + ((- pow(xt,3.) + 5.*pow(xt,2.) + 2.*xt)/(8.*pow((xt - 1.),3.)));
464 :
465 0 : double c7constmu = (626126./272277.)*pow(_etamu,(14./23.))
466 0 : - (56281./51730.)*pow(_etamu,(16./23.)) - (3./7.)*pow(_etamu,(6./23.))
467 0 : - (1./14.)*pow(_etamu,(-12./23.)) - .6494*pow(_etamu,.4086) - .038*pow(_etamu,-.423)
468 0 : - .0186*pow(_etamu,-.8994) - .0057*pow(_etamu,.1456);
469 :
470 0 : _c70mu = c7mWsm*pow(_etamu,(16./23.)) + (8./3.)*(pow(_etamu,(14./23.))
471 0 : -pow(_etamu,(16./23.)))*c8mWsm + c7constmu;
472 :
473 0 : double c8constmu = (313063./363036.)*pow(_etamu,(14./23.))
474 0 : -.9135*pow(_etamu,.4086) + .0873*pow(_etamu,-.423) - .0571*pow(_etamu,-.8994)
475 0 : + .0209*pow(_etamu,.1456);
476 :
477 0 : _c80mu = c8mWsm*pow(_etamu,(14./23.)) + c8constmu;
478 :
479 : //Compute the dilogarithm (PolyLog(2,x)) with the Simpson integrator
480 : //The dilogarithm is defined as: Li_2(x)=Int_0^x(-log(1.-z)/z)
481 : //however, Mathematica implements it as Sum[z^k/k^2,{k,1,Infinity}], so, althought the two
482 : //results are similar and both implemented in the program, we prefer to use the
483 : //one closer to the Mathematica implementation as identical to what used by the theorists.
484 :
485 : // EvtItgFunction *myDiLogFunc = new EvtItgFunction(&diLogFunc, 0., 1.-1./xt);
486 : //EvtItgAbsIntegrator *myDiLogSimp = new EvtItgSimpsonIntegrator(*myDiLogFunc, 1.0e-4, 50);
487 : //double li2 = myDiLogSimp->evaluate(1.0e-20,1.-1./xt);
488 :
489 0 : double li2=diLogMathematica(1.-1./xt);
490 :
491 0 : double c7mWsm1 = ( (-16. *pow(xt,4.) -122. *pow(xt,3.) + 80. *pow(xt,2.) -8. *xt)/
492 0 : (9. *pow((xt -1.),4.)) * li2 +
493 0 : (6. *pow(xt,4.) + 46. *pow(xt,3.) -28. *pow(xt,2.))/(3. *pow((xt-1.),5.)) *pow(log(xt),2.)
494 0 : + (-102. *pow(xt,5.) -588. *pow(xt,4.) -2262. *pow(xt,3.) + 3244. *pow(xt,2.) -1364. *xt
495 0 : + 208.)/(81. *pow((xt-1),5.)) *log(xt)
496 0 : + (1646. *pow(xt,4.) + 12205. *pow(xt,3.) -10740. *pow(xt,2.) + 2509. *xt -436.)/
497 0 : (486. *pow((xt-1),4.)) );
498 :
499 0 : double c8mWsm1 = ((-4. *pow(xt,4.) + 40. *pow(xt,3.) + 41. *pow(xt,2.) + xt)/
500 0 : (6. *pow((xt-1.),4.)) * li2
501 0 : + (-17. *pow(xt,3.) -31. *pow(xt,2.))/(2. *pow((xt-1.),5.) ) *pow(log(xt),2.)
502 0 : + (-210. *pow(xt,5.) + 1086. *pow(xt,4.) + 4893. *pow(xt,3.) + 2857. *pow(xt,2.)
503 0 : -1994. *xt + 280.)/(216. *pow((xt-1),5.)) *log(xt)
504 0 : + (737. *pow(xt,4.) -14102. *pow(xt,3.) -28209. *pow(xt,2.) + 610. *xt -508.)/
505 0 : (1296. *pow((xt-1),4.)) );
506 :
507 0 : double E1 = (xt *(18. -11. *xt -pow(xt,2.))/(12.*pow( (1. -xt),3.))
508 0 : + pow(xt,2.)* (15. -16. *xt + 4. *pow(xt,2.))/(6. *pow((1. -xt),4.)) *log(xt)
509 0 : -2./3. *log(xt) );
510 :
511 : double e1 = 4661194./816831.;
512 : double e2 = -8516./2217. ;
513 : double e3 = 0.;
514 : double e4 = 0.;
515 : double e5 = -1.9043;
516 : double e6 = -.1008;
517 : double e7 = .1216;
518 : double e8 = .0183;
519 :
520 : double f1 = -17.3023;
521 : double f2 = 8.5027;
522 : double f3 = 4.5508;
523 : double f4 = .7519;
524 : double f5 = 2.004;
525 : double f6 = .7476;
526 : double f7 = -.5385;
527 : double f8 = .0914;
528 :
529 : double g1 = 14.8088;
530 : double g2 = -10.809;
531 : double g3 = -.874;
532 : double g4 = .4218;
533 : double g5 = -2.9347;
534 : double g6 = .3971;
535 : double g7 = .1600;
536 : double g8 = .0225;
537 :
538 :
539 0 : double c71constmu = ((e1 *_etamu *E1 + f1 + g1 *_etamu) *pow(_etamu,(14./23.))
540 0 : + (e2 *_etamu *E1 + f2 + g2 *_etamu) *pow(_etamu,(16./23.))
541 0 : + (e3 *_etamu *E1 + f3 + g3 *_etamu) *pow(_etamu,(6./23.))
542 0 : + (e4 *_etamu *E1 + f4 + g4 *_etamu) *pow(_etamu,(-12./23.))
543 0 : + (e5 *_etamu *E1 + f5 + g5 *_etamu) *pow(_etamu,.4086)
544 0 : + (e6 *_etamu *E1 + f6 + g6 *_etamu) *pow(_etamu,(-.423))
545 0 : + (e7 *_etamu *E1 + f7 + g7 *_etamu) *pow(_etamu,(-.8994))
546 0 : + (e8 *_etamu *E1 + f8 + g8 *_etamu) *pow(_etamu,.1456 ));
547 :
548 0 : double c71pmu = ( ((297664./14283. *pow(_etamu,(16./23.))
549 0 : -7164416./357075. *pow(_etamu,(14./23.))
550 0 : + 256868./14283. *pow(_etamu,(37./23.)) - 6698884./357075. *pow(_etamu,(39./23.)))
551 0 : *(c8mWsm))
552 0 : + 37208./4761. *(pow(_etamu,(39./23.)) - pow(_etamu,(16./23.))) *(c7mWsm)
553 0 : + c71constmu );
554 :
555 0 : _c71mu = (_alphasmW/_alphasmu *(pow(_etamu,(16./23.))* c7mWsm1 + 8./3. *(pow(_etamu,(14./23.))
556 0 : - pow(_etamu,(16./23.)) ) *c8mWsm1 ) + c71pmu);
557 :
558 0 : _c7emmu = ((32./75. *pow(_etamu,(-9./23.)) - 40./69. *pow(_etamu,(-7./23.)) +
559 0 : 88./575. *pow(_etamu,(16./23.))) *c7mWsm + (-32./575. *pow(_etamu,(-9./23.)) +
560 0 : 32./1449. *pow(_etamu,(-7./23.)) + 640./1449.*pow(_etamu,(14./23.)) -
561 0 : 704./1725.*pow(_etamu,(16./23.)) ) *c8mWsm
562 0 : - 190./8073.*pow(_etamu,(-35./23.)) - 359./3105. *pow(_etamu,(-17./23.)) +
563 0 : 4276./121095. *pow(_etamu,(-12./23.)) + 350531./1009125.*pow(_etamu,(-9./23.))
564 0 : + 2./4347. *pow(_etamu,(-7./23.)) - 5956./15525. *pow(_etamu,(6./23.)) +
565 0 : 38380./169533. *pow(_etamu,(14./23.)) - 748./8625. *pow(_etamu,(16./23.)));
566 :
567 : // Wilson coefficients values as according to Kagan's program
568 : // _c2mu=1.10566;
569 : //_c70mu=-0.314292;
570 : // _c80mu=-0.148954;
571 : // _c71mu=0.480964;
572 : // _c7emmu=0.0323219;
573 :
574 0 : }
575 :
576 : void EvtBtoXsgammaKagan::CalcDelta() {
577 :
578 0 : 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.);
579 :
580 0 : double cDelta27 = ((_alphasmu/(2.*EvtConst::pi))*(_rer2 + _gam27*log(_mb/_mu)) - (_lam2/(9.*_z*pow(_mb,2.))))*_c2mu*_c70mu;
581 :
582 0 : double cDelta78 = (_alphasmu/(2.*EvtConst::pi))*(_rer8 + _gam87*log(_mb/_mu))*_c70mu*_c80mu;
583 :
584 0 : _cDeltatot = cDelta77 + cDelta27 + cDelta78 + (_alphasmu/(2.*EvtConst::pi))*_c71mu*_c70mu + (_alpha/_alphasmu)*(2.*_c7emmu*_c70mu - _kSLemmu*pow(_c70mu,2.));
585 :
586 0 : }
587 :
588 : double EvtBtoXsgammaKagan::Delta(double y, double alphasMu) {
589 :
590 : //Fix for singularity at endpoint
591 0 : if (y >= 1.0) y = 0.9999999999;
592 :
593 0 : return ( - 4.*(alphasMu/(3.*EvtConst::pi*(1. - y)))*(log(1. - y) + 7./4.)*
594 0 : exp(-2.*(alphasMu/(3.*EvtConst::pi))*(pow(log(1. - y),2) + (7./2.)*log(1. - y))));
595 :
596 : }
597 :
598 : double EvtBtoXsgammaKagan::s77(double y) {
599 :
600 : //Fix for singularity at endpoint
601 0 : if (y >= 1.0) y = 0.9999999999;
602 :
603 0 : return ((1./3.)*(7. + y - 2.*pow(y,2) - 2.*(1. + y)*log(1. - y)));
604 : }
605 :
606 : double EvtBtoXsgammaKagan::s88(double y, double mb, double ms) {
607 :
608 : //Fix for singularity at endpoint
609 0 : if (y >= 1.0) y = 0.9999999999;
610 :
611 0 : return ((1./27.)*((2.*(2. - 2.*y + pow(y,2))/y)*(log(1. - y) + 2.*log(mb/ms))
612 0 : - 2.*pow(y,2) - y - 8.*((1. - y)/y)));
613 : }
614 :
615 : double EvtBtoXsgammaKagan::s78(double y) {
616 :
617 : //Fix for singularity at endpoint
618 0 : if (y >= 1.0) y = 0.9999999999;
619 :
620 0 : return ((8./9.)*(((1. - y)/y)*log(1. - y) + 1. + (pow(y,2)/4.)));
621 : }
622 :
623 : double EvtBtoXsgammaKagan::ReG(double y) {
624 :
625 0 : if (y < 4.) return -2.*pow(atan(sqrt(y/(4. - y))),2.);
626 : else {
627 0 : return 2.*(pow(log((sqrt(y) + sqrt(y - 4.))/2.),2.)) - (1./2.)*pow(EvtConst::pi,2.);
628 : }
629 :
630 0 : }
631 :
632 : double EvtBtoXsgammaKagan::ImG(double y) {
633 :
634 0 : if (y < 4.) return 0.0;
635 : else {
636 0 : return (-2.*EvtConst::pi*log((sqrt(y) + sqrt(y - 4.))/2.));
637 : }
638 0 : }
639 :
640 : double EvtBtoXsgammaKagan::s22Func(double y, const std::vector<double> &coeffs) {
641 :
642 : //coeffs[0]=z
643 0 : return (1. - y)*((pow(coeffs[0],2.)/pow(y,2.))*(pow(ReG(y/coeffs[0]),2.) + pow(ImG(y/coeffs[0]),2.)) + (coeffs[0]/y)*ReG(y/coeffs[0]) + (1./4.));
644 :
645 : }
646 :
647 : double EvtBtoXsgammaKagan::s27Func(double y, const std::vector<double> &coeffs) {
648 :
649 : //coeffs[0] = z
650 0 : return (ReG(y/coeffs[0]) + y/(2.*coeffs[0]));
651 :
652 : }
653 :
654 : double EvtBtoXsgammaKagan::DeltaFermiFunc(double y, const std::vector<double> &coeffs1,
655 : const std::vector<double> &coeffs2, const std::vector<double> &coeffs3) {
656 :
657 : //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
658 : //coeffs2[2]=ymH, coeffs3[0]=DeltaCoeff (alphasmu)
659 :
660 0 : return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
661 0 : Delta((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0]);
662 :
663 : }
664 :
665 : double EvtBtoXsgammaKagan::s77FermiFunc(double y, const std::vector<double> &coeffs1,
666 : const std::vector<double> &coeffs2) {
667 :
668 : //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
669 : //coeffs2[2]=ymH
670 0 : return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
671 0 : s77((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y));
672 :
673 : }
674 :
675 : double EvtBtoXsgammaKagan::s88FermiFunc(double y, const std::vector<double> &coeffs1,
676 : const std::vector<double> &coeffs2, const std::vector<double> &coeffs3) {
677 :
678 : //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
679 : //coeffs2[2]=ymH, coeffs3=s88 coeffs
680 0 : return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
681 0 : s88((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y),coeffs3[0], coeffs3[1]);
682 :
683 : }
684 :
685 : double EvtBtoXsgammaKagan::s78FermiFunc(double y, const std::vector<double> &coeffs1,
686 : const std::vector<double> &coeffs2) {
687 :
688 : //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
689 : //coeffs2[2]=ymH
690 0 : return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
691 0 : s78((coeffs2[0]*coeffs2[2])/(coeffs2[1]+y));
692 :
693 : }
694 :
695 : double EvtBtoXsgammaKagan::sFermiFunc(double y, const std::vector<double> &coeffs1,
696 : const std::vector<double> &coeffs2, const std::vector<double> &coeffs3,
697 : const std::vector<double> &coeffs4) {
698 :
699 : //coeffs1=fermi function coeffs, coeffs2[0]=mB, coeffs2[1]=mb,
700 : //coeffs2[2]=ymH, coeffs3[0]=nIntervals in s22 or s27 array, coeffs3[1]=yMin,
701 : //coeffs3[2]=yMax, coeffs4=s22 or s27 array
702 0 : return FermiFunc(y,coeffs1)*(coeffs2[0]/(coeffs2[1]+y))*
703 0 : GetArrayVal(coeffs2[0]*coeffs2[2]/(coeffs2[1]+y), coeffs3[0], coeffs3[1], coeffs3[2], coeffs4);
704 :
705 0 : }
706 :
707 : double EvtBtoXsgammaKagan::Fz(double z) {
708 :
709 0 : return (1. -8.*z + 8.*pow(z,3.) - pow(z,4.) - 12.*pow(z,2.)*log(z));
710 : }
711 :
712 : double EvtBtoXsgammaKagan::GetArrayVal(double xp, double nInterval, double xMin, double xMax, std::vector<double> array) {
713 :
714 0 : double dx = (xMax - xMin)/nInterval;
715 0 : int bin1 = int(((xp-xMin)/(xMax - xMin))*nInterval);
716 :
717 0 : double x1 = double(bin1)*dx + xMin;
718 :
719 0 : if (xp == x1) return array[bin1];
720 :
721 : int bin2(0);
722 0 : if (xp > x1) {
723 0 : bin2 = bin1 + 1;
724 0 : }
725 0 : else if (xp < x1) {
726 0 : bin2 = bin1 - 1;
727 0 : }
728 :
729 0 : if (bin1 <= 0) {
730 : bin1=0;
731 : bin2 = 1;
732 0 : }
733 :
734 : //If xp is in the last bin, always interpolate between the last two bins
735 0 : if (bin1 == (int)nInterval){
736 : bin2 = (int)nInterval;
737 0 : bin1 = (int)nInterval - 1;
738 0 : x1 = double(bin1)*dx + xMin;
739 0 : }
740 :
741 0 : double x2 = double(bin2)*dx + xMin;
742 0 : double y1 = array[bin1];
743 :
744 0 : double y2 = array[bin2];
745 0 : double m = (y2 - y1)/(x2 - x1);
746 0 : double c = y1 - m*x1;
747 0 : double result = m*xp + c;
748 :
749 : return result;
750 :
751 0 : }
752 :
753 : double EvtBtoXsgammaKagan::FermiFunc(double y, const std::vector<double> &coeffs) {
754 :
755 : //Fermi shape functions :1=exponential, 2=gaussian, 3=roman
756 0 : if (int(coeffs[0]) == 1) return EvtBtoXsgammaFermiUtil::FermiExpFunc(y, coeffs);
757 0 : if (int(coeffs[0]) == 2) return EvtBtoXsgammaFermiUtil::FermiGaussFunc(y, coeffs);
758 0 : if (int(coeffs[0]) == 3) return EvtBtoXsgammaFermiUtil::FermiRomanFunc(y, coeffs);
759 0 : return 1.;
760 :
761 0 : }
762 :
763 : double EvtBtoXsgammaKagan::diLogFunc(double y) {
764 :
765 0 : return -log(fabs(1. - y))/y;
766 :
767 : }
768 :
769 :
770 : double EvtBtoXsgammaKagan::diLogMathematica(double y) {
771 :
772 : double li2(0);
773 0 : for(int i=1; i<1000; i++){ //the value 1000 should actually be Infinite...
774 0 : li2+=pow(y,i)/(i*i);
775 : }
776 0 : return li2;
777 : }
|
