00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <stdlib.h>
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtCPUtil.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenModels/EvtBTo3piCP.hh"
00029 #include "EvtGenBase/EvtId.hh"
00030 #include <string>
00031 #include "EvtGenBase/EvtConst.hh"
00032
00033 extern "C" {
00034 extern void evt3pions_(double *,int *,double *,
00035 double *,double *,double *,double *,
00036 double *,double *,double *);
00037 }
00038
00039
00040 EvtBTo3piCP::~EvtBTo3piCP() {}
00041
00042
00043 void EvtBTo3piCP::getName(std::string& model_name){
00044
00045 model_name="BTO3PI_CP";
00046
00047 }
00048
00049
00050 EvtDecayBase* EvtBTo3piCP::clone(){
00051
00052 return new EvtBTo3piCP;
00053
00054 }
00055
00056 void EvtBTo3piCP::init(){
00057
00058
00059 checkNArg(2);
00060 checkNDaug(3);
00061
00062 checkSpinParent(EvtSpinType::SCALAR);
00063
00064 checkSpinDaughter(0,EvtSpinType::SCALAR);
00065 checkSpinDaughter(1,EvtSpinType::SCALAR);
00066 checkSpinDaughter(2,EvtSpinType::SCALAR);
00067 }
00068
00069
00070
00071 void EvtBTo3piCP::initProbMax(){
00072
00073 setProbMax(1.5);
00074
00075 }
00076
00077 void EvtBTo3piCP::decay( EvtParticle *p){
00078
00079
00080 static EvtId B0=EvtPDL::getId("B0");
00081 static EvtId B0B=EvtPDL::getId("anti-B0");
00082
00083 double t;
00084 EvtId other_b;
00085
00086 EvtCPUtil::OtherB(p,t,other_b);
00087
00088 EvtParticle *pip,*pim,*pi0;
00089
00090 p->makeDaughters(getNDaug(),getDaugs());
00091
00092
00093 pip=p->getDaug(0);
00094 pim=p->getDaug(1);
00095 pi0=p->getDaug(2);
00096
00097 EvtVector4R p4[3];
00098
00099 double dm=getArg(0);
00100 double alpha=getArg(1);
00101 int iset;
00102
00103 static int first=1;
00104
00105 if (first==1) {
00106 iset=10000;
00107 first=0;
00108 }
00109 else{
00110 iset=0;
00111 }
00112
00113 double p4piplus[4],p4piminus[4],p4gamm1[4],p4gamm2[4];
00114
00115 double realA,imgA,realbarA,imgbarA;
00116
00117 evt3pions_(&alpha,&iset,p4piplus,p4piminus,p4gamm1,p4gamm2,
00118 &realA,&imgA,&realbarA,&imgbarA);
00119
00120 p4[0].set(p4piplus[3],p4piplus[0],p4piplus[1],p4piplus[2]);
00121 p4[1].set(p4piminus[3],p4piminus[0],p4piminus[1],p4piminus[2]);
00122 p4[2].set(p4gamm1[3]+p4gamm2[3],p4gamm1[0]+p4gamm2[0],
00123 p4gamm1[1]+p4gamm2[1],p4gamm1[2]+p4gamm2[2]);
00124
00125 if (pip->getId()==EvtPDL::getId("pi+")) {
00126 pip->init( getDaug(0), p4[0] );
00127 pim->init( getDaug(1), p4[1] );
00128 }
00129 else {
00130 pip->init( getDaug(0), p4[1] );
00131 pim->init( getDaug(1), p4[0] );
00132 }
00133
00134 pi0->init( getDaug(2), p4[2] );
00135
00136 EvtComplex amp;
00137
00138 EvtComplex A(realA,imgA);
00139 EvtComplex Abar(realbarA,imgbarA);
00140
00141 if (other_b==B0B){
00142 amp=A*cos(dm*t/(2*EvtConst::c))+
00143 EvtComplex(0.,1.)*Abar*sin(dm*t/(2*EvtConst::c));
00144 }
00145 if (other_b==B0){
00146 amp=Abar*cos(dm*t/(2*EvtConst::c))+
00147 EvtComplex(0.,1.)*A*sin(dm*t/(2*EvtConst::c));
00148 }
00149
00150 vertex(amp);
00151
00152 return ;
00153 }
00154