00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "EvtGenBase/EvtPatches.hh"
00017
00018
00019
00020
00021 #include <math.h>
00022 #include <assert.h>
00023 #include <stdio.h>
00024 #include "EvtGenBase/EvtPatches.hh"
00025 #include "EvtGenBase/EvtDalitzPlot.hh"
00026 #include "EvtGenBase/EvtTwoBodyVertex.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenBase/EvtDecayMode.hh"
00029
00030 using namespace EvtCyclic3;
00031
00032
00033 EvtDalitzPlot::EvtDalitzPlot()
00034 : _mA(0.), _mB(0.), _mC(0.), _bigM(0.),
00035 _ldel(0.), _rdel(0.)
00036 {}
00037
00038
00039 EvtDalitzPlot::EvtDalitzPlot(double mA, double mB, double mC, double bigM,
00040 double ldel, double rdel)
00041 : _mA(mA), _mB(mB), _mC(mC), _bigM(bigM),
00042 _ldel(ldel), _rdel(rdel)
00043 {
00044 sanityCheck();
00045 }
00046
00047 EvtDalitzPlot::EvtDalitzPlot(const EvtDecayMode& mode, double ldel, double rdel )
00048 {
00049 _mA = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(A)));
00050 _mB = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(B)));
00051 _mC = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(C)));
00052 _bigM = EvtPDL::getMeanMass(EvtPDL::getId(mode.mother()));
00053
00054 _ldel = ldel;
00055 _rdel = rdel;
00056
00057 sanityCheck();
00058 }
00059
00060
00061 EvtDalitzPlot::EvtDalitzPlot(const EvtDalitzPlot& other)
00062 : _mA(other._mA), _mB(other._mB), _mC(other._mC), _bigM(other._bigM),
00063 _ldel(other._ldel), _rdel(other._rdel)
00064 {}
00065
00066
00067 EvtDalitzPlot::~EvtDalitzPlot()
00068 {}
00069
00070
00071 bool EvtDalitzPlot::operator==(const EvtDalitzPlot& other) const
00072 {
00073 bool ret = false;
00074 if(_mA == other._mA &&
00075 _mB == other._mB &&
00076 _mC == other._mC &&
00077 _bigM == other._bigM) ret = true;
00078
00079 return ret;
00080 }
00081
00082
00083 const EvtDalitzPlot* EvtDalitzPlot::clone() const
00084 {
00085 return new EvtDalitzPlot(*this);
00086 }
00087
00088
00089 void EvtDalitzPlot::sanityCheck() const
00090 {
00091 if(_mA < 0 || _mB < 0 || _mC < 0 || _bigM <= 0 || _bigM - _mA - _mB - _mC < 0.) {
00092
00093 printf("Invalid Dalitz plot %f %f %f %f\n",_mA,_mB,_mC,_bigM);
00094 assert(0);
00095 }
00096 assert(_ldel <= 0.);
00097 assert(_rdel >= 0.);
00098 }
00099
00100
00101 double EvtDalitzPlot::m(Index i) const {
00102
00103 double m = _mA;
00104 if(i == B) m = _mB;
00105 else
00106 if(i == C) m = _mC;
00107
00108 return m;
00109 }
00110
00111
00112 double EvtDalitzPlot::sum() const
00113 {
00114 return _mA*_mA + _mB*_mB + _mC*_mC + _bigM*_bigM;
00115 }
00116
00117
00118 double EvtDalitzPlot::qAbsMin(Pair i) const
00119 {
00120 Index j = first(i);
00121 Index k = second(i);
00122
00123 return (m(j) + m(k))*(m(j) + m(k));
00124 }
00125
00126
00127 double EvtDalitzPlot::qAbsMax(Pair i) const
00128 {
00129 Index j = other(i);
00130 return (_bigM-m(j))*(_bigM-m(j));
00131 }
00132
00133
00134 double EvtDalitzPlot::qResAbsMin(EvtCyclic3::Pair i) const
00135 {
00136 return qAbsMin(i) - sum()/3.;
00137 }
00138
00139 double EvtDalitzPlot::qResAbsMax(EvtCyclic3::Pair i) const
00140 {
00141 return qAbsMax(i) - sum()/3.;
00142 }
00143
00144 double EvtDalitzPlot::qHelAbsMin(EvtCyclic3::Pair i) const
00145 {
00146 Pair j = next(i);
00147 Pair k = prev(i);
00148 return (qAbsMin(j) - qAbsMax(k))/2.;
00149 }
00150
00151 double EvtDalitzPlot::qHelAbsMax(EvtCyclic3::Pair i) const
00152 {
00153 Pair j = next(i);
00154 Pair k = prev(i);
00155 return (qAbsMax(j) - qAbsMin(k))/2.;
00156 }
00157
00158
00159 double EvtDalitzPlot::mAbsMin(Pair i) const
00160 {
00161 return sqrt(qAbsMin(i));
00162 }
00163
00164
00165 double EvtDalitzPlot::mAbsMax(Pair i) const
00166 {
00167 return sqrt(qAbsMax(i));
00168 }
00169
00170
00171
00172
00173 double EvtDalitzPlot::qMin(Pair i, Pair j, double q) const
00174 {
00175 if(i == j) return q;
00176
00177 else {
00178
00179
00180
00181
00182
00183
00184 Index k0 = common(i,j);
00185 Index k2 = other(j);
00186 Index k1 = other(k0,k2);
00187
00188
00189 EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
00190 double pk = jpair.p();
00191 double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
00192
00193
00194
00195 EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
00196 double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
00197 double pj = mother.p(EvtTwoBodyKine::A);
00198
00199
00200
00201 return (ek+ej)*(ek+ej) - (pk+pj)*(pk+pj);
00202 }
00203 }
00204
00205
00206
00207
00208 double EvtDalitzPlot::qMax(Pair i, Pair j, double q) const
00209 {
00210
00211 if(i == j) return q;
00212 else {
00213
00214
00215
00216
00217
00218
00219 Index k0 = common(i,j);
00220 Index k2 = other(j);
00221 Index k1 = other(k0,k2);
00222
00223
00224 EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
00225 double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
00226 double pk = jpair.p();
00227
00228
00229 EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
00230 double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
00231 double pj = mother.p(EvtTwoBodyKine::A);
00232
00233
00234
00235 return (ek+ej)*(ek+ej) - (pk-pj)*(pk-pj);
00236 }
00237 }
00238
00239
00240 double EvtDalitzPlot::getArea(int N, Pair i, Pair j) const
00241 {
00242
00243
00244
00245 double dh = (qAbsMax(i) - qAbsMin(i))/((double) N);
00246 double sum = 0;
00247
00248 int ii;
00249 for(ii=1;ii<N;ii++) {
00250
00251 double x = qAbsMin(i) + ii*dh;
00252 double dy = qMax(j,i,x) - qMin(j,i,x);
00253 sum += dy;
00254 }
00255
00256 return sum * dh;
00257 }
00258
00259
00260 double EvtDalitzPlot::cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
00261 {
00262 if(i1 == i2) return 1.;
00263
00264 double qmax = qMax(i1,i2,q2);
00265 double qmin = qMin(i1,i2,q2);
00266
00267 double cos = (qmax + qmin - 2*q1)/(qmax - qmin);
00268
00269 return cos;
00270 }
00271
00272
00273 double EvtDalitzPlot::e(Index i, Pair j, double q) const
00274 {
00275 if(i == other(j)) {
00276
00277
00278
00279 return (bigM()*bigM()-q-m(i)*m(i))/2/sqrt(q);
00280 }
00281 else {
00282
00283
00284
00285 Index k;
00286 if(first(j) == i) k = second(j);
00287 else k = first(j);
00288
00289 double e = (q + m(i)*m(i) - m(k)*m(k))/2/sqrt(q);
00290 return e;
00291 }
00292 }
00293
00294
00295 double EvtDalitzPlot::p(Index i, Pair j, double q) const
00296 {
00297 double en = e(i,j,q);
00298 double p2 = en*en - m(i)*m(i);
00299
00300 if(p2 < 0) {
00301 printf("Bad value of p2 %f %d %d %f %f\n",p2,i,j,en,m(i));
00302 assert(0);
00303 }
00304
00305 return sqrt(p2);
00306 }
00307
00308
00309 double EvtDalitzPlot::q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const
00310 {
00311 if(i1 == i2) return q2;
00312
00313 EvtCyclic3::Index f = first(i1);
00314 EvtCyclic3::Index s = second(i1);
00315 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;
00316 }
00317
00318
00319 double EvtDalitzPlot::jacobian(EvtCyclic3::Pair i, double q) const
00320 {
00321 return 2*p(first(i),i,q)*p(other(i),i,q);
00322 }
00323
00324
00325 EvtTwoBodyVertex EvtDalitzPlot::vD(Pair iRes, double m0, int L) const
00326 {
00327 return EvtTwoBodyVertex(m(first(iRes)),
00328 m(second(iRes)),m0,L);
00329 }
00330
00331
00332 EvtTwoBodyVertex EvtDalitzPlot::vB(Pair iRes, double m0, int L) const
00333 {
00334 return EvtTwoBodyVertex(m0,m(other(iRes)),bigM(),L);
00335 }
00336
00337
00338 void EvtDalitzPlot::print() const
00339 {
00340
00341 printf("Mass M %f\n",bigM());
00342 printf("Mass mA %f\n",_mA);
00343 printf("Mass mB %f\n",_mB);
00344 printf("Mass mC %f\n",_mC);
00345
00346 printf("Limits qAB %f : %f\n",qAbsMin(AB),qAbsMax(AB));
00347 printf("Limits qBC %f : %f\n",qAbsMin(BC),qAbsMax(BC));
00348 printf("Limits qCA %f : %f\n",qAbsMin(CA),qAbsMax(CA));
00349 printf("Sum q %f\n",sum());
00350 printf("Limit qsum %f : %f\n",qSumMin(),qSumMax());
00351 }
00352
00353