TEvtGen/EvtGenModels/EvtBcToNPi.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package. If you use all or part
   5 //      of it, please give an appropriate acknowledgement.
   6 //
   7 // Copyright Information: See EvtGen/COPYRIGHT
   8 //
   9 // Module: EvtGenModels/EvtBcToNPi.hh
  10 //
  11 // Description: General decay model for Bc -> V + npi and Bc -> P + npi
  12 //
  13 // Modification history:
  14 //
  15 //    A.Berezhnoy, A.Likhoded, A.Luchinsky  April 2011   Module created
  16 //
  17 //------------------------------------------------------------------------
  18
  19 #include "EvtGenBase/EvtPatches.hh"
  20
  21 #include "EvtGenBase/EvtPDL.hh"
  22 #include "EvtGenBase/EvtReport.hh"
  23 #include "EvtGenBase/EvtVector4C.hh"
  24 #include "EvtGenBase/EvtTensor4C.hh"
  25 #include "EvtGenBase/EvtIdSet.hh"
  26 #include "EvtGenBase/EvtSpinType.hh"
  27 #include "EvtGenBase/EvtScalarParticle.hh"
  28
  29 #include "EvtGenModels/EvtBcToNPi.hh"
  30
  31 #include <iostream>
  32 using std::endl;
  33
  34 EvtBcToNPi::EvtBcToNPi(bool printAuthorInfo) {
  35   nCall=0; maxAmp2=0;
  36   if (printAuthorInfo == true) {this->printAuthorInfo();}
  37 }
  38
  39 EvtBcToNPi::~EvtBcToNPi() {
  40 }
  41
  42 std::string EvtBcToNPi::getName(){
  43
  44   return "EvtBcToNPi";
  45
  46 }
  47
  48 EvtDecayBase* EvtBcToNPi::clone(){
  49
  50   return new EvtBcToNPi;
  51
  52 }
  53
  54
  55 void EvtBcToNPi::init(){
  56
  57         // check spins
  58         checkSpinParent(EvtSpinType::SCALAR);
  59
  60
  61         // the others are scalar
  62         for (int i=1; i<=(getNDaug()-1);i++) {
  63                 checkSpinDaughter(i,EvtSpinType::SCALAR);
  64         };
  65
  66         _beta=-0.108; _mRho=0.775; _gammaRho=0.149;
  67         _mRhopr=1.364; _gammaRhopr=0.400; _mA1=1.23; _gammaA1=0.4;
  68
  69         // read arguments
  70         if( EvtPDL::getSpinType(getDaug(0)) == EvtSpinType::VECTOR) {
  71                 checkNArg(10);
  72                 int n=0;
  73                 _maxProb=getArg(n++);
  74                 FA0_N=getArg(n++);
  75                 FA0_c1=getArg(n++);
  76                 FA0_c2=getArg(n++);
  77                 FAp_N=getArg(n++);
  78                 FAp_c1=getArg(n++);
  79                 FAp_c2=getArg(n++);
  80                 FV_N=getArg(n++);
  81                 FV_c1=getArg(n++);
  82                 FV_c2=getArg(n++);
  83                 FAm_N=0;
  84                 FAm_c1=0;
  85                 FAm_c2=0;
  86         }
  87         else if( EvtPDL::getSpinType(getDaug(0)) == EvtSpinType::SCALAR) {
  88                 checkNArg(4);
  89                 int n=0;
  90                 _maxProb=getArg(n++);
  91                 Fp_N=getArg(n++);
  92                 Fp_c1=getArg(n++);
  93                 Fp_c2=getArg(n++);
  94                 Fm_N=0;
  95                 Fm_c1=0;
  96                 Fm_c2=0;
  97         }
  98         else {
  99                 report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCPSINPI model" << endl;
 100                 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
 101                 for ( int id=0; id<(getNDaug()-1); id++ )
 102                         report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
 103                 return;
 104
 105         };
 106
 107         if(getNDaug()<2 || getNDaug()>4) {
 108                 report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCPSINPI model" << endl;
 109                 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
 110                 for ( int id=0; id<(getNDaug()-1); id++ )
 111                         report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
 112                 return;
 113         }
 114
 115 }
 116
 117 double EvtBcToNPi::_ee(double M, double m1, double m2) {
 118         return (M*M+m1*m1-m2*m2)/(2*M);
 119 };
 120
 121 double EvtBcToNPi::_pp(double M, double m1, double m2) {
 122         double __ee=_ee(M,m1,m2);
 123         return sqrt(__ee*__ee-m1*m1);
 124 };
 125
 126
 127 void EvtBcToNPi::initProbMax(){
 128         if(_maxProb>0.) setProbMax(_maxProb);
 129         else {
 130                 EvtId id=getParentId();
 131                 EvtScalarParticle *p=new EvtScalarParticle();
 132                 p->init(id, EvtPDL::getMass(id),0., 0., 0.);
 133                 p->setDiagonalSpinDensity();
 134                 // add daughters
 135                 p->makeDaughters(getNDaug(), getDaugs() );
 136
 137                 // fill the momenta
 138                 if(getNDaug()==2) {
 139                         double M=EvtPDL::getMass(id), m1=EvtPDL::getMass(getDaug(0)), m2=EvtPDL::getMass(getDaug(1));
 140                         double __pp=_pp(M,m1,m2);
 141                         p->getDaug(0)->setP4( EvtVector4R( _ee(M,m1,m2), 0., 0., __pp) );
 142                         p->getDaug(1)->setP4( EvtVector4R( _ee(M,m2,m1), 0., 0., -__pp) );
 143                 }
 144                 else if( getNDaug()==3) {
 145                         double M=EvtPDL::getMass(id),
 146                                 m1=EvtPDL::getMass(getDaug(0)),
 147                                 m2=EvtPDL::getMass(getDaug(1)),
 148                                 m3=EvtPDL::getMass(getDaug(2));
 149                         double __ppRho=_pp(M,m1,_mRho), __ppPi=_pp(_mRho,m2,m3);
 150                         p->getDaug(0)->setP4( EvtVector4R( _ee(M,m1,_mRho), 0., 0., __ppRho) );
 151                         EvtVector4R _pRho( _ee(M,_mRho,m1), 0., 0., -__ppRho);
 152                         EvtVector4R _p2( _ee(_mRho, m2, m3), 0., 0., __ppPi); _p2.applyBoostTo(_pRho);
 153                         EvtVector4R _p3( _ee(_mRho, m2, m3), 0., 0., -__ppPi); _p3.applyBoostTo(_pRho);
 154                         p->getDaug(1)->setP4(_p2);
 155                         p->getDaug(2)->setP4(_p3);
 156
 157                 }
 158                 else if( getNDaug()==4) {
 159                         double M=EvtPDL::getMass(id),
 160                                 m1=EvtPDL::getMass(getDaug(0)),
 161                                 m2=EvtPDL::getMass(getDaug(1)),
 162                                 m3=EvtPDL::getMass(getDaug(2)),
 163                                 m4=EvtPDL::getMass(getDaug(3));
 164                         if(M<m1+_mA1) return;
 165                         double   __ppA1=_pp(M,m1,_mA1),
 166                                  __ppRho=_pp(_mA1,_mRho,m4),
 167                                  __ppPi=_pp(_mRho, m2, m3);
 168                         p->getDaug(0)->setP4( EvtVector4R( _ee(M,m1,_mRho), 0., 0., __ppA1) );
 169                         EvtVector4R _pA1( _ee(M,_mA1,m1), 0., 0., -__ppA1);
 170                         EvtVector4R _pRho(_ee(_mA1, _mRho, m4), 0, 0, __ppRho);
 171                         _pRho.applyBoostTo(_pA1);
 172                         EvtVector4R _p4( _ee(_mA1, m4, _mRho), 0, 0, -__ppRho); _p4.applyBoostTo(_pA1);
 173                         p->getDaug(3)->setP4(_p4);
 174                         EvtVector4R _p2( _ee(_mRho, m2, m3), 0, 0, __ppPi); _p2.applyBoostTo(_pRho);
 175                         p->getDaug(1)->setP4(_p2);
 176                         EvtVector4R _p3( _ee(_mRho, m2, m3), 0, 0, -__ppPi); _p2.applyBoostTo(_pRho);
 177                         p->getDaug(2)->setP4(_p3);
 178                 };
 179
 180
 181                 _amp2.init(p->getId(),getNDaug(),getDaugs());
 182
 183                 decay(p);
 184
 185                 EvtSpinDensity rho=_amp2.getSpinDensity();
 186
 187                 double prob=p->getSpinDensityForward().normalizedProb(rho);
 188
 189                 if(prob>0) setProbMax(0.9*prob);
 190
 191
 192         };
 193 }
 194
 195 void EvtBcToNPi::decay( EvtParticle *root_particle ){
 196         ++nCall;
 197
 198         EvtIdSet thePis("pi+","pi-","pi0");
 199         EvtComplex I=EvtComplex(0.0, 1.0);
 200
 201
 202         root_particle->initializePhaseSpace(getNDaug(),getDaugs());
 203
 204         EvtVector4R
 205                 p(root_particle->mass(), 0., 0., 0.),                  // Bc momentum
 206                 k=root_particle->getDaug(0)->getP4(),                      // J/psi momenta
 207                 Q=p-k;
 208
 209         double Q2=Q.mass2();
 210
 211
 212 // check pi-mesons and calculate hadronic current
 213         EvtVector4C hardCur;
 214         bool foundHadCurr=false;
 215
 216         if ( getNDaug() == 2 ) // Bc -> psi pi+
 217         {
 218                 hardCur=Q;
 219                 foundHadCurr=true;
 220         }
 221         else if ( getNDaug() == 3 ) // Bc -> psi pi+ pi0
 222         {
 223                 EvtVector4R p1,p2;
 224                 p1=root_particle->getDaug(1)->getP4(),     // pi+ momenta
 225                 p2=root_particle->getDaug(2)->getP4(),    // pi0 momentum
 226                 hardCur=Fpi(p1,p2)*(p1-p2);
 227                 foundHadCurr=true;
 228         }
 229         else if( getNDaug()==4 ) // Bc -> psi pi+ pi pi
 230         {
 231                 int diffPi(0),samePi1(0),samePi2(0);
 232                 if ( getDaug( 1) == getDaug( 2) ) {diffPi= 3; samePi1= 1; samePi2= 2;}
 233                 if ( getDaug( 1) == getDaug( 3) ) {diffPi= 2; samePi1= 1; samePi2= 3;}
 234                 if ( getDaug( 2) == getDaug( 3) ) {diffPi= 1; samePi1= 2; samePi2= 3;}
 235
 236                 EvtVector4R p1=root_particle->getDaug(samePi1)->getP4();
 237                 EvtVector4R p2=root_particle->getDaug(samePi2)->getP4();
 238                 EvtVector4R p3=root_particle->getDaug(diffPi)->getP4();
 239
 240                 EvtComplex BA1;
 241                 double GA1=_gammaA1*pi3G(Q2,samePi1)/pi3G(_mA1*_mA1,samePi1);
 242                 EvtComplex denBA1(_mA1*_mA1 - Q.mass2(),-1.*_mA1*GA1);
 243                 BA1 = _mA1*_mA1 / denBA1;
 244
 245                 hardCur = BA1*( (p1-p3) - (Q*(Q*(p1-p3))/Q2)*Fpi(p2,p3) +
 246                                 (p2-p3) - (Q*(Q*(p2-p3))/Q2)*Fpi(p1,p3) );
 247                 foundHadCurr=true;
 248         }
 249
 250         if ( !foundHadCurr ) {
 251                 report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
 252                 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
 253                 int id;
 254                 for ( id=0; id<(getNDaug()-1); id++ )
 255                 report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
 256                 ::abort();
 257         };
 258
 259         EvtTensor4C H;
 260         double amp2=0.;
 261         if( root_particle->getDaug(0)->getSpinType() == EvtSpinType::VECTOR) {
 262                 double FA0=FA0_N*exp(FA0_c1*Q2 + FA0_c2*Q2*Q2);
 263                 double FAp=FAp_N*exp(FAp_c1*Q2 + FAp_c2*Q2*Q2);
 264                 double FAm=FAm_N*exp(FAm_c1*Q2 + FAm_c2*Q2*Q2);
 265                 double FV= FV_N* exp( FV_c1*Q2 + FV_c2*Q2*Q2 );
 266                 H=-FA0*EvtTensor4C::g()
 267                         -FAp*EvtGenFunctions::directProd(p,p+k)
 268                         +FAm*EvtGenFunctions::directProd(p,p-k)
 269                         +2*I*FV*dual(EvtGenFunctions::directProd(p,k));
 270                 EvtVector4C  Heps=H.cont2(hardCur);
 271
 272                 for(int i=0; i<4; i++) {
 273                         EvtVector4C  eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
 274                         EvtComplex amp=eps*Heps;
 275                         vertex(i,amp);
 276                         amp2+=pow( abs(amp),2);
 277                 }
 278         }
 279         else if( root_particle->getDaug(0)->getSpinType() == EvtSpinType::SCALAR) {
 280                 double Fp=Fp_N*exp(Fp_c1*Q2 + Fp_c2*Q2*Q2);
 281                 double Fm=Fm_N*exp(Fm_c1*Q2 + Fm_c2*Q2*Q2);
 282                 EvtVector4C H=Fp*(p+k)+Fm*(p-k);
 283                 EvtComplex amp=H*hardCur;
 284                 vertex(amp);
 285                 amp2+=pow( abs(amp),2);
 286         };
 287         if(amp2>maxAmp2) maxAmp2=amp2;
 288
 289   return ;
 290 }
 291
 292 EvtComplex EvtBcToNPi::Fpi( EvtVector4R q1, EvtVector4R q2) {
 293         double m1=q1.mass();
 294         double m2=q2.mass();
 295
 296         EvtVector4R Q = q1 + q2;
 297         double mQ2= Q*Q;
 298
 299         // momenta in the rho->pipi decay
 300         double dRho= _mRho*_mRho - m1*m1 - m2*m2;
 301         double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
 302
 303         double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
 304         double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
 305
 306         double dQ= mQ2 - m1*m1 - m2*m2;
 307         double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
 308
 309
 310         double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
 311         EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
 312         EvtComplex BRho= _mRho*_mRho / BRhoDem;
 313
 314         double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
 315         EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
 316         EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
 317
 318         return (BRho + _beta*BRhopr)/(1+_beta);
 319 }
 320
 321 double EvtBcToNPi::pi3G(double m2,int dupD) {
 322         double mPi= EvtPDL::getMeanMass(getDaug(dupD));
 323         if ( m2 > (_mRho+mPi) ) {
 324         return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
 325         }
 326         else {
 327         double t1=m2-9.0*mPi*mPi;
 328         return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
 329         }
 330 }
 331
 332 void EvtBcToNPi::printAuthorInfo() {
 333
 334   report(INFO,"EvtGen")<<"Defining EvtBcToNPi model: Bc -> V + npi and Bc -> P + npi decays\n"
 335                        <<"from A.V. Berezhnoy, A.K. Likhoded, A.V. Luchinsky: "
 336                        <<"Phys.Rev.D 82, 014012 (2010) and arXiV:1104.0808."<<endl;
 337
 338 }
 339