Line data Source code
1 : #include "EvtGenBase/EvtPatches.hh"
2 : /*******************************************************************************
3 : * Project: BaBar detector at the SLAC PEP-II B-factory
4 : * Package: EvtGenBase
5 : * File: $Id: EvtPto3PAmp.cpp,v 1.3 2009-03-16 15:44:04 robbep Exp $
6 : * Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
7 : *
8 : * Copyright (C) 2002 Caltech
9 : *
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 : *******************************************************************************/
14 :
15 : #include <assert.h>
16 : #include <math.h>
17 : #include <iostream>
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"
24 :
25 : using std::endl;
26 : using EvtCyclic3::Index;
27 : using EvtCyclic3::Pair;
28 :
29 :
30 : EvtPto3PAmp::EvtPto3PAmp(EvtDalitzPlot dp, Pair pairAng, Pair pairRes,
31 : EvtSpinType::spintype spin,
32 : const EvtPropagator& prop, NumType typeN)
33 0 : : EvtAmplitude<EvtDalitzPoint>(),
34 0 : _pairAng(pairAng), _pairRes(pairRes),
35 0 : _spin(spin),
36 0 : _typeN(typeN),
37 0 : _prop((EvtPropagator*) prop.clone()),
38 0 : _g0(prop.g0()),
39 0 : _min(0),
40 0 : _max(0),
41 0 : _vb(prop.m0(),dp.m(EvtCyclic3::other(pairRes)),dp.bigM(),spin),
42 0 : _vd(dp.m(EvtCyclic3::first(pairRes)),dp.m(EvtCyclic3::second(pairRes)),prop.m0(),spin)
43 0 : {}
44 :
45 :
46 :
47 : EvtPto3PAmp::EvtPto3PAmp(const EvtPto3PAmp& other)
48 0 : : EvtAmplitude<EvtDalitzPoint>(other),
49 0 : _pairAng(other._pairAng),
50 0 : _pairRes(other._pairRes),
51 0 : _spin(other._spin),
52 0 : _typeN(other._typeN),
53 0 : _prop( (other._prop) ? (EvtPropagator*) other._prop->clone() : 0),
54 0 : _g0(other._g0),
55 0 : _min(other._min),
56 0 : _max(other._max),
57 0 : _vb(other._vb), _vd(other._vd)
58 0 : {}
59 :
60 :
61 : EvtPto3PAmp::~EvtPto3PAmp()
62 0 : {
63 0 : if(_prop) delete _prop;
64 0 : }
65 :
66 :
67 : void EvtPto3PAmp::set_fd(double R)
68 : {
69 0 : _vd.set_f(R);
70 0 : }
71 :
72 : void EvtPto3PAmp::set_fb(double R)
73 : {
74 0 : _vb.set_f(R);
75 0 : }
76 :
77 :
78 : EvtComplex EvtPto3PAmp::amplitude(const EvtDalitzPoint& x) const
79 : {
80 0 : EvtComplex amp(1.0,0.0);
81 :
82 0 : double m = sqrt(x.q(_pairRes));
83 :
84 0 : if ((_max>0 && m > _max) || (_min>0 && m < _min)) return EvtComplex(0.0,0.0);
85 :
86 0 : EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
87 0 : x.m(EvtCyclic3::second(_pairRes)),m);
88 0 : EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
89 :
90 :
91 : // Compute mass-dependent width for relativistic propagators
92 :
93 0 : if(_typeN!=NBW && _typeN!=FLATTE) {
94 :
95 0 : _prop->set_g0(_g0*_vd.widthFactor(vd));
96 0 : }
97 :
98 :
99 : // Compute propagator
100 :
101 0 : amp *= evalPropagator(m);
102 :
103 :
104 :
105 : // Compute form-factors
106 :
107 0 : amp *= _vd.formFactor(vd);
108 0 : amp *= _vb.formFactor(vb);
109 :
110 0 : amp *= numerator(x);
111 :
112 0 : return amp;
113 0 : }
114 :
115 :
116 : EvtComplex EvtPto3PAmp::numerator(const EvtDalitzPoint& x) const
117 : {
118 0 : EvtComplex ret(0.,0.);
119 0 : double m = sqrt(x.q(_pairRes));
120 0 : EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
121 0 : x.m(EvtCyclic3::second(_pairRes)),m);
122 0 : EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
123 :
124 : // Non-relativistic Breit-Wigner
125 :
126 0 : if(NBW == _typeN) {
127 :
128 0 : ret = angDep(x);
129 0 : }
130 :
131 : // Standard relativistic Zemach propagator
132 :
133 0 : else if(RBW_ZEMACH == _typeN) {
134 :
135 0 : ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
136 0 : }
137 :
138 : // Kuehn-Santamaria normalization:
139 :
140 0 : else if(RBW_KUEHN == _typeN) {
141 :
142 0 : ret = _prop->m0()*_prop->m0() * angDep(x);
143 0 : }
144 :
145 :
146 : // CLEO amplitude is not factorizable
147 : //
148 : // The CLEO amplitude numerator is proportional to:
149 : //
150 : // m2_AC - m2_BC + (m2_D - m2_C)(m2_B - m2_A)/m2_0
151 : //
152 : // m2_AC = (eA + eC)^2 + (P - P_C cosTh(BC))^2
153 : // m2_BC = (eB + eC)^2 + (P + P_C cosTh(BC))^2
154 : //
155 : // The first term m2_AB-m2_BC is therefore a p-wave term
156 : // - 4PP_C cosTh(BC)
157 : // The second term is an s-wave, the amplitude
158 : // does not factorize!
159 : //
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.
164 :
165 :
166 0 : if(RBW_CLEO == _typeN || FLATTE == _typeN || GS == _typeN) {
167 :
168 0 : Index iA = EvtCyclic3::other(_pairAng); // A = other(BC)
169 0 : Index iB = EvtCyclic3::common(_pairRes,_pairAng); // B = common(AB,BC)
170 0 : Index iC = EvtCyclic3::other(_pairRes); // C = other(AB)
171 :
172 0 : double M = x.bigM();
173 0 : double mA = x.m(iA);
174 0 : double mB = x.m(iB);
175 0 : double mC = x.m(iC);
176 0 : double qAB = x.q(EvtCyclic3::combine(iA,iB));
177 0 : double qBC = x.q(EvtCyclic3::combine(iB,iC));
178 0 : double qCA = x.q(EvtCyclic3::combine(iC,iA));
179 :
180 : //double m0 = _prop->m0();
181 :
182 0 : if(_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
183 : else
184 0 : if(_spin == EvtSpinType::VECTOR) {
185 :
186 : //ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
187 0 : ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;
188 0 : }
189 : else
190 0 : if(_spin == EvtSpinType::TENSOR) {
191 :
192 : //double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
193 0 : 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 0 : double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/qAB;
197 0 : double x4 = mB*mB - mA*mA;
198 : //double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/m0/m0;
199 0 : double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/qAB;
200 0 : ret = (x1*x1 - 1./3.*x3*x5);
201 : }
202 0 : else assert(0);
203 0 : }
204 :
205 : return ret;
206 0 : }
207 :
208 :
209 : double EvtPto3PAmp::angDep(const EvtDalitzPoint& x) const
210 : {
211 : // Angular dependece for factorizable amplitudes
212 : // unphysical cosines indicate we are in big trouble
213 :
214 0 : double cosTh = x.cosTh(_pairAng,_pairRes);
215 0 : if(fabs(cosTh) > 1.) {
216 :
217 0 : report(Severity::Info,"EvtGen") << "cosTh " << cosTh << endl;
218 0 : assert(0);
219 : }
220 :
221 : // in units of half-spin
222 :
223 0 : return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
224 : }
225 :
|
