MatchboxXCombData.h

// -*- C++ -*-
//
// MatchboxXCombData.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_MatchboxXCombData_H
#define Herwig_MatchboxXCombData_H
//
// This is the declaration of the MatchboxXCombData class.
//
 
// work around a Boost 1.64 bug where ublas headers would fail otherwise
#include <boost/version.hpp>
#if (BOOST_VERSION / 100 >= 1064)
#include <boost/serialization/array_wrapper.hpp>
#endif
 
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/vector.hpp>
 
#include "ThePEG/MatrixElement/MEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/Models/StandardModel/StandardModel.h"
 
#include "ThePEG/Persistency/PersistentOStream.fh"
#include "ThePEG/Persistency/PersistentIStream.fh"
 
namespace Herwig {
 
  using namespace ThePEG;
 
  typedef boost::numeric::ublas::vector<Complex> CVector;
 
  typedef map<vector<int>,CVector> AmplitudeMap;
 
  typedef map<vector<int>,CVector>::iterator AmplitudeIterator;
 
  typedef map<vector<int>,CVector>::const_iterator AmplitudeConstIterator;
 
  class MatchboxXCombData {
 
  public:
 
    explicit MatchboxXCombData(tMEPtr newME);
 
    MatchboxXCombData();
 
    virtual ~MatchboxXCombData();
 
 
  public:
 
    void flushCaches();
 
  public:
 
    Ptr<MatchboxFactory>::tcptr factory() const;
 
    Ptr<MatchboxMEBase>::tptr matchboxME() const;
 
    Ptr<SubtractionDipole>::tptr subtractionDipole() const;
 
    const vector<int>& crossingMap() const { return theCrossingMap; }
 
    vector<int>& crossingMap() { return theCrossingMap; }
 
    const map<size_t,size_t>& amplitudeToColourMap() const { return theAmplitudeToColourMap; }
 
    map<size_t,size_t>& amplitudeToColourMap() { return theAmplitudeToColourMap; }
 
    const map<size_t,size_t>& colourToAmplitudeMap() const { return theColourToAmplitudeMap; }
 
    map<size_t,size_t>& colourToAmplitudeMap() { return theColourToAmplitudeMap; }
 
    double crossingSign() const { return theCrossingSign; }
 
    void crossingSign(double c) { theCrossingSign = c; }
 
    Energy2 lastRenormalizationScale() const { return theLastRenormalizationScale; }
 
    void lastRenormalizationScale(Energy2 lrs) { theLastRenormalizationScale = lrs; }
 
    const cPDVector& amplitudePartonData() const { return theAmplitudePartonData; }
 
    cPDVector& amplitudePartonData() { return theAmplitudePartonData; }
 
    const vector<Lorentz5Momentum>& amplitudeMomenta() const { return theAmplitudeMomenta; }
 
    vector<Lorentz5Momentum>& amplitudeMomenta() { return theAmplitudeMomenta; }
 
    bool calculateTreeAmplitudes() const { return theCalculateTreeAmplitudes; }
 
    const map<vector<int>,CVector>& lastAmplitudes() const { return theLastAmplitudes; }
 
    void haveTreeAmplitudes(bool f = true) { theCalculateTreeAmplitudes = !f; }
 
    map<vector<int>,CVector>& lastAmplitudes() { return theLastAmplitudes; }
 
    const map<vector<int>,CVector>& lastLargeNAmplitudes() const { return theLastLargeNAmplitudes; }
 
    map<vector<int>,CVector>& lastLargeNAmplitudes() { return theLastLargeNAmplitudes; }
 
    bool calculateOneLoopAmplitudes() const { return theCalculateOneLoopAmplitudes; }
 
    const map<vector<int>,CVector>& lastOneLoopAmplitudes() const { return theLastOneLoopAmplitudes; }
 
    void haveOneLoopAmplitudes(bool f = true) { theCalculateOneLoopAmplitudes = !f; }
 
    map<vector<int>,CVector>& lastOneLoopAmplitudes() { return theLastOneLoopAmplitudes; }
 
    bool calculateTreeME2() const { return theCalculateTreeME2; }
 
    double lastTreeME2() const { return theLastTreeME2; }
 
    void lastTreeME2(double v) { 
      theLastTreeME2 = v; theCalculateTreeME2 = false;
    }
 
    bool calculateLargeNME2() const { return theCalculateLargeNME2; }
 
    double lastLargeNME2() const { return theLastLargeNME2; }
 
    void lastLargeNME2(double v) { 
      theLastLargeNME2 = v; theCalculateLargeNME2 = false;
    }
 
    bool calculateOneLoopInterference() const { return theCalculateOneLoopInterference; }
 
    double lastOneLoopInterference() const { return theLastOneLoopInterference; }
 
    void lastOneLoopInterference(double v) { 
      theLastOneLoopInterference = v; theCalculateOneLoopInterference = false;
    }
 
    bool calculateOneLoopPoles() const { return theCalculateOneLoopPoles; }
 
    pair<double,double> lastOneLoopPoles() const { return theLastOneLoopPoles; }
 
    void lastOneLoopPoles(pair<double,double> v) { 
      theLastOneLoopPoles = v; theCalculateOneLoopPoles = false;
    }
 
    bool calculateColourCorrelator(pair<int,int> ij) const {
      if ( ij.first > ij.second )
        swap(ij.first,ij.second);
      map<pair<int,int>,bool>::const_iterator f =
        theCalculateColourCorrelators.find(ij);
      if ( f == theCalculateColourCorrelators.end() )
        return true;
      return f->second;
    }
 
    double lastColourCorrelator(pair<int,int> ij) const {
      if ( ij.first > ij.second )
        swap(ij.first,ij.second);
      map<pair<int,int>,double>::const_iterator v =
        theColourCorrelators.find(ij);
      if ( v == theColourCorrelators.end() )
        return 0.;
      return v->second;
    }
 
    void lastColourCorrelator(pair<int,int> ij, double v) {
      if ( ij.first > ij.second )
        swap(ij.first,ij.second);
      theColourCorrelators[ij] = v;
      theCalculateColourCorrelators[ij] = false;
    }
 
    bool calculateLargeNColourCorrelator(pair<int,int> ij) const {
      if ( ij.first > ij.second )
        swap(ij.first,ij.second);
      map<pair<int,int>,bool>::const_iterator f =
        theCalculateLargeNColourCorrelators.find(ij);
      if ( f == theCalculateLargeNColourCorrelators.end() )
        return true;
      return f->second;
    }
 
    double lastLargeNColourCorrelator(pair<int,int> ij) const {
      if ( ij.first > ij.second )
        swap(ij.first,ij.second);
      map<pair<int,int>,double>::const_iterator v =
        theLargeNColourCorrelators.find(ij);
      if ( v == theLargeNColourCorrelators.end() )
        return 0.;
      return v->second;
    }
 
    void lastLargeNColourCorrelator(pair<int,int> ij, double v) {
      if ( ij.first > ij.second )
        swap(ij.first,ij.second);
      theLargeNColourCorrelators[ij] = v;
      theCalculateLargeNColourCorrelators[ij] = false;
    }
 
    bool calculateColourSpinCorrelator(const pair<int,int>& ij) const {
      map<pair<int,int>,bool>::const_iterator f =
        theCalculateColourSpinCorrelators.find(ij);
      if ( f == theCalculateColourSpinCorrelators.end() )
        return true;
      return f->second;
    }
 
    Complex lastColourSpinCorrelator(const pair<int,int>& ij) const {
      map<pair<int,int>,Complex>::const_iterator v =
        theColourSpinCorrelators.find(ij);
      if ( v == theColourSpinCorrelators.end() )
        return 0.;
      return v->second;
    }
 
    void lastColourSpinCorrelator(const pair<int,int>& ij, Complex v) {
      theColourSpinCorrelators[ij] = v;
      theCalculateColourSpinCorrelators[ij] = false;
    }
 
    bool calculateSpinCorrelator(const pair<int,int>& ij) const {
      map<pair<int,int>,bool>::const_iterator f =
        theCalculateSpinCorrelators.find(ij);
      if ( f == theCalculateSpinCorrelators.end() )
        return true;
      return f->second;
    }
 
    Complex lastSpinCorrelator(const pair<int,int>& ij) const {
      map<pair<int,int>,Complex>::const_iterator v =
        theSpinCorrelators.find(ij);
      if ( v == theSpinCorrelators.end() )
        return 0.;
      return v->second;
    }
 
    void lastSpinCorrelator(const pair<int,int>& ij, Complex v) {
      theSpinCorrelators[ij] = v;
      theCalculateSpinCorrelators[ij] = false;
    }
 
    unsigned int nLight() const { return theNLight; }
 
    void nLight(unsigned int n) { theNLight = n; }
 
    vector<long> nLightJetVec() const { return theNLightJetVec; }
 
    void nLightJetVec(long n) { theNLightJetVec.push_back(n); }
 
    vector<long> nHeavyJetVec() const { return theNHeavyJetVec; }
 
    void nHeavyJetVec(long n) { theNHeavyJetVec.push_back(n); }
 
    vector<long> nLightProtonVec() const { return theNLightProtonVec; }
 
    void nLightProtonVec(long n) { theNLightProtonVec.push_back(n); }
 
    size_t colourBasisDim() const { return theColourBasisDim; }
 
    void colourBasisDim(size_t d) { theColourBasisDim = d; }
 
    int nDimPhasespace() const { return theNDimPhasespace; }
 
    void nDimPhasespace(int d) { theNDimPhasespace = d; }
 
    int nDimAmplitude() const { return theNDimAmplitude; }
 
    void nDimAmplitude(int d) { theNDimAmplitude = d; }
 
    int nDimInsertions() const { return theNDimInsertions; }
 
    void nDimInsertions(int d) { theNDimInsertions = d; }
 
    const vector<double>& amplitudeRandomNumbers() const { return theAmplitudeRandomNumbers; }
 
    vector<double>& amplitudeRandomNumbers() { return theAmplitudeRandomNumbers; }
 
    const vector<double>& insertionRandomNumbers() const { return theInsertionRandomNumbers; }
 
    vector<double>& insertionRandomNumbers() { return theInsertionRandomNumbers; }
 
    const map<int,double>& diagramWeights() const { return theDiagramWeights; }
 
    map<int,double>& diagramWeights() { return theDiagramWeights; }
 
    const set<pair<size_t,size_t> >& singularLimits() const { return theSingularLimits; }
 
    set<pair<size_t,size_t> >& singularLimits() { return theSingularLimits; }
 
    const set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() const { return theLastSingularLimit; }
 
    set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() { return theLastSingularLimit; }
 
    void hwStandardModel(Ptr<StandardModel>::tcptr sm) { theStandardModel = sm; }
 
    Ptr<StandardModel>::tcptr hwStandardModel() const { return theStandardModel; }
 
    double symmetryFactor() const { return theSymmetryFactor; }
 
    void symmetryFactor(double f) { theSymmetryFactor = f; }
 
    const vector<int>& olpId() const { return theOLPId; }
 
    void olpId(int pType, int id) {
      if ( theOLPId.empty() )
        theOLPId.resize(4,0);
      theOLPId[pType] = id;
    }
 
    void olpId(const vector<int>& id) { 
      theOLPId = id;
    }
 
    double* olpMomenta() { return theOLPMomenta; }
 
    void fillOLPMomenta(const vector<Lorentz5Momentum>& mm,
                        const cPDVector& mePartonData = cPDVector(),
                        const map<long,Energy>& reshuffleMap = map<long,Energy>());
 
    struct ReshuffleEquation {
 
      ReshuffleEquation(Energy xq,
                        cPDVector::const_iterator xdBegin,
                        cPDVector::const_iterator xdEnd,
                        vector<Lorentz5Momentum>::const_iterator xmBegin,
                        const map<long,Energy>* xreshuffleMap)
        : q(xq), dBegin(xdBegin), dEnd(xdEnd), mBegin(xmBegin),
          reshuffleMap(xreshuffleMap) {}
 
      typedef double ArgType;
      typedef double ValType;
 
      static double aUnit() { return 1.; }
      static double vUnit() { return 1.; }
 
      double operator() (double xi) const;
 
      Energy q;
      cPDVector::const_iterator dBegin;
      cPDVector::const_iterator dEnd;
      vector<Lorentz5Momentum>::const_iterator mBegin;
      const map<long,Energy>* reshuffleMap;
 
    };
 
    void reshuffle(vector<Lorentz5Momentum>& momenta,
                   const cPDVector& mePartonData,
                   const map<long,Energy>& reshuffleMap) const;
 
    int externalId() const { return theExternalId; }
 
    void externalId(int id) { theExternalId = id; }
 
    bool initialized() const { return theInitialized; }
 
    void isInitialized(bool is = true) { theInitialized = is; }
 
    const vector<double*>& externalMomenta() const { return theExternalMomenta; }
 
    void fillExternalMomenta(const vector<Lorentz5Momentum>&);
    
    const map<vector<int>,vector < complex<double> > >& heljamp() const { return theHelJamp; }
    
    const map<vector<int>,vector < complex<double> > >& helLNjamp() const { return theLNHelJamp; }
    
    void pushheljamp(const vector<int>& hel, const complex<double>& jamp) { theHelJamp[hel].push_back(jamp); }
    
    void clearheljamp() { theHelJamp.clear(); }
        
    void pushhelLNjamp(const vector<int>& hel, const complex<double>& jamp) { theLNHelJamp[hel].push_back(jamp); }
    
    void clearhelLNjamp() { theLNHelJamp.clear(); }
 
  public:
 
    void persistentOutput(PersistentOStream & os) const;
 
    void persistentInput(PersistentIStream & is, int version);
 
    static void putCVector(PersistentOStream&, const CVector&);
 
    static void getCVector(PersistentIStream&, CVector&);
 
    static void putAmplitudeMap(PersistentOStream&, const map<vector<int>,CVector>&);
 
    static void getAmplitudeMap(PersistentIStream&, map<vector<int>,CVector>&);
 
    static void Init();
 
  private:
 
    MatchboxXCombData & operator=(const MatchboxXCombData &) = delete;
 
  private:
 
    Ptr<MatchboxMEBase>::tptr theMatchboxME;
 
    Ptr<SubtractionDipole>::tptr theSubtractionDipole;
 
    vector<int> theCrossingMap;
 
    map<size_t,size_t> theAmplitudeToColourMap;
 
    map<size_t,size_t> theColourToAmplitudeMap;
 
    double theCrossingSign;
 
    Energy2 theLastRenormalizationScale;
 
    cPDVector theAmplitudePartonData;
 
    vector<Lorentz5Momentum> theAmplitudeMomenta;
 
    bool theCalculateTreeAmplitudes;
 
    map<vector<int>,CVector> theLastAmplitudes;
 
    map<vector<int>,CVector> theLastLargeNAmplitudes;
 
    bool theCalculateOneLoopAmplitudes;
 
    map<vector<int>,CVector> theLastOneLoopAmplitudes;
 
    bool theCalculateTreeME2;
 
    double theLastTreeME2;
 
    bool theCalculateLargeNME2;
 
    double theLastLargeNME2;
 
    bool theCalculateOneLoopInterference;
 
    double theLastOneLoopInterference;
 
    bool theCalculateOneLoopPoles;
 
    pair<double,double> theLastOneLoopPoles;
 
    map<pair<int,int>,bool> theCalculateColourCorrelators;
 
    map<pair<int,int>,double> theColourCorrelators;
 
    map<pair<int,int>,bool> theCalculateLargeNColourCorrelators;
 
    map<pair<int,int>,double> theLargeNColourCorrelators;
 
    map<pair<int,int>,bool> theCalculateColourSpinCorrelators;
 
    map<pair<int,int>,Complex> theColourSpinCorrelators;
 
    map<pair<int,int>,bool> theCalculateSpinCorrelators;
 
    map<pair<int,int>,Complex> theSpinCorrelators;
 
    static unsigned int theNLight;
 
    static vector<long> theNLightJetVec;
 
    static vector<long> theNHeavyJetVec;
 
    static vector<long> theNLightProtonVec;
 
    size_t theColourBasisDim;
 
    int theNDimPhasespace;
 
    int theNDimAmplitude;
 
    int theNDimInsertions;
 
    vector<double> theAmplitudeRandomNumbers;
 
    vector<double> theInsertionRandomNumbers;
 
    map<int,double> theDiagramWeights;
 
    set<pair<size_t,size_t> > theSingularLimits;
 
    set<pair<size_t,size_t> >::const_iterator theLastSingularLimit;
 
    Ptr<StandardModel>::tcptr theStandardModel;
 
    double theSymmetryFactor;
 
    vector<int> theOLPId;
 
    double* theOLPMomenta;
 
    bool filledOLPMomenta;
 
    int theExternalId;
 
    bool theInitialized;
 
    vector<double*> theExternalMomenta;
 
    bool filledExternalMomenta;
    
    map<vector<int>,vector < complex<double> > > theHelJamp;    
        
    map<vector<int>,vector < complex<double> > > theLNHelJamp;
 
  };
 
}
 
#endif /* Herwig_MatchboxXCombData_H */