Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtWnPi.cpp
CommitLineData
0ca57c2f 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
41using namespace std;
42
43// W+ -> pi_ current
44EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1) {
45 return q1;
46}
47
48// W+ -> pi+ pi0 current
49EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2) {
50 return BWr(q1+q2)*(q1-q2);
51}
52
53// W+ -> pi+ pi+ pi- current
54EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3) {
55 EvtVector4R Q=q1+q2+q3;
56 double Q2=Q.mass2();
57 return BWa(Q)*( (q1-q3) - (Q*(Q*(q1-q3))/Q2)*BWr(q2+q3) +
58 (q2-q3) - (Q*(Q*(q2-q3))/Q2)*BWr(q1+q3) );
59}
60
61// W+ -> pi+ pi+ pi- pi- pi+ current with symmetrization
62EvtVector4C 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 EvtVector4C V = JB(q1, q2, q3, q4, q5) + JB(q5, q2, q3, q4, q1) + JB(q1, q5, q3, q4, q2) +
66 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}
70
71
72// a1 -> pi+ pi+ pi- BW
73EvtComplex EvtWnPi::BWa(EvtVector4R q) {
74 double const _mA1=1.26, _GA1=0.4;
75 EvtComplex I(0,1);
76 double Q2 = q.mass2();
77 double GA1=_GA1*pi3G(Q2)/pi3G(_mA1*_mA1);
78 EvtComplex denBA1(_mA1*_mA1 - Q2,-1.*_mA1*GA1);
79 return _mA1*_mA1 / denBA1;
80}
81
82
83EvtComplex EvtWnPi::BWf(EvtVector4R q) {
84 double const mf=0.8, Gf=0.6;
85 EvtComplex I(0,1);
86 double Q2 = q.mass2();
87 return mf*mf/(mf*mf-Q2-I*mf*Gf);
88}
89
90EvtComplex EvtWnPi::BWr(EvtVector4R q) {
91 double _mRho = 0.775, _gammaRho=0.149, _mRhopr=1.364, _gammaRhopr=0.400, _beta=-0.108;
92 double m1=EvtPDL::getMeanMass(EvtPDL::getId("pi+")), m2=EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
93 double mQ2=q.mass2();
94
95 // momenta in the rho->pipi decay
96 double dRho= _mRho*_mRho - m1*m1 - m2*m2;
97 double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
98
99 double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
100 double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
101
102 double dQ= mQ2 - m1*m1 - m2*m2;
103 double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
104
105
106 double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
107 EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
108 EvtComplex BRho= _mRho*_mRho / BRhoDem;
109
110 double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
111 EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
112 EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
113
114 return (BRho + _beta*BRhopr)/(1+_beta);
115}
116
117double EvtWnPi::pi3G(double m2) {
118 double mPi = EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
119 double _mRho = 0.775;
120 if ( m2 > (_mRho+mPi) ) {
121 return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
122 }
123 else {
124 double t1=m2-9.0*mPi*mPi;
125 return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
126 };
127}
128
129EvtVector4C EvtWnPi::JB( EvtVector4R p1, EvtVector4R p2, EvtVector4R p3, EvtVector4R p4, EvtVector4R p5) {
130 EvtVector4R Qtot = p1+p2+p3+p4+p5, Qa=p1+p2+p3;
131 EvtTensor4C T= (1/Qtot.mass2())*EvtGenFunctions::directProd(Qtot,Qtot) - EvtTensor4C::g();
132 EvtVector4R V13 = Qa*( p2*(p1-p3) )/Qa.mass2() - (p1-p3);
133 EvtVector4R V23 = Qa*( p1*(p2-p3) )/Qa.mass2() - (p2-p3);
134 return BWa(Qtot)*BWa(Qa)*BWf(p4+p5)*(
135 T.cont1(V13)*BWr(p1+p3) + T.cont1(V23)*BWr(p2+p3)
136 );
137}