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