TEvtGen/EvtGenBase/EvtPto3PAmp.cxx

   1 #include "EvtGenBase/EvtPatches.hh"
   2 /*******************************************************************************
   3  * Project: BaBar detector at the SLAC PEP-II B-factory
   4  * Package: EvtGenBase
   5  *    File: $Id: EvtPto3PAmp.cc,v 1.21 2009/02/18 03:31:38 ryd Exp $
   6  *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
   7  *
   8  * Copyright (C) 2002 Caltech
   9  *
  10  * Modified: May 30, 2005, Denis Dujmic (ddujmic@slac.stanford.edu): flip sign
  11  *           for vector resonances (-i -> i at resonance mass).
  12  *           Introduce Flatte lineshape.
  13  *******************************************************************************/
  14
  15 #include <assert.h>
  16 #include <math.h>
  17 #include <iostream>
  18 #include "EvtGenBase/EvtComplex.hh"
  19 #include "EvtGenBase/EvtPto3PAmp.hh"
  20 #include "EvtGenBase/EvtDalitzCoord.hh"
  21 #include "EvtGenBase/EvtdFunction.hh"
  22 #include "EvtGenBase/EvtCyclic3.hh"
  23 #include "EvtGenBase/EvtReport.hh"
  24
  25 using std::endl;
  26 using EvtCyclic3::Index;
  27 using EvtCyclic3::Pair;
  28
  29
  30 EvtPto3PAmp::EvtPto3PAmp(EvtDalitzPlot dp, Pair pairAng, Pair pairRes,
  31                          EvtSpinType::spintype spin,
  32                          const EvtPropagator& prop, NumType typeN)
  33   : EvtAmplitude<EvtDalitzPoint>(),
  34   _pairAng(pairAng), _pairRes(pairRes),
  35   _spin(spin),
  36   _typeN(typeN),
  37   _prop((EvtPropagator*) prop.clone()),
  38   _g0(prop.g0()),
  39   _min(0),
  40   _max(0),
  41   _vb(prop.m0(),dp.m(EvtCyclic3::other(pairRes)),dp.bigM(),spin),
  42   _vd(dp.m(EvtCyclic3::first(pairRes)),dp.m(EvtCyclic3::second(pairRes)),prop.m0(),spin)
  43 {}
  44
  45
  46
  47 EvtPto3PAmp::EvtPto3PAmp(const EvtPto3PAmp& other)
  48   : EvtAmplitude<EvtDalitzPoint>(other),
  49   _pairAng(other._pairAng),
  50   _pairRes(other._pairRes),
  51   _spin(other._spin),
  52   _typeN(other._typeN),
  53   _prop( (other._prop) ? (EvtPropagator*) other._prop->clone() : 0),
  54   _g0(other._g0),
  55   _min(other._min),
  56   _max(other._max),
  57   _vb(other._vb), _vd(other._vd)
  58 {}
  59
  60
  61 EvtPto3PAmp::~EvtPto3PAmp()
  62 {
  63   if(_prop) delete _prop;
  64 }
  65
  66
  67 void EvtPto3PAmp::set_fd(double R)
  68 {
  69   _vd.set_f(R);
  70 }
  71
  72 void EvtPto3PAmp::set_fb(double R)
  73 {
  74   _vb.set_f(R);
  75 }
  76
  77
  78 EvtComplex EvtPto3PAmp::amplitude(const EvtDalitzPoint& x) const
  79 {
  80   EvtComplex amp(1.0,0.0);
  81
  82   double m = sqrt(x.q(_pairRes));
  83
  84   if ((_max>0 && m > _max) || (_min>0 && m < _min))  return EvtComplex(0.0,0.0);
  85
  86   EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
  87                     x.m(EvtCyclic3::second(_pairRes)),m);
  88   EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
  89
  90
  91   // Compute mass-dependent width for relativistic propagators
  92
  93   if(_typeN!=NBW && _typeN!=FLATTE) {
  94
  95     _prop->set_g0(_g0*_vd.widthFactor(vd));
  96   }
  97
  98
  99   // Compute propagator
 100
 101   amp *= evalPropagator(m);
 102
 103
 104
 105   // Compute form-factors
 106
 107   amp *= _vd.formFactor(vd);
 108   amp *= _vb.formFactor(vb);
 109
 110   amp *= numerator(x);
 111
 112   return amp;
 113 }
 114
 115
 116 EvtComplex EvtPto3PAmp::numerator(const EvtDalitzPoint& x) const
 117 {
 118   EvtComplex ret(0.,0.);
 119   double m = sqrt(x.q(_pairRes));
 120   EvtTwoBodyKine vd(x.m(EvtCyclic3::first(_pairRes)),
 121                     x.m(EvtCyclic3::second(_pairRes)),m);
 122   EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
 123
 124   // Non-relativistic Breit-Wigner
 125
 126   if(NBW == _typeN) {
 127
 128     ret = angDep(x);
 129   }
 130
 131   // Standard relativistic Zemach propagator
 132
 133   else if(RBW_ZEMACH == _typeN) {
 134
 135     ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
 136   }
 137
 138   // Kuehn-Santamaria normalization:
 139
 140   else if(RBW_KUEHN == _typeN) {
 141
 142     ret = _prop->m0()*_prop->m0() * angDep(x);
 143   }
 144
 145
 146   // CLEO amplitude is not factorizable
 147   //
 148   // The CLEO amplitude numerator is proportional to:
 149   //
 150   // m2_AC - m2_BC + (m2_D - m2_C)(m2_B - m2_A)/m2_0
 151   //
 152   // m2_AC = (eA + eC)^2 + (P - P_C cosTh(BC))^2
 153   // m2_BC = (eB + eC)^2 + (P + P_C cosTh(BC))^2
 154   //
 155   // The first term m2_AB-m2_BC is therefore a p-wave term
 156   // - 4PP_C cosTh(BC)
 157   // The second term is an s-wave, the amplitude
 158   // does not factorize!
 159   //
 160   // The first term is just Zemach. However, the sign is flipped!
 161   // Let's consistently use the convention in which the amplitude
 162   // is proportional to +cosTh(BC). In the CLEO expressions, I will
 163   // therefore exchange AB to get rid of the sign flip.
 164
 165
 166   if(RBW_CLEO == _typeN || FLATTE == _typeN || GS == _typeN) {
 167
 168     Index iA = EvtCyclic3::other(_pairAng);           // A = other(BC)
 169     Index iB = EvtCyclic3::common(_pairRes,_pairAng); // B = common(AB,BC)
 170     Index iC = EvtCyclic3::other(_pairRes);           // C = other(AB)
 171
 172     double M = x.bigM();
 173     double mA = x.m(iA);
 174     double mB = x.m(iB);
 175     double mC = x.m(iC);
 176     double qAB = x.q(EvtCyclic3::combine(iA,iB));
 177     double qBC = x.q(EvtCyclic3::combine(iB,iC));
 178     double qCA = x.q(EvtCyclic3::combine(iC,iA));
 179
 180     //double m0 = _prop->m0();
 181
 182     if(_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
 183     else
 184       if(_spin == EvtSpinType::VECTOR) {
 185
 186         //ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
 187         ret = qCA - qBC - (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;
 188       }
 189       else
 190         if(_spin == EvtSpinType::TENSOR) {
 191
 192           //double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/m0/m0;
 193           double x1 = qBC - qCA + (M*M - mC*mC)*(mA*mA - mB*mB)/qAB;
 194           double x2 = M*M - mC*mC;
 195           //double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/m0/m0;
 196           double x3 = qAB - 2*M*M - 2*mC*mC + x2*x2/qAB;
 197           double x4 = mB*mB - mA*mA;
 198           //double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/m0/m0;
 199           double x5 = qAB - 2*mB*mB - 2*mA*mA + x4*x4/qAB;
 200           ret = (x1*x1 - 1./3.*x3*x5);
 201         }
 202         else assert(0);
 203   }
 204
 205   return ret;
 206 }
 207
 208
 209 double EvtPto3PAmp::angDep(const EvtDalitzPoint& x) const
 210 {
 211   // Angular dependece for factorizable amplitudes
 212   // unphysical cosines indicate we are in big trouble
 213
 214   double cosTh = x.cosTh(_pairAng,_pairRes);
 215   if(fabs(cosTh) > 1.) {
 216
 217     report(INFO,"EvtGen") << "cosTh " << cosTh << endl;
 218     assert(0);
 219   }
 220
 221   // in units of half-spin
 222
 223   return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
 224 }
 225