TEvtGen/EvtGenBase/EvtResonance.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtResonance.cc
  12 //
  13 // Description: resonance-defining class
  14 //
  15 // Modification history:
  16 //
  17 //    NK        September 4, 1997      Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <math.h>
  23 #include "EvtGenBase/EvtVector4R.hh"
  24 #include "EvtGenBase/EvtKine.hh"
  25 #include "EvtGenBase/EvtComplex.hh"
  26 #include "EvtGenBase/EvtResonance.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtConst.hh"
  29 using std::endl;
  30
  31 EvtResonance::~EvtResonance(){}
  32
  33
  34 EvtResonance& EvtResonance::operator = ( const EvtResonance  &n)
  35 {
  36   if ( &n == this ) return *this;
  37   _p4_p = n._p4_p;
  38   _p4_d1 = n._p4_d1;
  39   _p4_d2 = n._p4_d2;
  40   _ampl = n._ampl;
  41   _theta = n._theta;
  42   _gamma = n._gamma;
  43   _spin = n._spin;
  44   _bwm = n._bwm;
  45    return  *this;
  46 }
  47
  48
  49 EvtResonance::EvtResonance(const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
  50                            const  EvtVector4R& p4_d2, double ampl,
  51                            double theta, double gamma, double bwm, int spin):
  52   _p4_p(p4_p),_p4_d1(p4_d1), _p4_d2(p4_d2),_ampl(ampl), _theta(theta),
  53   _gamma(gamma), _bwm(bwm), _spin(spin) {}
  54
  55 EvtComplex EvtResonance::resAmpl() {
  56
  57   double pi180inv = 1.0/EvtConst::radToDegrees;
  58
  59   EvtComplex ampl;
  60   //EvtVector4R  _p4_d3 = _p4_p-_p4_d1-_p4_d2;
  61
  62   //get cos of the angle between the daughters from their 4-momenta
  63   //and the 4-momentum of the parent
  64
  65   //in general, EvtDecayAngle(parent, part1+part2, part1) gives the angle
  66   //the missing particle (not listed in the arguments) makes
  67   //with part2 in the rest frame of both
  68   //listed particles (12)
  69
  70   //angle 3 makes with 2 in rest frame of 12 (CS3)
  71   double cos_phi_0 = EvtDecayAngle(_p4_p, _p4_d1+_p4_d2, _p4_d1);
  72   //angle 3 makes with 1 in 12 is, of course, -cos_phi_0
  73
  74   switch (_spin) {
  75
  76   case 0 :
  77     ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
  78           sqrt(_gamma/EvtConst::twoPi)*
  79           (1.0/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,0.5*_gamma))));
  80     break;
  81
  82   case 1 :
  83     ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
  84           sqrt(_gamma/EvtConst::twoPi)*
  85           (cos_phi_0/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,0.5*_gamma))));
  86     break;
  87
  88   case 2:
  89     ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
  90           sqrt(_gamma/EvtConst::twoPi)*
  91           ((1.5*cos_phi_0*cos_phi_0-0.5)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0, 0.5*_gamma))));
  92     break;
  93
  94   case 3:
  95     ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
  96           sqrt(_gamma/EvtConst::twoPi)*
  97           ((2.5*cos_phi_0*cos_phi_0*cos_phi_0-1.5*cos_phi_0)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0, 0.5*_gamma))));
  98     break;
  99
 100   default:
 101     report(DEBUG,"EvtGen") << "EvtGen: wrong spin in EvtResonance" << endl;
 102     ampl = EvtComplex(0.0);
 103     break;
 104
 105   }
 106
 107   return ampl;
 108 }
 109
 110 EvtComplex EvtResonance::relBrWig(int i) {
 111
 112 //this function returns relativistic Breit-Wigner amplitude
 113 //for a given resonance (for P-wave decays of scalars only at the moment!)
 114
 115     EvtComplex BW;
 116     EvtVector4R  _p4_d3 = _p4_p-_p4_d1-_p4_d2;
 117     EvtVector4R _p4_12 = _p4_d1 + _p4_d2;
 118
 119     double msq13 = (_p4_d1 + _p4_d3).mass2();
 120     double msq23 = (_p4_d2 + _p4_d3).mass2();
 121     double msqParent = _p4_p.mass2();
 122     double msq1 = _p4_d1.mass2();
 123     double msq2 = _p4_d2.mass2();
 124     double msq3 = _p4_d3.mass2();
 125
 126     double M;
 127
 128     double p2 = sqrt((_p4_12.mass2() - (_p4_d1.mass() + _p4_d2.mass())*(_p4_d1.mass() + _p4_d2.mass()))*(_p4_12.mass2() - (_p4_d1.mass() - _p4_d2.mass())*(_p4_d1.mass() - _p4_d2.mass())))/(2.0*_p4_12.mass());
 129
 130     double p2R = sqrt((_bwm*_bwm - (_p4_d1.mass() + _p4_d2.mass())*(_p4_d1.mass() + _p4_d2.mass()))*(_bwm*_bwm - (_p4_d1.mass() - _p4_d2.mass())*(_p4_d1.mass() - _p4_d2.mass())))/(2.0*_bwm);
 131
 132     double gam, R;
 133
 134     if (i == 1) {
 135
 136         R = 2.0/(0.197);
 137
 138     }
 139     else R = 5.0/(0.197);
 140
 141     gam = _gamma*(_bwm/_p4_12.mass())*(p2/p2R)*(p2/p2R)*(p2/p2R)*((1 + R*R*p2R*p2R)/(1 + R*R*p2*p2));
 142     M = (msq13 - msq23 - (msqParent - msq3)*(msq1 - msq2)/(_bwm*_bwm))*sqrt((1 + R*R*p2R*p2R)/(1 + R*R*p2*p2));
 143
 144     BW = sqrt(_gamma)*M/((_bwm*_bwm - _p4_12.mass2()) - EvtComplex(0.0,1.0)*gam*_bwm);
 145
 146     return BW;
 147
 148 }
 149
 150