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: EvtTVP.cc
12 : //
13 : // Description: Routine to calculate W -> (n pi) current
14 : // according to [Kuhn, Was, Acta.Phys.Polon B39 (2008) 147]
15 : //
16 : // Modification history:
17 : // AVL 6 July, 2012 Module created
18 : //
19 : //------------------------------------------------------------------------
20 : //
21 : #include "EvtGenBase/EvtPatches.hh"
22 : #include <iostream>
23 : #include <iomanip>
24 : #include <fstream>
25 : #include <ctype.h>
26 : #include <stdlib.h>
27 : #include <string.h>
28 : #include "EvtGenBase/EvtParticle.hh"
29 : #include "EvtGenBase/EvtPDL.hh"
30 : #include "EvtGenBase/EvtGenKine.hh"
31 : #include "EvtGenModels/EvtTauHadnu.hh"
32 : #include "EvtGenBase/EvtDiracSpinor.hh"
33 : #include "EvtGenBase/EvtReport.hh"
34 : #include "EvtGenBase/EvtVector4C.hh"
35 : #include "EvtGenBase/EvtTensor4C.hh"
36 : #include "EvtGenBase/EvtIdSet.hh"
37 : #include "EvtGenBase/EvtParser.hh"
38 :
39 : #include "EvtGenModels/EvtWnPi.hh"
40 :
41 : using namespace std;
42 :
43 : // W+ -> pi_ current
44 : EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1) {
45 0 : return q1;
46 : }
47 :
48 : // W+ -> pi+ pi0 current
49 : EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2) {
50 0 : return BWr(q1+q2)*(q1-q2);
51 : }
52 :
53 : // W+ -> pi+ pi+ pi- current
54 : EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3) {
55 0 : EvtVector4R Q=q1+q2+q3;
56 0 : double Q2=Q.mass2();
57 0 : return BWa(Q)*( (q1-q3) - (Q*(Q*(q1-q3))/Q2)*BWr(q2+q3) +
58 0 : (q2-q3) - (Q*(Q*(q2-q3))/Q2)*BWr(q1+q3) );
59 0 : }
60 :
61 : // W+ -> pi+ pi+ pi- pi- pi+ current with symmetrization
62 : EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3, EvtVector4R q4, EvtVector4R q5) {
63 : // double Q2 = Qtot*Qtot;
64 : // return q1-Qtot*(q1*Qtot)/Q2;
65 0 : EvtVector4C V = JB(q1, q2, q3, q4, q5) + JB(q5, q2, q3, q4, q1) + JB(q1, q5, q3, q4, q2) +
66 0 : JB(q1,q2,q4,q3,q5)+JB(q5,q2,q4,q3,q1)+JB(q1,q5,q4,q3,q2);
67 : // cout<<"BC2: Qtot="<<Qtot<<", V="<<V<<endl;
68 : return V;
69 0 : }
70 :
71 :
72 : // a1 -> pi+ pi+ pi- BW
73 : EvtComplex EvtWnPi::BWa(EvtVector4R q) {
74 : double const _mA1=1.26, _GA1=0.4;
75 0 : EvtComplex I(0,1);
76 0 : double Q2 = q.mass2();
77 0 : double GA1=_GA1*pi3G(Q2)/pi3G(_mA1*_mA1);
78 0 : EvtComplex denBA1(_mA1*_mA1 - Q2,-1.*_mA1*GA1);
79 0 : return _mA1*_mA1 / denBA1;
80 0 : }
81 :
82 :
83 : EvtComplex EvtWnPi::BWf(EvtVector4R q) {
84 : double const mf=0.8, Gf=0.6;
85 0 : EvtComplex I(0,1);
86 0 : double Q2 = q.mass2();
87 0 : return mf*mf/(mf*mf-Q2-I*mf*Gf);
88 0 : }
89 :
90 : EvtComplex EvtWnPi::BWr(EvtVector4R q) {
91 : double _mRho = 0.775, _gammaRho=0.149, _mRhopr=1.364, _gammaRhopr=0.400, _beta=-0.108;
92 0 : double m1=EvtPDL::getMeanMass(EvtPDL::getId("pi+")), m2=EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
93 0 : double mQ2=q.mass2();
94 :
95 : // momenta in the rho->pipi decay
96 0 : double dRho= _mRho*_mRho - m1*m1 - m2*m2;
97 0 : double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
98 :
99 0 : double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
100 0 : double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
101 :
102 0 : double dQ= mQ2 - m1*m1 - m2*m2;
103 0 : double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
104 :
105 :
106 0 : double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
107 0 : EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
108 0 : EvtComplex BRho= _mRho*_mRho / BRhoDem;
109 :
110 0 : double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
111 0 : EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
112 0 : EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
113 :
114 0 : return (BRho + _beta*BRhopr)/(1+_beta);
115 0 : }
116 :
117 : double EvtWnPi::pi3G(double m2) {
118 0 : double mPi = EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
119 : double _mRho = 0.775;
120 0 : if ( m2 > (_mRho+mPi) ) {
121 0 : return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
122 : }
123 : else {
124 0 : double t1=m2-9.0*mPi*mPi;
125 0 : return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
126 : };
127 0 : }
128 :
129 : EvtVector4C EvtWnPi::JB( EvtVector4R p1, EvtVector4R p2, EvtVector4R p3, EvtVector4R p4, EvtVector4R p5) {
130 0 : EvtVector4R Qtot = p1+p2+p3+p4+p5, Qa=p1+p2+p3;
131 0 : EvtTensor4C T= (1/Qtot.mass2())*EvtGenFunctions::directProd(Qtot,Qtot) - EvtTensor4C::g();
132 0 : EvtVector4R V13 = Qa*( p2*(p1-p3) )/Qa.mass2() - (p1-p3);
133 0 : EvtVector4R V23 = Qa*( p1*(p2-p3) )/Qa.mass2() - (p2-p3);
134 0 : return BWa(Qtot)*BWa(Qa)*BWf(p4+p5)*(
135 0 : T.cont1(V13)*BWr(p1+p3) + T.cont1(V23)*BWr(p2+p3)
136 : );
137 0 : }
|
