#include <EvtCPUtil.hh>
Static Public Member Functions | |
| void | fractB0CP (EvtComplex Af, EvtComplex Abarf, double deltam, double beta, double &fract) |
| void | fractB0nonCP (EvtComplex Af, EvtComplex Abarf, EvtComplex Afbar, EvtComplex Abarfbar, double deltam, double beta, int flip, double &fract) |
| void | incoherentMix (const EvtId id, double &t, int &mix) |
| void | OtherB (EvtParticle *p, double &t, EvtId &otherb, double probB0) |
| void | OtherB (EvtParticle *p, double &t, EvtId &otherb) |
|
||||||||||||||||||||||||
|
00041 {
00042 //this function returns the number of B0 tags for decays into CP-eigenstates
00043 //(the "probB0" in the new EvtOtherB)
00044
00045 //double gamma_B = EvtPDL::getWidth(B0);
00046 //double xd = deltam/gamma_B;
00047 //double xd = 0.65;
00048 double ratio = 1/(1 + 0.65*0.65);
00049
00050 EvtComplex rf, rbarf;
00051
00052 rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
00053 rbarf = EvtComplex(1.0)/rf;
00054
00055 double A2 = real(Af)*real(Af) + imag(Af)*imag(Af);
00056 double Abar2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
00057
00058 double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
00059 double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
00060
00061 fract = (Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio));
00062 return;
00063
00064 }
|
|
||||||||||||||||||||||||||||||||||||
|
00068 {
00069
00070 //this function returns the number of B0 tags for decays into non-CP eigenstates
00071 //(the "probB0" in the new EvtOtherB)
00072 //this needs more thought...
00073
00074 //double gamma_B = EvtPDL::getWidth(B0);
00075 //report(INFO,"EvtGen") << "gamma " << gamma_B<< endl;
00076 //double xd = deltam/gamma_B;
00077
00078 //why is the width of B0 0 in PDL??
00079
00080 double xd = 0.65;
00081 double gamma_B = deltam/xd;
00082 double IAf, IAfbar, IAbarf, IAbarfbar;
00083 EvtComplex rf, rfbar, rbarf, rbarfbar;
00084 double rf2, rfbar2, rbarf2, rbarfbar2;
00085 double Af2, Afbar2, Abarf2, Abarfbar2;
00086
00087 rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
00088 rfbar = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarfbar/Afbar;
00089 rbarf = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Af/Abarf;
00090 rbarfbar = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Afbar/Abarfbar;
00091
00092
00093 rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
00094 rfbar2 = real(rfbar)*real(rfbar) + imag(rfbar)*imag(rfbar);
00095 rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
00096 rbarfbar2 = real(rbarfbar)*real(rbarfbar) + imag(rbarfbar)*imag(rbarfbar);
00097
00098 Af2 = real(Af)*real(Af) + imag(Af)*imag(Af);
00099 Afbar2 = real(Afbar)*real(Afbar) + imag(Afbar)*imag(Afbar);
00100 Abarf2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
00101 Abarfbar2 = real(Abarfbar)*real(Abarfbar) + imag(Abarfbar)*imag(Abarfbar);
00102
00103
00104 //
00105 //IAf = integral(gamma(B0->f)), etc.
00106 //
00107
00108 IAf = (Af2/(2*gamma_B))*(1+rf2+(1-rf2)/(1+xd*xd));
00109 IAfbar = (Afbar2/(2*gamma_B))*(1+rfbar2+(1-rfbar2)/(1+xd*xd));
00110 IAbarf = (Abarf2/(2*gamma_B))*(1+rbarf2+(1-rbarf2)/(1+xd*xd));
00111 IAbarfbar = (Abarfbar2/(2*gamma_B))*(1+rbarfbar2+(1-rbarfbar2)/(1+xd*xd));
00112
00113 //flip specifies the relative fraction of fbar events
00114
00115 fract = IAbarf/(IAbarf+IAf) + flip*IAbarfbar/(IAfbar+IAbarfbar);
00116
00117
00118 return;
00119 }
|
|
||||||||||||||||
|
00293 {
00294
00295 int stdHepNum=EvtPDL::getStdHep(id);
00296 stdHepNum=abs(stdHepNum);
00297
00298 EvtId partId=EvtPDL::evtIdFromStdHep(stdHepNum);
00299
00300 std::string partName=EvtPDL::name(partId);
00301 std::string hname=partName+std::string("H");
00302 std::string lname=partName+std::string("L");
00303
00304 EvtId lId=EvtPDL::getId(lname);
00305 EvtId hId=EvtPDL::getId(hname);
00306
00307 double ctauL=EvtPDL::getctau(lId);
00308 double ctauH=EvtPDL::getctau(hId);
00309 double ctau=0.5*(ctauL+ctauH);
00310 double y=(ctauH-ctauL)/ctau;
00311
00312 //need to figure out how to get these parameters into the code...
00313
00314 std::string qoverpParmName=std::string("qoverp_incohMix_")+partName;
00315 std::string mdParmName=std::string("dm_incohMix_")+partName;
00316 int ierr;
00317 double qoverp=atof(EvtSymTable::Get(qoverpParmName,ierr).c_str());
00318 double x=atof(EvtSymTable::Get(mdParmName,ierr).c_str())*ctau/EvtConst::c;
00319 double fac;
00320
00321 if(id==partId){
00322 fac=1.0/(qoverp*qoverp);
00323 }
00324 else{
00325 fac=qoverp*qoverp;
00326 }
00327
00328 double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2+x*x-y*y));
00329
00330 int mixsign;
00331
00332 mixsign=(mixprob>EvtRandom::Flat(0.0,1.0))?-1:1;
00333
00334 double prob;
00335
00336 do{
00337 t=-log(EvtRandom::Flat())*ctauL;
00338 prob=1+exp(y*t/ctau)+mixsign*2*exp(0.5*y*t/ctau)*cos(x*t);
00339 }while(prob<4*EvtRandom::Flat());
00340
00341 mix=0;
00342
00343 if (mixsign==-1) mix=1;
00344
00345 return;
00346 }
|
|
||||||||||||||||||||
|
00121 {
00122
00123 //Can not call this recursively!!!
00124 static int entryCount=0;
00125 entryCount++;
00126
00127 //added by Lange Jan4,2000
00128 static EvtId B0B=EvtPDL::getId("anti-B0");
00129 static EvtId B0=EvtPDL::getId("B0");
00130
00131 int isB0=EvtRandom::Flat(0.0,1.0)<probB0;
00132
00133 int idaug;
00134
00135 p->setLifetime();
00136
00137 // now get the time between the decay of this B and the other B!
00138
00139 EvtParticle *parent=p->getParent();
00140
00141 EvtParticle *other;
00142
00143 if (parent==0) {
00144 //report(ERROR,"EvtGen") <<
00145 // "Warning CP violation with B having no parent!"<<endl;
00146 p->setLifetime();
00147 t=p->getLifetime();
00148 if (p->getId()==B0) otherb=B0B;
00149 if (p->getId()==B0B) otherb=B0;
00150 entryCount--;
00151 return;
00152 }
00153 else{
00154 if (parent->getDaug(0)!=p){
00155 other=parent->getDaug(0);
00156 idaug=0;
00157 }
00158 else{
00159 other=parent->getDaug(1);
00160 idaug=1;
00161 }
00162 }
00163
00164 if (parent != 0 ) {
00165
00166 //if (entryCount>1){
00167 // report(INFO,"EvtGen") << "Double CP decay:"<<entryCount<<endl;
00168 //}
00169
00170 //kludge!! Lange Mar21, 2003
00171 // if the other B is an alias... don't change the flavor..
00172 if ( other->getId().isAlias() ) {
00173 OtherB(p,t,otherb);
00174 entryCount--;
00175 return;
00176
00177 }
00178
00179 if (entryCount==1){
00180
00181 EvtVector4R p_init=other->getP4();
00182 //int decayed=other->getNDaug()>0;
00183 bool decayed = other->isDecayed();
00184
00185 other->deleteTree();
00186
00187 EvtScalarParticle* scalar_part;
00188
00189 scalar_part=new EvtScalarParticle;
00190 if (isB0) {
00191 scalar_part->init(B0,p_init);
00192 }
00193 else{
00194 scalar_part->init(B0B,p_init);
00195 }
00196 other=(EvtParticle *)scalar_part;
00197 // other->set_type(EvtSpinType::SCALAR);
00198 other->setDiagonalSpinDensity();
00199
00200 parent->insertDaugPtr(idaug,other);
00201
00202 if (decayed){
00203 //report(INFO,"EvtGen") << "In CP Util calling decay \n";
00204 other->decay();
00205 }
00206
00207 }
00208
00209 otherb=other->getId();
00210
00211 other->setLifetime();
00212 t=p->getLifetime()-other->getLifetime();
00213
00214 otherb = other->getId();
00215
00216 }
00217 else {
00218 report(INFO,"EvtGen") << "We have an error here!!!!"<<endl;
00219 otherb = EvtId(-1,-1);
00220 }
00221
00222 entryCount--;
00223 return ;
00224 }
|
|
||||||||||||||||
|
00228 {
00229
00230
00231 static EvtId BSB=EvtPDL::getId("anti-B_s0");
00232 static EvtId BS0=EvtPDL::getId("B_s0");
00233 static EvtId B0B=EvtPDL::getId("anti-B0");
00234 static EvtId B0=EvtPDL::getId("B0");
00235 static EvtId UPS4=EvtPDL::getId("Upsilon(4S)");
00236
00237 if (p->getId()==BS0||p->getId()==BSB){
00238 static double ctauL=EvtPDL::getctau(EvtPDL::getId("B_s0L"));
00239 static double ctauH=EvtPDL::getctau(EvtPDL::getId("B_s0H"));
00240 static double ctau=ctauL<ctauH?ctauH:ctauL;
00241 t=-log(EvtRandom::Flat())*ctau;
00242 EvtParticle* parent=p->getParent();
00243 if (parent!=0&&(parent->getId()==BS0||parent->getId()==BSB)){
00244 if (parent->getId()==BS0) otherb=BSB;
00245 if (parent->getId()==BSB) otherb=BS0;
00246 parent->setLifetime(t);
00247 return;
00248 }
00249 if (p->getId()==BS0) otherb=BSB;
00250 if (p->getId()==BSB) otherb=BS0;
00251 p->setLifetime(t);
00252 return;
00253 }
00254
00255
00256
00257
00258 p->setLifetime();
00259
00260
00261 // now get the time between the decay of this B and the other B!
00262
00263 EvtParticle *parent=p->getParent();
00264
00265 if (parent==0||parent->getId()!=UPS4) {
00266 //report(ERROR,"EvtGen") <<
00267 // "Warning CP violation with B having no parent!"<<endl;
00268 t=p->getLifetime();
00269 if (p->getId()==B0) otherb=B0B;
00270 if (p->getId()==B0B) otherb=B0;
00271 if (p->getId()==BS0) otherb=BSB;
00272 if (p->getId()==BSB) otherb=BS0;
00273 return;
00274 }
00275 else{
00276 if (parent->getDaug(0)!=p){
00277 otherb=parent->getDaug(0)->getId();
00278 parent->getDaug(0)->setLifetime();
00279 t=p->getLifetime()-parent->getDaug(0)->getLifetime();
00280 }
00281 else{
00282 otherb=parent->getDaug(1)->getId();
00283 parent->getDaug(1)->setLifetime();
00284 t=p->getLifetime()-parent->getDaug(1)->getLifetime();
00285 }
00286 }
00287
00288
00289 return ;
00290 }
|
1.3.9.1