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 <iostream>
00023 #include <math.h>
00024 #include <assert.h>
00025 #include "EvtGenBase/EvtComplex.hh"
00026 #include "EvtGenBase/EvtGammaMatrix.hh"
00027 #include "EvtGenBase/EvtDiracSpinor.hh"
00028 #include "EvtGenBase/EvtReport.hh"
00029 #include "EvtGenBase/EvtTensor4C.hh"
00030 #include "EvtGenBase/EvtVector4C.hh"
00031 using std::endl;
00032 using std::ostream;
00033
00034 EvtGammaMatrix::EvtGammaMatrix(){
00035 int i,j;
00036
00037 for(i=0;i<4;i++){
00038 for(j=0;j<4;j++){
00039 gamma[i][j]=EvtComplex(0.0,0.0);
00040 }
00041 }
00042 }
00043
00044 EvtGammaMatrix operator*(const EvtGammaMatrix& g, const EvtComplex& c)
00045 {
00046 return c*g;
00047 }
00048
00049
00050 EvtGammaMatrix operator*(const EvtComplex& c,const EvtGammaMatrix& g){
00051 int i,j;
00052
00053 EvtGammaMatrix temp;
00054
00055 for(i=0;i<4;i++){
00056 for(j=0;j<4;j++){
00057 temp.gamma[i][j]=g.gamma[i][j]*c;
00058 }
00059 }
00060
00061 return temp;
00062
00063 }
00064
00065
00066 ostream& operator<<(ostream& s, const EvtGammaMatrix& g){
00067
00068
00069 s<<"["<<g.gamma[0][0]<<","<<g.gamma[0][1]<<","<<g.gamma[0][2]<<","<<g.gamma[0][3]<<"]"<<endl;
00070 s<<"["<<g.gamma[1][0]<<","<<g.gamma[1][1]<<","<<g.gamma[1][2]<<","<<g.gamma[1][3]<<"]"<<endl;
00071 s<<"["<<g.gamma[2][0]<<","<<g.gamma[2][1]<<","<<g.gamma[2][2]<<","<<g.gamma[2][3]<<"]"<<endl;
00072 s<<"["<<g.gamma[3][0]<<","<<g.gamma[3][1]<<","<<g.gamma[3][2]<<","<<g.gamma[3][3]<<"]"<<endl;
00073
00074 return s;
00075
00076 }
00077
00078
00079
00080 EvtGammaMatrix::EvtGammaMatrix(const EvtGammaMatrix& gm){
00081 int i,j;
00082
00083 for(i=0;i<4;i++){
00084 for(j=0;j<4;j++){
00085 gamma[i][j]=gm.gamma[i][j];
00086 }
00087 }
00088 }
00089
00090 EvtGammaMatrix::~EvtGammaMatrix() {}
00091
00092 EvtGammaMatrix& EvtGammaMatrix::operator=(const EvtGammaMatrix& gm){
00093 int i,j;
00094
00095 for(i=0;i<4;i++){
00096 for(j=0;j<4;j++){
00097 gamma[i][j]=gm.gamma[i][j];
00098 }
00099 }
00100 return *this;
00101 }
00102
00103 void EvtGammaMatrix::init(){
00104 int i,j;
00105
00106 for(i=0;i<4;i++){
00107 for(j=0;j<4;j++){
00108 gamma[i][j]=EvtComplex(0.0,0.0);
00109 }
00110 }
00111 }
00112
00113 const EvtGammaMatrix& EvtGammaMatrix::va0(){
00114
00115 static EvtGammaMatrix g;
00116 static int first=1;
00117
00118 if (first){
00119 g.gamma[0][0]=EvtComplex(1.0,0.0);
00120 g.gamma[0][1]=EvtComplex(0.0,0.0);
00121 g.gamma[0][2]=EvtComplex(-1.0,0.0);
00122 g.gamma[0][3]=EvtComplex(0.0,0.0);
00123 g.gamma[1][0]=EvtComplex(0.0,0.0);
00124 g.gamma[1][1]=EvtComplex(1.0,0.0);
00125 g.gamma[1][2]=EvtComplex(0.0,0.0);
00126 g.gamma[1][3]=EvtComplex(-1.0,0.0);
00127 g.gamma[2][0]=EvtComplex(-1.0,0.0);
00128 g.gamma[2][1]=EvtComplex(0.0,0.0);
00129 g.gamma[2][2]=EvtComplex(1.0,0.0);
00130 g.gamma[2][3]=EvtComplex(0.0,0.0);
00131 g.gamma[3][0]=EvtComplex(0.0,0.0);
00132 g.gamma[3][1]=EvtComplex(-1.0,0.0);
00133 g.gamma[3][2]=EvtComplex(0.0,0.0);
00134 g.gamma[3][3]=EvtComplex(1.0,0.0);
00135 }
00136
00137 return g;
00138
00139 }
00140
00141
00142 const EvtGammaMatrix& EvtGammaMatrix::va1(){
00143
00144 static EvtGammaMatrix g;
00145 static int first=1;
00146
00147 if (first){
00148 g.gamma[0][0]=EvtComplex(0.0,0.0);
00149 g.gamma[0][1]=EvtComplex(-1.0,0.0);
00150 g.gamma[0][2]=EvtComplex(0.0,0.0);
00151 g.gamma[0][3]=EvtComplex(1.0,0.0);
00152 g.gamma[1][0]=EvtComplex(-1.0,0.0);
00153 g.gamma[1][1]=EvtComplex(0.0,0.0);
00154 g.gamma[1][2]=EvtComplex(1.0,0.0);
00155 g.gamma[1][3]=EvtComplex(0.0,0.0);
00156 g.gamma[2][0]=EvtComplex(0.0,0.0);
00157 g.gamma[2][1]=EvtComplex(1.0,0.0);
00158 g.gamma[2][2]=EvtComplex(0.0,0.0);
00159 g.gamma[2][3]=EvtComplex(-1.0,0.0);
00160 g.gamma[3][0]=EvtComplex(1.0,0.0);
00161 g.gamma[3][1]=EvtComplex(0.0,0.0);
00162 g.gamma[3][2]=EvtComplex(-1.0,0.0);
00163 g.gamma[3][3]=EvtComplex(0.0,0.0);
00164 }
00165
00166 return g;
00167
00168 }
00169
00170
00171
00172 const EvtGammaMatrix& EvtGammaMatrix::va2(){
00173
00174 static EvtGammaMatrix g;
00175 static int first=1;
00176
00177 if (first){
00178 g.gamma[0][0]=EvtComplex(0.0,0.0);
00179 g.gamma[0][1]=EvtComplex(0.0,1.0);
00180 g.gamma[0][2]=EvtComplex(0.0,0.0);
00181 g.gamma[0][3]=EvtComplex(0.0,-1.0);
00182 g.gamma[1][0]=EvtComplex(0.0,-1.0);
00183 g.gamma[1][1]=EvtComplex(0.0,0.0);
00184 g.gamma[1][2]=EvtComplex(0.0,1.0);
00185 g.gamma[1][3]=EvtComplex(0.0,0.0);
00186 g.gamma[2][0]=EvtComplex(0.0,0.0);
00187 g.gamma[2][1]=EvtComplex(0.0,-1.0);
00188 g.gamma[2][2]=EvtComplex(0.0,0.0);
00189 g.gamma[2][3]=EvtComplex(0.0,1.0);
00190 g.gamma[3][0]=EvtComplex(0.0,1.0);
00191 g.gamma[3][1]=EvtComplex(0.0,0.0);
00192 g.gamma[3][2]=EvtComplex(0.0,-1.0);
00193 g.gamma[3][3]=EvtComplex(0.0,0.0);
00194 }
00195
00196 return g;
00197
00198 }
00199
00200
00201
00202
00203 const EvtGammaMatrix& EvtGammaMatrix::va3(){
00204
00205 static EvtGammaMatrix g;
00206 static int first=1;
00207
00208 if (first){
00209 g.gamma[0][0]=EvtComplex(-1.0,0.0);
00210 g.gamma[0][1]=EvtComplex(0.0,0.0);
00211 g.gamma[0][2]=EvtComplex(1.0,0.0);
00212 g.gamma[0][3]=EvtComplex(0.0,0.0);
00213 g.gamma[1][0]=EvtComplex(0.0,0.0);
00214 g.gamma[1][1]=EvtComplex(1.0,0.0);
00215 g.gamma[1][2]=EvtComplex(0.0,0.0);
00216 g.gamma[1][3]=EvtComplex(-1.0,0.0);
00217 g.gamma[2][0]=EvtComplex(1.0,0.0);
00218 g.gamma[2][1]=EvtComplex(0.0,0.0);
00219 g.gamma[2][2]=EvtComplex(-1.0,0.0);
00220 g.gamma[2][3]=EvtComplex(0.0,0.0);
00221 g.gamma[3][0]=EvtComplex(0.0,0.0);
00222 g.gamma[3][1]=EvtComplex(-1.0,0.0);
00223 g.gamma[3][2]=EvtComplex(0.0,0.0);
00224 g.gamma[3][3]=EvtComplex(1.0,0.0);
00225 }
00226
00227 return g;
00228
00229 }
00230
00231
00232
00233
00234
00235 const EvtGammaMatrix& EvtGammaMatrix::g0(){
00236
00237 static EvtGammaMatrix g;
00238 static int first=1;
00239
00240 if (first){
00241
00242 first=0;
00243
00244 int i,j;
00245
00246 for(i=0;i<4;i++){
00247 for(j=0;j<4;j++){
00248 g.gamma[i][j]=EvtComplex(0.0,0.0);
00249 }
00250 }
00251
00252 g.gamma[0][0]=EvtComplex(1.0,0.0);
00253 g.gamma[1][1]=EvtComplex(1.0,0.0);
00254 g.gamma[2][2]=EvtComplex(-1.0,0.0);
00255 g.gamma[3][3]=EvtComplex(-1.0,0.0);
00256 }
00257
00258 return g;
00259
00260 }
00261
00262
00263
00264
00265 const EvtGammaMatrix& EvtGammaMatrix::g1(){
00266
00267 static EvtGammaMatrix g;
00268 static int first=1;
00269
00270 if (first){
00271 first=0;
00272 int i,j;
00273
00274 for(i=0;i<4;i++){
00275 for(j=0;j<4;j++){
00276 g.gamma[i][j]=EvtComplex(0.0,0.0);
00277 }
00278 }
00279
00280 g.gamma[0][3]=EvtComplex(1.0,0.0);
00281 g.gamma[1][2]=EvtComplex(1.0,0.0);
00282 g.gamma[2][1]=EvtComplex(-1.0,0.0);
00283 g.gamma[3][0]=EvtComplex(-1.0,0.0);
00284 }
00285
00286 return g;
00287
00288 }
00289
00290
00291
00292
00293 const EvtGammaMatrix& EvtGammaMatrix::g2(){
00294
00295 static EvtGammaMatrix g;
00296 static int first=1;
00297
00298 if (first){
00299 first=0;
00300 int i,j;
00301
00302 for(i=0;i<4;i++){
00303 for(j=0;j<4;j++){
00304 g.gamma[i][j]=EvtComplex(0.0,0.0);
00305 }
00306 }
00307
00308 g.gamma[0][3]=EvtComplex(0.0,-1.0);
00309 g.gamma[1][2]=EvtComplex(0.0,1.0);
00310 g.gamma[2][1]=EvtComplex(0.0,1.0);
00311 g.gamma[3][0]=EvtComplex(0.0,-1.0);
00312 }
00313
00314 return g;
00315
00316 }
00317
00318
00319
00320
00321
00322 const EvtGammaMatrix& EvtGammaMatrix::g3(){
00323
00324 static EvtGammaMatrix g;
00325 static int first=1;
00326
00327 if (first){
00328 first=0;
00329 int i,j;
00330
00331 for(i=0;i<4;i++){
00332 for(j=0;j<4;j++){
00333 g.gamma[i][j]=EvtComplex(0.0,0.0);
00334 }
00335 }
00336
00337 g.gamma[0][2]=EvtComplex(1.0,0.0);
00338 g.gamma[1][3]=EvtComplex(-1.0,0.0);
00339 g.gamma[2][0]=EvtComplex(-1.0,0.0);
00340 g.gamma[3][1]=EvtComplex(1.0,0.0);
00341 }
00342
00343 return g;
00344
00345 }
00346
00347
00348
00349
00350 const EvtGammaMatrix& EvtGammaMatrix::g5(){
00351
00352 static EvtGammaMatrix g;
00353 static int first=1;
00354
00355 if (first){
00356
00357 int i,j;
00358
00359 for(i=0;i<4;i++){
00360 for(j=0;j<4;j++){
00361 g.gamma[i][j]=EvtComplex(0.0,0.0);
00362 }
00363 }
00364
00365 g.gamma[0][2]=EvtComplex(1.0,0.0);
00366 g.gamma[1][3]=EvtComplex(1.0,0.0);
00367 g.gamma[2][0]=EvtComplex(1.0,0.0);
00368 g.gamma[3][1]=EvtComplex(1.0,0.0);
00369 }
00370
00371 return g;
00372
00373 }
00374
00375
00376
00377 const EvtGammaMatrix& EvtGammaMatrix::v0(){
00378
00379 static EvtGammaMatrix g;
00380 static int first=1;
00381
00382 if (first){
00383
00384 int i,j;
00385
00386 for(i=0;i<4;i++){
00387 for(j=0;j<4;j++){
00388 g.gamma[i][j]=EvtComplex(0.0,0.0);
00389 }
00390 }
00391
00392 g.gamma[0][0]=EvtComplex(1.0,0.0);
00393 g.gamma[1][1]=EvtComplex(1.0,0.0);
00394 g.gamma[2][2]=EvtComplex(1.0,0.0);
00395 g.gamma[3][3]=EvtComplex(1.0,0.0);
00396 }
00397
00398 return g;
00399
00400 }
00401
00402
00403
00404
00405
00406 const EvtGammaMatrix& EvtGammaMatrix::v1(){
00407
00408 static EvtGammaMatrix g;
00409 static int first=1;
00410
00411 if (first){
00412
00413 int i,j;
00414
00415 for(i=0;i<4;i++){
00416 for(j=0;j<4;j++){
00417 g.gamma[i][j]=EvtComplex(0.0,0.0);
00418 }
00419 }
00420
00421 g.gamma[0][3]=EvtComplex(1.0,0.0);
00422 g.gamma[1][2]=EvtComplex(1.0,0.0);
00423 g.gamma[2][1]=EvtComplex(1.0,0.0);
00424 g.gamma[3][0]=EvtComplex(1.0,0.0);
00425 }
00426
00427 return g;
00428
00429 }
00430
00431
00432
00433
00434 const EvtGammaMatrix& EvtGammaMatrix::v2(){
00435
00436 static EvtGammaMatrix g;
00437 static int first=1;
00438
00439 if (first){
00440
00441 int i,j;
00442
00443 for(i=0;i<4;i++){
00444 for(j=0;j<4;j++){
00445 g.gamma[i][j]=EvtComplex(0.0,0.0);
00446 }
00447 }
00448
00449 g.gamma[0][3]=EvtComplex(0.0,-1.0);
00450 g.gamma[1][2]=EvtComplex(0.0,1.0);
00451 g.gamma[2][1]=EvtComplex(0.0,-1.0);
00452 g.gamma[3][0]=EvtComplex(0.0,1.0);
00453 }
00454
00455 return g;
00456
00457 }
00458
00459
00460
00461
00462 const EvtGammaMatrix& EvtGammaMatrix::v3(){
00463
00464 static EvtGammaMatrix g;
00465 static int first=1;
00466
00467 if (first){
00468
00469 int i,j;
00470
00471 for(i=0;i<4;i++){
00472 for(j=0;j<4;j++){
00473 g.gamma[i][j]=EvtComplex(0.0,0.0);
00474 }
00475 }
00476
00477 g.gamma[0][2]=EvtComplex(1.0,0.0);
00478 g.gamma[1][3]=EvtComplex(-1.0,0.0);
00479 g.gamma[2][0]=EvtComplex(1.0,0.0);
00480 g.gamma[3][1]=EvtComplex(-1.0,0.0);
00481 }
00482
00483 return g;
00484
00485 }
00486
00487
00488
00489
00490
00491 const EvtGammaMatrix& EvtGammaMatrix::id(){
00492
00493 static EvtGammaMatrix g;
00494 static int first=1;
00495
00496 if (first){
00497
00498 int i,j;
00499
00500 for(i=0;i<4;i++){
00501 for(j=0;j<4;j++){
00502 g.gamma[i][j]=EvtComplex(0.0,0.0);
00503 }
00504 }
00505
00506 g.gamma[0][0]=EvtComplex(1.0,0.0);
00507 g.gamma[1][1]=EvtComplex(1.0,0.0);
00508 g.gamma[2][2]=EvtComplex(1.0,0.0);
00509 g.gamma[3][3]=EvtComplex(1.0,0.0);
00510 }
00511
00512 return g;
00513
00514 }
00515
00516
00517
00518
00519 EvtGammaMatrix& EvtGammaMatrix::operator+=(const EvtGammaMatrix &g){
00520
00521 int i,j;
00522
00523 for(i=0;i<4;i++){
00524 for(j=0;j<4;j++){
00525 gamma[i][j]+=g.gamma[i][j];
00526 }
00527 }
00528 return *this;
00529 }
00530
00531
00532
00533
00534
00535 EvtGammaMatrix& EvtGammaMatrix::operator-=(const EvtGammaMatrix &g){
00536
00537 int i,j;
00538
00539 for(i=0;i<4;i++){
00540 for(j=0;j<4;j++){
00541 gamma[i][j]-=g.gamma[i][j];
00542 }
00543 }
00544 return *this;
00545 }
00546
00547
00548
00549 EvtGammaMatrix& EvtGammaMatrix::operator*=(const EvtGammaMatrix &g){
00550
00551 int i,j,k;
00552 EvtGammaMatrix temp;
00553
00554 for(i=0;i<4;i++){
00555 for(j=0;j<4;j++){
00556 temp.gamma[i][j]=EvtComplex(0.0,0.0);
00557 for(k=0;k<4;k++){
00558 temp.gamma[i][j]+=gamma[i][k]*g.gamma[k][j];
00559 }
00560 }
00561 }
00562
00563 for(i=0;i<4;i++){
00564 for(j=0;j<4;j++){
00565 gamma[i][j]=temp.gamma[i][j];
00566 }
00567 }
00568
00569 return *this;
00570 }
00571
00572
00573 EvtDiracSpinor operator*(const EvtGammaMatrix& g,const EvtDiracSpinor& d){
00574
00575 int i,j;
00576 EvtDiracSpinor temp;
00577
00578 for(i=0;i<4;i++){
00579 temp.set_spinor(i,EvtComplex(0.0,0.0));
00580 for(j=0;j<4;j++){
00581 temp.set_spinor(i,temp.get_spinor(i)+g.gamma[i][j]*d.get_spinor(j));
00582 }
00583 }
00584
00585 return temp;
00586 }
00587
00588
00589 EvtComplex operator*(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
00590
00591 int i;
00592 EvtComplex temp;
00593
00594 temp=EvtComplex(0.0,0.0);
00595
00596 for(i=0;i<4;i++){
00597 temp+=::conj(d.get_spinor(i))*dp.get_spinor(i);
00598 }
00599 return temp;
00600 }
00601
00602
00603
00604 const EvtGammaMatrix& EvtGammaMatrix::sigmaUpper(unsigned int mu, unsigned int nu)
00605 {
00606 EvtGammaMatrix a, b;
00607 static const EvtTensor4C eta = EvtTensor4C::g();
00608 static EvtGammaMatrix sigma[4][4];
00609 static bool hasBeenCalled = false;
00610 if (!hasBeenCalled)
00611 {
00612 EvtComplex I(0, 1);
00613 for (int i=0; i<4; ++i)
00614 sigma[i][i].init();
00615
00616 EvtGammaMatrix s01 = I/2 * (g0()*g1() - g1()*g0());
00617 EvtGammaMatrix s02 = I/2 * (g0()*g2() - g2()*g0());
00618 EvtGammaMatrix s03 = I/2 * (g0()*g3() - g3()*g0());
00619 EvtGammaMatrix s12 = I/2 * (g1()*g2() - g2()*g1());
00620 EvtGammaMatrix s13 = I/2 * (g1()*g3() - g3()*g1());
00621 EvtGammaMatrix s23 = I/2 * (g2()*g3() - g3()*g2());
00622 sigma[0][1] = s01;
00623 sigma[1][0] = -1*s01;
00624 sigma[0][2] = s02;
00625 sigma[2][0] = -1*s02;
00626 sigma[0][3] = s03;
00627 sigma[3][0] = -1*s03;
00628 sigma[1][2] = s12;
00629 sigma[2][1] = -1*s12;
00630 sigma[1][3] = s13;
00631 sigma[3][1] = -1*s13;
00632 sigma[2][3] = s23;
00633 sigma[3][2] = -1*s23;
00634 }
00635 hasBeenCalled = true;
00636
00637 if (mu > 3 || nu > 3)
00638 {
00639 report(ERROR, "EvtSigmaTensor") << "Expected index between 0 and 3, but found " << nu << "!" << endl;
00640 assert(0);
00641 }
00642 return sigma[mu][nu];
00643
00644 }
00645
00646 const EvtGammaMatrix& EvtGammaMatrix::sigmaLower(unsigned int mu, unsigned int nu)
00647 {
00648 const EvtComplex I(0, 1);
00649 EvtGammaMatrix a, b;
00650 static EvtGammaMatrix sigma[4][4];
00651 static bool hasBeenCalled = false;
00652 static const EvtTensor4C eta = EvtTensor4C::g();
00653
00654 if (!hasBeenCalled)
00655 {
00656
00657 for (int i=0; i<4; ++i)
00658 {
00659 a = eta.get(i, 0)*g0() + eta.get(i, 1)*g1() + eta.get(i, 2)*g2() + eta.get(i, 3)*g3();
00660 for (int j=0; j<4; ++j)
00661 {
00662 b = eta.get(j, 0)*g0() + eta.get(j, 1)*g1() + eta.get(j, 2)*g2() + eta.get(j, 3)*g3();
00663 sigma[i][j] = I/2 * (a*b - b*a);
00664 }
00665 }
00666 }
00667 return sigma[mu][nu];
00668 }
00669
00670
00671 EvtGammaMatrix slash(const EvtVector4C& p)
00672 {
00673 return EvtGammaMatrix::g0()*p.get(0) + EvtGammaMatrix::g1()*p.get(1) + EvtGammaMatrix::g2()*p.get(2) + EvtGammaMatrix::g3()*p.get(3);
00674 }