TEvtGen/EvtGenModels/EvtSVVHelCPMix.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) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtSVVHelCPMix.cc
  12 //
  13 // Description: Routine to decay scalar -> 2 vectors
  14 //              by specifying the helicity amplitudes, taking appropriate
  15 //              weak phases into account to get mixing and CP violation through
  16 //              interference. Based on EvtSVVHelAmp. Particularly appropriate for
  17 //              Bs->J/Psi+Phi
  18 //
  19 // Modification history:
  20 //
  21 //    RYD       November 24, 1996       EvtSVVHelAmp Module
  22 //    CATMORE   March 2004              Modified to EvtSVVHelCPMix
  23 //
  24 // Model takes the following as user-specified arguments:
  25 //      deltaM, averageM - mass diference and average of light and heavy mass eigenstates (real scalars)
  26 //      gamma, deltagamma - average width and width difference of the l and h eigenstates (real scalars)
  27 //      delta1, delta2 - strong phases (real scalars)
  28 //      direct weak phase (real scalar) (for Bs->JPsiPhi this will be zero)
  29 //      weak mixing phase (real scalar) (this is equal to 2*arg(Vts Vtb) for Bs->JPsiPhi)
  30 //      Magnitudes of helicity amplitudes as in SVV_HELAMP
  31 // See Phys Rev D 34 p1404 - p1417 and chapters 5 and 7 of Physics Reports 370 p537-680 for more details
  32 //------------------------------------------------------------------------
  33 //
  34 #ifdef WIN32
  35   #pragma warning( disable : 4786 )
  36   // Disable anoying warning about symbol size
  37 #endif
  38 #include <iostream>
  39 #include <fstream>
  40 #include <stdlib.h>
  41 #include <ctype.h>
  42 #include "EvtGenBase/EvtParticle.hh"
  43 #include "EvtGenBase/EvtGenKine.hh"
  44 #include "EvtGenBase/EvtPDL.hh"
  45 #include "EvtGenBase/EvtVector4C.hh"
  46 #include "EvtGenBase/EvtTensor4C.hh"
  47 #include "EvtGenBase/EvtVector3C.hh"
  48 #include "EvtGenBase/EvtVector3R.hh"
  49 #include "EvtGenBase/EvtTensor3C.hh"
  50 #include "EvtGenBase/EvtReport.hh"
  51 #include "EvtGenModels/EvtSVVHelCPMix.hh"
  52 #include "EvtGenBase/EvtId.hh"
  53 #include <string>
  54
  55 EvtSVVHelCPMix::~EvtSVVHelCPMix() {}
  56
  57 std::string EvtSVVHelCPMix::getName(){
  58
  59   return "SVVHELCPMIX";
  60
  61 }
  62
  63
  64 EvtDecayBase* EvtSVVHelCPMix::clone(){
  65
  66   return new EvtSVVHelCPMix;
  67
  68 }
  69
  70 void EvtSVVHelCPMix::init(){
  71
  72   // check that there are 12 arguments
  73   checkNArg(12);
  74   checkNDaug(2);
  75
  76   checkSpinParent(EvtSpinType::SCALAR);
  77
  78   checkSpinDaughter(0,EvtSpinType::VECTOR);
  79   checkSpinDaughter(1,EvtSpinType::VECTOR);
  80
  81   hp = EvtComplex(getArg(0)*cos(getArg(1)),getArg(0)*sin(getArg(1)));
  82   h0 = EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
  83   hm = EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
  84   averageM = getArg(6);
  85   deltaM = getArg(7);
  86   gamma = getArg(8);
  87   deltagamma = getArg(9);
  88   weakmixingphase = EvtComplex(cos(getArg(10)),sin(getArg(10)));
  89   weakdirectphase = EvtComplex(cos(getArg(11)),sin(getArg(11)));
  90 }
  91
  92
  93 void EvtSVVHelCPMix::initProbMax(){
  94
  95   setProbMax(getArg(0)*getArg(0)+getArg(2)*getArg(2)+getArg(4)*getArg(4));
  96
  97 }
  98
  99
 100 void EvtSVVHelCPMix::decay( EvtParticle *p){
 101
 102   EvtParticle* parent = p;
 103   EvtAmp& amp = _amp2;
 104   EvtId n_v1 = getDaug(0);
 105   EvtId n_v2 = getDaug(1);
 106
 107   //  Routine to decay a vector into a vector and scalar.  Started
 108   //  by ryd on Oct 17, 1996.
 109   // Modified by J.Catmore to take account of CP-violation and mixing
 110
 111   int tndaug = 2;
 112   EvtId tdaug[2];
 113   EvtId Bs=EvtPDL::getId("B_s0");
 114   EvtId antiBs=EvtPDL::getId("anti-B_s0");
 115   tdaug[0] = n_v1;
 116   tdaug[1] = n_v2;
 117
 118 // Phase space and kinematics
 119
 120   parent->initializePhaseSpace(tndaug,tdaug);
 121
 122   EvtParticle *v1,*v2;
 123   v1 = parent->getDaug(0);
 124   v2 = parent->getDaug(1);
 125
 126   EvtVector4R momv1 = v1->getP4();
 127   EvtVector4R momv2 = v2->getP4();
 128
 129   EvtVector3R v1dir(momv1.get(1),momv1.get(2),momv1.get(3));
 130   v1dir=v1dir/v1dir.d3mag();
 131
 132 // Definition of quantities used in construction of complex amplitudes:
 133
 134 EvtTensor3C M;  // Tensor as defined in EvtGen manual, equ 117
 135 EvtComplex a,b,c; // Helicity amplitudes; EvtGen manual eqns 126-128, also see Phys Lett B 369 p144-150 eqn 15
 136 EvtComplex deltamu = EvtComplex(deltaM, -0.5*deltagamma); // See Phys Rev D 34 p1404
 137
 138 // conversion from times in mm/c to natural units [GeV]^-1
 139 double t = ((parent->getLifetime())/2.998e11)*6.58e-25;
 140
 141 // The following two quantities defined in Phys Rev D 34 p1404
 142 EvtComplex fplus =  EvtComplex(cos(averageM*t),-1.*sin(averageM*t))*exp(-(gamma/2.0)*t)*
 143     (cos(0.5*deltaM*t)*cosh(0.25*deltagamma*t)+EvtComplex(0.0,sin(0.5*deltaM*t)*sinh(0.25*deltagamma*t)));
 144 EvtComplex fminus = EvtComplex(cos(averageM*t), -1.*sin(averageM*t))*exp(-(gamma/2.0)*t)*EvtComplex(0.0,1.0)*
 145     (sin(0.5*deltaM*t)*cosh(0.25*deltagamma*t)-EvtComplex(0.0,1.0)*sinh(0.25*deltagamma*t)*cos(0.5*deltaM*t));
 146
 147 // See EvtGen manual pp 106-107
 148
 149   a=-0.5*(hp+hm);
 150   b=EvtComplex(0.0,0.5)*(hp-hm);
 151   c=(h0+0.5*(hp+hm));
 152
 153   M=a*EvtTensor3C::id()+
 154                 b*eps(v1dir)+
 155                 c*directProd(v1dir,v1dir);
 156
 157   EvtVector3C t0=M.cont1(v1->eps(0).vec().conj());
 158   EvtVector3C t1=M.cont1(v1->eps(1).vec().conj());
 159   EvtVector3C t2=M.cont1(v1->eps(2).vec().conj());
 160
 161   EvtVector3C eps0=v2->eps(0).vec().conj();
 162   EvtVector3C eps1=v2->eps(1).vec().conj();
 163   EvtVector3C eps2=v2->eps(2).vec().conj();
 164
 165 // We need two sets of equations, one for mesons which were in the Bs state at t=0, and another
 166 // for those which were in the antiBs state. Each equation consists of a sum of amplitudes - mod-squaring gives the interference terms.
 167
 168         EvtComplex amplSum00, amplSum01, amplSum02;
 169         EvtComplex amplSum10, amplSum11, amplSum12;
 170         EvtComplex amplSum20, amplSum21, amplSum22;
 171
 172   // First the Bs state:
 173
 174   if (parent->getId()==Bs) {
 175
 176         amplSum00 = (fplus*weakdirectphase*t0*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps0);
 177         amplSum01 = (fplus*weakdirectphase*t0*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps1);
 178         amplSum02 = (fplus*weakdirectphase*t0*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps2);
 179
 180         amplSum10 = (fplus*weakdirectphase*t1*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps0);
 181         amplSum11 = (fplus*weakdirectphase*t1*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps1);
 182         amplSum12 = (fplus*weakdirectphase*t1*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps2);
 183
 184         amplSum20 = (fplus*weakdirectphase*t2*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps0);
 185         amplSum21 = (fplus*weakdirectphase*t2*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps1);
 186         amplSum22 = (fplus*weakdirectphase*t2*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps2);
 187   }
 188
 189   // Now the anti-Bs state:
 190
 191   if (parent->getId()==antiBs) {
 192
 193         amplSum00 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps0) + (fplus*(1.0/weakdirectphase)*t0*eps0);
 194         amplSum01 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps1) + (fplus*(1.0/weakdirectphase)*t0*eps1);
 195         amplSum02 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps2) + (fplus*(1.0/weakdirectphase)*t0*eps2);
 196
 197         amplSum10 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps0) + (fplus*(1.0/weakdirectphase)*t1*eps0);
 198         amplSum11 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps1) + (fplus*(1.0/weakdirectphase)*t1*eps1);
 199         amplSum12 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps2) + (fplus*(1.0/weakdirectphase)*t1*eps2);
 200
 201         amplSum20 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps0) + (fplus*(1.0/weakdirectphase)*t2*eps0);
 202         amplSum21 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps1) + (fplus*(1.0/weakdirectphase)*t2*eps1);
 203         amplSum22 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps2) + (fplus*(1.0/weakdirectphase)*t2*eps2);
 204
 205   }
 206
 207   // Now set the amplitude
 208
 209   amp.vertex(0,0,amplSum00);
 210   report(INFO,"EvtGen") << "00: " << amplSum00 << std::endl;
 211   amp.vertex(0,1,amplSum01);
 212   report(INFO,"EvtGen") << "01: " << amplSum01 << std::endl;
 213   amp.vertex(0,2,amplSum02);
 214   report(INFO,"EvtGen") << "02: " << amplSum02 << std::endl;
 215
 216   amp.vertex(1,0,amplSum10);
 217   report(INFO,"EvtGen") << "10: " << amplSum10 << std::endl;
 218   amp.vertex(1,1,amplSum11);
 219   report(INFO,"EvtGen") << "11: " << amplSum11 << std::endl;
 220   amp.vertex(1,2,amplSum12);
 221   report(INFO,"EvtGen") << "12: " << amplSum12 << std::endl;
 222
 223   amp.vertex(2,0,amplSum20);
 224   report(INFO,"EvtGen") << "20: " << amplSum20 << std::endl;
 225   amp.vertex(2,1,amplSum21);
 226   report(INFO,"EvtGen") << "21: " << amplSum21 << std::endl;
 227   amp.vertex(2,2,amplSum22);
 228   report(INFO,"EvtGen") << "22: " << amplSum22 << std::endl;
 229
 230   return ;
 231
 232 }
 233
 234