Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtTauHadnu.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: 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   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