TEvtGen/EvtGenModels/EvtSSDCP.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   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.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 2001      Caltech
  10 //
  11 // Module: EvtSSDCP.cc
  12 //
  13 // Description: See EvtSSDCP.hh
  14 //
  15 // Modification history:
  16 //
  17 //    RYD       August 12, 2001        Module created
  18 //    F. Sandrelli, Fernando M-V  March 1, 2002     Debugged and added z parameter (CPT violation)
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtRandom.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtCPUtil.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include "EvtGenBase/EvtVector4C.hh"
  30 #include "EvtGenBase/EvtTensor4C.hh"
  31 #include "EvtGenModels/EvtSSDCP.hh"
  32 #include "EvtGenBase/EvtIncoherentMixing.hh"
  33 #include <string>
  34 #include "EvtGenBase/EvtConst.hh"
  35 using std::endl;
  36
  37 EvtSSDCP::~EvtSSDCP() {}
  38
  39 std::string EvtSSDCP::getName(){
  40
  41   return "SSD_CP";
  42
  43 }
  44
  45
  46 EvtDecayBase* EvtSSDCP::clone(){
  47
  48   return new EvtSSDCP;
  49
  50 }
  51
  52 void EvtSSDCP::init(){
  53
  54   // check that there are 8 or 12 or 14 arguments
  55
  56   checkNArg(14,12,8);
  57   checkNDaug(2);
  58
  59   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
  60   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
  61
  62   // Check it is a B0 or B0s
  63   if ( ( getParentId() != EvtPDL::getId( "B0" ) )
  64        && ( getParentId() != EvtPDL::getId( "anti-B0" ) )
  65        && ( getParentId() != EvtPDL::getId( "B_s0" ) )
  66        && ( getParentId() != EvtPDL::getId( "anti-B_s0" ) ) ) {
  67     report( ERROR , "EvtGen" ) << "EvtSSDCP only decays B0 and B0s"
  68                                << std::endl ;
  69     ::abort() ;
  70   }
  71
  72
  73   if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))||
  74        (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||(d2type==EvtSpinType::TENSOR)))||
  75        (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||(d1type==EvtSpinType::TENSOR)))
  76        ) {
  77     report(ERROR,"EvtGen") << "EvtSSDCP generator expected "
  78                            << "one of the daugters to be a scalar, the other either scalar, vector, or tensor, found:"
  79                            << EvtPDL::name(getDaug(0)).c_str()<<" and "<<EvtPDL::name(getDaug(1)).c_str()<<endl;
  80     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  81     ::abort();
  82   }
  83
  84   _dm=getArg(0)/EvtConst::c; //units of 1/mm
  85
  86   _dgog=getArg(1);
  87
  88
  89   _qoverp=getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
  90   _poverq=1.0/_qoverp;
  91
  92   _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5)));
  93
  94   _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7)));
  95
  96   if (getNArg()>=12){
  97     _eigenstate=false;
  98     _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
  99     _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
 100   }
 101   else{
 102     //I'm somewhat confused about this. For a CP eigenstate set the
 103     //amplitudes to the same. For a non CP eigenstate CPT invariance
 104     //is enforced. (ryd)
 105     if (
 106         (getDaug(0)==EvtPDL::chargeConj(getDaug(0))&&
 107          getDaug(1)==EvtPDL::chargeConj(getDaug(1)))||
 108         (getDaug(0)==EvtPDL::chargeConj(getDaug(1))&&
 109          getDaug(1)==EvtPDL::chargeConj(getDaug(0)))){
 110       _eigenstate=true;
 111     }else{
 112       _eigenstate=false;
 113       _A_fbar=conj(_Abar_f);
 114       _Abar_fbar=conj(_A_f);
 115     }
 116   }
 117
 118   //FS: new check for z
 119   if (getNArg()==14){ //FS Set _z parameter if provided else set it 0
 120     _z=EvtComplex(getArg(12),getArg(13));
 121   }
 122   else{
 123     _z=EvtComplex(0.0,0.0);
 124   }
 125
 126   // FS substituted next 2 lines...
 127
 128
 129   //
 130   //  _gamma=EvtPDL::getctau(EvtPDL::getId("B0"));  //units of 1/mm
 131   //_dgamma=_gamma*0.5*_dgog;
 132   //
 133   // ...with:
 134
 135   if ( ( getParentId() == EvtPDL::getId("B0") ) ||
 136        ( getParentId() == EvtPDL::getId("anti-B0") ) ) {
 137     _gamma=1./EvtPDL::getctau(EvtPDL::getId("B0")); //gamma/c (1/mm)
 138   }
 139   else{
 140     _gamma=1./EvtPDL::getctau(EvtPDL::getId("B_s0")) ;
 141   }
 142
 143   _dgamma=_gamma*_dgog;  //dgamma/c (1/mm)
 144
 145   if (verbose()){
 146     report(INFO,"EvtGen") << "SSD_CP will generate CP/CPT violation:"
 147                           << endl << endl
 148                           << "    " << EvtPDL::name(getParentId()).c_str() << " --> "
 149                           << EvtPDL::name(getDaug(0)).c_str() << " + "
 150                           << EvtPDL::name(getDaug(1)).c_str() << endl << endl
 151                           << "using parameters:" << endl << endl
 152                           << "  delta(m)  = " << _dm << " hbar/ps" << endl
 153                           << "dGamma      = "  << _dgamma <<" ps-1" <<endl
 154                           << "       q/p  = " << _qoverp << endl
 155                           << "        z  = " << _z << endl
 156                           << "       tau  = " << 1./_gamma << " ps" << endl;
 157   }
 158
 159 }
 160
 161 void EvtSSDCP::initProbMax() {
 162   double theProbMax =
 163     abs(_A_f) * abs(_A_f) +
 164     abs(_Abar_f) * abs(_Abar_f) +
 165     abs(_A_fbar) * abs(_A_fbar) +
 166     abs(_Abar_fbar) * abs(_Abar_fbar);
 167
 168   if (_eigenstate) theProbMax*=2;
 169
 170   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
 171   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
 172   if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10;
 173
 174   setProbMax(theProbMax);
 175 }
 176
 177 void EvtSSDCP::decay( EvtParticle *p){
 178
 179   static EvtId B0=EvtPDL::getId("B0");
 180   static EvtId B0B=EvtPDL::getId("anti-B0");
 181
 182   static EvtId B0s = EvtPDL::getId("B_s0");
 183   static EvtId B0Bs = EvtPDL::getId("anti-B_s0");
 184
 185   double t;
 186   EvtId other_b;
 187   EvtId daugs[2];
 188
 189   int flip=0;
 190   if (!_eigenstate){
 191     if (EvtRandom::Flat(0.0,1.0)<0.5) flip=1;
 192   }
 193
 194   if (!flip) {
 195     daugs[0]=getDaug(0);
 196     daugs[1]=getDaug(1);
 197   }
 198   else{
 199     daugs[0]=EvtPDL::chargeConj(getDaug(0));
 200     daugs[1]=EvtPDL::chargeConj(getDaug(1));
 201   }
 202
 203   EvtParticle *d;
 204   p->initializePhaseSpace(2, daugs);
 205
 206   EvtComplex amp;
 207
 208
 209   EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5); // t is c*Dt (mm)
 210 //  EvtIncoherentMixing::OtherB( p , t , other_b , 0.5 ) ;
 211
 212   //if (flip) t=-t;
 213
 214   //FS We assume DGamma=GammaLow-GammaHeavy and Dm=mHeavy-mLow
 215   EvtComplex expH=exp(-EvtComplex(-0.25*_dgamma*t,0.5*_dm*t));
 216   EvtComplex expL=exp(EvtComplex(-0.25*_dgamma*t,0.5*_dm*t));
 217   //FS Definition of gp and gm
 218   EvtComplex gp=0.5*(expL+expH);
 219   EvtComplex gm=0.5*(expL-expH);
 220   //FS Calculation os sqrt(1-z^2)
 221   EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
 222
 223   //EvtComplex BB=0.5*(expL+expH);                  // <B0|B0(t)>
 224   //EvtComplex barBB=_qoverp*0.5*(expL-expH);       // <B0bar|B0(t)>
 225   //EvtComplex BbarB=_poverq*0.5*(expL-expH);       // <B0|B0bar(t)>
 226   //EvtComplex barBbarB=BB;                         // <B0bar|B0bar(t)>
 227   //  FS redefinition of these guys... (See BAD #188 eq.35 for ref.)
 228   //  q/p is taken as in the BaBar Phys. Book (opposite sign wrt ref.)
 229   EvtComplex BB=gp+_z*gm;                 // <B0|B0(t)>
 230   EvtComplex barBB=sqz*_qoverp*gm;       // <B0bar|B0(t)>
 231   EvtComplex BbarB=sqz*_poverq*gm;       // <B0|B0bar(t)>
 232   EvtComplex barBbarB=gp-_z*gm;           // <B0bar|B0bar(t)>
 233
 234   if (!flip){
 235     if (other_b==B0B||other_b==B0Bs){
 236       //at t=0 we have a B0
 237       //report(INFO,"EvtGen") << "B0B"<<endl;
 238       amp=BB*_A_f+barBB*_Abar_f;
 239       //std::cout << "noflip B0B tag:"<<amp<<std::endl;
 240       //amp=0.0;
 241     }
 242     if (other_b==B0||other_b==B0s){
 243       //report(INFO,"EvtGen") << "B0"<<endl;
 244       amp=BbarB*_A_f+barBbarB*_Abar_f;
 245     }
 246   }else{
 247     if (other_b==B0||other_b==B0s){
 248       amp=BbarB*_A_fbar+barBbarB*_Abar_fbar;
 249       //std::cout << "flip B0 tag:"<<amp<<std::endl;
 250       //amp=0.0;
 251     }
 252     if (other_b==B0B||other_b==B0Bs){
 253       amp=BB*_A_fbar+barBB*_Abar_fbar;
 254     }
 255   }
 256
 257
 258   EvtVector4R p4_parent=p->getP4Restframe();
 259   double m_parent=p4_parent.mass();
 260
 261   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
 262
 263   EvtVector4R momv;
 264   EvtVector4R moms;
 265
 266   if (d2type==EvtSpinType::SCALAR){
 267     d2type=EvtPDL::getSpinType(getDaug(0));
 268     d= p->getDaug(0);
 269     momv = d->getP4();
 270     moms = p->getDaug(1)->getP4();
 271   }
 272   else{
 273     d= p->getDaug(1);
 274     momv = d->getP4();
 275     moms = p->getDaug(0)->getP4();
 276   }
 277
 278
 279
 280   if (d2type==EvtSpinType::SCALAR) {
 281     vertex(amp);
 282   }
 283
 284   if (d2type==EvtSpinType::VECTOR) {
 285
 286     double norm=momv.mass()/(momv.d3mag()*p->mass());
 287
 288     //std::cout << amp << " " << norm << " " << p4_parent << d->getP4()<< std::endl;
 289     //    std::cout << EvtPDL::name(d->getId()) << " " << EvtPDL::name(p->getDaug(0)->getId()) <<
 290     //  " 1and2 " << EvtPDL::name(p->getDaug(1)->getId()) << std::endl;
 291     //std::cout << d->eps(0) << std::endl;
 292     //std::cout << d->epsParent(0) << std::endl;
 293     vertex(0,amp*norm*p4_parent*(d->epsParent(0)));
 294     vertex(1,amp*norm*p4_parent*(d->epsParent(1)));
 295     vertex(2,amp*norm*p4_parent*(d->epsParent(2)));
 296
 297   }
 298
 299   if (d2type==EvtSpinType::TENSOR) {
 300
 301     double norm=d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag());
 302
 303
 304    vertex(0,amp*norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent);
 305    vertex(1,amp*norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent);
 306    vertex(2,amp*norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent);
 307    vertex(3,amp*norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent);
 308    vertex(4,amp*norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent);
 309
 310   }
 311
 312
 313   return ;
 314 }
 315
 316 std::string EvtSSDCP::getParamName(int i) {
 317   switch(i) {
 318   case 0:
 319     return "deltaM";
 320   case 1:
 321     return "deltaGammaOverGamma";
 322   case 2:
 323     return "qOverP";
 324   case 3:
 325     return "qOverPPhase";
 326   case 4:
 327     return "Af";
 328   case 5:
 329     return "AfPhase";
 330   case 6:
 331     return "Abarf";
 332   case 7:
 333     return "AbarfPhase";
 334   case 8:
 335     return "Afbar";
 336   case 9:
 337     return "AfbarPhase";
 338   case 10:
 339     return "Abarfbar";
 340   case 11:
 341     return "AbarfbarPhase";
 342   case 12:
 343     return "Z";
 344   case 13:
 345     return "ZPhase";
 346   default:
 347     return "";
 348   }
 349 }
 350
 351 std::string EvtSSDCP::getParamDefault(int i) {
 352   switch(i) {
 353   case 3:
 354     return "0.0";
 355   case 4:
 356     return "1.0";
 357   case 5:
 358     return "0.0";
 359   case 6:
 360     return "1.0";
 361   case 7:
 362     return "0.0";
 363   default:
 364     return "";
 365   }
 366 }