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 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtCPUtil.cc
12 : //
13 : // Description: Utilities needed for generation of CP violating
14 : // decays.
15 : //
16 : // Modification history:
17 : //
18 : // RYD March 24, 1998 Module created
19 : //
20 : // COWAN June 10, 2009 Added methods for getting dGamma(s)
21 : // and dm(s) using B(s)0H and B(s)0L.
22 : //------------------------------------------------------------------------
23 : //
24 : #include "EvtGenBase/EvtPatches.hh"
25 : #include "EvtGenBase/EvtPatches.hh"
26 : #include "EvtGenBase/EvtParticle.hh"
27 : #include "EvtGenBase/EvtScalarParticle.hh"
28 : #include "EvtGenBase/EvtRandom.hh"
29 : #include "EvtGenBase/EvtCPUtil.hh"
30 : #include "EvtGenBase/EvtPDL.hh"
31 : #include "EvtGenBase/EvtReport.hh"
32 : #include "EvtGenBase/EvtSymTable.hh"
33 : #include "EvtGenBase/EvtConst.hh"
34 : #include <stdio.h>
35 : #include <stdlib.h>
36 :
37 : #include <assert.h>
38 : using std::endl;
39 :
40 0 : EvtCPUtil::EvtCPUtil(int mixingType) {
41 0 : _enableFlip = false;
42 0 : _mixingType = mixingType;
43 0 : }
44 :
45 0 : EvtCPUtil::~EvtCPUtil() {
46 0 : }
47 :
48 : EvtCPUtil* EvtCPUtil::getInstance() {
49 :
50 : static EvtCPUtil* theCPUtil = 0;
51 :
52 0 : if (theCPUtil == 0) {
53 0 : theCPUtil = new EvtCPUtil(1);
54 0 : }
55 :
56 0 : return theCPUtil;
57 :
58 : }
59 :
60 : //added two functions for finding the fraction of B0 tags for decays into
61 : //both CP eigenstates and non-CP eigenstates -- NK, Jan. 27th, 1998
62 :
63 : void EvtCPUtil::fractB0CP(EvtComplex Af, EvtComplex Abarf,
64 : double /*deltam*/, double beta, double &fract) {
65 :
66 : //This function returns the number of B0 tags for decays into CP-eigenstates
67 : //(the "probB0" in the new EvtOtherB)
68 :
69 : //double gamma_B = EvtPDL::getWidth(B0);
70 : //double xd = deltam/gamma_B;
71 : //double xd = 0.65;
72 : double ratio = 1/(1 + 0.65*0.65);
73 :
74 0 : EvtComplex rf, rbarf;
75 :
76 0 : rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
77 0 : rbarf = EvtComplex(1.0)/rf;
78 :
79 0 : double A2 = real(Af)*real(Af) + imag(Af)*imag(Af);
80 0 : double Abar2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
81 :
82 0 : double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
83 0 : double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
84 :
85 0 : fract = (Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio));
86 : return;
87 :
88 0 : }
89 :
90 : void EvtCPUtil::fractB0nonCP(EvtComplex Af, EvtComplex Abarf,
91 : EvtComplex Afbar, EvtComplex Abarfbar,
92 : double deltam, double beta,
93 : int flip, double &fract) {
94 :
95 : //this function returns the number of B0 tags for decays into non-CP eigenstates
96 : //(the "probB0" in the new EvtOtherB)
97 : //this needs more thought...
98 :
99 : //double gamma_B = EvtPDL::getWidth(B0);
100 : //report(Severity::Info,"EvtGen") << "gamma " << gamma_B<< endl;
101 : //double xd = deltam/gamma_B;
102 :
103 : //why is the width of B0 0 in PDL??
104 :
105 : double xd = 0.65;
106 0 : double gamma_B = deltam/xd;
107 : double IAf, IAfbar, IAbarf, IAbarfbar;
108 0 : EvtComplex rf, rfbar, rbarf, rbarfbar;
109 : double rf2, rfbar2, rbarf2, rbarfbar2;
110 : double Af2, Afbar2, Abarf2, Abarfbar2;
111 :
112 0 : rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
113 0 : rfbar = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarfbar/Afbar;
114 0 : rbarf = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Af/Abarf;
115 0 : rbarfbar = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Afbar/Abarfbar;
116 :
117 :
118 0 : rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
119 0 : rfbar2 = real(rfbar)*real(rfbar) + imag(rfbar)*imag(rfbar);
120 0 : rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
121 0 : rbarfbar2 = real(rbarfbar)*real(rbarfbar) + imag(rbarfbar)*imag(rbarfbar);
122 :
123 0 : Af2 = real(Af)*real(Af) + imag(Af)*imag(Af);
124 0 : Afbar2 = real(Afbar)*real(Afbar) + imag(Afbar)*imag(Afbar);
125 0 : Abarf2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
126 0 : Abarfbar2 = real(Abarfbar)*real(Abarfbar) + imag(Abarfbar)*imag(Abarfbar);
127 :
128 :
129 : //
130 : //IAf = integral(gamma(B0->f)), etc.
131 : //
132 :
133 0 : IAf = (Af2/(2*gamma_B))*(1+rf2+(1-rf2)/(1+xd*xd));
134 0 : IAfbar = (Afbar2/(2*gamma_B))*(1+rfbar2+(1-rfbar2)/(1+xd*xd));
135 0 : IAbarf = (Abarf2/(2*gamma_B))*(1+rbarf2+(1-rbarf2)/(1+xd*xd));
136 0 : IAbarfbar = (Abarfbar2/(2*gamma_B))*(1+rbarfbar2+(1-rbarfbar2)/(1+xd*xd));
137 :
138 : //flip specifies the relative fraction of fbar events
139 :
140 0 : fract = IAbarf/(IAbarf+IAf) + flip*IAbarfbar/(IAfbar+IAbarfbar);
141 :
142 :
143 : return;
144 0 : }
145 :
146 : void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
147 :
148 0 : if (_mixingType == EvtCPUtil::Coherent) {
149 :
150 0 : OtherCoherentB(p, t, otherb, probB0);
151 :
152 0 : } else if (_mixingType == EvtCPUtil::Incoherent) {
153 :
154 0 : OtherIncoherentB(p, t, otherb, probB0);
155 :
156 0 : }
157 :
158 0 : }
159 :
160 : void EvtCPUtil::OtherCoherentB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
161 :
162 : //Can not call this recursively!!!
163 : static int entryCount=0;
164 0 : entryCount++;
165 :
166 : //added by Lange Jan4,2000
167 0 : static EvtId B0B=EvtPDL::getId("anti-B0");
168 0 : static EvtId B0=EvtPDL::getId("B0");
169 0 : static EvtId BSB=EvtPDL::getId("anti-B_s0");
170 0 : static EvtId BS=EvtPDL::getId("B_s0");
171 :
172 0 : static EvtId UPS4S=EvtPDL::getId("Upsilon(4S)");
173 :
174 0 : int isB0=EvtRandom::Flat(0.0,1.0)<probB0;
175 :
176 : int idaug;
177 :
178 0 : p->setLifetime();
179 :
180 : // now get the time between the decay of this B and the other B!
181 :
182 0 : EvtParticle *parent=p->getParent();
183 :
184 : EvtParticle *other;
185 :
186 : bool incoherentmix=false;
187 :
188 0 : if ((parent!=0)&&(parent->getId()==B0||
189 0 : parent->getId()==B0B||
190 0 : parent->getId()==BS||
191 0 : parent->getId()==BSB)) {
192 : incoherentmix=true;
193 0 : }
194 :
195 0 : if (incoherentmix) parent=parent->getParent();
196 :
197 0 : if (parent==0||parent->getId()!=UPS4S) {
198 : //Need to make this more general, but for now
199 : //assume no parent. If we have parent of B we
200 : //need to charge conj. full decay tree.
201 :
202 :
203 0 : if (parent!=0) {
204 0 : report(Severity::Info,"EvtGen") << "p="<<EvtPDL::name(p->getId())
205 0 : << " parent="<<EvtPDL::name(parent->getId())
206 0 : << endl;
207 0 : }
208 0 : assert(parent==0);
209 0 : p->setLifetime();
210 0 : t=p->getLifetime();
211 : bool needToChargeConj=false;
212 0 : if (p->getId()==B0B&&isB0) needToChargeConj=true;
213 0 : if (p->getId()==B0&&!isB0) needToChargeConj=true;
214 0 : if (p->getId()==BSB&&isB0) needToChargeConj=true;
215 0 : if (p->getId()==BS&&!isB0) needToChargeConj=true;
216 :
217 0 : if (needToChargeConj) {
218 0 : p->setId( EvtPDL::chargeConj(p->getId()));
219 0 : if (incoherentmix) {
220 0 : p->getDaug(0)->setId(EvtPDL::chargeConj(p->getDaug(0)->getId()));
221 0 : }
222 : }
223 0 : otherb=EvtPDL::chargeConj(p->getId());
224 :
225 0 : entryCount--;
226 : return;
227 : }
228 : else{
229 0 : if (parent->getDaug(0)!=p){
230 0 : other=parent->getDaug(0);
231 : idaug=0;
232 0 : }
233 : else{
234 0 : other=parent->getDaug(1);
235 : idaug=1;
236 : }
237 : }
238 :
239 0 : if (parent != 0 ) {
240 :
241 : //if (entryCount>1){
242 : // report(Severity::Info,"EvtGen") << "Double CP decay:"<<entryCount<<endl;
243 : //}
244 :
245 : //kludge!! Lange Mar21, 2003
246 : // if the other B is an alias... don't change the flavor..
247 0 : if ( other->getId().isAlias() ) {
248 0 : OtherB(p,t,otherb);
249 0 : entryCount--;
250 0 : return;
251 :
252 : }
253 :
254 0 : if (entryCount==1){
255 :
256 0 : EvtVector4R p_init=other->getP4();
257 : //int decayed=other->getNDaug()>0;
258 0 : bool decayed = other->isDecayed();
259 :
260 0 : other->deleteTree();
261 :
262 : EvtScalarParticle* scalar_part;
263 :
264 0 : scalar_part=new EvtScalarParticle;
265 0 : if (isB0) {
266 0 : scalar_part->init(B0,p_init);
267 0 : }
268 : else{
269 0 : scalar_part->init(B0B,p_init);
270 : }
271 0 : other=(EvtParticle *)scalar_part;
272 : // other->set_type(EvtSpinType::SCALAR);
273 0 : other->setDiagonalSpinDensity();
274 :
275 0 : parent->insertDaugPtr(idaug,other);
276 :
277 0 : if (decayed){
278 : //report(Severity::Info,"EvtGen") << "In CP Util calling decay \n";
279 0 : other->decay();
280 0 : }
281 :
282 0 : }
283 :
284 0 : otherb=other->getId();
285 :
286 0 : other->setLifetime();
287 0 : t=p->getLifetime()-other->getLifetime();
288 :
289 0 : otherb = other->getId();
290 :
291 0 : }
292 : else {
293 0 : report(Severity::Info,"EvtGen") << "We have an error here!!!!"<<endl;
294 0 : otherb = EvtId(-1,-1);
295 : }
296 :
297 0 : entryCount--;
298 0 : return ;
299 0 : }
300 :
301 : // ========================================================================
302 : bool EvtCPUtil::isBsMixed ( EvtParticle * p )
303 : {
304 0 : if ( ! ( p->getParent() ) ) return false ;
305 :
306 0 : static EvtId BS0=EvtPDL::getId("B_s0");
307 0 : static EvtId BSB=EvtPDL::getId("anti-B_s0");
308 :
309 0 : if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
310 :
311 0 : if ( ( p->getParent()->getId() == BS0 ) ||
312 0 : ( p->getParent()->getId() == BSB ) ) return true ;
313 :
314 0 : return false ;
315 0 : }
316 :
317 : // ========================================================================
318 : bool EvtCPUtil::isB0Mixed ( EvtParticle * p )
319 : {
320 0 : if ( ! ( p->getParent() ) ) return false ;
321 :
322 0 : static EvtId B0 =EvtPDL::getId("B0");
323 0 : static EvtId B0B=EvtPDL::getId("anti-B0");
324 :
325 0 : if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
326 :
327 0 : if ( ( p->getParent()->getId() == B0 ) ||
328 0 : ( p->getParent()->getId() == B0B ) ) return true ;
329 :
330 0 : return false ;
331 0 : }
332 : //============================================================================
333 : // Return the tag of the event (ie the anti-flavour of the produced
334 : // B meson). Flip the flavour of the event with probB probability
335 : //============================================================================
336 : void EvtCPUtil::OtherIncoherentB( EvtParticle * p ,
337 : double & t ,
338 : EvtId & otherb ,
339 : double probB )
340 : {
341 : //std::cout<<"New routine running"<<endl;
342 : //if(p->getId() == B0 || p->getId() == B0B)
343 : //added by liming Zhang
344 0 : enableFlip();
345 0 : if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
346 0 : p->getParent()->setLifetime() ;
347 0 : t = p->getParent()->getLifetime() ;
348 0 : }
349 : else {
350 0 : p->setLifetime() ;
351 0 : t = p->getLifetime() ;
352 : }
353 :
354 0 : if ( flipIsEnabled() ) {
355 : //std::cout << " liming << flipIsEnabled " << std::endl;
356 : // Flip the flavour of the particle with probability probB
357 0 : bool isFlipped = ( EvtRandom::Flat( 0. , 1. ) < probB ) ;
358 :
359 0 : if ( isFlipped ) {
360 0 : if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
361 0 : p->getParent()
362 0 : ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
363 0 : p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
364 0 : }
365 : else {
366 0 : p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
367 : }
368 : }
369 0 : }
370 :
371 0 : if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
372 : // if B has mixed, tag flavour is charge conjugate of parent of B-meson
373 0 : otherb = EvtPDL::chargeConj( p->getParent()->getId() ) ;
374 0 : }
375 : else {
376 : // else it is opposite flavour than this B hadron
377 0 : otherb = EvtPDL::chargeConj( p->getId() ) ;
378 : }
379 :
380 0 : return ;
381 : }
382 : //============================================================================
383 : void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb){
384 :
385 0 : static EvtId BSB=EvtPDL::getId("anti-B_s0");
386 0 : static EvtId BS0=EvtPDL::getId("B_s0");
387 0 : static EvtId B0B=EvtPDL::getId("anti-B0");
388 0 : static EvtId B0=EvtPDL::getId("B0");
389 0 : static EvtId D0B=EvtPDL::getId("anti-D0");
390 0 : static EvtId D0=EvtPDL::getId("D0");
391 0 : static EvtId UPS4=EvtPDL::getId("Upsilon(4S)");
392 :
393 0 : if (p->getId()==BS0||p->getId()==BSB){
394 0 : static double ctauL=EvtPDL::getctau(EvtPDL::getId("B_s0L"));
395 0 : static double ctauH=EvtPDL::getctau(EvtPDL::getId("B_s0H"));
396 0 : static double ctau=ctauL<ctauH?ctauH:ctauL;
397 0 : t=-log(EvtRandom::Flat())*ctau;
398 0 : EvtParticle* parent=p->getParent();
399 0 : if (parent!=0&&(parent->getId()==BS0||parent->getId()==BSB)){
400 0 : if (parent->getId()==BS0) otherb=BSB;
401 0 : if (parent->getId()==BSB) otherb=BS0;
402 0 : parent->setLifetime(t);
403 0 : return;
404 : }
405 0 : if (p->getId()==BS0) otherb=BSB;
406 0 : if (p->getId()==BSB) otherb=BS0;
407 0 : p->setLifetime(t);
408 0 : return;
409 : }
410 :
411 0 : if (p->getId()==D0||p->getId()==D0B){
412 0 : static double ctauL=EvtPDL::getctau(EvtPDL::getId("D0L"));
413 0 : static double ctauH=EvtPDL::getctau(EvtPDL::getId("D0H"));
414 0 : static double ctau=ctauL<ctauH?ctauH:ctauL;
415 0 : t=-log(EvtRandom::Flat())*ctau;
416 0 : EvtParticle* parent=p->getParent();
417 0 : if (parent!=0&&(parent->getId()==D0||parent->getId()==D0B)){
418 0 : if (parent->getId()==D0) otherb=D0B;
419 0 : if (parent->getId()==D0B) otherb=D0;
420 0 : parent->setLifetime(t);
421 0 : return;
422 : }
423 0 : if (p->getId()==D0) otherb=D0B;
424 0 : if (p->getId()==D0B) otherb=D0;
425 0 : p->setLifetime(t);
426 0 : return;
427 : }
428 :
429 0 : p->setLifetime();
430 :
431 : // now get the time between the decay of this B and the other B!
432 :
433 0 : EvtParticle *parent=p->getParent();
434 :
435 0 : if (parent==0||parent->getId()!=UPS4) {
436 : //report(Severity::Error,"EvtGen") <<
437 : // "Warning CP violation with B having no parent!"<<endl;
438 0 : t=p->getLifetime();
439 0 : if (p->getId()==B0) otherb=B0B;
440 0 : if (p->getId()==B0B) otherb=B0;
441 0 : if (p->getId()==BS0) otherb=BSB;
442 0 : if (p->getId()==BSB) otherb=BS0;
443 0 : return;
444 : }
445 : else{
446 0 : if (parent->getDaug(0)!=p){
447 0 : otherb=parent->getDaug(0)->getId();
448 0 : parent->getDaug(0)->setLifetime();
449 0 : t=p->getLifetime()-parent->getDaug(0)->getLifetime();
450 0 : }
451 : else{
452 0 : otherb=parent->getDaug(1)->getId();
453 0 : parent->getDaug(1)->setLifetime();
454 0 : t=p->getLifetime()-parent->getDaug(1)->getLifetime();
455 : }
456 : }
457 :
458 :
459 0 : return ;
460 0 : }
461 :
462 : // No CP violation is assumed
463 : void EvtCPUtil::incoherentMix(const EvtId id, double &t, int &mix){
464 :
465 0 : int stdHepNum=EvtPDL::getStdHep(id);
466 0 : stdHepNum=abs(stdHepNum);
467 :
468 0 : EvtId partId=EvtPDL::evtIdFromStdHep(stdHepNum);
469 :
470 0 : std::string partName=EvtPDL::name(partId);
471 0 : std::string hname=partName+std::string("H");
472 0 : std::string lname=partName+std::string("L");
473 :
474 0 : EvtId lId=EvtPDL::getId(lname);
475 0 : EvtId hId=EvtPDL::getId(hname);
476 :
477 0 : double ctauL=EvtPDL::getctau(lId);
478 0 : double ctauH=EvtPDL::getctau(hId);
479 :
480 : // Bug Fixed: Corrected the average as gamma is the relevent parameter
481 0 : double ctau=2.0*(ctauL*ctauH)/(ctauL+ctauH);
482 : //double ctau=0.5*(ctauL+ctauH);
483 :
484 : // Bug Fixed: ctau definition changed above
485 : //double y=(ctauH-ctauL)/(2*ctau);
486 0 : double y=(ctauH-ctauL)/(ctauH+ctauL);
487 :
488 : //deltam and qoverp defined in DECAY.DEC
489 :
490 0 : std::string qoverpParmName=std::string("qoverp_incohMix_")+partName;
491 0 : std::string mdParmName=std::string("dm_incohMix_")+partName;
492 0 : int ierr;
493 0 : double qoverp=atof(EvtSymTable::get(qoverpParmName,ierr).c_str());
494 0 : double x=atof(EvtSymTable::get(mdParmName,ierr).c_str())*ctau/EvtConst::c;
495 : double fac;
496 :
497 0 : if(id==partId){
498 0 : fac=1.0/(qoverp*qoverp);
499 0 : }
500 : else{
501 : fac=qoverp*qoverp;
502 : }
503 :
504 0 : double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2.0+x*x-y*y));
505 :
506 : int mixsign;
507 :
508 0 : mixsign=(mixprob>EvtRandom::Flat(0.0,1.0))?-1:1;
509 :
510 : double prob;
511 :
512 : // Find the longest of the two lifetimes
513 0 : double ctaulong = ctauL<=ctauH?ctauH:ctauL;
514 :
515 : // Bug fixed: Ensure cosine argument is dimensionless so /ctau
516 0 : do{
517 0 : t=-log(EvtRandom::Flat())*ctaulong;
518 0 : prob=1.0+exp(-2.0*fabs(y)*t/ctau)+mixsign*2.0*exp(-fabs(y)*t/ctau)*cos(x*t/ctau);
519 0 : }while(prob<4.0*EvtRandom::Flat());
520 :
521 0 : mix=0;
522 :
523 0 : if (mixsign==-1) mix=1;
524 :
525 : return;
526 0 : }
527 :
528 :
529 : double EvtCPUtil::getDeltaGamma(const EvtId id){
530 :
531 0 : int stdHepNum = EvtPDL::getStdHep(id);
532 0 : stdHepNum = abs(stdHepNum);
533 0 : EvtId partId = EvtPDL::evtIdFromStdHep(stdHepNum);
534 :
535 0 : std::string partName = EvtPDL::name(partId);
536 0 : std::string hname = partName + std::string("H");
537 0 : std::string lname = partName + std::string("L");
538 :
539 0 : EvtId lId = EvtPDL::getId(lname);
540 0 : EvtId hId = EvtPDL::getId(hname);
541 :
542 0 : double ctauL = EvtPDL::getctau(lId);
543 0 : double ctauH = EvtPDL::getctau(hId);
544 :
545 0 : double dGamma = (1/ctauL - 1/ctauH)*EvtConst::c;
546 : return dGamma;
547 0 : }
548 :
549 : double EvtCPUtil::getDeltaM(const EvtId id){
550 :
551 0 : int stdHepNum = EvtPDL::getStdHep(id);
552 0 : stdHepNum = abs(stdHepNum);
553 0 : EvtId partId = EvtPDL::evtIdFromStdHep(stdHepNum);
554 :
555 0 : std::string partName = EvtPDL::name(partId);
556 0 : std::string parmName = std::string("dm_incohMix_") + partName;
557 :
558 0 : int ierr;
559 0 : double dM = atof(EvtSymTable::get(parmName,ierr).c_str());
560 : return dM;
561 0 : }
562 :
563 0 : bool EvtCPUtil::flipIsEnabled() { return _enableFlip ; }
564 0 : void EvtCPUtil::enableFlip() { _enableFlip = true ; }
565 0 : void EvtCPUtil::disableFlip() { _enableFlip = false ; }
566 :
|
