TEvtGen/EvtGenBase/EvtResonance2.cpp

   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: EvtResonance2.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/EvtResonance2.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtConst.hh"
  29
  30 EvtResonance2::~EvtResonance2(){}
  31
  32 EvtResonance2& EvtResonance2::operator = ( const EvtResonance2  &n)
  33 {
  34   if ( &n == this ) return *this;
  35   _p4_p = n._p4_p;
  36   _p4_d1 = n._p4_d1;
  37   _p4_d2 = n._p4_d2;
  38   _ampl = n._ampl;
  39   _theta = n._theta;
  40   _gamma = n._gamma;
  41   _spin = n._spin;
  42   _bwm = n._bwm;
  43   _invmass_angdenom = n._invmass_angdenom;
  44    return  *this;
  45 }
  46
  47
  48 EvtResonance2::EvtResonance2(const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
  49                              const  EvtVector4R& p4_d2, double ampl,
  50                              double theta, double gamma, double bwm, int spin,
  51                              bool invmass_angdenom):
  52   _p4_p(p4_p),_p4_d1(p4_d1), _p4_d2(p4_d2),_ampl(ampl), _theta(theta),
  53   _gamma(gamma), _bwm(bwm), _spin(spin), _invmass_angdenom(invmass_angdenom) {}
  54
  55
  56 EvtComplex EvtResonance2::resAmpl() {
  57
  58   double pi180inv = 1.0/EvtConst::radToDegrees;
  59
  60   EvtComplex ampl;
  61   EvtVector4R  p4_d3 = _p4_p-_p4_d1-_p4_d2;
  62
  63   //get cos of the angle between the daughters from their 4-momenta
  64   //and the 4-momentum of the parent
  65
  66   //in general, EvtDecayAngle(parent, part1+part2, part1) gives the angle
  67   //the missing particle (not listed in the arguments) makes
  68   //with part2 in the rest frame of both
  69   //listed particles (12)
  70
  71   //angle 3 makes with 2 in rest frame of 12 (CS3)
  72   //double cos_phi_0 = EvtDecayAngle(_p4_p, _p4_d1+_p4_d2, _p4_d1);
  73   //angle 3 makes with 1 in 12 is, of course, -cos_phi_0
  74
  75   //first compute several quantities...follow CLEO preprint 00-23
  76
  77   double mAB=(_p4_d1+_p4_d2).mass();
  78   double mBC=(_p4_d2+p4_d3).mass();
  79   double mAC=(_p4_d1+p4_d3).mass();
  80   double mA=_p4_d1.mass();
  81   double mB=_p4_d2.mass();
  82   double mD=_p4_p.mass();
  83   double mC=p4_d3.mass();
  84
  85   double mR=_bwm;
  86   double gammaR=_gamma;
  87   double mdenom = _invmass_angdenom ? mAB : mR;
  88   double pAB=sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) -
  89                     mA*mA*mB*mB)/(mAB*mAB));
  90   double pR=sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) -
  91                    mA*mA*mB*mB)/(mR*mR));
  92
  93   double pD= (((mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0) -
  94                    mR*mR*mC*mC)/(mD*mD);
  95   if ( pD>0 ) { pD=sqrt(pD); } else {pD=0;}
  96   double pDAB=sqrt( (((mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0) -
  97                    mAB*mAB*mC*mC)/(mD*mD));
  98
  99
 100
 101   double fR=1;
 102   double fD=1;
 103   int power=0;
 104   switch (_spin) {
 105   case 0:
 106     fR=1.0;
 107     fD=1.0;
 108     power=1;
 109     break;
 110   case 1:
 111     fR=sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
 112     fD=sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
 113     power=3;
 114     break;
 115   case 2:
 116     fR = sqrt( (9+3*pow((1.5*pR),2)+pow((1.5*pR),4))/(9+3*pow((1.5*pAB),2)+pow((1.5*pAB),4)) );
 117     fD = sqrt( (9+3*pow((5.0*pD),2)+pow((5.0*pD),4))/(9+3*pow((5.0*pDAB),2)+pow((5.0*pDAB),4)) );
 118     power=5;
 119     break;
 120   default:
 121     report(INFO,"EvtGen") << "Incorrect spin in EvtResonance22.cc\n";
 122   }
 123
 124   double gammaAB= gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
 125   switch (_spin) {
 126   case 0:
 127     ampl=_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
 128           fR*fD/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB));
 129     break;
 130   case 1:
 131     ampl=_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
 132       (fR*fD*(mAC*mAC-mBC*mBC+((mD*mD-mC*mC)*(mB*mB-mA*mA)/(mdenom*mdenom)))/
 133        (mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB)));
 134     break;
 135   case 2:
 136     ampl=_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
 137       fR*fD/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB))*
 138       (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mdenom*mdenom)),2)-
 139        (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mdenom, 2))*
 140        (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mdenom,2)));
 141   break;
 142
 143   default:
 144     report(INFO,"EvtGen") << "Incorrect spin in EvtResonance22.cc\n";
 145   }
 146
 147   return ampl;
 148 }
 149
 150
 151