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 <stdlib.h>
00023 #include "EvtGenBase/EvtParticle.hh"
00024 #include "EvtGenBase/EvtGenKine.hh"
00025 #include "EvtGenBase/EvtCPUtil.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenModels/EvtBTo4piCP.hh"
00029 #include "EvtGenBase/EvtId.hh"
00030 #include "EvtGenBase/EvtConst.hh"
00031 #include <string>
00032
00033 EvtBTo4piCP::~EvtBTo4piCP() {}
00034
00035
00036 EvtComplex EvtAmpA2(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
00037 const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
00038
00039
00040 static EvtId A2M=EvtPDL::getId("a_2-");
00041 static EvtId RHO0=EvtPDL::getId("rho0");
00042
00043 EvtVector4R p4a2,p4rho,p4b;
00044
00045 p4rho=p4pi1+p4pi2;
00046
00047 p4a2=p4rho+p4pi3;
00048
00049 p4b=p4a2+p4pi4;
00050
00051 EvtVector4R p4b_a2,p4rho_a2,p4pi1_a2,p4a2_a2;
00052
00053 p4b_a2=boostTo(p4b,p4a2);
00054 p4rho_a2=boostTo(p4rho,p4a2);
00055 p4pi1_a2=boostTo(p4pi1,p4a2);
00056 p4a2_a2=boostTo(p4a2,p4a2);
00057
00058 EvtVector4R p4pi1_rho;
00059
00060 p4pi1_rho=boostTo(p4pi1_a2,p4rho_a2);
00061
00062 EvtVector4R vb,vrho,vpi,t;
00063
00064 vb=p4b_a2/p4b_a2.d3mag();
00065 vrho=p4rho_a2/p4rho_a2.d3mag();
00066 vpi=p4pi1_rho/p4pi1_rho.d3mag();
00067
00068 t.set(1.0,0.0,0.0,0.0);
00069
00070
00071 EvtComplex amp_a2;
00072
00073
00074
00075 double bwm_a2=EvtPDL::getMeanMass(A2M);
00076 double gamma_a2=EvtPDL::getWidth(A2M);
00077 double bwm_rho=EvtPDL::getMeanMass(RHO0);
00078 double gamma_rho=EvtPDL::getWidth(RHO0);
00079
00080 amp_a2=(sqrt(gamma_a2/EvtConst::twoPi)/
00081 ((p4a2).mass()-bwm_a2-EvtComplex(0.0,0.5*gamma_a2)))*
00082 (sqrt(gamma_rho/EvtConst::twoPi)/
00083 ((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
00084
00085 return amp_a2*
00086 (vb.get(1)*vrho.get(1)+vb.get(2)*vrho.get(2)+vb.get(3)*vrho.get(3))*
00087 (
00088 vpi.get(1)*(vb.get(2)*vrho.get(3)-vb.get(3)*vrho.get(2))+
00089 vpi.get(2)*(vb.get(3)*vrho.get(1)-vb.get(1)*vrho.get(3))+
00090 vpi.get(3)*(vb.get(1)*vrho.get(2)-vb.get(2)*vrho.get(1))
00091 );
00092
00093 }
00094
00095 EvtComplex EvtAmpA1(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
00096 const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
00097
00098
00099 static EvtId A1M=EvtPDL::getId("a_1-");
00100 static EvtId RHO0=EvtPDL::getId("rho0");
00101
00102 EvtVector4R p4a1,p4rho,p4b;
00103
00104 p4rho=p4pi1+p4pi2;
00105
00106 p4a1=p4rho+p4pi3;
00107
00108 p4b=p4a1+p4pi4;
00109
00110 EvtVector4R p4b_a1,p4rho_a1,p4pi1_a1,p4a1_a1;
00111
00112 p4b_a1=boostTo(p4b,p4a1);
00113 p4rho_a1=boostTo(p4rho,p4a1);
00114 p4pi1_a1=boostTo(p4pi1,p4a1);
00115 p4a1_a1=boostTo(p4a1,p4a1);
00116
00117 EvtVector4R p4pi1_rho;
00118
00119 p4pi1_rho=boostTo(p4pi1_a1,p4rho_a1);
00120
00121 EvtVector4R vb,vrho,vpi,t;
00122
00123 vb=p4b_a1/p4b_a1.d3mag();
00124 vrho=p4rho_a1/p4rho_a1.d3mag();
00125 vpi=p4pi1_rho/p4pi1_rho.d3mag();
00126
00127 t.set(1.0,0.0,0.0,0.0);
00128
00129 EvtComplex amp_a1;
00130
00131 double bwm_a1=EvtPDL::getMeanMass(A1M);
00132 double gamma_a1=EvtPDL::getWidth(A1M);
00133
00134
00135 double bwm_rho=EvtPDL::getMeanMass(RHO0);
00136 double gamma_rho=EvtPDL::getWidth(RHO0);
00137
00138 amp_a1=(sqrt(gamma_a1/EvtConst::twoPi)/
00139 ((p4a1).mass()-bwm_a1-EvtComplex(0.0,0.5*gamma_a1)))*
00140 (sqrt(gamma_rho/EvtConst::twoPi)/
00141 ((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
00142
00143 return amp_a1*
00144 (vb.get(1)*vpi.get(1)+vb.get(2)*vpi.get(2)+vb.get(3)*vpi.get(3));
00145
00146 }
00147
00148
00149 void EvtBTo4piCP::getName(std::string& model_name){
00150
00151 model_name="BTO4PI_CP";
00152
00153 }
00154
00155
00156 EvtDecayBase* EvtBTo4piCP::clone(){
00157
00158 return new EvtBTo4piCP;
00159
00160 }
00161
00162 void EvtBTo4piCP::init(){
00163
00164
00165 checkNArg(18);
00166 checkNDaug(4);
00167
00168 checkSpinParent(EvtSpinType::SCALAR);
00169
00170 checkSpinDaughter(0,EvtSpinType::SCALAR);
00171 checkSpinDaughter(1,EvtSpinType::SCALAR);
00172 checkSpinDaughter(2,EvtSpinType::SCALAR);
00173 checkSpinDaughter(3,EvtSpinType::SCALAR);
00174 }
00175
00176 void EvtBTo4piCP::decay( EvtParticle *p){
00177
00178
00179 static EvtId B0=EvtPDL::getId("B0");
00180 static EvtId B0B=EvtPDL::getId("anti-B0");
00181
00182
00183 double t;
00184 EvtId other_b;
00185
00186 EvtCPUtil::OtherB(p,t,other_b);
00187
00188 p->initializePhaseSpace(getNDaug(),getDaugs());
00189 EvtVector4R mom1 = p->getDaug(0)->getP4();
00190 EvtVector4R mom2 = p->getDaug(1)->getP4();
00191 EvtVector4R mom3 = p->getDaug(2)->getP4();
00192 EvtVector4R mom4 = p->getDaug(3)->getP4();
00193
00194
00195
00196
00197 EvtComplex amp;
00198
00199 EvtComplex A,Abar;
00200
00201
00202 EvtComplex A_a1p,Abar_a1p,A_a2p,Abar_a2p;
00203 EvtComplex A_a1m,Abar_a1m,A_a2m,Abar_a2m;
00204
00205 A_a1p=EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
00206 Abar_a1p=EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
00207
00208 A_a2p=EvtComplex(getArg(6)*cos(getArg(7)),getArg(6)*sin(getArg(7)));
00209 Abar_a2p=EvtComplex(getArg(8)*cos(getArg(9)),getArg(8)*sin(getArg(9)));
00210
00211 A_a1m=EvtComplex(getArg(10)*cos(getArg(11)),getArg(10)*sin(getArg(11)));
00212 Abar_a1m=EvtComplex(getArg(12)*cos(getArg(13)),getArg(12)*sin(getArg(13)));
00213
00214 A_a2m=EvtComplex(getArg(14)*cos(getArg(15)),getArg(14)*sin(getArg(15)));
00215 Abar_a2m=EvtComplex(getArg(16)*cos(getArg(17)),getArg(16)*sin(getArg(17)));
00216
00217 EvtComplex a2p_amp=EvtAmpA2(mom1,mom2,mom3,mom4)+
00218 EvtAmpA2(mom1,mom4,mom3,mom2)+
00219 EvtAmpA2(mom3,mom2,mom1,mom4)+
00220 EvtAmpA2(mom3,mom4,mom1,mom2);
00221
00222 EvtComplex a2m_amp=EvtAmpA2(mom2,mom3,mom4,mom1)+
00223 EvtAmpA2(mom2,mom1,mom4,mom3)+
00224 EvtAmpA2(mom4,mom3,mom2,mom1)+
00225 EvtAmpA2(mom4,mom1,mom2,mom3);
00226
00227 EvtComplex a1p_amp=EvtAmpA1(mom1,mom2,mom3,mom4)+
00228 EvtAmpA1(mom1,mom4,mom3,mom2)+
00229 EvtAmpA1(mom3,mom2,mom1,mom4)+
00230 EvtAmpA1(mom3,mom4,mom1,mom2);
00231
00232 EvtComplex a1m_amp=EvtAmpA1(mom2,mom3,mom4,mom1)+
00233 EvtAmpA1(mom2,mom1,mom4,mom3)+
00234 EvtAmpA1(mom4,mom3,mom2,mom1)+
00235 EvtAmpA1(mom4,mom1,mom2,mom3);
00236
00237
00238 A=A_a2p*a2p_amp+A_a1p*a1p_amp+
00239 A_a2m*a2m_amp+A_a1m*a1m_amp;
00240 Abar=Abar_a2p*a2p_amp+Abar_a1p*a1p_amp+
00241 Abar_a2m*a2m_amp+Abar_a1m*a1m_amp;
00242
00243
00244 if (other_b==B0B){
00245 amp=A*cos(getArg(1)*t/(2*EvtConst::c))+
00246 EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
00247 getArg(2)*EvtComplex(0.0,1.0)*Abar*sin(getArg(1)*t/(2*EvtConst::c));
00248 }
00249 if (other_b==B0){
00250 amp=A*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
00251 EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
00252 getArg(2)*Abar*cos(getArg(1)*t/(2*EvtConst::c));
00253 }
00254
00255 vertex(amp);
00256
00257 return ;
00258 }
00259