TEvtGen/EvtGenBase/EvtDalitzReso.cpp

   1 #include "EvtGenBase/EvtPatches.hh"
   2 /*****************************************************************************
   3  * Project: BaBar detector at the SLAC PEP-II B-factory
   4  * Package: EvtGenBase
   5  *    File: $Id: EvtDalitzReso.cpp,v 1.1 2009-03-16 16:47:51 robbep Exp $
   6  *
   7  * Description:
   8  *   Class to compute Dalitz amplitudes based on many models that cannot be
   9  *     handled with EvtResonance.
  10  *
  11  * Modification history:
  12  *   Jordi Garra Ticó     2008/07/03         File created
  13  *****************************************************************************/
  14
  15
  16 #include <assert.h>
  17 #include <cmath>
  18 #include <iostream>
  19
  20 #include <stdlib.h>
  21 #include "EvtGenBase/EvtParticle.hh"
  22 #include "EvtGenBase/EvtGenKine.hh"
  23 #include "EvtGenBase/EvtPDL.hh"
  24 #include "EvtGenBase/EvtReport.hh"
  25 #include "EvtGenBase/EvtMatrix.hh"
  26 #include "EvtGenBase/EvtDalitzReso.hh"
  27
  28 #include "EvtGenBase/EvtdFunction.hh"
  29 #include "EvtGenBase/EvtCyclic3.hh"
  30
  31 #define PRECISION ( 1.e-3 )
  32
  33 using EvtCyclic3::Index;
  34 using EvtCyclic3::Pair;
  35
  36
  37 // single Breit-Wigner
  38 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
  39                              EvtSpinType::spintype spin, double m0, double g0, NumType typeN, double f_b, double f_d)
  40   : _dp(dp),
  41     _pairAng(pairAng),
  42     _pairRes(pairRes),
  43     _spin(spin),
  44     _typeN(typeN),
  45     _m0(m0),_g0(g0),
  46     _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
  47     _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
  48     _g1(-1.),_g2(-1.),_coupling2(Undefined),
  49     _f_b(f_b), _f_d(f_d),
  50     _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
  51     _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.),_cutoff(-1.), _scaleByMOverQ(false),
  52     _alpha(0.)
  53 {
  54   _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
  55   _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
  56   _vb.set_f( _f_b ); // Default values for Blatt-Weisskopf factors are 0.0 and 1.5.
  57   _vd.set_f( _f_d );
  58   assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);  // single BW cannot be K-matrix
  59 }
  60
  61
  62 // Breit-Wigner with electromagnetic mass mixing
  63 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
  64                              EvtSpinType::spintype spin, double m0, double g0, NumType typeN,
  65                              double m0_mix, double g0_mix, double delta_mix, EvtComplex amp_mix)
  66   : _dp(dp),
  67     _pairAng(pairAng),
  68     _pairRes(pairRes),
  69     _spin(spin),
  70     _typeN(typeN),
  71     _m0(m0),_g0(g0),
  72     _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
  73     _m0_mix(m0_mix),_g0_mix(g0_mix),_delta_mix(delta_mix),_amp_mix(amp_mix),
  74     _g1(-1.),_g2(-1.),_coupling2(Undefined),
  75     _f_b(0.0), _f_d(1.5),
  76     _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
  77     _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.),_cutoff(-1.), _scaleByMOverQ(false),
  78     _alpha(0.)
  79 {
  80   _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
  81   _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
  82   _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
  83   _vd.set_f( 1.5 );
  84   // single BW (with electromagnetic mixing) cannot be K-matrix
  85   assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);
  86 }
  87
  88 // coupled Breit-Wigner
  89 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes,
  90                              EvtSpinType::spintype spin, double m0, NumType typeN, double g1, double g2, CouplingType coupling2)
  91   : _dp(dp),
  92     _pairAng(pairAng),
  93     _pairRes(pairRes),
  94     _spin(spin),
  95     _typeN(typeN),
  96     _m0(m0),_g0(-1.),
  97     _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
  98     _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
  99     _g1(g1),_g2(g2),_coupling2(coupling2),
 100     _f_b(0.0), _f_d(1.5),
 101     _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
 102     _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.),_cutoff(-1.), _scaleByMOverQ(false),
 103     _alpha(0.)
 104 {
 105   _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);
 106   _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
 107   _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
 108   _vd.set_f( 1.5 );
 109   assert(_coupling2 != Undefined);
 110   assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II); // coupled BW cannot be K-matrix
 111   assert(_typeN != LASS);     // coupled BW cannot be LASS
 112   assert(_typeN != NBW);      // for coupled BW, only relativistic BW
 113 }
 114
 115
 116 // K-Matrix (A&S)
 117 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes, std::string nameIndex, NumType typeN,
 118                              EvtComplex fr12prod, EvtComplex fr13prod, EvtComplex fr14prod, EvtComplex fr15prod, double s0prod)
 119   : _dp(dp),
 120     _pairRes(pairRes),
 121     _typeN(typeN),
 122     _m0(0.),_g0(0.),
 123     _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
 124     _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
 125     _g1(-1.),_g2(-1.),_coupling2(Undefined),
 126     _f_b(0.), _f_d(0.),
 127     _kmatrix_index(-1),_fr12prod(fr12prod),_fr13prod(fr13prod),_fr14prod(fr14prod),_fr15prod(fr15prod),_s0prod(s0prod),
 128     _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.),_cutoff(-1.), _scaleByMOverQ(false),
 129     _alpha(0.)
 130 {
 131   assert(_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II);
 132   _spin=EvtSpinType::SCALAR;
 133   if (nameIndex=="Pole1") _kmatrix_index=1;
 134   else if (nameIndex=="Pole2") _kmatrix_index=2;
 135   else if (nameIndex=="Pole3") _kmatrix_index=3;
 136   else if (nameIndex=="Pole4") _kmatrix_index=4;
 137   else if (nameIndex=="Pole5") _kmatrix_index=5;
 138   else if (nameIndex=="f11prod") _kmatrix_index=6;
 139   else assert(0);
 140 }
 141
 142
 143 // LASS parameterization
 144 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes,
 145                              double m0, double g0, double a, double r, double B, double phiB, double R, double phiR, double cutoff, bool scaleByMOverQ)
 146   : _dp(dp),
 147     _pairRes(pairRes),
 148     _typeN(LASS),
 149     _m0(m0),_g0(g0),
 150     _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
 151     _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
 152     _g1(-1.),_g2(-1.),_coupling2(Undefined),
 153     _f_b(0.0), _f_d(1.5),
 154     _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
 155     _a(a),_r(r),_Blass(B),_phiB(phiB),_R(R),_phiR(phiR), _cutoff(cutoff), _scaleByMOverQ(scaleByMOverQ),
 156     _alpha(0.)
 157 {
 158   _spin=EvtSpinType::SCALAR;
 159   _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
 160   _vd.set_f( 1.5 ); // Default values for Blatt-Weisskopf factors.
 161 }
 162
 163
 164 //Flatte
 165 EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, EvtCyclic3::Pair pairRes, double m0)
 166   : _dp(dp),
 167     _pairRes(pairRes),
 168     _typeN(FLATTE),
 169     _m0(m0), _g0(0.),
 170     _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
 171     _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
 172     _g1(-1.),_g2(-1.),_coupling2(Undefined),
 173     _f_b(0.), _f_d(0.),
 174     _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
 175     _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.),_cutoff(-1.), _scaleByMOverQ(false),
 176     _alpha(0.)
 177 {
 178   _spin=EvtSpinType::SCALAR;
 179 }
 180
 181
 182 EvtDalitzReso::EvtDalitzReso(const EvtDalitzReso& other)
 183   : _dp(other._dp),
 184     _pairAng(other._pairAng),
 185     _pairRes(other._pairRes),
 186     _spin(other._spin),
 187     _typeN(other._typeN),
 188     _m0(other._m0),_g0(other._g0),
 189     _vb(other._vb),_vd(other._vd),
 190     _massFirst(other._massFirst),_massSecond(other._massSecond),
 191     _m0_mix(other._m0_mix),_g0_mix(other._g0_mix),_delta_mix(other._delta_mix),_amp_mix(other._amp_mix),
 192     _g1(other._g1),_g2(other._g2),_coupling2(other._coupling2),
 193     _f_b(other._f_b), _f_d(other._f_d),
 194     _kmatrix_index(other._kmatrix_index),
 195     _fr12prod(other._fr12prod),_fr13prod(other._fr13prod),_fr14prod(other._fr14prod),_fr15prod(other._fr15prod),
 196     _s0prod(other._s0prod),
 197     _a(other._a),_r(other._r),_Blass(other._Blass),_phiB(other._phiB),_R(other._R),_phiR(other._phiR),_cutoff(other._cutoff), _scaleByMOverQ(other._scaleByMOverQ),
 198     _alpha(other._alpha),
 199     _flatteParams(other._flatteParams)
 200 {}
 201
 202
 203 EvtDalitzReso::~EvtDalitzReso()
 204 {}
 205
 206
 207 EvtComplex EvtDalitzReso::evaluate(const EvtDalitzPoint& x)
 208 {
 209   double m = sqrt(x.q(_pairRes));
 210
 211   if (_typeN==NON_RES)
 212     return EvtComplex(1.0,0.0);
 213
 214   if (_typeN==NON_RES_LIN)
 215     return m*m;
 216
 217   if (_typeN==NON_RES_EXP)
 218     return exp(-_alpha*m*m);
 219
 220   // do use always hash table (speed up fitting)
 221   if (_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II)
 222     return Fvector( m*m, _kmatrix_index );
 223
 224   if (_typeN==LASS)
 225     return lass(m*m);
 226
 227   if (_typeN==FLATTE)
 228     return flatte(m);
 229
 230   EvtComplex amp(1.0,0.0);
 231
 232   if (fabs(_dp.bigM() - x.bigM()) > 0.000001) {
 233     _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),x.bigM(),_spin);
 234     _vb.set_f(_f_b);
 235   }
 236   EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
 237   EvtTwoBodyKine vd(_massFirst,_massSecond,m);
 238
 239   EvtComplex prop(0,0);
 240   if (_typeN==NBW) {
 241     prop = propBreitWigner(_m0,_g0,m);
 242   } else if (_typeN==GAUSS_CLEO || _typeN==GAUSS_CLEO_ZEMACH) {
 243     prop = propGauss(_m0,_g0,m);
 244   } else {
 245     if (_coupling2==Undefined) {
 246       // single BW
 247       double g = (_g0<=0. || _vd.pD()<=0.)? -_g0 : _g0*_vd.widthFactor(vd);  // running width
 248       if (_typeN==GS_CLEO || _typeN==GS_CLEO_ZEMACH) {
 249         // Gounaris-Sakurai (GS)
 250         assert(_massFirst==_massSecond);
 251         prop = propGounarisSakurai(_m0,fabs(_g0),_vd.pD(),m,g,vd.p());
 252       } else {
 253         // standard relativistic BW
 254         prop = propBreitWignerRel(_m0,g,m);
 255       }
 256     } else {
 257       // coupled width BW
 258       EvtComplex G1,G2;
 259       switch (_coupling2) {
 260       case PicPic: {
 261         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 262         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
 263         G2 = _g2*_g2*psFactor(mPic,mPic,m);
 264         break;
 265       }
 266       case PizPiz: {
 267         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 268         static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
 269         G2 = _g2*_g2*psFactor(mPiz,mPiz,m);
 270         break;
 271       }
 272       case PiPi: {
 273         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 274         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
 275         static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
 276         G2 = _g2*_g2*psFactor(mPic,mPic,mPiz,mPiz,m);
 277         break;
 278       }
 279       case KcKc: {
 280         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 281         static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
 282         G2 = _g2*_g2*psFactor(mKc,mKc,m);
 283         break;
 284       }
 285       case KzKz: {
 286         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 287         static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
 288         G2 = _g2*_g2*psFactor(mKz,mKz,m);
 289         break;
 290       }
 291       case KK: {
 292         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 293         static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
 294         static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
 295         G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
 296         break;
 297       }
 298       case EtaPic: {
 299         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 300         static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
 301         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
 302         G2 = _g2*_g2*psFactor(mEta,mPic,m);
 303         break;
 304       }
 305       case EtaPiz: {
 306         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
 307         static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
 308         static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
 309         G2 = _g2*_g2*psFactor(mEta,mPiz,m);
 310         break;
 311       }
 312       case PicPicKK: {
 313         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
 314         //G1 = _g1*_g1*psFactor(mPic,mPic,m);
 315         G1 = _g1*psFactor(mPic,mPic,m);
 316         static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
 317         static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
 318         //G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
 319         G2 = _g2*psFactor(mKc,mKc,mKz,mKz,m);
 320         break;
 321       }
 322       default:
 323         std::cout << "EvtDalitzReso:evaluate(): PANIC, wrong coupling2 state." << std::endl;
 324         assert(0);
 325         break;
 326       }
 327       // calculate standard couple BW propagator
 328       if (_coupling2 != WA76)
 329         prop = _g1*propBreitWignerRelCoupled(_m0,G1,G2,m);
 330     }
 331   }
 332   amp *= prop;
 333
 334   // Compute form-factors (Blatt-Weisskopf penetration factor)
 335   amp *= _vb.formFactor(vb);
 336   amp *= _vd.formFactor(vd);
 337
 338   // Compute numerator (angular distribution)
 339   amp *= numerator(x,vb,vd);
 340
 341   // Compute electromagnetic mass mixing factor
 342   if (_m0_mix>0.) {
 343     EvtComplex prop_mix;
 344     if (_typeN==NBW) {
 345       prop_mix = propBreitWigner(_m0_mix,_g0_mix,m);
 346     } else {
 347       assert(_g1<0.); // running width only
 348       double g_mix = _g0_mix*_vd.widthFactor(vd);
 349       prop_mix = propBreitWignerRel(_m0_mix,g_mix,m);
 350     }
 351     amp *= mixFactor(prop,prop_mix);
 352   }
 353
 354   return amp;
 355 }
 356
 357
 358 EvtComplex EvtDalitzReso::psFactor(double & ma, double & mb, double& m)
 359 {
 360   if (m>(ma+mb)) {
 361     EvtTwoBodyKine vd(ma,mb,m);
 362     return EvtComplex(0,2*vd.p()/m);
 363   } else {
 364     // analytical continuation
 365     double s = m*m;
 366     double phaseFactor_analyticalCont = -0.5*(sqrt(4*ma*ma/s-1)+sqrt(4*mb*mb/s-1));
 367     return EvtComplex(phaseFactor_analyticalCont,0);
 368   }
 369 }
 370
 371
 372 EvtComplex EvtDalitzReso::psFactor(double & ma1,double & mb1, double & ma2, double & mb2, double& m)
 373 {
 374   return 0.5*(psFactor(ma1,mb1,m)+psFactor(ma2,mb2,m));
 375 }
 376
 377
 378 EvtComplex EvtDalitzReso::propGauss(const double& m0, const double& s0, const double& m)
 379 {
 380   // Gaussian
 381   double gauss = 1./sqrt(EvtConst::twoPi)/s0*exp(-(m-m0)*(m-m0)/2./(s0*s0));
 382   return EvtComplex(gauss,0.);
 383 }
 384
 385
 386 EvtComplex EvtDalitzReso::propBreitWigner(const double& m0, const double& g0, const double& m)
 387 {
 388   // non-relativistic BW
 389   return sqrt(g0/EvtConst::twoPi)/(m-m0-EvtComplex(0.0,g0/2.));
 390 }
 391
 392
 393 EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const double& g0, const double& m)
 394 {
 395   // relativistic BW with real width
 396   return 1./(m0*m0-m*m-EvtComplex(0.,m0*g0));
 397 }
 398
 399
 400
 401 EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const EvtComplex& g0, const double& m)
 402 {
 403   // relativistic BW with complex width
 404   return 1./(m0*m0-m*m-EvtComplex(0.,m0)*g0);
 405 }
 406
 407
 408 EvtComplex EvtDalitzReso::propBreitWignerRelCoupled(const double& m0, const EvtComplex& g1, const EvtComplex& g2, const double& m)
 409 {
 410   // relativistic coupled BW
 411   return 1./(m0*m0-m*m-(g1+g2));
 412 }
 413
 414 EvtComplex EvtDalitzReso::propGounarisSakurai(const double& m0, const double& g0, const double& k0,
 415                                             const double& m, const double& g, const double& k)
 416 {
 417   // Gounaris-Sakurai parameterization of pi+pi- P wave. PRD, Vol61, 112002. PRL, Vol21, 244.
 418   // Expressions taken from BAD637v4, after fixing the imaginary part of the BW denominator: i M_R Gamma_R(s) --> i sqrt(s) Gamma_R(s)
 419   return (1.+GS_d(m0,k0)*g0/m0)/(m0*m0-m*m-EvtComplex(0.,m*g)+GS_f(m0,g0,k0,m,k));
 420 }
 421
 422
 423 inline double EvtDalitzReso::GS_f(const double& m0, const double& g0, const double& k0, const double& m, const double& k)
 424 {
 425   // m: sqrt(s)
 426   // m0: nominal resonance mass
 427   // k: momentum of pion in resonance rest frame (at m)
 428   // k0: momentum of pion in resonance rest frame (at nominal resonance mass)
 429   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) );
 430 }
 431
 432 inline double EvtDalitzReso::GS_h(const double& m, const double& k)
 433 {return 2./EvtConst::pi*k/m*log((m+2.*k)/(2.*_massFirst)) ;}
 434
 435 inline double EvtDalitzReso::GS_dhods(const double& m0, const double& k0)
 436 {return GS_h(m0,k0)*( 0.125/(k0*k0) - 0.5/(m0*m0) ) + 0.5/(EvtConst::pi*m0*m0) ;}
 437
 438 inline double EvtDalitzReso::GS_d(const double& m0, const double& k0)
 439 {return 3./EvtConst::pi*_massFirst*_massFirst/(k0*k0)*log((m0+2.*k0)/(2.*_massFirst)) +
 440    m0/(2.*EvtConst::pi*k0) - _massFirst*_massFirst*m0/(EvtConst::pi*k0*k0*k0) ;}
 441
 442
 443 EvtComplex EvtDalitzReso::numerator(const EvtDalitzPoint& x, const EvtTwoBodyKine& vb, const EvtTwoBodyKine& vd)
 444 {
 445   EvtComplex ret(0.,0.);
 446
 447   // Non-relativistic Breit-Wigner
 448   if(NBW == _typeN) {
 449     ret = angDep(x);
 450   }
 451
 452   // Standard relativistic Zemach propagator
 453   else if(RBW_ZEMACH == _typeN) {
 454     ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
 455   }
 456
 457   // Standard relativistic Zemach propagator
 458   else if(RBW_ZEMACH2 == _typeN) {
 459     ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*_vb.phaseSpaceFactor(vb,EvtTwoBodyKine::AB)*angDep(x);
 460     if(_spin == EvtSpinType::VECTOR) {
 461       ret *= -4.;
 462     } else if(_spin == EvtSpinType::TENSOR) {
 463       ret *= 16./3.;
 464     } else if(_spin != EvtSpinType::SCALAR)
 465       assert(0);
 466   }
 467
 468   // Kuehn-Santamaria normalization:
 469   else if(RBW_KUEHN == _typeN) {
 470     ret = _m0*_m0 * angDep(x);
 471   }
 472
 473   // CLEO amplitude
 474   else if( ( RBW_CLEO        == _typeN ) || ( GS_CLEO           == _typeN ) ||
 475            ( RBW_CLEO_ZEMACH == _typeN ) || ( GS_CLEO_ZEMACH    == _typeN ) ||
 476            ( GAUSS_CLEO      == _typeN ) || ( GAUSS_CLEO_ZEMACH == _typeN)) {
 477
 478     Index iA = other(_pairAng);           // A = other(BC)
 479     Index iB = common(_pairRes,_pairAng); // B = common(AB,BC)
 480     Index iC = other(_pairRes);           // C = other(AB)
 481
 482     double M = x.bigM();
 483     double mA = x.m(iA);
 484     double mB = x.m(iB);
 485     double mC = x.m(iC);
 486     double qAB = x.q(combine(iA,iB));
 487     double qBC = x.q(combine(iB,iC));
 488     double qCA = x.q(combine(iC,iA));
 489
 490     double M2 = M*M;
 491     double m02 = ((RBW_CLEO_ZEMACH == _typeN)||(GS_CLEO_ZEMACH == _typeN)||(GAUSS_CLEO_ZEMACH == _typeN))?  qAB : _m0*_m0;
 492     double mA2 = mA*mA;
 493     double mB2 = mB*mB;
 494     double mC2 = mC*mC;
 495
 496     if (_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
 497     else if(_spin == EvtSpinType::VECTOR) {
 498       ret = qCA - qBC + (M2 - mC2)*(mB2 - mA2)/m02;
 499     } else if(_spin == EvtSpinType::TENSOR) {
 500       double x1 = qBC - qCA + (M2 - mC2)*(mA2 - mB2)/m02;
 501       double x2 = M2 - mC2;
 502       double x3 = qAB - 2*M2 - 2*mC2 + x2*x2/m02;
 503       double x4 = mA2 - mB2;
 504       double x5 = qAB - 2*mB2 - 2*mA2 + x4*x4/m02;
 505       ret = x1*x1 - x3*x5/3.;
 506     } else assert(0);
 507   }
 508
 509   return ret;
 510 }
 511
 512
 513 double EvtDalitzReso::angDep(const EvtDalitzPoint& x)
 514 {
 515   // Angular dependece for factorizable amplitudes
 516   // unphysical cosines indicate we are in big trouble
 517   double cosTh = x.cosTh(_pairAng,_pairRes);  // angle between common(reso,ang) and other(reso)
 518   if(fabs(cosTh) > 1.) {
 519     report(INFO,"EvtGen") << "cosTh " << cosTh << std::endl;
 520     assert(0);
 521   }
 522
 523   // in units of half-spin
 524   return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
 525 }
 526
 527
 528 EvtComplex EvtDalitzReso::mixFactor(EvtComplex prop, EvtComplex prop_mix)
 529 {
 530   double Delta = _delta_mix*(_m0+_m0_mix);
 531   return 1/(1-Delta*Delta*prop*prop_mix)*(1+_amp_mix*Delta*prop_mix);
 532 }
 533
 534
 535
 536 EvtComplex EvtDalitzReso::Fvector( double s, int index )
 537 {
 538   assert(index>=1 && index<=6);
 539
 540   //Define the complex coupling constant
 541   //The convection is as follow
 542   //i=0 --> pi+ pi-
 543   //i=1 --> KK
 544   //i=2 --> 4pi
 545   //i=3 --> eta eta
 546   //i=4 --> eta eta'
 547   //The first index is the resonace-pole index
 548
 549   double g[5][5]; // Coupling constants. The first index is the pole index. The second index is the decay channel
 550   double ma[5];   // Pole masses. The unit is in GeV
 551
 552   int solution = (_typeN==K_MATRIX)? 3 : (    (_typeN==K_MATRIX_I)? 1 : ( (_typeN==K_MATRIX_II)? 2 : 0 )    ) ;
 553   if (solution==0) { std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! " << std::endl; abort(); }
 554
 555   if (solution == 3 ) {
 556
 557     // coupling constants
 558     //pi+pi- channel
 559     g[0][0]=0.22889;
 560     g[1][0]=0.94128;
 561     g[2][0]=0.36856;
 562     g[3][0]=0.33650;
 563     g[4][0]=0.18171;
 564     //K+K- channel
 565     g[0][1]=-0.55377;
 566     g[1][1]=0.55095;
 567     g[2][1]=0.23888;
 568     g[3][1]=0.40907;
 569     g[4][1]=-0.17558;
 570     //4pi channel
 571     g[0][2]=0;
 572     g[1][2]=0;
 573     g[2][2]=0.55639;
 574     g[3][2]=0.85679;
 575     g[4][2]=-0.79658;
 576     //eta eta channel
 577     g[0][3]=-0.39899;
 578     g[1][3]=0.39065;
 579     g[2][3]=0.18340;
 580     g[3][3]=0.19906;
 581     g[4][3]=-0.00355;
 582     //eta eta' channel
 583     g[0][4]=-0.34639;
 584     g[1][4]=0.31503;
 585     g[2][4]=0.18681;
 586     g[3][4]=-0.00984;
 587     g[4][4]=0.22358;
 588
 589     // Pole masses
 590     ma[0]=0.651;
 591     ma[1]=1.20360;
 592     ma[2]=1.55817;
 593     ma[3]=1.21000;
 594     ma[4]=1.82206;
 595
 596   } else if (solution == 1) { // solnI.txt
 597
 598     // coupling constants
 599     //pi+pi- channel
 600     g[0][0]=0.31896;
 601     g[1][0]=0.85963;
 602     g[2][0]=0.47993;
 603     g[3][0]=0.45121;
 604     g[4][0]=0.39391;
 605     //K+K- channel
 606     g[0][1]=-0.49998;
 607     g[1][1]=0.52402;
 608     g[2][1]=0.40254;
 609     g[3][1]=0.42769;
 610     g[4][1]=-0.30860;
 611     //4pi channel
 612     g[0][2]=0;
 613     g[1][2]=0;
 614     g[2][2]=1.0;
 615     g[3][2]=1.15088;
 616     g[4][2]=0.33999;
 617     //eta eta channel
 618     g[0][3]=-0.21554;
 619     g[1][3]=0.38093;
 620     g[2][3]=0.21811;
 621     g[3][3]=0.22925;
 622     g[4][3]=0.06919;
 623     //eta eta' channel
 624     g[0][4]=-0.18294;
 625     g[1][4]=0.23788;
 626     g[2][4]=0.05454;
 627     g[3][4]=0.06444;
 628     g[4][4]=0.32620;
 629
 630     // Pole masses
 631     ma[0]=0.7369;
 632     ma[1]=1.24347;
 633     ma[2]=1.62681;
 634     ma[3]=1.21900;
 635     ma[4]=1.74932;
 636
 637   } else if (solution == 2) { // solnIIa.txt
 638
 639     // coupling constants
 640     //pi+pi- channel
 641     g[0][0]=0.26014;
 642     g[1][0]=0.95289;
 643     g[2][0]=0.46244;
 644     g[3][0]=0.41848;
 645     g[4][0]=0.01804;
 646     //K+K- channel
 647     g[0][1]=-0.57849;
 648     g[1][1]=0.55887;
 649     g[2][1]=0.31712;
 650     g[3][1]=0.49910;
 651     g[4][1]=-0.28430;
 652     //4pi channel
 653     g[0][2]=0;
 654     g[1][2]=0;
 655     g[2][2]=0.70340;
 656     g[3][2]=0.96819;
 657     g[4][2]=-0.90100;
 658     //eta eta channel
 659     g[0][3]=-0.32936;
 660     g[1][3]=0.39910;
 661     g[2][3]=0.22963;
 662     g[3][3]=0.24415;
 663     g[4][3]=-0.07252;
 664     //eta eta' channel
 665     g[0][4]=-0.30906;
 666     g[1][4]=0.31143;
 667     g[2][4]=0.19802;
 668     g[3][4]=-0.00522;
 669     g[4][4]=0.17097;
 670
 671     // Pole masses
 672     ma[0]=0.67460;
 673     ma[1]=1.21094;
 674     ma[2]=1.57896;
 675     ma[3]=1.21900;
 676     ma[4]=1.86602;
 677   }
 678
 679   //Now define the K-matrix pole
 680   double  rho1sq,rho2sq,rho4sq,rho5sq;
 681   EvtComplex rho[5];
 682   double f[5][5];
 683
 684   //Initalize the mass of the resonance
 685   double mpi=0.13957;
 686   double mK=0.493677;     //using charged K value
 687   double meta=0.54775;    //using PDG value
 688   double metap=0.95778;   //using PDG value
 689
 690   //Initialize the matrix to value zero
 691   EvtComplex K[5][5];
 692   for(int i=0;i<5;i++) {
 693     for(int j=0;j<5;j++) {
 694       K[i][j]=EvtComplex(0,0);
 695       f[i][j]=0;
 696     }
 697   }
 698
 699   //Input the _f[i][j] scattering data
 700   double s_scatt=0.0 ;
 701   if (solution == 3)
 702     s_scatt=-3.92637;
 703   else if (solution == 1)
 704     s_scatt= -5.0 ;
 705   else if (solution == 2)
 706     s_scatt= -5.0 ;
 707   double sa=1.0;
 708   double sa_0=-0.15;
 709   if (solution == 3) {
 710     f[0][0]=0.23399;  // f^scatt
 711     f[0][1]=0.15044;
 712     f[0][2]=-0.20545;
 713     f[0][3]=0.32825;
 714     f[0][4]=0.35412;
 715   }else if (solution == 1) {
 716     f[0][0]=0.04214;  // f^scatt
 717     f[0][1]=0.19865;
 718     f[0][2]=-0.63764;
 719     f[0][3]=0.44063;
 720     f[0][4]=0.36717;
 721   }else if (solution == 2) {
 722     f[0][0]=0.26447;  // f^scatt
 723     f[0][1]=0.10400;
 724     f[0][2]=-0.35445;
 725     f[0][3]=0.31596;
 726     f[0][4]=0.42483;
 727   }
 728   f[1][0]=f[0][1];
 729   f[2][0]=f[0][2];
 730   f[3][0]=f[0][3];
 731   f[4][0]=f[0][4];
 732
 733   //Now construct the phase-space factor
 734   //For eta-eta' there is no difference term
 735   rho1sq = 1. - pow( mpi + mpi, 2 ) / s;   //pi+ pi- phase factor
 736   if( rho1sq >= 0 )
 737     rho[ 0 ] = EvtComplex( sqrt( rho1sq ), 0 );
 738   else
 739     rho[ 0 ] = EvtComplex( 0, sqrt( -rho1sq ) );
 740
 741   rho2sq = 1. - pow( mK + mK, 2 ) / s;
 742   if( rho2sq >= 0 )
 743     rho[ 1 ] = EvtComplex( sqrt( rho2sq ), 0 );
 744   else
 745     rho[ 1 ] = EvtComplex( 0, sqrt( -rho2sq ) );
 746
 747   //using the A&S 4pi phase space Factor:
 748   //Shit, not continue
 749   if( s <= 1 )
 750     {
 751       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;
 752       double cont32 = sqrt(1.0-(16.0*mpi*mpi));
 753       rho[ 2 ] = EvtComplex( cont32 * real, 0 );
 754     }
 755   else
 756     rho[ 2 ] = EvtComplex( sqrt( 1. - 16. * mpi * mpi / s ), 0 );
 757
 758   rho4sq = 1. - pow( meta + meta, 2 ) / s;
 759   if( rho4sq >= 0 )
 760     rho[ 3 ] = EvtComplex( sqrt( rho4sq ), 0 );
 761   else
 762     rho[ 3 ] = EvtComplex( 0, sqrt( -rho4sq ) );
 763
 764   rho5sq = 1. - pow( meta + metap, 2 ) / s;
 765   if( rho5sq >= 0 )
 766     rho[ 4 ] = EvtComplex( sqrt( rho5sq ), 0 );
 767   else
 768     rho[ 4 ] = EvtComplex( 0, sqrt( -rho5sq ) );
 769
 770   double smallTerm = 1; // Factor to prevent divergences.
 771
 772   // Check if some pole may arise problems.
 773   for ( int pole = 0; pole < 5; pole++ )
 774     if ( fabs( pow( ma[ pole ], 2 ) - s ) < PRECISION )
 775       smallTerm = pow( ma[ pole ], 2 ) - s;
 776
 777   //now sum all the pole
 778   //equation (3) in the E791 K-matrix paper
 779   for(int i=0;i<5;i++) {
 780     for(int j=0;j<5;j++) {
 781       for (int pole_index=0;pole_index<5;pole_index++) {
 782         double A=g[pole_index][i]*g[pole_index][j];
 783         double B=ma[pole_index]*ma[pole_index]-s;
 784
 785         if ( fabs( B ) < PRECISION )
 786           K[ i ][ j ] += EvtComplex( A    , 0 );
 787         else
 788           K[ i ][ j ] += EvtComplex( A / B, 0 ) * smallTerm;
 789       }
 790     }
 791   }
 792
 793   //now add the SVT part
 794   for(int i=0;i<5;i++) {
 795     for(int j=0;j<5;j++) {
 796       double C=f[i][j]*(1.0-s_scatt);
 797       double D=(s-s_scatt);
 798       K[ i ][ j ] += EvtComplex( C / D, 0 ) * smallTerm;
 799     }
 800   }
 801
 802   //Fix the bug in the FOCUS paper
 803   //Include the Alder zero term:
 804   for(int i=0;i<5;i++) {
 805     for(int j=0;j<5;j++) {
 806       double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
 807       double F=(s-sa_0);
 808       K[ i ][ j ] *= EvtComplex(E/F,0);
 809     }
 810   }
 811
 812   //This is not correct!
 813   //(1-ipK) != (1-iKp)
 814   static EvtMatrix< EvtComplex > mat;
 815   mat.setRange( 5 ); // Try to do in only the first time. DEFINE ALLOCATION IN CONSTRUCTOR.
 816
 817   for ( int row = 0; row < 5; row++ )
 818     for ( int col = 0; col < 5; col++ )
 819       mat( row, col ) = ( row == col ) * smallTerm - EvtComplex( 0., 1. ) * K[ row ][ col ] * rho[ col ];
 820
 821
 822   EvtMatrix< EvtComplex >* matInverse = mat.inverse();  //The 1st row of the inverse matrix. This matrix is {(I-iKp)^-1}_0j
 823   vector< EvtComplex > U1j;
 824   for ( int j = 0; j < 5; j++ )
 825     U1j.push_back( (*matInverse)[ 0 ][ j ] );
 826
 827   delete matInverse;
 828
 829   //this calculates final F0 factor
 830   EvtComplex value( 0, 0 );
 831   if (index<=5) {
 832     //this calculates the beta_idx Factors
 833     for(int j=0;j<5;j++) {        // sum for 5 channel
 834       EvtComplex top    = U1j[j]*g[index-1][j];
 835       double     bottom = ma[index-1]*ma[index-1]-s;
 836
 837       if ( fabs( bottom ) < PRECISION )
 838         value += top;
 839       else
 840         value += top / bottom * smallTerm;
 841     }
 842   } else {
 843     //this calculates fprod Factors
 844     value += U1j[0];
 845     value += U1j[1]*_fr12prod;
 846     value += U1j[2]*_fr13prod;
 847     value += U1j[3]*_fr14prod;
 848     value += U1j[4]*_fr15prod;
 849
 850     value *= (1-_s0prod)/(s-_s0prod) * smallTerm;
 851   }
 852
 853   return value;
 854 }
 855
 856
 857 //replace Breit-Wigner with LASS
 858 EvtComplex EvtDalitzReso::lass(double s)
 859 {
 860   EvtTwoBodyKine vd(_massFirst,_massSecond, sqrt(s));
 861   double q = vd.p();
 862   double GammaM = _g0*_vd.widthFactor(vd);  // running width;
 863
 864   //calculate the background phase motion
 865   double cot_deltaB = 1.0/(_a*q) + 0.5*_r*q;
 866   double deltaB = atan( 1.0/cot_deltaB);
 867   double totalB = deltaB + _phiB ;
 868
 869   //calculate the resonant phase motion
 870   double deltaR = atan((_m0*GammaM/(_m0*_m0 - s)));
 871   double totalR = deltaR + _phiR ;
 872
 873   //sum them up
 874   EvtComplex  bkgB,resT;
 875   bkgB = EvtComplex(_Blass*sin(totalB),0)*EvtComplex(cos(totalB),sin(totalB));
 876   resT = EvtComplex(_R*sin(deltaR),0)*EvtComplex(cos(totalR),sin(totalR))*EvtComplex(cos(2*totalB),sin(2*totalB));
 877
 878   EvtComplex T;
 879   if(_cutoff>0 && sqrt(s)>_cutoff) T = resT;
 880   else T = bkgB + resT;
 881
 882   if(_scaleByMOverQ) T*=(sqrt(s)/q);
 883
 884   return T;
 885 }
 886
 887
 888 EvtComplex EvtDalitzReso::flatte(const double& m) {
 889
 890   EvtComplex w;
 891
 892   for (vector<EvtFlatteParam>::const_iterator param = _flatteParams.begin();
 893        param != _flatteParams.end();
 894        ++param) {
 895     double m1 = (*param).m1(); double m2 = (*param).m2();
 896     double g = (*param).g();
 897     w += (g*g*sqrtCplx((1-((m1-m2)*(m1-m2))/(m*m))*(1-((m1+m2)*(m1+m2))/(m*m))));
 898   }
 899
 900   EvtComplex denom = _m0*_m0 - m*m - EvtComplex(0,1)*w;
 901
 902   return EvtComplex(1.0,0.0)/denom;
 903 }