00001 #include "EvtGenBase/EvtPatches.hh"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <assert.h>
00012 #include <math.h>
00013 #include <iostream>
00014 #include "EvtGenBase/EvtComplex.hh"
00015 #include "EvtGenBase/EvtPto3PAmp.hh"
00016 #include "EvtGenBase/EvtDalitzCoord.hh"
00017 #include "EvtGenBase/EvtdFunction.hh"
00018 #include "EvtGenBase/EvtCyclic3.hh"
00019 using std::endl;
00020 using EvtCyclic3::Index;
00021 using EvtCyclic3::Pair;
00022
00023
00024 EvtPto3PAmp::EvtPto3PAmp(EvtDalitzPlot dp, Pair pairAng, Pair pairRes,
00025 EvtSpinType::spintype spin,
00026 const EvtPropagator& prop, NumType typeN)
00027 : EvtAmplitude<EvtDalitzPoint>(),
00028 _pairAng(pairAng), _pairRes(pairRes),
00029 _spin(spin),
00030 _typeN(typeN),
00031 _prop((EvtPropagator*) prop.clone()),
00032 _g0(prop.g0()),
00033 _vb(prop.m0(),dp.m(EvtCyclic3::other(pairRes)),dp.bigM(),spin),
00034 _vd(dp.m(EvtCyclic3::first(pairRes)),dp.m(EvtCyclic3::second(pairRes)),prop.m0(),spin)
00035 {}
00036
00037
00038
00039 EvtPto3PAmp::EvtPto3PAmp(const EvtPto3PAmp& other)
00040 : EvtAmplitude<EvtDalitzPoint>(other),
00041 _pairAng(other._pairAng),
00042 _pairRes(other._pairRes),
00043 _spin(other._spin),
00044 _typeN(other._typeN),
00045 _prop( (other._prop) ? (EvtPropagator*) other._prop->clone() : 0),
00046 _g0(other._g0),
00047 _vb(other._vb), _vd(other._vd)
00048 {}
00049
00050
00051 EvtPto3PAmp::~EvtPto3PAmp()
00052 {
00053 if(_prop) delete _prop;
00054 }
00055
00056
00057 void EvtPto3PAmp::set_fd(double R)
00058 {
00059 _vd.set_f(R);
00060 }
00061
00062 void EvtPto3PAmp::set_fb(double R)
00063 {
00064 _vb.set_f(R);
00065 }
00066
00067
00068 EvtComplex EvtPto3PAmp::amplitude(const EvtDalitzPoint& x) const
00069 {
00070 EvtComplex amp(1.0,0.0);
00071
00072 double m = sqrt(x.q(_pairRes));
00073 EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
00074 x.m(EvtCyclic3::second(_pairRes)),m);
00075 EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
00076
00077
00078
00079
00080 if(_typeN != NBW) {
00081
00082 _prop->set_g0(_g0*_vd.widthFactor(vd));
00083 }
00084
00085
00086
00087 amp *= _prop->evaluate(m);
00088
00089
00090
00091 amp *= _vd.formFactor(vd);
00092 amp *= _vb.formFactor(vb);
00093
00094 amp *= numerator(x);
00095
00096 return amp;
00097 }
00098
00099
00100 EvtComplex EvtPto3PAmp::numerator(const EvtDalitzPoint& x) const
00101 {
00102 EvtComplex ret(0.,0.);
00103 double m = sqrt(x.q(_pairRes));
00104 EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
00105 x.m(EvtCyclic3::second(_pairRes)),m);
00106 EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
00107
00108
00109
00110 if(NBW == _typeN) {
00111
00112 ret = angDep(x);
00113 }
00114
00115
00116
00117 else if(RBW_ZEMACH == _typeN) {
00118
00119 ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
00120 }
00121
00122
00123
00124 else if(RBW_KUEHN == _typeN) {
00125
00126 ret = _prop->m0()*_prop->m0() * angDep(x);
00127 }
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 if(RBW_CLEO == _typeN) {
00151
00152 Index iA = other(_pairAng);
00153 Index iB = common(_pairRes,_pairAng);
00154 Index iC = other(_pairRes);
00155
00156 double M = x.bigM();
00157 double mA = x.m(iA);
00158 double mB = x.m(iB);
00159 double mC = x.m(iC);
00160 double qAB = x.q(combine(iA,iB));
00161 double qBC = x.q(combine(iB,iC));
00162 double qCA = x.q(combine(iC,iA));
00163
00164 double m0 = _prop->m0();
00165
00166 if(_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
00167 else
00168 if(_spin == EvtSpinType::VECTOR) {
00169
00170 ret = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;;
00171 }
00172 else
00173 if(_spin == EvtSpinType::TENSOR) {
00174
00175 double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
00176 double x2 = M*M - mC*mC;
00177 double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/m0/m0;
00178 double x4 = mB*mB - mA*mA;
00179 double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/m0/m0;
00180 ret = (x1*x1 - 1./3.*x3*x5);
00181 }
00182 else assert(0);
00183 }
00184
00185 return ret;
00186 }
00187
00188
00189 double EvtPto3PAmp::angDep(const EvtDalitzPoint& x) const
00190 {
00191
00192
00193
00194 double cosTh = x.cosTh(_pairAng,_pairRes);
00195 if(fabs(cosTh) > 1.) {
00196
00197 report(INFO,"EvtGen") << "cosTh " << cosTh << endl;
00198 assert(0);
00199 }
00200
00201
00202
00203 return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
00204 }
00205