1 //--------------------------------------------------------------------------
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.
8 // Module: EvtBtoXsll.cc
10 // Description: Routine to generate non-resonant B -> Xs l+ l- decays.
11 // It generates a dilepton mass spectrum according to Kruger and Sehgal
12 // and then generates the two lepton momenta accoring to Ali et al.
13 // The resultant X_s particles may be decayed by JETSET.
15 // Modification history:
17 // Stephane Willocq Jan 17, 2001 Module created
18 // Stephane Willocq Jul 15, 2003 Input model parameters
20 //------------------------------------------------------------------------
22 #include "EvtGenBase/EvtPatches.hh"
25 #include "EvtGenBase/EvtRandom.hh"
26 #include "EvtGenBase/EvtParticle.hh"
27 #include "EvtGenBase/EvtGenKine.hh"
28 #include "EvtGenBase/EvtPDL.hh"
29 #include "EvtGenBase/EvtReport.hh"
30 #include "EvtGenModels/EvtbTosllAmp.hh"
31 #include "EvtGenModels/EvtBtoXsll.hh"
32 #include "EvtGenModels/EvtBtoXsllUtil.hh"
33 #include "EvtGenBase/EvtConst.hh"
34 #include "EvtGenBase/EvtId.hh"
37 EvtBtoXsll::~EvtBtoXsll() {
41 std::string EvtBtoXsll::getName(){
47 EvtDecayBase* EvtBtoXsll::clone(){
49 return new EvtBtoXsll;
54 void EvtBtoXsll::init(){
56 // check that there are no arguments
62 // Check that the two leptons are the same type
64 EvtId lepton1type = getDaug(1);
65 EvtId lepton2type = getDaug(2);
70 if ( lepton1type == EvtPDL::getId("e+") ||
71 lepton1type == EvtPDL::getId("e-") ) { etyp++;}
72 if ( lepton2type == EvtPDL::getId("e+") ||
73 lepton2type == EvtPDL::getId("e-") ) { etyp++;}
74 if ( lepton1type == EvtPDL::getId("mu+") ||
75 lepton1type == EvtPDL::getId("mu-") ) { mutyp++;}
76 if ( lepton2type == EvtPDL::getId("mu+") ||
77 lepton2type == EvtPDL::getId("mu-") ) { mutyp++;}
78 if ( lepton1type == EvtPDL::getId("tau+") ||
79 lepton1type == EvtPDL::getId("tau-") ) { tautyp++;}
80 if ( lepton2type == EvtPDL::getId("tau+") ||
81 lepton2type == EvtPDL::getId("tau-") ) { tautyp++;}
83 if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) {
85 report(ERROR,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
89 // Check that the second and third entries are leptons with positive
90 // and negative charge, respectively
94 if ( lepton1type == EvtPDL::getId("e+") ||
95 lepton1type == EvtPDL::getId("mu+") ||
96 lepton1type == EvtPDL::getId("tau+") ) { lpos++;}
97 if ( lepton2type == EvtPDL::getId("e-") ||
98 lepton2type == EvtPDL::getId("mu-") ||
99 lepton2type == EvtPDL::getId("tau-") ) { lneg++;}
101 if ( lpos != 1 || lneg != 1 ) {
103 report(ERROR,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n";
119 // spectator quark mass
121 // Fermi motion parameter
129 _calcprob = new EvtBtoXsllUtil;
131 double ml = EvtPDL::getMeanMass(getDaug(1));
133 // determine the maximum probability density from dGdsProb
138 double smin = 4.0 * ml * ml;
139 double smax = (_mb - _ms)*(_mb - _ms);
140 double probMax = -10000.0;
141 double sProbMax = -10.0;
142 double uProbMax = -10.0;
144 for (i=0;i<nsteps;i++)
146 s = smin + (i+0.002)*(smax - smin)/(double)nsteps;
147 double prob = _calcprob->dGdsProb(_mb, _ms, ml, s);
155 _dGdsProbMax = probMax;
158 report(INFO,"EvtGen") << "dGdsProbMax = " << probMax << " for s = " << sProbMax << endl;
161 // determine the maximum probability density from dGdsdupProb
166 for (i=0;i<nsteps;i++)
168 s = smin + (i+0.002)*(smax - smin)/(double)nsteps;
169 double umax = sqrt((s - (_mb + _ms)*(_mb + _ms)) * (s - (_mb - _ms)*(_mb - _ms)));
170 for (j=0;j<nsteps;j++)
172 double u = -umax + (j+0.002)*(2.0*umax)/(double)nsteps;
173 double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u);
183 _dGdsdupProbMax = 2.0*probMax;
186 report(INFO,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = " << sProbMax
187 << " and u = " << uProbMax << endl;
192 void EvtBtoXsll::initProbMax(){
198 void EvtBtoXsll::decay( EvtParticle *p ){
200 p->makeDaughters(getNDaug(),getDaugs());
202 EvtParticle* xhadron = p->getDaug(0);
203 EvtParticle* leptonp = p->getDaug(1);
204 EvtParticle* leptonn = p->getDaug(2);
208 findMasses( p, getNDaug(), getDaugs(), mass );
210 double mB = p->mass();
216 double xhadronMass = -999.0;
218 EvtVector4R p4xhadron;
219 EvtVector4R p4leptonp;
220 EvtVector4R p4leptonn;
222 // require the hadronic system has mass greater than that of a Kaon pion pair
224 // while (xhadronMass < 0.6333)
225 // the above minimum value of K+pi mass appears to be too close
226 // to threshold as far as JETSET is concerned
227 // (JETSET gets caught in an infinite loop)
228 // so we choose a lightly larger value for the threshold
229 while (xhadronMass < _mxmin)
233 // Apply Fermi motion and determine effective b-quark mass
235 // Old BaBar MC parameters
246 pb = _calcprob->FermiMomentum(_pf);
248 // effective b-quark mass
249 mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq);
250 if ( mb>0. && sqrt(mb)-_ms < 2.0*ml ) mb= -10.;
254 // cout << "b-quark momentum = " << pb << " mass = " << mb << endl;
256 // generate a dilepton invariant mass
259 double smin = 4.0 * ml * ml;
260 double smax = (mb - _ms)*(mb - _ms);
264 xbox = EvtRandom::Flat(smin, smax);
265 ybox = EvtRandom::Flat(_dGdsProbMax);
266 double prob= _calcprob->dGdsProb(mb, _ms, ml, xbox);
267 if ( !(prob>=0.0) && !(prob<=0.0)) {
268 // report(INFO,"EvtGen") << "nan from dGdsProb " << prob << " " << mb << " " << _ms << " " << ml << " " << xbox << std::endl;
270 if ( ybox < prob ) s=xbox;
273 // cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s)
274 // << " for s = " << s << endl;
276 // two-body decay of b quark at rest into s quark and dilepton pair:
279 EvtVector4R p4sdilep[2];
283 msdilep[1] = sqrt(s);
285 EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
287 // generate dilepton decay with the expected asymmetry: (ll) -> l+ l-
299 // (ll) -> l+ l- decay in dilepton rest frame
301 EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
303 // boost to b-quark rest frame
305 p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
306 p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
308 // compute kinematical variable u
310 EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
311 EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
313 double u = p4slp.mass2() - p4sln.mass2();
315 ybox = EvtRandom::Flat(_dGdsdupProbMax);
317 double prob = _calcprob->dGdsdupProb(mb, _ms, ml, s, u);
318 if ( !(prob>=0.0) && !(prob<=0.0)) {
319 report(INFO,"EvtGen") << "nan from dGdsProb " << prob << " " << mb << " " << _ms << " " << ml << " " << s << " " << u << std::endl;
321 if (prob > _dGdsdupProbMax && nmsg < 20)
323 report(INFO,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob
324 << " " << _dGdsdupProbMax
325 << " for s = " << s << " u = " << u << " mb = " << mb << endl;
331 // cout << "dGdsdupProb(s) = " << prob
332 // << " for u = " << u << endl;
337 // assign 4-momenta to valence quarks inside B meson in B rest frame
339 double phi = EvtRandom::Flat( EvtConst::twoPi );
340 double costh = EvtRandom::Flat( -1.0, 1.0 );
341 double sinth = sqrt(1.0 - costh*costh);
343 // b-quark four-momentum in B meson rest frame
345 EvtVector4R p4b(sqrt(mb*mb + pb*pb),
350 // B meson in its rest frame
352 // EvtVector4R p4B(mB, 0.0, 0.0, 0.0);
354 // boost B meson to b-quark rest frame
356 // p4B = boostTo(p4B, p4b);
358 // cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl;
360 // boost s, l+ and l- to B meson rest frame
362 // EvtVector4R p4s = boostTo(p4sdilep[0], p4B);
363 // p4leptonp = boostTo(p4ll[0], p4B);
364 // p4leptonn = boostTo(p4ll[1], p4B);
366 EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
367 p4leptonp = boostTo(p4ll[0], p4b);
368 p4leptonn = boostTo(p4ll[1], p4b);
370 // spectator quark in B meson rest frame
372 EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) );
374 // hadron system in B meson rest frame
376 p4xhadron = p4s + p4q;
377 xhadronMass = p4xhadron.mass();
379 // cout << "Xs mass = " << xhadronMass << " trial " << im << endl;
382 // initialize the decay products
384 xhadron->init(getDaug(0), p4xhadron);
386 // For B-bar mesons (i.e. containing a b quark) we have the normal
388 if ( p->getId() == EvtPDL::getId("anti-B0") ||
389 p->getId() == EvtPDL::getId("B-") )
391 leptonp->init(getDaug(1), p4leptonp);
392 leptonn->init(getDaug(2), p4leptonn);
394 // For B mesons (i.e. containing a b-bar quark) we need to flip the
395 // role of the positive and negative leptons in order to produce the
396 // correct forward-backward asymmetry between the two leptons
399 leptonp->init(getDaug(1), p4leptonn);
400 leptonn->init(getDaug(2), p4leptonp);