]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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 | // Copyright Information: | |
9 | // Copyright (C) 1998 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtVubHybrid.cc | |
12 | // | |
13 | // Description: Routine to decay a particle according to phase space. | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // Jochen Dingfelder February 1, 2005 Created Module as update of the | |
18 | // original module EvtVub by Sven Menke | |
19 | //--------------------------------------------------------------------------- | |
20 | // | |
21 | #include "EvtGenBase/EvtPatches.hh" | |
22 | #include <stdlib.h> | |
23 | #include "EvtGenBase/EvtParticle.hh" | |
24 | #include "EvtGenBase/EvtGenKine.hh" | |
25 | #include "EvtGenBase/EvtPDL.hh" | |
26 | #include "EvtGenBase/EvtReport.hh" | |
27 | #include "EvtGenModels/EvtVubHybrid.hh" | |
28 | #include "EvtGenBase/EvtVector4R.hh" | |
29 | #include "EvtGenModels/EvtPFermi.hh" | |
30 | #include "EvtGenModels/EvtVubdGamma.hh" | |
31 | #include "EvtGenBase/EvtRandom.hh" | |
32 | ||
33 | #include <string> | |
34 | #include <iostream> | |
35 | #include <fstream> | |
36 | using std::ifstream; | |
37 | using std::cout; | |
38 | using std::endl; | |
39 | ||
40 | ||
41 | // _noHybrid will be set TRUE if the DECAY.DEC file has no binning or weights | |
42 | // _storeQplus should alwasy be TRUE: writes out Fermi motion parameter | |
43 | ||
44 | EvtVubHybrid::EvtVubHybrid() | |
45 | : _noHybrid(false), _storeQplus(true), | |
46 | _mb(4.62), _a(2.27), _alphas(0.22), _dGMax(3.), | |
47 | _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0), | |
48 | _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0), | |
49 | _weights(0), _dGamma(0) | |
50 | {} | |
51 | ||
52 | EvtVubHybrid::~EvtVubHybrid() { | |
53 | delete _dGamma; | |
54 | delete [] _bins_mX; | |
55 | delete [] _bins_q2; | |
56 | delete [] _bins_El; | |
57 | delete [] _weights; | |
58 | ||
59 | } | |
60 | ||
61 | std::string EvtVubHybrid::getName(){ | |
62 | ||
63 | return "VUBHYBRID"; | |
64 | ||
65 | } | |
66 | ||
67 | EvtDecayBase* EvtVubHybrid::clone(){ | |
68 | ||
69 | return new EvtVubHybrid; | |
70 | ||
71 | } | |
72 | ||
73 | void EvtVubHybrid::init(){ | |
74 | ||
75 | // check that there are at least 3 arguments | |
76 | if (getNArg() < EvtVubHybrid::nParameters) { | |
77 | report(ERROR,"EvtVubHybrid") << "EvtVub generator expected " | |
78 | << "at least " << EvtVubHybrid::nParameters | |
79 | << " arguments but found: " << getNArg() | |
80 | << "\nWill terminate execution!"<<endl; | |
81 | ::abort(); | |
82 | } else if (getNArg() == EvtVubHybrid::nParameters) { | |
83 | report(WARNING,"EvtVubHybrid") << "EvtVub: generate B -> Xu l nu events " | |
84 | << "without using the hybrid reweighting." | |
85 | << endl; | |
86 | _noHybrid = true; | |
87 | } else if (getNArg() < EvtVubHybrid::nParameters+EvtVubHybrid::nVariables) { | |
88 | report(ERROR,"EvtVubHybrid") << "EvtVub could not read number of bins for " | |
89 | << "all variables used in the reweighting\n" | |
90 | << "Will terminate execution!" << endl; | |
91 | ::abort(); | |
92 | } | |
93 | ||
94 | // check that there are 3 daughters | |
95 | checkNDaug(3); | |
96 | ||
97 | // read minimum required parameters from decay.dec | |
98 | _mb = getArg(0); | |
99 | _a = getArg(1); | |
100 | _alphas = getArg(2); | |
101 | ||
102 | // the maximum dGamma*p2 value depends on alpha_s only: | |
103 | const double dGMax0 = 3.; | |
104 | _dGMax = 0.21344+8.905*_alphas; | |
105 | if ( _dGMax < dGMax0 ) _dGMax = dGMax0; | |
106 | ||
107 | // for the Fermi Motion we need a B-Meson mass - but it's not critical | |
108 | // to get an exact value; in order to stay in the phase space for | |
109 | // B+- and B0 use the smaller mass | |
110 | ||
111 | static double mB0 = EvtPDL::getMaxMass(EvtPDL::getId("B0")); | |
112 | static double mBP = EvtPDL::getMaxMass(EvtPDL::getId("B+")); | |
113 | static double mB = (mB0<mBP?mB0:mBP); | |
114 | ||
115 | const double xlow = -_mb; | |
116 | const double xhigh = mB-_mb; | |
117 | const int aSize = 10000; | |
118 | ||
119 | EvtPFermi pFermi(_a,mB,_mb); | |
120 | // pf is the cumulative distribution normalized to 1. | |
121 | _pf.resize(aSize); | |
122 | for(int i=0;i<aSize;i++){ | |
123 | double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow); | |
124 | if ( i== 0 ) | |
125 | _pf[i] = pFermi.getFPFermi(kplus); | |
126 | else | |
127 | _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus); | |
128 | } | |
129 | for (size_t index=0; index<_pf.size(); index++) { | |
130 | _pf[index]/=_pf[_pf.size()-1]; | |
131 | } | |
132 | ||
133 | _dGamma = new EvtVubdGamma(_alphas); | |
134 | ||
135 | if (_noHybrid) return; // Without hybrid weighting, nothing else to do | |
136 | ||
137 | _nbins_mX = abs((int)getArg(3)); | |
138 | _nbins_q2 = abs((int)getArg(4)); | |
139 | _nbins_El = abs((int)getArg(5)); | |
140 | ||
141 | int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables; | |
142 | ||
143 | _nbins = _nbins_mX*_nbins_q2*_nbins_El; // Binning of weight table | |
144 | ||
145 | int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins; | |
146 | ||
147 | if (getNArg() < expectArgs) { | |
148 | report(ERROR,"EvtVubHybrid") | |
149 | << " finds " << getNArg() << " arguments, expected " << expectArgs | |
150 | << ". Something is wrong with the tables of weights or thresholds." | |
151 | << "\nWill terminate execution!" << endl; | |
152 | ::abort(); | |
153 | } | |
154 | ||
155 | // read bin boundaries from decay.dec | |
156 | int i; | |
157 | ||
158 | _bins_mX = new double[_nbins_mX]; | |
159 | for (i = 0; i < _nbins_mX; i++,nextArg++) { | |
160 | _bins_mX[i] = getArg(nextArg); | |
161 | } | |
162 | _masscut = _bins_mX[0]; | |
163 | ||
164 | _bins_q2 = new double[_nbins_q2]; | |
165 | for (i = 0; i < _nbins_q2; i++,nextArg++) { | |
166 | _bins_q2[i] = getArg(nextArg); | |
167 | } | |
168 | ||
169 | _bins_El = new double[_nbins_El]; | |
170 | for (i = 0; i < _nbins_El; i++,nextArg++) { | |
171 | _bins_El[i] = getArg(nextArg); | |
172 | } | |
173 | ||
174 | // read in weights (and rescale to range 0..1) | |
175 | readWeights(nextArg); | |
176 | } | |
177 | ||
178 | void EvtVubHybrid::initProbMax(){ | |
179 | noProbMax(); | |
180 | } | |
181 | ||
182 | void EvtVubHybrid::decay( EvtParticle *p ){ | |
183 | ||
184 | int j; | |
185 | // B+ -> u-bar specflav l+ nu | |
186 | ||
187 | EvtParticle *xuhad, *lepton, *neutrino; | |
188 | EvtVector4R p4; | |
189 | // R. Faccini 21/02/03 | |
190 | // move the reweighting up , before also shooting the fermi distribution | |
191 | double x,z,p2; | |
192 | double sh=0.0; | |
193 | double mB,ml,xlow,xhigh,qplus; | |
194 | double El=0.0; | |
195 | double Eh=0.0; | |
196 | double kplus; | |
197 | double q2, mX; | |
198 | ||
199 | const double lp2epsilon=-10; | |
200 | bool rew(true); | |
201 | while(rew){ | |
202 | ||
203 | p->initializePhaseSpace(getNDaug(),getDaugs()); | |
204 | ||
205 | xuhad=p->getDaug(0); | |
206 | lepton=p->getDaug(1); | |
207 | neutrino=p->getDaug(2); | |
208 | ||
209 | mB = p->mass(); | |
210 | ml = lepton->mass(); | |
211 | ||
212 | xlow = -_mb; | |
213 | xhigh = mB-_mb; | |
214 | ||
215 | // Fermi motion does not need to be computed inside the | |
216 | // tryit loop as m_b in Gamma0 does not need to be replaced by (m_b+kplus). | |
217 | // The difference however should be of the Order (lambda/m_b)^2 which is | |
218 | // beyond the considered orders in the paper anyway ... | |
219 | ||
220 | // for alpha_S = 0 and a mass cut on X_u not all values of kplus are | |
221 | // possible. The maximum value is mB/2-_mb + sqrt(mB^2/4-_masscut^2) | |
222 | kplus = 2*xhigh; | |
223 | ||
224 | while( kplus >= xhigh || kplus <= xlow | |
225 | || (_alphas == 0 && kplus >= mB/2-_mb | |
226 | + sqrt(mB*mB/4-_masscut*_masscut))) { | |
227 | kplus = findPFermi(); //_pFermi->shoot(); | |
228 | kplus = xlow + kplus*(xhigh-xlow); | |
229 | } | |
230 | qplus = mB-_mb-kplus; | |
231 | if( (mB-qplus)/2.<=ml) continue; | |
232 | ||
233 | int tryit = 1; | |
234 | while (tryit) { | |
235 | ||
236 | x = EvtRandom::Flat(); | |
237 | z = EvtRandom::Flat(0,2); | |
238 | p2=EvtRandom::Flat(); | |
239 | p2 = pow(10,lp2epsilon*p2); | |
240 | ||
241 | El = x*(mB-qplus)/2; | |
242 | if ( El > ml && El < mB/2) { | |
243 | ||
244 | Eh = z*(mB-qplus)/2+qplus; | |
245 | if ( Eh > 0 && Eh < mB ) { | |
246 | ||
247 | sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus; | |
248 | if ( sh > _masscut*_masscut | |
249 | && mB*mB + sh - 2*mB*Eh > ml*ml) { | |
250 | ||
251 | double xran = EvtRandom::Flat(); | |
252 | ||
253 | double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2; | |
254 | ||
255 | if ( y > 1 ) report(WARNING,"EvtVubHybrid") <<"EvtVubHybrid decay probability > 1 found: " << y << endl; | |
256 | if ( y >= xran ) tryit = 0; | |
257 | } | |
258 | } | |
259 | } | |
260 | } | |
261 | ||
262 | // compute all kinematic variables needed for reweighting (J. Dingfelder) | |
263 | mX = sqrt(sh); | |
264 | q2 = mB*mB + sh - 2*mB*Eh; | |
265 | ||
266 | // Reweighting in bins of mX, q2, El (J. Dingfelder) | |
267 | if (_nbins>0) { | |
268 | double xran1 = EvtRandom::Flat(); | |
269 | double w = 1.0; | |
270 | if (!_noHybrid) w = getWeight(mX, q2, El); | |
271 | if ( w >= xran1 ) rew = false; | |
272 | } | |
273 | else { | |
274 | rew = false; | |
275 | } | |
276 | } | |
277 | ||
278 | // o.k. we have the three kineamtic variables | |
279 | // now calculate a flat cos Theta_H [-1,1] distribution of the | |
280 | // hadron flight direction w.r.t the B flight direction | |
281 | // because the B is a scalar and should decay isotropic. | |
282 | // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction | |
283 | // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the | |
284 | // W flight direction. | |
285 | ||
286 | double ctH = EvtRandom::Flat(-1,1); | |
287 | double phH = EvtRandom::Flat(0,2*M_PI); | |
288 | double phL = EvtRandom::Flat(0,2*M_PI); | |
289 | ||
290 | // now compute the four vectors in the B Meson restframe | |
291 | ||
292 | double ptmp,sttmp; | |
293 | // calculate the hadron 4 vector in the B Meson restframe | |
294 | ||
295 | sttmp = sqrt(1-ctH*ctH); | |
296 | ptmp = sqrt(Eh*Eh-sh); | |
297 | double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH}; | |
298 | p4.set(pHB[0],pHB[1],pHB[2],pHB[3]); | |
299 | xuhad->init( getDaug(0), p4); | |
300 | ||
301 | if (_storeQplus ) { | |
302 | // cludge to store the hidden parameter q+ with the decay; | |
303 | // the lifetime of the Xu is abused for this purpose. | |
304 | // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to | |
305 | // stay well below BaBars sensitivity we take q+/(10000 GeV) which | |
306 | // goes up to 0.0005 in the most extreme cases as ctau in mm. | |
307 | // To extract q+ back from the StdHepTrk its necessary to get | |
308 | // delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime() | |
309 | // where these pseudo calls refere to the StdHep time stored at | |
310 | // the production vertex in the lab for each particle. The boost | |
311 | // has to be reversed and the result is: | |
312 | // | |
313 | // q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu | |
314 | // | |
315 | xuhad->setLifetime(qplus/10000.); | |
316 | } | |
317 | ||
318 | // calculate the W 4 vector in the B Meson restrframe | |
319 | ||
320 | double apWB = ptmp; | |
321 | double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]}; | |
322 | ||
323 | // first go in the W restframe and calculate the lepton and | |
324 | // the neutrino in the W frame | |
325 | ||
326 | double mW2 = mB*mB + sh - 2*mB*Eh; | |
327 | double beta = ptmp/pWB[0]; | |
328 | double gamma = pWB[0]/sqrt(mW2); | |
329 | ||
330 | double pLW[4]; | |
331 | ||
332 | ptmp = (mW2-ml*ml)/2/sqrt(mW2); | |
333 | pLW[0] = sqrt(ml*ml + ptmp*ptmp); | |
334 | ||
335 | double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp; | |
336 | if ( ctL < -1 ) ctL = -1; | |
337 | if ( ctL > 1 ) ctL = 1; | |
338 | sttmp = sqrt(1-ctL*ctL); | |
339 | ||
340 | // eX' = eZ x eW | |
341 | double xW[3] = {-pWB[2],pWB[1],0}; | |
342 | // eZ' = eW | |
343 | double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB}; | |
344 | ||
345 | double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]); | |
346 | for (j=0;j<2;j++) | |
347 | xW[j] /= lx; | |
348 | ||
349 | // eY' = eZ' x eX' | |
350 | double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]}; | |
351 | double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]); | |
352 | for (j=0;j<3;j++) | |
353 | yW[j] /= ly; | |
354 | ||
355 | // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX' | |
356 | // + sin(Theta) * sin(Phi) * eY' | |
357 | // + cos(Theta) * eZ') | |
358 | for (j=0;j<3;j++) | |
359 | pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j] | |
360 | + sttmp*sin(phL)*ptmp*yW[j] | |
361 | + ctL *ptmp*zW[j]; | |
362 | ||
363 | double apLW = ptmp; | |
364 | // calculate the neutrino 4 vector in the W restframe | |
365 | //double pNW[4] = {sqrt(mW2)-pLW[0],-pLW[1],-pLW[2],-pLW[3]}; | |
366 | ||
367 | // boost them back in the B Meson restframe | |
368 | ||
369 | double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW; | |
370 | ||
371 | ptmp = sqrt(El*El-ml*ml); | |
372 | double ctLL = appLB/ptmp; | |
373 | ||
374 | if ( ctLL > 1 ) ctLL = 1; | |
375 | if ( ctLL < -1 ) ctLL = -1; | |
376 | ||
377 | double pLB[4] = {El,0,0,0}; | |
378 | double pNB[4] = {pWB[0]-El,0,0,0}; | |
379 | ||
380 | for (j=1;j<4;j++) { | |
381 | pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j]; | |
382 | pNB[j] = pWB[j] - pLB[j]; | |
383 | } | |
384 | ||
385 | p4.set(pLB[0],pLB[1],pLB[2],pLB[3]); | |
386 | lepton->init( getDaug(1), p4); | |
387 | ||
388 | p4.set(pNB[0],pNB[1],pNB[2],pNB[3]); | |
389 | neutrino->init( getDaug(2), p4); | |
390 | ||
391 | return ; | |
392 | } | |
393 | ||
394 | ||
395 | double EvtVubHybrid::findPFermi() { | |
396 | ||
397 | double ranNum=EvtRandom::Flat(); | |
398 | double oOverBins= 1.0/(float(_pf.size())); | |
399 | int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand | |
400 | int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand | |
401 | int middle; | |
402 | ||
403 | while (nBinsAbove > nBinsBelow+1) { | |
404 | middle = (nBinsAbove + nBinsBelow+1)>>1; | |
405 | if (ranNum >= _pf[middle]) { | |
406 | nBinsBelow = middle; | |
407 | } else { | |
408 | nBinsAbove = middle; | |
409 | } | |
410 | } | |
411 | ||
412 | double bSize = _pf[nBinsAbove] - _pf[nBinsBelow]; | |
413 | // binMeasure is always aProbFunc[nBinsBelow], | |
414 | ||
415 | if ( bSize == 0 ) { | |
416 | // rand lies right in a bin of measure 0. Simply return the center | |
417 | // of the range of that bin. (Any value between k/N and (k+1)/N is | |
418 | // equally good, in this rare case.) | |
419 | return (nBinsBelow + .5) * oOverBins; | |
420 | } | |
421 | ||
422 | double bFract = (ranNum - _pf[nBinsBelow]) / bSize; | |
423 | ||
424 | return (nBinsBelow + bFract) * oOverBins; | |
425 | ||
426 | } | |
427 | ||
428 | ||
429 | double EvtVubHybrid::getWeight(double mX, double q2, double El) { | |
430 | ||
431 | int ibin_mX = -1; | |
432 | int ibin_q2 = -1; | |
433 | int ibin_El = -1; | |
434 | ||
435 | for (int i = 0; i < _nbins_mX; i++) { | |
436 | if (mX >= _bins_mX[i]) ibin_mX = i; | |
437 | } | |
438 | for (int i = 0; i < _nbins_q2; i++) { | |
439 | if (q2 >= _bins_q2[i]) ibin_q2 = i; | |
440 | } | |
441 | for (int i = 0; i < _nbins_El; i++) { | |
442 | if (El >= _bins_El[i]) ibin_El = i; | |
443 | } | |
444 | int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2; | |
445 | ||
446 | if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) { | |
447 | report(ERROR,"EvtVubHybrid") << "Cannot determine hybrid weight " | |
448 | << "for this event " | |
449 | << "-> assign weight = 0" << endl; | |
450 | return 0.0; | |
451 | } | |
452 | ||
453 | return _weights[ibin]; | |
454 | } | |
455 | ||
456 | ||
457 | void EvtVubHybrid::readWeights(int startArg) { | |
458 | _weights = new double[_nbins]; | |
459 | ||
460 | double maxw = 0.0; | |
461 | for (int i = 0; i < _nbins; i++, startArg++) { | |
462 | _weights[i] = getArg(startArg); | |
463 | if (_weights[i] > maxw) maxw = _weights[i]; | |
464 | } | |
465 | ||
466 | if (maxw == 0) { | |
467 | report(ERROR,"EvtVubHybrid") << "EvtVub generator expected at least one " | |
468 | << " weight > 0, but found none! " | |
469 | << "Will terminate execution!"<<endl; | |
470 | ::abort(); | |
471 | } | |
472 | ||
473 | // rescale weights (to be in range 0..1) | |
474 | for (int i = 0; i < _nbins; i++) { | |
475 | _weights[i] /= maxw; | |
476 | } | |
477 | } |