TEvtGen/EvtGenModels/EvtWnPi.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) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtTVP.cc
  12 //
  13 // Description: Routine to calculate W -> (n pi) current
  14 //                      according to [Kuhn, Was, Acta.Phys.Polon B39 (2008) 147]
  15 //
  16 // Modification history:
  17 //      AVL     6 July, 2012    Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <iomanip>
  24 #include <fstream>
  25 #include <ctype.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28 #include "EvtGenBase/EvtParticle.hh"
  29 #include "EvtGenBase/EvtPDL.hh"
  30 #include "EvtGenBase/EvtGenKine.hh"
  31 #include "EvtGenModels/EvtTauHadnu.hh"
  32 #include "EvtGenBase/EvtDiracSpinor.hh"
  33 #include "EvtGenBase/EvtReport.hh"
  34 #include "EvtGenBase/EvtVector4C.hh"
  35 #include "EvtGenBase/EvtTensor4C.hh"
  36 #include "EvtGenBase/EvtIdSet.hh"
  37 #include "EvtGenBase/EvtParser.hh"
  38
  39 #include "EvtGenModels/EvtWnPi.hh"
  40
  41 using namespace std;
  42
  43 // W+ -> pi_ current
  44 EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1) {
  45   return q1;
  46 }
  47
  48 // W+ -> pi+ pi0 current
  49 EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2) {
  50   return BWr(q1+q2)*(q1-q2);
  51 }
  52
  53 // W+ -> pi+ pi+ pi- current
  54 EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3) {
  55  EvtVector4R Q=q1+q2+q3;
  56  double Q2=Q.mass2();
  57  return BWa(Q)*( (q1-q3) - (Q*(Q*(q1-q3))/Q2)*BWr(q2+q3) +
  58               (q2-q3) - (Q*(Q*(q2-q3))/Q2)*BWr(q1+q3) );
  59 }
  60
  61 // W+ -> pi+ pi+ pi- pi- pi+ current with symmetrization
  62 EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3, EvtVector4R q4, EvtVector4R q5) {
  63 //  double Q2 = Qtot*Qtot;
  64 //  return q1-Qtot*(q1*Qtot)/Q2;
  65   EvtVector4C V = JB(q1, q2, q3, q4, q5) + JB(q5, q2, q3, q4, q1) + JB(q1, q5, q3, q4, q2) +
  66                   JB(q1,q2,q4,q3,q5)+JB(q5,q2,q4,q3,q1)+JB(q1,q5,q4,q3,q2);
  67 //  cout<<"BC2: Qtot="<<Qtot<<", V="<<V<<endl;
  68   return V;
  69 }
  70
  71
  72 // a1 -> pi+ pi+ pi- BW
  73 EvtComplex EvtWnPi::BWa(EvtVector4R q) {
  74   double const _mA1=1.26, _GA1=0.4;
  75   EvtComplex I(0,1);
  76   double Q2 = q.mass2();
  77   double GA1=_GA1*pi3G(Q2)/pi3G(_mA1*_mA1);
  78   EvtComplex denBA1(_mA1*_mA1 - Q2,-1.*_mA1*GA1);
  79   return _mA1*_mA1 / denBA1;
  80 }
  81
  82
  83 EvtComplex EvtWnPi::BWf(EvtVector4R q) {
  84   double const mf=0.8, Gf=0.6;
  85   EvtComplex I(0,1);
  86   double Q2 = q.mass2();
  87   return mf*mf/(mf*mf-Q2-I*mf*Gf);
  88 }
  89
  90 EvtComplex EvtWnPi::BWr(EvtVector4R q) {
  91         double _mRho = 0.775, _gammaRho=0.149, _mRhopr=1.364,  _gammaRhopr=0.400, _beta=-0.108;
  92         double m1=EvtPDL::getMeanMass(EvtPDL::getId("pi+")), m2=EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
  93         double mQ2=q.mass2();
  94
  95         // momenta in the rho->pipi decay
  96         double dRho= _mRho*_mRho - m1*m1 - m2*m2;
  97         double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
  98
  99         double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
 100         double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
 101
 102         double dQ= mQ2 - m1*m1 - m2*m2;
 103         double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
 104
 105
 106         double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
 107         EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
 108         EvtComplex BRho= _mRho*_mRho / BRhoDem;
 109
 110         double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
 111         EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
 112         EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
 113
 114         return (BRho + _beta*BRhopr)/(1+_beta);
 115 }
 116
 117 double EvtWnPi::pi3G(double m2) {
 118  double mPi = EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
 119  double _mRho = 0.775;
 120   if ( m2 > (_mRho+mPi) ) {
 121     return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
 122   }
 123   else {
 124     double t1=m2-9.0*mPi*mPi;
 125     return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
 126   };
 127 }
 128
 129 EvtVector4C EvtWnPi::JB( EvtVector4R p1, EvtVector4R p2, EvtVector4R p3, EvtVector4R p4, EvtVector4R p5) {
 130   EvtVector4R Qtot = p1+p2+p3+p4+p5, Qa=p1+p2+p3;
 131   EvtTensor4C T= (1/Qtot.mass2())*EvtGenFunctions::directProd(Qtot,Qtot) - EvtTensor4C::g();
 132   EvtVector4R V13 = Qa*( p2*(p1-p3) )/Qa.mass2() - (p1-p3);
 133   EvtVector4R V23 = Qa*( p1*(p2-p3) )/Qa.mass2() - (p2-p3);
 134   return BWa(Qtot)*BWa(Qa)*BWf(p4+p5)*(
 135     T.cont1(V13)*BWr(p1+p3) + T.cont1(V23)*BWr(p2+p3)
 136   );
 137 }