EvtDiracSpinor.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: EvtDiracSpinor.cc
//
// Description:  Class to describe (EvtDiracParticle) spinors.
//
// Modification history:
//
//    DJL/RYD     September 25, 1996         Module created
//
//------------------------------------------------------------------------
//
#include "EvtDiracSpinor.hh"
#include "EvtComplex.hh"
#include "EvtGammaMatrix.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include "EvtTensor4C.hh"
#include "EvtVector4C.hh"
#include <assert.h>
#include <math.h>
using std::ostream;
 
EvtDiracSpinor::~EvtDiracSpinor() {}
 
EvtDiracSpinor::EvtDiracSpinor( const EvtComplex& sp0, const EvtComplex& sp1,
                                const EvtComplex& sp2, const EvtComplex& sp3 ) {
  set( sp0, sp1, sp2, sp3 );
}
 
void EvtDiracSpinor::set( const EvtComplex& sp0, const EvtComplex& sp1, const EvtComplex& sp2,
                          const EvtComplex& sp3 ) {
 
  spinor[0] = sp0;
  spinor[1] = sp1;
  spinor[2] = sp2;
  spinor[3] = sp3;
}
 
void EvtDiracSpinor::set_spinor( int i, const EvtComplex& sp ) { spinor[i] = sp; }
 
ostream& operator<<( ostream& s, const EvtDiracSpinor& sp ) {
 
  s << "[" << sp.spinor[0] << "," << sp.spinor[1] << "," << sp.spinor[2] << "," << sp.spinor[3]
    << "]";
  return s;
}
 
const EvtComplex& EvtDiracSpinor::get_spinor( int i ) const { return spinor[i]; }
 
EvtDiracSpinor rotateEuler( const EvtDiracSpinor& sp, double alpha, double beta,
                            double gamma ) {
 
  EvtDiracSpinor tmp( sp );
  tmp.applyRotateEuler( alpha, beta, gamma );
  return tmp;
}
 
EvtDiracSpinor boostTo( const EvtDiracSpinor& sp, const EvtVector4R p4 ) {
 
  EvtDiracSpinor tmp( sp );
  tmp.applyBoostTo( p4 );
  return tmp;
}
 
EvtDiracSpinor boostTo( const EvtDiracSpinor& sp, const EvtVector3R boost ) {
 
  EvtDiracSpinor tmp( sp );
  tmp.applyBoostTo( boost );
  return tmp;
}
 
void EvtDiracSpinor::applyBoostTo( const EvtVector4R& p4 ) {
 
  double e = p4.get( 0 );
 
  EvtVector3R boost( p4.get( 1 ) / e, p4.get( 2 ) / e, p4.get( 3 ) / e );
 
  applyBoostTo( boost );
 
  return;
}
 
void EvtDiracSpinor::applyBoostTo( const EvtVector3R& boost ) {
 
  double     bx, by, bz, gamma, b2, f1, f2;
  EvtComplex spinorp[4];
 
  bx = boost.get( 0 );
  by = boost.get( 1 );
  bz = boost.get( 2 );
  b2 = bx * bx + by * by + bz * bz;
 
  if ( b2 == 0.0 ) { return; }
 
  assert( b2 < 1.0 );
 
  gamma = 1.0 / sqrt( 1 - b2 );
 
  f1 = sqrt( ( gamma + 1.0 ) / 2.0 );
  f2 = f1 * gamma / ( gamma + 1.0 );
 
  spinorp[0] = f1 * spinor[0] + f2 * bz * spinor[2] + f2 * EvtComplex( bx, -by ) * spinor[3];
  spinorp[1] = f1 * spinor[1] + f2 * EvtComplex( bx, by ) * spinor[2] - f2 * bz * spinor[3];
  spinorp[2] = f2 * bz * spinor[0] + f2 * EvtComplex( bx, -by ) * spinor[1] + f1 * spinor[2];
  spinorp[3] = f2 * EvtComplex( bx, by ) * spinor[0] - f2 * bz * spinor[1] + f1 * spinor[3];
 
  spinor[0] = spinorp[0];
  spinor[1] = spinorp[1];
  spinor[2] = spinorp[2];
  spinor[3] = spinorp[3];
 
  return;
}
 
void EvtDiracSpinor::applyRotateEuler( double alpha, double beta, double gamma ) {
 
  EvtComplex retVal[4];
 
  double cb2   = cos( 0.5 * beta );
  double sb2   = sin( 0.5 * beta );
  double capg2 = cos( 0.5 * ( alpha + gamma ) );
  double camg2 = cos( 0.5 * ( alpha - gamma ) );
  double sapg2 = sin( 0.5 * ( alpha + gamma ) );
  double samg2 = sin( 0.5 * ( alpha - gamma ) );
 
  EvtComplex m11( cb2 * capg2, -cb2 * sapg2 );
  EvtComplex m12( -sb2 * camg2, sb2 * samg2 );
  EvtComplex m21( sb2 * camg2, sb2 * samg2 );
  EvtComplex m22( cb2 * capg2, cb2 * sapg2 );
 
  retVal[0] = m11 * spinor[0] + m12 * spinor[1];
  retVal[1] = m21 * spinor[0] + m22 * spinor[1];
  retVal[2] = m11 * spinor[2] + m12 * spinor[3];
  retVal[3] = m21 * spinor[2] + m22 * spinor[3];
 
  spinor[0] = retVal[0];
  spinor[1] = retVal[1];
  spinor[2] = retVal[2];
  spinor[3] = retVal[3];
 
  return;
}
 
EvtDiracSpinor EvtDiracSpinor::conj() const {
 
  EvtDiracSpinor sp;
 
  for ( int i = 0; i < 4; i++ ) sp.set_spinor( i, ::conj( spinor[i] ) );
 
  return sp;
}
 
EvtVector4C EvtLeptonVACurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) {
 
  // Old code; below is a new specialized code that does it more efficiently.
  // EvtGammaMatrix mat;
  // EvtVector4C temp;
  // mat.va0();
  // temp.set(0,d*(mat*dp));
  // mat.va1();
  // temp.set(1,d*(mat*dp));
  // mat.va2();
  // temp.set(2,d*(mat*dp));
  // mat.va3();
  // temp.set(3,d*(mat*dp));
  // return temp;
 
  EvtComplex u02 = ::conj( d.spinor[0] - d.spinor[2] );
  EvtComplex u13 = ::conj( d.spinor[1] - d.spinor[3] );
 
  EvtComplex v02 = dp.spinor[0] - dp.spinor[2];
  EvtComplex v13 = dp.spinor[1] - dp.spinor[3];
 
  EvtComplex a = u02 * v02;
  EvtComplex b = u13 * v13;
 
  EvtComplex c = u02 * v13;
  EvtComplex e = u13 * v02;
 
  return EvtVector4C( a + b, -( c + e ), EvtComplex( 0, 1 ) * ( c - e ), b - a );
}
 
EvtVector4C EvtLeptonVCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) {
 
  EvtVector4C temp;
 
  // no conjugate here; done in the multiplication
  // yes this is stupid and fooled me to for a long time (ryd)
 
  temp.set( 0, d * ( EvtGammaMatrix::v0() * dp ) );
  temp.set( 1, d * ( EvtGammaMatrix::v1() * dp ) );
  temp.set( 2, d * ( EvtGammaMatrix::v2() * dp ) );
  temp.set( 3, d * ( EvtGammaMatrix::v3() * dp ) );
 
  return temp;
}
 
EvtVector4C EvtLeptonACurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) {
 
  EvtVector4C temp;
 
  EvtGammaMatrix mat;
 
  // no conjugate here; done in the multiplication
  // yes this is stupid and fooled me to for a long time (ryd)
 
  mat = EvtGammaMatrix::v0() - EvtGammaMatrix::va0();
  temp.set( 0, d * ( mat * dp ) );
 
  mat = EvtGammaMatrix::v1() - EvtGammaMatrix::va1();
  temp.set( 1, d * ( mat * dp ) );
 
  mat = EvtGammaMatrix::v2() - EvtGammaMatrix::va2();
  temp.set( 2, d * ( mat * dp ) );
 
  mat = EvtGammaMatrix::v3() - EvtGammaMatrix::va3();
  temp.set( 3, d * ( mat * dp ) );
 
  return temp;
}
 
EvtComplex EvtLeptonSCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) {
 
  EvtComplex temp;
 
  // no conjugate here; done in the multiplication
  // yes this is stupid and fooled me to for a long time (ryd)
 
  temp = d * ( EvtGammaMatrix::g0() * dp );
 
  return temp;
}
 
EvtComplex EvtLeptonPCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) {
 
  EvtComplex temp;
 
  // no conjugate here; done in the multiplication
  // yes this is stupid and fooled me to for a long time (ryd)
  static EvtGammaMatrix m = EvtGammaMatrix::g0() * EvtGammaMatrix::g5();
  temp                    = d * ( m * dp );
 
  return temp;
}
 
EvtTensor4C EvtLeptonTCurrent( const EvtDiracSpinor& d, const EvtDiracSpinor& dp ) {
 
  EvtTensor4C temp;
  temp.zero();
  EvtComplex i2( 0, 0.5 );
 
  static EvtGammaMatrix mat01 =
      EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g1() -
                               EvtGammaMatrix::g1() * EvtGammaMatrix::g0() );
  static EvtGammaMatrix mat02 =
      EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g2() -
                               EvtGammaMatrix::g2() * EvtGammaMatrix::g0() );
  static EvtGammaMatrix mat03 =
      EvtGammaMatrix::g0() * ( EvtGammaMatrix::g0() * EvtGammaMatrix::g3() -
                               EvtGammaMatrix::g3() * EvtGammaMatrix::g0() );
  static EvtGammaMatrix mat12 =
      EvtGammaMatrix::g0() * ( EvtGammaMatrix::g1() * EvtGammaMatrix::g2() -
                               EvtGammaMatrix::g2() * EvtGammaMatrix::g1() );
  static EvtGammaMatrix mat13 =
      EvtGammaMatrix::g0() * ( EvtGammaMatrix::g1() * EvtGammaMatrix::g3() -
                               EvtGammaMatrix::g3() * EvtGammaMatrix::g1() );
  static EvtGammaMatrix mat23 =
      EvtGammaMatrix::g0() * ( EvtGammaMatrix::g2() * EvtGammaMatrix::g3() -
                               EvtGammaMatrix::g3() * EvtGammaMatrix::g2() );
 
  temp.set( 0, 1, i2 * ( d * ( mat01 * dp ) ) );
  temp.set( 1, 0, -temp.get( 0, 1 ) );
 
  temp.set( 0, 2, i2 * ( d * ( mat02 * dp ) ) );
  temp.set( 2, 0, -temp.get( 0, 2 ) );
 
  temp.set( 0, 3, i2 * ( d * ( mat03 * dp ) ) );
  temp.set( 3, 0, -temp.get( 0, 3 ) );
 
  temp.set( 1, 2, i2 * ( d * ( mat12 * dp ) ) );
  temp.set( 2, 1, -temp.get( 1, 2 ) );
 
  temp.set( 1, 3, i2 * ( d * ( mat13 * dp ) ) );
  temp.set( 3, 1, -temp.get( 1, 3 ) );
 
  temp.set( 2, 3, i2 * ( d * ( mat23 * dp ) ) );
  temp.set( 3, 2, -temp.get( 2, 3 ) );
 
  return temp;
}
 
EvtTensor4C EvtLeptonTg5Current( const EvtDiracSpinor& d,
                                 const EvtDiracSpinor& dp ) { // sigma^{mu,nu} gamma_5
 
  EvtTensor4C temp;
  temp.zero();
  EvtComplex i2( 0, 0.5 );
 
  static EvtGammaMatrix mat01 = EvtGammaMatrix::g0() *
                                ( EvtGammaMatrix::g0() * EvtGammaMatrix::g1() -
                                  EvtGammaMatrix::g1() * EvtGammaMatrix::g0() ) *
                                EvtGammaMatrix::g5();
  static EvtGammaMatrix mat02 = EvtGammaMatrix::g0() *
                                ( EvtGammaMatrix::g0() * EvtGammaMatrix::g2() -
                                  EvtGammaMatrix::g2() * EvtGammaMatrix::g0() ) *
                                EvtGammaMatrix::g5();
  static EvtGammaMatrix mat03 = EvtGammaMatrix::g0() *
                                ( EvtGammaMatrix::g0() * EvtGammaMatrix::g3() -
                                  EvtGammaMatrix::g3() * EvtGammaMatrix::g0() ) *
                                EvtGammaMatrix::g5();
  static EvtGammaMatrix mat12 = EvtGammaMatrix::g0() *
                                ( EvtGammaMatrix::g1() * EvtGammaMatrix::g2() -
                                  EvtGammaMatrix::g2() * EvtGammaMatrix::g1() ) *
                                EvtGammaMatrix::g5();
  static EvtGammaMatrix mat13 = EvtGammaMatrix::g0() *
                                ( EvtGammaMatrix::g1() * EvtGammaMatrix::g3() -
                                  EvtGammaMatrix::g3() * EvtGammaMatrix::g1() ) *
                                EvtGammaMatrix::g5();
  static EvtGammaMatrix mat23 = EvtGammaMatrix::g0() *
                                ( EvtGammaMatrix::g2() * EvtGammaMatrix::g3() -
                                  EvtGammaMatrix::g3() * EvtGammaMatrix::g2() ) *
                                EvtGammaMatrix::g5();
 
  temp.set( 0, 1, i2 * ( d * ( mat01 * dp ) ) );
  temp.set( 1, 0, -temp.get( 0, 1 ) );
 
  temp.set( 0, 2, i2 * ( d * ( mat02 * dp ) ) );
  temp.set( 2, 0, -temp.get( 0, 2 ) );
 
  temp.set( 0, 3, i2 * ( d * ( mat03 * dp ) ) );
  temp.set( 3, 0, -temp.get( 0, 3 ) );
 
  temp.set( 1, 2, i2 * ( d * ( mat12 * dp ) ) );
  temp.set( 2, 1, -temp.get( 1, 2 ) );
 
  temp.set( 1, 3, i2 * ( d * ( mat13 * dp ) ) );
  temp.set( 3, 1, -temp.get( 1, 3 ) );
 
  temp.set( 2, 3, i2 * ( d * ( mat23 * dp ) ) );
  temp.set( 3, 2, -temp.get( 2, 3 ) );
 
  return temp;
}
 
EvtDiracSpinor operator*( const EvtComplex& c, const EvtDiracSpinor& d ) {
  EvtDiracSpinor result;
  result.spinor[0] = c * d.spinor[0];
  result.spinor[1] = c * d.spinor[1];
  result.spinor[2] = c * d.spinor[2];
  result.spinor[3] = c * d.spinor[3];
  return result;
}
 
EvtDiracSpinor EvtDiracSpinor::adjoint() const {
  EvtDiracSpinor d  = this->conj(); // first conjugate, then multiply with gamma0
  EvtGammaMatrix g0 = EvtGammaMatrix::g0();
  EvtDiracSpinor result; // automatically initialized to 0
 
  for ( int i = 0; i < 4; ++i )
    for ( int j = 0; j < 4; ++j ) result.spinor[i] += d.spinor[j] * g0.gamma[i][j];
 
  return result;
}