00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "EvtGenBase/EvtPatches.hh"
00024 #include <stdlib.h>
00025 #include "EvtGenBase/EvtParticle.hh"
00026 #include "EvtGenBase/EvtGenKine.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenBase/EvtVector4C.hh"
00029 #include "EvtGenBase/EvtVector4R.hh"
00030 #include "EvtGenBase/EvtTensor4C.hh"
00031 #include "EvtGenBase/EvtReport.hh"
00032 #include "EvtGenBase/EvtHelSys.hh"
00033 #include "EvtJ2BB3.hh"
00034 #include "EvtGenBase/EvtdFunctionSingle.hh"
00035 #include "EvtGenBase/EvtRaritaSchwinger.hh"
00036 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
00037
00038 #include <string>
00039
00040 using std::endl;
00041 using std::cout;
00042 EvtJ2BB3::~EvtJ2BB3() {}
00043
00044 void EvtJ2BB3::getName(std::string& model_name){
00045
00046 model_name="J2BB3";
00047
00048 }
00049
00050
00051 EvtDecayBase* EvtJ2BB3::clone(){
00052
00053 return new EvtJ2BB3;
00054
00055 }
00056
00057 void EvtJ2BB3::init(){
00058
00059 checkNDaug(2);
00060 checkSpinParent(EvtSpinType::VECTOR);
00061 checkSpinDaughter(0,EvtSpinType::DIRAC);
00062 checkSpinDaughter(1,EvtSpinType::RARITASCHWINGER);
00063 }
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 void EvtJ2BB3::decay( EvtParticle *p){
00075 loop:
00076 p->initializePhaseSpace(getNDaug(),getDaugs());
00077
00078 EvtParticle *v,*s1;
00079 EvtVector4R pv,ps,ppr;
00080
00081 v =p->getDaug(0);
00082 s1=p->getDaug(1);
00083 pv=v->getP4();
00084 ps=s1->getP4();
00085 ppr=p->getP4();
00086
00087 EvtHelSys angles(ppr,pv);
00088 double theta =angles.getHelAng(1);
00089 double phi =angles.getHelAng(2);
00090 double gamma=0;
00091
00092 double m_b0 = EvtPDL::getMass(getDaug(0));
00093 double m_b1 = EvtPDL::getMass(getDaug(1));
00094
00095
00096 double m_M = EvtPDL::getMass(getParentId());
00097 double M2=pow(m_M,2.0);double b2_p=pow((m_b0+m_b1),2.0);
00098 double b2_m=pow((m_b1-m_b0),2.0);
00099
00100 if(M2-b2_p<=0) goto loop;
00101
00102 double P_c=sqrt((M2-b2_p)*(M2-b2_m)/(4.0*M2));
00103 double Q=M2-b2_p;
00104 double F1=-0.54/(m_b0*b2_p);
00105 double F2=-0.54*(m_b1-m_b0)/(m_b0*(m_b0+m_b1));
00106 double F3=-0.89/m_b0;
00107 double H1,H0,Hm1;
00108
00109
00110
00111 if(getNArg()>0)
00112 { alpha=getArg(0);
00113 H0=1.0;H1=sqrt((1+alpha)/(1-alpha));Hm1=H1;
00114 }
00115 else{
00116 H1=-1.15*P_c*m_M*(-M2*Q*F1-Q*F2+(M2-m_b0*(m_b0+m_b1))*F3)/(sqrt(Q)*m_b1);
00117 H0=-0.82*P_c*M2*(-(m_b1-m_b0)*Q*F1/m_b1-Q*F2/(m_b1*(m_b1-m_b0))+2.0*F3)/sqrt(Q);
00118 Hm1=2*P_c*m_M*(m_b0+m_b1)*F3/sqrt(Q);
00119 alpha=(H1*H1+Hm1*Hm1-2*H0*H0)/(H1*H1+Hm1*Hm1+2*H0*H0);
00120 }
00121
00122
00123
00124 vertex(0,0,0,Djmn(1, 1, 0,phi,theta,gamma)*H0);
00125 vertex(0,0,1,Djmn(1, 1,-1,phi,theta,gamma)*Hm1);
00126 vertex(0,0,2,0.0);
00127 vertex(0,0,3,Djmn(1, 1, 1,phi,theta,gamma)*H1);
00128 vertex(0,1,0,Djmn(1, 1,-1,phi,theta,gamma)*H1);
00129 vertex(0,1,1,0.0);
00130 vertex(0,1,2,Djmn(1, 1, 1,phi,theta,gamma)*Hm1);
00131 vertex(0,1,3,Djmn(1, 1, 0,phi,theta,gamma)*H0);
00132
00133 vertex(1,0,0,Djmn(1,-1, 0,phi,theta,gamma)*H0);
00134 vertex(1,0,1,Djmn(1,-1,-1,phi,theta,gamma)*Hm1);
00135 vertex(1,0,2,0.0);
00136 vertex(1,0,3,Djmn(1,-1, 1,phi,theta,gamma)*H1);
00137 vertex(1,1,0,Djmn(1,-1,-1,phi,theta,gamma)*H1);
00138 vertex(1,1,1,0.0);
00139 vertex(1,1,2,Djmn(1,-1, 1,phi,theta,gamma)*Hm1);
00140 vertex(1,1,3,Djmn(1,-1, 0,phi,theta,gamma)*H0);
00141
00142
00143 vertex(2,0,0,Djmn(1, 0, 0,phi,theta,gamma)*H0);
00144 vertex(2,0,1,Djmn(1, 0,-1,phi,theta,gamma)*Hm1);
00145 vertex(2,0,2,0.0);
00146 vertex(2,0,3,Djmn(1, 0, 1,phi,theta,gamma)*H1);
00147 vertex(2,1,0,Djmn(1, 0,-1,phi,theta,gamma)*H1);
00148 vertex(2,1,1,0.0);
00149 vertex(2,1,2,Djmn(1, 0, 1,phi,theta,gamma)*Hm1);
00150 vertex(2,1,3,Djmn(1, 0, 0,phi,theta,gamma)*H0);
00151
00152 return ;
00153
00154 }
00155
00156
00157
00158