00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <stdlib.h>
00022 #include "EvtGenBase/EvtParticle.hh"
00023 #include "EvtGenBase/EvtGenKine.hh"
00024 #include "EvtGenBase/EvtPDL.hh"
00025 #include "EvtGenBase/EvtVector4C.hh"
00026 #include "EvtGenBase/EvtVector4R.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtRandom.hh"
00029 #include "EvtGenModels/EvtVPHOtoVISR.hh"
00030 #include <string>
00031
00032 EvtVPHOtoVISR::~EvtVPHOtoVISR() {}
00033
00034 void EvtVPHOtoVISR::getName(std::string& model_name){
00035
00036 model_name="VPHOTOVISR";
00037
00038 }
00039
00040
00041 EvtDecayBase* EvtVPHOtoVISR::clone(){
00042
00043 return new EvtVPHOtoVISR;
00044
00045 }
00046
00047 void EvtVPHOtoVISR::init(){
00048
00049
00050 checkNArg(0,2);
00051
00052
00053 checkNDaug(2);
00054
00055
00056 checkSpinParent(EvtSpinType::VECTOR);
00057 checkSpinDaughter(0,EvtSpinType::VECTOR);
00058 checkSpinDaughter(1,EvtSpinType::PHOTON);
00059 }
00060
00061 void EvtVPHOtoVISR::initProbMax() {
00062
00063
00064
00065 }
00066
00067 void EvtVPHOtoVISR::decay( EvtParticle *p){
00068
00069
00070
00071
00072
00073 double w=p->mass();
00074 double s=w*w;
00075
00076 double L=2.0*log(w/0.000511);
00077 double alpha=1/137.0;
00078 double beta=(L-1)*2.0*alpha/EvtConst::pi;
00079
00080
00081
00082 assert(p->getDaug(0)->getDaug(0)!=0);
00083 double md=EvtPDL::getMeanMass(p->getDaug(0)->getDaug(0)->getId());
00084
00085 static double mD0=EvtPDL::getMeanMass(EvtPDL::getId("D0"));
00086 static double mDp=EvtPDL::getMeanMass(EvtPDL::getId("D+"));
00087
00088 double pgmax=(s-4.0*md*md)/(2.0*w);
00089
00090 assert(pgmax>0.0);
00091
00092 double pgz=0.99*pgmax*exp(log(EvtRandom::Flat(1.0))/beta);
00093
00094 if (EvtRandom::Flat(1.0)<0.5) pgz=-pgz;
00095
00096 double k=fabs(pgz);
00097
00098 EvtVector4R p4g(k,0.0,0.0,pgz);
00099
00100 EvtVector4R p4res=p->getP4Restframe()-p4g;
00101
00102 double mres=p4res.mass();
00103
00104 double ed=mres/2.0;
00105
00106 assert(ed>md);
00107
00108 double pd=sqrt(ed*ed-md*md);
00109
00110
00111
00112
00113 p->getDaug(0)->init(getDaug(0),p4res);
00114 p->getDaug(1)->init(getDaug(1),p4g);
00115
00116
00117 double sigma=beta*pow(2/w,beta)*(1+alpha*(1.5*L-2.0+EvtConst::pi*EvtConst::pi/3.0)/EvtConst::pi);
00118
00119 double m=EvtPDL::getMeanMass(p->getDaug(0)->getId());
00120 double Gamma=EvtPDL::getWidth(p->getDaug(0)->getId());
00121
00122
00123
00124 double p0=0.0;
00125 if (ed>mD0) p0=sqrt(ed*ed-mD0*mD0);
00126 double pp=0.0;
00127 if (ed>mDp) pp=sqrt(ed*ed-mDp*mDp);
00128
00129 double p0norm=sqrt(0.25*m*m-mD0*mD0);
00130 double ppnorm=sqrt(0.25*m*m-mDp*mDp);
00131
00132 double r0=12.7;
00133 double rp=12.7;
00134
00135 if (getNArg()==2){
00136 r0=getArg(0);
00137 rp=getArg(1);
00138 }
00139
00140 double GammaTot=Gamma*(pp*pp*pp/(1+pp*pp*rp*rp)+p0*p0*p0/(1+p0*p0*r0*r0))/
00141 (ppnorm*ppnorm*ppnorm/(1+ppnorm*ppnorm*rp*rp)+
00142 p0norm*p0norm*p0norm/(1+p0norm*p0norm*r0*r0));
00143
00144
00145 sigma*=pd*pd*pd/((mres-m)*(mres-m)+0.25*GammaTot*GammaTot);
00146
00147 assert(sigma>0.0);
00148
00149 static double sigmax=sigma;
00150
00151 if (sigma>sigmax){
00152 sigmax=sigma;
00153 }
00154
00155
00156
00157 static int count=0;
00158
00159 count++;
00160
00161
00162
00163
00164
00165
00166 double norm=sqrt(sigma);
00167
00168 vertex(0,0,0,norm*p->eps(0)*p->epsParent(0).conj());
00169 vertex(1,0,0,norm*p->eps(1)*p->epsParent(0).conj());
00170 vertex(2,0,0,norm*p->eps(2)*p->epsParent(0).conj());
00171
00172 vertex(0,1,0,norm*p->eps(0)*p->epsParent(1).conj());
00173 vertex(1,1,0,norm*p->eps(1)*p->epsParent(1).conj());
00174 vertex(2,1,0,norm*p->eps(2)*p->epsParent(1).conj());
00175
00176 vertex(0,2,0,norm*p->eps(0)*p->epsParent(2).conj());
00177 vertex(1,2,0,norm*p->eps(1)*p->epsParent(2).conj());
00178 vertex(2,2,0,norm*p->eps(2)*p->epsParent(2).conj());
00179
00180 vertex(0,0,1,norm*p->eps(0)*p->epsParent(0).conj());
00181 vertex(1,0,1,norm*p->eps(1)*p->epsParent(0).conj());
00182 vertex(2,0,1,norm*p->eps(2)*p->epsParent(0).conj());
00183
00184 vertex(0,1,1,norm*p->eps(0)*p->epsParent(1).conj());
00185 vertex(1,1,1,norm*p->eps(1)*p->epsParent(1).conj());
00186 vertex(2,1,1,norm*p->eps(2)*p->epsParent(1).conj());
00187
00188 vertex(0,2,1,norm*p->eps(0)*p->epsParent(2).conj());
00189 vertex(1,2,1,norm*p->eps(1)*p->epsParent(2).conj());
00190 vertex(2,2,1,norm*p->eps(2)*p->epsParent(2).conj());
00191
00192 return;
00193 }