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