00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <math.h>
00023 #include "EvtGenBase/EvtVector4R.hh"
00024 #include "EvtGenBase/EvtKine.hh"
00025 #include "EvtGenBase/EvtComplex.hh"
00026 #include "EvtGenBase/EvtResonance2.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtConst.hh"
00029
00030 EvtResonance2::~EvtResonance2(){}
00031
00032
00033
00034 EvtResonance2& EvtResonance2::operator = ( const EvtResonance2 &n)
00035 {
00036 if ( &n == this ) return *this;
00037 _p4_p = n._p4_p;
00038 _p4_d1 = n._p4_d1;
00039 _p4_d2 = n._p4_d2;
00040 _ampl = n._ampl;
00041 _theta = n._theta;
00042 _gamma = n._gamma;
00043 _spin = n._spin;
00044 _bwm = n._bwm;
00045 _invmass_angdenom = n._invmass_angdenom;
00046 return *this;
00047 }
00048
00049
00050
00051 EvtResonance2::EvtResonance2(const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
00052 const EvtVector4R& p4_d2, double ampl,
00053 double theta, double gamma, double bwm, int spin,
00054 bool invmass_angdenom):
00055 _p4_p(p4_p),_p4_d1(p4_d1), _p4_d2(p4_d2),_ampl(ampl), _theta(theta),
00056 _gamma(gamma), _bwm(bwm), _spin(spin), _invmass_angdenom(invmass_angdenom) {}
00057
00058
00059
00060 EvtComplex EvtResonance2::resAmpl() {
00061
00062 double pi180inv = 1.0/EvtConst::radToDegrees;
00063
00064 EvtComplex ampl;
00065 EvtVector4R p4_d3 = _p4_p-_p4_d1-_p4_d2;
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 double mAB=(_p4_d1+_p4_d2).mass();
00082 double mBC=(_p4_d2+p4_d3).mass();
00083 double mAC=(_p4_d1+p4_d3).mass();
00084 double mA=_p4_d1.mass();
00085 double mB=_p4_d2.mass();
00086 double mD=_p4_p.mass();
00087 double mC=p4_d3.mass();
00088
00089 double mR=_bwm;
00090 double gammaR=_gamma;
00091 double mdenom = _invmass_angdenom ? mAB : mR;
00092 double pAB=sqrt( (((mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0) -
00093 mA*mA*mB*mB)/(mAB*mAB));
00094 double pR=sqrt( (((mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0) -
00095 mA*mA*mB*mB)/(mR*mR));
00096
00097 double pD= (((mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0) -
00098 mR*mR*mC*mC)/(mD*mD);
00099 if ( pD>0 ) { pD=sqrt(pD); } else {pD=0;}
00100 double pDAB=sqrt( (((mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0) -
00101 mAB*mAB*mC*mC)/(mD*mD));
00102
00103
00104
00105
00106
00107 double fR=1;
00108 double fD=1;
00109 int power=0;
00110 switch (_spin) {
00111 case 0:
00112 fR=1.0;
00113 fD=1.0;
00114 power=1;
00115
00116 break;
00117 case 1:
00118 fR=sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
00119 fD=sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
00120
00121 power=3;
00122 break;
00123 case 2:
00124 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)) );
00125 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)) );
00126 power=5;
00127 break;
00128 default:
00129 report(INFO,"EvtGen") << "Incorrect spin in EvtResonance22.cc\n";
00130 }
00131
00132 double gammaAB= gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
00133
00134 switch (_spin) {
00135 case 0:
00136 ampl=_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00137 fR*fD/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB));
00138 break;
00139 case 1:
00140 ampl=_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00141 (fR*fD*(mAC*mAC-mBC*mBC+((mD*mD-mC*mC)*(mB*mB-mA*mA)/(mdenom*mdenom)))/
00142
00143 (mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB)));
00144 break;
00145 case 2:
00146 ampl=_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00147 fR*fD/(mR*mR-mAB*mAB-EvtComplex(0.0,mR*gammaAB))*
00148 (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mdenom*mdenom)),2)-
00149 (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mdenom, 2))*
00150 (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mdenom,2)));
00151 break;
00152
00153 default:
00154 report(INFO,"EvtGen") << "Incorrect spin in EvtResonance22.cc\n";
00155 }
00156
00157
00158 return ampl;
00159 }
00160
00161
00162