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