EvtResonance.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtResonance.cc
//
// Description: resonance-defining class
//
// Modification history:
//
//    NK        September 4, 1997      Module created
//
//------------------------------------------------------------------------
//
#include "EvtResonance.hh"
#include "EvtComplex.hh"
#include "EvtConst.hh"
#include "EvtKine.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include "EvtVector4R.hh"
#include <math.h>
using std::endl;
 
EvtResonance::~EvtResonance() {}
 
// operator
 
EvtResonance& EvtResonance::operator=( const EvtResonance& n ) {
  if ( &n == this ) return *this;
  _p4_p  = n._p4_p;
  _p4_d1 = n._p4_d1;
  _p4_d2 = n._p4_d2;
  _ampl  = n._ampl;
  _theta = n._theta;
  _gamma = n._gamma;
  _spin  = n._spin;
  _bwm   = n._bwm;
  return *this;
}
 
// constructor
 
EvtResonance::EvtResonance( const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
                            const EvtVector4R& p4_d2, double ampl, double theta, double gamma,
                            double bwm, int spin )
    : _p4_p( p4_p )
    , _p4_d1( p4_d1 )
    , _p4_d2( p4_d2 )
    , _ampl( ampl )
    , _theta( theta )
    , _gamma( gamma )
    , _bwm( bwm )
    , _spin( spin ) {}
 
// amplitude function
 
EvtComplex EvtResonance::resAmpl() {
 
  double pi180inv = 1.0 / EvtConst::radToDegrees;
 
  EvtComplex ampl;
  // EvtVector4R  _p4_d3 = _p4_p-_p4_d1-_p4_d2;
 
  // get cos of the angle between the daughters from their 4-momenta
  // and the 4-momentum of the parent
 
  // in general, EvtDecayAngle(parent, part1+part2, part1) gives the angle
  // the missing particle (not listed in the arguments) makes
  // with part2 in the rest frame of both
  // listed particles (12)
 
  // angle 3 makes with 2 in rest frame of 12 (CS3)
  double cos_phi_0 = EvtDecayAngle( _p4_p, _p4_d1 + _p4_d2, _p4_d1 );
  // angle 3 makes with 1 in 12 is, of course, -cos_phi_0
 
  switch ( _spin )
  {
 
  case 0:
    ampl =
        ( _ampl * EvtComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) *
          sqrt( _gamma / EvtConst::twoPi ) *
          ( 1.0 / ( ( _p4_d1 + _p4_d2 ).mass() - _bwm - EvtComplex( 0.0, 0.5 * _gamma ) ) ) );
    break;
 
  case 1:
    ampl = ( _ampl * EvtComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) *
             sqrt( _gamma / EvtConst::twoPi ) *
             ( cos_phi_0 /
               ( ( _p4_d1 + _p4_d2 ).mass() - _bwm - EvtComplex( 0.0, 0.5 * _gamma ) ) ) );
    break;
 
  case 2:
    ampl = ( _ampl * EvtComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) *
             sqrt( _gamma / EvtConst::twoPi ) *
             ( ( 1.5 * cos_phi_0 * cos_phi_0 - 0.5 ) /
               ( ( _p4_d1 + _p4_d2 ).mass() - _bwm - EvtComplex( 0.0, 0.5 * _gamma ) ) ) );
    break;
 
  case 3:
    ampl = ( _ampl * EvtComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) *
             sqrt( _gamma / EvtConst::twoPi ) *
             ( ( 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 ) ) ) );
    break;
 
  default:
    report( DEBUG, "EvtGen" ) << "EvtGen: wrong spin in EvtResonance" << endl;
    ampl = EvtComplex( 0.0 );
    break;
  }
 
  return ampl;
}
 
EvtComplex EvtResonance::relBrWig( int i ) {
 
  // this function returns relativistic Breit-Wigner amplitude
  // for a given resonance (for P-wave decays of scalars only at the moment!)
 
  EvtComplex  BW;
  EvtVector4R _p4_d3 = _p4_p - _p4_d1 - _p4_d2;
  EvtVector4R _p4_12 = _p4_d1 + _p4_d2;
 
  // EvtVector4R _p4_pm3 = _p4_p - _p4_d3;
  // EvtVector4R _p4_1m2 = _p4_d1 - _p4_d2;
 
  // double msq12 = _p4_12.mass2();
  double msq13     = ( _p4_d1 + _p4_d3 ).mass2();
  double msq23     = ( _p4_d2 + _p4_d3 ).mass2();
  double msqParent = _p4_p.mass2();
  double msq1      = _p4_d1.mass2();
  double msq2      = _p4_d2.mass2();
  double msq3      = _p4_d3.mass2();
 
  double M;
 
  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() );
 
  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 );
 
  double gam, R;
 
  if ( i == 1 )
  {
    // consider this the K resonance
 
    R = 2.0 / ( 0.197 );
  }
  else R = 5.0 / ( 0.197 );
 
  gam = _gamma * ( _bwm / _p4_12.mass() ) * ( p2 / p2R ) * ( p2 / p2R ) * ( p2 / p2R ) *
        ( ( 1 + R * R * p2R * p2R ) / ( 1 + R * R * p2 * p2 ) );
  M = ( msq13 - msq23 - ( msqParent - msq3 ) * ( msq1 - msq2 ) / ( _bwm * _bwm ) ) *
      sqrt( ( 1 + R * R * p2R * p2R ) / ( 1 + R * R * p2 * p2 ) );
  // M = (msq13 - msq23 - (msqParent - msq3)*(msq1 - msq2)/(_p4_12.mass2()))*sqrt((1 +
  // R*R*p2R*p2R)/(1 + R*R*p2*p2));
 
  BW = sqrt( _gamma ) * M /
       ( ( _bwm * _bwm - _p4_12.mass2() ) - EvtComplex( 0.0, 1.0 ) * gam * _bwm );
 
  return BW;
}