Matchbox/Phasespace/MatchboxPhasespace.h

0001 // -*- C++ -*-
0002 //
0003 // MatchboxPhasespace.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_MatchboxPhasespace_H
0010 #define HERWIG_MatchboxPhasespace_H
0011 //
0012 // This is the declaration of the MatchboxPhasespace class.
0013 //
0014
0015 #include "ThePEG/Handlers/StandardXComb.h"
0016 #include "ThePEG/Handlers/HandlerBase.h"
0017 #include "ThePEG/MatrixElement/Tree2toNDiagram.h"
0018 #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
0019 #include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.fh"
0020 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
0021 #include "Herwig/MatrixElement/Matchbox/Phasespace/PhasespaceCouplings.h"
0022
0023 namespace Herwig {
0024
0025 using namespace ThePEG;
0026
0027 /**
0028  * \ingroup Matchbox
0029  * \author Simon Platzer
0030  *
0031  * \brief Wrap around a vector of random numbers to behave as a stream
0032  * of those.
0033  */
0034 struct StreamingRnd {
0035
0036   /**
0037    * The random numbers
0038    */
0039   const double* numbers;
0040
0041   /**
0042    * The number of random numbers available.
0043    */
0044   size_t nRnd;
0045
0046   /**
0047    * Default constructor.
0048    */
0049   StreamingRnd()
0050     : numbers(0), nRnd(0) {}
0051
0052   /**
0053    * Construct from random numbers.
0054    */
0055   explicit StreamingRnd(const double* newNumbers,
0056             size_t n)
0057     : numbers(newNumbers), nRnd(n) {}
0058
0059   /**
0060    * Return next random number
0061    */
0062   inline double operator()() {
0063     assert(numbers && nRnd > 0);
0064     const double ret = numbers[0];
0065     ++numbers; --nRnd;
0066     return ret;
0067   }
0068
0069 };
0070
0071 /**
0072  * \ingroup Matchbox
0073  * \author Simon Platzer
0074  *
0075  * \brief MatchboxPhasespace defines an abstract interface to a phase
0076  * space generator.
0077  *
0078  */
0079 class MatchboxPhasespace:
0080     public HandlerBase,
0081     public LastXCombInfo<StandardXComb>,
0082     public LastMatchboxXCombInfo {
0083
0084 public:
0085
0086   /**
0087    * The default constructor.
0088    */
0089   MatchboxPhasespace();
0090
0091 public:
0092
0093   /**
0094    * Set the XComb object steering the Born matrix
0095    * element this class represents virtual corrections to.
0096    */
0097   virtual void setXComb(tStdXCombPtr xc) {
0098     theLastXComb = xc;
0099     lastMatchboxXComb(xc);
0100   }
0101
0102   /**
0103    * Return the factory object
0104    */
0105   Ptr<MatchboxFactory>::tcptr factory() const;
0106
0107   /**
0108    * Return the process data object
0109    */
0110   Ptr<ProcessData>::tptr processData() const;
0111
0112   /**
0113    * Generate a phase space point and return its weight.
0114    */
0115   virtual double generateKinematics(const double* r,
0116                     vector<Lorentz5Momentum>& momenta);
0117
0118   /**
0119    * Generate a phase space point and return its weight.
0120    */
0121   virtual double generateTwoToNKinematics(const double*,
0122                       vector<Lorentz5Momentum>& momenta) = 0;
0123
0124   /**
0125    * Generate a 2 -> 1 phase space point and return its weight.
0126    */
0127   virtual double generateTwoToOneKinematics(const double*,
0128                         vector<Lorentz5Momentum>& momenta);
0129
0130   /**
0131    * Return the number of random numbers required to produce a given
0132    * multiplicity final state.
0133    */
0134   virtual int nDim(const cPDVector&) const;
0135
0136   /**
0137    * Return the number of random numbers required to produce a given
0138    * multiplicity final state.
0139    */
0140   virtual int nDimPhasespace(int nFinal) const = 0;
0141
0142   /**
0143    * Return true, if this phasespace generator will generate incoming
0144    * partons itself.
0145    */
0146   virtual bool haveX1X2() const { return false; }
0147
0148   /**
0149    * Return true, if this phase space generator expects
0150    * the incoming partons in their center-of-mass system
0151    */
0152   virtual bool wantCMS() const { return true; }
0153
0154   /**
0155    * True, if mass generators should be used instead of fixed masses
0156    */
0157   bool useMassGenerators() const { return theUseMassGenerators; }
0158
0159   /**
0160    * Fill a diagram selector for the last phase space point.
0161    */
0162   virtual Selector<MEBase::DiagramIndex> selectDiagrams(const MEBase::DiagramVector&) const;
0163
0164   /**
0165    * Return the momentum and weight appropriate to the given timelike
0166    * branch of the diagram.
0167    */
0168   pair<double,Lorentz5Momentum> timeLikeWeight(const Tree2toNDiagram& diag,
0169                            int branch, double flatCut) const;
0170
0171   /**
0172    * Return the weight appropriate to the given spacelike branch of
0173    * the diagram.
0174    */
0175   double spaceLikeWeight(const Tree2toNDiagram& diag,
0176              const Lorentz5Momentum& incoming,
0177              int branch, double flatCut) const;
0178
0179   /**
0180    * Return the weight appropriate to the given diagram.
0181    */
0182   double diagramWeight(const Tree2toNDiagram& diag) const {
0183     assert( !diagramWeights().empty() );
0184     return diagramWeights().find(diag.id())->second;
0185   }
0186
0187   /**
0188    * Fill the diagram weights.
0189    */
0190   void fillDiagramWeights(double flatCut = 0.0);
0191
0192   /**
0193    * Clear the diagram weights.
0194    */
0195   void clearDiagramWeights() {
0196     diagramWeights().clear();
0197   }
0198
0199   /**
0200    * Clone this phase space generator.
0201    */
0202   Ptr<MatchboxPhasespace>::ptr cloneMe() const {
0203     return dynamic_ptr_cast<Ptr<MatchboxPhasespace>::ptr>(clone());
0204   }
0205
0206   /**
0207    * Clone the dependencies, using a given prefix.
0208    */
0209   virtual void cloneDependencies(const std::string& prefix = "");
0210
0211 public:
0212
0213   /**
0214    * Return true, if this phase space generator is invertible
0215    */
0216   virtual bool isInvertible() const { return false; }
0217
0218   /**
0219    * Invert the given phase space point to the random numbers which
0220    * would have generated it.
0221    */
0222   virtual double invertKinematics(const vector<Lorentz5Momentum>& momenta,
0223                   double* r) const;
0224
0225   /**
0226    * Invert the given phase space point to the random numbers which
0227    * would have generated it.
0228    */
0229   virtual double invertTwoToNKinematics(const vector<Lorentz5Momentum>&,
0230                     double*) const {
0231     return 0.;
0232   }
0233
0234   /**
0235    * Invert the given 2 -> 1 phase space point to the random numbers which
0236    * would have generated it.
0237    */
0238   virtual double invertTwoToOneKinematics(const vector<Lorentz5Momentum>&, double*) const;
0239
0240 public:
0241
0242   /**
0243    * Limit phasespace generation to a given collinear or soft limit.
0244    */
0245   void singularLimit(size_t i, size_t j) {
0246     if ( i > j )
0247       swap(i,j);
0248     singularLimits().insert(make_pair(i,j));
0249   }
0250
0251   /**
0252    * Return the last matched singular limit.
0253    */
0254   const pair<size_t,size_t>& lastSingularIndices() const {
0255     assert(lastSingularLimit() != singularLimits().end());
0256     return *lastSingularLimit();
0257   }
0258
0259   /**
0260    * Return true, if constraints on phasespace generation have been met.
0261    */
0262   bool matchConstraints(const vector<Lorentz5Momentum>& momenta);
0263
0264 protected:
0265
0266   /**
0267    * Set a coupling for the given vertex; the convention is that all
0268    * legs are outgoing, and all possible crossings will be taken care
0269    * of. If not set, coupling weights default to one.
0270    */
0271   void setCoupling(long a, long b, long c,
0272            double coupling, bool includeCrossings = true);
0273
0274 public:
0275
0276   /** @name Functions used by the persistent I/O system. */
0277   //@{
0278   /**
0279    * Function used to write out object persistently.
0280    * @param os the persistent output stream written to.
0281    */
0282   void persistentOutput(PersistentOStream & os) const;
0283
0284   /**
0285    * Function used to read in object persistently.
0286    * @param is the persistent input stream read from.
0287    * @param version the version number of the object when written.
0288    */
0289   void persistentInput(PersistentIStream & is, int version);
0290   //@}
0291
0292 public:
0293
0294   /**
0295    * The standard Init function used to initialize the interfaces.
0296    * Called exactly once for each class by the class description system
0297    * before the main function starts or
0298    * when this class is dynamically loaded.
0299    */
0300   static void Init();
0301
0302
0303 // If needed, insert declarations of virtual function defined in the
0304 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0305
0306
0307 private:
0308
0309   /**
0310    * A cutoff below which a region is considered singular.
0311    */
0312   Energy singularCutoff;
0313
0314   /**
0315    * True, if mass generators should be used instead of fixed masses
0316    */
0317   bool theUseMassGenerators;
0318
0319   /**
0320    * Couplings to be used in diagram weighting
0321    */
0322   Ptr<PhasespaceCouplings>::ptr theCouplings;
0323
0324   /**
0325    * Interface function to setcoupling
0326    */
0327   string doSetCoupling(string);
0328
0329   /**
0330    * Interface function to setcoupling
0331    */
0332   string doSetPhysicalCoupling(string);
0333
0334   /**
0335    * The first id in a range of id's meant to denote fictitious
0336    * 'ghost' particles to be used by the diagram generator
0337    * in loop induced processes.
0338    */
0339   int theLoopParticleIdMin;
0340
0341   /**
0342    * The last id in a range of id's meant to denote fictitious
0343    * 'ghost' particles to be used by the diagram generator
0344    * in loop induced processes.
0345    */
0346   int theLoopParticleIdMax;
0347
0348   /**
0349    * The assignment operator is private and must never be called.
0350    * In fact, it should not even be implemented.
0351    */
0352   MatchboxPhasespace & operator=(const MatchboxPhasespace &) = delete;
0353
0354 };
0355
0356 }
0357
0358 #endif /* HERWIG_MatchboxPhasespace_H */