DensityOperator.h

// -*- C++ -*-
//
// DensityOperator.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 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_DensityOperator_H
#define Herwig_DensityOperator_H
//
// This is the declaration of the DensityOperator class.
//
 
#include "ThePEG/Handlers/HandlerBase.h"
 
#include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h"
 
#include <tuple>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
 
namespace Herwig {
 
using namespace ThePEG;
 
typedef boost::numeric::ublas::vector<Complex> CVector;
 
 
 
class DensityOperator: public HandlerBase {
 
public:
 
  DensityOperator();
 
public:
  
  void clear();
 
  void prepare(const cPDVector&);
 
  void fill(const Ptr<MatchboxXComb>::ptr,
            const cPDVector&, const vector<Lorentz5Momentum>& momenta);
  
  void evolve(const map<pair<size_t,size_t>,Complex>& Vijk, 
              const cPDVector& before, 
              const cPDVector& after,
              const map<std::tuple<size_t,size_t,size_t>,map<size_t,size_t> >& emissionsMap,
              const bool splitAGluon,
              const bool initialGluonSplitting);
  
  double colourMatrixElementCorrection(const std::tuple<size_t,size_t,long>& ikemission,
                                       const cPDVector& particles);
 
  void colourConservation(const cPDVector& particles);
 
  Ptr<ColourBasis>::tptr colourBasis() { return theColourBasis; }
  
  const Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }
  
  void colourBasis(Ptr<ColourBasis>::ptr ptr) { theColourBasis = ptr; }
  
  const map<pair<vector<PDT::Colour>,pair<size_t,size_t> >,double>& correlatorMap() const {
    return theCorrelatorMap;
  }
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
 
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
 
private:
 
  double Nc;
 
  double TR;
 
  double colourNorm(const cPDPtr particle);
 
  matrix<Complex> prodSparseDense(const compressed_matrix<double>&,
                                  const matrix<Complex>&);
  matrix<Complex> prodDenseSparse(const matrix<Complex>&,
                                  const compressed_matrix<double>&);
 
  vector<Lorentz5Momentum> boostToRestFrame(const vector<Lorentz5Momentum>& momenta);
 
  bool compareMomentum(const Lorentz5Momentum& p, const Lorentz5Momentum& q);
  
  map<vector<PDT::Colour>,matrix<Complex> > theDensityOperatorMap;
  
  map<pair<vector<PDT::Colour>,pair<size_t,size_t> >,double> theCorrelatorMap;
 
  map<cPDVector, map<size_t,size_t> > theColourBasisToColourBasisMap;
  
  Ptr<ColourBasis>::ptr theColourBasis;
 
  DensityOperator & operator=(const DensityOperator &) = delete;
 
};
 
}
 
#endif /* Herwig_DensityOperator_H */