]>
Commit | Line | Data |
---|---|---|
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 | // 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 | EvtBtoXsll::~EvtBtoXsll() { | |
38 | delete _calcprob; | |
39 | } | |
40 | ||
41 | std::string EvtBtoXsll::getName(){ | |
42 | ||
43 | return "BTOXSLL"; | |
44 | ||
45 | } | |
46 | ||
47 | EvtDecayBase* EvtBtoXsll::clone(){ | |
48 | ||
49 | return new EvtBtoXsll; | |
50 | ||
51 | } | |
52 | ||
53 | ||
54 | void EvtBtoXsll::init(){ | |
55 | ||
56 | // check that there are no arguments | |
57 | ||
58 | checkNArg(0,4,5); | |
59 | ||
60 | checkNDaug(3); | |
61 | ||
62 | // Check that the two leptons are the same type | |
63 | ||
64 | EvtId lepton1type = getDaug(1); | |
65 | EvtId lepton2type = getDaug(2); | |
66 | ||
67 | int etyp = 0; | |
68 | int mutyp = 0; | |
69 | int tautyp = 0; | |
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++;} | |
82 | ||
83 | if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) { | |
84 | ||
85 | report(ERROR,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n"; | |
86 | ::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 | 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++;} | |
100 | ||
101 | if ( lpos != 1 || lneg != 1 ) { | |
102 | ||
103 | report(ERROR,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n"; | |
104 | ::abort(); | |
105 | } | |
106 | ||
107 | ||
108 | _mb=4.8; | |
109 | _ms=0.2; | |
110 | _mq=0.; | |
111 | _pf=0.41; | |
112 | _mxmin=1.1; | |
113 | if ( getNArg()==4) | |
114 | { | |
115 | // b-quark mass | |
116 | _mb = getArg(0); | |
117 | // s-quark mass | |
118 | _ms = getArg(1); | |
119 | // spectator quark mass | |
120 | _mq = getArg(2); | |
121 | // Fermi motion parameter | |
122 | _pf = getArg(3); | |
123 | } | |
124 | if ( getNArg()==5) | |
125 | { | |
126 | _mxmin = getArg(4); | |
127 | } | |
128 | ||
129 | _calcprob = new EvtBtoXsllUtil; | |
130 | ||
131 | 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 | 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; | |
143 | ||
144 | for (i=0;i<nsteps;i++) | |
145 | { | |
146 | s = smin + (i+0.002)*(smax - smin)/(double)nsteps; | |
147 | double prob = _calcprob->dGdsProb(_mb, _ms, ml, s); | |
148 | if (prob > probMax) | |
149 | { | |
150 | sProbMax = s; | |
151 | probMax = prob; | |
152 | } | |
153 | } | |
154 | ||
155 | _dGdsProbMax = probMax; | |
156 | ||
157 | if ( verbose() ) { | |
158 | report(INFO,"EvtGen") << "dGdsProbMax = " << probMax << " for s = " << sProbMax << endl; | |
159 | } | |
160 | ||
161 | // determine the maximum probability density from dGdsdupProb | |
162 | ||
163 | probMax = -10000.0; | |
164 | sProbMax = -10.0; | |
165 | ||
166 | for (i=0;i<nsteps;i++) | |
167 | { | |
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++) | |
171 | { | |
172 | double u = -umax + (j+0.002)*(2.0*umax)/(double)nsteps; | |
173 | double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u); | |
174 | if (prob > probMax) | |
175 | { | |
176 | sProbMax = s; | |
177 | uProbMax = u; | |
178 | probMax = prob; | |
179 | } | |
180 | } | |
181 | } | |
182 | ||
183 | _dGdsdupProbMax = 2.0*probMax; | |
184 | ||
185 | if ( verbose() ) { | |
186 | report(INFO,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = " << sProbMax | |
187 | << " and u = " << uProbMax << endl; | |
188 | } | |
189 | ||
190 | } | |
191 | ||
192 | void EvtBtoXsll::initProbMax(){ | |
193 | ||
194 | noProbMax(); | |
195 | ||
196 | } | |
197 | ||
198 | void EvtBtoXsll::decay( EvtParticle *p ){ | |
199 | ||
200 | p->makeDaughters(getNDaug(),getDaugs()); | |
201 | ||
202 | EvtParticle* xhadron = p->getDaug(0); | |
203 | EvtParticle* leptonp = p->getDaug(1); | |
204 | EvtParticle* leptonn = p->getDaug(2); | |
205 | ||
206 | double mass[3]; | |
207 | ||
208 | findMasses( p, getNDaug(), getDaugs(), mass ); | |
209 | ||
210 | double mB = p->mass(); | |
211 | double ml = mass[1]; | |
212 | double pb; | |
213 | ||
214 | int im = 0; | |
215 | static int nmsg = 0; | |
216 | double xhadronMass = -999.0; | |
217 | ||
218 | EvtVector4R p4xhadron; | |
219 | EvtVector4R p4leptonp; | |
220 | 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 | while (xhadronMass < _mxmin) | |
230 | { | |
231 | 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 | while (mb <= 0.0) | |
245 | { | |
246 | pb = _calcprob->FermiMomentum(_pf); | |
247 | ||
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.; | |
251 | } | |
252 | 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 | double smin = 4.0 * ml * ml; | |
260 | double smax = (mb - _ms)*(mb - _ms); | |
261 | ||
262 | while (s == 0.0) | |
263 | { | |
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; | |
269 | } | |
270 | 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 | EvtVector4R p4sdilep[2]; | |
280 | ||
281 | double msdilep[2]; | |
282 | msdilep[0] = _ms; | |
283 | msdilep[1] = sqrt(s); | |
284 | ||
285 | EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb); | |
286 | ||
287 | // generate dilepton decay with the expected asymmetry: (ll) -> l+ l- | |
288 | ||
289 | EvtVector4R p4ll[2]; | |
290 | ||
291 | double mll[2]; | |
292 | mll[0] = ml; | |
293 | mll[1] = ml; | |
294 | ||
295 | double tmp = 0.0; | |
296 | ||
297 | while (tmp == 0.0) | |
298 | { | |
299 | // (ll) -> l+ l- decay in dilepton rest frame | |
300 | ||
301 | EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]); | |
302 | ||
303 | // boost to b-quark rest frame | |
304 | ||
305 | p4ll[0] = boostTo(p4ll[0], p4sdilep[1]); | |
306 | p4ll[1] = boostTo(p4ll[1], p4sdilep[1]); | |
307 | ||
308 | // compute kinematical variable u | |
309 | ||
310 | EvtVector4R p4slp = p4sdilep[0] + p4ll[0]; | |
311 | EvtVector4R p4sln = p4sdilep[0] + p4ll[1]; | |
312 | ||
313 | double u = p4slp.mass2() - p4sln.mass2(); | |
314 | ||
315 | ybox = EvtRandom::Flat(_dGdsdupProbMax); | |
316 | ||
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; | |
320 | } | |
321 | if (prob > _dGdsdupProbMax && nmsg < 20) | |
322 | { | |
323 | report(INFO,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob | |
324 | << " " << _dGdsdupProbMax | |
325 | << " for s = " << s << " u = " << u << " mb = " << mb << endl; | |
326 | nmsg++; | |
327 | } | |
328 | if (ybox < prob) | |
329 | { | |
330 | tmp = 1.0; | |
331 | // cout << "dGdsdupProb(s) = " << prob | |
332 | // << " for u = " << u << endl; | |
333 | } | |
334 | } | |
335 | ||
336 | ||
337 | // assign 4-momenta to valence quarks inside B meson in B rest frame | |
338 | ||
339 | double phi = EvtRandom::Flat( EvtConst::twoPi ); | |
340 | double costh = EvtRandom::Flat( -1.0, 1.0 ); | |
341 | double sinth = sqrt(1.0 - costh*costh); | |
342 | ||
343 | // b-quark four-momentum in B meson rest frame | |
344 | ||
345 | EvtVector4R p4b(sqrt(mb*mb + pb*pb), | |
346 | pb*sinth*sin(phi), | |
347 | pb*sinth*cos(phi), | |
348 | 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 | EvtVector4R p4s = boostTo(p4sdilep[0], p4b); | |
367 | p4leptonp = boostTo(p4ll[0], p4b); | |
368 | p4leptonn = boostTo(p4ll[1], p4b); | |
369 | ||
370 | // spectator quark in B meson rest frame | |
371 | ||
372 | 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 | p4xhadron = p4s + p4q; | |
377 | xhadronMass = p4xhadron.mass(); | |
378 | ||
379 | // cout << "Xs mass = " << xhadronMass << " trial " << im << endl; | |
380 | } | |
381 | ||
382 | // initialize the decay products | |
383 | ||
384 | 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 | if ( p->getId() == EvtPDL::getId("anti-B0") || | |
389 | p->getId() == EvtPDL::getId("B-") ) | |
390 | { | |
391 | leptonp->init(getDaug(1), p4leptonp); | |
392 | leptonn->init(getDaug(2), p4leptonn); | |
393 | } | |
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 | leptonp->init(getDaug(1), p4leptonn); | |
400 | leptonn->init(getDaug(2), p4leptonp); | |
401 | } | |
402 | ||
403 | return ; | |
404 | } |