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 "EvtGenModels/EvtISGW2FF.hh"
00025 #include "EvtGenBase/EvtGenKine.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenModels/EvtISGW2FF.hh"
00028 #include "EvtGenBase/EvtTensor4C.hh"
00029 #include "EvtGenBase/EvtVector4C.hh"
00030 #include "EvtGenBase/EvtVector4R.hh"
00031 #include "EvtGenModels/EvtBHadronic.hh"
00032 #include "EvtGenBase/EvtReport.hh"
00033 #include <string>
00034 using std::endl;
00035
00036 EvtBHadronic::~EvtBHadronic() {}
00037
00038 void EvtBHadronic::getName(std::string& model_name){
00039
00040 model_name="BHADRONIC";
00041
00042 }
00043
00044
00045 EvtDecayBase* EvtBHadronic::clone(){
00046
00047 return new EvtBHadronic;
00048
00049 }
00050
00051
00052 void EvtBHadronic::init(){
00053
00054
00055 checkNArg(2);
00056
00057
00058 }
00059
00060 void EvtBHadronic::decay( EvtParticle *p){
00061
00062
00063 static EvtId B0=EvtPDL::getId("B0");
00064 static EvtId D0=EvtPDL::getId("D0");
00065 static EvtId DST0=EvtPDL::getId("D*0");
00066 static EvtId D3P10=EvtPDL::getId("D'_10");
00067 static EvtId D3P20=EvtPDL::getId("D_2*0");
00068 static EvtId D3P00=EvtPDL::getId("D_0*0");
00069 static EvtId D1P10=EvtPDL::getId("D_10");
00070
00071 static EvtISGW2FF ffmodel;
00072
00073 p->initializePhaseSpace(getNDaug(),getDaugs());
00074
00075
00076 EvtVector4R p4[MAX_DAUG];
00077 double m;
00078
00079 m = p->mass();
00080
00081
00082 int i,j;
00083
00084 for ( i=0; i<getNDaug(); i++) {
00085 p4[i]=p->getDaug(i)->getP4();
00086 }
00087
00088 int bcurrent,wcurrent;
00089 int nbcurrent=0;
00090 int nwcurrent=0;
00091
00092 bcurrent=(int)getArg(0);
00093 wcurrent=(int)getArg(1);
00094
00095 EvtVector4C jb[5],jw[5];
00096 EvtTensor4C g,tds;
00097
00098 EvtVector4R p4b;
00099 p4b.set(p->mass(),0.0,0.0,0.0);
00100
00101 EvtVector4R q;
00102 double q2;
00103
00104 EvtComplex ep_meson_bb[5];
00105 double f,gf,ap,am;
00106 double fp,fm;
00107 double hf,kf,bp,bm;
00108 EvtVector4R pp,pm;
00109 EvtVector4C ep_meson_b[5];
00110
00111 switch (bcurrent) {
00112
00113
00114 case 1:
00115 q=p4b-p4[0];
00116 q2=q*q;
00117 nbcurrent=1;
00118 ffmodel.getscalarff(B0,D0,EvtPDL::getMeanMass(D0),
00119 EvtPDL::getMeanMass(getDaugs()[1]),&fp,&fm);
00120 jb[0]=EvtVector4C(fp*(p4b+p4[0])-fm*q);
00121 break;
00122
00123 case 2:
00124 q=p4b-p4[0];
00125 q2=q*q;
00126 nbcurrent=3;
00127 ffmodel.getvectorff(B0,DST0,EvtPDL::getMeanMass(DST0),q2,&f,&gf,&ap,&am);
00128
00129 g.setdiag(1.0,-1.0,-1.0,-1.0);
00130 tds = -f*g
00131 -ap*(directProd(p4b,p4b)+directProd(p4b,p4[0]))
00132 -gf*EvtComplex(0.0,1.0)*dual(directProd(p4[0]+p4b,p4b-p4[0]))
00133 -am*((directProd(p4b,p4b)-directProd(p4b,p4[0])));
00134 jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
00135 jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
00136 jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
00137 break;
00138
00139 case 3:
00140 q=p4b-p4[0];
00141 q2=q*q;
00142 nbcurrent=3;
00143 ffmodel.getvectorff(B0,D3P10,EvtPDL::getMeanMass(D3P10),q2,&f,&gf,&ap,&am);
00144
00145 g.setdiag(1.0,-1.0,-1.0,-1.0);
00146 tds = -f*g
00147 -ap*(directProd(p4b,p4b)+directProd(p4b,p4[0]))
00148 -gf*EvtComplex(0.0,1.0)*dual(directProd(p4[0]+p4b,p4b-p4[0]))
00149 -am*((directProd(p4b,p4b)-directProd(p4b,p4[0])));
00150 jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
00151 jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
00152 jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
00153
00154 break;
00155
00156 case 4:
00157 q=p4b-p4[0];
00158 q2=q*q;
00159 nbcurrent=5;
00160 ffmodel.gettensorff(B0,D3P20,EvtPDL::getMeanMass(D3P20),q2,&hf,&kf,&bp,&bm);
00161 g.setdiag(1.0,-1.0,-1.0,-1.0);
00162
00163 ep_meson_b[0] = ((p->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
00164 ep_meson_b[1] = ((p->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
00165 ep_meson_b[2] = ((p->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
00166 ep_meson_b[3] = ((p->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
00167 ep_meson_b[4] = ((p->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
00168
00169 pp=p4b+p4[0];
00170 pm=p4b-p4[0];
00171
00172 ep_meson_bb[0]=ep_meson_b[0]*(p4b);
00173 ep_meson_bb[1]=ep_meson_b[1]*(p4b);
00174 ep_meson_bb[2]=ep_meson_b[2]*(p4b);
00175 ep_meson_bb[3]=ep_meson_b[3]*(p4b);
00176 ep_meson_bb[4]=ep_meson_b[4]*(p4b);
00177
00178 jb[0]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[0])
00179 -kf*ep_meson_b[0]
00180 -bp*ep_meson_bb[0]*pp-bm*ep_meson_bb[0]*pm;
00181
00182 jb[1]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[1])
00183 -kf*ep_meson_b[1]
00184 -bp*ep_meson_bb[1]*pp-bm*ep_meson_bb[1]*pm;
00185
00186 jb[2]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[2])
00187 -kf*ep_meson_b[2]
00188 -bp*ep_meson_bb[2]*pp-bm*ep_meson_bb[2]*pm;
00189
00190 jb[3]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[3])
00191 -kf*ep_meson_b[3]
00192 -bp*ep_meson_bb[3]*pp-bm*ep_meson_bb[3]*pm;
00193
00194 jb[4]=EvtComplex(0.0,(1.0)*hf)*dual(directProd(pp,pm)).cont2(ep_meson_b[4])
00195 -kf*ep_meson_b[4]
00196 -bp*ep_meson_bb[4]*pp-bm*ep_meson_bb[4]*pm;
00197 break;
00198
00199 case 5:
00200 q=p4b-p4[0];
00201 q2=q*q;
00202 double f,gf,ap,am;
00203 nbcurrent=3;
00204 ffmodel.getvectorff(B0,D1P10,EvtPDL::getMeanMass(D1P10),q2,&f,&gf,&ap,&am);
00205 g.setdiag(1.0,-1.0,-1.0,-1.0);
00206 tds = -f*g
00207 -ap*(directProd(p4b,p4b)+directProd(p4b,p4[0]))
00208 +gf*EvtComplex(0.0,1.0)*dual(directProd(p4[0]+p4b,p4b-p4[0]))
00209 -am*((directProd(p4b,p4b)-directProd(p4b,p4[0])));
00210 jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
00211 jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
00212 jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
00213 break;
00214
00215 case 6:
00216 q=p4b-p4[0];
00217 q2=q*q;
00218 nbcurrent=1;
00219 ffmodel.getscalarff(B0,D3P00,EvtPDL::getMeanMass(D3P00),q2,&fp,&fm);
00220 jb[0]=fp*(p4b+p4[0])+fm*q;
00221 break;
00222 default:
00223 report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown hadronic current."<<endl;
00224
00225 }
00226
00227 double norm;
00228
00229 switch (wcurrent) {
00230
00231
00232 case 1: case 3: case 4:
00233 nwcurrent=1;
00234 jw[0]=p4[getNDaug()-1];
00235 break;
00236
00237 case 2: case 5: case 6:
00238 nwcurrent=3;
00239 norm=1.0/sqrt(p4[1].get(0)*p4[1].get(0)/p4[1].mass2()-1.0);
00240 jw[0]=norm*p->getDaug(getNDaug()-1)->epsParent(0);
00241 jw[1]=norm*p->getDaug(getNDaug()-1)->epsParent(1);
00242 jw[2]=norm*p->getDaug(getNDaug()-1)->epsParent(2);
00243 break;
00244
00245
00246 default:
00247 report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown W current."<<endl;
00248
00249 }
00250
00251 if (nbcurrent==1&&nwcurrent==1){
00252 vertex(jb[0]*jw[0]);
00253 return;
00254 }
00255
00256 if (nbcurrent==1){
00257
00258
00259 for(j=0;j<nwcurrent;j++){
00260
00261
00262 vertex(j,jb[0]*jw[j]);
00263 }
00264
00265 return;
00266 }
00267
00268 if (nwcurrent==1){
00269 for(j=0;j<nbcurrent;j++){
00270 vertex(j,jb[j]*jw[0]);
00271 }
00272 return;
00273 }
00274
00275 for(j=0;j<nbcurrent;j++){
00276 for(i=0;i<nwcurrent;i++){
00277 vertex(j,i,jb[j]*jw[i]);
00278 }
00279 }
00280 return;
00281
00282 }
00283
00284
00285
00286
00287
00288
00289