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/EvtPDL.hh"
00026 #include "EvtGenBase/EvtReport.hh"
00027 #include "EvtGenModels/EvtGoityRoberts.hh"
00028 #include "EvtGenBase/EvtTensor4C.hh"
00029 #include "EvtGenBase/EvtDiracSpinor.hh"
00030 #include <string>
00031 #include "EvtGenBase/EvtVector4C.hh"
00032
00033 EvtGoityRoberts::~EvtGoityRoberts() {}
00034
00035 void EvtGoityRoberts::getName(std::string& model_name){
00036
00037 model_name="GOITY_ROBERTS";
00038
00039 }
00040
00041
00042 EvtDecayBase* EvtGoityRoberts::clone(){
00043
00044 return new EvtGoityRoberts;
00045
00046 }
00047
00048 void EvtGoityRoberts::init(){
00049
00050
00051 checkNArg(0);
00052 checkNDaug(4);
00053
00054 checkSpinParent(EvtSpinType::SCALAR);
00055 checkSpinDaughter(1,EvtSpinType::SCALAR);
00056 checkSpinDaughter(2,EvtSpinType::DIRAC);
00057 checkSpinDaughter(3,EvtSpinType::NEUTRINO);
00058
00059 }
00060
00061
00062 void EvtGoityRoberts::initProbMax() {
00063
00064 setProbMax( 3000.0);
00065 }
00066
00067 void EvtGoityRoberts::decay( EvtParticle *p){
00068
00069
00070 static EvtId DST0=EvtPDL::getId("D*0");
00071 static EvtId DSTB=EvtPDL::getId("anti-D*0");
00072 static EvtId DSTP=EvtPDL::getId("D*+");
00073 static EvtId DSTM=EvtPDL::getId("D*-");
00074 static EvtId D0=EvtPDL::getId("D0");
00075 static EvtId D0B=EvtPDL::getId("anti-D0");
00076 static EvtId DP=EvtPDL::getId("D+");
00077 static EvtId DM=EvtPDL::getId("D-");
00078
00079
00080
00081 EvtId meson=getDaug(0);
00082
00083 if (meson==DST0||meson==DSTP||meson==DSTM||meson==DSTB) {
00084 DecayBDstarpilnuGR(p,getDaug(0),getDaug(1),getDaug(2),getDaug(3));
00085 }
00086 else{
00087 if (meson==D0||meson==DP||meson==DM||meson==D0B) {
00088 DecayBDpilnuGR(p,getDaug(0),getDaug(1),getDaug(2),getDaug(3));
00089 }
00090 else{
00091 report(ERROR,"EvtGen") << "Wrong daugther in EvtGoityRoberts!\n";
00092 }
00093 }
00094 return ;
00095 }
00096
00097 void EvtGoityRoberts::DecayBDstarpilnuGR(EvtParticle *pb,EvtId ndstar,
00098 EvtId npion, EvtId nlep, EvtId nnu)
00099 {
00100
00101 pb->initializePhaseSpace(getNDaug(),getDaugs());
00102
00103
00104 static EvtId EM=EvtPDL::getId("e-");
00105 static EvtId EP=EvtPDL::getId("e+");
00106 static EvtId MUM=EvtPDL::getId("mu-");
00107 static EvtId MUP=EvtPDL::getId("mu+");
00108
00109 EvtParticle *dstar, *pion, *lepton, *neutrino;
00110
00111
00112 dstar=pb->getDaug(0);
00113 pion=pb->getDaug(1);
00114 lepton=pb->getDaug(2);
00115 neutrino=pb->getDaug(3);
00116
00117 EvtVector4C l1, l2, et0, et1, et2;
00118
00119 EvtVector4R v,vp,p4_pi;
00120 double w;
00121
00122 v.set(1.0,0.0,0.0,0.0);
00123 vp=(1.0/dstar->getP4().mass())*dstar->getP4();
00124 p4_pi=pion->getP4();
00125
00126 w=v*vp;
00127
00128 EvtTensor4C omega;
00129
00130 double mb=EvtPDL::getMeanMass(pb->getId());
00131 double md=EvtPDL::getMeanMass(ndstar);
00132
00133 EvtComplex dmb(0.0460,-0.5*0.00001);
00134 EvtComplex dmd(0.1421,-0.5*0.00006);
00135
00136
00137
00138 double g = 0.5;
00139 double alpha3 = 0.690;
00140 double alpha1 = -1.430;
00141 double alpha2 = -0.140;
00142 double f0 = 0.093;
00143
00144 EvtComplex dmt3(0.563,-0.5*0.191);
00145 EvtComplex dmt1(0.392,-0.5*1.040);
00146 EvtComplex dmt2(0.709,-0.5*0.405);
00147
00148 double betas=0.285;
00149 double betap=0.280;
00150 double betad=0.260;
00151 double betasp=betas*betas+betap*betap;
00152 double betasd=betas*betas+betad*betad;
00153
00154 double lambdabar=0.750;
00155
00156
00157 double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
00158 double xi1= -1.0*sqrt(2.0/3.0)*(
00159 lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
00160 exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
00161 double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
00162 pow((2*betas*betap/(betasp)),2.5)*
00163 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
00164 double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
00165 pow((2*betas*betad/(betasd)),3.5)*
00166 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
00167
00168
00169
00170 EvtComplex h1nr,h2nr,h3nr,f1nr,f2nr;
00171 EvtComplex f3nr,f4nr,f5nr,f6nr,knr,g1nr,g2nr,g3nr,g4nr,g5nr;
00172 EvtComplex h1r,h2r,h3r,f1r,f2r,f3r,f4r,f5r,f6r,kr,g1r,g2r,g3r,g4r,g5r;
00173 EvtComplex h1,h2,h3,f1,f2,f3,f4,f5,f6,k,g1,g2,g3,g4,g5;
00174
00175
00176 h1nr = -g*xi*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmb));
00177 h2nr = -g*xi/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmb));
00178 h3nr = -(g*xi/(f0*md))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
00179 -EvtComplex((1.0+w)/(p4_pi*vp),0.0));
00180
00181 f1nr = -(g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb) -
00182 1.0/(EvtComplex(p4_pi*vp,0.0)+dmd));
00183 f2nr = f1nr*mb/md;
00184 f3nr = EvtComplex(0.0);
00185 f4nr = EvtComplex(0.0);
00186 f5nr = (g*xi/(2*f0*mb*md))*(EvtComplex(1.0,0.0)
00187 +(p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmb));
00188 f6nr = (g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
00189 -EvtComplex(1.0/(p4_pi*vp),0.0));
00190
00191 knr = (g*xi/(2*f0))*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmb) +
00192 EvtComplex((p4_pi*(v-w*vp))/(p4_pi*vp),0.0));
00193
00194 g1nr = EvtComplex(0.0);
00195 g2nr = EvtComplex(0.0);
00196 g3nr = EvtComplex(0.0);
00197 g4nr = (g*xi)/(f0*md*EvtComplex(p4_pi*vp));
00198 g5nr = EvtComplex(0.0);
00199
00200
00201 h1r = -alpha1*rho1*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
00202 alpha2*rho2*(p4_pi*(v+2.0*w*v-vp))
00203 /(3*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
00204 alpha3*xi1*(p4_pi*v)/(f0*mb*md*EvtComplex(p4_pi*v,0.0)+dmt3);
00205 h2r = -alpha2*(1+w)*rho2/(3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
00206 alpha3*xi1/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
00207 h3r = alpha2*rho2*(1+w)/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
00208 alpha3*xi1/(f0*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
00209
00210 f1r = -alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
00211 alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
00212 f2r = f1r*mb/md;
00213 f3r = EvtComplex(0.0);
00214 f4r = EvtComplex(0.0);
00215 f5r = alpha1*rho1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
00216 alpha2*rho2*(p4_pi*(vp-v/3.0-2.0/3.0*w*v))/
00217 (2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
00218 alpha3*xi1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
00219 f6r = alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
00220 alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
00221
00222 kr = -alpha1*rho1*(w-1.0)*(p4_pi*v)/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt1)) -
00223 alpha2*rho2*(w-1.0)*(p4_pi*(vp-w*v))
00224 /(3*f0*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
00225 alpha3*xi1*(p4_pi*(vp-w*v))/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt3));
00226
00227 g1r = EvtComplex(0.0);
00228 g2r = EvtComplex(0.0);
00229 g3r = -g2r;
00230 g4r = 2.0*alpha2*rho2/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2));
00231 g5r = EvtComplex(0.0);
00232
00233
00234 h1 = h1nr + h1r;
00235 h2 = h2nr + h2r;
00236 h3 = h3nr + h3r;
00237
00238 f1 = f1nr + f1r;
00239 f2 = f2nr + f2r;
00240 f3 = f3nr + f3r;
00241 f4 = f4nr + f4r;
00242 f5 = f5nr + f5r;
00243 f6 = f6nr + f6r;
00244
00245 k = knr+kr;
00246
00247 g1 = g1nr + g1r;
00248 g2 = g2nr + g2r;
00249 g3 = g3nr + g3r;
00250 g4 = g4nr + g4r;
00251 g5 = g5nr + g5r;
00252
00253 EvtTensor4C g_metric;
00254 g_metric.setdiag(1.0,-1.0,-1.0,-1.0);
00255
00256 if (nlep==EM||nlep==MUM){
00257 omega=EvtComplex(0.0,0.5)*dual(h1*mb*md*directProd(v,vp)+
00258 h2*mb*directProd(v,p4_pi)+
00259 h3*md*directProd(vp,p4_pi))+
00260 f1*mb*directProd(v,p4_pi)+f2*md*directProd(vp,p4_pi)+
00261 f3*directProd(p4_pi,p4_pi)+f4*mb*mb*directProd(v,v)+
00262 f5*mb*md*directProd(vp,v)+f6*mb*directProd(p4_pi,v)+k*g_metric+
00263 EvtComplex(0.0,0.5)*directProd(dual(directProd(vp,p4_pi)).cont2(v),
00264 (g1*p4_pi+g2*mb*v))+
00265 EvtComplex(0.0,0.5)*directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
00266 dual(directProd(vp,p4_pi)).cont2(v));
00267
00268 l1=EvtLeptonVACurrent(lepton->spParent(0),neutrino->spParentNeutrino());
00269 l2=EvtLeptonVACurrent(lepton->spParent(1),neutrino->spParentNeutrino());
00270 }
00271 else{
00272 if (nlep==EP||nlep==MUP){
00273 omega=EvtComplex(0.0,-0.5)*dual(h1*mb*md*directProd(v,vp)+
00274 h2*mb*directProd(v,p4_pi)+
00275 h3*md*directProd(vp,p4_pi))+
00276 f1*mb*directProd(v,p4_pi)+f2*md*directProd(vp,p4_pi)+
00277 f3*directProd(p4_pi,p4_pi)+f4*mb*mb*directProd(v,v)+
00278 f5*mb*md*directProd(vp,v)+f6*mb*directProd(p4_pi,v)+k*g_metric+
00279 EvtComplex(0.0,-0.5)*directProd(dual(directProd(vp,p4_pi)).cont2(v),
00280 (g1*p4_pi+g2*mb*v))+
00281 EvtComplex(0.0,-0.5)*directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
00282 dual(directProd(vp,p4_pi)).cont2(v));
00283
00284 l1=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(0));
00285 l2=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(1));
00286 }
00287 else{
00288 report(DEBUG,"EvtGen") << "42387dfs878w wrong lepton number\n";
00289 }
00290 }
00291
00292 et0=omega.cont2( dstar->epsParent(0).conj() );
00293 et1=omega.cont2( dstar->epsParent(1).conj() );
00294 et2=omega.cont2( dstar->epsParent(2).conj() );
00295
00296 vertex(0,0,l1.cont(et0));
00297 vertex(0,1,l2.cont(et0));
00298
00299 vertex(1,0,l1.cont(et1));
00300 vertex(1,1,l2.cont(et1));
00301
00302 vertex(2,0,l1.cont(et2));
00303 vertex(2,1,l2.cont(et2));
00304
00305 return;
00306
00307 }
00308
00309 void EvtGoityRoberts::DecayBDpilnuGR(EvtParticle *pb,EvtId nd,
00310 EvtId npion, EvtId nlep, EvtId nnu)
00311
00312 {
00313
00314 static EvtId EM=EvtPDL::getId("e-");
00315 static EvtId EP=EvtPDL::getId("e+");
00316 static EvtId MUM=EvtPDL::getId("mu-");
00317 static EvtId MUP=EvtPDL::getId("mu+");
00318
00319 EvtParticle *d, *pion, *lepton, *neutrino;
00320
00321 pb->initializePhaseSpace(getNDaug(),getDaugs());
00322 d=pb->getDaug(0);
00323 pion=pb->getDaug(1);
00324 lepton=pb->getDaug(2);
00325 neutrino=pb->getDaug(3);
00326
00327 EvtVector4C l1, l2, et0, et1, et2;
00328
00329 EvtVector4R v,vp,p4_pi;
00330 double w;
00331
00332 v.set(1.0,0.0,0.0,0.0);
00333 vp=(1.0/d->getP4().mass())*d->getP4();
00334 p4_pi=pion->getP4();
00335 w=v*vp;
00336
00337 double mb=EvtPDL::getMeanMass(pb->getId());
00338 double md=EvtPDL::getMeanMass(nd);
00339 EvtComplex dmb(0.0460,-0.5*0.00001);
00340
00341
00342
00343
00344 double g = 0.5;
00345 double alpha3 = 0.690;
00346 double alpha1 = -1.430;
00347 double alpha2 = -0.140;
00348 double f0=0.093;
00349
00350 EvtComplex dmt3(0.563,-0.5*0.191);
00351 EvtComplex dmt1(0.392,-0.5*1.040);
00352 EvtComplex dmt2(0.709,-0.5*0.405);
00353
00354 double betas=0.285;
00355 double betap=0.280;
00356 double betad=0.260;
00357 double betasp=betas*betas+betap*betap;
00358 double betasd=betas*betas+betad*betad;
00359
00360 double lambdabar=0.750;
00361
00362
00363 double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
00364 double xi1= -1.0*sqrt(2.0/3.0)*(lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
00365 exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
00366 double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
00367 pow((2*betas*betap/(betasp)),2.5)*
00368 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
00369 double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
00370 pow((2*betas*betad/(betasd)),3.5)*
00371 exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
00372
00373 EvtComplex h,a1,a2,a3;
00374 EvtComplex hnr,a1nr,a2nr,a3nr;
00375 EvtComplex hr,a1r,a2r,a3r;
00376
00377
00378 hnr = g*xi*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb*md);
00379 a1nr= -1.0*g*xi*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0);
00380 a2nr= g*xi*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb);
00381 a3nr=EvtComplex(0.0,0.0);
00382
00383
00384 hr = alpha2*rho2*(w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0*mb*md) +
00385 alpha3*xi1*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb*md);
00386 a1r= -1.0*alpha2*rho2*(w*w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0) -
00387 alpha3*xi1*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0);
00388 a2r= alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*mb) +
00389 alpha2*rho2*(0.5*p4_pi*(w*vp-v)+p4_pi*(vp-w*v))/
00390 (3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
00391 alpha3*xi1*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb);
00392 a3r= -1.0*alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*md) -
00393 alpha2*rho2*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmt2))/(2*f0*md);
00394
00395
00396 h=hnr+hr;
00397 a1=a1nr+a1r;
00398 a2=a2nr+a2r;
00399 a3=a3nr+a3r;
00400
00401 EvtVector4C omega;
00402
00403 if ( nlep==EM|| nlep==MUM ) {
00404 omega=EvtComplex(0.0,-1.0)*h*mb*md*dual(directProd(vp,p4_pi)).cont2(v)+
00405 a1*p4_pi+a2*mb*v+a3*md*vp;
00406 l1=EvtLeptonVACurrent(
00407 lepton->spParent(0),neutrino->spParentNeutrino());
00408 l2=EvtLeptonVACurrent(
00409 lepton->spParent(1),neutrino->spParentNeutrino());
00410 }
00411 else{
00412 if ( nlep==EP|| nlep==MUP ) {
00413 omega=EvtComplex(0.0,1.0)*h*mb*md*dual(directProd(vp,p4_pi)).cont2(v)+
00414 a1*p4_pi+a2*mb*v+a3*md*vp;
00415 l1=EvtLeptonVACurrent(
00416 neutrino->spParentNeutrino(),lepton->spParent(0));
00417 l2=EvtLeptonVACurrent(
00418 neutrino->spParentNeutrino(),lepton->spParent(1));
00419 }
00420 else{
00421 report(ERROR,"EvtGen") << "42387dfs878w wrong lepton number\n";
00422 }
00423 }
00424
00425 vertex(0,l1*omega);
00426 vertex(1,l2*omega);
00427
00428 return;
00429
00430 }
00431