Merger.h

 
  //
  //
  //
#ifndef HERWIG_Merger_H
#define HERWIG_Merger_H
  //
  //
#include "MergingFactory.fh"
#include "Node.fh"
 
 
 
#include "ThePEG/Handlers/HandlerBase.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
  //#include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
 
#include "Herwig/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h"
 
#include "ThePEG/Cuts/JetFinder.h"
#include "ThePEG/Cuts/Cuts.h"
 
 
 
namespace Herwig {
  
  using namespace ThePEG;
  
  
  class Merger;
  
  ThePEG_DECLARE_POINTERS(Merger , MergerPtr );
  
  typedef vector<NodePtr> NodePtrVec;
    //definition of a history step
  struct HistoryStep {
    NodePtr node;
    double weight;
    Energy scale;
  };
  
  typedef vector< HistoryStep > History;
  
  typedef multimap<DipoleIndex, Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap2;
  
  class Merger: public MergerBase {
    
    friend class MergingFactory;
    friend class Node;
    
  public:
    
      // define the ME region for a particle vector.
    bool matrixElementRegion(PVector incoming, 
                             PVector outgoing, 
                             Energy winnerScale = ZERO, 
                             Energy cutscale = ZERO)const;
    Energy mergingScale()const{return theMergePt;}
    Energy mergePt()const {return theMergePt;}
    int M()const;
    int N()const;
    int N0()const{return theN0;}
    CrossSection MergingDSigDR();
    void setXComb( tStdXCombPtr );
    void setKinematics();
    void clearKinematics();
    bool generateKinematics( const double * );
    void flushCaches();
    size_t maxLegs() const {return theCurrentMaxLegs;}
    void setME(MatchboxMEBasePtr me){
      theCurrentME=me;
      assert(theFirstNodeMap.count(theCurrentME));
      theCurrentNode=theFirstNodeMap[theCurrentME];
    }
 
  protected:
    void N0(int n){ theN0=n;}
    Ptr<ColourBasis>::ptr largeNBasis()const{return theLargeNBasis;}
    void smearMergePt(){
        const double factor = 1. + (-1. + 2.*UseRandom::rnd() ) * smear();
        theMergePt = factor * centralMergePt();
    }
    int openZBoundaries()const{return DSH()->showerPhaseSpaceOption();}
    MatchboxMEBasePtr currentME() const { return theCurrentME; }
    NodePtr currentNode() const { return theCurrentNode; }
    
  private:
    double singlesudakov(Dipole, Energy, Energy, pair<bool, bool>);
    bool   dosudakov(NodePtr Born, Energy  running, Energy next, double& sudakov0_n);
    void   cleanup(NodePtr);
    bool   isProjectorStage( NodePtr , int )const;
    Energy CKKW_StartScale(NodePtr) const;
    void   CKKW_PrepareSudakov(NodePtr, Energy);
    double Nf(Energy scale)const{return DSH()->Nf(scale);}
    MergingFactoryPtr treefactory() const;
    map<MatchboxMEBasePtr, NodePtr> firstNodeMap() const ;
    void mergePt(Energy x) {theMergePt = x;}
    Energy centralMergePt() const {return theCentralMergePt;}
    
  private:
    CrossSection MergingDSigDRBornStandard();
    CrossSection MergingDSigDRVirtualStandard();
    CrossSection MergingDSigDRRealStandard();
    CrossSection MergingDSigDRRealAllAbove();
    CrossSection MergingDSigDRRealBelowSubReal();
    CrossSection MergingDSigDRRealBelowSubInt();
    size_t maxLegsLO() const {return N0()+N();}
    CrossSection TreedSigDR(Energy startscale);
    Energy fillProjector(int);
    void   fillHistory(Energy, NodePtr, NodePtr );
    double pdfratio(NodePtr, Energy, Energy, int, bool fromIsME, bool toIsME);
    double pdfReweight();
    double alphaReweight(bool nocmw=false);
    size_t maxLegsNLO()const {return N0()+M();}
    CrossSection LoopdSigDR(Energy startscale );
    double sumPdfReweightExpansion()const;
    double sumAlphaSReweightExpansion()const;
    double sumFillHistoryExpansion();
    double alphasExpansion( Energy next, Energy fixedScale)const;
    double pdfExpansion(NodePtr, int, Energy, Energy, double, int, Energy)const;
    bool   doHistExpansion(NodePtr Born, Energy  running, Energy next, Energy fixedScale, double& HistExpansion);
    double singleHistExpansion(Dipole, Energy, Energy, Energy, pair<bool, bool>);
      //alpha_s as given in the shower
    double as(Energy q)const{return DSH()->as(q);}
      // set the pointer to the Mergingfactory.
    void setFactory(MergingFactoryPtr f){theTreeFactory=f;}
      // set the pointer to the DipoleShower.
    void setDipoleShower(DipoleShowerHandlerPtr dsh){theDipoleShowerHandler=dsh;}
      //return the dipole shower handler
    DipoleShowerHandlerPtr DSH(){return theDipoleShowerHandler;}
      //return the const dipole shower handler
    cDipoleShowerHandlerPtr DSH()const{return theDipoleShowerHandler;}
    void firstNodeMap(MatchboxMEBasePtr, NodePtr);
    int chooseHistory()const {return theChooseHistory;}
    double smear()const{return theSmearing;}
    void largeNBasis(Ptr<ColourBasis>::ptr x){theLargeNBasis=x;}
    bool notOnlyMulti()const;
    double cmwAlphaS(Energy q)const;
    void debugVirt(double, double, double, double, CrossSection,
                   double, double, double, NodePtr,CrossSection) const;
    void debugReal( string, double, CrossSection, CrossSection) const;
    
  private:
    
    unsigned int theShowerExpansionWeights = 2;
    unsigned int theCMWScheme = 0;
    bool projected = true;
    bool isUnitarized = true;
    bool isNLOUnitarized = true;
    int theChooseHistory = 0;
    int theN0 = 0;
    int theOnlyN = -1;
    int theCurrentMaxLegs = -1;
    double theSmearing = 0.;
    Energy theIRSafePT = 1_GeV;
    Energy theMergePt = 4_GeV;
    Energy theCentralMergePt = 4_GeV;
    double theRealSubtractionRatio=3.;
    History history;
    Ptr<ColourBasis>::ptr theLargeNBasis;
    NodePtr theCurrentNode;
    MatchboxMEBasePtr theCurrentME;
    Ptr<FFLightTildeKinematics>::ptr FFLTK = new_ptr( FFLightTildeKinematics() );
    Ptr<FILightTildeKinematics>::ptr FILTK = new_ptr( FILightTildeKinematics() );
    Ptr<IFLightTildeKinematics>::ptr IFLTK = new_ptr( IFLightTildeKinematics() );
    Ptr<IILightTildeKinematics>::ptr IILTK = new_ptr( IILightTildeKinematics() );
    Ptr<FFMassiveTildeKinematics>::ptr FFMTK = new_ptr( FFMassiveTildeKinematics() );
    Ptr<FIMassiveTildeKinematics>::ptr FIMTK = new_ptr( FIMassiveTildeKinematics() );
    Ptr<IFMassiveTildeKinematics>::ptr IFMTK = new_ptr( IFMassiveTildeKinematics() );
      //pointer to the shower handler
    DipoleShowerHandlerPtr theDipoleShowerHandler;
    MergingFactoryPtr theTreeFactory;
    map<MatchboxMEBasePtr, NodePtr> theFirstNodeMap;
    map<NodePtr, CrossSection> theHighMeWeightMap;
    
  protected:
    
    virtual void doinit();
    
    
    
  public:
    
    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:
    
    
    
    Merger & operator=(const Merger &) = delete;
    
  };
  
}
 
 
 
#endif /* HERWIG_Merger_H */