MixingMatrix.h

// -*- C++ -*-
//
// MixingMatrix.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MixingMatrix_H
#define HERWIG_MixingMatrix_H
 
//
// This is the declaration of the MixingMatrix class.
//
 
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/Config/Complex.h"
#include "MixingMatrix.fh"
#include "ThePEG/Utilities/Exception.h"
 
namespace Herwig {
using namespace ThePEG;
 
/*@name Some convenient typedefs. */
 
  typedef vector<vector<Complex> > CMatrix;
  
  struct MixingElement{
 
    MixingElement(unsigned int irow, unsigned int icol,Complex ivalue) 
      : row(irow), col(icol), value(ivalue) {}
 
    unsigned int row;
 
    unsigned int col;
 
    Complex value;
  };
 
  typedef vector<MixingElement> MixingVector;
 
  typedef pair<unsigned int, unsigned int> MatrixSize;
 
class MixingMatrix: public Interfaced {
 
  class MixingMatrixError : public Exception {};
 
public:
 
  MixingMatrix(const CMatrix & mix,const vector<long> & ids) : 
    mixingMatrix_(mix),ids_(ids), size_(make_pair(mix.size(),mix[0].size())) 
  {}
 
  MixingMatrix(unsigned int row, unsigned int col) :
    mixingMatrix_(row,vector<Complex>(col,Complex(0.,0.))), size_(row,col)
  {}
 
  MixingMatrix() {}
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
  void setMatrix(const CMatrix & mixing)  {
    mixingMatrix_ = mixing;
    size_ = make_pair(mixing.size(),mixing[0].size());
  }
 
  CMatrix getMatrix() const {return mixingMatrix_;}
 
  void setIds(const vector<long> & mixingCodes)  {
    if(mixingCodes.size() != size_.first) {
      ostringstream codes;
      for(unsigned int ix=0;ix<mixingCodes.size();++ix)
        codes << mixingCodes[ix] << " ";
      throw MixingMatrixError() << "MixingMatrix::setIds() - The number "
                                << "of PDG codes (" << mixingCodes.size()
                                << ") does not match the size of the "
                                << "matrix (" << size_.first
                                << ")"
                                << "Ids are " << codes.str() 
                                << Exception::warning;
      return;
    }
    ids_ = mixingCodes;
  }
  
  const vector<long> & getIds() const {return ids_;}
  
  void adjustPhase(long id);
    
  const Complex operator()(unsigned int row, unsigned int col) const {
    return mixingMatrix_.at(row).at(col);
  }
 
  Complex & operator()(unsigned int row, unsigned int col) {
    return mixingMatrix_.at(row).at(col);
  }
  
  void addCode(long id) {
    if(ids_.size() >= size_.first) {
      throw MixingMatrixError() << "MixingMatrix::addCode() - Trying to add a"
                                << "PDG code but the vector already contains the "
                                << "same number as the matrix size " 
                                << Exception::warning;
      return;
    }
    ids_.push_back(id);
  }
  
  MatrixSize size() const {return size_;}
  
protected:
 
  virtual IBPtr clone() const {return new_ptr(*this);}
 
  virtual IBPtr fullclone() const {return new_ptr(*this);}
 
private:
 
  MixingMatrix & operator=(const MixingMatrix &) = delete;
 
  CMatrix mixingMatrix_;
 
  vector<long> ids_;
 
  pair<unsigned int,unsigned int> size_;
  
  friend ostream & operator<<(ostream & os,const MixingMatrix & mix);
};
 
  ostream & operator<<(ostream &,const MixingMatrix &);
}
 
#endif /* HERWIG_MixingMatrix_H */