Line data Source code
1 : //-----------------------------------------------------------------------
2 : // File and Version Information:
3 : // $Id: EvtDalitzPlot.cpp,v 1.3 2009-03-16 15:53:27 robbep Exp $
4 : //
5 : // Environment:
6 : // This software is part of the EvtGen package developed jointly
7 : // for the BaBar and CLEO collaborations. If you use all or part
8 : // of it, please give an appropriate acknowledgement.
9 : //
10 : // Copyright Information:
11 : // Copyright (C) 1998 Caltech, UCSB
12 : //
13 : // Module creator:
14 : // Alexei Dvoretskii, Caltech, 2001-2002.
15 : //-----------------------------------------------------------------------
16 : #include "EvtGenBase/EvtPatches.hh"
17 :
18 : // Global 3-body Dalitz decay kinematics as defined by the mass
19 : // of the mother and the daughters. Spins are not considered.
20 :
21 : #include <math.h>
22 : #include <assert.h>
23 : #include <stdio.h>
24 : #include "EvtGenBase/EvtPatches.hh"
25 : #include "EvtGenBase/EvtDalitzPlot.hh"
26 : #include "EvtGenBase/EvtTwoBodyVertex.hh"
27 : #include "EvtGenBase/EvtPDL.hh"
28 : #include "EvtGenBase/EvtDecayMode.hh"
29 :
30 : using namespace EvtCyclic3;
31 :
32 :
33 : EvtDalitzPlot::EvtDalitzPlot()
34 0 : : _mA(0.), _mB(0.), _mC(0.), _bigM(0.),
35 0 : _ldel(0.), _rdel(0.)
36 0 : {}
37 :
38 :
39 : EvtDalitzPlot::EvtDalitzPlot(double mA, double mB, double mC, double bigM,
40 : double ldel, double rdel)
41 0 : : _mA(mA), _mB(mB), _mC(mC), _bigM(bigM),
42 0 : _ldel(ldel), _rdel(rdel)
43 0 : {
44 0 : sanityCheck();
45 0 : }
46 :
47 : EvtDalitzPlot::EvtDalitzPlot(const EvtDecayMode& mode, double ldel, double rdel )
48 0 : {
49 0 : _mA = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(A)));
50 0 : _mB = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(B)));
51 0 : _mC = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(C)));
52 0 : _bigM = EvtPDL::getMeanMass(EvtPDL::getId(mode.mother()));
53 :
54 0 : _ldel = ldel;
55 0 : _rdel = rdel;
56 :
57 0 : sanityCheck();
58 0 : }
59 :
60 :
61 : EvtDalitzPlot::EvtDalitzPlot(const EvtDalitzPlot& other)
62 0 : : _mA(other._mA), _mB(other._mB), _mC(other._mC), _bigM(other._bigM),
63 0 : _ldel(other._ldel), _rdel(other._rdel)
64 0 : {}
65 :
66 :
67 : EvtDalitzPlot::~EvtDalitzPlot()
68 0 : {}
69 :
70 :
71 : bool EvtDalitzPlot::operator==(const EvtDalitzPlot& other) const
72 : {
73 : bool ret = false;
74 0 : if(_mA == other._mA &&
75 0 : _mB == other._mB &&
76 0 : _mC == other._mC &&
77 0 : _bigM == other._bigM) ret = true;
78 :
79 0 : return ret;
80 : }
81 :
82 :
83 : const EvtDalitzPlot* EvtDalitzPlot::clone() const
84 : {
85 0 : return new EvtDalitzPlot(*this);
86 : }
87 :
88 :
89 : void EvtDalitzPlot::sanityCheck() const
90 : {
91 0 : if(_mA < 0 || _mB < 0 || _mC < 0 || _bigM <= 0 || _bigM - _mA - _mB - _mC < 0.) {
92 :
93 0 : printf("Invalid Dalitz plot %f %f %f %f\n",_mA,_mB,_mC,_bigM);
94 0 : assert(0);
95 : }
96 0 : assert(_ldel <= 0.);
97 0 : assert(_rdel >= 0.);
98 0 : }
99 :
100 :
101 : double EvtDalitzPlot::m(Index i) const {
102 :
103 0 : double m = _mA;
104 0 : if(i == B) m = _mB;
105 : else
106 0 : if(i == C) m = _mC;
107 :
108 0 : return m;
109 : }
110 :
111 :
112 : double EvtDalitzPlot::sum() const
113 : {
114 0 : return _mA*_mA + _mB*_mB + _mC*_mC + _bigM*_bigM;
115 : }
116 :
117 :
118 : double EvtDalitzPlot::qAbsMin(Pair i) const
119 : {
120 0 : Index j = first(i);
121 0 : Index k = second(i);
122 :
123 0 : return (m(j) + m(k))*(m(j) + m(k));
124 : }
125 :
126 :
127 : double EvtDalitzPlot::qAbsMax(Pair i) const
128 : {
129 0 : Index j = other(i);
130 0 : return (_bigM-m(j))*(_bigM-m(j));
131 : }
132 :
133 :
134 : double EvtDalitzPlot::qResAbsMin(EvtCyclic3::Pair i) const
135 : {
136 0 : return qAbsMin(i) - sum()/3.;
137 : }
138 :
139 : double EvtDalitzPlot::qResAbsMax(EvtCyclic3::Pair i) const
140 : {
141 0 : return qAbsMax(i) - sum()/3.;
142 : }
143 :
144 : double EvtDalitzPlot::qHelAbsMin(EvtCyclic3::Pair i) const
145 : {
146 0 : Pair j = next(i);
147 0 : Pair k = prev(i);
148 0 : return (qAbsMin(j) - qAbsMax(k))/2.;
149 : }
150 :
151 : double EvtDalitzPlot::qHelAbsMax(EvtCyclic3::Pair i) const
152 : {
153 0 : Pair j = next(i);
154 0 : Pair k = prev(i);
155 0 : return (qAbsMax(j) - qAbsMin(k))/2.;
156 : }
157 :
158 :
159 : double EvtDalitzPlot::mAbsMin(Pair i) const
160 : {
161 0 : return sqrt(qAbsMin(i));
162 : }
163 :
164 :
165 : double EvtDalitzPlot::mAbsMax(Pair i) const
166 : {
167 0 : return sqrt(qAbsMax(i));
168 : }
169 :
170 :
171 : // parallel
172 :
173 : double EvtDalitzPlot::qMin(Pair i, Pair j, double q) const
174 : {
175 0 : if(i == j) return q;
176 :
177 : else {
178 :
179 : // Particle pair j defines the rest-frame
180 : // 0 - particle common to r.f. and angle calculations
181 : // 1 - particle belonging to r.f. but not angle
182 : // 2 - particle not belonging to r.f.
183 :
184 0 : Index k0 = common(i,j);
185 0 : Index k2 = other(j);
186 0 : Index k1 = other(k0,k2);
187 :
188 : // Energy, momentum of particle common to rest-frame and angle
189 0 : EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
190 0 : double pk = jpair.p();
191 0 : double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
192 :
193 :
194 : // Energy and momentum of the other particle
195 0 : EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
196 0 : double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
197 0 : double pj = mother.p(EvtTwoBodyKine::A);
198 :
199 :
200 : // See PDG 34.4.3.1
201 0 : return (ek+ej)*(ek+ej) - (pk+pj)*(pk+pj);
202 0 : }
203 0 : }
204 :
205 :
206 : // antiparallel
207 :
208 : double EvtDalitzPlot::qMax(Pair i, Pair j, double q) const
209 : {
210 :
211 0 : if(i == j) return q;
212 : else {
213 :
214 : // Particle pair j defines the rest-frame
215 : // 0 - particle common to r.f. and angle calculations
216 : // 1 - particle belonging to r.f. but not angle
217 : // 2 - particle not belonging to r.f.
218 :
219 0 : Index k0 = common(i,j);
220 0 : Index k2 = other(j);
221 0 : Index k1 = other(k0,k2);
222 :
223 : // Energy, momentum of particle common to rest-frame and angle
224 0 : EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
225 0 : double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
226 0 : double pk = jpair.p();
227 :
228 : // Energy and momentum of the other particle
229 0 : EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
230 0 : double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
231 0 : double pj = mother.p(EvtTwoBodyKine::A);
232 :
233 :
234 : // See PDG 34.4.3.1
235 0 : return (ek+ej)*(ek+ej) - (pk-pj)*(pk-pj);
236 0 : }
237 0 : }
238 :
239 :
240 : double EvtDalitzPlot::getArea(int N, Pair i, Pair j) const
241 : {
242 : // Trapezoidal integral over qi. qj can be calculated.
243 : // The first and the last point are zero, so they are not counted
244 :
245 0 : double dh = (qAbsMax(i) - qAbsMin(i))/((double) N);
246 : double sum = 0;
247 :
248 : int ii;
249 0 : for(ii=1;ii<N;ii++) {
250 :
251 0 : double x = qAbsMin(i) + ii*dh;
252 0 : double dy = qMax(j,i,x) - qMin(j,i,x);
253 0 : sum += dy;
254 : }
255 :
256 0 : return sum * dh;
257 : }
258 :
259 :
260 : double EvtDalitzPlot::cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
261 : {
262 0 : if(i1 == i2) return 1.;
263 :
264 0 : double qmax = qMax(i1,i2,q2);
265 0 : double qmin = qMin(i1,i2,q2);
266 :
267 0 : double cos = (qmax + qmin - 2*q1)/(qmax - qmin);
268 :
269 : return cos;
270 0 : }
271 :
272 :
273 : double EvtDalitzPlot::e(Index i, Pair j, double q) const
274 : {
275 0 : if(i == other(j)) {
276 :
277 : // i does not belong to pair j
278 :
279 0 : return (bigM()*bigM()-q-m(i)*m(i))/2/sqrt(q);
280 : }
281 : else {
282 :
283 : // i and k make pair j
284 :
285 : Index k;
286 0 : if(first(j) == i) k = second(j);
287 0 : else k = first(j);
288 :
289 0 : double e = (q + m(i)*m(i) - m(k)*m(k))/2/sqrt(q);
290 : return e;
291 : }
292 0 : }
293 :
294 :
295 : double EvtDalitzPlot::p(Index i, Pair j, double q) const
296 : {
297 0 : double en = e(i,j,q);
298 0 : double p2 = en*en - m(i)*m(i);
299 :
300 0 : if(p2 < 0) {
301 0 : printf("Bad value of p2 %f %d %d %f %f\n",p2,i,j,en,m(i));
302 0 : assert(0);
303 : }
304 :
305 0 : return sqrt(p2);
306 : }
307 :
308 :
309 : double EvtDalitzPlot::q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const
310 : {
311 0 : if(i1 == i2) return q2;
312 :
313 0 : EvtCyclic3::Index f = first(i1);
314 0 : EvtCyclic3::Index s = second(i1);
315 0 : return m(f)*m(f) + m(s)*m(s) + 2*e(f,i2,q2)*e(s,i2,q2) - 2*p(f,i2,q2)*p(s,i2,q2)*cosTh;
316 0 : }
317 :
318 :
319 : double EvtDalitzPlot::jacobian(EvtCyclic3::Pair i, double q) const
320 : {
321 0 : return 2*p(first(i),i,q)*p(other(i),i,q); // J(BC) = 2pA*pB = 2pA*pC
322 : }
323 :
324 :
325 : EvtTwoBodyVertex EvtDalitzPlot::vD(Pair iRes, double m0, int L) const
326 : {
327 0 : return EvtTwoBodyVertex(m(first(iRes)),
328 0 : m(second(iRes)),m0,L);
329 : }
330 :
331 :
332 : EvtTwoBodyVertex EvtDalitzPlot::vB(Pair iRes, double m0, int L) const
333 : {
334 0 : return EvtTwoBodyVertex(m0,m(other(iRes)),bigM(),L);
335 : }
336 :
337 :
338 : void EvtDalitzPlot::print() const
339 : {
340 : // masses
341 0 : printf("Mass M %f\n",bigM());
342 0 : printf("Mass mA %f\n",_mA);
343 0 : printf("Mass mB %f\n",_mB);
344 0 : printf("Mass mC %f\n",_mC);
345 : // limits
346 0 : printf("Limits qAB %f : %f\n",qAbsMin(AB),qAbsMax(AB));
347 0 : printf("Limits qBC %f : %f\n",qAbsMin(BC),qAbsMax(BC));
348 0 : printf("Limits qCA %f : %f\n",qAbsMin(CA),qAbsMax(CA));
349 0 : printf("Sum q %f\n",sum());
350 0 : printf("Limit qsum %f : %f\n",qSumMin(),qSumMax());
351 0 : }
352 :
353 :
|
