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: EvtDalitzPoint.cpp,v 1.3 2009-03-16 15:53:27 robbep 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 0 : : _mA(-1.), _mB(-1.), _mC(-1.), _qAB(-1.), _qBC(-1.), _qCA(-1.)
20 0 : {}
21 :
22 : EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC, double qAB, double qBC, double qCA)
23 0 : : _mA(mA), _mB(mB), _mC(mC), _qAB(qAB), _qBC(qBC), _qCA(qCA)
24 0 : {}
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 0 : : _mA(mA), _mB(mB), _mC(mC)
32 0 : {
33 0 : double qi = qres + qsum/3.;
34 0 : double qj = -qres/2. + qhel + qsum/3.;
35 0 : double qk = -qres/2. - qhel + qsum/3.;
36 :
37 0 : if(i == AB) { _qAB = qi; _qBC = qj; _qCA = qk; }
38 0 : else if(i == BC) { _qAB = qk; _qBC = qi; _qCA = qj; }
39 0 : else if(i == CA) { _qAB = qj; _qBC = qk; _qCA = qi; }
40 0 : }
41 :
42 : EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPlot& dp, const EvtDalitzCoord& x)
43 0 : : _mA(dp.m(A)), _mB(dp.m(B)), _mC(dp.m(C))
44 0 : {
45 0 : if(x.pair1() == AB) _qAB = x.q1();
46 : else
47 0 : if(x.pair2() == AB) _qAB = x.q2();
48 0 : else _qAB = dp.sum() - x.q1() - x.q2();
49 :
50 0 : if(x.pair1() == BC) _qBC = x.q1();
51 : else
52 0 : if(x.pair2() == BC) _qBC = x.q2();
53 0 : else _qBC = dp.sum() - x.q1() - x.q2();
54 :
55 0 : if(x.pair1() == CA) _qCA = x.q1();
56 : else
57 0 : if(x.pair2() == CA) _qCA = x.q2();
58 0 : else _qCA = dp.sum() - x.q1() - x.q2();
59 :
60 0 : }
61 :
62 : EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPoint& other)
63 0 : : _mA(other._mA), _mB(other._mB), _mC(other._mC),
64 0 : _qAB(other._qAB), _qBC(other._qBC), _qCA(other._qCA)
65 0 : {}
66 :
67 : EvtDalitzPoint::~EvtDalitzPoint()
68 0 : {}
69 :
70 : double EvtDalitzPoint::q(EvtCyclic3::Pair i) const
71 : {
72 0 : double ret = _qAB;
73 0 : if(BC == i) ret = _qBC;
74 : else
75 0 : if(CA == i) ret = _qCA;
76 :
77 0 : return ret;
78 : }
79 :
80 : double EvtDalitzPoint::m(EvtCyclic3::Index i) const
81 : {
82 0 : double ret = _mA;
83 0 : if(B == i) ret = _mB;
84 : else
85 0 : if(C == i) ret = _mC;
86 :
87 0 : return ret;
88 : }
89 :
90 : // Zemach variables
91 :
92 : double EvtDalitzPoint::qres(EvtCyclic3::Pair i) const
93 : {
94 0 : return (2.*q(i) - q(EvtCyclic3::prev(i)) - q(EvtCyclic3::next(i)))/3.;
95 : }
96 : double EvtDalitzPoint::qhel(EvtCyclic3::Pair i) const
97 : {
98 0 : Pair j = next(i);
99 0 : Pair k = prev(i);
100 0 : return (q(j) - q(k))/2.;
101 : }
102 : double EvtDalitzPoint::qsum() const
103 : {
104 0 : return _qAB + _qBC + _qCA;
105 : }
106 :
107 :
108 : double EvtDalitzPoint::qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
109 : {
110 0 : EvtDalitzPlot dp = getDalitzPlot();
111 0 : return dp.qMin(i,j,q(j));
112 0 : }
113 :
114 : double EvtDalitzPoint::qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
115 : {
116 0 : EvtDalitzPlot dp = getDalitzPlot();
117 0 : return dp.qMax(i,j,q(j));
118 0 : }
119 :
120 : double EvtDalitzPoint::pp(EvtCyclic3::Index i, EvtCyclic3::Index j) const
121 : {
122 0 : if(i == j) return m(i)*m(i);
123 0 : else return (q(combine(i,j)) - m(i)*m(i) - m(j)*m(j))/2.;
124 0 : }
125 :
126 : double EvtDalitzPoint::e(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
127 : {
128 0 : EvtDalitzPlot dp = getDalitzPlot();
129 0 : return dp.e(i,j,q(j));
130 0 : }
131 :
132 : double EvtDalitzPoint::p(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
133 : {
134 0 : EvtDalitzPlot dp = getDalitzPlot();
135 0 : return dp.p(i,j,q(j));
136 0 : }
137 :
138 : double EvtDalitzPoint::cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const
139 : {
140 0 : EvtDalitzPlot dp = getDalitzPlot();
141 0 : return dp.cosTh(pairAng,q(pairAng),pairRes,q(pairRes));
142 0 : }
143 :
144 :
145 : EvtDalitzCoord EvtDalitzPoint::getDalitzPoint(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
146 : {
147 0 : return EvtDalitzCoord(i,q(i),j,q(j));
148 : }
149 :
150 :
151 : EvtDalitzPlot EvtDalitzPoint::getDalitzPlot() const
152 : {
153 0 : return EvtDalitzPlot(_mA,_mB,_mC,bigM());
154 : }
155 :
156 : bool EvtDalitzPoint::isValid() const
157 : {
158 : // Check masses
159 :
160 0 : double M = bigM();
161 0 : if(_mA < 0 || _mB < 0 || _mC < 0 || M <= 0) return false;
162 0 : if(M < _mA + _mB + _mC) return false;
163 :
164 : // Check that first coordinate is within absolute limits
165 :
166 : bool inside = false;
167 0 : EvtDalitzPlot dp = getDalitzPlot();
168 :
169 0 : if(dp.qAbsMin(AB) <= _qAB && _qAB <= dp.qAbsMax(AB))
170 0 : if(qMin(BC,AB) <= _qBC && _qBC <= qMax(BC,AB))
171 0 : inside = true;
172 :
173 0 : return inside;
174 0 : };
175 :
176 : double EvtDalitzPoint::bigM() const
177 : {
178 0 : return sqrt(_qAB+_qBC+_qCA - _mA*_mA - _mB*_mB - _mC*_mC);
179 : }
180 :
181 :
182 : void EvtDalitzPoint::print() const
183 : {
184 0 : getDalitzPlot().print();
185 0 : printf("%f %f %f\n",_qAB,_qBC,_qCA);
186 0 : }
187 :
188 :
|
