TEvtGen/EvtGenModels/EvtVSSBMixCPT.cxx

   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) 2002      INFN-Pisa
  10 //
  11 // Module: EvtVSSBMixCPT.cc
  12 //
  13 // Description:
  14 //    Routine to decay vector-> scalar scalar with coherent BB-like mixing
  15 //                              including CPT effects
  16 //    Based on VSSBMIX
  17 //
  18 // Modification history:
  19 //
  20 //    F. Sandrelli, Fernando M-V March 03, 2002
  21 //
  22 //------------------------------------------------------------------------
  23 //
  24 #include "EvtGenBase/EvtPatches.hh"
  25 #include <stdlib.h>
  26 #include "EvtGenBase/EvtConst.hh"
  27 #include "EvtGenBase/EvtParticle.hh"
  28 #include "EvtGenBase/EvtGenKine.hh"
  29 #include "EvtGenBase/EvtPDL.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include "EvtGenBase/EvtVector4C.hh"
  32 #include "EvtGenModels/EvtVSSBMixCPT.hh"
  33 #include "EvtGenBase/EvtId.hh"
  34 #include <string>
  35 #include "EvtGenBase/EvtRandom.hh"
  36 using std::endl;
  37
  38 EvtVSSBMixCPT::~EvtVSSBMixCPT() {}
  39
  40 std::string EvtVSSBMixCPT::getName(){
  41   return "VSS_BMIX";
  42 }
  43
  44
  45 EvtDecayBase* EvtVSSBMixCPT::clone(){
  46   return new EvtVSSBMixCPT;
  47 }
  48
  49 void EvtVSSBMixCPT::init(){
  50
  51   if ( getNArg()>4) checkNArg(14,12,8);
  52
  53   if (getNArg()<1) {
  54     report(ERROR,"EvtGen") << "EvtVSSBMix generator expected "
  55                            << " at least 1 argument (deltam) but found:"<<getNArg()<<endl;
  56     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  57     ::abort();
  58   }
  59   // check that we are asked to produced exactly 2 daughters
  60   //4 are allowed if they are aliased..
  61   checkNDaug(2,4);
  62
  63   if ( getNDaug()==4) {
  64     if ( getDaug(0)!=getDaug(2)||getDaug(1)!=getDaug(3)){
  65       report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator allows "
  66                              << " 4 daughters only if 1=3 and 2=4"
  67                              << " (but 3 and 4 are aliased "<<endl;
  68       report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  69       ::abort();
  70     }
  71   }
  72   // check that we are asked to decay a vector particle into a pair
  73   // of scalar particles
  74
  75   checkSpinParent(EvtSpinType::VECTOR);
  76
  77   checkSpinDaughter(0,EvtSpinType::SCALAR);
  78   checkSpinDaughter(1,EvtSpinType::SCALAR);
  79
  80   // check that our daughter particles are charge conjugates of each other
  81   if(!(EvtPDL::chargeConj(getDaug(0)) == getDaug(1))) {
  82     report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
  83                            << "to be charge conjugate." << endl
  84                            << "  Found " << EvtPDL::name(getDaug(0)).c_str() << " and "
  85                            << EvtPDL::name(getDaug(1)).c_str() << endl;
  86     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  87     ::abort();
  88   }
  89   // check that both daughter particles have the same lifetime
  90   if(EvtPDL::getctau(getDaug(0)) != EvtPDL::getctau(getDaug(1))) {
  91     report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
  92                            << "to have the same lifetime." << endl
  93                            << "  Found ctau = "
  94                            << EvtPDL::getctau(getDaug(0)) << " mm and "
  95                            << EvtPDL::getctau(getDaug(1)) << " mm" << endl;
  96     report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  97     ::abort();
  98   }
  99   // precompute quantities that will be used to generate events
 100   // and print out a summary of parameters for this decay
 101
 102   // mixing frequency in hbar/mm
 103   _freq= getArg(0)/EvtConst::c;
 104
 105   // deltaG
 106   double gamma= 1/EvtPDL::getctau(getDaug(0));  // gamma/c (1/mm)
 107   _dGamma=0.0;
 108   double dgog=0.0;
 109   if ( getNArg() > 1 ) {
 110     dgog=getArg(1);
 111     _dGamma=dgog*gamma;
 112   }
 113   // q/p
 114   _qoverp = EvtComplex(1.0,0.0);
 115   if ( getNArg() > 2){
 116     _qoverp = EvtComplex(getArg(2),0.0);
 117   }
 118   if ( getNArg() > 3) {
 119     _qoverp = getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
 120   }
 121   _poverq=1.0/_qoverp;
 122
 123   // decay amplitudes
 124   _A_f=EvtComplex(1.0,0.0);
 125   _Abar_f=EvtComplex(0.0,0.0);
 126   _A_fbar=_Abar_f;  // CPT conservation
 127   _Abar_fbar=_A_f;  // CPT conservation
 128   if ( getNArg() > 4){
 129     _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5)));    // this allows for DCSD
 130     _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7))); // this allows for DCSD
 131     if ( getNArg() > 8 ){
 132       // CPT violation in decay
 133       _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
 134       _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
 135     } else {
 136       // CPT conservation in decay
 137       _A_fbar=_Abar_f;
 138       _Abar_fbar=_A_f;
 139     }
 140   }
 141
 142   // CPT violation in mixing
 143   _z = EvtComplex(0.0,0.0);
 144   if ( getNArg() > 12 ){
 145     _z = EvtComplex(getArg(12),getArg(13));
 146   }
 147
 148
 149   // some printout
 150   double tau= 1e12*EvtPDL::getctau(getDaug(0))/EvtConst::c; // in ps
 151   double dm= 1e-12*getArg(0); // B0/anti-B0 mass difference in hbar/ps
 152   double x= dm*tau;
 153   double y= dgog*0.5; //y=dgamma/(2*gamma)
 154   double qop2 = abs(_qoverp*_qoverp);
 155   _chib0_b0bar=qop2*(x*x+y*y)/(qop2*(x*x+y*y)+2+x*x-y*y);  // does not include CPT in mixing
 156   _chib0bar_b0=(1/qop2)*(x*x+y*y)/((1/qop2)*(x*x+y*y)+2+x*x-y*y); // does not include CPT in mixing
 157
 158   if ( verbose() ) {
 159     report(INFO,"EvtGen") << "VSS_BMIXCPT will generate mixing and CPT/CP effects in mixing:"
 160                           << endl << endl
 161                           << "    " << EvtPDL::name(getParentId()).c_str() << " --> "
 162                           << EvtPDL::name(getDaug(0)).c_str() << " + "
 163                           << EvtPDL::name(getDaug(1)).c_str() << endl << endl
 164                           << "using parameters:" << endl << endl
 165                           << "  delta(m)  = " << dm << " hbar/ps" << endl
 166                           << "  _freq     = " << _freq << " hbar/mm" << endl
 167                           << "  dgog      = "  << dgog <<endl
 168                           << "  dGamma    = "  << _dGamma <<" hbar/mm" <<endl
 169                           << "  q/p       = " << _qoverp << endl
 170                           << "  z         = " << _z << endl
 171                           << "  tau       = " << tau << " ps" << endl
 172                           << "  x         = " << x << endl
 173                           << " chi(B0->B0bar) = " << _chib0_b0bar << endl
 174                           << " chi(B0bar->B0) = " << _chib0bar_b0 << endl
 175                           << " Af         = " << _A_f << endl
 176                           << " Abarf      = " << _Abar_f << endl
 177                           << " Afbar      = " << _A_fbar << endl
 178                           << " Abarfbar   = " << _Abar_fbar << endl
 179                           << endl;
 180   }
 181 }
 182
 183 void EvtVSSBMixCPT::initProbMax(){
 184   // this value is ok for reasonable values of all the parameters
 185   setProbMax(4.0);
 186 }
 187
 188 void EvtVSSBMixCPT::decay( EvtParticle *p ){
 189
 190   static EvtId B0=EvtPDL::getId("B0");
 191   static EvtId B0B=EvtPDL::getId("anti-B0");
 192
 193   // generate a final state according to phase space
 194
 195   double rndm= EvtRandom::random();
 196
 197   if ( getNDaug()==4) {
 198     EvtId tempDaug[2];
 199
 200     if ( rndm < 0.5 ) { tempDaug[0]=getDaug(0);  tempDaug[1]=getDaug(3); }
 201     else{ tempDaug[0]=getDaug(2);  tempDaug[1]=getDaug(1); }
 202
 203     p->initializePhaseSpace(2,tempDaug);
 204   }
 205   else{ //nominal case.
 206     p->initializePhaseSpace(2,getDaugs());
 207   }
 208
 209   EvtParticle *s1,*s2;
 210
 211   s1 = p->getDaug(0);
 212   s2 = p->getDaug(1);
 213   //delete any daughters - if there are daughters, they
 214   //are from the initialization and will be redone later
 215   if ( s1->getNDaug() > 0 ) { s1->deleteDaughters();}
 216   if ( s2->getNDaug() > 0 ) { s2->deleteDaughters();}
 217
 218   EvtVector4R p1= s1->getP4();
 219   EvtVector4R p2= s2->getP4();
 220
 221   // throw a random number to decide if this final state should be mixed
 222   rndm= EvtRandom::random();
 223   int mixed= (rndm < 0.5) ? 1 : 0;
 224
 225   // if this decay is mixed, choose one of the 2 possible final states
 226   // with equal probability (re-using the same random number)
 227   if(mixed==1) {
 228     EvtId mixedId= (rndm < 0.25) ? getDaug(0) : getDaug(1);
 229     EvtId mixedId2= mixedId;
 230     if (getNDaug()==4&&rndm<0.25)  mixedId2=getDaug(2);
 231     if (getNDaug()==4&&rndm>0.25)  mixedId2=getDaug(3);
 232     s1->init(mixedId, p1);
 233     s2->init(mixedId2, p2);
 234   }
 235
 236
 237   // if this decay is unmixed, choose one of the 2 possible final states
 238   // with equal probability (re-using the same random number)
 239   if(mixed==0) {
 240     EvtId unmixedId = (rndm < 0.75) ? getDaug(0) : getDaug(1);
 241     EvtId unmixedId2= (rndm < 0.75) ? getDaug(1) : getDaug(0);
 242     if (getNDaug()==4&&rndm<0.75)  unmixedId2=getDaug(3);
 243     if (getNDaug()==4&&rndm>0.75)  unmixedId2=getDaug(2);
 244     s1->init(unmixedId, p1);
 245     s2->init(unmixedId2, p2);
 246   }
 247
 248   // choose a decay time for each final state particle using the
 249   // lifetime (which must be the same for both particles) in pdt.table
 250   // and calculate the lifetime difference for this event
 251   s1->setLifetime();
 252   s2->setLifetime();
 253   double dct= s1->getLifetime() - s2->getLifetime(); // in mm
 254
 255   // Convention: _dGamma=GammaLight-GammaHeavy
 256   EvtComplex exp1(-0.25*_dGamma*dct,0.5*_freq*dct);
 257
 258   /*
 259   //Find the flavor of the B that decayed first.
 260   EvtId firstDec = (dct > 0 ) ? s2->getId() : s1->getId();
 261
 262   //This tags the flavor of the other particle at that time.
 263   EvtId stateAtDeltaTeq0 = ( firstDec==B0 ) ? B0B : B0;
 264   */
 265   EvtId stateAtDeltaTeq0 = (s2->getId()==B0) ? B0B : B0;
 266
 267   // calculate the oscillation amplitude, based on wether this event is mixed or not
 268   EvtComplex osc_amp;
 269
 270   //define some useful functions: (see BAD #188 eq. 39 for ref.)
 271   EvtComplex gp=0.5*(exp(-1.0*exp1)+exp(exp1));
 272   EvtComplex gm=0.5*(exp(-1.0*exp1)-exp(exp1));
 273   EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
 274
 275   EvtComplex BB=gp+_z*gm;                // <B0|B0(t)>
 276   EvtComplex barBB=-sqz*_qoverp*gm;      // <B0bar|B0(t)>
 277   EvtComplex BbarB=-sqz*_poverq*gm;      // <B0|B0bar(t)>
 278   EvtComplex barBbarB=gp-_z*gm;          // <B0bar|B0bar(t)>
 279
 280   //
 281   if ( !mixed&&stateAtDeltaTeq0==B0 ) {
 282     osc_amp= BB*_A_f+barBB*_Abar_f;
 283   }
 284   if ( !mixed&&stateAtDeltaTeq0==B0B ) {
 285     osc_amp= barBbarB*_Abar_fbar+BbarB*_A_fbar;
 286   }
 287
 288   if ( mixed&&stateAtDeltaTeq0==B0 ) {
 289     osc_amp=barBB*_Abar_fbar+BB*_A_fbar;
 290   }
 291   if ( mixed&&stateAtDeltaTeq0==B0B ) {
 292     osc_amp=BbarB*_A_f+barBbarB*_Abar_f;
 293   }
 294
 295   // store the amplitudes for each parent spin basis state
 296   double norm=1.0/p1.d3mag();
 297   vertex(0,norm*osc_amp*p1*(p->eps(0)));
 298   vertex(1,norm*osc_amp*p1*(p->eps(1)));
 299   vertex(2,norm*osc_amp*p1*(p->eps(2)));
 300
 301   return ;
 302 }
 303
 304
 305
 306
 307
 308
 309