00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <stdlib.h>
00024 #include <iostream>
00025 #include <string>
00026 #include "EvtGenBase/EvtParticle.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenBase/EvtGenKine.hh"
00029 #include "EvtGenModels/EvtTauHadnu.hh"
00030 #include "EvtGenBase/EvtDiracSpinor.hh"
00031 #include "EvtGenBase/EvtReport.hh"
00032 #include "EvtGenBase/EvtVector4C.hh"
00033 #include "EvtGenBase/EvtIdSet.hh"
00034 using std::endl;
00035
00036 EvtTauHadnu::~EvtTauHadnu() {}
00037
00038 void EvtTauHadnu::getName(std::string& model_name){
00039
00040 model_name="TAUHADNU";
00041
00042 }
00043
00044
00045 EvtDecayBase* EvtTauHadnu::clone(){
00046
00047 return new EvtTauHadnu;
00048
00049 }
00050
00051 void EvtTauHadnu::init() {
00052
00053
00054
00055 checkSpinParent(EvtSpinType::DIRAC);
00056
00057
00058 checkSpinDaughter(getNDaug()-1,EvtSpinType::NEUTRINO);
00059
00060 int i;
00061 for ( i=0; i<(getNDaug()-1);i++) {
00062 checkSpinDaughter(i,EvtSpinType::SCALAR);
00063 }
00064
00065 bool validndaug=false;
00066
00067 if ( getNDaug()==4 ) {
00068
00069 validndaug=true;
00070 checkNArg(7);
00071 _beta = getArg(0);
00072 _mRho = getArg(1);
00073 _gammaRho = getArg(2);
00074 _mRhopr = getArg(3);
00075 _gammaRhopr = getArg(4);
00076 _mA1 = getArg(5);
00077 _gammaA1 = getArg(6);
00078 }
00079 if ( getNDaug()==3 ) {
00080
00081 validndaug=true;
00082 checkNArg(5);
00083 _beta = getArg(0);
00084 _mRho = getArg(1);
00085 _gammaRho = getArg(2);
00086 _mRhopr = getArg(3);
00087 _gammaRhopr = getArg(4);
00088 }
00089 if ( getNDaug()==2 ) {
00090
00091 validndaug=true;
00092 checkNArg(0);
00093 }
00094
00095 if ( !validndaug ) {
00096 report(ERROR,"EvtGen") << "Have not yet implemented this final state in TAUHADNU model" << endl;
00097 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
00098 int id;
00099 for ( id=0; id<(getNDaug()-1); id++ )
00100 report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
00101
00102 }
00103
00104 }
00105
00106 void EvtTauHadnu::initProbMax(){
00107
00108 if ( getNDaug()==2 ) setProbMax(90.0);
00109 if ( getNDaug()==3 ) setProbMax(2500.0);
00110 if ( getNDaug()==4 ) setProbMax(5000.0);
00111
00112 }
00113
00114 void EvtTauHadnu::decay(EvtParticle *p){
00115
00116 static EvtId TAUM=EvtPDL::getId("tau-");
00117
00118 EvtIdSet thePis("pi+","pi-","pi0");
00119 EvtIdSet theKs("K+","K-");
00120
00121 p->initializePhaseSpace(getNDaug(),getDaugs());
00122
00123 EvtParticle *nut;
00124 nut = p->getDaug(getNDaug()-1);
00125
00126
00127 EvtVector4C tau1, tau2;
00128
00129 if (p->getId()==TAUM) {
00130 tau1=EvtLeptonVACurrent(nut->spParentNeutrino(),p->sp(0));
00131 tau2=EvtLeptonVACurrent(nut->spParentNeutrino(),p->sp(1));
00132 }
00133 else{
00134 tau1=EvtLeptonVACurrent(p->sp(0),nut->spParentNeutrino());
00135 tau2=EvtLeptonVACurrent(p->sp(1),nut->spParentNeutrino());
00136 }
00137
00138 EvtVector4C hadCurr;
00139 bool foundHadCurr=false;
00140 if ( getNDaug() == 2 ) {
00141 hadCurr = p->getDaug(0)->getP4();
00142 foundHadCurr=true;
00143 }
00144 if ( getNDaug() == 3 ) {
00145
00146
00147 if ( thePis.contains(getDaug(0)) &&
00148 thePis.contains(getDaug(1)) ) {
00149
00150 EvtVector4R q1 = p->getDaug(0)->getP4();
00151 EvtVector4R q2 = p->getDaug(1)->getP4();
00152
00153 hadCurr = Fpi(q1,q2)*(q1-q2);
00154
00155 foundHadCurr = true;
00156 }
00157
00158 }
00159 if ( getNDaug() == 4 ) {
00160 if ( thePis.contains(getDaug(0)) &&
00161 thePis.contains(getDaug(1)) &&
00162 thePis.contains(getDaug(2)) ) {
00163 foundHadCurr = true;
00164
00165
00166
00167 int diffPi(0),samePi1(0),samePi2(0);
00168 if ( getDaug(0) == getDaug(1) ) {diffPi=2; samePi1=0; samePi2=1;}
00169 if ( getDaug(0) == getDaug(2) ) {diffPi=1; samePi1=0; samePi2=2;}
00170 if ( getDaug(1) == getDaug(2) ) {diffPi=0; samePi1=1; samePi2=2;}
00171
00172 EvtVector4R q1=p->getDaug(samePi1)->getP4();
00173 EvtVector4R q2=p->getDaug(samePi2)->getP4();
00174 EvtVector4R q3=p->getDaug(diffPi)->getP4();
00175
00176 EvtVector4R Q=q1+q2+q3;
00177 double qMass2=Q.mass2();
00178
00179 double GA1=_gammaA1*pi3G(Q.mass2(),samePi1)/pi3G(_mA1*_mA1,samePi1);
00180
00181 EvtComplex denBA1(_mA1*_mA1 - Q.mass2(),-1.*_mA1*GA1);
00182 EvtComplex BA1 = _mA1*_mA1 / denBA1;
00183
00184 hadCurr = BA1*( (q1-q3) - (Q*(Q*(q1-q3))/qMass2)*Fpi(q2,q3) +
00185 (q2-q3) - (Q*(Q*(q2-q3))/qMass2)*Fpi(q1,q3) );
00186
00187
00188 }
00189
00190
00191 }
00192
00193
00194
00195 if ( !foundHadCurr ) {
00196 report(ERROR,"EvtGen") << "Have not yet implemented this final state in TAUHADNU model" << endl;
00197 report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
00198 int id;
00199 for ( id=0; id<(getNDaug()-1); id++ )
00200 report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
00201
00202 }
00203
00204
00205 vertex(0,tau1*hadCurr);
00206 vertex(1,tau2*hadCurr);
00207
00208 return;
00209
00210 }
00211
00212
00213 double EvtTauHadnu::pi3G(double m2,int dupD) {
00214 double mPi= EvtPDL::getMeanMass(getDaug(dupD));
00215 if ( m2 > (_mRho+mPi) ) {
00216 return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
00217 }
00218 else {
00219 double t1=m2-9.0*mPi*mPi;
00220 return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
00221 }
00222 return 0.;
00223 }
00224
00225 EvtComplex EvtTauHadnu::Fpi( EvtVector4R q1, EvtVector4R q2) {
00226
00227 double m1=q1.mass();
00228 double m2=q2.mass();
00229
00230 EvtVector4R Q = q1 + q2;
00231 double mQ2= Q*Q;
00232
00233
00234 double dRho= _mRho*_mRho - m1*m1 - m2*m2;
00235 double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
00236
00237 double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
00238 double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
00239
00240 double dQ= mQ2 - m1*m1 - m2*m2;
00241 double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
00242
00243
00244 double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
00245 EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
00246 EvtComplex BRho= _mRho*_mRho / BRhoDem;
00247
00248 double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
00249 EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
00250 EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
00251
00252 return (BRho + _beta*BRhopr)/(1+_beta);
00253 }
00254
00255
00256
00257
00258