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 : // Module: EvtBtoXsll.cc
9 : //
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.
14 : //
15 : // Modification history:
16 : //
17 : // Stephane Willocq Jan 17, 2001 Module created
18 : // Stephane Willocq Jul 15, 2003 Input model parameters
19 : //
20 : //------------------------------------------------------------------------
21 : //
22 : #include "EvtGenBase/EvtPatches.hh"
23 :
24 : #include <stdlib.h>
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"
35 : using std::endl;
36 :
37 0 : EvtBtoXsll::~EvtBtoXsll() {
38 0 : delete _calcprob;
39 0 : }
40 :
41 : std::string EvtBtoXsll::getName(){
42 :
43 0 : return "BTOXSLL";
44 :
45 : }
46 :
47 : EvtDecayBase* EvtBtoXsll::clone(){
48 :
49 0 : return new EvtBtoXsll;
50 :
51 0 : }
52 :
53 :
54 : void EvtBtoXsll::init(){
55 :
56 : // check that there are no arguments
57 :
58 0 : checkNArg(0,4,5);
59 :
60 0 : checkNDaug(3);
61 :
62 : // Check that the two leptons are the same type
63 :
64 0 : EvtId lepton1type = getDaug(1);
65 0 : EvtId lepton2type = getDaug(2);
66 :
67 : int etyp = 0;
68 : int mutyp = 0;
69 : int tautyp = 0;
70 0 : if ( lepton1type == EvtPDL::getId("e+") ||
71 0 : lepton1type == EvtPDL::getId("e-") ) { etyp++;}
72 0 : if ( lepton2type == EvtPDL::getId("e+") ||
73 0 : lepton2type == EvtPDL::getId("e-") ) { etyp++;}
74 0 : if ( lepton1type == EvtPDL::getId("mu+") ||
75 0 : lepton1type == EvtPDL::getId("mu-") ) { mutyp++;}
76 0 : if ( lepton2type == EvtPDL::getId("mu+") ||
77 0 : lepton2type == EvtPDL::getId("mu-") ) { mutyp++;}
78 0 : if ( lepton1type == EvtPDL::getId("tau+") ||
79 0 : lepton1type == EvtPDL::getId("tau-") ) { tautyp++;}
80 0 : if ( lepton2type == EvtPDL::getId("tau+") ||
81 0 : lepton2type == EvtPDL::getId("tau-") ) { tautyp++;}
82 :
83 0 : if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) {
84 :
85 0 : report(Severity::Error,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
86 0 : ::abort();
87 : }
88 :
89 : // Check that the second and third entries are leptons with positive
90 : // and negative charge, respectively
91 :
92 : int lpos = 0;
93 : int lneg = 0;
94 0 : if ( lepton1type == EvtPDL::getId("e+") ||
95 0 : lepton1type == EvtPDL::getId("mu+") ||
96 0 : lepton1type == EvtPDL::getId("tau+") ) { lpos++;}
97 0 : if ( lepton2type == EvtPDL::getId("e-") ||
98 0 : lepton2type == EvtPDL::getId("mu-") ||
99 0 : lepton2type == EvtPDL::getId("tau-") ) { lneg++;}
100 :
101 0 : if ( lpos != 1 || lneg != 1 ) {
102 :
103 0 : report(Severity::Error,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n";
104 0 : ::abort();
105 : }
106 :
107 :
108 0 : _mb=4.8;
109 0 : _ms=0.2;
110 0 : _mq=0.;
111 0 : _pf=0.41;
112 0 : _mxmin=1.1;
113 0 : if ( getNArg()==4)
114 : {
115 : // b-quark mass
116 0 : _mb = getArg(0);
117 : // s-quark mass
118 0 : _ms = getArg(1);
119 : // spectator quark mass
120 0 : _mq = getArg(2);
121 : // Fermi motion parameter
122 0 : _pf = getArg(3);
123 0 : }
124 0 : if ( getNArg()==5)
125 : {
126 0 : _mxmin = getArg(4);
127 0 : }
128 :
129 0 : _calcprob = new EvtBtoXsllUtil;
130 :
131 0 : double ml = EvtPDL::getMeanMass(getDaug(1));
132 :
133 : // determine the maximum probability density from dGdsProb
134 :
135 : int i, j;
136 : int nsteps = 100;
137 : double s = 0.0;
138 0 : double smin = 4.0 * ml * ml;
139 0 : double smax = (_mb - _ms)*(_mb - _ms);
140 : double probMax = -10000.0;
141 : double sProbMax = -10.0;
142 : double uProbMax = -10.0;
143 :
144 0 : for (i=0;i<nsteps;i++)
145 : {
146 0 : s = smin + (i+0.002)*(smax - smin)/(double)nsteps;
147 0 : double prob = _calcprob->dGdsProb(_mb, _ms, ml, s);
148 0 : if (prob > probMax)
149 : {
150 : sProbMax = s;
151 : probMax = prob;
152 0 : }
153 : }
154 :
155 0 : _dGdsProbMax = probMax;
156 :
157 0 : if ( verbose() ) {
158 0 : report(Severity::Info,"EvtGen") << "dGdsProbMax = " << probMax << " for s = " << sProbMax << endl;
159 0 : }
160 :
161 : // determine the maximum probability density from dGdsdupProb
162 :
163 : probMax = -10000.0;
164 : sProbMax = -10.0;
165 :
166 0 : for (i=0;i<nsteps;i++)
167 : {
168 0 : s = smin + (i+0.002)*(smax - smin)/(double)nsteps;
169 0 : double umax = sqrt((s - (_mb + _ms)*(_mb + _ms)) * (s - (_mb - _ms)*(_mb - _ms)));
170 0 : for (j=0;j<nsteps;j++)
171 : {
172 0 : double u = -umax + (j+0.002)*(2.0*umax)/(double)nsteps;
173 0 : double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u);
174 0 : if (prob > probMax)
175 : {
176 : sProbMax = s;
177 : uProbMax = u;
178 : probMax = prob;
179 0 : }
180 : }
181 : }
182 :
183 0 : _dGdsdupProbMax = 2.0*probMax;
184 :
185 0 : if ( verbose() ) {
186 0 : report(Severity::Info,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = " << sProbMax
187 0 : << " and u = " << uProbMax << endl;
188 0 : }
189 :
190 0 : }
191 :
192 : void EvtBtoXsll::initProbMax(){
193 :
194 0 : noProbMax();
195 :
196 0 : }
197 :
198 : void EvtBtoXsll::decay( EvtParticle *p ){
199 :
200 0 : p->makeDaughters(getNDaug(),getDaugs());
201 :
202 0 : EvtParticle* xhadron = p->getDaug(0);
203 0 : EvtParticle* leptonp = p->getDaug(1);
204 0 : EvtParticle* leptonn = p->getDaug(2);
205 :
206 0 : double mass[3];
207 :
208 0 : findMasses( p, getNDaug(), getDaugs(), mass );
209 :
210 0 : double mB = p->mass();
211 0 : double ml = mass[1];
212 : double pb;
213 :
214 : int im = 0;
215 : static int nmsg = 0;
216 : double xhadronMass = -999.0;
217 :
218 0 : EvtVector4R p4xhadron;
219 0 : EvtVector4R p4leptonp;
220 0 : EvtVector4R p4leptonn;
221 :
222 : // require the hadronic system has mass greater than that of a Kaon pion pair
223 :
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 0 : while (xhadronMass < _mxmin)
230 : {
231 0 : im++;
232 :
233 : // Apply Fermi motion and determine effective b-quark mass
234 :
235 : // Old BaBar MC parameters
236 : // double pf = 0.25;
237 : // double ms = 0.2;
238 : // double mq = 0.3;
239 :
240 : double mb = 0.0;
241 :
242 : double xbox, ybox;
243 :
244 0 : while (mb <= 0.0)
245 : {
246 0 : pb = _calcprob->FermiMomentum(_pf);
247 :
248 : // effective b-quark mass
249 0 : mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq);
250 0 : if ( mb>0. && sqrt(mb)-_ms < 2.0*ml ) mb= -10.;
251 : }
252 0 : mb = sqrt(mb);
253 :
254 : // cout << "b-quark momentum = " << pb << " mass = " << mb << endl;
255 :
256 : // generate a dilepton invariant mass
257 :
258 : double s = 0.0;
259 0 : double smin = 4.0 * ml * ml;
260 0 : double smax = (mb - _ms)*(mb - _ms);
261 :
262 0 : while (s == 0.0)
263 : {
264 0 : xbox = EvtRandom::Flat(smin, smax);
265 0 : ybox = EvtRandom::Flat(_dGdsProbMax);
266 0 : double prob= _calcprob->dGdsProb(mb, _ms, ml, xbox);
267 0 : if ( !(prob>=0.0) && !(prob<=0.0)) {
268 : // report(Severity::Info,"EvtGen") << "nan from dGdsProb " << prob << " " << mb << " " << _ms << " " << ml << " " << xbox << std::endl;
269 : }
270 0 : if ( ybox < prob ) s=xbox;
271 : }
272 :
273 : // cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s)
274 : // << " for s = " << s << endl;
275 :
276 : // two-body decay of b quark at rest into s quark and dilepton pair:
277 : // b -> s (ll)
278 :
279 0 : EvtVector4R p4sdilep[2];
280 :
281 0 : double msdilep[2];
282 0 : msdilep[0] = _ms;
283 0 : msdilep[1] = sqrt(s);
284 :
285 0 : EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
286 :
287 : // generate dilepton decay with the expected asymmetry: (ll) -> l+ l-
288 :
289 0 : EvtVector4R p4ll[2];
290 :
291 0 : double mll[2];
292 0 : mll[0] = ml;
293 0 : mll[1] = ml;
294 :
295 : double tmp = 0.0;
296 :
297 0 : while (tmp == 0.0)
298 : {
299 : // (ll) -> l+ l- decay in dilepton rest frame
300 :
301 0 : EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
302 :
303 : // boost to b-quark rest frame
304 :
305 0 : p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
306 0 : p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
307 :
308 : // compute kinematical variable u
309 :
310 0 : EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
311 0 : EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
312 :
313 0 : double u = p4slp.mass2() - p4sln.mass2();
314 :
315 0 : ybox = EvtRandom::Flat(_dGdsdupProbMax);
316 :
317 0 : double prob = _calcprob->dGdsdupProb(mb, _ms, ml, s, u);
318 0 : if ( !(prob>=0.0) && !(prob<=0.0)) {
319 0 : report(Severity::Info,"EvtGen") << "nan from dGdsProb " << prob << " " << mb << " " << _ms << " " << ml << " " << s << " " << u << std::endl;
320 0 : }
321 0 : if (prob > _dGdsdupProbMax && nmsg < 20)
322 : {
323 0 : report(Severity::Info,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob
324 0 : << " " << _dGdsdupProbMax
325 0 : << " for s = " << s << " u = " << u << " mb = " << mb << endl;
326 0 : nmsg++;
327 0 : }
328 0 : if (ybox < prob)
329 : {
330 : tmp = 1.0;
331 : // cout << "dGdsdupProb(s) = " << prob
332 : // << " for u = " << u << endl;
333 0 : }
334 0 : }
335 :
336 :
337 : // assign 4-momenta to valence quarks inside B meson in B rest frame
338 :
339 0 : double phi = EvtRandom::Flat( EvtConst::twoPi );
340 0 : double costh = EvtRandom::Flat( -1.0, 1.0 );
341 0 : double sinth = sqrt(1.0 - costh*costh);
342 :
343 : // b-quark four-momentum in B meson rest frame
344 :
345 0 : EvtVector4R p4b(sqrt(mb*mb + pb*pb),
346 0 : pb*sinth*sin(phi),
347 0 : pb*sinth*cos(phi),
348 0 : pb*costh);
349 :
350 : // B meson in its rest frame
351 : //
352 : // EvtVector4R p4B(mB, 0.0, 0.0, 0.0);
353 : //
354 : // boost B meson to b-quark rest frame
355 : //
356 : // p4B = boostTo(p4B, p4b);
357 : //
358 : // cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl;
359 :
360 : // boost s, l+ and l- to B meson rest frame
361 :
362 : // EvtVector4R p4s = boostTo(p4sdilep[0], p4B);
363 : // p4leptonp = boostTo(p4ll[0], p4B);
364 : // p4leptonn = boostTo(p4ll[1], p4B);
365 :
366 0 : EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
367 0 : p4leptonp = boostTo(p4ll[0], p4b);
368 0 : p4leptonn = boostTo(p4ll[1], p4b);
369 :
370 : // spectator quark in B meson rest frame
371 :
372 0 : EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) );
373 :
374 : // hadron system in B meson rest frame
375 :
376 0 : p4xhadron = p4s + p4q;
377 0 : xhadronMass = p4xhadron.mass();
378 :
379 : // cout << "Xs mass = " << xhadronMass << " trial " << im << endl;
380 0 : }
381 :
382 : // initialize the decay products
383 :
384 0 : xhadron->init(getDaug(0), p4xhadron);
385 :
386 : // For B-bar mesons (i.e. containing a b quark) we have the normal
387 : // order of leptons
388 0 : if ( p->getId() == EvtPDL::getId("anti-B0") ||
389 0 : p->getId() == EvtPDL::getId("B-") )
390 : {
391 0 : leptonp->init(getDaug(1), p4leptonp);
392 0 : leptonn->init(getDaug(2), p4leptonn);
393 0 : }
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
397 : else
398 : {
399 0 : leptonp->init(getDaug(1), p4leptonn);
400 0 : leptonn->init(getDaug(2), p4leptonp);
401 : }
402 :
403 : return ;
404 0 : }
|
