00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <stdlib.h>
00023 #include "EvtGenBase/EvtParticle.hh"
00024 #include "EvtGenBase/EvtRandom.hh"
00025 #include "EvtGenBase/EvtGenKine.hh"
00026 #include "EvtGenBase/EvtCPUtil.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenBase/EvtReport.hh"
00029 #include "EvtGenBase/EvtVector4C.hh"
00030 #include "EvtGenBase/EvtTensor4C.hh"
00031 #include "EvtGenModels/EvtSSDCP.hh"
00032 #include <string>
00033 #include "EvtGenBase/EvtConst.hh"
00034 using std::endl;
00035
00036 EvtSSDCP::~EvtSSDCP() {}
00037
00038 void EvtSSDCP::getName(std::string& model_name){
00039
00040 model_name="SSD_CP";
00041
00042 }
00043
00044
00045 EvtDecayBase* EvtSSDCP::clone(){
00046
00047 return new EvtSSDCP;
00048
00049 }
00050
00051 void EvtSSDCP::init(){
00052
00053
00054
00055 checkNArg(14,12,8);
00056 checkNDaug(2);
00057
00058 EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
00059 EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))||
00072 (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||(d2type==EvtSpinType::TENSOR)))||
00073 (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||(d1type==EvtSpinType::TENSOR)))
00074 ) {
00075 report(ERROR,"EvtGen") << "EvtSSDCP generator expected "
00076 << "one of the daugters to be a scalar, the other either scalar, vector, or tensor, found:"
00077 << EvtPDL::name(getDaug(0)).c_str()<<" and "<<EvtPDL::name(getDaug(1)).c_str()<<endl;
00078 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00079 ::abort();
00080 }
00081
00082 _dm=getArg(0)/EvtConst::c;
00083
00084 _dgog=getArg(1);
00085
00086
00087 _qoverp=getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
00088 _poverq=1.0/_qoverp;
00089
00090 _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5)));
00091
00092 _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7)));
00093
00094 if (getNArg()>=12){
00095 _eigenstate=false;
00096 _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
00097 _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
00098 }
00099 else{
00100
00101
00102
00103 if (
00104 (getDaug(0)==EvtPDL::chargeConj(getDaug(0))&&
00105 getDaug(1)==EvtPDL::chargeConj(getDaug(1)))||
00106 (getDaug(0)==EvtPDL::chargeConj(getDaug(1))&&
00107 getDaug(1)==EvtPDL::chargeConj(getDaug(0)))){
00108 _eigenstate=true;
00109 }else{
00110 _eigenstate=false;
00111 _A_fbar=conj(_Abar_f);
00112 _Abar_fbar=conj(_A_f);
00113 }
00114 }
00115
00116
00117 if (getNArg()==14){
00118 _z=EvtComplex(getArg(12),getArg(13));
00119 }
00120 else{
00121 _z=EvtComplex(0.0,0.0);
00122 }
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 _gamma=1/EvtPDL::getctau(EvtPDL::getId("B0"));
00134 _dgamma=_gamma*_dgog;
00135
00136 if (verbose()){
00137 report(INFO,"EvtGen") << "SSD_CP will generate CP/CPT violation:"
00138 << endl << endl
00139 << " " << EvtPDL::name(getParentId()).c_str() << " --> "
00140 << EvtPDL::name(getDaug(0)).c_str() << " + "
00141 << EvtPDL::name(getDaug(1)).c_str() << endl << endl
00142 << "using parameters:" << endl << endl
00143 << " delta(m) = " << _dm << " hbar/ps" << endl
00144 << "dGamma = " << _dgamma <<" ps-1" <<endl
00145 << " q/p = " << _qoverp << endl
00146 << " z = " << _z << endl
00147 << " tau = " << 1./_gamma << " ps" << endl;
00148 }
00149
00150 }
00151
00152 void EvtSSDCP::initProbMax() {
00153 double theProbMax =
00154 abs(_A_f) * abs(_A_f) +
00155 abs(_Abar_f) * abs(_Abar_f) +
00156 abs(_A_fbar) * abs(_A_fbar) +
00157 abs(_Abar_fbar) * abs(_Abar_fbar);
00158
00159 if (_eigenstate) theProbMax*=2;
00160
00161 EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
00162 EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
00163 if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10;
00164
00165 setProbMax(theProbMax);
00166 }
00167
00168 void EvtSSDCP::decay( EvtParticle *p){
00169
00170 static EvtId B0=EvtPDL::getId("B0");
00171 static EvtId B0B=EvtPDL::getId("anti-B0");
00172
00173 double t;
00174 EvtId other_b;
00175 EvtId daugs[2];
00176
00177 int flip=0;
00178 if (!_eigenstate){
00179 if (EvtRandom::Flat(0.0,1.0)<0.5) flip=1;
00180 }
00181
00182 if (!flip) {
00183 daugs[0]=getDaug(0);
00184 daugs[1]=getDaug(1);
00185 }
00186 else{
00187 daugs[0]=EvtPDL::chargeConj(getDaug(0));
00188 daugs[1]=EvtPDL::chargeConj(getDaug(1));
00189 }
00190
00191 EvtParticle *d;
00192 p->initializePhaseSpace(2, daugs);
00193
00194 EvtComplex amp;
00195
00196
00197 EvtCPUtil::OtherB(p,t,other_b,0.5);
00198
00199
00200
00201
00202 EvtComplex expL=exp(-EvtComplex(-0.25*_dgamma*t,0.5*_dm*t));
00203 EvtComplex expH=exp(EvtComplex(-0.25*_dgamma*t,0.5*_dm*t));
00204
00205 EvtComplex gp=0.5*(expL+expH);
00206 EvtComplex gm=0.5*(expL-expH);
00207
00208 EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
00209
00210
00211
00212
00213
00214
00215
00216 EvtComplex BB=gp+_z*gm;
00217 EvtComplex barBB=sqz*_qoverp*gm;
00218 EvtComplex BbarB=sqz*_poverq*gm;
00219 EvtComplex barBbarB=gp-_z*gm;
00220
00221 if (!flip){
00222 if (other_b==B0B){
00223
00224
00225 amp=BB*_A_f+BbarB*_Abar_f;
00226
00227
00228 }
00229 if (other_b==B0){
00230
00231 amp=barBB*_A_f+barBbarB*_Abar_f;
00232 }
00233 }else{
00234 if (other_b==B0){
00235 amp=barBB*_A_fbar+barBbarB*_Abar_fbar;
00236
00237
00238 }
00239 if (other_b==B0B){
00240 amp=BB*_A_fbar+BbarB*_Abar_fbar;
00241 }
00242 }
00243
00244
00245 EvtVector4R p4_parent=p->getP4Restframe();
00246 double m_parent=p4_parent.mass();
00247
00248 EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
00249
00250 EvtVector4R momv;
00251 EvtVector4R moms;
00252
00253 if (d2type==EvtSpinType::SCALAR){
00254 d2type=EvtPDL::getSpinType(getDaug(0));
00255 d= p->getDaug(0);
00256 momv = d->getP4();
00257 moms = p->getDaug(1)->getP4();
00258 }
00259 else{
00260 d= p->getDaug(1);
00261 momv = d->getP4();
00262 moms = p->getDaug(0)->getP4();
00263 }
00264
00265
00266
00267 if (d2type==EvtSpinType::SCALAR) {
00268 vertex(amp);
00269 }
00270
00271 if (d2type==EvtSpinType::VECTOR) {
00272
00273 double norm=momv.mass()/(momv.d3mag()*p->mass());
00274
00275 vertex(0,amp*norm*p4_parent*(d->epsParent(0)));
00276 vertex(1,amp*norm*p4_parent*(d->epsParent(1)));
00277 vertex(2,amp*norm*p4_parent*(d->epsParent(2)));
00278
00279 }
00280
00281 if (d2type==EvtSpinType::TENSOR) {
00282
00283 double norm=d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag());
00284
00285
00286 vertex(0,amp*norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent);
00287 vertex(1,amp*norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent);
00288 vertex(2,amp*norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent);
00289 vertex(3,amp*norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent);
00290 vertex(4,amp*norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent);
00291
00292 }
00293
00294
00295 return ;
00296 }
00297