00001 #include "EvtGenBase/EvtPatches.hh"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <assert.h>
00017 #include <cmath>
00018 #include <iostream>
00019
00020 #include <stdlib.h>
00021 #include "EvtGenBase/EvtParticle.hh"
00022 #include "EvtGenBase/EvtGenKine.hh"
00023 #include "EvtGenBase/EvtPDL.hh"
00024 #include "EvtGenBase/EvtReport.hh"
00025 #include "EvtGenBase/EvtMatrix.hh"
00026 #include "EvtGenBase/EvtDalitzReso.hh"
00027
00028 #include "EvtGenBase/EvtdFunction.hh"
00029 #include "EvtGenBase/EvtCyclic3.hh"
00030
00031 #define PRECISION ( 1.e-3 )
00032
00033 using EvtCyclic3::Index;
00034 using EvtCyclic3::Pair;
00035
00036
00037
00038 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
00039 EvtSpinType::spintype spin, double m0, double g0, NumType typeN)
00040 : _dp(dp),
00041 _pairAng(pairAng),
00042 _pairRes(pairRes),
00043 _spin(spin),
00044 _typeN(typeN),
00045 _m0(m0),_g0(g0),
00046 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
00047 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
00048 _g1(-1.),_g2(-1.),_coupling2(Undefined),
00049 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
00050 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
00051 {
00052 _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
00053 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
00054 _vb.set_f( 0.0 );
00055 _vd.set_f( 1.5 );
00056 assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);
00057 }
00058
00059
00060
00061 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
00062 EvtSpinType::spintype spin, double m0, double g0, NumType typeN,
00063 double m0_mix, double g0_mix, double delta_mix, EvtComplex amp_mix)
00064 : _dp(dp),
00065 _pairAng(pairAng),
00066 _pairRes(pairRes),
00067 _spin(spin),
00068 _typeN(typeN),
00069 _m0(m0),_g0(g0),
00070 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
00071 _m0_mix(m0_mix),_g0_mix(g0_mix),_delta_mix(delta_mix),_amp_mix(amp_mix),
00072 _g1(-1.),_g2(-1.),_coupling2(Undefined),
00073 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
00074 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
00075 {
00076 _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
00077 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
00078 _vb.set_f( 0.0 );
00079 _vd.set_f( 1.5 );
00080
00081 assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);
00082 }
00083
00084
00085 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
00086 EvtSpinType::spintype spin, double m0, NumType typeN, double g1, double g2, CouplingType coupling2)
00087 : _dp(dp),
00088 _pairAng(pairAng),
00089 _pairRes(pairRes),
00090 _spin(spin),
00091 _typeN(typeN),
00092 _m0(m0),_g0(-1.),
00093 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
00094 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
00095 _g1(g1),_g2(g2),_coupling2(coupling2),
00096 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
00097 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
00098 {
00099 _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
00100 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
00101 _vb.set_f( 0.0 );
00102 _vd.set_f( 1.5 );
00103 assert(_coupling2 != Undefined);
00104 assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);
00105 assert(_typeN != LASS);
00106 assert(_typeN != NBW);
00107 }
00108
00109
00110
00111 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes, std::string nameIndex, NumType typeN,
00112 EvtComplex fr12prod, EvtComplex fr13prod, EvtComplex fr14prod, EvtComplex fr15prod, double s0prod)
00113 : _dp(dp),
00114 _pairRes(pairRes),
00115 _typeN(typeN),
00116 _m0(0.),_g0(0.),
00117 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
00118 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
00119 _g1(-1.),_g2(-1.),_coupling2(Undefined),
00120 _kmatrix_index(-1),_fr12prod(fr12prod),_fr13prod(fr13prod),_fr14prod(fr14prod),_fr15prod(fr15prod),_s0prod(s0prod),
00121 _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
00122 {
00123 assert(_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II);
00124 _spin=EvtSpinType::SCALAR;
00125 if (nameIndex=="Pole1") _kmatrix_index=1;
00126 else if (nameIndex=="Pole2") _kmatrix_index=2;
00127 else if (nameIndex=="Pole3") _kmatrix_index=3;
00128 else if (nameIndex=="Pole4") _kmatrix_index=4;
00129 else if (nameIndex=="Pole5") _kmatrix_index=5;
00130 else if (nameIndex=="f11prod") _kmatrix_index=6;
00131 else assert(0);
00132 }
00133
00134
00135
00136 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes,
00137 double m0, double g0, double a, double r, double B, double phiB, double R, double phiR)
00138 : _dp(dp),
00139 _pairRes(pairRes),
00140 _typeN(LASS),
00141 _m0(m0),_g0(g0),
00142 _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
00143 _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
00144 _g1(-1.),_g2(-1.),_coupling2(Undefined),
00145 _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
00146 _a(a),_r(r),_Blass(B),_phiB(phiB),_R(R),_phiR(phiR)
00147 {
00148 _spin=EvtSpinType::SCALAR;
00149 _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
00150 _vd.set_f( 1.5 );
00151 }
00152
00153
00154
00155 EvtDalitzReso::EvtDalitzReso(const EvtDalitzReso& other)
00156 : _dp(other._dp),
00157 _pairAng(other._pairAng),
00158 _pairRes(other._pairRes),
00159 _spin(other._spin),
00160 _typeN(other._typeN),
00161 _m0(other._m0),_g0(other._g0),
00162 _vb(other._vb),_vd(other._vd),
00163 _massFirst(other._massFirst),_massSecond(other._massSecond),
00164 _m0_mix(other._m0_mix),_g0_mix(other._g0_mix),_delta_mix(other._delta_mix),_amp_mix(other._amp_mix),
00165 _g1(other._g1),_g2(other._g2),_coupling2(other._coupling2),
00166 _kmatrix_index(other._kmatrix_index),
00167 _fr12prod(other._fr12prod),_fr13prod(other._fr13prod),_fr14prod(other._fr14prod),_fr15prod(other._fr15prod),
00168 _s0prod(other._s0prod),
00169 _a(other._a),_r(other._r),_Blass(other._Blass),_phiB(other._phiB),_R(other._R),_phiR(other._phiR)
00170 {}
00171
00172
00173 EvtDalitzReso::~EvtDalitzReso()
00174 {}
00175
00176
00177 EvtComplex EvtDalitzReso::evaluate(const EvtDalitzPoint& x)
00178 {
00179 double m = sqrt(x.q(_pairRes));
00180
00181
00182 if (_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II)
00183 return Fvector( m*m, _kmatrix_index );
00184
00185 if (_typeN==LASS)
00186 return lass(m*m);
00187
00188 EvtComplex amp(1.0,0.0);
00189
00190 if (_dp.bigM() != x.bigM()) _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),x.bigM(),_spin);
00191 EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
00192 EvtTwoBodyKine vd(_massFirst,_massSecond,m);
00193
00194 EvtComplex prop(0,0);
00195 if (_typeN==NBW) {
00196 prop = propBreitWigner(_m0,_g0,m);
00197 } else if (_typeN==GAUSS_CLEO || _typeN==GAUSS_CLEO_ZEMACH) {
00198 prop = propGauss(_m0,_g0,m);
00199 } else {
00200 if (_coupling2==Undefined) {
00201
00202 double g = (_g0<=0. || _vd.pD()<=0.)? -_g0 : _g0*_vd.widthFactor(vd);
00203 if (_typeN==GS_CLEO || _typeN==GS_CLEO_ZEMACH) {
00204
00205 assert(_massFirst==_massSecond);
00206 prop = propGounarisSakurai(_m0,fabs(_g0),_vd.pD(),m,g,vd.p());
00207 } else {
00208
00209 prop = propBreitWignerRel(_m0,g,m);
00210 }
00211 } else {
00212
00213 EvtComplex G1,G2;
00214 switch (_coupling2) {
00215 case PicPic: {
00216 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00217 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00218 G2 = _g2*_g2*psFactor(mPic,mPic,m);
00219 break;
00220 }
00221 case PizPiz: {
00222 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00223 static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
00224 G2 = _g2*_g2*psFactor(mPiz,mPiz,m);
00225 break;
00226 }
00227 case PiPi: {
00228 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00229 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00230 static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
00231 G2 = _g2*_g2*psFactor(mPic,mPic,mPiz,mPiz,m);
00232 break;
00233 }
00234 case KcKc: {
00235 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00236 static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
00237 G2 = _g2*_g2*psFactor(mKc,mKc,m);
00238 break;
00239 }
00240 case KzKz: {
00241 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00242 static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
00243 G2 = _g2*_g2*psFactor(mKz,mKz,m);
00244 break;
00245 }
00246 case KK: {
00247 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00248 static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
00249 static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
00250 G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
00251 break;
00252 }
00253 case EtaPic: {
00254 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00255 static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
00256 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00257 G2 = _g2*_g2*psFactor(mEta,mPic,m);
00258 break;
00259 }
00260 case EtaPiz: {
00261 G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00262 static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
00263 static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
00264 G2 = _g2*_g2*psFactor(mEta,mPiz,m);
00265 break;
00266 }
00267 case PicPicKK: {
00268 static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00269
00270 G1 = _g1*psFactor(mPic,mPic,m);
00271 static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
00272 static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
00273
00274 G2 = _g2*psFactor(mKc,mKc,mKz,mKz,m);
00275 break;
00276 }
00277 default:
00278 std::cout << "EvtDalitzReso:evaluate(): PANIC, wrong coupling2 state." << std::endl;
00279 assert(0);
00280 break;
00281 }
00282
00283 if (_coupling2 != WA76)
00284 prop = _g1*propBreitWignerRelCoupled(_m0,G1,G2,m);
00285 }
00286 }
00287 amp *= prop;
00288
00289
00290 amp *= _vb.formFactor(vb);
00291 amp *= _vd.formFactor(vd);
00292
00293
00294 amp *= numerator(x,vb,vd);
00295
00296
00297 if (_m0_mix>0.) {
00298 EvtComplex prop_mix;
00299 if (_typeN==NBW) {
00300 prop_mix = propBreitWigner(_m0_mix,_g0_mix,m);
00301 } else {
00302 assert(_g1<0.);
00303 double g_mix = _g0_mix*_vd.widthFactor(vd);
00304 prop_mix = propBreitWignerRel(_m0_mix,g_mix,m);
00305 }
00306 amp *= mixFactor(prop,prop_mix);
00307 }
00308
00309 return amp;
00310 }
00311
00312
00313 EvtComplex EvtDalitzReso::psFactor(double & ma, double & mb, double& m)
00314 {
00315 if (m>(ma+mb)) {
00316 EvtTwoBodyKine vd(ma,mb,m);
00317 return EvtComplex(0,2*vd.p()/m);
00318 } else {
00319
00320 double s = m*m;
00321 double phaseFactor_analyticalCont = -0.5*(sqrt(4*ma*ma/s-1)+sqrt(4*mb*mb/s-1));
00322 return EvtComplex(phaseFactor_analyticalCont,0);
00323 }
00324 }
00325
00326
00327 EvtComplex EvtDalitzReso::psFactor(double & ma1,double & mb1, double & ma2, double & mb2, double& m)
00328 {
00329 return 0.5*(psFactor(ma1,mb1,m)+psFactor(ma2,mb2,m));
00330 }
00331
00332
00333 EvtComplex EvtDalitzReso::propGauss(const double& m0, const double& s0, const double& m)
00334 {
00335
00336 double gauss = 1./sqrt(EvtConst::twoPi)/s0*exp(-(m-m0)*(m-m0)/2./(s0*s0));
00337 return EvtComplex(gauss,0.);
00338 }
00339
00340
00341 EvtComplex EvtDalitzReso::propBreitWigner(const double& m0, const double& g0, const double& m)
00342 {
00343
00344 return sqrt(g0/EvtConst::twoPi)/(m-m0-EvtComplex(0.0,g0/2.));
00345 }
00346
00347
00348 EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const double& g0, const double& m)
00349 {
00350
00351 return 1./(m0*m0-m*m-EvtComplex(0.,m0*g0));
00352 }
00353
00354
00355
00356 EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const EvtComplex& g0, const double& m)
00357 {
00358
00359 return 1./(m0*m0-m*m-EvtComplex(0.,m0)*g0);
00360 }
00361
00362
00363 EvtComplex EvtDalitzReso::propBreitWignerRelCoupled(const double& m0, const EvtComplex& g1, const EvtComplex& g2, const double& m)
00364 {
00365
00366 return 1./(m0*m0-m*m-(g1+g2));
00367 }
00368
00369 EvtComplex EvtDalitzReso::propGounarisSakurai(const double& m0, const double& g0, const double& k0,
00370 const double& m, const double& g, const double& k)
00371 {
00372
00373
00374 return (1.+GS_d(m0,k0)*g0/m0)/(m0*m0-m*m-EvtComplex(0.,m*g)+GS_f(m0,g0,k0,m,k));
00375 }
00376
00377
00378 inline double EvtDalitzReso::GS_f(const double& m0, const double& g0, const double& k0, const double& m, const double& k)
00379 {
00380
00381
00382
00383
00384 return g0*m0*m0/(k0*k0*k0)*( k*k*(GS_h(m,k)-GS_h(m0,k0)) + (m0*m0-m*m)*k0*k0*GS_dhods(m0,k0) );
00385 }
00386
00387 inline double EvtDalitzReso::GS_h(const double& m, const double& k)
00388 {return 2./EvtConst::pi*k/m*log((m+2.*k)/(2.*_massFirst)) ;}
00389
00390 inline double EvtDalitzReso::GS_dhods(const double& m0, const double& k0)
00391 {return GS_h(m0,k0)*( 0.125/(k0*k0) - 0.5/(m0*m0) ) + 0.5/(EvtConst::pi*m0*m0) ;}
00392
00393 inline double EvtDalitzReso::GS_d(const double& m0, const double& k0)
00394 {return 3./EvtConst::pi*_massFirst*_massFirst/(k0*k0)*log((m0+2.*k0)/(2.*_massFirst)) +
00395 m0/(2.*EvtConst::pi*k0) - _massFirst*_massFirst*m0/(EvtConst::pi*k0*k0*k0) ;}
00396
00397
00398 EvtComplex EvtDalitzReso::numerator(const EvtDalitzPoint& x, const EvtTwoBodyKine& vb, const EvtTwoBodyKine& vd)
00399 {
00400 EvtComplex ret(0.,0.);
00401
00402
00403 if(NBW == _typeN) {
00404 ret = angDep(x);
00405 }
00406
00407
00408 else if(RBW_ZEMACH == _typeN) {
00409 ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
00410 }
00411
00412
00413 else if(RBW_ZEMACH2 == _typeN) {
00414 ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*_vb.phaseSpaceFactor(vb,EvtTwoBodyKine::AB)*angDep(x);
00415 if(_spin == EvtSpinType::VECTOR) {
00416 ret *= -4.;
00417 } else if(_spin == EvtSpinType::TENSOR) {
00418 ret *= 16./3.;
00419 } else if(_spin != EvtSpinType::SCALAR)
00420 assert(0);
00421 }
00422
00423
00424 else if(RBW_KUEHN == _typeN) {
00425 ret = _m0*_m0 * angDep(x);
00426 }
00427
00428
00429 else if( ( RBW_CLEO == _typeN ) || ( GS_CLEO == _typeN ) ||
00430 ( RBW_CLEO_ZEMACH == _typeN ) || ( GS_CLEO_ZEMACH == _typeN ) ||
00431 ( GAUSS_CLEO == _typeN ) || ( GAUSS_CLEO_ZEMACH == _typeN)) {
00432
00433 Index iA = other(_pairAng);
00434 Index iB = common(_pairRes,_pairAng);
00435 Index iC = other(_pairRes);
00436
00437 double M = x.bigM();
00438 double mA = x.m(iA);
00439 double mB = x.m(iB);
00440 double mC = x.m(iC);
00441 double qAB = x.q(combine(iA,iB));
00442 double qBC = x.q(combine(iB,iC));
00443 double qCA = x.q(combine(iC,iA));
00444
00445 double M2 = M*M;
00446 double m02 = ((RBW_CLEO_ZEMACH == _typeN)||(GS_CLEO_ZEMACH == _typeN)||(GAUSS_CLEO_ZEMACH == _typeN))? qAB : _m0*_m0;
00447 double mA2 = mA*mA;
00448 double mB2 = mB*mB;
00449 double mC2 = mC*mC;
00450
00451 if (_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
00452 else if(_spin == EvtSpinType::VECTOR) {
00453 ret = qCA - qBC + (M2 - mC2)*(mB2 - mA2)/m02;
00454 } else if(_spin == EvtSpinType::TENSOR) {
00455 double x1 = qBC - qCA + (M2 - mC2)*(mA2 - mB2)/m02;
00456 double x2 = M2 - mC2;
00457 double x3 = qAB - 2*M2 - 2*mC2 + x2*x2/m02;
00458 double x4 = mA2 - mB2;
00459 double x5 = qAB - 2*mB2 - 2*mA2 + x4*x4/m02;
00460 ret = x1*x1 - x3*x5/3.;
00461 } else assert(0);
00462 }
00463
00464 return ret;
00465 }
00466
00467
00468 double EvtDalitzReso::angDep(const EvtDalitzPoint& x)
00469 {
00470
00471
00472 double cosTh = x.cosTh(_pairAng,_pairRes);
00473 if(fabs(cosTh) > 1.) {
00474 report(INFO,"EvtGen") << "cosTh " << cosTh << std::endl;
00475 assert(0);
00476 }
00477
00478
00479 return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
00480 }
00481
00482
00483 EvtComplex EvtDalitzReso::mixFactor(EvtComplex prop, EvtComplex prop_mix)
00484 {
00485 double Delta = _delta_mix*(_m0+_m0_mix);
00486 return 1/(1-Delta*Delta*prop*prop_mix)*(1+_amp_mix*Delta*prop_mix);
00487 }
00488
00489
00490
00491 EvtComplex EvtDalitzReso::Fvector( double s, int index )
00492 {
00493 assert(index>=1 && index<=6);
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504 double g[5][5];
00505 double ma[5];
00506
00507 int solution = (_typeN==K_MATRIX)? 3 : ( (_typeN==K_MATRIX_I)? 1 : ( (_typeN==K_MATRIX_II)? 2 : 0 ) ) ;
00508 if (solution==0) { std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! " << std::endl; abort(); }
00509
00510 if (solution == 3 ) {
00511
00512
00513
00514 g[0][0]=0.22889;
00515 g[1][0]=0.94128;
00516 g[2][0]=0.36856;
00517 g[3][0]=0.33650;
00518 g[4][0]=0.18171;
00519
00520 g[0][1]=-0.55377;
00521 g[1][1]=0.55095;
00522 g[2][1]=0.23888;
00523 g[3][1]=0.40907;
00524 g[4][1]=-0.17558;
00525
00526 g[0][2]=0;
00527 g[1][2]=0;
00528 g[2][2]=0.55639;
00529 g[3][2]=0.85679;
00530 g[4][2]=-0.79658;
00531
00532 g[0][3]=-0.39899;
00533 g[1][3]=0.39065;
00534 g[2][3]=0.18340;
00535 g[3][3]=0.19906;
00536 g[4][3]=-0.00355;
00537
00538 g[0][4]=-0.34639;
00539 g[1][4]=0.31503;
00540 g[2][4]=0.18681;
00541 g[3][4]=-0.00984;
00542 g[4][4]=0.22358;
00543
00544
00545 ma[0]=0.651;
00546 ma[1]=1.20360;
00547 ma[2]=1.55817;
00548 ma[3]=1.21000;
00549 ma[4]=1.82206;
00550
00551 } else if (solution == 1) {
00552
00553
00554
00555 g[0][0]=0.31896;
00556 g[1][0]=0.85963;
00557 g[2][0]=0.47993;
00558 g[3][0]=0.45121;
00559 g[4][0]=0.39391;
00560
00561 g[0][1]=-0.49998;
00562 g[1][1]=0.52402;
00563 g[2][1]=0.40254;
00564 g[3][1]=0.42769;
00565 g[4][1]=-0.30860;
00566
00567 g[0][2]=0;
00568 g[1][2]=0;
00569 g[2][2]=1.0;
00570 g[3][2]=1.15088;
00571 g[4][2]=0.33999;
00572
00573 g[0][3]=-0.21554;
00574 g[1][3]=0.38093;
00575 g[2][3]=0.21811;
00576 g[3][3]=0.22925;
00577 g[4][3]=0.06919;
00578
00579 g[0][4]=-0.18294;
00580 g[1][4]=0.23788;
00581 g[2][4]=0.05454;
00582 g[3][4]=0.06444;
00583 g[4][4]=0.32620;
00584
00585
00586 ma[0]=0.7369;
00587 ma[1]=1.24347;
00588 ma[2]=1.62681;
00589 ma[3]=1.21900;
00590 ma[4]=1.74932;
00591
00592 } else if (solution == 2) {
00593
00594
00595
00596 g[0][0]=0.26014;
00597 g[1][0]=0.95289;
00598 g[2][0]=0.46244;
00599 g[3][0]=0.41848;
00600 g[4][0]=0.01804;
00601
00602 g[0][1]=-0.57849;
00603 g[1][1]=0.55887;
00604 g[2][1]=0.31712;
00605 g[3][1]=0.49910;
00606 g[4][1]=-0.28430;
00607
00608 g[0][2]=0;
00609 g[1][2]=0;
00610 g[2][2]=0.70340;
00611 g[3][2]=0.96819;
00612 g[4][2]=-0.90100;
00613
00614 g[0][3]=-0.32936;
00615 g[1][3]=0.39910;
00616 g[2][3]=0.22963;
00617 g[3][3]=0.24415;
00618 g[4][3]=-0.07252;
00619
00620 g[0][4]=-0.30906;
00621 g[1][4]=0.31143;
00622 g[2][4]=0.19802;
00623 g[3][4]=-0.00522;
00624 g[4][4]=0.17097;
00625
00626
00627 ma[0]=0.67460;
00628 ma[1]=1.21094;
00629 ma[2]=1.57896;
00630 ma[3]=1.21900;
00631 ma[4]=1.86602;
00632 }
00633
00634
00635 double rho1sq,rho2sq,rho4sq,rho5sq;
00636 EvtComplex rho[5];
00637 double f[5][5];
00638
00639
00640 double mpi=0.13957;
00641 double mK=0.493677;
00642 double meta=0.54775;
00643 double metap=0.95778;
00644
00645
00646 EvtComplex K[5][5];
00647 for(int i=0;i<5;i++) {
00648 for(int j=0;j<5;j++) {
00649 K[i][j]=EvtComplex(0,0);
00650 f[i][j]=0;
00651 }
00652 }
00653
00654
00655 double s_scatt=0.0 ;
00656 if (solution == 3)
00657 s_scatt=-3.92637;
00658 else if (solution == 1)
00659 s_scatt= -5.0 ;
00660 else if (solution == 2)
00661 s_scatt= -5.0 ;
00662 double sa=1.0;
00663 double sa_0=-0.15;
00664 if (solution == 3) {
00665 f[0][0]=0.23399;
00666 f[0][1]=0.15044;
00667 f[0][2]=-0.20545;
00668 f[0][3]=0.32825;
00669 f[0][4]=0.35412;
00670 }else if (solution == 1) {
00671 f[0][0]=0.04214;
00672 f[0][1]=0.19865;
00673 f[0][2]=-0.63764;
00674 f[0][3]=0.44063;
00675 f[0][4]=0.36717;
00676 }else if (solution == 2) {
00677 f[0][0]=0.26447;
00678 f[0][1]=0.10400;
00679 f[0][2]=-0.35445;
00680 f[0][3]=0.31596;
00681 f[0][4]=0.42483;
00682 }
00683 f[1][0]=f[0][1];
00684 f[2][0]=f[0][2];
00685 f[3][0]=f[0][3];
00686 f[4][0]=f[0][4];
00687
00688
00689
00690 rho1sq = 1. - pow( mpi + mpi, 2 ) / s;
00691 if( rho1sq >= 0 )
00692 rho[ 0 ] = EvtComplex( sqrt( rho1sq ), 0 );
00693 else
00694 rho[ 0 ] = EvtComplex( 0, sqrt( -rho1sq ) );
00695
00696 rho2sq = 1. - pow( mK + mK, 2 ) / s;
00697 if( rho2sq >= 0 )
00698 rho[ 1 ] = EvtComplex( sqrt( rho2sq ), 0 );
00699 else
00700 rho[ 1 ] = EvtComplex( 0, sqrt( -rho2sq ) );
00701
00702
00703
00704 if( s <= 1 )
00705 {
00706 double real = 1.2274 + .00370909 / ( s * s ) - .111203 / s - 6.39017 * s + 16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
00707 double cont32 = sqrt(1.0-(16.0*mpi*mpi));
00708 rho[ 2 ] = EvtComplex( cont32 * real, 0 );
00709 }
00710 else
00711 rho[ 2 ] = EvtComplex( sqrt( 1. - 16. * mpi * mpi / s ), 0 );
00712
00713 rho4sq = 1. - pow( meta + meta, 2 ) / s;
00714 if( rho4sq >= 0 )
00715 rho[ 3 ] = EvtComplex( sqrt( rho4sq ), 0 );
00716 else
00717 rho[ 3 ] = EvtComplex( 0, sqrt( -rho4sq ) );
00718
00719 rho5sq = 1. - pow( meta + metap, 2 ) / s;
00720 if( rho5sq >= 0 )
00721 rho[ 4 ] = EvtComplex( sqrt( rho5sq ), 0 );
00722 else
00723 rho[ 4 ] = EvtComplex( 0, sqrt( -rho5sq ) );
00724
00725 double smallTerm = 1;
00726
00727
00728 for ( int pole = 0; pole < 5; pole++ )
00729 if ( fabs( pow( ma[ pole ], 2 ) - s ) < PRECISION )
00730 smallTerm = pow( ma[ pole ], 2 ) - s;
00731
00732
00733
00734 for(int i=0;i<5;i++) {
00735 for(int j=0;j<5;j++) {
00736 for (int pole_index=0;pole_index<5;pole_index++) {
00737 double A=g[pole_index][i]*g[pole_index][j];
00738 double B=ma[pole_index]*ma[pole_index]-s;
00739
00740 if ( fabs( B ) < PRECISION )
00741 K[ i ][ j ] += EvtComplex( A , 0 );
00742 else
00743 K[ i ][ j ] += EvtComplex( A / B, 0 ) * smallTerm;
00744 }
00745 }
00746 }
00747
00748
00749 for(int i=0;i<5;i++) {
00750 for(int j=0;j<5;j++) {
00751 double C=f[i][j]*(1.0-s_scatt);
00752 double D=(s-s_scatt);
00753 K[ i ][ j ] += EvtComplex( C / D, 0 ) * smallTerm;
00754 }
00755 }
00756
00757
00758
00759 for(int i=0;i<5;i++) {
00760 for(int j=0;j<5;j++) {
00761 double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
00762 double F=(s-sa_0);
00763 K[ i ][ j ] *= EvtComplex(E/F,0);
00764 }
00765 }
00766
00767
00768
00769 static EvtMatrix< EvtComplex > mat;
00770 mat.setRange( 5 );
00771
00772 for ( int row = 0; row < 5; row++ )
00773 for ( int col = 0; col < 5; col++ )
00774 mat( row, col ) = ( row == col ) * smallTerm - EvtComplex( 0., 1. ) * K[ row ][ col ] * rho[ col ];
00775
00776
00777 EvtMatrix< EvtComplex >* matInverse = mat.inverse();
00778 vector< EvtComplex > U1j;
00779 for ( int j = 0; j < 5; j++ )
00780 U1j.push_back( (*matInverse)[ 0 ][ j ] );
00781
00782 delete matInverse;
00783
00784
00785 EvtComplex value( 0, 0 );
00786 if (index<=5) {
00787
00788 for(int j=0;j<5;j++) {
00789 EvtComplex top = U1j[j]*g[index-1][j];
00790 double bottom = ma[index-1]*ma[index-1]-s;
00791
00792 if ( fabs( bottom ) < PRECISION )
00793 value += top;
00794 else
00795 value += top / bottom * smallTerm;
00796 }
00797 } else {
00798
00799 value += U1j[0];
00800 value += U1j[1]*_fr12prod;
00801 value += U1j[2]*_fr13prod;
00802 value += U1j[3]*_fr14prod;
00803 value += U1j[4]*_fr15prod;
00804
00805 value *= (1-_s0prod)/(s-_s0prod) * smallTerm;
00806 }
00807
00808 return value;
00809 }
00810
00811
00812
00813 EvtComplex EvtDalitzReso::lass(double s)
00814 {
00815 EvtTwoBodyKine vd(_massFirst,_massSecond, sqrt(s));
00816 double q = vd.p();
00817 double GammaM = _g0*_vd.widthFactor(vd);
00818
00819
00820 double cot_deltaB = 1.0/(_a*q) + 0.5*_r*q;
00821 double deltaB = atan( 1.0/cot_deltaB);
00822 double totalB = deltaB + _phiB ;
00823
00824
00825 double deltaR = atan((_m0*GammaM/(_m0*_m0 - s)));
00826 double totalR = deltaR + _phiR ;
00827
00828
00829 EvtComplex bkgB,resT;
00830 bkgB = EvtComplex(_Blass*sin(totalB),0)*EvtComplex(cos(totalB),sin(totalB));
00831 resT = EvtComplex(_R*sin(deltaR),0)*EvtComplex(cos(totalR),sin(totalR))*EvtComplex(cos(2*totalB),sin(2*totalB));
00832 EvtComplex T = bkgB + resT;
00833
00834 return T;
00835 }
00836
00837
00838
00839