EvtTensor4C.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: EvtTensor4C.cc
//
// Description: Implementation of tensor particles.
//
// Modification history:
//
//    DJL/RYD   September 25,1996           Module created
//
//------------------------------------------------------------------------
//
#include "EvtTensor4C.hh"
#include "EvtComplex.hh"
#include "EvtPatches.hh"
#include "EvtVector4C.hh"
#include <assert.h>
#include <iostream>
#include <math.h>
using std::endl;
using std::ostream;
 
EvtTensor4C::EvtTensor4C( const EvtTensor4C& t1 ) {
 
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] = t1.t[i][j]; }
  }
}
 
EvtTensor4C::~EvtTensor4C() {}
 
const EvtTensor4C& EvtTensor4C::g() {
 
  static EvtTensor4C g_metric( 1.0, -1.0, -1.0, -1.0 );
 
  return g_metric;
}
 
EvtTensor4C& EvtTensor4C::operator=( const EvtTensor4C& t1 ) {
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] = t1.t[i][j]; }
  }
  return *this;
}
 
EvtTensor4C EvtTensor4C::conj() const {
  EvtTensor4C temp;
 
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { temp.set( j, i, ::conj( t[i][j] ) ); }
  }
  return temp;
}
 
EvtTensor4C rotateEuler( const EvtTensor4C& rs, double alpha, double beta, double gamma ) {
 
  EvtTensor4C tmp( rs );
  tmp.applyRotateEuler( alpha, beta, gamma );
  return tmp;
}
 
EvtTensor4C boostTo( const EvtTensor4C& rs, const EvtVector4R p4 ) {
 
  EvtTensor4C tmp( rs );
  tmp.applyBoostTo( p4 );
  return tmp;
}
 
EvtTensor4C boostTo( const EvtTensor4C& rs, const EvtVector3R boost ) {
 
  EvtTensor4C tmp( rs );
  tmp.applyBoostTo( boost );
  return tmp;
}
 
void EvtTensor4C::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 EvtTensor4C::applyBoostTo( const EvtVector3R& boost ) {
 
  double     bx, by, bz, gamma, b2;
  double     lambda[4][4];
  EvtComplex tt[4][4];
 
  bx = boost.get( 0 );
  by = boost.get( 1 );
  bz = boost.get( 2 );
 
  double bxx = bx * bx;
  double byy = by * by;
  double bzz = bz * bz;
 
  b2 = bxx + byy + bzz;
 
  if ( b2 == 0.0 ) { return; }
 
  assert( b2 < 1.0 );
 
  gamma = 1.0 / sqrt( 1 - b2 );
 
  int i, j, k;
 
  if ( b2 == 0.0 ) { return; }
 
  lambda[0][0] = gamma;
  lambda[0][1] = gamma * bx;
  lambda[1][0] = gamma * bx;
  lambda[0][2] = gamma * by;
  lambda[2][0] = gamma * by;
  lambda[0][3] = gamma * bz;
  lambda[3][0] = gamma * bz;
 
  lambda[1][1] = 1.0 + ( gamma - 1.0 ) * bx * bx / b2;
  lambda[2][2] = 1.0 + ( gamma - 1.0 ) * by * by / b2;
  lambda[3][3] = 1.0 + ( gamma - 1.0 ) * bz * bz / b2;
 
  lambda[1][2] = ( gamma - 1.0 ) * bx * by / b2;
  lambda[2][1] = ( gamma - 1.0 ) * bx * by / b2;
 
  lambda[1][3] = ( gamma - 1.0 ) * bx * bz / b2;
  lambda[3][1] = ( gamma - 1.0 ) * bx * bz / b2;
 
  lambda[3][2] = ( gamma - 1.0 ) * bz * by / b2;
  lambda[2][3] = ( gamma - 1.0 ) * bz * by / b2;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ )
    {
      tt[i][j] = EvtComplex( 0.0 );
      for ( k = 0; k < 4; k++ ) { tt[i][j] = tt[i][j] + lambda[j][k] * t[i][k]; }
    }
  }
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ )
    {
      t[i][j] = EvtComplex( 0.0 );
      for ( k = 0; k < 4; k++ ) { t[i][j] = t[i][j] + lambda[i][k] * tt[k][j]; }
    }
  }
}
 
void EvtTensor4C::zero() {
  int i, j;
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] = EvtComplex( 0.0, 0.0 ); }
  }
}
 
ostream& operator<<( ostream& s, const EvtTensor4C& t ) {
 
  int i, j;
  s << endl;
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { s << t.t[i][j]; }
    s << endl;
  }
  return s;
}
 
void EvtTensor4C::setdiag( double g00, double g11, double g22, double g33 ) {
  t[0][0] = EvtComplex( g00 );
  t[1][1] = EvtComplex( g11 );
  t[2][2] = EvtComplex( g22 );
  t[3][3] = EvtComplex( g33 );
  t[0][1] = EvtComplex( 0.0 );
  t[0][2] = EvtComplex( 0.0 );
  t[0][3] = EvtComplex( 0.0 );
  t[1][0] = EvtComplex( 0.0 );
  t[1][2] = EvtComplex( 0.0 );
  t[1][3] = EvtComplex( 0.0 );
  t[2][0] = EvtComplex( 0.0 );
  t[2][1] = EvtComplex( 0.0 );
  t[2][3] = EvtComplex( 0.0 );
  t[3][0] = EvtComplex( 0.0 );
  t[3][1] = EvtComplex( 0.0 );
  t[3][2] = EvtComplex( 0.0 );
}
 
EvtTensor4C& EvtTensor4C::operator+=( const EvtTensor4C& t2 ) {
 
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] += t2.get( i, j ); }
  }
  return *this;
}
 
EvtTensor4C& EvtTensor4C::operator-=( const EvtTensor4C& t2 ) {
 
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] -= t2.get( i, j ); }
  }
  return *this;
}
 
EvtTensor4C& EvtTensor4C::operator*=( const EvtComplex& c ) {
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] *= c; }
  }
  return *this;
}
 
EvtTensor4C operator*( const EvtTensor4C& t1, const EvtComplex& c ) {
 
  return EvtTensor4C( t1 ) *= c;
}
 
EvtTensor4C operator*( const EvtComplex& c, const EvtTensor4C& t1 ) {
 
  return EvtTensor4C( t1 ) *= c;
}
 
EvtTensor4C& EvtTensor4C::operator*=( double d ) {
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] *= EvtComplex( d, 0.0 ); }
  }
  return *this;
}
 
EvtTensor4C operator*( const EvtTensor4C& t1, double d ) {
 
  return EvtTensor4C( t1 ) *= EvtComplex( d, 0.0 );
}
 
EvtTensor4C operator*( double d, const EvtTensor4C& t1 ) {
 
  return EvtTensor4C( t1 ) *= EvtComplex( d, 0.0 );
}
 
EvtComplex cont( const EvtTensor4C& t1, const EvtTensor4C& t2 ) {
 
  EvtComplex sum( 0.0, 0.0 );
  int        i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { sum += t1.t[i][j] * t2.t[i][j]; }
  }
 
  return sum;
}
 
EvtTensor4C directProd( const EvtVector4C& c1, const EvtVector4C& c2 ) {
  EvtTensor4C temp;
  int         i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { temp.set( i, j, c1.get( i ) * c2.get( j ) ); }
  }
  return temp;
}
 
EvtTensor4C directProd( const EvtVector4C& c1, const EvtVector4R& c2 ) {
  EvtTensor4C temp;
  int         i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { temp.set( i, j, c1.get( i ) * c2.get( j ) ); }
  }
  return temp;
}
 
EvtTensor4C directProd( const EvtVector4R& c1, const EvtVector4R& c2 ) {
 
  EvtTensor4C temp;
  int         i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { temp.t[i][j] = EvtComplex( c1.get( i ) * c2.get( j ), 0.0 ); }
  }
  return temp;
}
 
EvtTensor4C& EvtTensor4C::addDirProd( const EvtVector4R& p1, const EvtVector4R& p2 ) {
 
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { t[i][j] += p1.get( i ) * p2.get( j ); }
  }
  return *this;
}
 
EvtTensor4C dual( const EvtTensor4C& t2 ) {
 
  EvtTensor4C temp;
 
  temp.set( 0, 0, EvtComplex( 0.0, 0.0 ) );
  temp.set( 1, 1, EvtComplex( 0.0, 0.0 ) );
  temp.set( 2, 2, EvtComplex( 0.0, 0.0 ) );
  temp.set( 3, 3, EvtComplex( 0.0, 0.0 ) );
 
  temp.set( 0, 1, t2.get( 3, 2 ) - t2.get( 2, 3 ) );
  temp.set( 0, 2, -t2.get( 3, 1 ) + t2.get( 1, 3 ) );
  temp.set( 0, 3, t2.get( 2, 1 ) - t2.get( 1, 2 ) );
 
  temp.set( 1, 2, -t2.get( 3, 0 ) + t2.get( 0, 3 ) );
  temp.set( 1, 3, t2.get( 2, 0 ) - t2.get( 0, 2 ) );
 
  temp.set( 2, 3, -t2.get( 1, 0 ) + t2.get( 0, 1 ) );
 
  temp.set( 1, 0, -temp.get( 0, 1 ) );
  temp.set( 2, 0, -temp.get( 0, 2 ) );
  temp.set( 3, 0, -temp.get( 0, 3 ) );
 
  temp.set( 2, 1, -temp.get( 1, 2 ) );
  temp.set( 3, 1, -temp.get( 1, 3 ) );
 
  temp.set( 3, 2, -temp.get( 2, 3 ) );
 
  return temp;
}
 
EvtTensor4C conj( const EvtTensor4C& t2 ) {
  EvtTensor4C temp;
 
  int i, j;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ ) { temp.set( i, j, ::conj( ( t2.get( i, j ) ) ) ); }
  }
 
  return temp;
}
 
EvtTensor4C cont22( const EvtTensor4C& t1, const EvtTensor4C& t2 ) {
  EvtTensor4C temp;
 
  int        i, j;
  EvtComplex c;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ )
    {
      c = t1.get( i, 0 ) * t2.get( j, 0 ) - t1.get( i, 1 ) * t2.get( j, 1 ) -
          t1.get( i, 2 ) * t2.get( j, 2 ) - t1.get( i, 3 ) * t2.get( j, 3 );
      temp.set( i, j, c );
    }
  }
 
  return temp;
}
 
EvtTensor4C cont11( const EvtTensor4C& t1, const EvtTensor4C& t2 ) {
  EvtTensor4C temp;
 
  int        i, j;
  EvtComplex c;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ )
    {
      c = t1.get( 0, i ) * t2.get( 0, j ) - t1.get( 1, i ) * t2.get( 1, j ) -
          t1.get( 2, i ) * t2.get( 2, j ) - t1.get( 3, i ) * t2.get( 3, j );
      temp.set( i, j, c );
    }
  }
 
  return temp;
}
 
EvtVector4C EvtTensor4C::cont1( const EvtVector4C& v4 ) const {
  EvtVector4C temp;
 
  int i;
 
  for ( i = 0; i < 4; i++ )
  {
    temp.set( i, t[0][i] * v4.get( 0 ) - t[1][i] * v4.get( 1 ) - t[2][i] * v4.get( 2 ) -
                     t[3][i] * v4.get( 3 ) );
  }
 
  return temp;
}
 
EvtVector4C EvtTensor4C::cont2( const EvtVector4C& v4 ) const {
  EvtVector4C temp;
 
  int i;
 
  for ( i = 0; i < 4; i++ )
  {
    temp.set( i, t[i][0] * v4.get( 0 ) - t[i][1] * v4.get( 1 ) - t[i][2] * v4.get( 2 ) -
                     t[i][3] * v4.get( 3 ) );
  }
 
  return temp;
}
 
EvtVector4C EvtTensor4C::cont1( const EvtVector4R& v4 ) const {
  EvtVector4C temp;
 
  int i;
 
  for ( i = 0; i < 4; i++ )
  {
    temp.set( i, t[0][i] * v4.get( 0 ) - t[1][i] * v4.get( 1 ) - t[2][i] * v4.get( 2 ) -
                     t[3][i] * v4.get( 3 ) );
  }
 
  return temp;
}
 
EvtVector4C EvtTensor4C::cont2( const EvtVector4R& v4 ) const {
  EvtVector4C temp;
 
  int i;
 
  for ( i = 0; i < 4; i++ )
  {
    temp.set( i, t[i][0] * v4.get( 0 ) - t[i][1] * v4.get( 1 ) - t[i][2] * v4.get( 2 ) -
                     t[i][3] * v4.get( 3 ) );
  }
 
  return temp;
}
 
void EvtTensor4C::applyRotateEuler( double phi, double theta, double ksi ) {
 
  EvtComplex tt[4][4];
  double     sp, st, sk, cp, ct, ck;
  double     lambda[4][4];
 
  sp = sin( phi );
  st = sin( theta );
  sk = sin( ksi );
  cp = cos( phi );
  ct = cos( theta );
  ck = cos( ksi );
 
  lambda[0][0] = 1.0;
  lambda[0][1] = 0.0;
  lambda[1][0] = 0.0;
  lambda[0][2] = 0.0;
  lambda[2][0] = 0.0;
  lambda[0][3] = 0.0;
  lambda[3][0] = 0.0;
 
  lambda[1][1] = ck * ct * cp - sk * sp;
  lambda[1][2] = -sk * ct * cp - ck * sp;
  lambda[1][3] = st * cp;
 
  lambda[2][1] = ck * ct * sp + sk * cp;
  lambda[2][2] = -sk * ct * sp + ck * cp;
  lambda[2][3] = st * sp;
 
  lambda[3][1] = -ck * st;
  lambda[3][2] = sk * st;
  lambda[3][3] = ct;
 
  int i, j, k;
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ )
    {
      tt[i][j] = EvtComplex( 0.0 );
      for ( k = 0; k < 4; k++ ) { tt[i][j] += lambda[j][k] * t[i][k]; }
    }
  }
 
  for ( i = 0; i < 4; i++ )
  {
    for ( j = 0; j < 4; j++ )
    {
      t[i][j] = EvtComplex( 0.0 );
      for ( k = 0; k < 4; k++ ) { t[i][j] += lambda[i][k] * tt[k][j]; }
    }
  }
}