TEvtGen/EvtGen/EvtGenModels/EvtBBScalar.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) 2003      Caltech
  10 //
  11 // Module: EvtGen/EvtBBScalar
  12 //
  13 // Description:Implementation of the decay B- -> lambda p_bar pi according to
  14 // hep-ph/0204185, hep-ph/0211240
  15 // This model is intended to be applicable to all decays of the type B-> baryon baryon scalar
  16 //
  17 // Modification history:
  18 //
  19 //    Jan Strube     March 24, 2006         Module created
  20 //
  21 //------------------------------------------------------------------------
  22 #include "EvtGenBase/EvtPatches.hh"
  23
  24 #include "EvtGenModels/EvtBBScalar.hh"
  25 #include "EvtGenBase/EvtGammaMatrix.hh"
  26 #include "EvtGenBase/EvtDiracSpinor.hh"
  27 #include "EvtGenBase/EvtSpinType.hh"
  28 #include "EvtGenBase/EvtTensor4C.hh"
  29 #include <cmath>
  30
  31 using namespace std;
  32
  33 const float pi = 3.14159;
  34 const EvtComplex EvtBBScalar::I = EvtComplex(0, 1);
  35 const EvtComplex EvtBBScalar::V_ub = EvtComplex(3.67e-3*cos(60/180*pi), 3.67e-3*cos(60/180*pi));
  36 const EvtComplex EvtBBScalar::V_us_star = EvtComplex(0.22, 0);
  37 const EvtComplex EvtBBScalar::a1 = EvtComplex(1.05, 0);
  38 const EvtComplex EvtBBScalar::V_tb = EvtComplex(0.99915, 0);
  39 const EvtComplex EvtBBScalar::V_ts_star = EvtComplex(-0.04029-0.000813*cos(60/180*pi), -0.000813*cos(60/180*pi));
  40 const EvtComplex EvtBBScalar::a4 = EvtComplex(-387.3e-4, -121e-4);
  41 const EvtComplex EvtBBScalar::a6 = EvtComplex(-555.3e-4, -121e-4);
  42 const double EvtBBScalar::x[] = {420.96, -10485.50, 100639.97, -433916.61, 613780.15};
  43 const double EvtBBScalar::y[] = {292.62, -735.73};
  44 const double EvtBBScalar::m_s = 0.120;
  45 const double EvtBBScalar::m_u = 0.029 * 0.120;
  46 const double EvtBBScalar::m_b = 4.88;
  47
  48
  49 EvtBBScalar::EvtBBScalar()
  50     : EvtDecayAmp()
  51     , _massRatio(0)
  52     , _baryonMassSum(0)
  53 {
  54     FormFactor dummy;
  55     dummy.value = 0.36;
  56     dummy.sigma1 = 0.43;
  57     dummy.sigma2 = 0.0;
  58     dummy.mV = 5.42;
  59     _f1Map.insert(make_pair(string("K"), dummy));
  60     dummy.sigma1 = 0.70;
  61     dummy.sigma2 = 0.27;
  62     _f0Map.insert(make_pair(string("K"), dummy));
  63     dummy.value = 0.29;
  64     dummy.sigma1 = 0.48;
  65     dummy.sigma2 = 0.0;
  66     dummy.mV = 5.32;
  67     _f1Map.insert(make_pair(string("pi"), dummy));
  68     dummy.sigma1 = 0.76;
  69     dummy.sigma2 = 0.28;
  70     _f0Map.insert(make_pair(string("pi"), dummy));
  71 }
  72
  73
  74
  75 std::string EvtBBScalar::getName(){
  76     return "B_TO_2BARYON_SCALAR";
  77 }
  78
  79 EvtDecayBase* EvtBBScalar::clone(){
  80     return new EvtBBScalar;
  81 }
  82
  83
  84 void EvtBBScalar::setKnownBaryonTypes(const EvtId& baryon) {
  85     int baryonId = EvtPDL::getStdHep(baryon);
  86     if (EvtPDL::getStdHep(EvtPDL::getId("Lambda0")) == baryonId
  87      or EvtPDL::getStdHep(EvtPDL::getId("anti-Lambda0")) == baryonId ) {
  88         _baryonCombination.set(Lambda);
  89     } else if (EvtPDL::getStdHep(EvtPDL::getId("p+")) == baryonId
  90             or EvtPDL::getStdHep(EvtPDL::getId("anti-p-")) == baryonId ) {
  91         _baryonCombination.set(Proton);
  92     } else if (EvtPDL::getStdHep(EvtPDL::getId("n0")) == baryonId
  93             or EvtPDL::getStdHep(EvtPDL::getId("anti-n0")) == baryonId) {
  94         _baryonCombination.set(Neutron);
  95     } else if (EvtPDL::getStdHep(EvtPDL::getId("Sigma0")) == baryonId
  96             or EvtPDL::getStdHep(EvtPDL::getId("anti-Sigma0")) == baryonId ) {
  97         _baryonCombination.set(Sigma0);
  98     } else if (EvtPDL::getStdHep(EvtPDL::getId("Sigma-")) == baryonId
  99             or EvtPDL::getStdHep(EvtPDL::getId("anti-Sigma+")) == baryonId ) {
 100         _baryonCombination.set(Sigma_minus);
 101     } else if (EvtPDL::getStdHep(EvtPDL::getId("Xi0")) == baryonId
 102             or EvtPDL::getStdHep(EvtPDL::getId("anti-Xi0")) == baryonId) {
 103         _baryonCombination.set(Xi0);
 104     } else if (EvtPDL::getStdHep(EvtPDL::getId("Xi-")) == baryonId
 105             or EvtPDL::getStdHep(EvtPDL::getId("anti-Xi+")) == baryonId) {
 106         _baryonCombination.set(Xi_minus);
 107     } else {
 108         report(ERROR, "EvtGen")
 109             << "EvtBBScalar::init: Don't know what to do with this type as the first or second baryon\n";
 110         exit(2);
 111     }
 112 }
 113
 114 double EvtBBScalar::baryonF1F2(double t) const {
 115     // check for known form factors for combination of baryons
 116     if (_baryonCombination.test(Lambda) and _baryonCombination.test(Proton)) {
 117         return -sqrt(1.5) * G_p(t);
 118     } else if (_baryonCombination.test(Sigma0) and _baryonCombination.test(Proton)) {
 119         return -sqrt(0.5) * (G_p(t) + 2* G_n(t));
 120     } else if (_baryonCombination.test(Sigma_minus) and _baryonCombination.test(Neutron)) {
 121         return -G_p(t) - 2* G_n(t);
 122     } else if (_baryonCombination.test(Xi0) and _baryonCombination.test(Sigma_minus)) {
 123         return G_p(t) - G_n(t);
 124     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Sigma0)) {
 125         return sqrt(0.5) * (G_p(t) - G_n(t));
 126     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Lambda)) {
 127         return sqrt(1.5) * (G_p(t) + G_n(t));
 128     } else {
 129         report(ERROR, "EvtGen")
 130                 << "EvtBBScalar::baryonF1F2: Don't know what to do with this type as the first or second baryon\n";
 131         exit(2);
 132     }
 133 }
 134
 135 double EvtBBScalar::formFactorFit(double t, const vector<double>& params) const {
 136     static const double gamma = 2.148;
 137     static const double Lambda_0 = 0.3;
 138     double result = 0;
 139     for (size_t i=0; i<params.size(); ++i) {
 140         result += params[i]/pow(t, static_cast<int>(i+1));
 141     }
 142     return result * pow(log(t/pow(Lambda_0, 2)), -gamma);
 143 }
 144
 145
 146 double EvtBBScalar::G_p(double t) const {
 147     const vector<double> v_x(x, x+5);
 148     return formFactorFit(t, v_x);
 149 }
 150
 151 double EvtBBScalar::G_n(double t) const {
 152     const vector<double> v_y(y, y+2);
 153     return -formFactorFit(t, v_y);
 154 }
 155
 156 double EvtBBScalar::baryon_gA(double t) const {
 157     // check for known form factors for combination of baryons
 158     if (_baryonCombination.test(Lambda) and _baryonCombination.test(Proton)) {
 159         return -1/sqrt(6.) * (D_A(t) + 3*F_A(t));
 160     } else if (_baryonCombination.test(Sigma0) and _baryonCombination.test(Proton)) {
 161         return 1/sqrt(2.) * (D_A(t) - F_A(t));
 162     } else if (_baryonCombination.test(Sigma_minus) and _baryonCombination.test(Neutron)) {
 163         return D_A(t) - F_A(t);
 164     } else if (_baryonCombination.test(Xi0) and _baryonCombination.test(Sigma_minus)) {
 165         return D_A(t) + F_A(t);
 166     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Sigma0)) {
 167         return 1/sqrt(2.) * (D_A(t) + F_A(t));
 168     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Lambda)) {
 169         return -1 / sqrt(6.) * (D_A(t) - 3*F_A(t));
 170     } else {
 171         report(ERROR, "EvtGen")
 172                 << "EvtBBScalar::baryon_gA: Don't know what to do with this type as the first or second baryon\n";
 173         exit(2);
 174     }
 175 }
 176
 177 double EvtBBScalar::baryon_gP(double t) const {
 178     // check for known form factors for combination of baryons
 179     if (_baryonCombination.test(Lambda) and _baryonCombination.test(Proton)) {
 180         return -1/sqrt(6.) * (D_P(t) + 3*F_P(t));
 181     } else if (_baryonCombination.test(Sigma0) and _baryonCombination.test(Proton)) {
 182         return 1/sqrt(2.) * (D_P(t) - F_P(t));
 183     } else if (_baryonCombination.test(Sigma_minus) and _baryonCombination.test(Neutron)) {
 184         return D_P(t) - F_P(t);
 185     } else if (_baryonCombination.test(Xi0) and _baryonCombination.test(Sigma_minus)) {
 186         return D_P(t) + F_P(t);
 187     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Sigma0)) {
 188         return 1/sqrt(2.) * (D_P(t) + F_P(t));
 189     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Lambda)) {
 190         return -1 / sqrt(6.) * (D_P(t) - 3*F_P(t));
 191     } else {
 192         report(ERROR, "EvtGen")
 193                 << "EvtBBScalar::baryon_gP: Don't know what to do with this type as the first or second baryon\n";
 194         exit(2);
 195     }
 196 }
 197
 198 double EvtBBScalar::baryon_fS(double t) const {
 199     // check for known form factors for combination of baryons
 200     if (_baryonCombination.test(Lambda) and _baryonCombination.test(Proton)) {
 201         return -1/sqrt(6.) * (D_S(t) + 3*F_S(t));
 202     } else if (_baryonCombination.test(Sigma0) and _baryonCombination.test(Proton)) {
 203         return 1/sqrt(2.) * (D_S(t) - F_S(t));
 204     } else if (_baryonCombination.test(Sigma_minus) and _baryonCombination.test(Neutron)) {
 205         return D_S(t) - F_S(t);
 206     } else if (_baryonCombination.test(Xi0) and _baryonCombination.test(Sigma_minus)) {
 207         return D_S(t) + F_S(t);
 208     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Sigma0)) {
 209         return 1/sqrt(2.) * (D_S(t) + F_S(t));
 210     } else if (_baryonCombination.test(Xi_minus) and _baryonCombination.test(Lambda)) {
 211         return -1 / sqrt(6.) * (D_S(t) - 3*F_S(t));
 212     } else {
 213         report(ERROR, "EvtGen")
 214                 << "EvtBBScalar::baryon_fS: Don't know what to do with this type as the first or second baryon\n";
 215         exit(2);
 216     }
 217 }
 218
 219 double EvtBBScalar::D_A(double t) const {
 220     const double d_tilde[] = {x[0]-1.5*y[0], -478};
 221     const vector<double> v_d_tilde(d_tilde, d_tilde+2);
 222     return formFactorFit(t, v_d_tilde);
 223 }
 224
 225 double EvtBBScalar::F_A(double t) const {
 226     const double f_tilde[] = {2./3*x[0]+0.5*y[0], -478};
 227     const vector<double> v_f_tilde(f_tilde, f_tilde+2);
 228     return formFactorFit(t, v_f_tilde);
 229 }
 230
 231 double EvtBBScalar::D_P(double t) const {
 232     const double d_bar[] = {1.5*y[0]* _massRatio, /*-952*/0};
 233     const vector<double> v_d_bar(d_bar, d_bar+2);
 234     return formFactorFit(t, v_d_bar);
 235 }
 236
 237 double EvtBBScalar::F_P(double t) const {
 238     const double f_bar[] = {(x[0]-0.5*y[0]) * _massRatio, /*-952*/0};
 239     const vector<double> v_f_bar(f_bar, f_bar+2);
 240     return formFactorFit(t, v_f_bar);
 241 }
 242
 243 double EvtBBScalar::D_S(double t) const {
 244     return -1.5 * _massRatio * G_n(t);
 245 }
 246
 247 double EvtBBScalar::F_S(double t) const {
 248     return (G_p(t) + 0.5*G_n(t)) * _massRatio;
 249 }
 250
 251 double EvtBBScalar::baryon_hA(double t) const {
 252     return (1/_massRatio*baryon_gP(t)-baryon_gA(t))*pow(_baryonMassSum, 2)/t;
 253 }
 254
 255
 256 void EvtBBScalar::init() {
 257     // no arguments, daughter lambda p_bar pi
 258     // charge conservation is checked by base class
 259     checkNArg(0);
 260     checkNDaug(3);
 261     checkSpinParent(EvtSpinType::SCALAR);
 262     checkSpinDaughter(0, EvtSpinType::DIRAC);
 263     checkSpinDaughter(1, EvtSpinType::DIRAC);
 264     checkSpinDaughter(2, EvtSpinType::SCALAR);
 265     EvtId baryon1 = getDaug(0);
 266     EvtId baryon2 = getDaug(1);
 267     EvtId scalar = getDaug(2);
 268     int scalarId = EvtPDL::getStdHep(scalar);
 269
 270     // Different form factors for the B-pi or B-K transition.
 271     if (   scalarId == EvtPDL::getStdHep(EvtPDL::getId("pi+"))
 272         or scalarId == EvtPDL::getStdHep(EvtPDL::getId("pi-"))
 273         or scalarId == EvtPDL::getStdHep(EvtPDL::getId("pi0"))) {
 274         _scalarType = "pi";
 275     } else if (scalarId == EvtPDL::getStdHep(EvtPDL::getId("K+"))
 276         or scalarId == EvtPDL::getStdHep(EvtPDL::getId("K-"))
 277         or scalarId == EvtPDL::getStdHep(EvtPDL::getId("K0"))
 278         or scalarId == EvtPDL::getStdHep(EvtPDL::getId("anti-K0"))) {
 279         _scalarType = "K";
 280     } else {
 281         report(ERROR, "EvtGen")
 282             << "EvtBBScalar::init: Can only deal with Kaons or pions as the third particle\n"
 283                 << "\tFound: " << scalarId << endl;
 284         exit(2);
 285     }
 286     // check for known particles
 287     setKnownBaryonTypes(baryon1);
 288     setKnownBaryonTypes(baryon2);
 289     double mass1 = EvtPDL::getMass(baryon1);
 290     double mass2 = EvtPDL::getMass(baryon2);
 291     // This whole model deals only with baryons that differ in s-u
 292     if (mass1 > mass2)
 293         _massRatio = (mass1-mass2) / (m_s-m_u);
 294     else
 295         _massRatio = (mass2-mass1) / (m_s-m_u);
 296     _baryonMassSum = mass1 + mass2;
 297 }
 298
 299
 300 // initialize phasespace and calculate the amplitude
 301 void EvtBBScalar::decay(EvtParticle* p) {
 302     p->initializePhaseSpace(getNDaug(), getDaugs());
 303     EvtVector4R B_Momentum = p->getP4Lab();
 304     EvtDiracParticle* theLambda = dynamic_cast<EvtDiracParticle*>(p->getDaug(0));
 305     EvtDiracParticle* theAntiP = dynamic_cast<EvtDiracParticle*>(p->getDaug(1));
 306     EvtScalarParticle* theScalar = dynamic_cast<EvtScalarParticle*>(p->getDaug(2));
 307     EvtVector4R scalarMomentum = theScalar->getP4Lab();
 308
 309     // The amplitude consists of three matrix elements. These will be calculated one by one here.
 310
 311     // loop over all possible spin states
 312     for (int i=0; i<2; ++i) {
 313     EvtDiracSpinor lambdaPol = theLambda->spParent(i);
 314         for (int j=0; j<2; ++j)  {
 315             EvtDiracSpinor antiP_Pol = theAntiP->spParent(j);
 316             EvtVector4C theAmplitudePartA = amp_A(B_Momentum, scalarMomentum);
 317             EvtComplex amplitude;
 318             for (int index=0; index<4; ++index) {
 319                 amplitude += theAmplitudePartA.get(index)
 320                         * ( const_B*amp_B(theLambda, lambdaPol, theAntiP, antiP_Pol, index)
 321                           + const_C*amp_C(theLambda, lambdaPol, theAntiP, antiP_Pol, index) );
 322             }
 323             vertex(i, j, amplitude);
 324         }
 325     }
 326 }
 327
 328 void EvtBBScalar::initProbMax()
 329 {
 330     // setProbMax(1);
 331     setProbMax(0.2); // found by trial and error
 332 }
 333
 334 // Form factor f1 for B-pi transition
 335 double EvtBBScalar::B_pi_f1(double t) const
 336 {
 337     FormFactor f = _f1Map[_scalarType];
 338     double mv2 = f.mV*f.mV;
 339     return f.value / ((1-t/mv2) * (1-f.sigma1*t/mv2+f.sigma2*t*t/mv2/mv2));
 340 }
 341
 342 // Form factor f0 for B-pi transition
 343 double EvtBBScalar::B_pi_f0(double t) const
 344 {
 345     FormFactor f = _f0Map[_scalarType];
 346     double mv2 = f.mV*f.mV;
 347     return f.value / (1 - f.sigma1*t/mv2 + f.sigma2*t*t/mv2/mv2);
 348 }
 349
 350 // constants of the B and C parts of the amplitude
 351 const EvtComplex EvtBBScalar::const_B = V_ub*V_us_star*a1 - V_tb*V_ts_star*a4;
 352 const EvtComplex EvtBBScalar::const_C = 2*a6*V_tb*V_ts_star;
 353
 354 // part A of the amplitude, see hep-ph/0204185
 355 const EvtVector4C
 356 EvtBBScalar::amp_A(const EvtVector4R& p4B, const EvtVector4R& p4Scalar)
 357 {
 358     double mB2 = p4B.mass2();
 359     double mScalar2 = p4Scalar.mass2();
 360     double t = (p4B-p4Scalar).mass2();
 361     return ((p4B+p4Scalar) - (mB2-mScalar2)/t * (p4B-p4Scalar)) * B_pi_f1(t)
 362          + (mB2-mScalar2)/t * (p4B-p4Scalar) * B_pi_f0(t);
 363 }
 364
 365 // part B of the amplitude, Vector and Axial Vector parts
 366 const EvtComplex
 367 EvtBBScalar::amp_B(const EvtDiracParticle* baryon1, const EvtDiracSpinor& b1Pol
 368                  , const EvtDiracParticle* baryon2, const EvtDiracSpinor& b2Pol
 369                  , int index)
 370 {
 371     return amp_B_vectorPart(baryon1, b1Pol, baryon2, b2Pol, index)
 372          - amp_B_axialPart(baryon1, b1Pol, baryon2, b2Pol, index);
 373 }
 374
 375
 376 const EvtComplex
 377 EvtBBScalar::amp_B_vectorPart(const EvtDiracParticle* baryon1, const EvtDiracSpinor& b1Pol
 378                             , const EvtDiracParticle* baryon2, const EvtDiracSpinor& b2Pol
 379                             , int index)
 380 {
 381     double t = (baryon1->getP4Lab() + baryon2->getP4Lab()).mass2();
 382     EvtGammaMatrix gamma;
 383     for (int i=0; i<4; ++i) {
 384         gamma += EvtTensor4C::g().get(index, i) * EvtGammaMatrix::g(i);
 385     }
 386     // The F2 contribution that is written out in the paper is neglected here.
 387     // see hep-ph/0204185
 388     EvtDiracSpinor A = EvtComplex(baryonF1F2(t))*b2Pol ;
 389     EvtDiracSpinor Adjb1Pol = b1Pol.adjoint() ;
 390     EvtDiracSpinor gammaA = gamma * A ;
 391     return Adjb1Pol * gammaA ;
 392     //    return b1Pol.adjoint()*(gamma*(EvtComplex(baryonF1F2(t))*b2Pol));
 393 }
 394
 395 const EvtComplex
 396 EvtBBScalar::amp_B_axialPart(const EvtDiracParticle* baryon1, const EvtDiracSpinor& b1Pol
 397                            , const EvtDiracParticle* baryon2, const EvtDiracSpinor& b2Pol
 398                            , int index)
 399 {
 400     EvtGammaMatrix gamma;
 401     for (int i=0; i<4; ++i) {
 402         gamma += EvtTensor4C::g().get(index, i) * EvtGammaMatrix::g(i);
 403     }
 404     double t = (baryon1->getP4Lab() + baryon2->getP4Lab()).mass2();
 405     double mSum = baryon1->mass() + baryon2->mass();
 406     EvtVector4C momentum_upper = (baryon1->getP4Lab()+baryon2->getP4Lab());
 407     EvtVector4C momentum;
 408     for (int mu=0; mu<0; ++mu) {
 409         EvtComplex dummy;
 410         for (int i=0; i<4; ++i) {
 411             dummy += EvtTensor4C::g().get(index, i)*momentum_upper.get(i);
 412         }
 413         momentum.set(mu, dummy);
 414     }
 415     return b1Pol.adjoint() * (((baryon_gA(t) * gamma +
 416                               EvtGammaMatrix::id()*baryon_hA(t)/mSum*momentum.get(index))
 417             * EvtGammaMatrix::g5()) * b2Pol);
 418 }
 419
 420
 421 // part C of the amplitude, Scalar and Pseudoscalar parts
 422 const EvtComplex
 423 EvtBBScalar::amp_C(const EvtDiracParticle* baryon1, const EvtDiracSpinor& b1Pol
 424                  , const EvtDiracParticle* baryon2, const EvtDiracSpinor& b2Pol
 425                  , int index)
 426 {
 427     EvtVector4C baryonSumP4_upper = baryon1->getP4Lab() + baryon2->getP4Lab();
 428     EvtVector4C baryonSumP4;
 429     for (int mu=0; mu<4; ++mu) {
 430         EvtComplex dummy;
 431         for (int i=0; i<4; ++i) {
 432             dummy += EvtTensor4C::g().get(mu, i) * baryonSumP4_upper.get(i);
 433         }
 434         baryonSumP4.set(mu, dummy);
 435     }
 436     double t = (baryon1->getP4Lab() + baryon2->getP4Lab()).mass2();
 437     return baryonSumP4.get(index)/(m_b-m_u)*(amp_C_scalarPart(b1Pol, b2Pol, t) + amp_C_pseudoscalarPart(b1Pol, b2Pol, t));
 438 }
 439
 440
 441 const EvtComplex
 442 EvtBBScalar::amp_C_scalarPart(const EvtDiracSpinor& b1Pol, const EvtDiracSpinor& b2Pol, double t)
 443 {
 444     return baryon_fS(t) * b1Pol.adjoint()*b2Pol;
 445 }
 446
 447 const EvtComplex
 448 EvtBBScalar::amp_C_pseudoscalarPart(const EvtDiracSpinor& b1Pol, const EvtDiracSpinor& b2Pol, double t)
 449 {
 450     return baryon_gP(t) * b1Pol.adjoint()*(EvtGammaMatrix::g5()*b2Pol);
 451 }
 452