BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtGammaMatrix.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtGammaMatrix.cc
//
// Description: Make gamma matrices availible for the calc. of amplitudes, etc.
//
// Modification history:
//
//    DJL/RYD     September 25, 1996         Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include <math.h>
#include <assert.h>
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtGammaMatrix.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
using std::endl;
using std::ostream;
 
EvtGammaMatrix::EvtGammaMatrix(){
  int i,j;
 
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      gamma[i][j]=EvtComplex(0.0,0.0);
    }
  }
}
 
EvtGammaMatrix operator*(const EvtGammaMatrix& g, const EvtComplex& c)
{
    return c*g;
}
 
 
EvtGammaMatrix operator*(const EvtComplex& c,const EvtGammaMatrix& g){
  int i,j;
 
  EvtGammaMatrix temp;
 
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      temp.gamma[i][j]=g.gamma[i][j]*c;
    }
  }
 
  return temp;
 
}
 
 
ostream& operator<<(ostream& s, const EvtGammaMatrix& g){
 
 
  s<<"["<<g.gamma[0][0]<<","<<g.gamma[0][1]<<","<<g.gamma[0][2]<<","<<g.gamma[0][3]<<"]"<<endl;
  s<<"["<<g.gamma[1][0]<<","<<g.gamma[1][1]<<","<<g.gamma[1][2]<<","<<g.gamma[1][3]<<"]"<<endl;
  s<<"["<<g.gamma[2][0]<<","<<g.gamma[2][1]<<","<<g.gamma[2][2]<<","<<g.gamma[2][3]<<"]"<<endl;
  s<<"["<<g.gamma[3][0]<<","<<g.gamma[3][1]<<","<<g.gamma[3][2]<<","<<g.gamma[3][3]<<"]"<<endl;
 
  return s;
 
}
 
 
 
EvtGammaMatrix::EvtGammaMatrix(const EvtGammaMatrix& gm){
  int i,j;
  
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      gamma[i][j]=gm.gamma[i][j];
    }
  }
}
 
EvtGammaMatrix::~EvtGammaMatrix() {}
 
EvtGammaMatrix& EvtGammaMatrix::operator=(const EvtGammaMatrix& gm){
  int i,j;
  
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      gamma[i][j]=gm.gamma[i][j];
    }
  }
  return *this;
}
 
void EvtGammaMatrix::init(){
  int i,j;
  
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      gamma[i][j]=EvtComplex(0.0,0.0);
    }
  }
}
 
const EvtGammaMatrix& EvtGammaMatrix::va0(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    g.gamma[0][0]=EvtComplex(1.0,0.0);
    g.gamma[0][1]=EvtComplex(0.0,0.0);
    g.gamma[0][2]=EvtComplex(-1.0,0.0);
    g.gamma[0][3]=EvtComplex(0.0,0.0);
    g.gamma[1][0]=EvtComplex(0.0,0.0);
    g.gamma[1][1]=EvtComplex(1.0,0.0);
    g.gamma[1][2]=EvtComplex(0.0,0.0);
    g.gamma[1][3]=EvtComplex(-1.0,0.0);
    g.gamma[2][0]=EvtComplex(-1.0,0.0);
    g.gamma[2][1]=EvtComplex(0.0,0.0);
    g.gamma[2][2]=EvtComplex(1.0,0.0);
    g.gamma[2][3]=EvtComplex(0.0,0.0);
    g.gamma[3][0]=EvtComplex(0.0,0.0);
    g.gamma[3][1]=EvtComplex(-1.0,0.0);
    g.gamma[3][2]=EvtComplex(0.0,0.0);
    g.gamma[3][3]=EvtComplex(1.0,0.0);
  }
 
  return g;
 
}
 
 
const EvtGammaMatrix& EvtGammaMatrix::va1(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    g.gamma[0][0]=EvtComplex(0.0,0.0);
    g.gamma[0][1]=EvtComplex(-1.0,0.0);
    g.gamma[0][2]=EvtComplex(0.0,0.0);
    g.gamma[0][3]=EvtComplex(1.0,0.0);
    g.gamma[1][0]=EvtComplex(-1.0,0.0);
    g.gamma[1][1]=EvtComplex(0.0,0.0);
    g.gamma[1][2]=EvtComplex(1.0,0.0);
    g.gamma[1][3]=EvtComplex(0.0,0.0);
    g.gamma[2][0]=EvtComplex(0.0,0.0);
    g.gamma[2][1]=EvtComplex(1.0,0.0);
    g.gamma[2][2]=EvtComplex(0.0,0.0);
    g.gamma[2][3]=EvtComplex(-1.0,0.0);
    g.gamma[3][0]=EvtComplex(1.0,0.0);
    g.gamma[3][1]=EvtComplex(0.0,0.0);
    g.gamma[3][2]=EvtComplex(-1.0,0.0);
    g.gamma[3][3]=EvtComplex(0.0,0.0);
  }
 
  return g;
 
}
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::va2(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    g.gamma[0][0]=EvtComplex(0.0,0.0);
    g.gamma[0][1]=EvtComplex(0.0,1.0);
    g.gamma[0][2]=EvtComplex(0.0,0.0);
    g.gamma[0][3]=EvtComplex(0.0,-1.0);
    g.gamma[1][0]=EvtComplex(0.0,-1.0);
    g.gamma[1][1]=EvtComplex(0.0,0.0);
    g.gamma[1][2]=EvtComplex(0.0,1.0);
    g.gamma[1][3]=EvtComplex(0.0,0.0);
    g.gamma[2][0]=EvtComplex(0.0,0.0);
    g.gamma[2][1]=EvtComplex(0.0,-1.0);
    g.gamma[2][2]=EvtComplex(0.0,0.0);
    g.gamma[2][3]=EvtComplex(0.0,1.0);
    g.gamma[3][0]=EvtComplex(0.0,1.0);
    g.gamma[3][1]=EvtComplex(0.0,0.0);
    g.gamma[3][2]=EvtComplex(0.0,-1.0);
    g.gamma[3][3]=EvtComplex(0.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::va3(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    g.gamma[0][0]=EvtComplex(-1.0,0.0);
    g.gamma[0][1]=EvtComplex(0.0,0.0);
    g.gamma[0][2]=EvtComplex(1.0,0.0);
    g.gamma[0][3]=EvtComplex(0.0,0.0);
    g.gamma[1][0]=EvtComplex(0.0,0.0);
    g.gamma[1][1]=EvtComplex(1.0,0.0);
    g.gamma[1][2]=EvtComplex(0.0,0.0);
    g.gamma[1][3]=EvtComplex(-1.0,0.0);
    g.gamma[2][0]=EvtComplex(1.0,0.0);
    g.gamma[2][1]=EvtComplex(0.0,0.0);
    g.gamma[2][2]=EvtComplex(-1.0,0.0);
    g.gamma[2][3]=EvtComplex(0.0,0.0);
    g.gamma[3][0]=EvtComplex(0.0,0.0);
    g.gamma[3][1]=EvtComplex(-1.0,0.0);
    g.gamma[3][2]=EvtComplex(0.0,0.0);
    g.gamma[3][3]=EvtComplex(1.0,0.0);
  }
 
  return g;
  
}
 
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::g0(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
 
    first=0;
 
    int i,j;
  
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][0]=EvtComplex(1.0,0.0);
    g.gamma[1][1]=EvtComplex(1.0,0.0);
    g.gamma[2][2]=EvtComplex(-1.0,0.0);
    g.gamma[3][3]=EvtComplex(-1.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::g1(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    first=0;
    int i,j;
    
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][3]=EvtComplex(1.0,0.0);
    g.gamma[1][2]=EvtComplex(1.0,0.0);
    g.gamma[2][1]=EvtComplex(-1.0,0.0);
    g.gamma[3][0]=EvtComplex(-1.0,0.0);
  }
  
  return g;
 
}
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::g2(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    first=0;
    int i,j;
    
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][3]=EvtComplex(0.0,-1.0);
    g.gamma[1][2]=EvtComplex(0.0,1.0);
    g.gamma[2][1]=EvtComplex(0.0,1.0);
    g.gamma[3][0]=EvtComplex(0.0,-1.0);
  }
  
  return g;
 
}
 
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::g3(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
    first=0;
    int i,j;
    
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][2]=EvtComplex(1.0,0.0);
    g.gamma[1][3]=EvtComplex(-1.0,0.0);
    g.gamma[2][0]=EvtComplex(-1.0,0.0);
    g.gamma[3][1]=EvtComplex(1.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::g5(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
 
    int i,j;
    
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][2]=EvtComplex(1.0,0.0);
    g.gamma[1][3]=EvtComplex(1.0,0.0);
    g.gamma[2][0]=EvtComplex(1.0,0.0);
    g.gamma[3][1]=EvtComplex(1.0,0.0);
  }
 
  return g;
 
}
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::v0(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
 
    int i,j;
  
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][0]=EvtComplex(1.0,0.0);
    g.gamma[1][1]=EvtComplex(1.0,0.0);
    g.gamma[2][2]=EvtComplex(1.0,0.0);
    g.gamma[3][3]=EvtComplex(1.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::v1(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
 
    int i,j;
    
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][3]=EvtComplex(1.0,0.0);
    g.gamma[1][2]=EvtComplex(1.0,0.0);
    g.gamma[2][1]=EvtComplex(1.0,0.0);
    g.gamma[3][0]=EvtComplex(1.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::v2(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
  
    int i,j;
 
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][3]=EvtComplex(0.0,-1.0);
    g.gamma[1][2]=EvtComplex(0.0,1.0);
    g.gamma[2][1]=EvtComplex(0.0,-1.0);
    g.gamma[3][0]=EvtComplex(0.0,1.0);
  }
 
  return g;
 
}
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::v3(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
 
    int i,j;
  
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][2]=EvtComplex(1.0,0.0);
    g.gamma[1][3]=EvtComplex(-1.0,0.0);
    g.gamma[2][0]=EvtComplex(1.0,0.0);
    g.gamma[3][1]=EvtComplex(-1.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
 
const EvtGammaMatrix& EvtGammaMatrix::id(){
 
  static EvtGammaMatrix g;
  static int first=1;
 
  if (first){
 
    int i,j;
    
    for(i=0;i<4;i++){
      for(j=0;j<4;j++){
    g.gamma[i][j]=EvtComplex(0.0,0.0);
      }
    }
    
    g.gamma[0][0]=EvtComplex(1.0,0.0);
    g.gamma[1][1]=EvtComplex(1.0,0.0);
    g.gamma[2][2]=EvtComplex(1.0,0.0);
    g.gamma[3][3]=EvtComplex(1.0,0.0);
  }
 
  return g;
 
}
 
 
 
 
EvtGammaMatrix& EvtGammaMatrix::operator+=(const EvtGammaMatrix &g){
 
  int i,j;
  
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      gamma[i][j]+=g.gamma[i][j];
    }
  }
  return *this;
}
 
 
 
 
 
EvtGammaMatrix& EvtGammaMatrix::operator-=(const EvtGammaMatrix &g){
 
  int i,j;
  
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      gamma[i][j]-=g.gamma[i][j];
    }
  }
  return *this;
}
 
 
 
EvtGammaMatrix& EvtGammaMatrix::operator*=(const EvtGammaMatrix &g){
 
  int i,j,k;
  EvtGammaMatrix temp;
 
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      temp.gamma[i][j]=EvtComplex(0.0,0.0);
      for(k=0;k<4;k++){
    temp.gamma[i][j]+=gamma[i][k]*g.gamma[k][j];
      }
    }
  }
 
  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
       gamma[i][j]=temp.gamma[i][j];
    }
  }
 
  return *this;
}
 
 
EvtDiracSpinor operator*(const EvtGammaMatrix& g,const EvtDiracSpinor& d){
 
  int i,j;
  EvtDiracSpinor temp;
  
   for(i=0;i<4;i++){
     temp.set_spinor(i,EvtComplex(0.0,0.0));
     for(j=0;j<4;j++){
       temp.set_spinor(i,temp.get_spinor(i)+g.gamma[i][j]*d.get_spinor(j));      
     }
   }
   
   return temp;
}  
 
 
EvtComplex operator*(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 
  int i;
  EvtComplex temp;
  
  temp=EvtComplex(0.0,0.0);
  
  for(i=0;i<4;i++){
    temp+=::conj(d.get_spinor(i))*dp.get_spinor(i);
  }
  return temp;
}  
 
 
// upper index
const EvtGammaMatrix& EvtGammaMatrix::sigmaUpper(unsigned int mu, unsigned int nu)
{
    EvtGammaMatrix a, b;
    static const EvtTensor4C eta = EvtTensor4C::g(); //metric
    static EvtGammaMatrix sigma[4][4];
    static bool hasBeenCalled = false;
    if (!hasBeenCalled)
    {
        EvtComplex I(0, 1);
        for (int i=0; i<4; ++i)
            sigma[i][i].init(); // set to 0
        
        EvtGammaMatrix s01 = I/2 * (g0()*g1() - g1()*g0());
        EvtGammaMatrix s02 = I/2 * (g0()*g2() - g2()*g0());
        EvtGammaMatrix s03 = I/2 * (g0()*g3() - g3()*g0());
        EvtGammaMatrix s12 = I/2 * (g1()*g2() - g2()*g1());
        EvtGammaMatrix s13 = I/2 * (g1()*g3() - g3()*g1());
        EvtGammaMatrix s23 = I/2 * (g2()*g3() - g3()*g2());
        sigma[0][1] = s01;
        sigma[1][0] = -1*s01;
        sigma[0][2] = s02;
        sigma[2][0] = -1*s02;
        sigma[0][3] = s03;
        sigma[3][0] = -1*s03;
        sigma[1][2] = s12;
        sigma[2][1] = -1*s12;
        sigma[1][3] = s13;
        sigma[3][1] = -1*s13;
        sigma[2][3] = s23;
        sigma[3][2] = -1*s23;
    }
    hasBeenCalled = true;
        
    if (mu > 3 || nu > 3)
    {
        report(ERROR, "EvtSigmaTensor") << "Expected index between 0 and 3, but found " << nu << "!" << endl;
        assert(0);
    }
    return sigma[mu][nu];
    
}
 
const EvtGammaMatrix& EvtGammaMatrix::sigmaLower(unsigned int mu, unsigned int nu)
{
    const EvtComplex I(0, 1);
    EvtGammaMatrix a, b;
    static EvtGammaMatrix sigma[4][4];
    static bool hasBeenCalled = false;
    static const EvtTensor4C eta = EvtTensor4C::g();
    
    if (!hasBeenCalled) // has to be initialized only at the first call
    {
        // lower index
        for (int i=0; i<4; ++i)
        {
            a = eta.get(i, 0)*g0() + eta.get(i, 1)*g1() + eta.get(i, 2)*g2() + eta.get(i, 3)*g3();
            for (int j=0; j<4; ++j)
            {
                b = eta.get(j, 0)*g0() + eta.get(j, 1)*g1() + eta.get(j, 2)*g2() + eta.get(j, 3)*g3();
                sigma[i][j] = I/2 * (a*b - b*a);
            }
        }
    }
    return sigma[mu][nu];    
}
 
 
EvtGammaMatrix slash(const EvtVector4C& p)
{
    return EvtGammaMatrix::g0()*p.get(0) + EvtGammaMatrix::g1()*p.get(1) + EvtGammaMatrix::g2()*p.get(2) + EvtGammaMatrix::g3()*p.get(3);
}