ColourBasis.h

// -*- C++ -*-
//
// ColourBasis.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_ColourBasis_H
#define HERWIG_ColourBasis_H
//
// This is the declaration of the ColourBasis class.
//
 
#include "ThePEG/Handlers/HandlerBase.h"
 
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/MatrixElement/MEBase.h"
 
#include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
 
#include <iterator>
#include <tuple>
 
namespace Herwig {
 
using std::iterator_traits;
using std::distance;
 
using namespace ThePEG;
 
using boost::numeric::ublas::matrix;
using boost::numeric::ublas::symmetric_matrix;
using boost::numeric::ublas::compressed_matrix;
using boost::numeric::ublas::upper;
 
class ColourBasis: public HandlerBase {
 
public:
 
  ColourBasis();
 
  virtual ~ColourBasis();
 
public:
 
  Ptr<MatchboxFactory>::tptr factory() const;
 
  Ptr<ColourBasis>::ptr cloneMe() const {
    return dynamic_ptr_cast<Ptr<ColourBasis>::ptr>(clone());
  }
 
  virtual void clear();
 
  size_t prepare(const cPDVector&, bool);
 
  size_t prepare(const MEBase::DiagramVector&, bool);
 
  const map<cPDVector,map<size_t,size_t> >& indexMap() const { return theIndexMap; }
 
  virtual map<size_t,vector<vector<size_t> > > basisList(const vector<PDT::Colour>&) const {
    return map<size_t,vector<vector<size_t> > >();
  }
 
  const string& orderingString(const cPDVector& sub, 
                               const map<size_t,size_t>& colourToAmplitude,
                               size_t tensorId);
 
  const set<vector<size_t> >& ordering(const cPDVector& sub, 
                                       const map<size_t,size_t>& colourToAmplitude,
                                       size_t tensorId, size_t shift = 0);
 
  double me2(const cPDVector&, const map<vector<int>,CVector>&) const;
 
  double interference(const cPDVector&, 
                      const map<vector<int>,CVector>&,
                      const map<vector<int>,CVector>&) const;
 
  double colourCorrelatedME2(const pair<size_t,size_t>&,
                             const cPDVector&, 
                             const map<vector<int>,CVector>&) const;
 
  Complex interference(const cPDVector&, 
                       const CVector&, const CVector&) const;
 
  Complex colourCorrelatedInterference(const pair<size_t,size_t>&,
                                       const cPDVector&, 
                                       const CVector&, const CVector&) const;
 
  double me2(const cPDVector&, const matrix<Complex>&) const;
 
  double colourCorrelatedME2(const pair<size_t,size_t>&,
                             const cPDVector&, 
                             const matrix<Complex>&) const;
 
  const symmetric_matrix<double,upper>& scalarProducts(const cPDVector&) const;
 
  const symmetric_matrix<double,upper>& correlator(const cPDVector&,
                                                   const pair<size_t,size_t>&) const;
 
  virtual bool haveColourFlows() const { return false; }
 
  Selector<const ColourLines *> colourGeometries(tcDiagPtr diag,
                                                 const map<vector<int>,CVector>& amps);
 
  size_t tensorIdFromFlow(tcDiagPtr diag, const ColourLines * cl);
 
  struct matchRep {
    PDT::Colour m;
    matchRep(PDT::Colour n)
      : m(n) {}
    bool operator()(PDT::Colour c) const {
      return c == m;
    }
  };
 
  virtual bool largeN() const { return theLargeN; }
 
  void doLargeN(bool yes = true) { theLargeN = yes; }
 
  vector<PDT::Colour> projectColour(const cPDVector&) const;
 
  virtual vector<PDT::Colour> normalOrder(const vector<PDT::Colour>&) const;
 
  vector<PDT::Colour> normalOrderMap(const cPDVector& sub);
 
  const vector<PDT::Colour>& normalOrderedLegs(const cPDVector& sub) const;
 
  const std::tuple<vector<PDT::Colour>,vector<PDT::Colour>,
                   size_t,size_t,size_t,map<size_t,size_t> >& 
  normalOrderEmissionMap(const cPDVector& subFrom,
                         const cPDVector& subTo,
                         size_t ij, size_t i, size_t j,
                         const map<size_t,size_t>& emissionMap);
 
  const pair<compressed_matrix<double>,vector<pair<size_t,size_t> > >&
  charge(const cPDVector& subFrom,
         const cPDVector& subTo,
         size_t ij, size_t i, size_t j,
         const map<size_t,size_t>& emissionMap);
 
  string file(const vector<PDT::Colour>&) const;
 
  void chargeProduct(const compressed_matrix<double>& ti,
                     const vector<pair<size_t,size_t> >& tiNonZero,
                     const symmetric_matrix<double,upper>& X,
                     const compressed_matrix<double>& tj,
                     const vector<pair<size_t,size_t> >& tjNonZero,
                     symmetric_matrix<double,upper>& result) const;
 
  void chargeProductAdd(const compressed_matrix<double>& ti,
                        const vector<pair<size_t,size_t> >& tiNonZero,
                        const matrix<Complex>& X,
                        const compressed_matrix<double>& tj,
                        const vector<pair<size_t,size_t> >& tjNonZero,
                        matrix<Complex>& result,
                        double factor = 1.) const;
 
public:
 
  static list<pair<int,bool> > colouredPath(pair<int,bool> a, pair<int,bool> b,
                                            Ptr<Tree2toNDiagram>::tcptr);
 
  static list<list<list<pair<int,bool> > > > colourFlows(Ptr<Tree2toNDiagram>::tcptr);
 
  static string cfstring(const list<list<pair<int,bool> > >&);
 
  virtual map<size_t,size_t> indexChange(const vector<PDT::Colour>&,
                                         const size_t,
                                         const map<size_t,size_t>&) const {
    map<size_t,size_t> aMap;
    return aMap;
  }
 
 
protected:
 
  virtual size_t prepareBasis(const vector<PDT::Colour>&) = 0;
 
  virtual double scalarProduct(size_t a, size_t b,
                               const vector<PDT::Colour>& abBasis) const = 0;
 
  virtual double tMatrixElement(size_t i, size_t a, size_t b,
                                const vector<PDT::Colour>& aBasis,
                                const vector<PDT::Colour>& bBasis,
                                size_t k, size_t l,
                                const map<size_t,size_t>& dict) const = 0;
 
  virtual bool colourConnected(const cPDVector&,
                               const vector<PDT::Colour>&,
                               const pair<int,bool>&, 
                               const pair<int,bool>&, 
                               size_t) const;
 
  virtual bool colourConnected(const vector<PDT::Colour>&,
                               int, int, size_t) const {
    return false;
  }
 
  vector<string> makeFlows(Ptr<Tree2toNDiagram>::tcptr, size_t) const;
 
  map<Ptr<Tree2toNDiagram>::tcptr,vector<ColourLines*> >&
  colourLineMap();
 
  void updateColourLines(Ptr<Tree2toNDiagram>::tcptr);
 
  
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
 
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
 
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
  virtual void dofinish();
 
private:
 
  typedef map<vector<PDT::Colour>,symmetric_matrix<double,upper> >
  ScalarProductMap;
 
  typedef map<vector<PDT::Colour>,map<pair<size_t,size_t>,symmetric_matrix<double,upper> > >
  CorrelatorMap;
 
  typedef map<std::tuple<vector<PDT::Colour>,vector<PDT::Colour>,
                         size_t,size_t,size_t,map<size_t,size_t> >,
              pair<compressed_matrix<double>,vector<pair<size_t,size_t> > > > TSMap;
 
  bool theLargeN;
 
  map<cPDVector,vector<PDT::Colour> > theNormalOrderedLegs;
 
  map<cPDVector,map<size_t,size_t> > theIndexMap;
 
  map<std::tuple<cPDVector,cPDVector,
                 size_t,size_t,size_t,map<size_t,size_t> >,
      std::tuple<vector<PDT::Colour>,vector<PDT::Colour>,
                 size_t,size_t,size_t,map<size_t,size_t> > >
  theEmissionMaps;
 
  ScalarProductMap theScalarProducts;
 
  CorrelatorMap theCorrelators;
 
  TSMap theCharges;
 
  map<Ptr<Tree2toNDiagram>::tcptr,vector<string> > theFlowMap;
 
  map<Ptr<Tree2toNDiagram>::tcptr,vector<ColourLines*> > theColourLineMap;
 
  map<cPDVector,map<size_t,string> > theOrderingStringIdentifiers;
 
  map<cPDVector,map<size_t,set<vector<size_t> > > > theOrderingIdentifiers;
 
  void writeBasis(const string& prefix = "") const;
 
  void readBasis();
 
  bool readBasis(const vector<PDT::Colour>&);
 
  virtual void readBasisDetails(const vector<PDT::Colour>&) {}
 
  void write(const symmetric_matrix<double,upper>&, ostream&) const;
 
  void read(symmetric_matrix<double,upper>&, istream&);
 
  void write(const compressed_matrix<double>&, ostream&,
             const vector<pair<size_t,size_t> >&) const;
 
  void read(compressed_matrix<double>&, istream&,
            vector<pair<size_t,size_t> >&);
 
  bool didRead;
 
  mutable bool didWrite;
 
  matrix<double> tmp;
 
  string theSearchPath;
 
  ColourBasis & operator=(const ColourBasis &) = delete;
 
};
 
}
 
#endif /* HERWIG_ColourBasis_H */