TEvtGen/EvtGenBase/EvtDalitzPlot.cxx

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //      $Id: EvtDalitzPlot.cc,v 1.23 2004/12/21 19:58:42 ryd 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   : _mA(0.), _mB(0.), _mC(0.), _bigM(0.),
  35     _ldel(0.), _rdel(0.)
  36 {}
  37
  38
  39 EvtDalitzPlot::EvtDalitzPlot(double mA, double mB, double mC, double bigM,
  40                              double ldel, double rdel)
  41   : _mA(mA), _mB(mB), _mC(mC), _bigM(bigM),
  42     _ldel(ldel), _rdel(rdel)
  43 {
  44   sanityCheck();
  45 }
  46
  47 EvtDalitzPlot::EvtDalitzPlot(const EvtDecayMode& mode, double ldel, double rdel )
  48 {
  49   _mA = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(A)));
  50   _mB = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(B)));
  51   _mC = EvtPDL::getMeanMass(EvtPDL::getId(mode.dau(C)));
  52   _bigM = EvtPDL::getMeanMass(EvtPDL::getId(mode.mother()));
  53
  54   _ldel = ldel;
  55   _rdel = rdel;
  56
  57   sanityCheck();
  58 }
  59
  60
  61 EvtDalitzPlot::EvtDalitzPlot(const EvtDalitzPlot& other)
  62   : _mA(other._mA), _mB(other._mB), _mC(other._mC), _bigM(other._bigM),
  63     _ldel(other._ldel), _rdel(other._rdel)
  64 {}
  65
  66
  67 EvtDalitzPlot::~EvtDalitzPlot()
  68 {}
  69
  70
  71 bool EvtDalitzPlot::operator==(const EvtDalitzPlot& other) const
  72 {
  73   bool ret = false;
  74   if(_mA == other._mA &&
  75      _mB == other._mB &&
  76      _mC == other._mC &&
  77      _bigM == other._bigM) ret = true;
  78
  79   return ret;
  80 }
  81
  82
  83 const EvtDalitzPlot* EvtDalitzPlot::clone() const
  84 {
  85   return new EvtDalitzPlot(*this);
  86 }
  87
  88
  89 void EvtDalitzPlot::sanityCheck() const
  90 {
  91   if(_mA < 0 || _mB < 0 || _mC < 0 || _bigM <= 0 || _bigM - _mA - _mB - _mC < 0.) {
  92
  93     printf("Invalid Dalitz plot %f %f %f %f\n",_mA,_mB,_mC,_bigM);
  94     assert(0);
  95   }
  96   assert(_ldel <= 0.);
  97   assert(_rdel >= 0.);
  98 }
  99
 100
 101 double EvtDalitzPlot::m(Index i) const {
 102
 103   double m = _mA;
 104   if(i == B) m = _mB;
 105   else
 106     if(i == C) m = _mC;
 107
 108   return m;
 109 }
 110
 111
 112 double EvtDalitzPlot::sum() const
 113 {
 114   return _mA*_mA + _mB*_mB + _mC*_mC + _bigM*_bigM;
 115 }
 116
 117
 118 double EvtDalitzPlot::qAbsMin(Pair i) const
 119 {
 120   Index j = first(i);
 121   Index k = second(i);
 122
 123   return (m(j) + m(k))*(m(j) + m(k));
 124 }
 125
 126
 127 double EvtDalitzPlot::qAbsMax(Pair i) const
 128 {
 129   Index j = other(i);
 130   return (_bigM-m(j))*(_bigM-m(j));
 131 }
 132
 133
 134 double EvtDalitzPlot::qResAbsMin(EvtCyclic3::Pair i) const
 135 {
 136   return qAbsMin(i) - sum()/3.;
 137 }
 138
 139 double EvtDalitzPlot::qResAbsMax(EvtCyclic3::Pair i) const
 140 {
 141   return  qAbsMax(i) - sum()/3.;
 142 }
 143
 144 double EvtDalitzPlot::qHelAbsMin(EvtCyclic3::Pair i) const
 145 {
 146   Pair j = next(i);
 147   Pair k = prev(i);
 148   return  (qAbsMin(j) - qAbsMax(k))/2.;
 149 }
 150
 151 double EvtDalitzPlot::qHelAbsMax(EvtCyclic3::Pair i) const
 152 {
 153   Pair j = next(i);
 154   Pair k = prev(i);
 155   return  (qAbsMax(j) - qAbsMin(k))/2.;
 156 }
 157
 158
 159 double EvtDalitzPlot::mAbsMin(Pair i) const
 160 {
 161   return sqrt(qAbsMin(i));
 162 }
 163
 164
 165 double EvtDalitzPlot::mAbsMax(Pair i) const
 166 {
 167   return sqrt(qAbsMax(i));
 168 }
 169
 170
 171 // parallel
 172
 173 double EvtDalitzPlot::qMin(Pair i, Pair j, double q) const
 174 {
 175   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     Index k0 = common(i,j);
 185     Index k2 = other(j);
 186     Index k1 = other(k0,k2);
 187
 188     // Energy, momentum of particle common to rest-frame and angle
 189     EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
 190     double pk = jpair.p();
 191     double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
 192
 193
 194     // Energy and momentum of the other particle
 195     EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
 196     double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
 197     double pj = mother.p(EvtTwoBodyKine::A);
 198
 199
 200     // See PDG 34.4.3.1
 201     return (ek+ej)*(ek+ej) - (pk+pj)*(pk+pj);
 202   }
 203 }
 204
 205
 206 // antiparallel
 207
 208 double EvtDalitzPlot::qMax(Pair i, Pair j, double q) const
 209 {
 210
 211   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     Index k0 = common(i,j);
 220     Index k2 = other(j);
 221     Index k1 = other(k0,k2);
 222
 223     // Energy, momentum of particle common to rest-frame and angle
 224     EvtTwoBodyKine jpair(m(k0),m(k1),sqrt(q));
 225     double ek = jpair.e(EvtTwoBodyKine::A,EvtTwoBodyKine::AB);
 226     double pk = jpair.p();
 227
 228     // Energy and momentum of the other particle
 229     EvtTwoBodyKine mother(sqrt(q),m(k2),bigM());
 230     double ej = mother.e(EvtTwoBodyKine::B,EvtTwoBodyKine::A);
 231     double pj = mother.p(EvtTwoBodyKine::A);
 232
 233
 234     // See PDG 34.4.3.1
 235     return (ek+ej)*(ek+ej) - (pk-pj)*(pk-pj);
 236   }
 237 }
 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   double dh = (qAbsMax(i) - qAbsMin(i))/((double) N);
 246   double sum = 0;
 247
 248   int ii;
 249   for(ii=1;ii<N;ii++) {
 250
 251     double x = qAbsMin(i) + ii*dh;
 252     double dy = qMax(j,i,x) - qMin(j,i,x);
 253     sum += dy;
 254   }
 255
 256   return sum * dh;
 257 }
 258
 259
 260 double EvtDalitzPlot::cosTh(EvtCyclic3::Pair i1, double q1, EvtCyclic3::Pair i2, double q2) const
 261 {
 262   if(i1 == i2) return 1.;
 263
 264   double qmax = qMax(i1,i2,q2);
 265   double qmin = qMin(i1,i2,q2);
 266
 267   double cos = (qmax + qmin - 2*q1)/(qmax - qmin);
 268
 269   return cos;
 270 }
 271
 272
 273 double EvtDalitzPlot::e(Index i, Pair j, double q) const
 274 {
 275   if(i == other(j)) {
 276
 277     // i does not belong to pair j
 278
 279     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     if(first(j) == i) k = second(j);
 287     else k = first(j);
 288
 289     double e = (q + m(i)*m(i) - m(k)*m(k))/2/sqrt(q);
 290     return e;
 291   }
 292 }
 293
 294
 295 double EvtDalitzPlot::p(Index i, Pair j, double q) const
 296 {
 297   double en = e(i,j,q);
 298   double p2 = en*en - m(i)*m(i);
 299
 300   if(p2 < 0) {
 301     printf("Bad value of p2 %f %d %d %f %f\n",p2,i,j,en,m(i));
 302     assert(0);
 303   }
 304
 305   return sqrt(p2);
 306 }
 307
 308
 309 double EvtDalitzPlot::q(EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2, double q2) const
 310 {
 311   if(i1 == i2) return q2;
 312
 313   EvtCyclic3::Index f = first(i1);
 314   EvtCyclic3::Index s = second(i1);
 315   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 }
 317
 318
 319 double EvtDalitzPlot::jacobian(EvtCyclic3::Pair i, double q) const
 320 {
 321   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   return EvtTwoBodyVertex(m(first(iRes)),
 328                           m(second(iRes)),m0,L);
 329 }
 330
 331
 332 EvtTwoBodyVertex EvtDalitzPlot::vB(Pair iRes, double m0, int L) const
 333 {
 334   return EvtTwoBodyVertex(m0,m(other(iRes)),bigM(),L);
 335 }
 336
 337
 338 void EvtDalitzPlot::print() const
 339 {
 340   // masses
 341   printf("Mass  M    %f\n",bigM());
 342   printf("Mass mA    %f\n",_mA);
 343   printf("Mass mB    %f\n",_mB);
 344   printf("Mass mC    %f\n",_mC);
 345   // limits
 346   printf("Limits qAB %f : %f\n",qAbsMin(AB),qAbsMax(AB));
 347   printf("Limits qBC %f : %f\n",qAbsMin(BC),qAbsMax(BC));
 348   printf("Limits qCA %f : %f\n",qAbsMin(CA),qAbsMax(CA));
 349   printf("Sum q       %f\n",sum());
 350   printf("Limit qsum  %f : %f\n",qSumMin(),qSumMax());
 351 }
 352
 353