TEvtGen/EvtGenBase/EvtGammaMatrix.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtGammaMatrix.cc
  12 //
  13 // Description: Make gamma matrices availible for the calc. of amplitudes, etc.
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     September 25, 1996         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <math.h>
  24 #include <assert.h>
  25 #include "EvtGenBase/EvtComplex.hh"
  26 #include "EvtGenBase/EvtGammaMatrix.hh"
  27 #include "EvtGenBase/EvtDiracSpinor.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include "EvtGenBase/EvtTensor4C.hh"
  30 #include "EvtGenBase/EvtVector4C.hh"
  31 #include <stdlib.h>
  32 using std::endl;
  33 using std::ostream;
  34
  35 EvtGammaMatrix::EvtGammaMatrix(){
  36   int i,j;
  37
  38   static EvtComplex zero(0.0,0.0);
  39
  40   for(i=0;i<4;i++){
  41     for(j=0;j<4;j++){
  42       _gamma[i][j]=zero;
  43     }
  44   }
  45 }
  46
  47 EvtGammaMatrix operator*(const EvtGammaMatrix& g, const EvtComplex& c)
  48 {
  49     return c*g;
  50 }
  51
  52
  53 EvtGammaMatrix operator*(const EvtComplex& c,const EvtGammaMatrix& g){
  54   int i,j;
  55
  56   EvtGammaMatrix temp;
  57
  58   for(i=0;i<4;i++){
  59     for(j=0;j<4;j++){
  60       temp._gamma[i][j]=g._gamma[i][j]*c;
  61     }
  62   }
  63
  64   return temp;
  65
  66 }
  67
  68
  69 ostream& operator<<(ostream& s, const EvtGammaMatrix& g){
  70
  71
  72   s<<"["<<g._gamma[0][0]<<","<<g._gamma[0][1]<<","<<g._gamma[0][2]<<","<<g._gamma[0][3]<<"]"<<endl;
  73   s<<"["<<g._gamma[1][0]<<","<<g._gamma[1][1]<<","<<g._gamma[1][2]<<","<<g._gamma[1][3]<<"]"<<endl;
  74   s<<"["<<g._gamma[2][0]<<","<<g._gamma[2][1]<<","<<g._gamma[2][2]<<","<<g._gamma[2][3]<<"]"<<endl;
  75   s<<"["<<g._gamma[3][0]<<","<<g._gamma[3][1]<<","<<g._gamma[3][2]<<","<<g._gamma[3][3]<<"]"<<endl;
  76
  77   return s;
  78
  79 }
  80
  81
  82
  83 EvtGammaMatrix::EvtGammaMatrix(const EvtGammaMatrix& gm){
  84   int i,j;
  85
  86   for(i=0;i<4;i++){
  87     for(j=0;j<4;j++){
  88       _gamma[i][j]=gm._gamma[i][j];
  89     }
  90   }
  91 }
  92
  93 EvtGammaMatrix::~EvtGammaMatrix() {}
  94
  95 EvtGammaMatrix& EvtGammaMatrix::operator=(const EvtGammaMatrix& gm){
  96   int i,j;
  97
  98   for(i=0;i<4;i++){
  99     for(j=0;j<4;j++){
 100       _gamma[i][j]=gm._gamma[i][j];
 101     }
 102   }
 103   return *this;
 104 }
 105
 106 void EvtGammaMatrix::init(){
 107   int i,j;
 108
 109   static EvtComplex zero(0.0,0.0);
 110
 111   for(i=0;i<4;i++){
 112     for(j=0;j<4;j++){
 113       _gamma[i][j]=zero;
 114     }
 115   }
 116 }
 117
 118 const EvtGammaMatrix& EvtGammaMatrix::va0(){
 119
 120   static EvtGammaMatrix g;
 121   static int first=1;
 122
 123   if (first){
 124     first = 0;
 125     g._gamma[0][0]=EvtComplex(1.0,0.0);
 126     g._gamma[0][1]=EvtComplex(0.0,0.0);
 127     g._gamma[0][2]=EvtComplex(-1.0,0.0);
 128     g._gamma[0][3]=EvtComplex(0.0,0.0);
 129     g._gamma[1][0]=EvtComplex(0.0,0.0);
 130     g._gamma[1][1]=EvtComplex(1.0,0.0);
 131     g._gamma[1][2]=EvtComplex(0.0,0.0);
 132     g._gamma[1][3]=EvtComplex(-1.0,0.0);
 133     g._gamma[2][0]=EvtComplex(-1.0,0.0);
 134     g._gamma[2][1]=EvtComplex(0.0,0.0);
 135     g._gamma[2][2]=EvtComplex(1.0,0.0);
 136     g._gamma[2][3]=EvtComplex(0.0,0.0);
 137     g._gamma[3][0]=EvtComplex(0.0,0.0);
 138     g._gamma[3][1]=EvtComplex(-1.0,0.0);
 139     g._gamma[3][2]=EvtComplex(0.0,0.0);
 140     g._gamma[3][3]=EvtComplex(1.0,0.0);
 141   }
 142
 143   return g;
 144
 145 }
 146
 147
 148 const EvtGammaMatrix& EvtGammaMatrix::va1(){
 149
 150   static EvtGammaMatrix g;
 151   static int first=1;
 152
 153   if (first){
 154     first = 0;
 155     g._gamma[0][0]=EvtComplex(0.0,0.0);
 156     g._gamma[0][1]=EvtComplex(-1.0,0.0);
 157     g._gamma[0][2]=EvtComplex(0.0,0.0);
 158     g._gamma[0][3]=EvtComplex(1.0,0.0);
 159     g._gamma[1][0]=EvtComplex(-1.0,0.0);
 160     g._gamma[1][1]=EvtComplex(0.0,0.0);
 161     g._gamma[1][2]=EvtComplex(1.0,0.0);
 162     g._gamma[1][3]=EvtComplex(0.0,0.0);
 163     g._gamma[2][0]=EvtComplex(0.0,0.0);
 164     g._gamma[2][1]=EvtComplex(1.0,0.0);
 165     g._gamma[2][2]=EvtComplex(0.0,0.0);
 166     g._gamma[2][3]=EvtComplex(-1.0,0.0);
 167     g._gamma[3][0]=EvtComplex(1.0,0.0);
 168     g._gamma[3][1]=EvtComplex(0.0,0.0);
 169     g._gamma[3][2]=EvtComplex(-1.0,0.0);
 170     g._gamma[3][3]=EvtComplex(0.0,0.0);
 171   }
 172
 173   return g;
 174
 175 }
 176
 177
 178
 179 const EvtGammaMatrix& EvtGammaMatrix::va2(){
 180
 181   static EvtGammaMatrix g;
 182   static int first=1;
 183
 184   if (first){
 185     first = 0;
 186     g._gamma[0][0]=EvtComplex(0.0,0.0);
 187     g._gamma[0][1]=EvtComplex(0.0,1.0);
 188     g._gamma[0][2]=EvtComplex(0.0,0.0);
 189     g._gamma[0][3]=EvtComplex(0.0,-1.0);
 190     g._gamma[1][0]=EvtComplex(0.0,-1.0);
 191     g._gamma[1][1]=EvtComplex(0.0,0.0);
 192     g._gamma[1][2]=EvtComplex(0.0,1.0);
 193     g._gamma[1][3]=EvtComplex(0.0,0.0);
 194     g._gamma[2][0]=EvtComplex(0.0,0.0);
 195     g._gamma[2][1]=EvtComplex(0.0,-1.0);
 196     g._gamma[2][2]=EvtComplex(0.0,0.0);
 197     g._gamma[2][3]=EvtComplex(0.0,1.0);
 198     g._gamma[3][0]=EvtComplex(0.0,1.0);
 199     g._gamma[3][1]=EvtComplex(0.0,0.0);
 200     g._gamma[3][2]=EvtComplex(0.0,-1.0);
 201     g._gamma[3][3]=EvtComplex(0.0,0.0);
 202   }
 203
 204   return g;
 205
 206 }
 207
 208
 209
 210
 211 const EvtGammaMatrix& EvtGammaMatrix::va3(){
 212
 213   static EvtGammaMatrix g;
 214   static int first=1;
 215
 216   if (first){
 217     first = 0;
 218     g._gamma[0][0]=EvtComplex(-1.0,0.0);
 219     g._gamma[0][1]=EvtComplex(0.0,0.0);
 220     g._gamma[0][2]=EvtComplex(1.0,0.0);
 221     g._gamma[0][3]=EvtComplex(0.0,0.0);
 222     g._gamma[1][0]=EvtComplex(0.0,0.0);
 223     g._gamma[1][1]=EvtComplex(1.0,0.0);
 224     g._gamma[1][2]=EvtComplex(0.0,0.0);
 225     g._gamma[1][3]=EvtComplex(-1.0,0.0);
 226     g._gamma[2][0]=EvtComplex(1.0,0.0);
 227     g._gamma[2][1]=EvtComplex(0.0,0.0);
 228     g._gamma[2][2]=EvtComplex(-1.0,0.0);
 229     g._gamma[2][3]=EvtComplex(0.0,0.0);
 230     g._gamma[3][0]=EvtComplex(0.0,0.0);
 231     g._gamma[3][1]=EvtComplex(-1.0,0.0);
 232     g._gamma[3][2]=EvtComplex(0.0,0.0);
 233     g._gamma[3][3]=EvtComplex(1.0,0.0);
 234   }
 235
 236   return g;
 237
 238 }
 239
 240
 241
 242
 243
 244 const EvtGammaMatrix& EvtGammaMatrix::g0(){
 245
 246   static EvtGammaMatrix g;
 247   static int first=1;
 248
 249   if (first){
 250
 251     first=0;
 252
 253     int i,j;
 254
 255     for(i=0;i<4;i++){
 256       for(j=0;j<4;j++){
 257         g._gamma[i][j]=EvtComplex(0.0,0.0);
 258       }
 259     }
 260
 261     g._gamma[0][0]=EvtComplex(1.0,0.0);
 262     g._gamma[1][1]=EvtComplex(1.0,0.0);
 263     g._gamma[2][2]=EvtComplex(-1.0,0.0);
 264     g._gamma[3][3]=EvtComplex(-1.0,0.0);
 265   }
 266
 267   return g;
 268
 269 }
 270
 271
 272
 273
 274 const EvtGammaMatrix& EvtGammaMatrix::g1(){
 275
 276   static EvtGammaMatrix g;
 277   static int first=1;
 278
 279   if (first){
 280     first=0;
 281     int i,j;
 282
 283     for(i=0;i<4;i++){
 284       for(j=0;j<4;j++){
 285         g._gamma[i][j]=EvtComplex(0.0,0.0);
 286       }
 287     }
 288
 289     g._gamma[0][3]=EvtComplex(1.0,0.0);
 290     g._gamma[1][2]=EvtComplex(1.0,0.0);
 291     g._gamma[2][1]=EvtComplex(-1.0,0.0);
 292     g._gamma[3][0]=EvtComplex(-1.0,0.0);
 293   }
 294
 295   return g;
 296
 297 }
 298
 299
 300
 301
 302 const EvtGammaMatrix& EvtGammaMatrix::g2(){
 303
 304   static EvtGammaMatrix g;
 305   static int first=1;
 306
 307   if (first){
 308     first=0;
 309     int i,j;
 310
 311     for(i=0;i<4;i++){
 312       for(j=0;j<4;j++){
 313         g._gamma[i][j]=EvtComplex(0.0,0.0);
 314       }
 315     }
 316
 317     g._gamma[0][3]=EvtComplex(0.0,-1.0);
 318     g._gamma[1][2]=EvtComplex(0.0,1.0);
 319     g._gamma[2][1]=EvtComplex(0.0,1.0);
 320     g._gamma[3][0]=EvtComplex(0.0,-1.0);
 321   }
 322
 323   return g;
 324
 325 }
 326
 327
 328
 329
 330
 331 const EvtGammaMatrix& EvtGammaMatrix::g3(){
 332
 333   static EvtGammaMatrix g;
 334   static int first=1;
 335
 336   if (first){
 337     first=0;
 338     int i,j;
 339
 340     for(i=0;i<4;i++){
 341       for(j=0;j<4;j++){
 342         g._gamma[i][j]=EvtComplex(0.0,0.0);
 343       }
 344     }
 345
 346     g._gamma[0][2]=EvtComplex(1.0,0.0);
 347     g._gamma[1][3]=EvtComplex(-1.0,0.0);
 348     g._gamma[2][0]=EvtComplex(-1.0,0.0);
 349     g._gamma[3][1]=EvtComplex(1.0,0.0);
 350   }
 351
 352   return g;
 353
 354 }
 355
 356
 357
 358
 359 const EvtGammaMatrix& EvtGammaMatrix::g5(){
 360
 361   static EvtGammaMatrix g;
 362   static int first=1;
 363
 364   if (first){
 365     first = 0;
 366     int i,j;
 367
 368     for(i=0;i<4;i++){
 369       for(j=0;j<4;j++){
 370         g._gamma[i][j]=EvtComplex(0.0,0.0);
 371       }
 372     }
 373
 374     g._gamma[0][2]=EvtComplex(1.0,0.0);
 375     g._gamma[1][3]=EvtComplex(1.0,0.0);
 376     g._gamma[2][0]=EvtComplex(1.0,0.0);
 377     g._gamma[3][1]=EvtComplex(1.0,0.0);
 378   }
 379
 380   return g;
 381
 382 }
 383
 384
 385
 386
 387 const EvtGammaMatrix& EvtGammaMatrix::g(int index) {
 388   switch (index) {
 389   case 0:
 390     return g0();
 391   case 1:
 392     return g1();
 393   case 2:
 394     return g2();
 395   case 3:
 396     return g3();
 397   case 5:
 398     return g5();
 399   default:
 400     report(ERROR, "EvtGen") << "Invalid index for four vector: " << index << endl;
 401     exit(-2);
 402   }
 403 }
 404
 405
 406
 407 const EvtGammaMatrix& EvtGammaMatrix::v0(){
 408
 409   static EvtGammaMatrix g;
 410   static int first=1;
 411
 412   if (first){
 413     first = 0;
 414     int i,j;
 415
 416     for(i=0;i<4;i++){
 417       for(j=0;j<4;j++){
 418         g._gamma[i][j]=EvtComplex(0.0,0.0);
 419       }
 420     }
 421
 422     g._gamma[0][0]=EvtComplex(1.0,0.0);
 423     g._gamma[1][1]=EvtComplex(1.0,0.0);
 424     g._gamma[2][2]=EvtComplex(1.0,0.0);
 425     g._gamma[3][3]=EvtComplex(1.0,0.0);
 426   }
 427
 428   return g;
 429
 430 }
 431
 432
 433
 434
 435
 436 const EvtGammaMatrix& EvtGammaMatrix::v1(){
 437
 438   static EvtGammaMatrix g;
 439   static int first=1;
 440
 441   if (first){
 442     first = 0;
 443     int i,j;
 444
 445     for(i=0;i<4;i++){
 446       for(j=0;j<4;j++){
 447         g._gamma[i][j]=EvtComplex(0.0,0.0);
 448       }
 449     }
 450
 451     g._gamma[0][3]=EvtComplex(1.0,0.0);
 452     g._gamma[1][2]=EvtComplex(1.0,0.0);
 453     g._gamma[2][1]=EvtComplex(1.0,0.0);
 454     g._gamma[3][0]=EvtComplex(1.0,0.0);
 455   }
 456
 457   return g;
 458
 459 }
 460
 461
 462
 463
 464 const EvtGammaMatrix& EvtGammaMatrix::v2(){
 465
 466   static EvtGammaMatrix g;
 467   static int first=1;
 468
 469   if (first){
 470     first = 0;
 471     int i,j;
 472
 473     for(i=0;i<4;i++){
 474       for(j=0;j<4;j++){
 475         g._gamma[i][j]=EvtComplex(0.0,0.0);
 476       }
 477     }
 478
 479     g._gamma[0][3]=EvtComplex(0.0,-1.0);
 480     g._gamma[1][2]=EvtComplex(0.0,1.0);
 481     g._gamma[2][1]=EvtComplex(0.0,-1.0);
 482     g._gamma[3][0]=EvtComplex(0.0,1.0);
 483   }
 484
 485   return g;
 486
 487 }
 488
 489
 490
 491
 492 const EvtGammaMatrix& EvtGammaMatrix::v3(){
 493
 494   static EvtGammaMatrix g;
 495   static int first=1;
 496
 497   if (first){
 498     first = 0;
 499     int i,j;
 500
 501     for(i=0;i<4;i++){
 502       for(j=0;j<4;j++){
 503         g._gamma[i][j]=EvtComplex(0.0,0.0);
 504       }
 505     }
 506
 507     g._gamma[0][2]=EvtComplex(1.0,0.0);
 508     g._gamma[1][3]=EvtComplex(-1.0,0.0);
 509     g._gamma[2][0]=EvtComplex(1.0,0.0);
 510     g._gamma[3][1]=EvtComplex(-1.0,0.0);
 511   }
 512
 513   return g;
 514
 515 }
 516
 517
 518
 519
 520
 521 const EvtGammaMatrix& EvtGammaMatrix::id(){
 522
 523   static EvtGammaMatrix g;
 524   static int first=1;
 525
 526   if (first){
 527     first = 0;
 528     int i,j;
 529
 530     for(i=0;i<4;i++){
 531       for(j=0;j<4;j++){
 532         g._gamma[i][j]=EvtComplex(0.0,0.0);
 533       }
 534     }
 535
 536     g._gamma[0][0]=EvtComplex(1.0,0.0);
 537     g._gamma[1][1]=EvtComplex(1.0,0.0);
 538     g._gamma[2][2]=EvtComplex(1.0,0.0);
 539     g._gamma[3][3]=EvtComplex(1.0,0.0);
 540   }
 541
 542   return g;
 543
 544 }
 545
 546
 547
 548
 549 EvtGammaMatrix& EvtGammaMatrix::operator+=(const EvtGammaMatrix &g){
 550
 551   int i,j;
 552
 553   for(i=0;i<4;i++){
 554     for(j=0;j<4;j++){
 555       _gamma[i][j]+=g._gamma[i][j];
 556     }
 557   }
 558   return *this;
 559 }
 560
 561
 562
 563
 564
 565 EvtGammaMatrix& EvtGammaMatrix::operator-=(const EvtGammaMatrix &g){
 566
 567   int i,j;
 568
 569   for(i=0;i<4;i++){
 570     for(j=0;j<4;j++){
 571       _gamma[i][j]-=g._gamma[i][j];
 572     }
 573   }
 574   return *this;
 575 }
 576
 577
 578
 579 EvtGammaMatrix& EvtGammaMatrix::operator*=(const EvtGammaMatrix &g){
 580
 581   int i,j,k;
 582   EvtGammaMatrix temp;
 583
 584   for(i=0;i<4;i++){
 585     for(j=0;j<4;j++){
 586       temp._gamma[i][j]=EvtComplex(0.0,0.0);
 587       for(k=0;k<4;k++){
 588         temp._gamma[i][j]+=_gamma[i][k]*g._gamma[k][j];
 589       }
 590     }
 591   }
 592
 593   for(i=0;i<4;i++){
 594     for(j=0;j<4;j++){
 595        _gamma[i][j]=temp._gamma[i][j];
 596     }
 597   }
 598
 599   return *this;
 600 }
 601
 602
 603 EvtDiracSpinor operator*(const EvtGammaMatrix& g,const EvtDiracSpinor& d){
 604
 605   int i,j;
 606   EvtDiracSpinor temp;
 607
 608    for(i=0;i<4;i++){
 609      temp.set_spinor(i,EvtComplex(0.0,0.0));
 610      for(j=0;j<4;j++){
 611        temp.set_spinor(i,temp.get_spinor(i)+g._gamma[i][j]*d.get_spinor(j));
 612      }
 613    }
 614
 615    return temp;
 616 }
 617
 618 // upper index
 619 const EvtGammaMatrix& EvtGammaMatrix::sigmaUpper(unsigned int mu, unsigned int nu)
 620 {
 621     EvtGammaMatrix a, b;
 622     static const EvtTensor4C eta = EvtTensor4C::g(); //metric
 623     static EvtGammaMatrix sigma[4][4];
 624     static bool hasBeenCalled = false;
 625     if (!hasBeenCalled)
 626     {
 627         EvtComplex I(0, 1);
 628         for (int i=0; i<4; ++i)
 629             sigma[i][i].init(); // set to 0
 630
 631         EvtGammaMatrix s01 = I/2 * (g0()*g1() - g1()*g0());
 632         EvtGammaMatrix s02 = I/2 * (g0()*g2() - g2()*g0());
 633         EvtGammaMatrix s03 = I/2 * (g0()*g3() - g3()*g0());
 634         EvtGammaMatrix s12 = I/2 * (g1()*g2() - g2()*g1());
 635         EvtGammaMatrix s13 = I/2 * (g1()*g3() - g3()*g1());
 636         EvtGammaMatrix s23 = I/2 * (g2()*g3() - g3()*g2());
 637         sigma[0][1] = s01;
 638         sigma[1][0] = -1*s01;
 639         sigma[0][2] = s02;
 640         sigma[2][0] = -1*s02;
 641         sigma[0][3] = s03;
 642         sigma[3][0] = -1*s03;
 643         sigma[1][2] = s12;
 644         sigma[2][1] = -1*s12;
 645         sigma[1][3] = s13;
 646         sigma[3][1] = -1*s13;
 647         sigma[2][3] = s23;
 648         sigma[3][2] = -1*s23;
 649     }
 650     hasBeenCalled = true;
 651
 652     if (mu > 3 || nu > 3)
 653     {
 654         report(ERROR, "EvtSigmaTensor") << "Expected index between 0 and 3, but found " << nu << "!" << endl;
 655         assert(0);
 656     }
 657     return sigma[mu][nu];
 658
 659 }
 660
 661 const EvtGammaMatrix& EvtGammaMatrix::sigmaLower(unsigned int mu, unsigned int nu)
 662 {
 663     const EvtComplex I(0, 1);
 664     EvtGammaMatrix a, b;
 665     static EvtGammaMatrix sigma[4][4];
 666     static bool hasBeenCalled = false;
 667     static const EvtTensor4C eta = EvtTensor4C::g();
 668
 669     if (!hasBeenCalled) // has to be initialized only at the first call
 670     {
 671         // lower index
 672         for (int i=0; i<4; ++i)
 673         {
 674             a = eta.get(i, 0)*g0() + eta.get(i, 1)*g1() + eta.get(i, 2)*g2() + eta.get(i, 3)*g3();
 675             for (int j=0; j<4; ++j)
 676             {
 677                 b = eta.get(j, 0)*g0() + eta.get(j, 1)*g1() + eta.get(j, 2)*g2() + eta.get(j, 3)*g3();
 678                 sigma[i][j] = I/2 * (a*b - b*a);
 679             }
 680         }
 681     }
 682     return sigma[mu][nu];
 683 }
 684
 685
 686 EvtGammaMatrix EvtGenFunctions::slash(const EvtVector4C& p)
 687 {
 688     return EvtGammaMatrix::g0()*p.get(0) +
 689       EvtGammaMatrix::g1()*p.get(1) +
 690       EvtGammaMatrix::g2()*p.get(2) +
 691       EvtGammaMatrix::g3()*p.get(3);
 692 }
 693
 694 EvtGammaMatrix EvtGenFunctions::slash(const EvtVector4R& p)
 695 {
 696   return EvtGammaMatrix::g0()*p.get(0) +
 697     EvtGammaMatrix::g1()*p.get(1) +
 698     EvtGammaMatrix::g2()*p.get(2) +
 699     EvtGammaMatrix::g3()*p.get(3);
 700 }