00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "EvtGenBase/EvtPatches.hh"
00025 #include <stdlib.h>
00026 #include "EvtGenBase/EvtConst.hh"
00027 #include "EvtGenBase/EvtParticle.hh"
00028 #include "EvtGenBase/EvtGenKine.hh"
00029 #include "EvtGenBase/EvtPDL.hh"
00030 #include "EvtGenBase/EvtReport.hh"
00031 #include "EvtGenBase/EvtVector4C.hh"
00032 #include "EvtGenModels/EvtVSSBMixCPT.hh"
00033 #include "EvtGenBase/EvtId.hh"
00034 #include <string>
00035 #include "EvtGenBase/EvtRandom.hh"
00036 using std::endl;
00037
00038 EvtVSSBMixCPT::~EvtVSSBMixCPT() {}
00039
00040 void EvtVSSBMixCPT::getName(std::string& model_name){
00041 model_name="VSS_BMIX";
00042 }
00043
00044
00045 EvtDecayBase* EvtVSSBMixCPT::clone(){
00046 return new EvtVSSBMixCPT;
00047 }
00048
00049 void EvtVSSBMixCPT::init(){
00050
00051
00052 if ( getNArg()>3) checkNArg(14,12,8);
00053
00054 if (getNArg()<1) {
00055 report(ERROR,"EvtGen") << "EvtVSSBMix generator expected "
00056 << " at least 1 argument (deltam) but found:"<<getNArg()<<endl;
00057 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00058 ::abort();
00059 }
00060
00061
00062 checkNDaug(2,4);
00063
00064 if ( getNDaug()==4) {
00065 if ( getDaug(0)!=getDaug(2)||getDaug(1)!=getDaug(3)){
00066 report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator allows "
00067 << " 4 daughters only if 1=3 and 2=4"
00068 << " (but 3 and 4 are aliased "<<endl;
00069 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00070 ::abort();
00071 }
00072 }
00073
00074
00075
00076 checkSpinParent(EvtSpinType::VECTOR);
00077
00078 checkSpinDaughter(0,EvtSpinType::SCALAR);
00079 checkSpinDaughter(1,EvtSpinType::SCALAR);
00080
00081
00082 if(!(EvtPDL::chargeConj(getDaug(0)) == getDaug(1))) {
00083 report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
00084 << "to be charge conjugate." << endl
00085 << " Found " << EvtPDL::name(getDaug(0)).c_str() << " and "
00086 << EvtPDL::name(getDaug(1)).c_str() << endl;
00087 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00088 ::abort();
00089 }
00090
00091 if(EvtPDL::getctau(getDaug(0)) != EvtPDL::getctau(getDaug(1))) {
00092 report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
00093 << "to have the same lifetime." << endl
00094 << " Found ctau = "
00095 << EvtPDL::getctau(getDaug(0)) << " mm and "
00096 << EvtPDL::getctau(getDaug(1)) << " mm" << endl;
00097 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00098 ::abort();
00099 }
00100
00101
00102
00103
00104 _freq= getArg(0)/EvtConst::c;
00105
00106
00107 double gamma= 1/EvtPDL::getctau(getDaug(0));
00108 _dGamma=0.0;
00109 double dgog=0.0;
00110 if ( getNArg() > 1 ) {
00111 dgog=getArg(1);
00112 _dGamma=dgog*gamma;
00113 }
00114
00115 _qoverp = EvtComplex(1.0,0.0);
00116 if ( getNArg() == 3){
00117 _qoverp = EvtComplex(getArg(2),0.0);
00118 }
00119 if ( getNArg() > 3) {
00120 _qoverp = getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
00121 }
00122 _poverq=1.0/_qoverp;
00123
00124
00125 _A_f=EvtComplex(1.0,0.0);
00126 _Abar_f=EvtComplex(0.0,0.0);
00127 _A_fbar=_Abar_f;
00128 _Abar_fbar=_A_f;
00129 if ( getNArg() > 3){
00130 _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5)));
00131 _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7)));
00132 if ( getNArg() > 8 ){
00133
00134 _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
00135 _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
00136 } else {
00137
00138 _A_fbar=_Abar_f;
00139 _Abar_fbar=_A_f;
00140 }
00141 }
00142
00143
00144 _z = EvtComplex(0.0,0.0);
00145 if ( getNArg() > 12 ){
00146 _z = EvtComplex(getArg(12),getArg(13));
00147 }
00148
00149
00150
00151 double tau= 1e12*EvtPDL::getctau(getDaug(0))/EvtConst::c;
00152 double dm= 1e-12*getArg(0);
00153 double x= dm*tau;
00154 double y= dgog*0.5;
00155 double qop2 = abs(_qoverp*_qoverp);
00156 _chib0_b0bar=qop2*(x*x+y*y)/(qop2*(x*x+y*y)+2+x*x-y*y);
00157 _chib0bar_b0=(1/qop2)*(x*x+y*y)/((1/qop2)*(x*x+y*y)+2+x*x-y*y);
00158
00159 report(INFO,"EvtGen") << "VSS_BMIXCPT will generate mixing and CPT/CP effects in mixing:"
00160 << endl << endl
00161 << " " << EvtPDL::name(getParentId()).c_str() << " --> "
00162 << EvtPDL::name(getDaug(0)).c_str() << " + "
00163 << EvtPDL::name(getDaug(1)).c_str() << endl << endl
00164 << "using parameters:" << endl << endl
00165 << " delta(m) = " << dm << " hbar/ps" << endl
00166 << " _freq = " << _freq << " hbar/mm" << endl
00167 << " dgog = " << dgog <<endl
00168 << " dGamma = " << _dGamma <<" hbar/mm" <<endl
00169 << " q/p = " << _qoverp << endl
00170 << " z = " << _z << endl
00171 << " tau = " << tau << " ps" << endl
00172 << " x = " << x << endl
00173 << " chi(B0->B0bar) = " << _chib0_b0bar << endl
00174 << " chi(B0bar->B0) = " << _chib0bar_b0 << endl
00175 << " Af = " << _A_f << endl
00176 << " Abarf = " << _Abar_f << endl
00177 << " Afbar = " << _A_fbar << endl
00178 << " Abarfbar = " << _Abar_fbar << endl
00179 << endl;
00180 }
00181
00182 void EvtVSSBMixCPT::initProbMax(){
00183
00184 setProbMax(4.0);
00185 }
00186
00187 void EvtVSSBMixCPT::decay( EvtParticle *p ){
00188
00189 static EvtId B0=EvtPDL::getId("B0");
00190 static EvtId B0B=EvtPDL::getId("anti-B0");
00191
00192
00193
00194 double rndm= EvtRandom::random();
00195
00196 if ( getNDaug()==4) {
00197 EvtId tempDaug[2];
00198
00199 if ( rndm < 0.5 ) { tempDaug[0]=getDaug(0); tempDaug[1]=getDaug(3); }
00200 else{ tempDaug[0]=getDaug(2); tempDaug[1]=getDaug(1); }
00201
00202 p->initializePhaseSpace(2,tempDaug);
00203 }
00204 else{
00205 p->initializePhaseSpace(2,getDaugs());
00206 }
00207
00208 EvtParticle *s1,*s2;
00209
00210 s1 = p->getDaug(0);
00211 s2 = p->getDaug(1);
00212
00213
00214 if ( s1->getNDaug() > 0 ) { s1->deleteDaughters();}
00215 if ( s2->getNDaug() > 0 ) { s2->deleteDaughters();}
00216
00217 EvtVector4R p1= s1->getP4();
00218 EvtVector4R p2= s2->getP4();
00219
00220
00221 rndm= EvtRandom::random();
00222 int mixed= (rndm < 0.5) ? 1 : 0;
00223
00224
00225
00226 if(mixed==1) {
00227 EvtId mixedId= (rndm < 0.25) ? getDaug(0) : getDaug(1);
00228 EvtId mixedId2= mixedId;
00229 if (getNDaug()==4&&rndm<0.25) mixedId2=getDaug(2);
00230 if (getNDaug()==4&&rndm>0.25) mixedId2=getDaug(3);
00231 s1->init(mixedId, p1);
00232 s2->init(mixedId2, p2);
00233 }
00234
00235
00236
00237
00238 if(mixed==0) {
00239 EvtId unmixedId = (rndm < 0.75) ? getDaug(0) : getDaug(1);
00240 EvtId unmixedId2= (rndm < 0.75) ? getDaug(1) : getDaug(0);
00241 if (getNDaug()==4&&rndm<0.75) unmixedId2=getDaug(3);
00242 if (getNDaug()==4&&rndm>0.75) unmixedId2=getDaug(2);
00243 s1->init(unmixedId, p1);
00244 s2->init(unmixedId2, p2);
00245 }
00246
00247
00248
00249
00250 s1->setLifetime();
00251 s2->setLifetime();
00252 double dct= s1->getLifetime() - s2->getLifetime();
00253
00254
00255 EvtComplex exp1(-0.25*_dGamma*dct,0.5*_freq*dct);
00256
00257
00258
00259
00260
00261
00262
00263
00264 EvtId stateAtDeltaTeq0 = (s2->getId()==B0) ? B0B : B0;
00265
00266
00267 EvtComplex osc_amp;
00268
00269
00270 EvtComplex gp=0.5*(exp(-1.0*exp1)+exp(exp1));
00271 EvtComplex gm=0.5*(exp(-1.0*exp1)-exp(exp1));
00272 EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
00273
00274 EvtComplex BB=gp+_z*gm;
00275 EvtComplex barBB=-sqz*_qoverp*gm;
00276 EvtComplex BbarB=-sqz*_poverq*gm;
00277 EvtComplex barBbarB=gp-_z*gm;
00278
00279
00280 if ( !mixed&&stateAtDeltaTeq0==B0 ) {
00281 osc_amp= BB*_A_f+barBB*_Abar_f;
00282 }
00283 if ( !mixed&&stateAtDeltaTeq0==B0B ) {
00284 osc_amp= barBbarB*_Abar_fbar+BbarB*_A_fbar;
00285 }
00286
00287 if ( mixed&&stateAtDeltaTeq0==B0 ) {
00288 osc_amp=barBB*_Abar_fbar+BB*_A_fbar;
00289 }
00290 if ( mixed&&stateAtDeltaTeq0==B0B ) {
00291 osc_amp=BbarB*_A_f+barBbarB*_Abar_f;
00292 }
00293
00294
00295 double norm=1.0/p1.d3mag();
00296 vertex(0,norm*osc_amp*p1*(p->eps(0)));
00297 vertex(1,norm*osc_amp*p1*(p->eps(1)));
00298 vertex(2,norm*osc_amp*p1*(p->eps(2)));
00299
00300 return ;
00301 }
00302
00303
00304
00305
00306
00307
00308