CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtBtoXsll Class Reference

#include <EvtBtoXsll.hh>

+ Inheritance diagram for EvtBtoXsll:

Public Member Functions

 EvtBtoXsll ()
 
virtual ~EvtBtoXsll ()
 
void getName (std::string &name)
 
EvtDecayBaseclone ()
 
void initProbMax ()
 
void init ()
 
void decay (EvtParticle *p)
 
- Public Member Functions inherited from EvtDecayIncoherent
void makeDecay (EvtParticle *p)
 
virtual ~EvtDecayIncoherent ()
 
void setDaughterSpinDensity (int daughter)
 
int isDaughterSpinDensitySet (int daughter)
 
- Public Member Functions inherited from EvtDecayBase
virtual void getName (std::string &name)=0
 
virtual void decay (EvtParticle *p)=0
 
virtual void makeDecay (EvtParticle *p)=0
 
virtual EvtDecayBaseclone ()=0
 
virtual void init ()
 
virtual void initProbMax ()
 
virtual std::string commandName ()
 
virtual void command (std::string cmd)
 
double getProbMax (double prob)
 
double resetProbMax (double prob)
 
 EvtDecayBase ()
 
virtual ~EvtDecayBase ()
 
virtual bool matchingDecay (const EvtDecayBase &other) const
 
EvtId getParentId ()
 
double getBranchingFraction ()
 
void disableCheckQ ()
 
void checkQ ()
 
int getNDaug ()
 
EvtIdgetDaugs ()
 
EvtId getDaug (int i)
 
int getNArg ()
 
int getPHOTOS ()
 
void setPHOTOS ()
 
void setVerbose ()
 
void setSummary ()
 
double * getArgs ()
 
std::string * getArgsStr ()
 
double getArg (int j)
 
std::string getArgStr (int j)
 
std::string getModelName ()
 
int getDSum ()
 
int summary ()
 
int verbose ()
 
void saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
 
void printSummary ()
 
void setProbMax (double prbmx)
 
void noProbMax ()
 
void checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
 
void checkNDaug (int d1, int d2=-1)
 
void checkSpinParent (EvtSpinType::spintype sp)
 
void checkSpinDaughter (int d1, EvtSpinType::spintype sp)
 
virtual int nRealDaughters ()
 

Additional Inherited Members

- Static Public Member Functions inherited from EvtDecayBase
static void findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
 
static void findMass (EvtParticle *p)
 
static double findMaxMass (EvtParticle *p)
 
- Protected Member Functions inherited from EvtDecayBase
bool daugsDecayedByParentModel ()
 
- Protected Attributes inherited from EvtDecayBase
bool _daugsDecayedByParentModel
 

Detailed Description

Definition at line 28 of file EvtBtoXsll.hh.

Constructor & Destructor Documentation

◆ EvtBtoXsll()

EvtBtoXsll::EvtBtoXsll ( )
inline

Definition at line 32 of file EvtBtoXsll.hh.

32{}

Referenced by clone().

◆ ~EvtBtoXsll()

EvtBtoXsll::~EvtBtoXsll ( )
virtual

Definition at line 37 of file EvtBtoXsll.cc.

37{}

Member Function Documentation

◆ clone()

EvtDecayBase * EvtBtoXsll::clone ( )
virtual

Implements EvtDecayBase.

Definition at line 45 of file EvtBtoXsll.cc.

45 {
46
47 return new EvtBtoXsll;
48
49}

◆ decay()

void EvtBtoXsll::decay ( EvtParticle p)
virtual

Implements EvtDecayBase.

Definition at line 196 of file EvtBtoXsll.cc.

196 {
197
199
200 EvtParticle* xhadron = p->getDaug(0);
201 EvtParticle* leptonp = p->getDaug(1);
202 EvtParticle* leptonn = p->getDaug(2);
203
204 double mass[3];
205
206 findMasses( p, getNDaug(), getDaugs(), mass );
207
208 double mB = p->mass();
209 double ml = mass[1];
210 double pb;
211
212 int im = 0;
213 static int nmsg = 0;
214 double xhadronMass = -999.0;
215
216 EvtVector4R p4xhadron;
217 EvtVector4R p4leptonp;
218 EvtVector4R p4leptonn;
219
220 // require the hadronic system has mass greater than that of a Kaon pion pair
221
222 // while (xhadronMass < 0.6333)
223 // the above minimum value of K+pi mass appears to be too close
224 // to threshold as far as JETSET is concerned
225 // (JETSET gets caught in an infinite loop)
226 // so we choose a lightly larger value for the threshold
227 while (xhadronMass < _mxmin)
228 {
229 im++;
230
231 // Apply Fermi motion and determine effective b-quark mass
232
233 // Old BaBar MC parameters
234 // double pf = 0.25;
235 // double ms = 0.2;
236 // double mq = 0.3;
237
238 double mb = 0.0;
239
240 double xbox, ybox;
241
242 while (mb <= 0.0)
243 {
244 pb = _calcprob->FermiMomentum(_pf);
245
246 // effective b-quark mass
247 mb = mB*mB + _mq*_mq - 2.0*mB*sqrt(pb*pb + _mq*_mq);
248 }
249 mb = sqrt(mb);
250
251 // cout << "b-quark momentum = " << pb << " mass = " << mb << endl;
252
253 // generate a dilepton invariant mass
254
255 double s = 0.0;
256 double smin = 4.0 * ml * ml;
257 double smax = (mb - _ms)*(mb - _ms);
258
259 while (s == 0.0)
260 {
261 xbox = EvtRandom::Flat(smin, smax);
262 ybox = EvtRandom::Flat(_dGdsProbMax);
263 if (ybox < _calcprob->dGdsProb(mb, _ms, ml, xbox)) { s = xbox;}
264 }
265
266 // cout << "dGdsProb(s) = " << _calcprob->dGdsProb(mb, _ms, ml, s)
267 // << " for s = " << s << endl;
268
269 // two-body decay of b quark at rest into s quark and dilepton pair:
270 // b -> s (ll)
271
272 EvtVector4R p4sdilep[2];
273
274 double msdilep[2];
275 msdilep[0] = _ms;
276 msdilep[1] = sqrt(s);
277
278 EvtGenKine::PhaseSpace(2, msdilep, p4sdilep, mb);
279
280 // generate dilepton decay with the expected asymmetry: (ll) -> l+ l-
281
282 EvtVector4R p4ll[2];
283
284 double mll[2];
285 mll[0] = ml;
286 mll[1] = ml;
287
288 double tmp = 0.0;
289
290 while (tmp == 0.0)
291 {
292 // (ll) -> l+ l- decay in dilepton rest frame
293
294 EvtGenKine::PhaseSpace(2, mll, p4ll, msdilep[1]);
295
296 // boost to b-quark rest frame
297
298 p4ll[0] = boostTo(p4ll[0], p4sdilep[1]);
299 p4ll[1] = boostTo(p4ll[1], p4sdilep[1]);
300
301 // compute kinematical variable u
302
303 EvtVector4R p4slp = p4sdilep[0] + p4ll[0];
304 EvtVector4R p4sln = p4sdilep[0] + p4ll[1];
305
306 double u = p4slp.mass2() - p4sln.mass2();
307
308 ybox = EvtRandom::Flat(_dGdsdupProbMax);
309
310 double prob = _calcprob->dGdsdupProb(mb, _ms, ml, s, u);
311 if (prob > _dGdsdupProbMax && nmsg < 20)
312 {
313 report(INFO,"EvtGen") << "d2gdsdup GT d2gdsdup_max:" << prob
314 << " " << _dGdsdupProbMax
315 << " for s = " << s << " u = " << u << " mb = " << mb << endl;
316 nmsg++;
317 }
318 if (ybox < prob)
319 {
320 tmp = 1.0;
321 // cout << "dGdsdupProb(s) = " << prob
322 // << " for u = " << u << endl;
323 }
324 }
325
326 // assign 4-momenta to valence quarks inside B meson in B rest frame
327
328 double phi = EvtRandom::Flat( EvtConst::twoPi );
329 double costh = EvtRandom::Flat( -1.0, 1.0 );
330 double sinth = sqrt(1.0 - costh*costh);
331
332 // b-quark four-momentum in B meson rest frame
333
334 EvtVector4R p4b(sqrt(mb*mb + pb*pb),
335 pb*sinth*sin(phi),
336 pb*sinth*cos(phi),
337 pb*costh);
338
339 // B meson in its rest frame
340 //
341 // EvtVector4R p4B(mB, 0.0, 0.0, 0.0);
342 //
343 // boost B meson to b-quark rest frame
344 //
345 // p4B = boostTo(p4B, p4b);
346 //
347 // cout << " B meson mass in b-quark rest frame = " << p4B.mass() << endl;
348
349 // boost s, l+ and l- to B meson rest frame
350
351 // EvtVector4R p4s = boostTo(p4sdilep[0], p4B);
352 // p4leptonp = boostTo(p4ll[0], p4B);
353 // p4leptonn = boostTo(p4ll[1], p4B);
354
355 EvtVector4R p4s = boostTo(p4sdilep[0], p4b);
356 p4leptonp = boostTo(p4ll[0], p4b);
357 p4leptonn = boostTo(p4ll[1], p4b);
358
359 // spectator quark in B meson rest frame
360
361 EvtVector4R p4q( sqrt(pb*pb + _mq*_mq), -p4b.get(1), -p4b.get(2), -p4b.get(3) );
362
363 // hadron system in B meson rest frame
364
365 p4xhadron = p4s + p4q;
366 xhadronMass = p4xhadron.mass();
367
368 // cout << "Xs mass = " << xhadronMass << " trial " << im << endl;
369 }
370
371 // initialize the decay products
372
373 xhadron->init(getDaug(0), p4xhadron);
374
375 // For B-bar mesons (i.e. containing a b quark) we have the normal
376 // order of leptons
377 if ( p->getId() == EvtPDL::getId("anti-B0") ||
378 p->getId() == EvtPDL::getId("B-") )
379 {
380 leptonp->init(getDaug(1), p4leptonp);
381 leptonn->init(getDaug(2), p4leptonn);
382 }
383 // For B mesons (i.e. containing a b-bar quark) we need to flip the
384 // role of the positive and negative leptons in order to produce the
385 // correct forward-backward asymmetry between the two leptons
386 else
387 {
388 leptonp->init(getDaug(1), p4leptonn);
389 leptonn->init(getDaug(2), p4leptonp);
390 }
391
392 return ;
393}
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
double mass
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
@ INFO
Definition: EvtReport.hh:52
XmlRpcServer s
Definition: HelloServer.cpp:11
double FermiMomentum(double pf)
double dGdsdupProb(double mb, double ms, double ml, double s, double u)
static const double twoPi
Definition: EvtConst.hh:29
static void findMasses(EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
static double PhaseSpace(int ndaug, double mass[30], EvtVector4R p4[30], double mp)
Definition: EvtGenKine.cc:50
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
double mass() const
Definition: EvtParticle.cc:127
static double Flat()
Definition: EvtRandom.cc:73
double mass() const
Definition: EvtVector4R.cc:39
double mass2() const
Definition: EvtVector4R.hh:116

◆ getName()

void EvtBtoXsll::getName ( std::string &  name)
virtual

Implements EvtDecayBase.

Definition at line 39 of file EvtBtoXsll.cc.

39 {
40
41 model_name="BTOXSLL";
42
43}

◆ init()

void EvtBtoXsll::init ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 52 of file EvtBtoXsll.cc.

52 {
53
54 // check that there are no arguments
55
56 checkNArg(0,4,5);
57
58 checkNDaug(3);
59
60 // Check that the two leptons are the same type
61
62 EvtId lepton1type = getDaug(1);
63 EvtId lepton2type = getDaug(2);
64
65 int etyp = 0;
66 int mutyp = 0;
67 int tautyp = 0;
68 if ( lepton1type == EvtPDL::getId("e+") ||
69 lepton1type == EvtPDL::getId("e-") ) { etyp++;}
70 if ( lepton2type == EvtPDL::getId("e+") ||
71 lepton2type == EvtPDL::getId("e-") ) { etyp++;}
72 if ( lepton1type == EvtPDL::getId("mu+") ||
73 lepton1type == EvtPDL::getId("mu-") ) { mutyp++;}
74 if ( lepton2type == EvtPDL::getId("mu+") ||
75 lepton2type == EvtPDL::getId("mu-") ) { mutyp++;}
76 if ( lepton1type == EvtPDL::getId("tau+") ||
77 lepton1type == EvtPDL::getId("tau-") ) { tautyp++;}
78 if ( lepton2type == EvtPDL::getId("tau+") ||
79 lepton2type == EvtPDL::getId("tau-") ) { tautyp++;}
80
81 if ( etyp != 2 && mutyp != 2 && tautyp != 2 ) {
82
83 report(ERROR,"EvtGen") << "Expect two leptons of the same type in EvtBtoXsll.cc\n";
84 ::abort();
85 }
86
87 // Check that the second and third entries are leptons with positive
88 // and negative charge, respectively
89
90 int lpos = 0;
91 int lneg = 0;
92 if ( lepton1type == EvtPDL::getId("e+") ||
93 lepton1type == EvtPDL::getId("mu+") ||
94 lepton1type == EvtPDL::getId("tau+") ) { lpos++;}
95 if ( lepton2type == EvtPDL::getId("e-") ||
96 lepton2type == EvtPDL::getId("mu-") ||
97 lepton2type == EvtPDL::getId("tau-") ) { lneg++;}
98
99 if ( lpos != 1 || lneg != 1 ) {
100
101 report(ERROR,"EvtGen") << "Expect 2nd and 3rd particles to be positive and negative leptons in EvtBtoXsll.cc\n";
102 ::abort();
103 }
104
105
106 _mb=4.8;
107 _ms=0.2;
108 _mq=0.;
109 _pf=0.41;
110 _mxmin=1.1;
111 if ( getNArg()==4)
112 {
113 // b-quark mass
114 _mb = getArg(0);
115 // s-quark mass
116 _ms = getArg(1);
117 // spectator quark mass
118 _mq = getArg(2);
119 // Fermi motion parameter
120 _pf = getArg(3);
121 }
122 if ( getNArg()==5)
123 {
124 _mxmin = getArg(4);
125 }
126
127 _calcprob = new EvtBtoXsllUtil;
128
129 double ml = EvtPDL::getMeanMass(getDaug(1));
130
131 // determine the maximum probability density from dGdsProb
132
133 int i, j;
134 int nsteps = 100;
135 double s = 0.0;
136 double smin = 4.0 * ml * ml;
137 double smax = (_mb - _ms)*(_mb - _ms);
138 double probMax = -10000.0;
139 double sProbMax = -10.0;
140 double uProbMax = -10.0;
141
142 for (i=0;i<nsteps;i++)
143 {
144 s = smin + (i+0.0005)*(smax - smin)/(double)nsteps;
145 double prob = _calcprob->dGdsProb(_mb, _ms, ml, s);
146 if (prob > probMax)
147 {
148 sProbMax = s;
149 probMax = prob;
150 }
151 }
152
153 _dGdsProbMax = probMax;
154
155 if ( verbose() ) {
156 report(INFO,"EvtGen") << "dGdsProbMax = " << probMax << " for s = " << sProbMax << endl;
157 }
158
159 // determine the maximum probability density from dGdsdupProb
160
161 probMax = -10000.0;
162 sProbMax = -10.0;
163
164 for (i=0;i<nsteps;i++)
165 {
166 s = smin + (i+0.0005)*(smax - smin)/(double)nsteps;
167 double umax = sqrt((s - (_mb + _ms)*(_mb + _ms)) * (s - (_mb - _ms)*(_mb - _ms)));
168 for (j=0;j<nsteps;j++)
169 {
170 double u = -umax + (j+0.0005)*(2.0*umax)/(double)nsteps;
171 double prob = _calcprob->dGdsdupProb(_mb, _ms, ml, s, u);
172 if (prob > probMax)
173 {
174 sProbMax = s;
175 uProbMax = u;
176 probMax = prob;
177 }
178 }
179 }
180
181 _dGdsdupProbMax = probMax;
182
183 if ( verbose() ) {
184 report(INFO,"EvtGen") << "dGdsdupProbMax = " << probMax << " for s = " << sProbMax
185 << " and u = " << uProbMax << endl;
186 }
187
188}
@ ERROR
Definition: EvtReport.hh:49
double dGdsProb(double mb, double ms, double ml, double s)
double getArg(int j)
void checkNDaug(int d1, int d2=-1)
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtId.hh:27
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45

◆ initProbMax()

void EvtBtoXsll::initProbMax ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 190 of file EvtBtoXsll.cc.

190 {
191
192 noProbMax();
193
194}
void noProbMax()

The documentation for this class was generated from the following files: