ThePEG/EventRecord/RhoDMatrix.h

0001 // -*- C++ -*-
0002 //
0003 // RhoDMatrix.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 2003-2019 Peter Richardson, Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_RhoDMatrix_H
0010 #define ThePEG_RhoDMatrix_H
0011 // This is the declaration of the RhoDMatrix class.
0012
0013 #include "ThePEG/PDT/PDT.h"
0014 #include "ThePEG/Helicity/HelicityDefinitions.h"
0015 #include <cassert>
0016 #include <array>
0017
0018 namespace ThePEG {
0019
0020 /**
0021  *  The RhoDMatrix class is designed to implement the storage of the
0022  *  rho and D matrices which are required for the spin correlation
0023  *  algorithm.  The matrix stores the spin as 2s+1.
0024  *
0025  * @author Peter Richardson
0026  *
0027  */
0028 class RhoDMatrix {
0029
0030 public:
0031
0032   /** @name Standard constructors and destructors. */
0033   //@{
0034   /**
0035    * Default constructor with undefined spin.
0036    */
0037   RhoDMatrix() = default;
0038
0039   /**
0040    * Standard constructor giving the spin as 2s+1. The matrix starts out averaged,
0041    * unless the second argument is false, when it is zeroed.
0042    */
0043   RhoDMatrix(PDT::Spin inspin, bool average = true)
0044   : _spin(inspin), _ispin(abs(int(inspin))) {
0045     assert(_ispin <= MAXSPIN);
0046     // initialize to average
0047     if ( average )
0048         for(size_t ix=0; ix<_ispin; ++ix)
0049         _matrix[ix][ix] = 1./_ispin;
0050   }
0051   //@}
0052
0053 public:
0054
0055   /** @name Access matrix elements. */
0056   //@{
0057   /**
0058    * Return an element of the matrix.
0059    */
0060   Complex operator() (size_t ix, size_t iy) const {
0061     assert(ix < _ispin);
0062     assert(iy < _ispin);
0063     return _matrix[ix][iy];
0064   }
0065
0066   /**
0067    * Set an element of the matrix.
0068    */
0069   Complex & operator() (size_t ix, size_t iy) {
0070     assert(ix < _ispin);
0071     assert(iy < _ispin);
0072     return _matrix[ix][iy];
0073   }
0074
0075   /**
0076    * renormalise the matrix so it has unit trace
0077    */
0078   void normalize() {
0079 #ifndef NDEBUG
0080     static const double epsa=1e-40, epsb=1e-10;
0081 #endif
0082     Complex norm = 0.;
0083     for(size_t ix=0; ix<_ispin; ++ix)
0084       norm += _matrix[ix][ix];
0085     assert(norm.real() > epsa);
0086     assert(norm.imag()/norm.real() < epsb);
0087     double invnorm = 1./norm.real();
0088     for(size_t ix=0; ix<_ispin; ++ix)
0089       for(size_t iy=0; iy<_ispin; ++iy)
0090     _matrix[ix][iy]*=invnorm;
0091   }
0092   //@}
0093
0094   /**
0095    *  Reset
0096    */
0097   void reset(bool average = true) {
0098     for(size_t ix=0; ix<_ispin; ++ix)
0099       for(size_t iy=0; iy<_ispin; ++iy)
0100     _matrix[ix][iy]=0.;
0101     if ( average )
0102       for(size_t ix=0; ix<_ispin; ++ix)
0103     _matrix[ix][ix] = 1./_ispin;
0104   }
0105
0106   /** @name Access the spin. */
0107   //@{
0108
0109   /**
0110    * Get the spin. The spin is returned as 2J+1 in units of hbar/2.
0111    */
0112   PDT::Spin iSpin() const { return _spin; }
0113   //@}
0114
0115   /**
0116    * Output the spin density matrix for debugging purposes.
0117    */
0118   friend ostream & operator<<(ostream & os, const RhoDMatrix & rd);
0119
0120 private:
0121
0122   /**
0123    * 2s+1 for the particle.
0124    */
0125   PDT::Spin _spin;
0126
0127   /**
0128    *  Storage of 2s+1 for speed.
0129    */
0130   size_t _ispin;
0131
0132   /**
0133    * Spin matrix size
0134    */
0135   enum { MAXSPIN = 7 };
0136
0137   /**
0138    * Storage for the matrix allowing up to spin 2 particles.
0139    */
0140   // Deliberately not using vector<> to avoid calls to 'new'
0141   // from this commonly used class.
0142   std::array<std::array<Complex,MAXSPIN>,MAXSPIN> _matrix;
0143
0144 };
0145
0146 /** Output operator */
0147 inline ostream & operator<<(ostream & os, const RhoDMatrix & rd) {
0148   for (size_t ix = 0; ix < rd._ispin; ++ix) {
0149     for (size_t iy = 0; iy < rd._ispin; ++iy)
0150       os << rd._matrix[ix][iy] << "  ";
0151     os << '\n';
0152   }
0153   return os;
0154 }
0155
0156 }
0157
0158 #endif /* ThePEG_RhoDMatrix_H */