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 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtCPUtil.cc
13 // Description: Utilities needed for generation of CP violating
16 // Modification history:
18 // RYD March 24, 1998 Module created
20 //------------------------------------------------------------------------
22 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtPatches.hh"
24 #include "EvtGenBase/EvtParticle.hh"
25 #include "EvtGenBase/EvtScalarParticle.hh"
26 #include "EvtGenBase/EvtRandom.hh"
27 #include "EvtGenBase/EvtCPUtil.hh"
28 #include "EvtGenBase/EvtPDL.hh"
29 #include "EvtGenBase/EvtReport.hh"
30 #include "EvtGenBase/EvtSymTable.hh"
31 #include "EvtGenBase/EvtConst.hh"
39 //added two functions for finding the fraction of B0 tags for decays into
40 //both CP eigenstates and non-CP eigenstates -- NK, Jan. 27th, 1998
42 void EvtCPUtil::fractB0CP(EvtComplex Af, EvtComplex Abarf,
43 double deltam, double beta, double &fract) {
44 //this function returns the number of B0 tags for decays into CP-eigenstates
45 //(the "probB0" in the new EvtOtherB)
47 //double gamma_B = EvtPDL::getWidth(B0);
48 //double xd = deltam/gamma_B;
50 double ratio = 1/(1 + 0.65*0.65);
54 rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
55 rbarf = EvtComplex(1.0)/rf;
57 double A2 = real(Af)*real(Af) + imag(Af)*imag(Af);
58 double Abar2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
60 double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
61 double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
63 fract = (Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio));
67 void EvtCPUtil::fractB0nonCP(EvtComplex Af, EvtComplex Abarf,
68 EvtComplex Afbar, EvtComplex Abarfbar,
69 double deltam, double beta,
70 int flip, double &fract) {
72 //this function returns the number of B0 tags for decays into non-CP eigenstates
73 //(the "probB0" in the new EvtOtherB)
74 //this needs more thought...
76 //double gamma_B = EvtPDL::getWidth(B0);
77 //report(INFO,"EvtGen") << "gamma " << gamma_B<< endl;
78 //double xd = deltam/gamma_B;
80 //why is the width of B0 0 in PDL??
83 double gamma_B = deltam/xd;
84 double IAf, IAfbar, IAbarf, IAbarfbar;
85 EvtComplex rf, rfbar, rbarf, rbarfbar;
86 double rf2, rfbar2, rbarf2, rbarfbar2;
87 double Af2, Afbar2, Abarf2, Abarfbar2;
89 rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
90 rfbar = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarfbar/Afbar;
91 rbarf = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Af/Abarf;
92 rbarfbar = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Afbar/Abarfbar;
95 rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
96 rfbar2 = real(rfbar)*real(rfbar) + imag(rfbar)*imag(rfbar);
97 rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
98 rbarfbar2 = real(rbarfbar)*real(rbarfbar) + imag(rbarfbar)*imag(rbarfbar);
100 Af2 = real(Af)*real(Af) + imag(Af)*imag(Af);
101 Afbar2 = real(Afbar)*real(Afbar) + imag(Afbar)*imag(Afbar);
102 Abarf2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
103 Abarfbar2 = real(Abarfbar)*real(Abarfbar) + imag(Abarfbar)*imag(Abarfbar);
107 //IAf = integral(gamma(B0->f)), etc.
110 IAf = (Af2/(2*gamma_B))*(1+rf2+(1-rf2)/(1+xd*xd));
111 IAfbar = (Afbar2/(2*gamma_B))*(1+rfbar2+(1-rfbar2)/(1+xd*xd));
112 IAbarf = (Abarf2/(2*gamma_B))*(1+rbarf2+(1-rbarf2)/(1+xd*xd));
113 IAbarfbar = (Abarfbar2/(2*gamma_B))*(1+rbarfbar2+(1-rbarfbar2)/(1+xd*xd));
115 //flip specifies the relative fraction of fbar events
117 fract = IAbarf/(IAbarf+IAf) + flip*IAbarfbar/(IAfbar+IAbarfbar);
123 void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
125 //Can not call this recursively!!!
126 static int entryCount=0;
129 //added by Lange Jan4,2000
130 static EvtId B0B=EvtPDL::getId("anti-B0");
131 static EvtId B0=EvtPDL::getId("B0");
132 static EvtId BSB=EvtPDL::getId("anti-B_s0");
133 static EvtId BS=EvtPDL::getId("B_s0");
135 static EvtId UPS4S=EvtPDL::getId("Upsilon(4S)");
137 int isB0=EvtRandom::Flat(0.0,1.0)<probB0;
143 // now get the time between the decay of this B and the other B!
145 EvtParticle *parent=p->getParent();
149 bool incoherentmix=false;
151 if ((parent!=0)&&(parent->getId()==B0||
152 parent->getId()==B0B||
153 parent->getId()==BS||
154 parent->getId()==BSB)) {
158 if (incoherentmix) parent=parent->getParent();
160 if (parent==0||parent->getId()!=UPS4S) {
161 //Need to make this more general, but for now
162 //assume no parent. If we have parent of B we
163 //need to charge conj. full decay tree.
167 report(INFO,"EvtGen") << "p="<<EvtPDL::name(p->getId())
168 << " parent="<<EvtPDL::name(parent->getId())
174 bool needToChargeConj=false;
175 if (p->getId()==B0B&&isB0) needToChargeConj=true;
176 if (p->getId()==B0&&!isB0) needToChargeConj=true;
177 if (p->getId()==BSB&&isB0) needToChargeConj=true;
178 if (p->getId()==BS&&!isB0) needToChargeConj=true;
180 if (needToChargeConj) {
181 p->setId( EvtPDL::chargeConj(p->getId()));
183 p->getDaug(0)->setId(EvtPDL::chargeConj(p->getDaug(0)->getId()));
186 otherb=EvtPDL::chargeConj(p->getId());
192 if (parent->getDaug(0)!=p){
193 other=parent->getDaug(0);
197 other=parent->getDaug(1);
205 // report(INFO,"EvtGen") << "Double CP decay:"<<entryCount<<endl;
208 //kludge!! Lange Mar21, 2003
209 // if the other B is an alias... don't change the flavor..
210 if ( other->getId().isAlias() ) {
219 EvtVector4R p_init=other->getP4();
220 //int decayed=other->getNDaug()>0;
221 bool decayed = other->isDecayed();
225 EvtScalarParticle* scalar_part;
227 scalar_part=new EvtScalarParticle;
229 scalar_part->init(B0,p_init);
232 scalar_part->init(B0B,p_init);
234 other=(EvtParticle *)scalar_part;
235 // other->set_type(EvtSpinType::SCALAR);
236 other->setDiagonalSpinDensity();
238 parent->insertDaugPtr(idaug,other);
241 //report(INFO,"EvtGen") << "In CP Util calling decay \n";
247 otherb=other->getId();
249 other->setLifetime();
250 t=p->getLifetime()-other->getLifetime();
252 otherb = other->getId();
256 report(INFO,"EvtGen") << "We have an error here!!!!"<<endl;
257 otherb = EvtId(-1,-1);
266 void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb){
269 static EvtId BSB=EvtPDL::getId("anti-B_s0");
270 static EvtId BS0=EvtPDL::getId("B_s0");
271 static EvtId B0B=EvtPDL::getId("anti-B0");
272 static EvtId B0=EvtPDL::getId("B0");
273 static EvtId D0B=EvtPDL::getId("anti-D0");
274 static EvtId D0=EvtPDL::getId("D0");
275 static EvtId UPS4=EvtPDL::getId("Upsilon(4S)");
277 if (p->getId()==BS0||p->getId()==BSB){
278 static double ctauL=EvtPDL::getctau(EvtPDL::getId("B_s0L"));
279 static double ctauH=EvtPDL::getctau(EvtPDL::getId("B_s0H"));
280 static double ctau=ctauL<ctauH?ctauH:ctauL;
281 t=-log(EvtRandom::Flat())*ctau;
282 EvtParticle* parent=p->getParent();
283 if (parent!=0&&(parent->getId()==BS0||parent->getId()==BSB)){
284 if (parent->getId()==BS0) otherb=BSB;
285 if (parent->getId()==BSB) otherb=BS0;
286 parent->setLifetime(t);
289 if (p->getId()==BS0) otherb=BSB;
290 if (p->getId()==BSB) otherb=BS0;
295 if (p->getId()==D0||p->getId()==D0B){
296 static double ctauL=EvtPDL::getctau(EvtPDL::getId("D0L"));
297 static double ctauH=EvtPDL::getctau(EvtPDL::getId("D0H"));
298 static double ctau=ctauL<ctauH?ctauH:ctauL;
299 t=-log(EvtRandom::Flat())*ctau;
300 EvtParticle* parent=p->getParent();
301 if (parent!=0&&(parent->getId()==D0||parent->getId()==D0B)){
302 if (parent->getId()==D0) otherb=D0B;
303 if (parent->getId()==D0B) otherb=D0;
304 parent->setLifetime(t);
307 if (p->getId()==D0) otherb=D0B;
308 if (p->getId()==D0B) otherb=D0;
319 // now get the time between the decay of this B and the other B!
321 EvtParticle *parent=p->getParent();
323 if (parent==0||parent->getId()!=UPS4) {
324 //report(ERROR,"EvtGen") <<
325 // "Warning CP violation with B having no parent!"<<endl;
327 if (p->getId()==B0) otherb=B0B;
328 if (p->getId()==B0B) otherb=B0;
329 if (p->getId()==BS0) otherb=BSB;
330 if (p->getId()==BSB) otherb=BS0;
334 if (parent->getDaug(0)!=p){
335 otherb=parent->getDaug(0)->getId();
336 parent->getDaug(0)->setLifetime();
337 t=p->getLifetime()-parent->getDaug(0)->getLifetime();
340 otherb=parent->getDaug(1)->getId();
341 parent->getDaug(1)->setLifetime();
342 t=p->getLifetime()-parent->getDaug(1)->getLifetime();
351 void EvtCPUtil::incoherentMix(const EvtId id, double &t, int &mix){
353 int stdHepNum=EvtPDL::getStdHep(id);
354 stdHepNum=abs(stdHepNum);
356 EvtId partId=EvtPDL::evtIdFromStdHep(stdHepNum);
358 std::string partName=EvtPDL::name(partId);
359 std::string hname=partName+std::string("H");
360 std::string lname=partName+std::string("L");
362 EvtId lId=EvtPDL::getId(lname);
363 EvtId hId=EvtPDL::getId(hname);
365 double ctauL=EvtPDL::getctau(lId);
366 double ctauH=EvtPDL::getctau(hId);
367 double ctau=0.5*(ctauL+ctauH);
368 double y=(ctauH-ctauL)/ctau;
370 //need to figure out how to get these parameters into the code...
372 std::string qoverpParmName=std::string("qoverp_incohMix_")+partName;
373 std::string mdParmName=std::string("dm_incohMix_")+partName;
375 double qoverp=atof(EvtSymTable::get(qoverpParmName,ierr).c_str());
376 double x=atof(EvtSymTable::get(mdParmName,ierr).c_str())*ctau/EvtConst::c;
380 fac=1.0/(qoverp*qoverp);
386 double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2+x*x-y*y));
390 mixsign=(mixprob>EvtRandom::Flat(0.0,1.0))?-1:1;
395 t=-log(EvtRandom::Flat())*ctauL;
396 prob=1.0+exp(2.0*y*t/ctau)+mixsign*2.0*exp(y*t/ctau)*cos(x*t);
397 }while(prob<4*EvtRandom::Flat());
401 if (mixsign==-1) mix=1;