TEvtGen/EvtGen/EvtGenModels/EvtBcVNpi.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: EvtBcVNpi.cc
  12 //
  13 // Description: Module to implement Bc -> psi + (n pi) decays
  14 //
  15 // Modification history:
  16 //
  17 //    AVL     July 6, 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/EvtBcVNpi.hh"
  40 #include "EvtGenModels/EvtWnPi.hh"
  41
  42
  43
  44 EvtBcVNpi::~EvtBcVNpi() {
  45 //   cout<<"BcVNpi::destructor : nCall = "<<nCall<<" getProbMax(-1) = "<<getProbMax(-1)<<endl;
  46
  47 }
  48
  49 std::string EvtBcVNpi::getName(){ return "BC_VNPI";}
  50
  51 EvtDecayBase* EvtBcVNpi::clone() {  return new EvtBcVNpi;}
  52
  53 //======================================================
  54 void EvtBcVNpi::init(){
  55     //cout<<"BcVNpi::init()"<<endl;
  56
  57     checkNArg(1);
  58     checkSpinParent(EvtSpinType::SCALAR);
  59     checkSpinDaughter(0,EvtSpinType::VECTOR);
  60     for (int i=1; i<=(getNDaug()-1);i++) {
  61       checkSpinDaughter(i,EvtSpinType::SCALAR);
  62     };
  63
  64     if(getNDaug()<2 || getNDaug()>6) {
  65       report(ERROR,"EvtGen") << "Have not yet implemented this final state in BcVNpi model" << endl;
  66       report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
  67       for ( int id=0; id<(getNDaug()-1); id++ )
  68         report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
  69       return;
  70     }
  71
  72
  73 //     for(int i=0; i<getNDaug(); i++)
  74 //       cout<<"BcVNpi::init \t\t daughter "<<i<<" : "<<getDaug(i).getId()<<"   "<<EvtPDL::name(getDaug(i)).c_str()<<endl;
  75
  76    idVector = getDaug(0).getId();
  77     whichfit = int(getArg(0)+0.1);
  78 //     cout<<"BcVNpi: whichfit ="<<whichfit<<"  idVector="<<idVector<<endl;
  79     ffmodel = new EvtBCVFF(idVector,whichfit);
  80
  81     wcurr = new EvtWnPi();
  82
  83     nCall = 0;
  84 }
  85
  86 //======================================================
  87 void EvtBcVNpi::initProbMax() {
  88 //     cout<<"BcVNpi::initProbMax()"<<endl;
  89     if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==6) setProbMax(720000.);
  90     else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==6) setProbMax(471817.);
  91     else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==4) setProbMax(42000.);
  92     else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==4) setProbMax(16000.);
  93
  94     else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==4) setProbMax(1200.);
  95     else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==4) setProbMax(2600.);
  96     else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==6) setProbMax(40000.);
  97     else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==6) setProbMax(30000.);
  98 }
  99
 100 //======================================================
 101 void EvtBcVNpi::decay( EvtParticle *root_particle ) {
 102    ++nCall;
 103 //     cout<<"BcVNpi::decay()"<<endl;
 104     root_particle->initializePhaseSpace(getNDaug(),getDaugs());
 105
 106     EvtVector4R
 107             p4b(root_particle->mass(), 0., 0., 0.),                  // Bc momentum
 108             p4meson=root_particle->getDaug(0)->getP4(),                    // J/psi momenta
 109             Q=p4b-p4meson;
 110     double Q2=Q.mass2();
 111
 112
 113 // check pi-mesons and calculate hadronic current
 114     EvtVector4C hardCur;
 115 //     bool foundHadCurr=false;
 116     if( getNDaug() == 2) {
 117       hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() );
 118 //       foundHadCurr=true;
 119     }
 120     else if( getNDaug() == 3) {
 121       hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
 122                                  root_particle->getDaug(2)->getP4()
 123                                );
 124 //       foundHadCurr=true;
 125     }
 126     else if( getNDaug() == 4) {
 127       hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
 128                                  root_particle->getDaug(2)->getP4(),
 129                                  root_particle->getDaug(3)->getP4()
 130                                );
 131 //       foundHadCurr=true;
 132     }
 133     else if( getNDaug() == 6) // Bc -> psi pi+ pi+ pi- pi- pi+ from [Kuhn, Was, hep-ph/0602162
 134     {
 135
 136                 hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(),
 137                                           root_particle->getDaug(2)->getP4(),
 138                                           root_particle->getDaug(3)->getP4(),
 139                                           root_particle->getDaug(4)->getP4(),
 140                                           root_particle->getDaug(5)->getP4()
 141                                  );
 142 //              foundHadCurr=true;
 143     }
 144     else {
 145             report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
 146             report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
 147             int id;
 148             for ( id=0; id<(getNDaug()-1); id++ )
 149             report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
 150             ::abort();
 151     };
 152
 153 // calculate Bc -> V W form-factors
 154         double a1f, a2f, vf, a0f;
 155         double m_meson = root_particle->getDaug(0)->mass();
 156         double m_b = root_particle->mass();
 157         ffmodel->getvectorff(root_particle->getId(),
 158                                 root_particle->getDaug(0)->getId(),
 159                                 Q2,
 160                                 m_meson,
 161                                 &a1f,
 162                                 &a2f,
 163                                 &vf,
 164                                 &a0f);
 165         double a3f = ((m_b+m_meson)/(2.0*m_meson))*a1f -
 166               ((m_b-m_meson)/(2.0*m_meson))*a2f;
 167
 168 // calculate Bc -> V W current
 169         EvtTensor4C H;
 170         H = a1f*(m_b+m_meson)*EvtTensor4C::g();
 171         H.addDirProd((-a2f/(m_b+m_meson))*p4b,p4b+p4meson);
 172         H+=EvtComplex(0.0,vf/(m_b+m_meson))*dual(EvtGenFunctions::directProd(p4meson+p4b,p4b-p4meson));
 173         H.addDirProd((a0f-a3f)*2.0*(m_meson/Q2)*p4b,p4b-p4meson);
 174         EvtVector4C  Heps=H.cont2(hardCur);
 175
 176         for(int i=0; i<4; i++) {
 177                 EvtVector4C  eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
 178                 EvtComplex amp=eps*Heps;
 179                 vertex(i,amp);
 180         };
 181
 182 }
 183