]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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 |