Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtVub.cpp
CommitLineData
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: 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"
33using std::endl;
34
35EvtVub::~EvtVub() {
36 if (_dGamma) delete _dGamma;
37 if (_masses) delete [] _masses;
38 if (_weights) delete [] _weights;
39}
40
41std::string EvtVub::getName(){
42
43 return "VUB";
44
45}
46
47EvtDecayBase* EvtVub::clone(){
48
49 return new EvtVub;
50
51}
52
53
54void EvtVub::init(){
55
56 // check that there are at least 6 arguments
57
58 if (getNArg()<6) {
59
60 report(ERROR,"EvtGen") << "EvtVub generator expected "
61 << " at least 6 arguments (mb,a,alpha_s,Nbins,m1,w1,...) but found: "
62 <<getNArg()<<endl;
63 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
64 ::abort();
65
66 }
67
68 _mb = getArg(0);
69 _a = getArg(1);
70 _alphas = getArg(2);
71 _nbins = abs((int)getArg(3));
72 _storeQplus = (getArg(3)<0?1:0);
73 _masses = new double[_nbins];
74 _weights = new double[_nbins];
75
76 if (getNArg()-4 != 2*_nbins) {
77 report(ERROR,"EvtGen") << "EvtVub generator expected "
78 << _nbins << " masses and weights but found: "
79 <<(getNArg()-4)/2 <<endl;
80 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
81 ::abort();
82 }
83 int i,j = 4;
84 double maxw = 0.;
85 for (i=0;i<_nbins;i++) {
86 _masses[i] = getArg(j++);
87 if (i>0 && _masses[i] <= _masses[i-1]) {
88 report(ERROR,"EvtGen") << "EvtVub generator expected "
89 << " mass bins in ascending order!"
90 << "Will terminate execution!"<<endl;
91 ::abort();
92 }
93 _weights[i] = getArg(j++);
94 if (_weights[i] < 0) {
95 report(ERROR,"EvtGen") << "EvtVub generator expected "
96 << " weights >= 0, but found: "
97 <<_weights[i] <<endl;
98 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
99 ::abort();
100 }
101 if ( _weights[i] > maxw ) maxw = _weights[i];
102 }
103 if (maxw == 0) {
104 report(ERROR,"EvtGen") << "EvtVub generator expected at least one "
105 << " weight > 0, but found none! "
106 << "Will terminate execution!"<<endl;
107 ::abort();
108 }
109 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 _dGMax = 0.21344+8.905*_alphas;
115 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 EvtId BP=EvtPDL::getId("B+");
122 EvtId B0=EvtPDL::getId("B0");
123
124 double mB0 = EvtPDL::getMaxMass(B0);
125 double mBP = EvtPDL::getMaxMass(BP);
126
127 double mB = (mB0<mBP?mB0:mBP);
128
129 const double xlow = -_mb;
130 const double xhigh = mB-_mb;
131 const int aSize = 10000;
132
133 EvtPFermi pFermi(_a,mB,_mb);
134 // pf is the cumulative distribution
135 // normalized to 1.
136 _pf.resize(aSize);
137 for(i=0;i<aSize;i++){
138 double kplus = xlow + (double)(i+0.5)/((double)aSize)*(xhigh-xlow);
139 if ( i== 0 )
140 _pf[i] = pFermi.getFPFermi(kplus);
141 else
142 _pf[i] = _pf[i-1] + pFermi.getFPFermi(kplus);
143 }
144 for (size_t index=0; index<_pf.size(); index++) {
145 _pf[index]/=_pf[_pf.size()-1];
146 }
147
148 // static EvtHepRandomEngine myEngine;
149
150 // _pFermi = new RandGeneral(myEngine,pf,aSize,0);
151 _dGamma = new EvtVubdGamma(_alphas);
152
153 // check that there are 3 daughters
154 checkNDaug(3);
155}
156
157void EvtVub::initProbMax(){
158
159 noProbMax();
160
161}
162
163void EvtVub::decay( EvtParticle *p ){
164
165 int j;
166 // B+ -> u-bar specflav l+ nu
167
168 EvtParticle *xuhad, *lepton, *neutrino;
169 EvtVector4R p4;
170 // R. Faccini 21/02/03
171 // move the reweighting up , before also shooting the fermi distribution
172 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 while(rew){
181
182 p->initializePhaseSpace(getNDaug(),getDaugs());
183
184 xuhad=p->getDaug(0);
185 lepton=p->getDaug(1);
186 neutrino=p->getDaug(2);
187
188 mB = p->mass();
189 ml = lepton->mass();
190
191 xlow = -_mb;
192 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 kplus = 2*xhigh;
203
204 while( kplus >= xhigh || kplus <= xlow
205 || (_alphas == 0 && kplus >= mB/2-_mb
206 + sqrt(mB*mB/4-_masses[0]*_masses[0]))) {
207 kplus = findPFermi(); //_pFermi->shoot();
208 kplus = xlow + kplus*(xhigh-xlow);
209 }
210 qplus = mB-_mb-kplus;
211 if( (mB-qplus)/2.<=ml)continue;
212
213 int tryit = 1;
214 while (tryit) {
215
216 x = EvtRandom::Flat();
217 z = EvtRandom::Flat(0,2);
218 p2=EvtRandom::Flat();
219 p2 = pow(10.0,lp2epsilon*p2);
220
221 El = x*(mB-qplus)/2;
222 if ( El > ml && El < mB/2) {
223
224 Eh = z*(mB-qplus)/2+qplus;
225 if ( Eh > 0 && Eh < mB ) {
226
227 sh = p2*pow(mB-qplus,2)+2*qplus*(Eh-qplus)+qplus*qplus;
228 if ( sh > _masses[0]*_masses[0]
229 && mB*mB + sh - 2*mB*Eh > ml*ml) {
230
231 double xran = EvtRandom::Flat();
232
233 double y = _dGamma->getdGdxdzdp(x,z,p2)/_dGMax*p2;
234
235 if ( y > 1 ) report(WARNING,"EvtGen")<<"EvtVub decay probability > 1 found: " << y << endl;
236 if ( y >= xran ) tryit = 0;
237 }
238 }
239 }
240 }
241 // reweight the Mx distribution
242 if(_nbins>0){
243 double xran1 = EvtRandom::Flat();
244 double m = sqrt(sh);j=0;
245 while ( j < _nbins && m > _masses[j] ) j++;
246 double w = _weights[j-1];
247 if ( w >= xran1 ) rew = false;
248 } 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 double ctH = EvtRandom::Flat(-1,1);
263 double phH = EvtRandom::Flat(0,2*EvtConst::pi);
264 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 sttmp = sqrt(1-ctH*ctH);
272 ptmp = sqrt(Eh*Eh-sh);
273 double pHB[4] = {Eh,ptmp*sttmp*cos(phH),ptmp*sttmp*sin(phH),ptmp*ctH};
274 p4.set(pHB[0],pHB[1],pHB[2],pHB[3]);
275 xuhad->init( getDaug(0), p4);
276
277 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 xuhad->setLifetime(qplus/10000.);
292 }
293
294 // calculate the W 4 vector in the B Meson restrframe
295
296 double apWB = ptmp;
297 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 double mW2 = mB*mB + sh - 2*mB*Eh;
303 double beta = ptmp/pWB[0];
304 double gamma = pWB[0]/sqrt(mW2);
305
306 double pLW[4];
307
308 ptmp = (mW2-ml*ml)/2/sqrt(mW2);
309 pLW[0] = sqrt(ml*ml + ptmp*ptmp);
310
311 double ctL = (El - gamma*pLW[0])/beta/gamma/ptmp;
312 if ( ctL < -1 ) ctL = -1;
313 if ( ctL > 1 ) ctL = 1;
314 sttmp = sqrt(1-ctL*ctL);
315
316 // eX' = eZ x eW
317 double xW[3] = {-pWB[2],pWB[1],0};
318 // eZ' = eW
319 double zW[3] = {pWB[1]/apWB,pWB[2]/apWB,pWB[3]/apWB};
320
321 double lx = sqrt(xW[0]*xW[0]+xW[1]*xW[1]);
322 for (j=0;j<2;j++)
323 xW[j] /= lx;
324
325 // eY' = eZ' x eX'
326 double yW[3] = {-pWB[1]*pWB[3],-pWB[2]*pWB[3],pWB[1]*pWB[1]+pWB[2]*pWB[2]};
327 double ly = sqrt(yW[0]*yW[0]+yW[1]*yW[1]+yW[2]*yW[2]);
328 for (j=0;j<3;j++)
329 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 for (j=0;j<3;j++)
335 pLW[j+1] = sttmp*cos(phL)*ptmp*xW[j]
336 + sttmp*sin(phL)*ptmp*yW[j]
337 + ctL *ptmp*zW[j];
338
339 double apLW = ptmp;
340
341 // boost them back in the B Meson restframe
342 double appLB = beta*gamma*pLW[0] + gamma*ctL*apLW;
343
344 ptmp = sqrt(El*El-ml*ml);
345 double ctLL = appLB/ptmp;
346
347 if ( ctLL > 1 ) ctLL = 1;
348 if ( ctLL < -1 ) ctLL = -1;
349
350 double pLB[4] = {El,0,0,0};
351 double pNB[4] = {pWB[0]-El,0,0,0};
352
353 for (j=1;j<4;j++) {
354 pLB[j] = pLW[j] + (ctLL*ptmp - ctL*apLW)/apWB*pWB[j];
355 pNB[j] = pWB[j] - pLB[j];
356 }
357
358 p4.set(pLB[0],pLB[1],pLB[2],pLB[3]);
359 lepton->init( getDaug(1), p4);
360
361 p4.set(pNB[0],pNB[1],pNB[2],pNB[3]);
362 neutrino->init( getDaug(2), p4);
363
364 return ;
365}
366
367
368double EvtVub::findPFermi() {
369
370 double ranNum=EvtRandom::Flat();
371 double oOverBins= 1.0/(float(_pf.size()));
372 int nBinsBelow = 0; // largest k such that I[k] is known to be <= rand
373 int nBinsAbove = _pf.size(); // largest k such that I[k] is known to be > rand
374 int middle;
375
376 while (nBinsAbove > nBinsBelow+1) {
377 middle = (nBinsAbove + nBinsBelow+1)>>1;
378 if (ranNum >= _pf[middle]) {
379 nBinsBelow = middle;
380 } else {
381 nBinsAbove = middle;
382 }
383 }
384
385 double bSize = _pf[nBinsAbove] - _pf[nBinsBelow];
386 // binMeasure is always aProbFunc[nBinsBelow],
387
388 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 return (nBinsBelow + .5) * oOverBins;
393 }
394
395 double bFract = (ranNum - _pf[nBinsBelow]) / bSize;
396
397 return (nBinsBelow + bFract) * oOverBins;
398
399}