RhoDMatrix.h

// -*- C++ -*-
//
// RhoDMatrix.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2019 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_RhoDMatrix_H
#define ThePEG_RhoDMatrix_H
// This is the declaration of the RhoDMatrix class.
 
#include "ThePEG/PDT/PDT.h"
#include "ThePEG/Helicity/HelicityDefinitions.h"
#include <cassert>
#include <array>
 
namespace ThePEG {
 
class RhoDMatrix {
 
public:
 
  RhoDMatrix() = default;
 
  RhoDMatrix(PDT::Spin inspin, bool average = true) 
  : _spin(inspin), _ispin(abs(int(inspin))) {
    assert(_ispin <= MAXSPIN);
    // initialize to average
    if ( average )
        for(size_t ix=0; ix<_ispin; ++ix)
            _matrix[ix][ix] = 1./_ispin;
  }
 
public:
 
  Complex operator() (size_t ix, size_t iy) const {
    assert(ix < _ispin);
    assert(iy < _ispin);
    return _matrix[ix][iy];
  }
 
  Complex & operator() (size_t ix, size_t iy) {
    assert(ix < _ispin);
    assert(iy < _ispin);
    return _matrix[ix][iy];
  }
 
  void normalize() {
#ifndef NDEBUG
    static const double epsa=1e-40, epsb=1e-10;
#endif
    Complex norm = 0.;
    for(size_t ix=0; ix<_ispin; ++ix) 
      norm += _matrix[ix][ix];
    assert(norm.real() > epsa);
    assert(norm.imag()/norm.real() < epsb);
    double invnorm = 1./norm.real();
    for(size_t ix=0; ix<_ispin; ++ix)
      for(size_t iy=0; iy<_ispin; ++iy) 
        _matrix[ix][iy]*=invnorm;
  }
  
  void reset(bool average = true) {
    for(size_t ix=0; ix<_ispin; ++ix)
      for(size_t iy=0; iy<_ispin; ++iy)
        _matrix[ix][iy]=0.;
    if ( average )
      for(size_t ix=0; ix<_ispin; ++ix)
        _matrix[ix][ix] = 1./_ispin;
  }
    
 
  PDT::Spin iSpin() const { return _spin; }
 
  friend ostream & operator<<(ostream & os, const RhoDMatrix & rd);
 
private:
 
  PDT::Spin _spin;
 
  size_t _ispin;
 
  enum { MAXSPIN = 7 };
 
  // Deliberately not using vector<> to avoid calls to 'new' 
  // from this commonly used class.
  std::array<std::array<Complex,MAXSPIN>,MAXSPIN> _matrix;
 
};
 
inline ostream & operator<<(ostream & os, const RhoDMatrix & rd) {
  for (size_t ix = 0; ix < rd._ispin; ++ix) {
    for (size_t iy = 0; iy < rd._ispin; ++iy)
      os << rd._matrix[ix][iy] << "  ";
    os << '\n';
  }
  return os;
}
 
}
 
#endif /* ThePEG_RhoDMatrix_H */