RhoDMatrix.h

// -*- C++ -*-
//
// RhoDMatrix.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2007 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 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>

namespace ThePEG {
namespace Helicity {

class RhoDMatrix {

public:

  RhoDMatrix() : _spin(), _ispin() {}

  RhoDMatrix(PDT::Spin inspin, bool average = true) 
  : _spin(inspin), _ispin(abs(int(inspin))) {
    assert(_ispin <= MAXSPIN);
    // initialize to average
    for(size_t ix=0; ix<_ispin; ++ix)
      for(size_t iy=0; iy<_ispin; ++iy) 
        _matrix[ix][iy] = (average && ix==iy) ? 1./_ispin : 0.;
  }

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-30, 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;
  }


  PDT::Spin iSpin() const { return _spin; }

  friend ostream & operator<<(ostream & os, const RhoDMatrix & rd);

private:

  PDT::Spin _spin;

  size_t _ispin;

  enum { MAXSPIN = 5 };

  // Deliberately not using vector<> to avoid calls to 'new' 
  // from this commonly used class.
  Complex _matrix[MAXSPIN][MAXSPIN];

};

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 */

---