MatchboxMEBase.h

// -*- C++ -*-
//
// MatchboxMEBase.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_MatchboxMEBase_H
#define HERWIG_MatchboxMEBase_H
//
// This is the declaration of the MatchboxMEBase class.
//
 
#include "ThePEG/MatrixElement/MEBase.h"
#include "Herwig/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
#include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxReweightBase.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
 
namespace Herwig {
 
using namespace ThePEG;
 
class MatchboxMEBase: 
    public MEBase, public LastMatchboxXCombInfo {
 
public:
 
  MatchboxMEBase();
 
public:
 
  Ptr<MatchboxFactory>::tptr factory() const;
 
 
  const Process& subProcess() const { return theSubprocess; }
 
  Process& subProcess() { return theSubprocess; }
 
  Ptr<Tree2toNGenerator>::tptr diagramGenerator() const;
 
  Ptr<ProcessData>::tptr processData() const;
 
  virtual bool noMirror () const { return true; }
 
  virtual void getDiagrams() const;
  using MEBase::getDiagrams;
 
  virtual Selector<DiagramIndex> diagrams(const DiagramVector &) const;
  using MEBase::diagrams;
 
  virtual Selector<const ColourLines *>
  colourGeometries(tcDiagPtr diag) const;
 
  virtual bool canFillRhoMatrix() const { 
    if ( matchboxAmplitude() )
      return matchboxAmplitude()->canFillRhoMatrix();
    return false;
  }
 
  virtual void constructVertex(tSubProPtr) {}
 
  virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);
 
  virtual unsigned int orderInAlphaS() const;
  using MEBase::orderInAlphaS;
 
  virtual unsigned int orderInAlphaEW() const;
  using MEBase::orderInAlphaEW;  
 
  virtual bool hasInitialAverage() const { 
    return matchboxAmplitude() ? matchboxAmplitude()->hasInitialAverage() : false;
  }
 
  virtual bool hasFinalStateSymmetry() const { 
    return matchboxAmplitude() ? matchboxAmplitude()->hasFinalStateSymmetry() : false; 
  }
 
 
  virtual unsigned int getNLight() const;
 
  virtual vector<long> getNLightJetVec() const;
 
  virtual vector<long> getNHeavyJetVec() const;
 
  virtual vector<long> getNLightProtonVec() const;
 
  virtual bool isOLPTree() const { 
    return matchboxAmplitude() ? matchboxAmplitude()->isOLPTree() : false;
  }
 
  virtual bool isOLPLoop() const { 
    return matchboxAmplitude() ? matchboxAmplitude()->isOLPLoop() : false;
  }
 
  virtual bool needsOLPCorrelators() const { 
    return matchboxAmplitude() ? matchboxAmplitude()->needsOLPCorrelators() : true;
  }
 
  const vector<int>& olpProcess() const { return theOLPProcess; }
 
  void olpProcess(int pType, int id) { 
    if ( theOLPProcess.empty() )
      theOLPProcess.resize(5,0);
    theOLPProcess[pType] = id;
  }
 
  bool noCorrelations() const {
    return theNoCorrelations;
  }
 
  void needsNoCorrelations() {
    theNoCorrelations = true;
  }
 
  void needsCorrelations() {
    theNoCorrelations = false;
  }
 
 
 
  Ptr<MatchboxPhasespace>::tptr phasespace() const { return thePhasespace; }
 
  void phasespace(Ptr<MatchboxPhasespace>::ptr ps) { thePhasespace = ps; }
 
  virtual void setXComb(tStdXCombPtr xc);
 
  virtual bool keepRandomNumbers() const { return true; }
 
  double generateIncomingPartons(const double* r1, const double* r2);
 
  virtual bool generateKinematics(const double * r);
 
  virtual void setKinematics();
 
  virtual void clearKinematics();
 
  virtual int nDim() const;
 
  virtual int nDimBorn() const;
 
  virtual bool haveX1X2() const { 
    return 
      (phasespace() ? phasespace()->haveX1X2() : false) ||
      diagrams().front()->partons().size() == 3;
  }
 
  virtual bool wantCMS() const { 
    return 
      (phasespace() ? phasespace()->wantCMS() : true) &&
      diagrams().front()->partons().size() != 3; }
 
  const vector<Lorentz5Momentum>& lastMEMomenta() const { return meMomenta(); }
 
  vector<Lorentz5Momentum>& lastMEMomenta() { return meMomenta(); }
  
  
  int legsize() const {return int(meMomenta().size());}
 
 
 
  void scaleChoice(Ptr<MatchboxScaleChoice>::ptr sc) { theScaleChoice = sc; }
 
  Ptr<MatchboxScaleChoice>::tptr scaleChoice() const { return theScaleChoice; }
 
  void setScale(Energy2 ren=ZERO,Energy2 fac=ZERO) const;
 
  virtual bool hasRunningAlphaS() const { 
    if ( matchboxAmplitude() )
      return matchboxAmplitude()->hasRunningAlphaS();
    return false;
  }
 
  virtual bool hasRunningAlphaEW() const {
    if ( matchboxAmplitude() )
      return matchboxAmplitude()->hasRunningAlphaEW();
    return false;
  }
 
  virtual Energy2 scale() const { return lastScale(); }
 
  virtual Energy2 factorizationScale() const;
 
  virtual double factorizationScaleFactor() const;
      
      
  virtual double facFac() const{return factorizationScaleFactor();}
 
  virtual Energy2 renormalizationScale() const;
 
  virtual double renormalizationScaleFactor() const;
      
      
  virtual double renFac() const{return renormalizationScaleFactor();}
 
  virtual Energy2 renormalizationScaleQED() const;
 
  virtual Energy2 showerScale() const;
 
  virtual void setVetoScales(tSubProPtr) const;
 
  bool fixedCouplings() const;
 
  bool fixedQEDCouplings() const;
 
  virtual double alphaS() const { return lastAlphaS(); }
 
  virtual double alphaEM() const { return lastAlphaEM(); }
 
  virtual bool havePDFWeight1() const;
 
  virtual bool havePDFWeight2() const;
 
  void getPDFWeight(Energy2 factorizationScale = ZERO) const;
 
  double pdf1(Energy2 factorizationScale = ZERO,
              double xEx = 1., double xFactor = 1.) const;
 
  double pdf2(Energy2 factorizationScale = ZERO,
              double xEx = 1., double xFactor = 1.) const;
 
 
 
  Ptr<MatchboxAmplitude>::tptr matchboxAmplitude() const { return theAmplitude; }
 
  void matchboxAmplitude(Ptr<MatchboxAmplitude>::ptr amp) { theAmplitude = amp; }
 
  virtual double me2() const;
 
  virtual double largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const;
 
  virtual double finalStateSymmetry() const;
 
  double me2Norm(unsigned int addAlphaS = 0) const;
 
  virtual CrossSection dSigHatDR() const;
 
  CrossSection prefactor()const;
 
  CrossSection dSigHatDRB() const ;
      
  CrossSection dSigHatDRV() const ;
      
  CrossSection dSigHatDRI() const ;
      
  CrossSection dSigHatDRAlphaDiff(double alpha) const ;
      
 
 
 
  virtual double oneLoopInterference() const;
 
  virtual bool haveOneLoop() const;
 
  virtual bool onlyOneLoop() const;
 
  virtual bool isDRbar() const;
 
  virtual bool isDR() const;
 
  virtual bool isCS() const;
 
  virtual bool isBDK() const;
 
  virtual bool isExpanded() const;
 
  virtual Energy2 mu2() const;
 
  virtual double oneLoopDoublePole() const;
 
  virtual double oneLoopSinglePole() const;
 
  bool checkPoles() const;
 
  struct AccuracyHistogram {
 
    double lower;
 
    double upper;
 
    map<double,double> bins;
 
    unsigned long sameSign;
 
    unsigned long oppositeSign;
 
    unsigned long nans;
 
    unsigned long overflow;
 
    unsigned long underflow;
 
    AccuracyHistogram(double low = -40.,
                      double up = 0.,
                      unsigned int nbins = 80);
 
    void book(double a, double b);
 
    void dump(const std::string& folder, const std::string& prefix,
              const cPDVector& proc) const;
 
    void persistentOutput(PersistentOStream&) const;
 
    void persistentInput(PersistentIStream&);
 
  };
 
  void logPoles() const;
 
  const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const {
    return theVirtuals;
  }
 
  vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() {
    return theVirtuals;
  }
 
  void doOneLoop() { theOneLoop = true; }
      
  void noOneLoop() { theOneLoop = false; }
 
  bool oneLoop() const { return theOneLoop; }
 
  void doOneLoopNoBorn() { theOneLoop = true; theOneLoopNoBorn = true; }
      
  void noOneLoopNoBorn() { theOneLoop = false; theOneLoopNoBorn = false; }
 
  bool oneLoopNoBorn() const { return theOneLoopNoBorn || onlyOneLoop(); }
 
  void doOneLoopNoLoops() { theOneLoop = true; theOneLoopNoLoops = true; }
 
  bool oneLoopNoLoops() const { return theOneLoopNoLoops; }
 
 
 
  vector<Ptr<SubtractionDipole>::ptr> 
  getDipoles(const vector<Ptr<SubtractionDipole>::ptr>&,
             const vector<Ptr<MatchboxMEBase>::ptr>&,bool slim=false) const;
 
  
  virtual bool noDipole(int,int,int) const { return false; }
 
  virtual bool noDipole(int,int) const { return false; }
 
  virtual double colourCorrelatedME2(pair<int,int>) const;
 
  virtual double largeNColourCorrelatedME2(pair<int,int> ij,
                                           Ptr<ColourBasis>::tptr largeNBasis) const;
 
  virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
                                         const SpinCorrelationTensor& c) const;
 
  virtual double spinCorrelatedME2(pair<int,int> emitterSpectator,
                                   const SpinCorrelationTensor& c) const;
 
 
 
  virtual void flushCaches();
 
  bool verbose() const;
 
  bool initVerbose() const;
 
  void print(ostream&) const;
 
  virtual void printLastEvent(ostream&) const;
 
  void logGenerateKinematics(const double * r) const;
 
  void logSetScale() const;
 
  void logPDFWeight() const;
 
  void logME2() const;
 
  void logDSigHatDR() const;
 
 
 
  void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }
 
  const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }
 
  vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }
 
  const MergerBasePtr merger() const;
    
  MergerBasePtr merger() ;
      
  void merger(MergerBasePtr v);
  
 
 
  virtual bool preInitialize() const { return true; }
 
  Ptr<MatchboxMEBase>::ptr cloneMe() const {
    return dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(clone());
  }
 
  void cloneDependencies(const std::string& prefix = "",bool slim = false );
 
  void prepareXComb(MatchboxXCombData&) const;
 
  virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
                                tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
                                tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
                                const PBPair & newPartonBins, tCutsPtr newCuts,
                                const DiagramVector & newDiagrams, bool mir,
                                const PartonPairVec& allPBins,
                                tStdXCombPtr newHead = tStdXCombPtr(),
                                tMEPtr newME = tMEPtr());
 
  virtual StdXCombPtr makeXComb(tStdXCombPtr newHead,
                                const PBPair & newPartonBins,
                                const DiagramVector & newDiagrams,
                                tMEPtr newME = tMEPtr());
 
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
protected:
 
  virtual void doinit();
 
  virtual void doinitrun();
 
  virtual void dofinish();
 
private:
 
  Ptr<MatchboxPhasespace>::ptr thePhasespace;
 
  Ptr<MatchboxAmplitude>::ptr theAmplitude;
 
  Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
 
  vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
 
  vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
 
private:
 
  Process theSubprocess;
 
  bool theOneLoop;
 
  bool theOneLoopNoBorn;
 
  bool theOneLoopNoLoops;
 
  vector<int> theOLPProcess;
 
  mutable map<cPDVector,AccuracyHistogram> epsilonSquarePoleHistograms;
 
  mutable map<cPDVector,AccuracyHistogram> epsilonPoleHistograms;
 
  bool theNoCorrelations;
 
  mutable pair<bool,bool> theHavePDFs;
 
  mutable bool checkedPDFs;
  
  MergerBasePtr theMerger;
 
 
private:
 
  MatchboxMEBase & operator=(const MatchboxMEBase &) = delete;
 
};
 
inline PersistentOStream& operator<<(PersistentOStream& os,
                                     const MatchboxMEBase::AccuracyHistogram& h) {
  h.persistentOutput(os);
  return os;
}
 
inline PersistentIStream& operator>>(PersistentIStream& is,
                                     MatchboxMEBase::AccuracyHistogram& h) {
  h.persistentInput(is);
  return is;
}
 
}
 
#endif /* HERWIG_MatchboxMEBase_H */