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/EvtRandom.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenModels/EvtVectorIsr.hh"
00029 #include "EvtGenBase/EvtReport.hh"
00030 #include "EvtGenBase/EvtConst.hh"
00031 #include <string>
00032 #include "EvtGenBase/EvtVector4C.hh"
00033
00034 EvtVectorIsr::~EvtVectorIsr() {}
00035
00036 void EvtVectorIsr::getName(std::string& model_name){
00037
00038 model_name="VECTORISR";
00039
00040 }
00041
00042 EvtDecayBase* EvtVectorIsr::clone(){
00043
00044 return new EvtVectorIsr;
00045
00046 }
00047
00048 void EvtVectorIsr::init(){
00049
00050
00051 checkNArg(2);
00052 checkNDaug(2);
00053
00054 checkSpinParent(EvtSpinType::VECTOR);
00055 checkSpinDaughter(0,EvtSpinType::VECTOR);
00056 checkSpinDaughter(1,EvtSpinType::PHOTON);
00057
00058
00059
00060 csfrmn=getArg(0);
00061 csbkmn=getArg(1);
00062
00063
00064 }
00065
00066 void EvtVectorIsr::initProbMax(){
00067
00068 noProbMax();
00069
00070 }
00071
00072 void EvtVectorIsr::decay( EvtParticle *p ){
00073
00074 EvtParticle *phi;
00075 EvtParticle *gamma;
00076
00077
00078 double electMass=EvtPDL::getMeanMass(EvtPDL::getId("e-"));
00079
00080
00081
00082
00083
00084 p->initializePhaseSpace(getNDaug(),getDaugs());
00085 phi=p->getDaug(0);
00086 gamma=p->getDaug(1);
00087
00088 double wcm=p->mass();
00089 double beta=2.*electMass/wcm;
00090 beta=sqrt(1. - beta*beta);
00091
00092
00093 double pg=(wcm*wcm-phi->mass()*phi->mass())/(2*wcm);
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 double ymax=log((1.+beta*csfrmn)/(1.-beta*csfrmn));
00106 double ymin=log((1.-beta*csbkmn)/(1.+beta*csbkmn));
00107
00108
00109 double y=(ymax-ymin)*EvtRandom::Flat(0.0,1.0) + ymin;
00110 double cs=exp(y);
00111 cs=(cs - 1.)/(cs + 1.)/beta;
00112 double sn=sqrt(1-cs*cs);
00113
00114 double fi=EvtRandom::Flat(EvtConst::twoPi);
00115
00116
00117 EvtVector4R p4phi(sqrt(phi->mass()*phi->mass()+pg*pg),pg*sn*cos(fi),
00118 pg*sn*sin(fi),pg*cs);
00119
00120 EvtVector4R p4gamma(pg,-p4phi.get(1),-p4phi.get(2),-p4phi.get(3));
00121
00122
00123 phi->init( getDaug(0),p4phi);
00124 gamma->init( getDaug(1),p4gamma);
00125
00126
00127
00128
00129
00130
00131 EvtVector4C phi0=phi->epsParent(0);
00132 EvtVector4C phi1=phi->epsParent(1);
00133 EvtVector4C phi2=phi->epsParent(2);
00134
00135 EvtVector4C gamma0=gamma->epsParentPhoton(0);
00136 EvtVector4C gamma1=gamma->epsParentPhoton(1);
00137
00138 EvtComplex r1p=phi0*gamma0;
00139 EvtComplex r2p=phi1*gamma0;
00140 EvtComplex r3p=phi2*gamma0;
00141
00142
00143 EvtComplex r1m=phi0*gamma1;
00144 EvtComplex r2m=phi1*gamma1;
00145 EvtComplex r3m=phi2*gamma1;
00146
00147 EvtComplex rho33=r3p*conj(r3p)+r3m*conj(r3m);
00148 EvtComplex rho22=r2p*conj(r2p)+r2m*conj(r2m);
00149 EvtComplex rho11=r1p*conj(r1p)+r1m*conj(r1m);
00150
00151 EvtComplex rho13=r3p*conj(r1p)+r3m*conj(r1m);
00152 EvtComplex rho12=r2p*conj(r1p)+r2m*conj(r1m);
00153 EvtComplex rho23=r3p*conj(r2p)+r3m*conj(r2m);
00154
00155 EvtComplex rho31=conj(rho13);
00156 EvtComplex rho32=conj(rho23);
00157 EvtComplex rho21=conj(rho12);
00158
00159
00160 EvtSpinDensity rho;
00161 rho.SetDim(3);
00162
00163 rho.Set(0,0,rho11);
00164 rho.Set(0,1,rho12);
00165 rho.Set(0,2,rho13);
00166 rho.Set(1,0,rho21);
00167 rho.Set(1,1,rho22);
00168 rho.Set(1,2,rho23);
00169 rho.Set(2,0,rho31);
00170 rho.Set(2,1,rho32);
00171 rho.Set(2,2,rho33);
00172
00173
00174 setDaughterSpinDensity(0);
00175 phi->setSpinDensityForward(rho);
00176
00177
00178 return ;
00179 }
00180