00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include "EvtGenBase/EvtPatches.hh"
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtScalarParticle.hh"
00026 #include "EvtGenBase/EvtRandom.hh"
00027 #include "EvtGenBase/EvtCPUtil.hh"
00028 #include "EvtGenBase/EvtPDL.hh"
00029 #include "EvtGenBase/EvtReport.hh"
00030 #include "EvtGenBase/EvtSymTable.hh"
00031 #include "EvtGenBase/EvtConst.hh"
00032
00033 #include <assert.h>
00034 #include <stdlib.h>
00035 using std::endl;
00036
00037
00038
00039
00040
00041 void EvtCPUtil::fractB0CP(EvtComplex Af, EvtComplex Abarf,
00042 double deltam, double beta, double &fract) {
00043
00044
00045
00046
00047
00048
00049 double ratio = 1/(1 + 0.65*0.65);
00050
00051 EvtComplex rf, rbarf;
00052
00053 rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
00054 rbarf = EvtComplex(1.0)/rf;
00055
00056 double A2 = real(Af)*real(Af) + imag(Af)*imag(Af);
00057 double Abar2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
00058
00059 double rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
00060 double rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
00061
00062 fract = (Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio))/(Abar2*(1+ rbarf2 + (1 - rbarf2)*ratio) + A2*(1+ rf2 + (1 - rf2)*ratio));
00063 return;
00064
00065 }
00066 void EvtCPUtil::fractB0nonCP(EvtComplex Af, EvtComplex Abarf,
00067 EvtComplex Afbar, EvtComplex Abarfbar,
00068 double deltam, double beta,
00069 int flip, double &fract) {
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 double xd = 0.65;
00082 double gamma_B = deltam/xd;
00083 double IAf, IAfbar, IAbarf, IAbarfbar;
00084 EvtComplex rf, rfbar, rbarf, rbarfbar;
00085 double rf2, rfbar2, rbarf2, rbarfbar2;
00086 double Af2, Afbar2, Abarf2, Abarfbar2;
00087
00088 rf = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarf/Af;
00089 rfbar = EvtComplex(cos(2.0*beta),sin(2.0*beta))*Abarfbar/Afbar;
00090 rbarf = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Af/Abarf;
00091 rbarfbar = EvtComplex(cos(-2.0*beta),sin(-2.0*beta))*Afbar/Abarfbar;
00092
00093
00094 rf2 = real(rf)*real(rf) + imag(rf)*imag(rf);
00095 rfbar2 = real(rfbar)*real(rfbar) + imag(rfbar)*imag(rfbar);
00096 rbarf2 = real(rbarf)*real(rbarf) + imag(rbarf)*imag(rbarf);
00097 rbarfbar2 = real(rbarfbar)*real(rbarfbar) + imag(rbarfbar)*imag(rbarfbar);
00098
00099 Af2 = real(Af)*real(Af) + imag(Af)*imag(Af);
00100 Afbar2 = real(Afbar)*real(Afbar) + imag(Afbar)*imag(Afbar);
00101 Abarf2 = real(Abarf)*real(Abarf) + imag(Abarf)*imag(Abarf);
00102 Abarfbar2 = real(Abarfbar)*real(Abarfbar) + imag(Abarfbar)*imag(Abarfbar);
00103
00104
00105
00106
00107
00108
00109 IAf = (Af2/(2*gamma_B))*(1+rf2+(1-rf2)/(1+xd*xd));
00110 IAfbar = (Afbar2/(2*gamma_B))*(1+rfbar2+(1-rfbar2)/(1+xd*xd));
00111 IAbarf = (Abarf2/(2*gamma_B))*(1+rbarf2+(1-rbarf2)/(1+xd*xd));
00112 IAbarfbar = (Abarfbar2/(2*gamma_B))*(1+rbarfbar2+(1-rbarfbar2)/(1+xd*xd));
00113
00114
00115
00116 fract = IAbarf/(IAbarf+IAf) + flip*IAbarfbar/(IAfbar+IAbarfbar);
00117
00118
00119 return;
00120 }
00121
00122 void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb, double probB0){
00123
00124
00125 static int entryCount=0;
00126 entryCount++;
00127
00128
00129 static EvtId B0B=EvtPDL::getId("anti-B0");
00130 static EvtId B0=EvtPDL::getId("B0");
00131
00132 int isB0=EvtRandom::Flat(0.0,1.0)<probB0;
00133
00134 int idaug;
00135
00136 p->setLifetime();
00137
00138
00139
00140 EvtParticle *parent=p->getParent();
00141
00142 EvtParticle *other;
00143
00144 if (parent==0) {
00145
00146
00147 p->setLifetime();
00148 t=p->getLifetime();
00149 if (p->getId()==B0) otherb=B0B;
00150 if (p->getId()==B0B) otherb=B0;
00151 entryCount--;
00152 return;
00153 }
00154 else{
00155 if (parent->getDaug(0)!=p){
00156 other=parent->getDaug(0);
00157 idaug=0;
00158 }
00159 else{
00160 other=parent->getDaug(1);
00161 idaug=1;
00162 }
00163 }
00164
00165 if (parent != 0 ) {
00166
00167
00168
00169
00170
00171
00172
00173 if ( other->getId().isAlias() ) {
00174 OtherB(p,t,otherb);
00175 entryCount--;
00176 return;
00177
00178 }
00179
00180 if (entryCount==1){
00181
00182 EvtVector4R p_init=other->getP4();
00183
00184 bool decayed = other->isDecayed();
00185
00186 other->deleteTree();
00187
00188 EvtScalarParticle* scalar_part;
00189
00190 scalar_part=new EvtScalarParticle;
00191 if (isB0) {
00192 scalar_part->init(B0,p_init);
00193 }
00194 else{
00195 scalar_part->init(B0B,p_init);
00196 }
00197 other=(EvtParticle *)scalar_part;
00198
00199 other->setDiagonalSpinDensity();
00200
00201 parent->insertDaugPtr(idaug,other);
00202
00203 if (decayed){
00204
00205 other->decay();
00206 }
00207
00208 }
00209
00210 otherb=other->getId();
00211
00212 other->setLifetime();
00213 t=p->getLifetime()-other->getLifetime();
00214
00215 otherb = other->getId();
00216
00217 }
00218 else {
00219 report(INFO,"EvtGen") << "We have an error here!!!!"<<endl;
00220 otherb = EvtId(-1,-1);
00221 }
00222
00223 entryCount--;
00224 return ;
00225 }
00226
00227
00228
00229 void EvtCPUtil::OtherB( EvtParticle *p,double &t, EvtId &otherb){
00230
00231
00232 static EvtId BSB=EvtPDL::getId("anti-B_s0");
00233 static EvtId BS0=EvtPDL::getId("B_s0");
00234 static EvtId B0B=EvtPDL::getId("anti-B0");
00235 static EvtId B0=EvtPDL::getId("B0");
00236 static EvtId UPS4=EvtPDL::getId("Upsilon(4S)");
00237
00238 if (p->getId()==BS0||p->getId()==BSB){
00239 static double ctauL=EvtPDL::getctau(EvtPDL::getId("B_s0L"));
00240 static double ctauH=EvtPDL::getctau(EvtPDL::getId("B_s0H"));
00241 static double ctau=ctauL<ctauH?ctauH:ctauL;
00242 t=-log(EvtRandom::Flat())*ctau;
00243 EvtParticle* parent=p->getParent();
00244 if (parent!=0&&(parent->getId()==BS0||parent->getId()==BSB)){
00245 if (parent->getId()==BS0) otherb=BSB;
00246 if (parent->getId()==BSB) otherb=BS0;
00247 parent->setLifetime(t);
00248 return;
00249 }
00250 if (p->getId()==BS0) otherb=BSB;
00251 if (p->getId()==BSB) otherb=BS0;
00252 p->setLifetime(t);
00253 return;
00254 }
00255
00256
00257
00258
00259 p->setLifetime();
00260
00261
00262
00263
00264 EvtParticle *parent=p->getParent();
00265
00266 if (parent==0||parent->getId()!=UPS4) {
00267
00268
00269 t=p->getLifetime();
00270 if (p->getId()==B0) otherb=B0B;
00271 if (p->getId()==B0B) otherb=B0;
00272 if (p->getId()==BS0) otherb=BSB;
00273 if (p->getId()==BSB) otherb=BS0;
00274 return;
00275 }
00276 else{
00277 if (parent->getDaug(0)!=p){
00278 otherb=parent->getDaug(0)->getId();
00279 parent->getDaug(0)->setLifetime();
00280 t=p->getLifetime()-parent->getDaug(0)->getLifetime();
00281 }
00282 else{
00283 otherb=parent->getDaug(1)->getId();
00284 parent->getDaug(1)->setLifetime();
00285 t=p->getLifetime()-parent->getDaug(1)->getLifetime();
00286 }
00287 }
00288
00289
00290 return ;
00291 }
00292
00293
00294 void EvtCPUtil::incoherentMix(const EvtId id, double &t, int &mix){
00295
00296 int stdHepNum=EvtPDL::getStdHep(id);
00297 stdHepNum=abs(stdHepNum);
00298
00299 EvtId partId=EvtPDL::evtIdFromStdHep(stdHepNum);
00300
00301 std::string partName=EvtPDL::name(partId);
00302 std::string hname=partName+std::string("H");
00303 std::string lname=partName+std::string("L");
00304
00305 EvtId lId=EvtPDL::getId(lname);
00306 EvtId hId=EvtPDL::getId(hname);
00307
00308 double ctauL=EvtPDL::getctau(lId);
00309 double ctauH=EvtPDL::getctau(hId);
00310 double ctau=0.5*(ctauL+ctauH);
00311 double y=(ctauH-ctauL)/ctau;
00312
00313
00314
00315 std::string qoverpParmName=std::string("qoverp_incohMix_")+partName;
00316 std::string mdParmName=std::string("dm_incohMix_")+partName;
00317 int ierr;
00318 double qoverp=atof(EvtSymTable::Get(qoverpParmName,ierr).c_str());
00319 double x=atof(EvtSymTable::Get(mdParmName,ierr).c_str())*ctau/EvtConst::c;
00320 double fac;
00321
00322 if(id==partId){
00323 fac=1.0/(qoverp*qoverp);
00324 }
00325 else{
00326 fac=qoverp*qoverp;
00327 }
00328
00329 double mixprob=(x*x+y*y)/(x*x+y*y+fac*(2+x*x-y*y));
00330
00331 int mixsign;
00332
00333 mixsign=(mixprob>EvtRandom::Flat(0.0,1.0))?-1:1;
00334
00335 double prob;
00336
00337 do{
00338 t=-log(EvtRandom::Flat())*ctauL;
00339 prob=1+exp(y*t/ctau)+mixsign*2*exp(0.5*y*t/ctau)*cos(x*t);
00340 }while(prob<4*EvtRandom::Flat());
00341
00342 mix=0;
00343
00344 if (mixsign==-1) mix=1;
00345
00346 return;
00347 }
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357