SubtractionDipole.h

// -*- C++ -*-
//
// SubtractionDipole.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_SubtractionDipole_H
#define HERWIG_SubtractionDipole_H
//
// This is the declaration of the SubtractionDipole class.
//
 
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.fh"
#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.fh"
 
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
 
namespace Herwig {
 
using namespace ThePEG;
 
class SubtractionDipole: 
    public MEBase,  public LastMatchboxXCombInfo {
 
public:
 
  SubtractionDipole();
 
  virtual ~SubtractionDipole();
 
public:
 
  Ptr<MatchboxFactory>::tptr factory() const;
 
 
  struct MergeInfo {
 
    int emitter;
 
    Ptr<Tree2toNDiagram>::ptr diagram;
 
    map<int,int> mergeLegs;
 
  };
 
  virtual bool canHandleEmitter(const cPDVector& partons, int emitter) const = 0;
 
  virtual bool canHandleSplitting(const cPDVector& partons, int emitter, int emission) const = 0;
 
  virtual bool canHandleSpectator(const cPDVector& partons, int spectator) const = 0;
 
  virtual bool canHandle(const cPDVector& partons,
                         int emitter, int emission, int spectator) const = 0;
 
  virtual bool isSymmetric() const { return false; }
 
  virtual bool apply() const { return theApply; }
 
  void clearBookkeeping();
 
  void setupBookkeeping(const map<Ptr<DiagramBase>::ptr,MergeInfo>& mergeInfo,bool slim);
 
  void subtractionBookkeeping();
 
  void splittingBookkeeping();
 
  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());
 
  StdXCombPtr makeBornXComb(tStdXCombPtr realXC);
 
  vector<StdXCombPtr> makeRealXCombs(tStdXCombPtr bornXC);
 
  bool empty() const { return theSplittingMap.empty()&&theMergingMap.empty(); }
 
  int realEmitter() const { return theRealEmitter; }
 
  void realEmitter(int id) { theRealEmitter = id; }
 
  int realEmission() const { return theRealEmission; }
 
  void realEmission(int id) { theRealEmission = id; }
 
  int realSpectator() const { return theRealSpectator; }
 
  void realSpectator(int id) { theRealSpectator = id; }
 
  int bornEmitter() const { return theBornEmitter; }
 
  void bornEmitter(int id) { theBornEmitter = id; }
 
  int bornSpectator() const { return theBornSpectator; }
 
  void bornSpectator(int id) { theBornSpectator = id; }
 
  typedef pair<pair<cPDVector,int>,pair<int,int> > RealEmissionKey;
 
  static RealEmissionKey realEmissionKey(const cPDVector& proc, 
                                         int em, int emm, int sp) {
    return make_pair(make_pair(proc,emm),make_pair(em,sp));
  }
 
  static const cPDVector& process(const RealEmissionKey& key) {
    return key.first.first;
  }
 
  static int emission(const RealEmissionKey& key) {
    return key.first.second;
  }
 
  static int emitter(const RealEmissionKey& key) {
    return key.second.first;
  }
 
  static int spectator(const RealEmissionKey& key) {
    return key.second.second;
  }
 
  typedef pair<cPDVector,pair<int,int> > UnderlyingBornKey;
 
  static UnderlyingBornKey underlyingBornKey(const cPDVector& proc, 
                                             int em, int sp) {
    return make_pair(proc,make_pair(em,sp));
  }
 
  static const cPDVector& process(const UnderlyingBornKey& key) {
    return key.first;
  }
 
  static int emitter(const UnderlyingBornKey& key) {
    return key.second.first;
  }
 
  static int spectator(const UnderlyingBornKey& key) {
    return key.second.second;
  }
 
  typedef pair<RealEmissionKey,map<int,int> > RealEmissionInfo;
 
  typedef pair<UnderlyingBornKey,map<int,int> > UnderlyingBornInfo;
 
  const map<RealEmissionKey,UnderlyingBornInfo>& mergingMap() const { return theMergingMap; }
 
  const multimap<UnderlyingBornKey,RealEmissionInfo>& splittingMap() const { return theSplittingMap; }
 
  const DiagramVector& underlyingBornDiagrams(const cPDVector& real) const;
 
  tcDiagPtr underlyingBornDiagram(tcDiagPtr realDiag) const;
 
  const DiagramVector& realEmissionDiagrams(const cPDVector& born) const;
 
  tcDiagPtr realEmissionDiagram(tcDiagPtr bornDiag) const;
 
  virtual void getDiagrams() const;
 
  virtual bool noMirror () const { return true; }
 
  virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
 
  virtual Selector<const ColourLines *>
  colourGeometries(tcDiagPtr diag) const;
 
  virtual const ColourLines &
  selectColourGeometry(tcDiagPtr diag) const;
 
  virtual unsigned int orderInAlphaS() const { return realEmissionME()->orderInAlphaS(); }
 
  virtual unsigned int orderInAlphaEW() const { return underlyingBornME()->orderInAlphaEW(); }
 
 
 
  virtual void setXComb(tStdXCombPtr xc);
 
  virtual void setKinematics();
 
  virtual bool generateKinematics(const double * r);
 
  virtual int nDim() const;
 
  virtual bool wantCMS () const { return realEmissionME()->wantCMS(); }
 
  virtual void clearKinematics();
 
  virtual bool ignoreCuts() const { return theIgnoreCuts; }
 
  void doIgnoreCuts(bool is = true) { theIgnoreCuts = is; }
 
 
 
  Ptr<TildeKinematics>::tcptr tildeKinematics() const { return theTildeKinematics; }
 
  void tildeKinematics(Ptr<TildeKinematics>::tptr);
 
  virtual bool generateTildeKinematics();
 
  Ptr<InvertedTildeKinematics>::tcptr invertedTildeKinematics() const { return theInvertedTildeKinematics; }
 
  void invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr);
 
  virtual int nDimRadiation() const;
 
  virtual bool generateRadiationKinematics(const double *);
 
  void ptCut(Energy cut);
 
  Energy lastDipoleScale() const {
    return splitting() ? theLastSplittingScale : theLastSubtractionScale;
  }
 
  Energy lastPt() const {
    return splitting() ? theLastSplittingPt : theLastSubtractionPt;
  }
 
  double lastZ() const {
    return splitting() ? theLastSplittingZ : theLastSubtractionZ;
  }
 
  bool splitting() const { return theSplitting; }
 
  void doSplitting() { theSplitting = true; }
 
  void doSubtraction() { theSplitting = false; }
 
  const vector<double>& subtractionParameters() const { return theSubtractionParameters; }
 
  vector<double>& subtractionParameters() { return theSubtractionParameters; }
 
  Energy showerHardScale() const { return theShowerHardScale; }
 
  void showerHardScale(Energy s) { theShowerHardScale = s; }
 
  Energy showerScale() const { return theShowerScale; }
 
  void showerScale(Energy s) { theShowerScale = s; }
 
  const vector<double>& showerParameters() const { return theShowerParameters; }
 
  vector<double>& showerParameters() { return theShowerParameters; }
 
  bool isInShowerPhasespace() const { return theIsInShowerPhasespace; }
 
  void isInShowerPhasespace(bool yes) { theIsInShowerPhasespace = yes; }
 
  bool isAboveCutoff() const { return theIsAboveCutoff; }
 
  void isAboveCutoff(bool yes) { theIsAboveCutoff = yes; }
 
 
 
  bool realEmissionScales() const { return theRealEmissionScales; }
 
  void doRealEmissionScales(bool on = true) { theRealEmissionScales = on; }
 
  virtual Energy2 scale() const { 
    return realEmissionScales() ? 
      realEmissionME()->scale() :
      underlyingBornME()->scale();
  }
 
  virtual double alphaS() const { 
    return realEmissionScales() ? 
      realEmissionME()->alphaS() :
      underlyingBornME()->alphaS();
  }
 
  virtual double alphaEM() const { 
    return realEmissionScales() ? 
      realEmissionME()->alphaEM() :
      underlyingBornME()->alphaEM();
  }
 
  virtual bool havePDFWeight1() const { return realEmissionME()->havePDFWeight1(); }
 
  virtual bool havePDFWeight2() const { return realEmissionME()->havePDFWeight2(); }
 
  virtual int samplingZ() const {return 4;}
 
 
  Ptr<MatchboxMEBase>::tcptr realEmissionME() const { 
    return theRealEmissionME;
  }
 
  Ptr<MatchboxMEBase>::tptr realEmissionME() { 
    return theRealEmissionME;
  }
 
  void realEmissionME(Ptr<MatchboxMEBase>::tptr me) { theRealEmissionME = me; }
 
  Ptr<MatchboxMEBase>::tcptr underlyingBornME() const { 
    return theUnderlyingBornME;
  }
 
  Ptr<MatchboxMEBase>::tptr underlyingBornME() { 
    return theUnderlyingBornME;
  }
 
  void underlyingBornME(Ptr<MatchboxMEBase>::tptr me) { theUnderlyingBornME = me; }
 
  void partnerDipoles(const vector<Ptr<SubtractionDipole>::tptr>& p) {
    thePartners = p;
  }
 
  const vector<Ptr<SubtractionDipole>::tptr>& partnerDipoles() const {
    return thePartners;
  }
 
  virtual double me2Avg(double ccme2) const = 0;
 
  virtual CrossSection dSigHatDR(Energy2 factorizationScale) const;
 
  virtual CrossSection dSigHatDR() const { return dSigHatDR(ZERO); }
 
      
  CrossSection prefactor(Energy2 factorizationScale)const;
      
  CrossSection ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;
 
  CrossSection dip(Energy2 factorizationScale) const;
 
  
 
  pair<CrossSection,CrossSection> dipandPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;
 
 
 
  void showerApproximation(Ptr<ShowerApproximation>::tptr app) {
    theShowerApproximation = app;
  }
 
  Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
 
  void doRealShowerSubtraction() { theRealShowerSubtraction = true; }
 
  bool realShowerSubtraction() const { return theRealShowerSubtraction; }
 
  void doVirtualShowerSubtraction() { theVirtualShowerSubtraction = true; }
 
  bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
 
  void doLoopSimSubtraction() { theLoopSimSubtraction = true; }
 
  bool loopSimSubtraction() const { return theLoopSimSubtraction; }
 
 
 
  virtual void flushCaches();
 
  void doTestSubtraction() { theSubtractionTest = true; }
 
  bool testSubtraction() const { return theSubtractionTest; }
 
  bool verbose() const { return realEmissionME()->verbose() || underlyingBornME()->verbose(); }
 
  void print(ostream&) const;
 
  virtual void printLastEvent(ostream&) const;
 
  void logGenerateTildeKinematics() const;
 
  void logGenerateRadiationKinematics(const double * r) const;
 
  void logME2() const;
 
  void logDSigHatDR(double effectiveJac) 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; }
 
 
 
  Ptr<SubtractionDipole>::ptr cloneMe() const {
    return dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(clone());
  }
 
  void cloneDependencies(const std::string& prefix = "", bool slim=false);
 
 
 
  virtual void constructVertex(tSubProPtr sub);
 
  virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);
 
  virtual void generateSubCollision(SubProcess & sub);
      
   double alpha() const;
      
   /*
    * True if phase space point is above the alpha cut for this dipole.
    */
      
   bool aboveAlpha() const;
 
 
public:
 
  void persistentOutput(PersistentOStream & os) const;
 
  void persistentInput(PersistentIStream & is, int version);
 
  static void Init();
 
protected:
 
 
  virtual void doinit();
 
  virtual void doinitrun();
 
private:
 
  bool theSplitting;
 
  bool theApply;
 
  bool theSubtractionTest;
 
  bool theIgnoreCuts;
 
  Ptr<MatchboxMEBase>::ptr theRealEmissionME;
 
  Ptr<MatchboxMEBase>::ptr theUnderlyingBornME;
 
  vector<Ptr<SubtractionDipole>::tptr> thePartners;
 
  Ptr<TildeKinematics>::ptr theTildeKinematics;
 
  Ptr<InvertedTildeKinematics>::ptr theInvertedTildeKinematics;
 
  vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
 
  int theRealEmitter;
 
  int theRealEmission;
 
  int theRealSpectator;
 
  vector<double> theSubtractionParameters;
 
  map<RealEmissionKey,UnderlyingBornInfo> theMergingMap;
 
  multimap<UnderlyingBornKey,RealEmissionInfo> theSplittingMap;
 
  map<cPDVector,pair<int,int> > theIndexMap;
 
  map<cPDVector,DiagramVector> theUnderlyingBornDiagrams;
 
  map<cPDVector,DiagramVector> theRealEmissionDiagrams;
 
  map<tcDiagPtr,tcDiagPtr> theBornToRealDiagrams;
 
  map<tcDiagPtr,tcDiagPtr> theRealToBornDiagrams;
 
  RealEmissionKey lastRealEmissionKey;
 
  UnderlyingBornKey lastUnderlyingBornKey;
 
  multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator lastRealEmissionInfo;
 
  int theBornEmitter;
 
  int theBornSpectator;
 
  Energy theLastSubtractionScale;
 
  Energy theLastSplittingScale;
 
  Energy theLastSubtractionPt;
 
  Energy theLastSplittingPt;
 
  double theLastSubtractionZ;
 
  double theLastSplittingZ;
 
  Ptr<ShowerApproximation>::ptr theShowerApproximation;
 
  bool theRealShowerSubtraction;
 
  bool theVirtualShowerSubtraction;
 
  bool theLoopSimSubtraction;
 
  bool theRealEmissionScales;
 
  Energy theShowerHardScale;
 
  Energy theShowerScale;
 
  vector<double> theShowerParameters;
 
  bool theIsInShowerPhasespace;
 
  bool theIsAboveCutoff;
 
private:
 
  SubtractionDipole & operator=(const SubtractionDipole &) = delete;
 
};
 
}
 
#endif /* HERWIG_SubtractionDipole_H */