00001 #include "EvtGenBase/EvtPatches.hh"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "EvtGenBase/EvtPatches.hh"
00012 #include <assert.h>
00013 #include <math.h>
00014 #include <stdio.h>
00015 #include "EvtGenBase/EvtDalitzPoint.hh"
00016 using namespace EvtCyclic3;
00017
00018 EvtDalitzPoint::EvtDalitzPoint()
00019 : _mA(-1.), _mB(-1.), _mC(-1.), _qAB(-1.), _qBC(-1.), _qCA(-1.)
00020 {}
00021
00022 EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC, double qAB, double qBC, double qCA)
00023 : _mA(mA), _mB(mB), _mC(mC), _qAB(qAB), _qBC(qBC), _qCA(qCA)
00024 {}
00025
00026
00027
00028 EvtDalitzPoint::EvtDalitzPoint(double mA, double mB, double mC,
00029 EvtCyclic3::Pair i,
00030 double qres, double qhel, double qsum)
00031 : _mA(mA), _mB(mB), _mC(mC)
00032 {
00033 double qi = qres + qsum/3.;
00034 double qj = -qres/2. + qhel + qsum/3.;
00035 double qk = -qres/2. - qhel + qsum/3.;
00036
00037 if(i == AB) { _qAB = qi; _qBC = qj; _qCA = qk; }
00038 else if(i == BC) { _qAB = qk; _qBC = qi; _qCA = qj; }
00039 else if(i == CA) { _qAB = qj; _qBC = qk; _qCA = qi; }
00040 }
00041
00042 EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPlot& dp, const EvtDalitzCoord& x)
00043 : _mA(dp.m(A)), _mB(dp.m(B)), _mC(dp.m(C))
00044 {
00045 if(x.pair1() == AB) _qAB = x.q1();
00046 else
00047 if(x.pair2() == AB) _qAB = x.q2();
00048 else _qAB = dp.sum() - x.q1() - x.q2();
00049
00050 if(x.pair1() == BC) _qBC = x.q1();
00051 else
00052 if(x.pair2() == BC) _qBC = x.q2();
00053 else _qBC = dp.sum() - x.q1() - x.q2();
00054
00055 if(x.pair1() == CA) _qCA = x.q1();
00056 else
00057 if(x.pair2() == CA) _qCA = x.q2();
00058 else _qCA = dp.sum() - x.q1() - x.q2();
00059
00060 }
00061
00062 EvtDalitzPoint::EvtDalitzPoint(const EvtDalitzPoint& other)
00063 : _mA(other._mA), _mB(other._mB), _mC(other._mC),
00064 _qAB(other._qAB), _qBC(other._qBC), _qCA(other._qCA)
00065 {}
00066
00067 EvtDalitzPoint::~EvtDalitzPoint()
00068 {}
00069
00070 double EvtDalitzPoint::q(EvtCyclic3::Pair i) const
00071 {
00072 double ret = _qAB;
00073 if(BC == i) ret = _qBC;
00074 else
00075 if(CA == i) ret = _qCA;
00076
00077 return ret;
00078 }
00079
00080 double EvtDalitzPoint::m(EvtCyclic3::Index i) const
00081 {
00082 double ret = _mA;
00083 if(B == i) ret = _mB;
00084 else
00085 if(C == i) ret = _mC;
00086
00087 return ret;
00088 }
00089
00090
00091
00092 double EvtDalitzPoint::qres(EvtCyclic3::Pair i) const
00093 {
00094 return (2.*q(i) - q(EvtCyclic3::prev(i)) - q(EvtCyclic3::next(i)))/3.;
00095 }
00096 double EvtDalitzPoint::qhel(EvtCyclic3::Pair i) const
00097 {
00098 Pair j = next(i);
00099 Pair k = prev(i);
00100 return (q(j) - q(k))/2.;
00101 }
00102 double EvtDalitzPoint::qsum() const
00103 {
00104 return _qAB + _qBC + _qCA;
00105 }
00106
00107
00108 double EvtDalitzPoint::qMin(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
00109 {
00110 EvtDalitzPlot dp = getDalitzPlot();
00111 return dp.qMin(i,j,q(j));
00112 }
00113
00114 double EvtDalitzPoint::qMax(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
00115 {
00116 EvtDalitzPlot dp = getDalitzPlot();
00117 return dp.qMax(i,j,q(j));
00118 }
00119
00120 double EvtDalitzPoint::pp(EvtCyclic3::Index i, EvtCyclic3::Index j) const
00121 {
00122 if(i == j) return m(i)*m(i);
00123 else return (q(combine(i,j)) - m(i)*m(i) - m(j)*m(j))/2.;
00124 }
00125
00126 double EvtDalitzPoint::e(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
00127 {
00128 EvtDalitzPlot dp = getDalitzPlot();
00129 return dp.e(i,j,q(j));
00130 }
00131
00132 double EvtDalitzPoint::p(EvtCyclic3::Index i, EvtCyclic3::Pair j) const
00133 {
00134 EvtDalitzPlot dp = getDalitzPlot();
00135 return dp.p(i,j,q(j));
00136 }
00137
00138 double EvtDalitzPoint::cosTh(EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes) const
00139 {
00140 EvtDalitzPlot dp = getDalitzPlot();
00141 return dp.cosTh(pairAng,q(pairAng),pairRes,q(pairRes));
00142 }
00143
00144
00145 EvtDalitzCoord EvtDalitzPoint::getDalitzPoint(EvtCyclic3::Pair i, EvtCyclic3::Pair j) const
00146 {
00147 return EvtDalitzCoord(i,q(i),j,q(j));
00148 }
00149
00150
00151 EvtDalitzPlot EvtDalitzPoint::getDalitzPlot() const
00152 {
00153 return EvtDalitzPlot(_mA,_mB,_mC,bigM());
00154 }
00155
00156 bool EvtDalitzPoint::isValid() const
00157 {
00158
00159
00160 double M = bigM();
00161 if(_mA < 0 || _mB < 0 || _mC < 0 || M <= 0) return false;
00162 if(M < _mA + _mB + _mC) return false;
00163
00164
00165
00166 bool inside = false;
00167 EvtDalitzPlot dp = getDalitzPlot();
00168
00169 if(dp.qAbsMin(AB) <= _qAB && _qAB <= dp.qAbsMax(AB))
00170 if(qMin(BC,AB) <= _qBC && _qBC <= qMax(BC,AB))
00171 inside = true;
00172
00173 return inside;
00174 }
00175
00176 double EvtDalitzPoint::bigM() const
00177 {
00178 return sqrt(_qAB+_qBC+_qCA - _mA*_mA - _mB*_mB - _mC*_mC);
00179 }
00180
00181
00182 void EvtDalitzPoint::print() const
00183 {
00184 getDalitzPlot().print();
00185 printf("%f %f %f\n",_qAB,_qBC,_qCA);
00186 }
00187
00188