]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | #include "EvtGenBase/EvtPatches.hh" |
2 | /******************************************************************************* | |
3 | * Project: BaBar detector at the SLAC PEP-II B-factory | |
4 | * Package: EvtGenBase | |
5 | * File: $Id: EvtDalitzPoint.cc,v 1.14 2004/12/21 19:58:42 ryd Exp $ | |
6 | * Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002 | |
7 | * | |
8 | * Copyright (C) 2002 Caltech | |
9 | *******************************************************************************/ | |
10 | ||
11 | #include "EvtGenBase/EvtPatches.hh" | |
12 | #include <assert.h> | |
13 | #include <math.h> | |
14 | #include <stdio.h> | |
15 | #include "EvtGenBase/EvtDalitzPoint.hh" | |
16 | using namespace EvtCyclic3; | |
17 | ||
18 | EvtDalitzPoint::EvtDalitzPoint() | |
19 | : _mA(-1.), _mB(-1.), _mC(-1.), _qAB(-1.), _qBC(-1.), _qCA(-1.) | |
20 | {} | |
21 | ||
22 | EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC, double qAB, double qBC, double qCA) | |
23 | : _mA(mA), _mB(mB), _mC(mC), _qAB(qAB), _qBC(qBC), _qCA(qCA) | |
24 | {} | |
25 | ||
26 | // Constructor from Zemach coordinates | |
27 | ||
28 | EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC, | |
29 | EvtCyclic3::Pair i, | |
30 | double qres, double qhel, double qsum) | |
31 | : _mA(mA), _mB(mB), _mC(mC) | |
32 | { | |
33 | double qi = qres + qsum/3.; | |
34 | double qj = -qres/2. + qhel + qsum/3.; | |
35 | double qk = -qres/2. - qhel + qsum/3.; | |
36 | ||
37 | if(i == AB) { _qAB = qi; _qBC = qj; _qCA = qk; } | |
38 | else if(i == BC) { _qAB = qk; _qBC = qi; _qCA = qj; } | |
39 | else if(i == CA) { _qAB = qj; _qBC = qk; _qCA = qi; } | |
40 | } | |
41 | ||
42 | EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPlot& dp, const EvtDalitzCoord& x) | |
43 | : _mA(dp.m(A)), _mB(dp.m(B)), _mC(dp.m(C)) | |
44 | { | |
45 | if(x.pair1() == AB) _qAB = x.q1(); | |
46 | else | |
47 | if(x.pair2() == AB) _qAB = x.q2(); | |
48 | else _qAB = dp.sum() - x.q1() - x.q2(); | |
49 | ||
50 | if(x.pair1() == BC) _qBC = x.q1(); | |
51 | else | |
52 | if(x.pair2() == BC) _qBC = x.q2(); | |
53 | else _qBC = dp.sum() - x.q1() - x.q2(); | |
54 | ||
55 | if(x.pair1() == CA) _qCA = x.q1(); | |
56 | else | |
57 | if(x.pair2() == CA) _qCA = x.q2(); | |
58 | else _qCA = dp.sum() - x.q1() - x.q2(); | |
59 | ||
60 | } | |
61 | ||
62 | EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPoint& other) | |
63 | : _mA(other._mA), _mB(other._mB), _mC(other._mC), | |
64 | _qAB(other._qAB), _qBC(other._qBC), _qCA(other._qCA) | |
65 | {} | |
66 | ||
67 | EvtDalitzPoint::~EvtDalitzPoint() | |
68 | {} | |
69 | ||
70 | double EvtDalitzPoint::q(EvtCyclic3::Pair i) const | |
71 | { | |
72 | double ret = _qAB; | |
73 | if(BC == i) ret = _qBC; | |
74 | else | |
75 | if(CA == i) ret = _qCA; | |
76 | ||
77 | return ret; | |
78 | } | |
79 | ||
80 | double EvtDalitzPoint::m(EvtCyclic3::Index i) const | |
81 | { | |
82 | double ret = _mA; | |
83 | if(B == i) ret = _mB; | |
84 | else | |
85 | if(C == i) ret = _mC; | |
86 | ||
87 | return ret; | |
88 | } | |
89 | ||
90 | // Zemach variables | |
91 | ||
92 | double EvtDalitzPoint::qres(EvtCyclic3::Pair i) const | |
93 | { | |
94 | return (2.*q(i) - q(EvtCyclic3::prev(i)) - q(EvtCyclic3::next(i)))/3.; | |
95 | } | |
96 | double EvtDalitzPoint::qhel(EvtCyclic3::Pair i) const | |
97 | { | |
98 | Pair j = next(i); | |
99 | Pair k = prev(i); | |
100 | return (q(j) - q(k))/2.; | |
101 | } | |
102 | double EvtDalitzPoint::qsum() const | |
103 | { | |
104 | return _qAB + _qBC + _qCA; | |
105 | } | |
106 | ||
107 | ||
108 | double EvtDalitzPoint::qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const | |
109 | { | |
110 | EvtDalitzPlot dp = getDalitzPlot(); | |
111 | return dp.qMin(i,j,q(j)); | |
112 | } | |
113 | ||
114 | double EvtDalitzPoint::qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const | |
115 | { | |
116 | EvtDalitzPlot dp = getDalitzPlot(); | |
117 | return dp.qMax(i,j,q(j)); | |
118 | } | |
119 | ||
120 | double EvtDalitzPoint::pp(EvtCyclic3::Index i, EvtCyclic3::Index j) const | |
121 | { | |
122 | if(i == j) return m(i)*m(i); | |
123 | else return (q(combine(i,j)) - m(i)*m(i) - m(j)*m(j))/2.; | |
124 | } | |
125 | ||
126 | double EvtDalitzPoint::e(EvtCyclic3::Index i, EvtCyclic3::Pair j) const | |
127 | { | |
128 | EvtDalitzPlot dp = getDalitzPlot(); | |
129 | return dp.e(i,j,q(j)); | |
130 | } | |
131 | ||
132 | double EvtDalitzPoint::p(EvtCyclic3::Index i, EvtCyclic3::Pair j) const | |
133 | { | |
134 | EvtDalitzPlot dp = getDalitzPlot(); | |
135 | return dp.p(i,j,q(j)); | |
136 | } | |
137 | ||
138 | double EvtDalitzPoint::cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const | |
139 | { | |
140 | EvtDalitzPlot dp = getDalitzPlot(); | |
141 | return dp.cosTh(pairAng,q(pairAng),pairRes,q(pairRes)); | |
142 | } | |
143 | ||
144 | ||
145 | EvtDalitzCoord EvtDalitzPoint::getDalitzPoint(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const | |
146 | { | |
147 | return EvtDalitzCoord(i,q(i),j,q(j)); | |
148 | } | |
149 | ||
150 | ||
151 | EvtDalitzPlot EvtDalitzPoint::getDalitzPlot() const | |
152 | { | |
153 | return EvtDalitzPlot(_mA,_mB,_mC,bigM()); | |
154 | } | |
155 | ||
156 | bool EvtDalitzPoint::isValid() const | |
157 | { | |
158 | // Check masses | |
159 | ||
160 | double M = bigM(); | |
161 | if(_mA < 0 || _mB < 0 || _mC < 0 || M <= 0) return false; | |
162 | if(M < _mA + _mB + _mC) return false; | |
163 | ||
164 | // Check that first coordinate is within absolute limits | |
165 | ||
166 | bool inside = false; | |
167 | EvtDalitzPlot dp = getDalitzPlot(); | |
168 | ||
169 | if(dp.qAbsMin(AB) <= _qAB && _qAB <= dp.qAbsMax(AB)) | |
170 | if(qMin(BC,AB) <= _qBC && _qBC <= qMax(BC,AB)) | |
171 | inside = true; | |
172 | ||
173 | return inside; | |
174 | } | |
175 | ||
176 | double EvtDalitzPoint::bigM() const | |
177 | { | |
178 | return sqrt(_qAB+_qBC+_qCA - _mA*_mA - _mB*_mB - _mC*_mC); | |
179 | } | |
180 | ||
181 | ||
182 | void EvtDalitzPoint::print() const | |
183 | { | |
184 | getDalitzPlot().print(); | |
185 | printf("%f %f %f\n",_qAB,_qBC,_qCA); | |
186 | } | |
187 | ||
188 |