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 <iostream>
00024 #include <math.h>
00025 #include "EvtGenBase/EvtComplex.hh"
00026 #include "EvtGenBase/EvtPhotonParticle.hh"
00027 #include "EvtGenBase/EvtVector4C.hh"
00028 #include "EvtGenBase/EvtReport.hh"
00029 using std::endl;
00030
00031 EvtPhotonParticle::~EvtPhotonParticle(){}
00032
00033 void EvtPhotonParticle::init(EvtId part_n,const EvtVector4R& p4){
00034
00035 init(part_n,p4.get(0),p4.get(1),
00036 p4.get(2),p4.get(3));
00037 }
00038
00039 void EvtPhotonParticle::init(EvtId part_n,double e,double px,double py,double pz){
00040
00041 _validP4=true;
00042 setp(e,px,py,pz);
00043 setpart_num(part_n);
00044
00045 setLifetime();
00046
00047
00048 _evalBasis=0;
00049
00050 }
00051
00052
00053 EvtVector4C EvtPhotonParticle::epsParentPhoton(int i){
00054
00055 if (!_evalBasis){
00056
00057 _evalBasis=1;
00058 eps1.set(EvtComplex(0.0,0.0), EvtComplex(-1.0/sqrt(2.0),0.0),
00059 EvtComplex(0.0,-1.0/sqrt(2.0)), EvtComplex(0.0,0.0));
00060 eps2.set(EvtComplex(0.0,0.0), EvtComplex(1.0/sqrt(2.0),0.0),
00061 EvtComplex(0.0,-1.0/sqrt(2.0)), EvtComplex(0.0,0.0));
00062
00063
00064
00065
00066 double phi,theta;
00067
00068 EvtVector4R p=this->getP4();
00069
00070 double px=p.get(1);
00071 double py=p.get(2);
00072 double pz=p.get(3);
00073
00074 phi = atan2(py,px);
00075 theta = acos(pz/sqrt(px*px+py*py+pz*pz));
00076 eps1.applyRotateEuler(phi,theta,-phi);
00077 eps2.applyRotateEuler(phi,theta,-phi);
00078
00079 }
00080
00081
00082 EvtVector4C temp;
00083
00084 switch(i) {
00085
00086 case 0:
00087 temp=eps1;
00088 break;
00089 case 1:
00090 temp=eps2;
00091 break;
00092 default:
00093 report(ERROR,"EvtGen") << "EvtPhotonParticle.cc: Asked "
00094 << "for state:"<<i<<endl;
00095 ::abort();
00096 break;
00097 }
00098
00099 return temp;
00100 }
00101
00102 EvtVector4C EvtPhotonParticle::epsPhoton(int i){
00103
00104 report(ERROR,"EvtGen") << "EvtPhotonParticle.cc: Can not get "
00105 << "state in photons restframe."<<endl;;
00106 ::abort();
00107 return EvtVector4C();
00108
00109 }
00110
00111
00112 EvtSpinDensity EvtPhotonParticle::rotateToHelicityBasis() const {
00113
00114 EvtVector4C eplus(0.0,-1.0/sqrt(2.0),
00115 EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00116 EvtVector4C eminus(0.0,1.0/sqrt(2.0),
00117 EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00118
00119
00120
00121 EvtVector4C e1=((EvtParticle*)this)->epsParentPhoton(0);
00122 EvtVector4C e2=((EvtParticle*)this)->epsParentPhoton(1);
00123
00124
00125 EvtSpinDensity R;
00126 R.SetDim(2);
00127
00128 R.Set(0,0,(eplus.conj())*e1);
00129 R.Set(0,1,(eplus.conj())*e2);
00130
00131 R.Set(1,0,(eminus.conj())*e1);
00132 R.Set(1,1,(eminus.conj())*e2);
00133
00134 return R;
00135
00136 }
00137
00138
00139 EvtSpinDensity EvtPhotonParticle::rotateToHelicityBasis(double alpha,
00140 double beta,
00141 double gamma) const{
00142
00143 EvtVector4C eplus(0.0,-1.0/sqrt(2.0),
00144 EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00145 EvtVector4C eminus(0.0,1.0/sqrt(2.0),
00146 EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
00147
00148 eplus.applyRotateEuler(alpha,beta,gamma);
00149 eminus.applyRotateEuler(alpha,beta,gamma);
00150
00151
00152
00153
00154 EvtVector4C e1=((EvtParticle*)this)->epsParentPhoton(0);
00155 EvtVector4C e2=((EvtParticle*)this)->epsParentPhoton(1);
00156
00157 EvtSpinDensity R;
00158 R.SetDim(2);
00159
00160 R.Set(0,0,(eplus.conj())*e1);
00161 R.Set(0,1,(eplus.conj())*e2);
00162
00163 R.Set(1,0,(eminus.conj())*e1);
00164 R.Set(1,1,(eminus.conj())*e2);
00165
00166 return R;
00167
00168 }
00169
00170