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: EvtVub.cc
12 : //
13 : // Description: Routine to decay a particle according th phase space
14 : //
15 : // Modification history:
16 : //
17 : // Sven Menke January 17, 2001 Module created
18 : //
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/EvtVub.hh"
28 : #include <string>
29 : #include "EvtGenBase/EvtVector4R.hh"
30 : #include "EvtGenModels/EvtPFermi.hh"
31 : #include "EvtGenModels/EvtVubdGamma.hh"
32 : #include "EvtGenBase/EvtRandom.hh"
33 : using std::endl;
34 :
35 0 : EvtVub::~EvtVub() {
36 0 : if (_dGamma) delete _dGamma;
37 0 : if (_masses) delete [] _masses;
38 0 : if (_weights) delete [] _weights;
39 0 : }
40 :
41 : std::string EvtVub::getName(){
42 :
43 0 : return "VUB";
44 :
45 : }
46 :
47 : EvtDecayBase* EvtVub::clone(){
48 :
49 0 : return new EvtVub;
50 :
51 0 : }
52 :
53 :
54 : void EvtVub::init(){
55 :
56 : // check that there are at least 6 arguments
57 :
58 0 : if (getNArg()<6) {
59 :
60 0 : report(Severity::Error,"EvtGen") << "EvtVub generator expected "
61 0 : << " at least 6 arguments (mb,a,alpha_s,Nbins,m1,w1,...) but found: "
62 0 : <<getNArg()<<endl;
63 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
64 0 : ::abort();
65 :
66 : }
67 :
68 0 : _mb = getArg(0);
69 0 : _a = getArg(1);
70 0 : _alphas = getArg(2);
71 0 : _nbins = abs((int)getArg(3));
72 0 : _storeQplus = (getArg(3)<0?1:0);
73 0 : _masses = new double[_nbins];
74 0 : _weights = new double[_nbins];
75 :
76 0 : if (getNArg()-4 != 2*_nbins) {
77 0 : report(Severity::Error,"EvtGen") << "EvtVub generator expected "
78 0 : << _nbins << " masses and weights but found: "
79 0 : <<(getNArg()-4)/2 <<endl;
80 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
81 0 : ::abort();
82 : }
83 : int i,j = 4;
84 : double maxw = 0.;
85 0 : for (i=0;i<_nbins;i++) {
86 0 : _masses[i] = getArg(j++);
87 0 : if (i>0 && _masses[i] <= _masses[i-1]) {
88 0 : report(Severity::Error,"EvtGen") << "EvtVub generator expected "
89 0 : << " mass bins in ascending order!"
90 0 : << "Will terminate execution!"<<endl;
91 0 : ::abort();
92 : }
93 0 : _weights[i] = getArg(j++);
94 0 : if (_weights[i] < 0) {
95 0 : report(Severity::Error,"EvtGen") << "EvtVub generator expected "
96 0 : << " weights >= 0, but found: "
97 0 : <<_weights[i] <<endl;
98 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
99 0 : ::abort();
100 : }
101 0 : if ( _weights[i] > maxw ) maxw = _weights[i];
102 : }
103 0 : if (maxw == 0) {
104 0 : report(Severity::Error,"EvtGen") << "EvtVub generator expected at least one "
105 0 : << " weight > 0, but found none! "
106 0 : << "Will terminate execution!"<<endl;
107 0 : ::abort();
108 : }
109 0 : for (i=0;i<_nbins;i++) _weights[i]/=maxw;
110 :
111 : // the maximum dGamma*p2 value depends on alpha_s only:
112 :
113 : const double dGMax0 = 3.;
114 0 : _dGMax = 0.21344+8.905*_alphas;
115 0 : if ( _dGMax < dGMax0 ) _dGMax = dGMax0;
116 :
117 : // for the Fermi Motion we need a B-Meson mass - but it's not critical
118 : // to get an exact value; in order to stay in the phase space for
119 : // B+- and B0 use the smaller mass
120 :
121 0 : EvtId BP=EvtPDL::getId("B+");
122 0 : EvtId B0=EvtPDL::getId("B0");
123 :
124 0 : double mB0 = EvtPDL::getMaxMass(B0);
125 0 : double mBP = EvtPDL::getMaxMass(BP);
126 :
127 0 : double mB = (mB0<mBP?mB0:mBP);
128 :
129 0 : const double xlow = -_mb;
130 0 : const double xhigh = mB-_mb;
131 : const int aSize = 10000;
132 :
133 0 : EvtPFermi pFermi(_a,mB,_mb);
134 : // pf is the cumulative distribution
135 : // normalized to 1.
136 0 : _pf.resize(aSize);
137 0 : for(i=0;i<aSize;i++){
138 0 : double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
139 0 : if ( i== 0 )
140 0 : _pf[i] = pFermi.getFPFermi(kplus);
141 : else
142 0 : _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
143 0 : }
144 0 : for (size_t index=0; index<_pf.size(); index++) {
145 0 : _pf[index]/=_pf[_pf.size()-1];
146 : }
147 :
148 : // static EvtHepRandomEngine myEngine;
149 :
150 : // _pFermi = new RandGeneral(myEngine,pf,aSize,0);
151 0 : _dGamma = new EvtVubdGamma(_alphas);
152 :
153 : // check that there are 3 daughters
154 0 : checkNDaug(3);
155 0 : }
156 :
157 : void EvtVub::initProbMax(){
158 :
159 0 : noProbMax();
160 :
161 0 : }
162 :
163 : void EvtVub::decay( EvtParticle *p ){
164 :
165 : int j;
166 : // B+ -> u-bar specflav l+ nu
167 :
168 : EvtParticle *xuhad, *lepton, *neutrino;
169 0 : EvtVector4R p4;
170 : // R. Faccini 21/02/03
171 : // move the reweighting up , before also shooting the fermi distribution
172 0 : double x,z,p2;
173 : double sh=0.0;
174 : double mB,ml,xlow,xhigh,qplus;
175 : double El=0.0;
176 : double Eh=0.0;
177 : double kplus;
178 : const double lp2epsilon=-10;
179 : bool rew(true);
180 0 : while(rew){
181 :
182 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
183 :
184 0 : xuhad=p->getDaug(0);
185 0 : lepton=p->getDaug(1);
186 0 : neutrino=p->getDaug(2);
187 :
188 0 : mB = p->mass();
189 0 : ml = lepton->mass();
190 :
191 0 : xlow = -_mb;
192 0 : xhigh = mB-_mb;
193 :
194 :
195 : // Fermi motion does not need to be computed inside the
196 : // tryit loop as m_b in Gamma0 does not need to be replaced by (m_b+kplus).
197 : // The difference however should be of the Order (lambda/m_b)^2 which is
198 : // beyond the considered orders in the paper anyway ...
199 :
200 : // for alpha_S = 0 and a mass cut on X_u not all values of kplus are
201 : // possible. The maximum value is mB/2-_mb + sqrt(mB^2/4-_masses[0]^2)
202 0 : kplus = 2*xhigh;
203 :
204 0 : while( kplus >= xhigh || kplus <= xlow
205 0 : || (_alphas == 0 && kplus >= mB/2-_mb
206 0 : + sqrt(mB*mB/4-_masses[0]*_masses[0]))) {
207 0 : kplus = findPFermi(); //_pFermi->shoot();
208 0 : kplus = xlow + kplus*(xhigh-xlow);
209 : }
210 0 : qplus = mB-_mb-kplus;
211 0 : if( (mB-qplus)/2.<=ml)continue;
212 :
213 : int tryit = 1;
214 0 : while (tryit) {
215 :
216 0 : x = EvtRandom::Flat();
217 0 : z = EvtRandom::Flat(0,2);
218 0 : p2=EvtRandom::Flat();
219 0 : p2 = pow(10.0,lp2epsilon*p2);
220 :
221 0 : El = x*(mB-qplus)/2;
222 0 : if ( El > ml && El < mB/2) {
223 :
224 0 : Eh = z*(mB-qplus)/2+qplus;
225 0 : if ( Eh > 0 && Eh < mB ) {
226 :
227 0 : sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
228 0 : if ( sh > _masses[0]*_masses[0]
229 0 : && mB*mB + sh - 2*mB*Eh > ml*ml) {
230 :
231 0 : double xran = EvtRandom::Flat();
232 :
233 0 : double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
234 :
235 0 : if ( y > 1 ) report(Severity::Warning,"EvtGen")<<"EvtVub decay probability > 1 found: " << y << endl;
236 0 : if ( y >= xran ) tryit = 0;
237 0 : }
238 : }
239 : }
240 : }
241 : // reweight the Mx distribution
242 0 : if(_nbins>0){
243 0 : double xran1 = EvtRandom::Flat();
244 0 : double m = sqrt(sh);j=0;
245 0 : while ( j < _nbins && m > _masses[j] ) j++;
246 0 : double w = _weights[j-1];
247 0 : if ( w >= xran1 ) rew = false;
248 0 : } else {
249 : rew = false;
250 : }
251 :
252 : }
253 :
254 : // o.k. we have the three kineamtic variables
255 : // now calculate a flat cos Theta_H [-1,1] distribution of the
256 : // hadron flight direction w.r.t the B flight direction
257 : // because the B is a scalar and should decay isotropic.
258 : // Then chose a flat Phi_H [0,2Pi] w.r.t the B flight direction
259 : // and and a flat Phi_L [0,2Pi] in the W restframe w.r.t the
260 : // W flight direction.
261 :
262 0 : double ctH = EvtRandom::Flat(-1,1);
263 0 : double phH = EvtRandom::Flat(0,2*EvtConst::pi);
264 0 : double phL = EvtRandom::Flat(0,2*EvtConst::pi);
265 :
266 : // now compute the four vectors in the B Meson restframe
267 :
268 : double ptmp,sttmp;
269 : // calculate the hadron 4 vector in the B Meson restframe
270 :
271 0 : sttmp = sqrt(1-ctH*ctH);
272 0 : ptmp = sqrt(Eh*Eh-sh);
273 0 : double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
274 0 : p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
275 0 : xuhad->init( getDaug(0), p4);
276 :
277 0 : if (_storeQplus ) {
278 : // cludge to store the hidden parameter q+ with the decay;
279 : // the lifetime of the Xu is abused for this purpose.
280 : // tau = 1 ps corresponds to ctau = 0.3 mm -> in order to
281 : // stay well below BaBars sensitivity we take q+/(10000 GeV) which
282 : // goes up to 0.0005 in the most extreme cases as ctau in mm.
283 : // To extract q+ back from the StdHepTrk its necessary to get
284 : // delta_ctau = Xu->anyDaughter->getVertexTime()-Xu->getVertexTime()
285 : // where these pseudo calls refere to the StdHep time stored at
286 : // the production vertex in the lab for each particle. The boost
287 : // has to be reversed and the result is:
288 : //
289 : // q+ = delta_ctau * 10000 GeV/mm * Mass_Xu/Energy_Xu
290 : //
291 0 : xuhad->setLifetime(qplus/10000.);
292 0 : }
293 :
294 : // calculate the W 4 vector in the B Meson restrframe
295 :
296 : double apWB = ptmp;
297 0 : double pWB[4] = {mB-Eh,-pHB[1],-pHB[2],-pHB[3]};
298 :
299 : // first go in the W restframe and calculate the lepton and
300 : // the neutrino in the W frame
301 :
302 0 : double mW2 = mB*mB + sh - 2*mB*Eh;
303 0 : double beta = ptmp/pWB[0];
304 0 : double gamma = pWB[0]/sqrt(mW2);
305 :
306 0 : double pLW[4];
307 :
308 0 : ptmp = (mW2-ml*ml)/2/sqrt(mW2);
309 0 : pLW[0] = sqrt(ml*ml + ptmp*ptmp);
310 :
311 0 : double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
312 0 : if ( ctL < -1 ) ctL = -1;
313 0 : if ( ctL > 1 ) ctL = 1;
314 0 : sttmp = sqrt(1-ctL*ctL);
315 :
316 : // eX' = eZ x eW
317 0 : double xW[3] = {-pWB[2],pWB[1],0};
318 : // eZ' = eW
319 0 : double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
320 :
321 0 : double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
322 0 : for (j=0;j<2;j++)
323 0 : xW[j] /= lx;
324 :
325 : // eY' = eZ' x eX'
326 0 : double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
327 0 : double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
328 0 : for (j=0;j<3;j++)
329 0 : yW[j] /= ly;
330 :
331 : // p_lep = |p_lep| * ( sin(Theta) * cos(Phi) * eX'
332 : // + sin(Theta) * sin(Phi) * eY'
333 : // + cos(Theta) * eZ')
334 0 : for (j=0;j<3;j++)
335 0 : pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
336 0 : + sttmp*sin(phL)*ptmp*yW[j]
337 0 : + ctL *ptmp*zW[j];
338 :
339 : double apLW = ptmp;
340 :
341 : // boost them back in the B Meson restframe
342 0 : double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
343 :
344 0 : ptmp = sqrt(El*El-ml*ml);
345 0 : double ctLL = appLB/ptmp;
346 :
347 0 : if ( ctLL > 1 ) ctLL = 1;
348 0 : if ( ctLL < -1 ) ctLL = -1;
349 :
350 0 : double pLB[4] = {El,0,0,0};
351 0 : double pNB[4] = {pWB[0]-El,0,0,0};
352 :
353 0 : for (j=1;j<4;j++) {
354 0 : pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
355 0 : pNB[j] = pWB[j] - pLB[j];
356 : }
357 :
358 0 : p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
359 0 : lepton->init( getDaug(1), p4);
360 :
361 0 : p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
362 0 : neutrino->init( getDaug(2), p4);
363 :
364 : return ;
365 0 : }
366 :
367 :
368 : double EvtVub::findPFermi() {
369 :
370 0 : double ranNum=EvtRandom::Flat();
371 0 : double oOverBins= 1.0/(float(_pf.size()));
372 : int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
373 0 : int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
374 : int middle;
375 :
376 0 : while (nBinsAbove > nBinsBelow+1) {
377 0 : middle = (nBinsAbove + nBinsBelow+1)>>1;
378 0 : if (ranNum >= _pf[middle]) {
379 : nBinsBelow = middle;
380 0 : } else {
381 : nBinsAbove = middle;
382 : }
383 : }
384 :
385 0 : double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
386 : // binMeasure is always aProbFunc[nBinsBelow],
387 :
388 0 : if ( bSize == 0 ) {
389 : // rand lies right in a bin of measure 0. Simply return the center
390 : // of the range of that bin. (Any value between k/N and (k+1)/N is
391 : // equally good, in this rare case.)
392 0 : return (nBinsBelow + .5) * oOverBins;
393 : }
394 :
395 0 : double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
396 :
397 0 : return (nBinsBelow + bFract) * oOverBins;
398 :
399 0 : }
|
