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 <math.h>
00023 #include <assert.h>
00024 #include "EvtGenBase/EvtDiracSpinor.hh"
00025 #include "EvtGenBase/EvtGammaMatrix.hh"
00026 #include "EvtGenBase/EvtComplex.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtVector4C.hh"
00029 #include "EvtGenBase/EvtTensor4C.hh"
00030 using std::ostream;
00031
00032
00033 EvtDiracSpinor::~EvtDiracSpinor(){}
00034
00035 EvtDiracSpinor::EvtDiracSpinor(const EvtComplex& sp0,const EvtComplex& sp1,
00036 const EvtComplex& sp2,const EvtComplex& sp3){
00037 set(sp0,sp1,sp2,sp3);
00038 }
00039
00040 void EvtDiracSpinor::set(const EvtComplex& sp0,const EvtComplex& sp1,
00041 const EvtComplex& sp2,const EvtComplex& sp3){
00042
00043 spinor[0]=sp0;spinor[1]=sp1;spinor[2]=sp2;spinor[3]=sp3;
00044 }
00045
00046 void EvtDiracSpinor::set_spinor(int i,const EvtComplex& sp){
00047
00048 spinor[i]=sp;
00049 }
00050
00051 ostream& operator<<(ostream& s, const EvtDiracSpinor& sp){
00052
00053 s <<"["<<sp.spinor[0]<<","<<sp.spinor[1]<<","
00054 <<sp.spinor[2]<<","<<sp.spinor[3]<<"]";
00055 return s;
00056
00057 }
00058
00059
00060 const EvtComplex& EvtDiracSpinor::get_spinor(int i) const {
00061
00062 return spinor[i];
00063
00064 }
00065
00066 EvtDiracSpinor rotateEuler(const EvtDiracSpinor& sp,
00067 double alpha,double beta,double gamma){
00068
00069 EvtDiracSpinor tmp(sp);
00070 tmp.applyRotateEuler(alpha,beta,gamma);
00071 return tmp;
00072
00073 }
00074
00075 EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
00076 const EvtVector4R p4){
00077
00078 EvtDiracSpinor tmp(sp);
00079 tmp.applyBoostTo(p4);
00080 return tmp;
00081
00082 }
00083
00084 EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
00085 const EvtVector3R boost){
00086
00087 EvtDiracSpinor tmp(sp);
00088 tmp.applyBoostTo(boost);
00089 return tmp;
00090
00091 }
00092
00093 void EvtDiracSpinor::applyBoostTo(const EvtVector4R& p4){
00094
00095 double e=p4.get(0);
00096
00097 EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
00098
00099 applyBoostTo(boost);
00100
00101 return;
00102
00103 }
00104
00105
00106
00107 void EvtDiracSpinor::applyBoostTo(const EvtVector3R& boost) {
00108
00109 double bx,by,bz,gamma,b2,f1,f2;
00110 EvtComplex spinorp[4];
00111
00112 bx=boost.get(0);
00113 by=boost.get(1);
00114 bz=boost.get(2);
00115 b2=bx*bx+by*by+bz*bz;
00116
00117 if (b2==0.0){
00118 return;
00119 }
00120
00121 assert(b2<1.0);
00122
00123 gamma=1.0/sqrt(1-b2);
00124
00125 f1=sqrt((gamma+1.0)/2.0);
00126 f2=f1*gamma/(gamma+1.0);
00127
00128 spinorp[0]=f1*spinor[0]+f2*bz*spinor[2]+
00129 f2*EvtComplex(bx,-by)*spinor[3];
00130 spinorp[1]=f1*spinor[1]+f2*EvtComplex(bx,by)*spinor[2]-
00131 f2*bz*spinor[3];
00132 spinorp[2]=f2*bz*spinor[0]+f2*EvtComplex(bx,-by)*spinor[1]+
00133 f1*spinor[2];
00134 spinorp[3]=f2*EvtComplex(bx,by)*spinor[0]-
00135 f2*bz*spinor[1]+f1*spinor[3];
00136
00137 spinor[0]=spinorp[0];
00138 spinor[1]=spinorp[1];
00139 spinor[2]=spinorp[2];
00140 spinor[3]=spinorp[3];
00141
00142 return;
00143 }
00144
00145 void EvtDiracSpinor::applyRotateEuler(double alpha,double beta,
00146 double gamma) {
00147
00148 EvtComplex retVal[4];
00149
00150 double cb2=cos(0.5*beta);
00151 double sb2=sin(0.5*beta);
00152 double capg2=cos(0.5*(alpha+gamma));
00153 double camg2=cos(0.5*(alpha-gamma));
00154 double sapg2=sin(0.5*(alpha+gamma));
00155 double samg2=sin(0.5*(alpha-gamma));
00156
00157 EvtComplex m11(cb2*capg2,-cb2*sapg2);
00158 EvtComplex m12(-sb2*camg2,sb2*samg2);
00159 EvtComplex m21(sb2*camg2,sb2*samg2);
00160 EvtComplex m22(cb2*capg2,cb2*sapg2);
00161
00162 retVal[0]=m11*spinor[0]+m12*spinor[1];
00163 retVal[1]=m21*spinor[0]+m22*spinor[1];
00164 retVal[2]=m11*spinor[2]+m12*spinor[3];
00165 retVal[3]=m21*spinor[2]+m22*spinor[3];
00166
00167 spinor[0]=retVal[0];
00168 spinor[1]=retVal[1];
00169 spinor[2]=retVal[2];
00170 spinor[3]=retVal[3];
00171
00172 return;
00173 }
00174
00175
00176
00177 EvtDiracSpinor EvtDiracSpinor::conj() const {
00178
00179 EvtDiracSpinor sp;
00180
00181 for ( int i=0; i<4; i++)
00182 sp.set_spinor(i,::conj(spinor[i]));
00183
00184 return sp;
00185 }
00186
00187 EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 EvtComplex u02=::conj(d.spinor[0]-d.spinor[2]);
00204 EvtComplex u13=::conj(d.spinor[1]-d.spinor[3]);
00205
00206 EvtComplex v02=dp.spinor[0]-dp.spinor[2];
00207 EvtComplex v13=dp.spinor[1]-dp.spinor[3];
00208
00209 EvtComplex a=u02*v02;
00210 EvtComplex b=u13*v13;
00211
00212 EvtComplex c=u02*v13;
00213 EvtComplex e=u13*v02;
00214
00215 return EvtVector4C(a+b,-(c+e),EvtComplex(0,1)*(c-e),b-a);
00216
00217
00218 }
00219
00220 EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00221
00222 EvtVector4C temp;
00223
00224
00225
00226
00227 temp.set(0,d*(EvtGammaMatrix::v0()*dp));
00228 temp.set(1,d*(EvtGammaMatrix::v1()*dp));
00229 temp.set(2,d*(EvtGammaMatrix::v2()*dp));
00230 temp.set(3,d*(EvtGammaMatrix::v3()*dp));
00231
00232 return temp;
00233 }
00234
00235
00236 EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00237
00238 EvtVector4C temp;
00239
00240 EvtGammaMatrix mat;
00241
00242
00243
00244
00245 mat = EvtGammaMatrix::v0()-EvtGammaMatrix::va0();
00246 temp.set(0,d*(mat*dp));
00247
00248 mat = EvtGammaMatrix::v1()-EvtGammaMatrix::va1();
00249 temp.set(1,d*(mat*dp));
00250
00251 mat = EvtGammaMatrix::v2()-EvtGammaMatrix::va2();
00252 temp.set(2,d*(mat*dp));
00253
00254 mat = EvtGammaMatrix::v3()-EvtGammaMatrix::va3();
00255 temp.set(3,d*(mat*dp));
00256
00257 return temp;
00258 }
00259
00260 EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00261
00262 EvtComplex temp;
00263
00264
00265
00266
00267 temp=d*(EvtGammaMatrix::g0()*dp);
00268
00269 return temp;
00270 }
00271
00272 EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00273
00274 EvtComplex temp;
00275
00276
00277
00278 static EvtGammaMatrix m=EvtGammaMatrix::g0()*EvtGammaMatrix::g5();
00279 temp=d*(m*dp);
00280
00281 return temp;
00282 }
00283
00284 EvtTensor4C EvtLeptonTCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00285
00286 EvtTensor4C temp;
00287 temp.zero();
00288 EvtComplex i2(0,0.5);
00289
00290 static EvtGammaMatrix mat01=EvtGammaMatrix::g0()*
00291 (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()-
00292 EvtGammaMatrix::g1()*EvtGammaMatrix::g0());
00293 static EvtGammaMatrix mat02=EvtGammaMatrix::g0()*
00294 (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()-
00295 EvtGammaMatrix::g2()*EvtGammaMatrix::g0());
00296 static EvtGammaMatrix mat03=EvtGammaMatrix::g0()*
00297 (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()-
00298 EvtGammaMatrix::g3()*EvtGammaMatrix::g0());
00299 static EvtGammaMatrix mat12=EvtGammaMatrix::g0()*
00300 (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()-
00301 EvtGammaMatrix::g2()*EvtGammaMatrix::g1());
00302 static EvtGammaMatrix mat13=EvtGammaMatrix::g0()*
00303 (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()-
00304 EvtGammaMatrix::g3()*EvtGammaMatrix::g1());
00305 static EvtGammaMatrix mat23=EvtGammaMatrix::g0()*
00306 (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()-
00307 EvtGammaMatrix::g3()*EvtGammaMatrix::g2());
00308
00309
00310 temp.set(0,1,i2*(d*(mat01*dp)));
00311 temp.set(1,0,-temp.get(0,1));
00312
00313 temp.set(0,2,i2*(d*(mat02*dp)));
00314 temp.set(2,0,-temp.get(0,2));
00315
00316 temp.set(0,3,i2*(d*(mat03*dp)));
00317 temp.set(3,0,-temp.get(0,3));
00318
00319 temp.set(1,2,i2*(d*(mat12*dp)));
00320 temp.set(2,1,-temp.get(1,2));
00321
00322 temp.set(1,3,i2*(d*(mat13*dp)));
00323 temp.set(3,1,-temp.get(1,3));
00324
00325 temp.set(2,3,i2*(d*(mat23*dp)));
00326 temp.set(3,2,-temp.get(2,3));
00327
00328 return temp;
00329 }
00330
00331 EvtTensor4C EvtLeptonTg5Current(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00332
00333 EvtTensor4C temp;
00334 temp.zero();
00335 EvtComplex i2(0,0.5);
00336
00337 static EvtGammaMatrix mat01=EvtGammaMatrix::g0()*
00338 (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()-
00339 EvtGammaMatrix::g1()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5();
00340 static EvtGammaMatrix mat02=EvtGammaMatrix::g0()*
00341 (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()-
00342 EvtGammaMatrix::g2()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5();
00343 static EvtGammaMatrix mat03=EvtGammaMatrix::g0()*
00344 (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()-
00345 EvtGammaMatrix::g3()*EvtGammaMatrix::g0())*EvtGammaMatrix::g5();
00346 static EvtGammaMatrix mat12=EvtGammaMatrix::g0()*
00347 (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()-
00348 EvtGammaMatrix::g2()*EvtGammaMatrix::g1())*EvtGammaMatrix::g5();
00349 static EvtGammaMatrix mat13=EvtGammaMatrix::g0()*
00350 (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()-
00351 EvtGammaMatrix::g3()*EvtGammaMatrix::g1())*EvtGammaMatrix::g5();
00352 static EvtGammaMatrix mat23=EvtGammaMatrix::g0()*
00353 (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()-
00354 EvtGammaMatrix::g3()*EvtGammaMatrix::g2())*EvtGammaMatrix::g5();
00355
00356
00357 temp.set(0,1,i2*(d*(mat01*dp)));
00358 temp.set(1,0,-temp.get(0,1));
00359
00360 temp.set(0,2,i2*(d*(mat02*dp)));
00361 temp.set(2,0,-temp.get(0,2));
00362
00363 temp.set(0,3,i2*(d*(mat03*dp)));
00364 temp.set(3,0,-temp.get(0,3));
00365
00366 temp.set(1,2,i2*(d*(mat12*dp)));
00367 temp.set(2,1,-temp.get(1,2));
00368
00369 temp.set(1,3,i2*(d*(mat13*dp)));
00370 temp.set(3,1,-temp.get(1,3));
00371
00372 temp.set(2,3,i2*(d*(mat23*dp)));
00373 temp.set(3,2,-temp.get(2,3));
00374
00375 return temp;
00376 }
00377
00378 EvtDiracSpinor operator*(const EvtComplex& c, const EvtDiracSpinor& d) {
00379 EvtDiracSpinor result;
00380 result.spinor[0] = c*d.spinor[0];
00381 result.spinor[1] = c*d.spinor[1];
00382 result.spinor[2] = c*d.spinor[2];
00383 result.spinor[3] = c*d.spinor[3];
00384 return result;
00385 }
00386
00387 EvtDiracSpinor EvtDiracSpinor::adjoint() const
00388 {
00389 EvtDiracSpinor d = this->conj();
00390 EvtGammaMatrix g0 = EvtGammaMatrix::g0();
00391 EvtDiracSpinor result;
00392
00393 for (int i=0; i<4; ++i)
00394 for (int j=0; j<4; ++j)
00395 result.spinor[i] += d.spinor[j] * g0.gamma[i][j];
00396
00397 return result;
00398 }