1 #include "EvtGenBase/EvtPatches.hh"
2 /*******************************************************************************
3 * Project: BaBar detector at the SLAC PEP-II B-factory
5 * File: $Id: EvtPto3PAmp.cc,v 1.21 2009/02/18 03:31:38 ryd Exp $
6 * Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
8 * Copyright (C) 2002 Caltech
10 * Modified: May 30, 2005, Denis Dujmic (ddujmic@slac.stanford.edu): flip sign
11 * for vector resonances (-i -> i at resonance mass).
12 * Introduce Flatte lineshape.
13 *******************************************************************************/
18 #include "EvtGenBase/EvtComplex.hh"
19 #include "EvtGenBase/EvtPto3PAmp.hh"
20 #include "EvtGenBase/EvtDalitzCoord.hh"
21 #include "EvtGenBase/EvtdFunction.hh"
22 #include "EvtGenBase/EvtCyclic3.hh"
23 #include "EvtGenBase/EvtReport.hh"
26 using EvtCyclic3::Index;
27 using EvtCyclic3::Pair;
30 EvtPto3PAmp::EvtPto3PAmp(EvtDalitzPlot dp, Pair pairAng, Pair pairRes,
31 EvtSpinType::spintype spin,
32 const EvtPropagator& prop, NumType typeN)
33 : EvtAmplitude<EvtDalitzPoint>(),
34 _pairAng(pairAng), _pairRes(pairRes),
37 _prop((EvtPropagator*) prop.clone()),
41 _vb(prop.m0(),dp.m(EvtCyclic3::other(pairRes)),dp.bigM(),spin),
42 _vd(dp.m(EvtCyclic3::first(pairRes)),dp.m(EvtCyclic3::second(pairRes)),prop.m0(),spin)
47 EvtPto3PAmp::EvtPto3PAmp(const EvtPto3PAmp& other)
48 : EvtAmplitude<EvtDalitzPoint>(other),
49 _pairAng(other._pairAng),
50 _pairRes(other._pairRes),
53 _prop( (other._prop) ? (EvtPropagator*) other._prop->clone() : 0),
57 _vb(other._vb), _vd(other._vd)
61 EvtPto3PAmp::~EvtPto3PAmp()
63 if(_prop) delete _prop;
67 void EvtPto3PAmp::set_fd(double R)
72 void EvtPto3PAmp::set_fb(double R)
78 EvtComplex EvtPto3PAmp::amplitude(const EvtDalitzPoint& x) const
80 EvtComplex amp(1.0,0.0);
82 double m = sqrt(x.q(_pairRes));
84 if ((_max>0 && m > _max) || (_min>0 && m < _min)) return EvtComplex(0.0,0.0);
86 EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
87 x.m(EvtCyclic3::second(_pairRes)),m);
88 EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
91 // Compute mass-dependent width for relativistic propagators
93 if(_typeN!=NBW && _typeN!=FLATTE) {
95 _prop->set_g0(_g0*_vd.widthFactor(vd));
101 amp *= evalPropagator(m);
105 // Compute form-factors
107 amp *= _vd.formFactor(vd);
108 amp *= _vb.formFactor(vb);
116 EvtComplex EvtPto3PAmp::numerator(const EvtDalitzPoint& x) const
118 EvtComplex ret(0.,0.);
119 double m = sqrt(x.q(_pairRes));
120 EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
121 x.m(EvtCyclic3::second(_pairRes)),m);
122 EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
124 // Non-relativistic Breit-Wigner
131 // Standard relativistic Zemach propagator
133 else if(RBW_ZEMACH == _typeN) {
135 ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
138 // Kuehn-Santamaria normalization:
140 else if(RBW_KUEHN == _typeN) {
142 ret = _prop->m0()*_prop->m0() * angDep(x);
146 // CLEO amplitude is not factorizable
148 // The CLEO amplitude numerator is proportional to:
150 // m2_AC - m2_BC + (m2_D - m2_C)(m2_B - m2_A)/m2_0
152 // m2_AC = (eA + eC)^2 + (P - P_C cosTh(BC))^2
153 // m2_BC = (eB + eC)^2 + (P + P_C cosTh(BC))^2
155 // The first term m2_AB-m2_BC is therefore a p-wave term
157 // The second term is an s-wave, the amplitude
158 // does not factorize!
160 // The first term is just Zemach. However, the sign is flipped!
161 // Let's consistently use the convention in which the amplitude
162 // is proportional to +cosTh(BC). In the CLEO expressions, I will
163 // therefore exchange AB to get rid of the sign flip.
166 if(RBW_CLEO == _typeN || FLATTE == _typeN || GS == _typeN) {
168 Index iA = EvtCyclic3::other(_pairAng); // A = other(BC)
169 Index iB = EvtCyclic3::common(_pairRes,_pairAng); // B = common(AB,BC)
170 Index iC = EvtCyclic3::other(_pairRes); // C = other(AB)
176 double qAB = x.q(EvtCyclic3::combine(iA,iB));
177 double qBC = x.q(EvtCyclic3::combine(iB,iC));
178 double qCA = x.q(EvtCyclic3::combine(iC,iA));
180 //double m0 = _prop->m0();
182 if(_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
184 if(_spin == EvtSpinType::VECTOR) {
186 //ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
187 ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;
190 if(_spin == EvtSpinType::TENSOR) {
192 //double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
193 double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;
194 double x2 = M*M - mC*mC;
195 //double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/m0/m0;
196 double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/qAB;
197 double x4 = mB*mB - mA*mA;
198 //double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/m0/m0;
199 double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/qAB;
200 ret = (x1*x1 - 1./3.*x3*x5);
209 double EvtPto3PAmp::angDep(const EvtDalitzPoint& x) const
211 // Angular dependece for factorizable amplitudes
212 // unphysical cosines indicate we are in big trouble
214 double cosTh = x.cosTh(_pairAng,_pairRes);
215 if(fabs(cosTh) > 1.) {
217 report(INFO,"EvtGen") << "cosTh " << cosTh << endl;
221 // in units of half-spin
223 return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));