Dipole/Merging/Merger.h

0001   /// -*- C++ -*-
0002   //
0003   /// Merger.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004   /// Copyright (C) 2002-2019 The Herwig Collaboration
0005   //
0006   /// Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007   /// Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008   //
0009 #ifndef HERWIG_Merger_H
0010 #define HERWIG_Merger_H
0011   //
0012   /// This is the declaration of the Merger class.
0013   //
0014 #include "MergingFactory.fh"
0015 #include "Node.fh"
0016
0017
0018
0019 #include "ThePEG/Handlers/HandlerBase.h"
0020 #include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
0021   //#include "Herwig/Shower/Dipole/Base/DipoleSplittingGenerator.h"
0022 #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
0023
0024 #include "Herwig/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
0025 #include "Herwig/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
0026 #include "Herwig/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
0027 #include "Herwig/MatrixElement/Matchbox/Phasespace/IFMassiveTildeKinematics.h"
0028 #include "Herwig/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
0029 #include "Herwig/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
0030 #include "Herwig/MatrixElement/Matchbox/Phasespace/FIMassiveTildeKinematics.h"
0031
0032 #include "ThePEG/Cuts/JetFinder.h"
0033 #include "ThePEG/Cuts/Cuts.h"
0034
0035
0036
0037 namespace Herwig {
0038
0039   using namespace ThePEG;
0040
0041
0042   class Merger;
0043
0044   ThePEG_DECLARE_POINTERS(Merger , MergerPtr );
0045
0046   typedef vector<NodePtr> NodePtrVec;
0047     //definition of a history step
0048   struct HistoryStep {
0049       /// containing the full information
0050     NodePtr node;
0051       /// current sudakov weight of the history
0052     double weight;
0053       /// current scale of the history
0054     Energy scale;
0055   };
0056
0057   typedef vector< HistoryStep > History;
0058
0059   typedef multimap<DipoleIndex, Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap2;
0060
0061   /**
0062    * \ingroup DipoleShower
0063    * \author Johannes Bellm
0064    *
0065    * \brief This class is responsible for the handling of the merging process
0066    * after the setup stage, performed by the merging factory.
0067    * The class is inherited from the MergerBase class that is visible
0068    * (and the interface) to the shower and matrix elements.
0069    *
0070    * The main functions of the Merger are:
0071    * - matrixElementRegion:
0072    *           Given vectors of incoming and outgoing particles, the
0073    *           function defines the matrix element region.
0074    *           This function is used by the shower after an emission
0075    *           that produces a multiplicity still handeld by the merging
0076    *           and the Matrix element to determine the region to calculate the
0077    *           full ME contribution.
0078    * - mergingScale:
0079    *           This is the possibly smeared version of the merging pt
0080    *           (centralMergePt) given as input to the Merger. This scaleis used
0081    *           in the definitions of ME regions.
0082    * - N and M:
0083    *           are the number of additional legs of the highest LO or NLO (virtual)
0084    *           contribution.
0085    * - MergingDSigDR:
0086    *           The main function to call the calculation of the underlying
0087    *           process and the history reweights.
0088    *           Here it is also decided if the process should be
0089    *           clustered/projected to a process with less legs.
0090    *
0091    *
0092    * @see \ref MergerInterfaces "The interfaces"
0093    * defined for Merger.
0094    */
0095   class Merger: public MergerBase {
0096
0097     friend class MergingFactory;
0098     friend class Node;
0099
0100   public:
0101
0102       // define the ME region for a particle vector.
0103     bool matrixElementRegion(PVector incoming,
0104                              PVector outgoing,
0105                              Energy winnerScale = ZERO,
0106                              Energy cutscale = ZERO)const;
0107       /// return the current merging scale,
0108       /// gets smeared around the central merging scale in generate kinematics.
0109     Energy mergingScale()const{return theMergePt;}
0110       /// return the current merging pt (should be unified with mergingScale)
0111     Energy mergePt()const {return theMergePt;}
0112       /// legsize of highest process with NLO corrections
0113     int M()const;
0114       /// legsize of the highest LO merged process
0115     int N()const;
0116       /// legsize of the production process
0117     int N0()const{return theN0;}
0118       /// cross section of as given by the merging
0119     CrossSection MergingDSigDR();
0120       /// ***** virtual functions of the base class ****///
0121       /// set the current xcomb, called from ME
0122     void setXComb( tStdXCombPtr );
0123       /// set kinematics, called from ME
0124     void setKinematics();
0125       ///  clear kinematics, called from ME
0126     void clearKinematics();
0127       /// generate kinematics, called from ME
0128     bool generateKinematics( const double * );
0129       /// generate kinematics, called from ME
0130     void flushCaches();
0131       /// return the current maximum legs, the shower should veto
0132     size_t maxLegs() const {return theCurrentMaxLegs;}
0133       /// set the current ME
0134     void setME(MatchboxMEBasePtr me){
0135       theCurrentME=me;
0136       assert(theFirstNodeMap.count(theCurrentME));
0137       theCurrentNode=theFirstNodeMap[theCurrentME];
0138     }
0139
0140   protected:
0141       /// the merging factory needs to set the legsize of the production process
0142     void N0(int n){ theN0=n;}
0143       /// return the large-N basis (TODO: implement check if born ME works with the choice)
0144     Ptr<ColourBasis>::ptr largeNBasis()const{return theLargeNBasis;}
0145       /// smear the merging pt
0146     void smearMergePt(){
0147         const double factor = 1. + (-1. + 2.*UseRandom::rnd() ) * smear();
0148         theMergePt = factor * centralMergePt();
0149     }
0150       /// true if the phase space for initial emissions should not be restricted in z.
0151     int openZBoundaries()const{return DSH()->showerPhaseSpaceOption();}
0152       /// return the current ME
0153     MatchboxMEBasePtr currentME() const { return theCurrentME; }
0154       /// return the current Node
0155     NodePtr currentNode() const { return theCurrentNode; }
0156
0157   private:
0158       /// calculate a single sudakov step for a given dipole
0159     double singlesudakov(Dipole, Energy, Energy, pair<bool, bool>);
0160       /// calculate the sudakov supression for a clusternode between
0161       /// the current running scale and next scale
0162     bool   dosudakov(NodePtr Born, Energy  running, Energy next, double& sudakov0_n);
0163       /// cleanup
0164     void   cleanup(NodePtr);
0165       /// return true if the cluster node has the matching number of
0166       /// legs to the current projector stage
0167     bool   isProjectorStage( NodePtr , int )const;
0168       /**
0169        * Calculate the staring scale:
0170        * if Node is part of the production process, calculate according to the
0171        * scale choice object in the merging scale objekt, else
0172        * return max(scale as scalechoice , min(Mass(i, j)))
0173        */
0174     Energy CKKW_StartScale(NodePtr) const;
0175       /// prepare the sudakov calculation
0176     void   CKKW_PrepareSudakov(NodePtr, Energy);
0177       /// number of active flavours as given by the shower
0178     double Nf(Energy scale)const{return DSH()->Nf(scale);}
0179       /// pointer to the factory
0180     MergingFactoryPtr treefactory() const;
0181       /// map from ME to first clusternode
0182     map<MatchboxMEBasePtr, NodePtr> firstNodeMap() const ;
0183       /// set the current merging pt, smeared in generate kinematics
0184     void mergePt(Energy x) {theMergePt = x;}
0185       /// return the central merging pt
0186     Energy centralMergePt() const {return theCentralMergePt;}
0187
0188   private:
0189       /// calculate the history weighted born cross section
0190     CrossSection MergingDSigDRBornStandard();
0191       /// calculate the history weighted virtual contribution
0192     CrossSection MergingDSigDRVirtualStandard();
0193       /**
0194        * calculate the history weighted real contribution
0195        * splitted into 3 differnt contibutions
0196        */
0197     CrossSection MergingDSigDRRealStandard();
0198       /// calculate the history weighted real contribution
0199       /// all dipoles above:
0200       /// N*(R rnd(i)-Dip_i) history_i U(\phi^n_i)
0201     CrossSection MergingDSigDRRealAllAbove();
0202       /// calculate the history weighted real contribution
0203       /// not all dipoles above:
0204       /// (R - sum PS_i) history_rnd U(\phi^n+1)
0205     CrossSection MergingDSigDRRealBelowSubReal();
0206       /// calculate the history weighted real contribution
0207       /// not all dipoles above:
0208       /// rnd(i)-> N*(PS_i - Dip_i) history_i U(\phi^n_i)
0209     CrossSection MergingDSigDRRealBelowSubInt();
0210       /// max legssize the shower should veto for LO
0211     size_t maxLegsLO() const {return N0()+N();}
0212       /// Calculate the LO partonic cross section.
0213     CrossSection TreedSigDR(Energy startscale);
0214       /// fill the projecting xcomb
0215     Energy fillProjector(int);
0216       /// fill the history, including calculation of sudakov supression
0217     void   fillHistory(Energy, NodePtr, NodePtr );
0218       /// calculate the pdf ratio for the given clusternode
0219     double pdfratio(NodePtr, Energy, Energy, int, bool fromIsME, bool toIsME);
0220       /// return the pdf-ratio reweight for the history
0221     double pdfReweight();
0222       /// return the alpha_s reweight for the history
0223     double alphaReweight(bool nocmw=false);
0224       /// max legssize the shower should veto for NLO
0225     size_t maxLegsNLO()const {return N0()+M();}
0226       /// calculate the virtual contribution.
0227     CrossSection LoopdSigDR(Energy startscale );
0228       /// calculate alpha_s expansion of the pdf-ratios
0229     double sumPdfReweightExpansion()const;
0230       /// calculate alpha_s expansion of the alpha_s-ratios, including K_g
0231     double sumAlphaSReweightExpansion()const;
0232       /// calculate alpha_s expansion of the sudakov exponents
0233     double sumFillHistoryExpansion();
0234       /// calculate alpha_s expansion of the single step alpha_s-ratio, including K_g
0235     double alphasExpansion( Energy next, Energy fixedScale)const;
0236       /// calculate alpha_s expansion of the single step pdf-ratio
0237     double pdfExpansion(NodePtr, int, Energy, Energy, double, int, Energy)const;
0238       /// calculate alpha_s expansion of the single step sudakov exponent
0239     bool   doHistExpansion(NodePtr Born, Energy  running, Energy next, Energy fixedScale, double& HistExpansion);
0240       /// calculate alpha_s expansion of the single dipole sudakov exponent
0241     double singleHistExpansion(Dipole, Energy, Energy, Energy, pair<bool, bool>);
0242       //alpha_s as given in the shower
0243     double as(Energy q)const{return DSH()->as(q);}
0244       // set the pointer to the Mergingfactory.
0245     void setFactory(MergingFactoryPtr f){theTreeFactory=f;}
0246       // set the pointer to the DipoleShower.
0247     void setDipoleShower(DipoleShowerHandlerPtr dsh){theDipoleShowerHandler=dsh;}
0248       //return the dipole shower handler
0249     DipoleShowerHandlerPtr DSH(){return theDipoleShowerHandler;}
0250       //return the const dipole shower handler
0251     cDipoleShowerHandlerPtr DSH()const{return theDipoleShowerHandler;}
0252       /// insert map from ME to first clusternode
0253     void firstNodeMap(MatchboxMEBasePtr, NodePtr);
0254       /// history choice: weighted history choice
0255     int chooseHistory()const {return theChooseHistory;}
0256       /// the smearing factor for the merging scale
0257     double smear()const{return theSmearing;}
0258       /// return the large-N colour basis
0259     void largeNBasis(Ptr<ColourBasis>::ptr x){theLargeNBasis=x;}
0260       /// helper function to check the only multi condition.
0261     bool notOnlyMulti()const;
0262       /// Calculate the CMW AlphaS
0263     double cmwAlphaS(Energy q)const;
0264       /// debug output for virtual
0265     void debugVirt(double, double, double, double, CrossSection,
0266                    double, double, double, NodePtr,CrossSection) const;
0267       /// debug output for reals
0268     void debugReal( string, double, CrossSection, CrossSection) const;
0269
0270   private:
0271
0272       /// calculate the history expansion
0273     unsigned int theShowerExpansionWeights = 2;
0274       /// use CMW scheme
0275     unsigned int theCMWScheme = 0;
0276       /// true if current point should be projected
0277     bool projected = true;
0278       /// true if LO cross sections should be unitarised
0279     bool isUnitarized = true;
0280       /// true if NLO contributions should be unitarised
0281     bool isNLOUnitarized = true;
0282       /// history weight choice
0283     int theChooseHistory = 0;
0284       /// legsize of production process
0285     int theN0 = 0;
0286       /// calculate only the N particle contribution
0287     int theOnlyN = -1;
0288       /// the current maxlegs (either LO or NLO maxlegs)
0289     int theCurrentMaxLegs = -1;
0290       /// smearing factor for merging scale
0291     double theSmearing = 0.;
0292       /// cutoff for real emission contribution
0293     Energy theIRSafePT = 1_GeV;
0294       /// current merging scale
0295     Energy theMergePt = 4_GeV;
0296       /// central merging scale
0297     Energy theCentralMergePt = 4_GeV;
0298       /// below mergingscale/theRealSubtractionRatio the dipoles are used to subtract.
0299       /// above the shower approximation is in use.
0300     double theRealSubtractionRatio=3.;
0301       /// current cluster histoy including sudakov weights
0302     History history;
0303       /// pointer to the large-N basis
0304     Ptr<ColourBasis>::ptr theLargeNBasis;
0305       /// current Node
0306     NodePtr theCurrentNode;
0307       /// current ME
0308     MatchboxMEBasePtr theCurrentME;
0309       /// Tilde kinematics pointers, only to use lastPt(emitter, emission, spectator)
0310     Ptr<FFLightTildeKinematics>::ptr FFLTK = new_ptr( FFLightTildeKinematics() );
0311     Ptr<FILightTildeKinematics>::ptr FILTK = new_ptr( FILightTildeKinematics() );
0312     Ptr<IFLightTildeKinematics>::ptr IFLTK = new_ptr( IFLightTildeKinematics() );
0313     Ptr<IILightTildeKinematics>::ptr IILTK = new_ptr( IILightTildeKinematics() );
0314     Ptr<FFMassiveTildeKinematics>::ptr FFMTK = new_ptr( FFMassiveTildeKinematics() );
0315     Ptr<FIMassiveTildeKinematics>::ptr FIMTK = new_ptr( FIMassiveTildeKinematics() );
0316     Ptr<IFMassiveTildeKinematics>::ptr IFMTK = new_ptr( IFMassiveTildeKinematics() );
0317       //pointer to the shower handler
0318     DipoleShowerHandlerPtr theDipoleShowerHandler;
0319       /// pointer to the MergingFactory
0320     MergingFactoryPtr theTreeFactory;
0321       /// map from ME to first Node
0322     map<MatchboxMEBasePtr, NodePtr> theFirstNodeMap;
0323       /// map from ME to highest ME weight so far
0324     map<NodePtr, CrossSection> theHighMeWeightMap;
0325
0326   protected:
0327
0328     /** @name Standard Interfaced functions. */
0329       //@{
0330     /**
0331      * Initialize this object after the setup phase before saving an
0332      * EventGenerator to disk.
0333      * @throws InitException if object could not be initialized properly.
0334      */
0335     virtual void doinit();
0336
0337       //@}
0338
0339
0340   public:
0341
0342     /** @name Functions used by the persistent I/O system. */
0343       //@{
0344     /**
0345      * Function used to write out object persistently.
0346      * @param os the persistent output stream written to.
0347      */
0348     void persistentOutput(PersistentOStream & os) const;
0349
0350     /**
0351      * Function used to read in object persistently.
0352      * @param is the persistent input stream read from.
0353      * @param version the version number of the object when written.
0354      */
0355     void persistentInput(PersistentIStream & is, int version);
0356       //@}
0357
0358     /**
0359      * The standard Init function used to initialize the interfaces.
0360      * Called exactly once for each class by the class description system
0361      * before the main function starts or
0362      * when this class is dynamically loaded.
0363      */
0364     static void Init();
0365
0366   protected:
0367
0368     /** @name Clone Methods. */
0369       //@{
0370     /**
0371      * Make a simple clone of this object.
0372      * @return a pointer to the new object.
0373      */
0374     virtual IBPtr clone() const;
0375
0376     /** Make a clone of this object, possibly modifying the cloned object
0377      * to make it sane.
0378      * @return a pointer to the new object.
0379      */
0380     virtual IBPtr fullclone() const;
0381       //@}
0382
0383
0384       // If needed, insert declarations of virtual function defined in the
0385       // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0386
0387
0388   private:
0389
0390
0391
0392     /**
0393      * The assignment operator is private and must never be called.
0394      * In fact, it should not even be implemented.
0395      */
0396     Merger & operator=(const Merger &) = delete;
0397
0398   };
0399
0400 }
0401
0402
0403
0404 #endif /* HERWIG_Merger_H */