TEvtGen/EvtGenModels/EvtLNuGamma.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) 2001      Caltech
  10 //
  11 // Module: EvtLNuGamma.cc
  12 //
  13 // Description: B+ -> l+ nu gamma
  14 //
  15 //
  16 // Modification history:
  17 //
  18 //    Edward Chen   April 24, 2001           Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include <stdlib.h>
  24 #include <iostream>
  25 #include <string>
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtGenKine.hh"
  29 #include "EvtGenModels/EvtLNuGamma.hh"
  30 #include "EvtGenBase/EvtDiracSpinor.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32 #include "EvtGenBase/EvtComplex.hh"
  33 #include "EvtGenBase/EvtVector4C.hh"
  34 #include "EvtGenBase/EvtVector4R.hh"
  35 #include "EvtGenBase/EvtTensor4C.hh"
  36
  37 EvtLNuGamma::EvtLNuGamma(){
  38   _fafvzero = false;
  39 }
  40
  41 EvtLNuGamma::~EvtLNuGamma() {}
  42
  43 std::string EvtLNuGamma::getName(){
  44
  45   return "LNUGAMMA";
  46
  47 }
  48
  49
  50 EvtDecayBase* EvtLNuGamma::clone(){
  51
  52   return new EvtLNuGamma;
  53
  54 }
  55
  56 void EvtLNuGamma::init(){
  57
  58   // check that there are 3 or 4 arguments
  59   checkNArg(3,4);
  60   checkNDaug(3);
  61
  62   if (getNArg() == 4){
  63     //      Argv[3] is a flag set to 0 if abs(f_a/f_v) is 1
  64     //       and not set to 0 if f_a/f_v is set to 0.
  65     if (getArg(3) > 0){
  66       _fafvzero = true;
  67     }
  68     else{
  69       _fafvzero = false;
  70     }
  71   }
  72   else{
  73     _fafvzero = false;
  74   }
  75
  76   checkSpinParent(EvtSpinType::SCALAR);
  77
  78   checkSpinDaughter(0,EvtSpinType::DIRAC);
  79   checkSpinDaughter(1,EvtSpinType::NEUTRINO);
  80   checkSpinDaughter(2,EvtSpinType::PHOTON);
  81 }
  82
  83 void EvtLNuGamma::initProbMax(){
  84
  85
  86   setProbMax(7000.0);
  87
  88 }
  89
  90 void EvtLNuGamma::decay(EvtParticle *p){
  91
  92   static EvtId BM=EvtPDL::getId("B-");
  93   static EvtId DM=EvtPDL::getId("D-");
  94   p->initializePhaseSpace(getNDaug(),getDaugs());
  95
  96   EvtComplex myI(0,1);
  97
  98   EvtParticle *lept, *neut,*phot;
  99   lept=p->getDaug(0);
 100   neut=p->getDaug(1);
 101   phot=p->getDaug(2);
 102
 103   EvtVector4C lept1,lept2,photon1,photon2;
 104
 105   if (p->getId()==BM||p->getId()==DM) {
 106     lept1=EvtLeptonVACurrent(lept->spParent(0),neut->spParentNeutrino());
 107     lept2=EvtLeptonVACurrent(lept->spParent(1),neut->spParentNeutrino());
 108   }
 109   else{
 110     lept1=EvtLeptonVACurrent(neut->spParentNeutrino(),lept->spParent(0));
 111     lept2=EvtLeptonVACurrent(neut->spParentNeutrino(),lept->spParent(1));
 112   }
 113
 114
 115
 116   EvtVector4R photp = phot->getP4();  // Photon 4-momentum in parent rest frame
 117   double photE = photp.get(0);  // Photon energy in parent rest frame
 118
 119   EvtVector4C photone1 = phot->epsParentPhoton(0).conj();
 120   EvtVector4C photone2 = phot->epsParentPhoton(1).conj();
 121
 122   EvtVector4R parVelocity(1,0,0,0);  // Parent velocity in parent rest-frame
 123
 124   double fv,fa;
 125
 126   fv = getFormFactor(photE);
 127   if (_fafvzero){
 128     fa = 0.0;
 129   }
 130   else if (p->getId()==BM||p->getId()==DM) {
 131     fa = - fv;
 132   }
 133   else{
 134     fa = fv;
 135   }
 136
 137   EvtVector4C temp1a = dual(directProd(parVelocity,photp)).cont2(photone1);
 138   EvtVector4C temp2a = dual(directProd(parVelocity,photp)).cont2(photone2);
 139
 140   EvtVector4C temp1b = (photone1)*(parVelocity*photp);
 141   EvtVector4C temp1c = (photp)*(photone1*parVelocity);
 142
 143   EvtVector4C temp2b = (photone2)*(parVelocity*photp);
 144   EvtVector4C temp2c = (photp)*(photone2*parVelocity);
 145
 146   photon1 = (temp1a*fv) + (myI*fa*(temp1b - temp1c));
 147   photon2 = (temp2a*fv) + (myI*fa*(temp2b - temp2c));
 148
 149   vertex(0,0,lept1.cont(photon1));
 150   vertex(0,1,lept1.cont(photon2));
 151   vertex(1,0,lept2.cont(photon1));
 152   vertex(1,1,lept2.cont(photon2));
 153
 154   return;
 155
 156 }
 157
 158 double EvtLNuGamma::getFormFactor(double photonEnergy){
 159   // Arg[0] = photon mass cutoff (GeV)
 160   // Arg[1] = R (GeV^(-1))
 161   // Arg[2] = m_b (GeV)
 162   // Using Korchemsky et al. Phy Rev D 61 (2000) 114510
 163   // Up to a constant
 164
 165   double formFactor = 0;
 166   double qu = 2./3.;
 167   double qb = -1./3.;
 168
 169   if (photonEnergy > getArg(0)){
 170     formFactor = (1/photonEnergy) * ((qu*getArg(1)) - (qb/getArg(2)));
 171   }
 172   return formFactor;
 173 }