TEvtGen/EvtGenModels/EvtXPsiGamma.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: EvtXPsiGamma.cc
  12 //
  13 // Description: Routine to implement radiative decay X3872(2-+) -> J/psi gamma
  14 //      according to [F. Brazzi et al, arXiv:1103.3155
  15 //
  16 // Modification history:
  17 //
  18 //    May, 7, 2012        Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtTensorParticle.hh"
  26 #include "EvtGenBase/EvtGenKine.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include "EvtGenBase/EvtVector4C.hh"
  30 #include "EvtGenBase/EvtTensor4C.hh"
  31
  32 #include "EvtGenModels/EvtXPsiGamma.hh"
  33
  34 #include <string>
  35 #include <iostream>
  36
  37 using namespace std;
  38
  39 EvtXPsiGamma::~EvtXPsiGamma() {
  40 /*  cout<<"(* AVL EvtXPsiGamma::destructor getProbMax(-1) = "<<getProbMax(-1)<<" *)"<<endl;
  41   cout<<"(* AVL EvtXPsiGamma::destructor "<<ncall<<" calls *)"<<endl;*/
  42 }
  43
  44 std::string EvtXPsiGamma::getName(){
  45   return "X38722-+_PSI_GAMMA";
  46 }
  47
  48
  49 EvtDecayBase* EvtXPsiGamma::clone(){
  50 //  cout<<" (* AVL: === EvtXPsiGamma::clone() ============ *)"<<endl;
  51   return new EvtXPsiGamma;
  52
  53 }
  54
  55 EvtComplex EvtXPsiGamma::fT2(EvtVector4R p, EvtVector4R q , EvtTensor4C epsPI, EvtVector4C epsEps, EvtVector4C epsEta) {
  56 // T2 term from [Bazi](10)
  57   EvtTensor4C epsPQ = EvtGenFunctions::directProd(q,p); // e_{mu nu a b} p^a q^b;
  58   epsPQ = dual(epsPQ);
  59
  60   EvtVector4C tmp1 = epsPI.cont1(epsEps);
  61   EvtVector4C tmp2=epsPQ.cont1(tmp1);
  62   EvtComplex T2 = tmp2*epsEta; // epa^a pi_{a mu} e_{mu nu rho si} p_nu q_rho eta_si
  63
  64   tmp1 = epsPI.cont1(epsEta);
  65   tmp2=epsPQ.cont1(tmp1);
  66   T2+=tmp2*epsEps;   // T2 - eta^a pi_{a mu} e_{mu nu rho si} q_nu p_rhi eps_si
  67
  68   return T2;
  69 }
  70
  71 EvtComplex EvtXPsiGamma::fT3(EvtVector4R p, EvtVector4R q , EvtTensor4C epsPI, EvtVector4C epsEps, EvtVector4C epsEta) {
  72 // T3 term from [Bazi](11)
  73   EvtVector4R   Q = p-q, P = p+q;
  74   EvtVector4C tmp1 = epsPI.cont1(Q); // Q_a pi_{a mu}
  75   EvtTensor4C tmp3 = dual(EvtGenFunctions::directProd(P,epsEps)); // e_{mu nu rho si} P^rho eps^si
  76   EvtVector4C tmp4 = tmp3.cont1(tmp1);
  77   EvtComplex T3 = tmp4*epsEta; // Q_a pi_{a mu} e_{mu nu rho si} P^rho eps_si eta_nu
  78   return T3;
  79 }
  80
  81
  82 void EvtXPsiGamma::decay( EvtParticle *root ){
  83   ncall++;
  84   root -> initializePhaseSpace(getNDaug(),getDaugs());
  85
  86   double gOmega = 1.58, gPOmega = -0.74;  // X -> omega psi couplings from table II
  87   double gRho = 1.58, gPRho = -0.74;  // X -> omega psi couplings from table II
  88   double fRho=0.121, mRho2 = 0.770*0.770, fOmega=0.036, mOmega2 = 0.782*0.782;
  89
  90   EvtComplex amp;
  91
  92   if(_ID0 == EvtPDL::getId("gamma") ) {
  93     for(int iPsi = 0; iPsi < 4; iPsi++) {
  94       for(int iGamma = 0; iGamma < 1; iGamma++) {
  95         for(int iChi = 0; iChi<4; iChi++) {
  96
  97           EvtComplex T2 = fT2(
  98                               root->getDaug(1)->getP4(),
  99                               root->getDaug(0)->getP4(),
 100                               root->epsTensor(iChi),
 101                               root->getDaug(1)->epsParent(iPsi).conj(),
 102                               root->getDaug(0)->epsParentPhoton(iGamma).conj()
 103                               );
 104           EvtComplex T3 = fT3(
 105                               root->getDaug(1)->getP4(),
 106                               root->getDaug(0)->getP4(),
 107                               root->epsTensor(iChi),
 108                               root->getDaug(1)->epsParent(iPsi).conj(),
 109                               root->getDaug(0)->epsParentPhoton(iGamma).conj()
 110                               );
 111           amp = (fOmega/mOmega2*gOmega+fRho/mRho2*gRho)*T2
 112             + (fOmega/mOmega2*gPOmega+fRho/mRho2*gPRho)*T3;
 113           vertex(iChi, iGamma, iPsi, amp);
 114         };};};
 115   }
 116   else if(_ID0 == EvtPDL::getId("omega") ) {
 117     for(int iPsi = 0; iPsi < 4; iPsi++) {
 118       for(int iGamma = 0; iGamma < 4; iGamma++) {
 119         for(int iChi = 0; iChi<4; iChi++) {
 120
 121           EvtComplex T2 = fT2(
 122                               root->getDaug(1)->getP4(),
 123                               root->getDaug(0)->getP4(),
 124                               root->epsTensor(iChi),
 125                               root->getDaug(1)->epsParent(iPsi).conj(),
 126                               root->getDaug(0)->epsParent(iGamma).conj()
 127                               );
 128           EvtComplex T3 = fT3(
 129                               root->getDaug(1)->getP4(),
 130                               root->getDaug(0)->getP4(),
 131                               root->epsTensor(iChi),
 132                               root->getDaug(1)->epsParent(iPsi).conj(),
 133                               root->getDaug(0)->epsParent(iGamma).conj()
 134                               );
 135           //      cout << "AVL:: omega"<<endl;
 136           amp = gOmega*T2 + gPOmega*T3;
 137           vertex(iChi, iGamma, iPsi, amp);
 138         };};};
 139   }
 140   else if(_ID0 == EvtPDL::getId("rho0") ) {
 141     for(int iPsi = 0; iPsi < 4; iPsi++) {
 142       for(int iGamma = 0; iGamma < 4; iGamma++) {
 143         for(int iChi = 0; iChi<4; iChi++) {
 144
 145           EvtComplex T2 = fT2(
 146                               root->getDaug(1)->getP4(),
 147                               root->getDaug(0)->getP4(),
 148                               root->epsTensor(iChi),
 149                               root->getDaug(1)->epsParent(iPsi).conj(),
 150                               root->getDaug(0)->epsParent(iGamma).conj()
 151                               );
 152           EvtComplex T3 = fT3(
 153                               root->getDaug(1)->getP4(),
 154                               root->getDaug(0)->getP4(),
 155                               root->epsTensor(iChi),
 156                               root->getDaug(1)->epsParent(iPsi).conj(),
 157                               root->getDaug(0)->epsParent(iGamma).conj()
 158                               );
 159           //      cout << "AVL:: rho"<<endl;
 160           amp = gRho*T2 + gPRho*T3;
 161           vertex(iChi, iGamma, iPsi, amp);
 162         };};};
 163   }
 164   else {
 165     cout<<"AVL:: Not realized yet"<<endl;
 166   };
 167
 168 }
 169
 170
 171 void EvtXPsiGamma::init(){
 172 //  cout<<" (* AVL: ==== EvtXPsiGamma::init() ============ *)"<<endl;
 173
 174   ncall = 0;
 175
 176   checkNArg(0);
 177   checkNDaug(2);
 178
 179
 180   checkSpinParent(EvtSpinType::TENSOR);
 181
 182 //  checkSpinDaughter(0,EvtSpinType::PHOTON);
 183   checkSpinDaughter(1,EvtSpinType::VECTOR);
 184
 185   _ID0 = getDaug(0);
 186 /*  if(_ID0 == EvtPDL::getId("gamma") ) {
 187     cout << "AVL:: gamma"<<endl;
 188   }
 189   else if(_ID0 == EvtPDL::getId("omega") ) {
 190     cout << "AVL:: omega"<<endl;
 191   }
 192   else if(_ID0 == EvtPDL::getId("rho0") ) {
 193     cout << "AVL:: rho"<<endl;
 194   };
 195 */
 196
 197 }
 198
 199 void EvtXPsiGamma::initProbMax() {
 200   if(_ID0 == EvtPDL::getId("gamma") )  setProbMax(2.400);
 201   else if(_ID0 == EvtPDL::getId("omega") )  setProbMax(16.);
 202   else if(_ID0 == EvtPDL::getId("rho0") )  setProbMax(70.);
 203 };
 204
 205