EvtSpinAmp.cc

Go to the documentation of this file.

#include "EvtSpinAmp.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include <cstdlib>
 
using std::endl;
 
std::ostream& operator<<( std::ostream& os, const EvtSpinAmp& amp ) {
  vector<int> index = amp.iterinit();
 
  os << ":";
  do {
    os << "<";
    for ( int i = 0; i < index.size() - 1; ++i ) { os << index[i]; }
    os << index[index.size() - 1] << ">" << amp( index ) << ":";
  } while ( amp.iterate( index ) );
 
  return os;
}
 
EvtSpinAmp operator*( const EvtComplex& real, const EvtSpinAmp& cont ) {
  EvtSpinAmp ret( cont );
 
  for ( int i = 0; i < ret._elem.size(); ++i ) { ret._elem[i] *= real; }
 
  return ret;
}
 
EvtSpinAmp operator*( const EvtSpinAmp& cont, const EvtComplex& real ) { return real * cont; }
 
EvtSpinAmp operator/( const EvtSpinAmp& cont, const EvtComplex& real ) {
  EvtSpinAmp ret( cont );
 
  for ( int i = 0; i < ret._elem.size(); ++i ) { ret._elem[i] /= real; }
 
  return ret;
}
 
vector<int> EvtSpinAmp::calctwospin( const vector<EvtSpinType::spintype>& type ) const {
  vector<int> twospin;
 
  for ( int i = 0; i < type.size(); ++i )
  { twospin.push_back( EvtSpinType::getSpin2( type[i] ) ); }
 
  return twospin;
}
 
EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type ) {
  int num  = 1;
  _type    = type;
  _twospin = calctwospin( type );
 
  for ( int i = 0; i < _twospin.size(); ++i ) num *= _twospin[i] + 1;
 
  _elem = vector<EvtComplex>( num );
}
 
EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type, const EvtComplex& val ) {
  int num  = 1;
  _type    = type;
  _twospin = calctwospin( type );
 
  for ( int i = 0; i < _twospin.size(); ++i ) num *= _twospin[i] + 1;
 
  _elem = vector<EvtComplex>( num, val );
}
 
EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type,
                        const vector<EvtComplex>&            elem ) {
  int num = 1;
 
  _type    = type;
  _twospin = calctwospin( type );
  _elem    = elem;
 
  for ( int i = 0; i < _twospin.size(); ++i ) num *= _twospin[i] + 1;
 
  if ( _elem.size() != num )
  {
    report( ERROR, "EvtGen" ) << "Wrong number of elements input:" << _elem.size() << " vs. "
                              << num << endl;
    ::abort();
  }
}
 
EvtSpinAmp::EvtSpinAmp( const EvtSpinAmp& copy ) {
  _twospin = copy._twospin;
  _elem    = copy._elem;
  _type    = copy._type;
}
 
void EvtSpinAmp::checktwospin( const vector<int>& twospin ) const {
  if ( _twospin == twospin ) return;
 
  report( ERROR, "EvtGen" ) << "Dimension or order of tensors being operated on does not match"
                            << endl;
  ::abort();
}
 
void EvtSpinAmp::checkindexargs( const vector<int>& index ) const {
  if ( index.size() == 0 )
  {
    report( ERROR, "EvtGen" ) << "EvtSpinAmp can't handle no indices" << endl;
    ::abort();
  }
 
  if ( index.size() != _twospin.size() )
  {
    report( ERROR, "EvtGen" ) << "Rank of EvtSpinAmp index does not match: " << _twospin.size()
                              << " expected " << index.size() << " input." << endl;
    ::abort();
  }
 
  for ( int i = 0; i < _twospin.size(); ++i )
  {
    if ( _twospin[i] >= abs( index[i] ) && _twospin[i] % 2 == index[i] % 2 ) continue;
    report( ERROR, "EvtGen" ) << "EvtSpinAmp index out of range" << endl;
    report( ERROR, "EvtGen" ) << " Index: ";
    for ( int j = 0; j < _twospin.size(); ++j ) report( ERROR, " " ) << _twospin[j];
 
    report( ERROR, " " ) << endl << " Index: ";
    for ( int j = 0; j < index.size(); ++j ) report( ERROR, " " ) << index[j];
    report( ERROR, " " ) << endl;
    ::abort();
  }
}
 
int EvtSpinAmp::findtrueindex( const vector<int>& index ) const {
  int trueindex = 0;
 
  for ( int i = index.size() - 1; i > 0; --i )
  {
    trueindex += ( index[i] + _twospin[i] ) / 2;
    trueindex *= _twospin[i - 1] + 1;
  }
 
  trueindex += ( index[0] + _twospin[0] ) / 2;
 
  return trueindex;
}
 
EvtComplex& EvtSpinAmp::operator()( const vector<int>& index ) {
  checkindexargs( index );
 
  int trueindex = findtrueindex( index );
  if ( trueindex >= _elem.size() )
  {
    report( ERROR, "EvtGen" ) << "indexing error " << trueindex << " " << _elem.size() << endl;
    for ( int i = 0; i < _twospin.size(); ++i ) { report( ERROR, "" ) << _twospin[i] << " "; }
    report( ERROR, "" ) << endl;
 
    for ( int i = 0; i < index.size(); ++i ) { report( ERROR, "" ) << index[i] << " "; }
    report( ERROR, "" ) << endl;
 
    ::abort();
  }
 
  return _elem[trueindex];
}
 
const EvtComplex& EvtSpinAmp::operator()( const vector<int>& index ) const {
  checkindexargs( index );
 
  int trueindex = findtrueindex( index );
  if ( trueindex >= _elem.size() )
  {
    report( ERROR, "EvtGen" ) << "indexing error " << trueindex << " " << _elem.size() << endl;
    for ( int i = 0; i < _twospin.size(); ++i ) { report( ERROR, "" ) << _twospin[i] << " "; }
    report( ERROR, "" ) << endl;
 
    for ( int i = 0; i < index.size(); ++i ) { report( ERROR, "" ) << index[i] << " "; }
    report( ERROR, "" ) << endl;
 
    ::abort();
  }
 
  return _elem[trueindex];
}
 
EvtComplex& EvtSpinAmp::operator()( int i, ... ) {
  va_list     ap;
  vector<int> index( _twospin.size() );
 
  va_start( ap, i );
 
  index[0] = i;
  for ( int n = 1; n < _twospin.size(); ++n ) index[n] = va_arg( ap, int );
 
  va_end( ap );
 
  return ( *this )( index );
}
 
const EvtComplex& EvtSpinAmp::operator()( int i, ... ) const {
  vector<int> index( _twospin.size() );
  va_list     ap;
 
  va_start( ap, i );
 
  index[0] = i;
  for ( int n = 1; n < _twospin.size(); ++n ) index[n] = va_arg( ap, int );
 
  va_end( ap );
 
  return ( *this )( index );
}
 
EvtSpinAmp& EvtSpinAmp::operator=( const EvtSpinAmp& cont ) {
  _twospin = cont._twospin;
  _elem    = cont._elem;
  _type    = cont._type;
 
  return *this;
}
 
EvtSpinAmp EvtSpinAmp::operator+( const EvtSpinAmp& cont ) const {
  checktwospin( cont._twospin );
 
  EvtSpinAmp ret( cont );
  for ( int i = 0; i < ret._elem.size(); ++i ) { ret._elem[i] += _elem[i]; }
 
  return ret;
}
 
EvtSpinAmp& EvtSpinAmp::operator+=( const EvtSpinAmp& cont ) {
  checktwospin( cont._twospin );
 
  for ( int i = 0; i < _elem.size(); ++i ) _elem[i] += cont._elem[i];
 
  return *this;
}
 
EvtSpinAmp EvtSpinAmp::operator-( const EvtSpinAmp& cont ) const {
  checktwospin( cont._twospin );
 
  EvtSpinAmp ret( *this );
  for ( int i = 0; i < ret._elem.size(); ++i ) ret._elem[i] -= cont._elem[i];
 
  return ret;
}
 
EvtSpinAmp& EvtSpinAmp::operator-=( const EvtSpinAmp& cont ) {
  checktwospin( cont._twospin );
 
  for ( int i = 0; i < _elem.size(); ++i ) _elem[i] -= cont._elem[i];
 
  return *this;
}
 
// amp = amp1 * amp2
EvtSpinAmp EvtSpinAmp::operator*( const EvtSpinAmp& amp2 ) const {
  vector<int> index( rank() + amp2.rank() );
  vector<int> index1( rank() ), index2( amp2.rank() );
  EvtSpinAmp  amp;
 
  amp._twospin = _twospin;
  amp._type    = _type;
 
  for ( int i = 0; i < amp2._twospin.size(); ++i )
  {
    amp._twospin.push_back( amp2._twospin[i] );
    amp._type.push_back( amp2._type[i] );
  }
 
  amp._elem = vector<EvtComplex>( _elem.size() * amp2._elem.size() );
 
  for ( int i = 0; i < index1.size(); ++i ) index[i] = index1[i] = -_twospin[i];
 
  for ( int i = 0; i < index2.size(); ++i ) index[i + rank()] = index2[i] = -amp2._twospin[i];
 
  while ( true )
  {
    amp( index ) = ( *this )(index1)*amp2( index2 );
    if ( !amp.iterate( index ) ) break;
 
    for ( int i = 0; i < index1.size(); ++i ) index1[i] = index[i];
 
    for ( int i = 0; i < index2.size(); ++i ) index2[i] = index[i + rank()];
  }
 
  return amp;
}
 
EvtSpinAmp& EvtSpinAmp::operator*=( const EvtSpinAmp& cont ) {
  EvtSpinAmp ret = ( *this ) * cont;
  *this          = ret;
  return *this;
}
 
EvtSpinAmp& EvtSpinAmp::operator*=( const EvtComplex& real ) {
  for ( int i = 0; i < _elem.size(); ++i ) _elem[i] *= real;
 
  return *this;
}
 
EvtSpinAmp& EvtSpinAmp::operator/=( const EvtComplex& real ) {
  for ( int i = 0; i < _elem.size(); ++i ) _elem[i] /= real;
 
  return *this;
}
 
vector<int> EvtSpinAmp::iterinit() const {
  vector<int> init( _twospin.size() );
 
  for ( int i = 0; i < _twospin.size(); ++i ) init[i] = -_twospin[i];
 
  return init;
}
 
bool EvtSpinAmp::iterate( vector<int>& index ) const {
  int last = _twospin.size() - 1;
 
  index[0] += 2;
  for ( int j = 0; j < last; ++j )
  {
    if ( index[j] > _twospin[j] )
    {
      index[j] = -_twospin[j];
      index[j + 1] += 2;
    }
  }
 
  return abs( index[last] ) <= _twospin[last];
}
 
// Test whether a particular index is an allowed one (specifically to deal with
// photons and possibly neutrinos)
bool EvtSpinAmp::allowed( const vector<int>& index ) const {
  if ( index.size() != _type.size() )
  {
    report( ERROR, "EvtGen" ) << "Wrong dimensino index input to allowed." << endl;
    ::abort();
  }
 
  for ( int i = 0; i < index.size(); ++i )
  {
    switch ( _type[i] )
    {
    case EvtSpinType::PHOTON:
      if ( abs( index[i] ) != 2 ) return false;
      break;
    case EvtSpinType::NEUTRINO:
      report( ERROR, "EvtGen" ) << "EvMultibody currently cannot handle neutrinos." << endl;
      ::abort();
    default: break;
    }
  }
 
  return true;
}
 
bool EvtSpinAmp::iterateallowed( vector<int>& index ) const {
  while ( true )
  {
    if ( !iterate( index ) ) return false;
    if ( allowed( index ) ) return true;
  }
}
 
vector<int> EvtSpinAmp::iterallowedinit() const {
  vector<int> init = iterinit();
  while ( !allowed( init ) ) { iterate( init ); }
 
  return init;
}
 
void EvtSpinAmp::intcont( int a, int b ) {
  int newrank = rank() - 2;
  if ( newrank <= 0 )
  {
    report( ERROR, "EvtGen" ) << "EvtSpinAmp can't handle no indices" << endl;
    ::abort();
  }
 
  if ( _twospin[a] != _twospin[b] )
  {
    report( ERROR, "EvtGen" ) << "Contaction called on indices of different dimension" << endl;
    report( ERROR, "EvtGen" ) << "Called on " << _twospin[a] << " and " << _twospin[b] << endl;
    ::abort();
  }
 
  vector<int>                   newtwospin( newrank );
  vector<EvtSpinType::spintype> newtype( newrank );
 
  for ( int i = 0, j = 0; i < _twospin.size(); ++i )
  {
    if ( i == a || i == b ) continue;
 
    newtwospin[j] = _twospin[i];
    newtype[j]    = _type[i];
    ++j;
  }
 
  EvtSpinAmp  newamp( newtype );
  vector<int> index( rank() ), newindex = newamp.iterinit();
 
  for ( int i = 0; i < newrank; ++i ) newindex[i] = -newtwospin[i];
 
  while ( true )
  {
 
    for ( int i = 0, j = 0; i < rank(); ++i )
    {
      if ( i == a || i == b ) continue;
      index[i] = newindex[j];
      ++j;
    }
 
    index[b] = index[a] = -_twospin[a];
    newamp( newindex )  = ( *this )( index );
    for ( int i = -_twospin[a] + 2; i <= _twospin[a]; i += 2 )
    {
      index[b] = index[a] = i;
      newamp( newindex ) += ( *this )( index );
    }
 
    if ( !newamp.iterate( newindex ) ) break;
  }
 
  *this = newamp;
}
 
// In A.extcont(B), a is the index in A and b is the index in B - note that this
// routine can be extremely improved!
void EvtSpinAmp::extcont( const EvtSpinAmp& cont, int a, int b ) {
  EvtSpinAmp ret = ( *this ) * cont;
  ret.intcont( a, rank() + b );
 
  *this = ret;
}