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 : // 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 0 : EvtVubHybrid::EvtVubHybrid()
45 0 : : _noHybrid(false), _storeQplus(true),
46 0 : _mb(4.62), _a(2.27), _alphas(0.22), _dGMax(3.),
47 0 : _nbins_mX(0), _nbins_q2(0), _nbins_El(0), _nbins(0),
48 0 : _masscut(0.28), _bins_mX(0), _bins_q2(0), _bins_El(0),
49 0 : _weights(0), _dGamma(0)
50 0 : {}
51 :
52 0 : EvtVubHybrid::~EvtVubHybrid() {
53 0 : delete _dGamma;
54 0 : delete [] _bins_mX;
55 0 : delete [] _bins_q2;
56 0 : delete [] _bins_El;
57 0 : delete [] _weights;
58 :
59 0 : }
60 :
61 : std::string EvtVubHybrid::getName(){
62 :
63 0 : return "VUBHYBRID";
64 :
65 : }
66 :
67 : EvtDecayBase* EvtVubHybrid::clone(){
68 :
69 0 : return new EvtVubHybrid;
70 :
71 0 : }
72 :
73 : void EvtVubHybrid::init(){
74 :
75 : // check that there are at least 3 arguments
76 0 : if (getNArg() < EvtVubHybrid::nParameters) {
77 0 : report(Severity::Error,"EvtVubHybrid") << "EvtVub generator expected "
78 0 : << "at least " << EvtVubHybrid::nParameters
79 0 : << " arguments but found: " << getNArg()
80 0 : << "\nWill terminate execution!"<<endl;
81 0 : ::abort();
82 0 : } else if (getNArg() == EvtVubHybrid::nParameters) {
83 0 : report(Severity::Warning,"EvtVubHybrid") << "EvtVub: generate B -> Xu l nu events "
84 0 : << "without using the hybrid reweighting."
85 0 : << endl;
86 0 : _noHybrid = true;
87 0 : } else if (getNArg() < EvtVubHybrid::nParameters+EvtVubHybrid::nVariables) {
88 0 : report(Severity::Error,"EvtVubHybrid") << "EvtVub could not read number of bins for "
89 0 : << "all variables used in the reweighting\n"
90 0 : << "Will terminate execution!" << endl;
91 0 : ::abort();
92 : }
93 :
94 : // check that there are 3 daughters
95 0 : checkNDaug(3);
96 :
97 : // read minimum required parameters from decay.dec
98 0 : _mb = getArg(0);
99 0 : _a = getArg(1);
100 0 : _alphas = getArg(2);
101 :
102 : // the maximum dGamma*p2 value depends on alpha_s only:
103 : const double dGMax0 = 3.;
104 0 : _dGMax = 0.21344+8.905*_alphas;
105 0 : 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 0 : static double mB0 = EvtPDL::getMaxMass(EvtPDL::getId("B0"));
112 0 : static double mBP = EvtPDL::getMaxMass(EvtPDL::getId("B+"));
113 0 : static double mB = (mB0<mBP?mB0:mBP);
114 :
115 0 : const double xlow = -_mb;
116 0 : const double xhigh = mB-_mb;
117 : const int aSize = 10000;
118 :
119 0 : EvtPFermi pFermi(_a,mB,_mb);
120 : // pf is the cumulative distribution normalized to 1.
121 0 : _pf.resize(aSize);
122 0 : for(int i=0;i<aSize;i++){
123 0 : double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
124 0 : if ( i== 0 )
125 0 : _pf[i] = pFermi.getFPFermi(kplus);
126 : else
127 0 : _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
128 0 : }
129 0 : for (size_t index=0; index<_pf.size(); index++) {
130 0 : _pf[index]/=_pf[_pf.size()-1];
131 : }
132 :
133 0 : _dGamma = new EvtVubdGamma(_alphas);
134 :
135 0 : if (_noHybrid) return; // Without hybrid weighting, nothing else to do
136 :
137 0 : _nbins_mX = abs((int)getArg(3));
138 0 : _nbins_q2 = abs((int)getArg(4));
139 0 : _nbins_El = abs((int)getArg(5));
140 :
141 : int nextArg = EvtVubHybrid::nParameters + EvtVubHybrid::nVariables;
142 :
143 0 : _nbins = _nbins_mX*_nbins_q2*_nbins_El; // Binning of weight table
144 :
145 0 : int expectArgs = nextArg + _nbins_mX +_nbins_q2 + _nbins_El + _nbins;
146 :
147 0 : if (getNArg() < expectArgs) {
148 0 : report(Severity::Error,"EvtVubHybrid")
149 0 : << " finds " << getNArg() << " arguments, expected " << expectArgs
150 0 : << ". Something is wrong with the tables of weights or thresholds."
151 0 : << "\nWill terminate execution!" << endl;
152 0 : ::abort();
153 : }
154 :
155 : // read bin boundaries from decay.dec
156 : int i;
157 :
158 0 : _bins_mX = new double[_nbins_mX];
159 0 : for (i = 0; i < _nbins_mX; i++,nextArg++) {
160 0 : _bins_mX[i] = getArg(nextArg);
161 : }
162 0 : _masscut = _bins_mX[0];
163 :
164 0 : _bins_q2 = new double[_nbins_q2];
165 0 : for (i = 0; i < _nbins_q2; i++,nextArg++) {
166 0 : _bins_q2[i] = getArg(nextArg);
167 : }
168 :
169 0 : _bins_El = new double[_nbins_El];
170 0 : for (i = 0; i < _nbins_El; i++,nextArg++) {
171 0 : _bins_El[i] = getArg(nextArg);
172 : }
173 :
174 : // read in weights (and rescale to range 0..1)
175 0 : readWeights(nextArg);
176 0 : }
177 :
178 : void EvtVubHybrid::initProbMax(){
179 0 : noProbMax();
180 0 : }
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 0 : EvtVector4R p4;
189 : // R. Faccini 21/02/03
190 : // move the reweighting up , before also shooting the fermi distribution
191 0 : 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 0 : while(rew){
202 :
203 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
204 :
205 0 : xuhad=p->getDaug(0);
206 0 : lepton=p->getDaug(1);
207 0 : neutrino=p->getDaug(2);
208 :
209 0 : mB = p->mass();
210 0 : ml = lepton->mass();
211 :
212 0 : xlow = -_mb;
213 0 : 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 0 : kplus = 2*xhigh;
223 :
224 0 : while( kplus >= xhigh || kplus <= xlow
225 0 : || (_alphas == 0 && kplus >= mB/2-_mb
226 0 : + sqrt(mB*mB/4-_masscut*_masscut))) {
227 0 : kplus = findPFermi(); //_pFermi->shoot();
228 0 : kplus = xlow + kplus*(xhigh-xlow);
229 : }
230 0 : qplus = mB-_mb-kplus;
231 0 : if( (mB-qplus)/2.<=ml) continue;
232 :
233 : int tryit = 1;
234 0 : while (tryit) {
235 :
236 0 : x = EvtRandom::Flat();
237 0 : z = EvtRandom::Flat(0,2);
238 0 : p2=EvtRandom::Flat();
239 0 : p2 = pow(10,lp2epsilon*p2);
240 :
241 0 : El = x*(mB-qplus)/2;
242 0 : if ( El > ml && El < mB/2) {
243 :
244 0 : Eh = z*(mB-qplus)/2+qplus;
245 0 : if ( Eh > 0 && Eh < mB ) {
246 :
247 0 : sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
248 0 : if ( sh > _masscut*_masscut
249 0 : && mB*mB + sh - 2*mB*Eh > ml*ml) {
250 :
251 0 : double xran = EvtRandom::Flat();
252 :
253 0 : double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
254 :
255 0 : if ( y > 1 ) report(Severity::Warning,"EvtVubHybrid") <<"EvtVubHybrid decay probability > 1 found: " << y << endl;
256 0 : if ( y >= xran ) tryit = 0;
257 0 : }
258 : }
259 : }
260 : }
261 :
262 : // compute all kinematic variables needed for reweighting (J. Dingfelder)
263 0 : mX = sqrt(sh);
264 0 : q2 = mB*mB + sh - 2*mB*Eh;
265 :
266 : // Reweighting in bins of mX, q2, El (J. Dingfelder)
267 0 : if (_nbins>0) {
268 0 : double xran1 = EvtRandom::Flat();
269 : double w = 1.0;
270 0 : if (!_noHybrid) w = getWeight(mX, q2, El);
271 0 : if ( w >= xran1 ) rew = false;
272 0 : }
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 0 : double ctH = EvtRandom::Flat(-1,1);
287 0 : double phH = EvtRandom::Flat(0,2*M_PI);
288 0 : 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 0 : sttmp = sqrt(1-ctH*ctH);
296 0 : ptmp = sqrt(Eh*Eh-sh);
297 0 : double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
298 0 : p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
299 0 : xuhad->init( getDaug(0), p4);
300 :
301 0 : 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 0 : xuhad->setLifetime(qplus/10000.);
316 0 : }
317 :
318 : // calculate the W 4 vector in the B Meson restrframe
319 :
320 : double apWB = ptmp;
321 0 : 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 0 : double mW2 = mB*mB + sh - 2*mB*Eh;
327 0 : double beta = ptmp/pWB[0];
328 0 : double gamma = pWB[0]/sqrt(mW2);
329 :
330 0 : double pLW[4];
331 :
332 0 : ptmp = (mW2-ml*ml)/2/sqrt(mW2);
333 0 : pLW[0] = sqrt(ml*ml + ptmp*ptmp);
334 :
335 0 : double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
336 0 : if ( ctL < -1 ) ctL = -1;
337 0 : if ( ctL > 1 ) ctL = 1;
338 0 : sttmp = sqrt(1-ctL*ctL);
339 :
340 : // eX' = eZ x eW
341 0 : double xW[3] = {-pWB[2],pWB[1],0};
342 : // eZ' = eW
343 0 : double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
344 :
345 0 : double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
346 0 : for (j=0;j<2;j++)
347 0 : xW[j] /= lx;
348 :
349 : // eY' = eZ' x eX'
350 0 : double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
351 0 : double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
352 0 : for (j=0;j<3;j++)
353 0 : 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 0 : for (j=0;j<3;j++)
359 0 : pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
360 0 : + sttmp*sin(phL)*ptmp*yW[j]
361 0 : + 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 0 : double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
370 :
371 0 : ptmp = sqrt(El*El-ml*ml);
372 0 : double ctLL = appLB/ptmp;
373 :
374 0 : if ( ctLL > 1 ) ctLL = 1;
375 0 : if ( ctLL < -1 ) ctLL = -1;
376 :
377 0 : double pLB[4] = {El,0,0,0};
378 0 : double pNB[4] = {pWB[0]-El,0,0,0};
379 :
380 0 : for (j=1;j<4;j++) {
381 0 : pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
382 0 : pNB[j] = pWB[j] - pLB[j];
383 : }
384 :
385 0 : p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
386 0 : lepton->init( getDaug(1), p4);
387 :
388 0 : p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
389 0 : neutrino->init( getDaug(2), p4);
390 :
391 : return ;
392 0 : }
393 :
394 :
395 : double EvtVubHybrid::findPFermi() {
396 :
397 0 : double ranNum=EvtRandom::Flat();
398 0 : double oOverBins= 1.0/(float(_pf.size()));
399 : int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
400 0 : int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
401 : int middle;
402 :
403 0 : while (nBinsAbove > nBinsBelow+1) {
404 0 : middle = (nBinsAbove + nBinsBelow+1)>>1;
405 0 : if (ranNum >= _pf[middle]) {
406 : nBinsBelow = middle;
407 0 : } else {
408 : nBinsAbove = middle;
409 : }
410 : }
411 :
412 0 : double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
413 : // binMeasure is always aProbFunc[nBinsBelow],
414 :
415 0 : 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 0 : return (nBinsBelow + .5) * oOverBins;
420 : }
421 :
422 0 : double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
423 :
424 0 : return (nBinsBelow + bFract) * oOverBins;
425 :
426 0 : }
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 0 : for (int i = 0; i < _nbins_mX; i++) {
436 0 : if (mX >= _bins_mX[i]) ibin_mX = i;
437 : }
438 0 : for (int i = 0; i < _nbins_q2; i++) {
439 0 : if (q2 >= _bins_q2[i]) ibin_q2 = i;
440 : }
441 0 : for (int i = 0; i < _nbins_El; i++) {
442 0 : if (El >= _bins_El[i]) ibin_El = i;
443 : }
444 0 : int ibin = ibin_mX + ibin_q2*_nbins_mX + ibin_El*_nbins_mX*_nbins_q2;
445 :
446 0 : if ( (ibin_mX < 0) || (ibin_q2 < 0) || (ibin_El < 0) ) {
447 0 : report(Severity::Error,"EvtVubHybrid") << "Cannot determine hybrid weight "
448 0 : << "for this event "
449 0 : << "-> assign weight = 0" << endl;
450 0 : return 0.0;
451 : }
452 :
453 0 : return _weights[ibin];
454 0 : }
455 :
456 :
457 : void EvtVubHybrid::readWeights(int startArg) {
458 0 : _weights = new double[_nbins];
459 :
460 : double maxw = 0.0;
461 0 : for (int i = 0; i < _nbins; i++, startArg++) {
462 0 : _weights[i] = getArg(startArg);
463 0 : if (_weights[i] > maxw) maxw = _weights[i];
464 : }
465 :
466 0 : if (maxw == 0) {
467 0 : report(Severity::Error,"EvtVubHybrid") << "EvtVub generator expected at least one "
468 0 : << " weight > 0, but found none! "
469 0 : << "Will terminate execution!"<<endl;
470 0 : ::abort();
471 : }
472 :
473 : // rescale weights (to be in range 0..1)
474 0 : for (int i = 0; i < _nbins; i++) {
475 0 : _weights[i] /= maxw;
476 : }
477 0 : }
|
