TEvtGen/EvtGenBase/EvtDalitzPoint.cxx

   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