TEvtGen/EvtGenModels/EvtTauHadnu.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) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtTauHadnu.cc
  12 //
  13 // Description: The leptonic decay of the tau meson.
  14 //              E.g., tau- -> e- nueb nut
  15 //
  16 // Modification history:
  17 //
  18 //    RYD       January 17, 1997       Module created
  19 //
  20 //------------------------------------------------------------------------
  21 //
  22 #include <stdlib.h>
  23 #include <iostream>
  24 #include <string>
  25 #include "EvtGenBase/EvtParticle.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtGenKine.hh"
  28 #include "EvtGenModels/EvtTauHadnu.hh"
  29 #include "EvtGenBase/EvtDiracSpinor.hh"
  30 #include "EvtGenBase/EvtReport.hh"
  31 #include "EvtGenBase/EvtVector4C.hh"
  32 #include "EvtGenBase/EvtIdSet.hh"
  33
  34 using namespace std;
  35
  36 EvtTauHadnu::~EvtTauHadnu() {}
  37
  38 std::string EvtTauHadnu::getName(){
  39
  40   return "TAUHADNU";
  41
  42 }
  43
  44
  45 EvtDecayBase* EvtTauHadnu::clone(){
  46
  47   return new EvtTauHadnu;
  48
  49 }
  50
  51 void EvtTauHadnu::init() {
  52
  53   // check that there are 0 arguments
  54
  55   checkSpinParent(EvtSpinType::DIRAC);
  56
  57   //the last daughter should be a neutrino
  58   checkSpinDaughter(getNDaug()-1,EvtSpinType::NEUTRINO);
  59
  60   int i;
  61   for ( i=0; i<(getNDaug()-1);i++) {
  62     checkSpinDaughter(i,EvtSpinType::SCALAR);
  63   }
  64
  65   bool validndaug=false;
  66
  67   if ( getNDaug()==4 ) {
  68     //pipinu
  69     validndaug=true;
  70     checkNArg(7);
  71     _beta = getArg(0);
  72     _mRho = getArg(1);
  73     _gammaRho = getArg(2);
  74     _mRhopr = getArg(3);
  75     _gammaRhopr = getArg(4);
  76     _mA1 = getArg(5);
  77     _gammaA1 = getArg(6);
  78   }
  79   if ( getNDaug()==3 ) {
  80     //pipinu
  81     validndaug=true;
  82     checkNArg(5);
  83     _beta = getArg(0);
  84     _mRho = getArg(1);
  85     _gammaRho = getArg(2);
  86     _mRhopr = getArg(3);
  87     _gammaRhopr = getArg(4);
  88   }
  89   if ( getNDaug()==2 ) {
  90     //pipinu
  91     validndaug=true;
  92     checkNArg(0);
  93   }
  94
  95   if ( !validndaug ) {
  96     report(ERROR,"EvtGen") << "Have not yet implemented this final state in TAUHADNUKS model" << endl;
  97     report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
  98     int id;
  99     for ( id=0; id<(getNDaug()-1); id++ )
 100       report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
 101
 102   }
 103
 104 }
 105
 106 void EvtTauHadnu::initProbMax(){
 107
 108   if ( getNDaug()==2 )  setProbMax(90.0);
 109   if ( getNDaug()==3 )  setProbMax(2500.0);
 110   if ( getNDaug()==4 )  setProbMax(30000.0);
 111
 112 }
 113
 114 void EvtTauHadnu::decay(EvtParticle *p){
 115
 116   static EvtId TAUM=EvtPDL::getId("tau-");
 117
 118   EvtIdSet thePis("pi+","pi-","pi0");
 119   EvtIdSet theKs("K+","K-");
 120
 121   p->initializePhaseSpace(getNDaug(),getDaugs());
 122
 123   EvtParticle *nut;
 124   nut = p->getDaug(getNDaug()-1);
 125   EvtVector4R momscalar = p->getDaug(0)->getP4();
 126
 127   //get the leptonic current
 128   EvtVector4C tau1, tau2;
 129
 130   if (p->getId()==TAUM) {
 131     tau1=EvtLeptonVACurrent(nut->spParentNeutrino(),p->sp(0));
 132     tau2=EvtLeptonVACurrent(nut->spParentNeutrino(),p->sp(1));
 133   }
 134   else{
 135     tau1=EvtLeptonVACurrent(p->sp(0),nut->spParentNeutrino());
 136     tau2=EvtLeptonVACurrent(p->sp(1),nut->spParentNeutrino());
 137   }
 138
 139   EvtVector4C hadCurr;
 140   bool foundHadCurr=false;
 141   if ( getNDaug() == 2 ) {
 142     hadCurr = p->getDaug(0)->getP4();
 143     foundHadCurr=true;
 144   }
 145   if ( getNDaug() == 3 ) {
 146
 147     //pi pi0 nu with rho and rhopr resonance
 148     if ( thePis.contains(getDaug(0)) &&
 149          thePis.contains(getDaug(1)) ) {
 150
 151       EvtVector4R q1 = p->getDaug(0)->getP4();
 152       EvtVector4R q2 = p->getDaug(1)->getP4();
 153       double m1 = q1.mass();
 154       double m2 = q2.mass();
 155
 156       hadCurr = Fpi( (q1+q2).mass2(), m1, m2 )  * (q1-q2);
 157
 158       foundHadCurr = true;
 159     }
 160
 161   }
 162   if ( getNDaug() == 4 ) {
 163     if ( thePis.contains(getDaug(0)) &&
 164          thePis.contains(getDaug(1)) &&
 165          thePis.contains(getDaug(2)) ) {
 166       //figure out which is the different charged pi
 167       //want it to be q3
 168
 169       int diffPi(0),samePi1(0),samePi2(0);
 170       if ( getDaug(0) == getDaug(1) ) {diffPi=2; samePi1=0; samePi2=1;}
 171       if ( getDaug(0) == getDaug(2) ) {diffPi=1; samePi1=0; samePi2=2;}
 172       if ( getDaug(1) == getDaug(2) ) {diffPi=0; samePi1=1; samePi2=2;}
 173
 174       EvtVector4R q1=p->getDaug(samePi1)->getP4();
 175       EvtVector4R q2=p->getDaug(samePi2)->getP4();
 176       EvtVector4R q3=p->getDaug(diffPi)->getP4();
 177
 178       double m1 = q1.mass();
 179       double m2 = q2.mass();
 180       double m3 = q3.mass();
 181
 182       EvtVector4R Q = q1 + q2 + q3;
 183       double Q2 = Q.mass2();
 184       double _mA12 = _mA1*_mA1;
 185
 186       double _gammaA1X = _gammaA1 * gFunc( Q2, samePi1 )/gFunc( _mA12, samePi1 );
 187
 188       EvtComplex denBW_A1( _mA12 - Q2, -1.*_mA1*_gammaA1X );
 189       EvtComplex BW_A1 = _mA12 / denBW_A1;
 190
 191       hadCurr = BW_A1 * ( ((q1-q3)-(Q*(Q*(q1-q3))/Q2)) * Fpi( (q1+q3).mass2(), m1, m3) +
 192                           ((q2-q3)-(Q*(Q*(q2-q3))/Q2)) * Fpi( (q2+q3).mass2(), m2, m3) );
 193
 194       foundHadCurr = true;
 195
 196     }
 197
 198
 199   }
 200
 201
 202
 203   if ( !foundHadCurr ) {
 204     report(ERROR,"EvtGen") << "Have not yet implemented this final state in TAUHADNUKS model" << endl;
 205     report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
 206     int id;
 207     for ( id=0; id<(getNDaug()-1); id++ )
 208       report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
 209
 210   }
 211
 212
 213   vertex(0,tau1*hadCurr);
 214   vertex(1,tau2*hadCurr);
 215
 216
 217
 218   return;
 219
 220 }
 221
 222 double EvtTauHadnu::gFunc(double Q2, int dupD) {
 223
 224   double mpi= EvtPDL::getMeanMass(getDaug(dupD));
 225   double mpi2 = pow( mpi,2.);
 226   if ( Q2 < pow(_mRho + mpi, 2.) ) {
 227     double arg = Q2-9.*mpi2;
 228     return 4.1 * pow(arg,3.) * (1. - 3.3*arg + 5.8*pow(arg,2.));
 229   }
 230   else
 231     return Q2 * (1.623 + 10.38/Q2 - 9.32/pow(Q2,2.) + 0.65/pow(Q2,3.));
 232 }
 233
 234 EvtComplex EvtTauHadnu::Fpi( double s, double xm1, double xm2 ) {
 235
 236   EvtComplex BW_rho = BW( s, _mRho, _gammaRho, xm1, xm2 );
 237   EvtComplex BW_rhopr = BW( s, _mRhopr, _gammaRhopr, xm1, xm2 );
 238
 239
 240   return (BW_rho + _beta*BW_rhopr)/(1.+_beta);
 241 }
 242
 243 EvtComplex EvtTauHadnu::BW( double s, double m, double gamma, double xm1, double xm2 ) {
 244
 245   double m2 = pow( m, 2.);
 246
 247   if ( s > pow( xm1+xm2, 2.) ) {
 248     double qs = sqrt( fabs( (s-pow(xm1+xm2,2.)) * (s-pow(xm1-xm2,2.)) ) ) / sqrt(s);
 249     double qm = sqrt( fabs( (m2-pow(xm1+xm2,2.)) * (m2-pow(xm1-xm2,2.)) ) ) / m;
 250
 251     gamma *= m2/s * pow( qs/qm, 3.);
 252   }
 253   else
 254     gamma = 0.;
 255
 256   EvtComplex denBW( m2 - s, -1.* sqrt(s) * gamma );
 257
 258
 259   return m2 / denBW;
 260 }
 261
 262
 263
 264
 265