Matchbox/Base/SubtractedME.h

0001 // -*- C++ -*-
0002 //
0003 // SubtractedME.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_SubtractedME_H
0010 #define HERWIG_SubtractedME_H
0011 //
0012 // This is the declaration of the SubtractedME class.
0013 //
0014
0015 #include "ThePEG/MatrixElement/MEGroup.h"
0016 #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
0017 #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
0018 #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
0019 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
0020
0021 namespace Herwig {
0022
0023 using namespace ThePEG;
0024
0025 /**
0026  * \ingroup Matchbox
0027  * \author Simon Platzer
0028  *
0029  * \brief SubtractedME represents a subtracted real emission matrix element.
0030  *
0031  * @see \ref SubtractedMEInterfaces "The interfaces"
0032  * defined for SubtractedME.
0033  */
0034 class SubtractedME:
0035     public MEGroup,
0036     public LastMatchboxXCombInfo {
0037
0038 public:
0039
0040   /**
0041    * The default constructor.
0042    */
0043   SubtractedME();
0044
0045 public:
0046
0047   /**
0048    * Return the factory which produced this matrix element
0049    */
0050   Ptr<MatchboxFactory>::tcptr factory() const;
0051
0052   /** @name Phasespace and subprocess information */
0053   //@{
0054
0055   /**
0056    * For the given event generation setup return an xcomb object
0057    * appropriate to this matrix element.
0058    */
0059   virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
0060                 tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
0061                 tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
0062                 const PBPair & newPartonBins, tCutsPtr newCuts,
0063                 const DiagramVector & newDiagrams, bool mir,
0064                 const PartonPairVec& allPBins,
0065                 tStdXCombPtr newHead = tStdXCombPtr(),
0066                 tMEPtr newME = tMEPtr());
0067
0068   /**
0069    * Set the XComb object to be used in the next call to
0070    * generateKinematics() and dSigHatDR().
0071    */
0072   virtual void setXComb(tStdXCombPtr);
0073
0074   /**
0075    * Return true, if the same additional random numbers
0076    * should be presented to any of the dependent
0077    * matrix elements.
0078    */
0079   virtual bool uniformAdditional() const { return true; }
0080
0081   /**
0082    * Return true, if the XComb steering this matrix element
0083    * should keep track of the random numbers used to generate
0084    * the last phase space point
0085    */
0086   virtual bool keepRandomNumbers() const { return true; }
0087
0088   /**
0089    * Given a process from the head matrix element,
0090    * return a list of diagrams which should be considered for
0091    * the given dependent matrix element.
0092    */
0093   virtual MEBase::DiagramVector dependentDiagrams(const cPDVector& proc,
0094                           tMEPtr depME) const;
0095
0096   /**
0097    * Return true, if SubProcessGroups should be
0098    * setup from this MEGroup. If not, a single SubProcess
0099    * is constructed from the data provided by the
0100    * head matrix element.
0101    */
0102   virtual bool subProcessGroups() const;
0103
0104   /**
0105    * Return true, if one of the dependent subprocesses should be
0106    * constructed in place of the one driven by the head matrix element
0107    * or a full subprocess group.
0108    */
0109   virtual bool selectDependentSubProcess() const { return false; }
0110
0111   /**
0112    * Fill the projectors object of xcombs to choose subprocesses
0113    * different than the one currently integrated.
0114    */
0115   virtual void fillProjectors();
0116
0117   /**
0118    * Return true, if projectors will be used
0119    */
0120   virtual bool willProject() const {
0121     return virtualShowerSubtraction() || loopSimSubtraction();
0122   }
0123
0124   /**
0125    * Return true, if this MEGroup will reweight the contributing cross
0126    * sections.
0127    */
0128   virtual bool groupReweighted() const {
0129     return showerApproximation();
0130   }
0131
0132   /**
0133    * Reweight the head cross section
0134    */
0135   virtual double reweightHead(const vector<tStdXCombPtr>&);
0136
0137   /**
0138    * Reweight the dependent cross section
0139    */
0140   virtual double reweightDependent(tStdXCombPtr, const vector<tStdXCombPtr>&);
0141
0142   /**
0143    * Switch on or off that scales should be calculated from real emission kinematics
0144    */
0145   void doRealEmissionScales();
0146
0147   //@}
0148
0149   /** @name Methods relevant to matching */
0150   //@{
0151
0152   /**
0153    * Inform this matrix element that a new phase space
0154    * point is about to be generated, so all caches should
0155    * be flushed.
0156    */
0157   virtual void flushCaches() {
0158     MEGroup::flushCaches();
0159     if ( showerApproximation() )
0160       showerApproximation()->resetBelowCutoff();
0161   }
0162
0163   /**
0164    * Return the shower approximation.
0165    */
0166   Ptr<ShowerApproximation>::tptr showerApproximation() const;
0167
0168   /**
0169    * Indicate that the shower real emission contribution should be subtracted.
0170    */
0171   void doRealShowerSubtraction();
0172
0173   /**
0174    * Return true, if the shower real emission contribution should be subtracted.
0175    */
0176   bool realShowerSubtraction() const { return theRealShowerSubtraction; }
0177
0178   /**
0179    * Indicate that the shower virtual contribution should be subtracted.
0180    */
0181   void doVirtualShowerSubtraction();
0182
0183   /**
0184    * Return true, if the shower virtual contribution should be subtracted.
0185    */
0186   bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
0187
0188   /**
0189    * Indicate that the loopsim matched virtual contribution should be subtracted.
0190    */
0191   void doLoopSimSubtraction();
0192
0193   /**
0194    * Return true, if the loopsim matched virtual contribution should be subtracted.
0195    */
0196   bool loopSimSubtraction() const { return theLoopSimSubtraction; }
0197
0198   /**
0199    * Return true, if this configuration of cross sections should not
0200    * be included due to their relative magnitude. Arguments are head
0201    * cross section and dependent cross section, including all
0202    * reweights.
0203    */
0204   virtual bool discard(const CrossSection&, const CrossSection&) const { return false; }
0205
0206   //@}
0207
0208   /** @name Matrix element and dipole information */
0209   //@{
0210
0211   /**
0212    * Return the subtraction dipoles.
0213    */
0214   vector<Ptr<SubtractionDipole>::ptr> dipoles();
0215
0216   /**
0217    * Return the underlying born matrix elements.
0218    */
0219   const vector<Ptr<MatchboxMEBase>::ptr>& borns() const;
0220
0221   /**
0222    * Access the underlying born matrix elements,
0223    * overriding the ones contained in the factory object.
0224    */
0225   void setBorns(const vector<Ptr<MatchboxMEBase>::ptr>& newBorns) { theBorns = newBorns; }
0226
0227   /**
0228    * Build up dipoles needed.
0229    */
0230   void getDipoles();
0231
0232   /**
0233    * Clone all dipoles.
0234    */
0235   void cloneDipoles(const string& prefix = "");
0236
0237   /**
0238    * Clone the real emission matrix element.
0239    */
0240   void cloneRealME(const string& prefix = "");
0241
0242   /**
0243    * Clone all dependencies.
0244    */
0245   void cloneDependencies(const string& prefix = "") {
0246     cloneDipoles(prefix);
0247     cloneRealME(prefix);
0248   }
0249
0250   /**
0251    * Return all dipoles matching the given Born process
0252    */
0253   vector<Ptr<SubtractionDipole>::ptr> splitDipoles(const cPDVector&);
0254
0255   //@}
0256
0257   /** @name Veto scale settings */
0258   //@{
0259   /**
0260    * Set veto scales on the particles at the given
0261    * SubProcess which has been generated using this
0262    * matrix element.
0263    */
0264   virtual void setVetoScales(tSubProPtr) const;
0265   //@}
0266
0267   /** @name Diagnostic information */
0268   //@{
0269
0270   /**
0271    * Dump the setup to an ostream
0272    */
0273   void print(ostream&) const;
0274
0275   /**
0276    * Collect information on the last evaluated phasespace
0277    * point for verification or debugging purposes. This
0278    * only called, if the StdXCombGroup did accumulate
0279    * a non-zero cross section from this ME group.
0280    */
0281   virtual void lastEventStatistics();
0282
0283   /**
0284    * Print debug information on the last event
0285    */
0286   void printLastEvent(ostream&) const;
0287
0288   /**
0289    * Check the subtraction for the last event
0290    */
0291   void lastEventSubtraction();
0292
0293   /**
0294    * Return true, if verbose
0295    */
0296   bool verbose() const;
0297
0298   /**
0299    * Return true, if verbose
0300    */
0301   bool initVerbose() const;
0302
0303   //@}
0304
0305   /** @name Setup of Subtracted ME objects */
0306   //@{
0307
0308   /**
0309    * Return true if this object needs to be initialized before all
0310    * other objects (except those for which this function also returns
0311    * true).  This default version always returns false, but subclasses
0312    * may override it to return true.
0313    */
0314   virtual bool preInitialize() const { return true; }
0315
0316   /**
0317    * Simple envelope histogram to keep track of subtraction
0318    */
0319   struct SubtractionHistogram {
0320
0321     /**
0322      * The lower bound
0323      */
0324     double lower;
0325
0326     /**
0327      * The bins, indexed by upper bound.
0328      */
0329     map<double,pair<double,double> > bins;
0330
0331     /**
0332      * Constructor
0333      */
0334     SubtractionHistogram(double low = 0.001,
0335              double up = 10.,
0336              unsigned int nbins = 100);
0337
0338     /**
0339      * Book an event.
0340      */
0341     void book(double inv, double diff) {
0342       map<double,pair<double,double> >::iterator b =
0343     bins.upper_bound(inv);
0344       if ( b == bins.end() ) return;
0345       b->second.first = min(b->second.first,diff);
0346       b->second.second = max(b->second.second,diff);
0347     }
0348
0349     /**
0350      * Write to file given name and invariant.
0351      */
0352     void dump(const std::string& prefix,
0353         const int& plottype,
0354         const bool& scatterplot,
0355           const cPDVector& proc,
0356           int i, int j) const;
0357
0358     /**
0359      * Write to persistent ostream
0360      */
0361     void persistentOutput(PersistentOStream&) const;
0362
0363     /**
0364      * Read from persistent istream
0365      */
0366     void persistentInput(PersistentIStream&);
0367
0368   };
0369
0370   //@}
0371
0372 public:
0373
0374   /** @name Functions used by the persistent I/O system. */
0375   //@{
0376   /**
0377    * Function used to write out object persistently.
0378    * @param os the persistent output stream written to.
0379    */
0380   void persistentOutput(PersistentOStream & os) const;
0381
0382   /**
0383    * Function used to read in object persistently.
0384    * @param is the persistent input stream read from.
0385    * @param version the version number of the object when written.
0386    */
0387   void persistentInput(PersistentIStream & is, int version);
0388   //@}
0389
0390   /**
0391    * The standard Init function used to initialize the interfaces.
0392    * Called exactly once for each class by the class description system
0393    * before the main function starts or
0394    * when this class is dynamically loaded.
0395    */
0396   static void Init();
0397
0398 protected:
0399
0400   /** @name Clone Methods. */
0401   //@{
0402   /**
0403    * Make a simple clone of this object.
0404    * @return a pointer to the new object.
0405    */
0406   virtual IBPtr clone() const;
0407
0408   /** Make a clone of this object, possibly modifying the cloned object
0409    * to make it sane.
0410    * @return a pointer to the new object.
0411    */
0412   virtual IBPtr fullclone() const;
0413   //@}
0414
0415
0416 // If needed, insert declarations of virtual function defined in the
0417 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0418
0419
0420 protected:
0421
0422   /** @name Standard Interfaced functions. */
0423   //@{
0424   /**
0425    * Initialize this object after the setup phase before saving an
0426    * EventGenerator to disk.
0427    * @throws InitException if object could not be initialized properly.
0428    */
0429   virtual void doinit();
0430
0431   /**
0432    * Initialize this object. Called in the run phase just before
0433    * a run begins.
0434    */
0435   virtual void doinitrun();
0436
0437   /**
0438    * Finalize this object. Called in the run phase just after a
0439    * run has ended. Used eg. to write out statistics.
0440    */
0441   virtual void dofinish();
0442
0443   //@}
0444
0445 private:
0446
0447   /**
0448    * The underlying born matrix elements, overriding the ones
0449    * contained in the factory object.
0450    */
0451   vector<Ptr<MatchboxMEBase>::ptr> theBorns;
0452
0453   /**
0454    * Pointer to the head real emission ME casted to a MatchboxMEBase
0455    * object.
0456    */
0457   Ptr<MatchboxMEBase>::ptr theReal;
0458
0459   /**
0460    * Define the key for the collinear subtraction data.
0461    */
0462   typedef pair<cPDVector,pair<size_t, size_t> > CollinearSubtractionIndex;
0463
0464   /**
0465    * subtraction data for collinear limits.
0466    */
0467   map<CollinearSubtractionIndex,SubtractionHistogram> collinearHistograms;
0468
0469   /**
0470    * names of files to which subtraction data is written for all phase space points individually
0471    */
0472   map<CollinearSubtractionIndex,string> fnamesCollinearSubtraction;
0473
0474   /**
0475    * Define the key for the soft subtraction data.
0476    */
0477   typedef pair<cPDVector,size_t> SoftSubtractionIndex;
0478
0479   /**
0480    * subtraction data for soft limits.
0481    */
0482   map<SoftSubtractionIndex,SubtractionHistogram> softHistograms;
0483
0484   /**
0485    * names of files to which subtraction data is written for all phase space points individually
0486    */
0487   map<SoftSubtractionIndex,string> fnamesSoftSubtraction;
0488
0489   /**
0490    * True, if the shower real emission contribution should be subtracted.
0491    */
0492   bool theRealShowerSubtraction;
0493
0494   /**
0495    * True, if the shower virtual contribution should be subtracted.
0496    */
0497   bool theVirtualShowerSubtraction;
0498
0499   /**
0500    * True, if the loopsim matched virtual contribution should be subtracted.
0501    */
0502   bool theLoopSimSubtraction;
0503
0504 private:
0505
0506   /**
0507    * The assignment operator is private and must never be called.
0508    * In fact, it should not even be implemented.
0509    */
0510   SubtractedME & operator=(const SubtractedME &) = delete;
0511
0512 };
0513
0514 inline PersistentOStream& operator<<(PersistentOStream& os,
0515                      const SubtractedME::SubtractionHistogram& h) {
0516   h.persistentOutput(os);
0517   return os;
0518 }
0519
0520 inline PersistentIStream& operator>>(PersistentIStream& is,
0521                      SubtractedME::SubtractionHistogram& h) {
0522   h.persistentInput(is);
0523   return is;
0524 }
0525
0526 }
0527
0528 #endif /* HERWIG_SubtractedME_H */