TEvtGen/EvtGenBase/EvtPropFlatte.cxx

   1 #include "EvtGenBase/EvtPatches.hh"
   2 /*******************************************************************************
   3  * Project: BaBar detector at the SLAC PEP-II B-factory
   4  * Package: EvtGenBase
   5  * Author : D. Dujmic, J. Thompson
   6  *
   7  * Copyright (C) 2005 SLAC
   8  *******************************************************************************/
   9
  10 #include <math.h>
  11 #include "EvtGenBase/EvtPropFlatte.hh"
  12
  13 #include <iostream>
  14 using std::cout;
  15 using std::endl;
  16
  17 EvtPropFlatte::EvtPropFlatte(double m0,
  18                              double g0, double m0a, double m0b,
  19                              double g1, double m1a, double m1b) :
  20   EvtPropagator( m0, g0),
  21   _m0a(m0a),
  22   _m0b(m0b),
  23   _g1 (g1),
  24   _m1a(m1a),
  25   _m1b(m1b)
  26 {}
  27
  28
  29 EvtPropFlatte::EvtPropFlatte(const EvtPropFlatte& other) :
  30   EvtPropagator(other),
  31   _m0a (other._m0a),
  32   _m0b (other._m0b),
  33   _g1  (other._g1),
  34   _m1a (other._m1a),
  35   _m1b (other._m1b)
  36 {}
  37
  38
  39 EvtPropFlatte::~EvtPropFlatte()
  40 {}
  41
  42
  43 EvtAmplitude<EvtPoint1D>* EvtPropFlatte::clone() const
  44 {
  45   return new EvtPropFlatte(*this);
  46 }
  47
  48
  49
  50 EvtComplex EvtPropFlatte::amplitude(const EvtPoint1D& x) const
  51 {
  52
  53   /*
  54
  55   Use BES parameterization:
  56
  57                         1.
  58       -----------------------------------------
  59        m0^2 - m^2 - i*m0*( g1*rho1 + g2*rho2 )
  60
  61
  62   Resonance mass: m0
  63   Channel1: m0a, m0b, g0
  64   Channel2: m1a, m1b, g1
  65
  66   where breakup momenta q's are:
  67
  68           E0a = (m^2 + m0a^2 - m0b^2) / 2m
  69           q0  = sqrt( E0a^2 - m0a^2 )
  70
  71           E1a = (m^2 + m1a^2 - m1b^2) / 2m
  72           q1  = sqrt( E1a^2 - m1a^2 )
  73
  74
  75   */
  76
  77
  78
  79   double s = x.value()*x.value();
  80   double m = x.value();
  81
  82   double E0a  = 0.5 * (s + _m0a*_m0a - _m0b*_m0b) / m;
  83   double qSq0 = E0a*E0a - _m0a*_m0a;
  84
  85   double E1a  = 0.5 * (s + _m1a*_m1a - _m1b*_m1b) / m;
  86   double qSq1 = E1a*E1a - _m1a*_m1a;
  87
  88   EvtComplex gamma0 = qSq0 >= 0 ?  EvtComplex(  _g0 * sqrt(qSq0), 0)  : EvtComplex( 0, _g0 * sqrt(-qSq0) );
  89   EvtComplex gamma1 = qSq1 >= 0 ?  EvtComplex(  _g1 * sqrt(qSq1), 0)  : EvtComplex( 0, _g1 * sqrt(-qSq1) );
  90
  91   EvtComplex gamma = gamma0 + gamma1;
  92
  93   EvtComplex a = 1.0/( _m0*_m0 - s - EvtComplex(0.0,2*_m0/m)*gamma  );
  94
  95   return a;
  96 }
  97