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/EvtResonance.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtConst.hh"
00029 using std::endl;
00030
00031 EvtResonance::~EvtResonance(){}
00032
00033
00034
00035 EvtResonance& EvtResonance::operator = ( const EvtResonance &n)
00036 {
00037 if ( &n == this ) return *this;
00038 _p4_p = n._p4_p;
00039 _p4_d1 = n._p4_d1;
00040 _p4_d2 = n._p4_d2;
00041 _ampl = n._ampl;
00042 _theta = n._theta;
00043 _gamma = n._gamma;
00044 _spin = n._spin;
00045 _bwm = n._bwm;
00046 return *this;
00047 }
00048
00049
00050
00051 EvtResonance::EvtResonance(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 _p4_p(p4_p),_p4_d1(p4_d1), _p4_d2(p4_d2),_ampl(ampl), _theta(theta),
00055 _gamma(gamma), _bwm(bwm), _spin(spin) {}
00056
00057
00058
00059 EvtComplex EvtResonance::resAmpl() {
00060
00061 double pi180inv = 1.0/EvtConst::radToDegrees;
00062
00063 EvtComplex ampl;
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 double cos_phi_0 = EvtDecayAngle(_p4_p, _p4_d1+_p4_d2, _p4_d1);
00076
00077
00078 switch (_spin) {
00079
00080 case 0 :
00081 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00082 sqrt(_gamma/EvtConst::twoPi)*
00083 (1.0/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,0.5*_gamma))));
00084 break;
00085
00086 case 1 :
00087 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00088 sqrt(_gamma/EvtConst::twoPi)*
00089 (cos_phi_0/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,0.5*_gamma))));
00090 break;
00091
00092 case 2:
00093 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00094 sqrt(_gamma/EvtConst::twoPi)*
00095 ((1.5*cos_phi_0*cos_phi_0-0.5)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0, 0.5*_gamma))));
00096 break;
00097
00098 case 3:
00099 ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00100 sqrt(_gamma/EvtConst::twoPi)*
00101 ((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))));
00102 break;
00103
00104 default:
00105 report(DEBUG,"EvtGen") << "EvtGen: wrong spin in EvtResonance" << endl;
00106 ampl = EvtComplex(0.0);
00107 break;
00108
00109 }
00110
00111 return ampl;
00112 }
00113
00114 EvtComplex EvtResonance::relBrWig(int i) {
00115
00116
00117
00118
00119 EvtComplex BW;
00120 EvtVector4R _p4_d3 = _p4_p-_p4_d1-_p4_d2;
00121 EvtVector4R _p4_12 = _p4_d1 + _p4_d2;
00122
00123
00124
00125
00126
00127 double msq13 = (_p4_d1 + _p4_d3).mass2();
00128 double msq23 = (_p4_d2 + _p4_d3).mass2();
00129 double msqParent = _p4_p.mass2();
00130 double msq1 = _p4_d1.mass2();
00131 double msq2 = _p4_d2.mass2();
00132 double msq3 = _p4_d3.mass2();
00133
00134 double M;
00135
00136 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());
00137
00138 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);
00139
00140 double gam, R;
00141
00142 if (i == 1) {
00143
00144
00145 R = 2.0/(0.197);
00146
00147 }
00148 else R = 5.0/(0.197);
00149
00150 gam = _gamma*(_bwm/_p4_12.mass())*(p2/p2R)*(p2/p2R)*(p2/p2R)*((1 + R*R*p2R*p2R)/(1 + R*R*p2*p2));
00151 M = (msq13 - msq23 - (msqParent - msq3)*(msq1 - msq2)/(_bwm*_bwm))*sqrt((1 + R*R*p2R*p2R)/(1 + R*R*p2*p2));
00152
00153
00154 BW = sqrt(_gamma)*M/((_bwm*_bwm - _p4_12.mass2()) - EvtComplex(0.0,1.0)*gam*_bwm);
00155
00156 return BW;
00157
00158 }
00159
00160