TEvtGen/EvtGenModels/EvtPhiDalitz.cxx

   1 #include "EvtGenBase/EvtPatches.hh"
   2
   3 #include <stdlib.h>
   4 #include <math.h>
   5 #include "EvtGenBase/EvtVector4R.hh"
   6 #include "EvtGenBase/EvtParticle.hh"
   7 #include "EvtGenBase/EvtGenKine.hh"
   8 #include "EvtGenBase/EvtPDL.hh"
   9 #include "EvtGenModels/EvtPhiDalitz.hh"
  10 #include "EvtGenBase/EvtReport.hh"
  11 #include <string>
  12
  13 // Implementation of KLOE measurement
  14 // PL B561: 55-60 (2003) + Erratum B609:449-450 (2005)
  15 // or hep-ex/0303016v2
  16
  17
  18 EvtPhiDalitz::~EvtPhiDalitz() {}
  19
  20 std::string EvtPhiDalitz::getName(){
  21
  22   return "PHI_DALITZ";
  23
  24 }
  25
  26
  27 EvtDecayBase* EvtPhiDalitz::clone(){
  28
  29   return new EvtPhiDalitz;
  30
  31 }
  32
  33 void EvtPhiDalitz::init(){
  34
  35   // check that there are 0 arguments
  36   checkNArg(0);
  37   checkNDaug(3);
  38
  39   checkSpinParent(EvtSpinType::VECTOR);
  40
  41   checkSpinDaughter(0,EvtSpinType::SCALAR);
  42   checkSpinDaughter(1,EvtSpinType::SCALAR);
  43   checkSpinDaughter(2,EvtSpinType::SCALAR);
  44
  45   _mRho=0.7758;
  46   _gRho=0.1439;
  47   _aD=0.78;
  48   _phiD=-2.47;
  49   _aOmega=0.0071;
  50   _phiOmega=-0.22;
  51
  52   _locPip=-1;
  53   _locPim=-1;
  54   _locPi0=-1;
  55
  56   for ( int i=0; i<3; i++) {
  57      if ( getDaug(i) == EvtPDL::getId("pi+")) _locPip=i;
  58      if ( getDaug(i) == EvtPDL::getId("pi-")) _locPim=i;
  59      if ( getDaug(i) == EvtPDL::getId("pi0")) _locPi0=i;
  60   }
  61   if ( _locPip == -1 || _locPim == -1 || _locPi0 == -1 ) {
  62     report(ERROR,"EvtGen") << getModelName() << "generator expects daughters to be pi+ pi- pi0\n";
  63     report(ERROR,"EvtGen") << "Found " << EvtPDL::name(getDaug(0)) << " "
  64                            << EvtPDL::name(getDaug(1)) << " "
  65                            << EvtPDL::name(getDaug(2)) << std::endl;
  66
  67   }
  68
  69
  70 }
  71
  72
  73
  74
  75 void EvtPhiDalitz::decay( EvtParticle *p){
  76
  77   EvtId PIP=EvtPDL::getId("pi+");
  78   EvtId PIM=EvtPDL::getId("pi-");
  79   EvtId PIZ=EvtPDL::getId("pi0");
  80   EvtId OMEGA=EvtPDL::getId("omega");
  81
  82   p->initializePhaseSpace(getNDaug(),getDaugs());
  83
  84   EvtVector4R Ppip = p->getDaug(_locPip)->getP4();
  85   EvtVector4R Ppim = p->getDaug(_locPim)->getP4();
  86   EvtVector4R Ppi0 = p->getDaug(_locPi0)->getP4();
  87   EvtVector4R Qp = (Ppim + Ppi0);
  88   EvtVector4R Qm = (Ppip + Ppi0);
  89   EvtVector4R Q0 = (Ppip + Ppim);
  90   double m2_pip = pow(EvtPDL::getMeanMass(PIP),2);
  91   double m2_pim = pow(EvtPDL::getMeanMass(PIM),2);
  92   double m2_pi0 = pow(EvtPDL::getMeanMass(PIZ),2);
  93   double M2rhop = pow(_mRho,2);
  94   double M2rhom = pow(_mRho,2);
  95   double M2rho0 = pow(_mRho,2);
  96   double M2omega = pow(EvtPDL::getMeanMass(OMEGA),2);
  97
  98   double Wrhop = _gRho;
  99   double Wrhom = _gRho;
 100   double Wrho0 = _gRho;
 101   double Womega = EvtPDL::getWidth(OMEGA);
 102
 103   EvtComplex Atot(0,0);
 104
 105   //Rho+ Risonance Amplitude
 106   double Gp = Wrhop*pow(((Qp.mass2()-m2_pim-m2_pi0)/2-M2rhop/4)/(M2rhop/4-(m2_pim+m2_pi0)/2),3/2)*(M2rhop/Qp.mass2());
 107   EvtComplex Drhop((Qp.mass2()-M2rhop),Qp.mass()*Gp);
 108   EvtComplex A1(M2rhop/Drhop);
 109
 110   //Rho- Risonance Amplitude
 111   double Gm = Wrhom*pow(((Qm.mass2()-m2_pip-m2_pi0)/2-M2rhom/4)/(M2rhom/4-(m2_pip+m2_pi0)/2),3/2)*(M2rhom/Qm.mass2());
 112   EvtComplex Drhom((Qm.mass2()-M2rhom),Qm.mass()*Gm);
 113   EvtComplex A2(M2rhom/Drhom);
 114
 115   //Rho0 Risonance Amplitude
 116   double G0 = Wrho0*pow(((Q0.mass2()-m2_pip-m2_pim)/2-M2rho0/4)/(M2rho0/4-(m2_pip+m2_pim)/2),3/2)*(M2rho0/Q0.mass2());
 117   EvtComplex Drho0((Q0.mass2()-M2rho0),Q0.mass()*G0);
 118   EvtComplex A3(M2rho0/Drho0);
 119
 120   //Omega Risonance Amplitude
 121   EvtComplex OmegaPhase(0,_phiOmega);
 122   EvtComplex DOmega((Q0.mass2()-M2omega),Q0.mass()*Womega);
 123   EvtComplex A4(_aOmega*M2omega*exp(OmegaPhase)/DOmega);
 124
 125   //Direct Decay Amplitude
 126   EvtComplex DirPhase(0,_phiD);
 127   EvtComplex A5(_aD*exp(DirPhase));
 128
 129   Atot=A1+A2+A3+A4+A5;
 130
 131   vertex(0,Atot);
 132   vertex(1,Atot);
 133   vertex(2,Atot);
 134
 135   return ;
 136
 137 }
 138
 139
 140
 141
 142
 143
 144