Matchbox/Phasespace/FlatInvertiblePhasespace.h

0001 // -*- C++ -*-
0002 //
0003 // FlatInvertiblePhasespace.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_FlatInvertiblePhasespace_H
0010 #define Herwig_FlatInvertiblePhasespace_H
0011 //
0012 // This is the declaration of the FlatInvertiblePhasespace class.
0013 //
0014
0015 #include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
0016
0017 namespace Herwig {
0018
0019 using namespace ThePEG;
0020
0021 /**
0022  * \ingroup Matchbox
0023  * \author Simon Platzer
0024  *
0025  * \brief FlatInvertiblePhasespace implements flat, invertible phase space generation.
0026  *
0027  */
0028 class FlatInvertiblePhasespace: public MatchboxPhasespace {
0029
0030 public:
0031
0032   /**
0033    * Generate a phase space point and return its weight.
0034    */
0035   virtual double generateTwoToNKinematics(const double*,
0036                       vector<Lorentz5Momentum>& momenta);
0037
0038   /**
0039    * Return the number of random numbers required to produce a given
0040    * multiplicity final state.
0041    */
0042   virtual int nDimPhasespace(int nFinal) const {
0043     if ( nFinal == 1 )
0044       return 1;
0045     return 3*nFinal - 4;
0046   }
0047
0048 public:
0049
0050   /**
0051    * Return true, if this phase space generator is invertible
0052    */
0053   virtual bool isInvertible() const { return true; }
0054
0055   /**
0056    * Invert the given phase space point to the random numbers which
0057    * would have generated it.
0058    */
0059   virtual double invertTwoToNKinematics(const vector<Lorentz5Momentum>& momenta,
0060                     double* r) const {
0061     return invertKinematics(momenta,(momenta[0]+momenta[1]).m(),r);
0062   }
0063
0064 private:
0065
0066   /**
0067    * Solve v = (n+2) * u^(n+1) - (n+1) * u^(n+2) for u
0068    */
0069   double bisect(double v, double n,
0070         double target = -16., double maxLevel = 80.) const;
0071
0072   /**
0073    * Return rho
0074    */
0075   double rho(Energy M, Energy N, Energy m) const {
0076     return sqrt((sqr(M)-sqr(N+m))*(sqr(M)-sqr(N-m)))/(8.*sqr(M));
0077   }
0078
0079   /**
0080    * Generate intermediate masses for a massless final state
0081    */
0082   double generateIntermediates(vector<Energy>& K,
0083                    const double* r) const;
0084
0085   /**
0086    * Invert intermediate masses for a massless final state
0087    */
0088   double invertIntermediates(const vector<Energy>& K,
0089                  double* r) const;
0090
0091   /**
0092    * Generate intermediate masses for a massive final state
0093    */
0094   double generateIntermediates(vector<Energy>& M,
0095                    const vector<Energy>& m,
0096                    const double* r) const;
0097
0098   /**
0099    * Invert intermediate masses for a massive final state
0100    */
0101   double invertIntermediates(const vector<Energy>& M,
0102                  const vector<Energy>& m,
0103                  double* r) const;
0104
0105   /**
0106    * Generate momenta in the CMS
0107    */
0108   double generateKinematics(vector<Lorentz5Momentum>& P,
0109                 Energy Ecm,
0110                 const double* r) const;
0111
0112   /**
0113    * Invert momenta in the CMS
0114    */
0115   double invertKinematics(const vector<Lorentz5Momentum>& P,
0116               Energy Ecm,
0117               double* r) const;
0118
0119   /**
0120    * Return the appropriate phase space weight,
0121    * Eq. 11 in 1308.2922
0122    * with the factor (2 pi)^4/(2 pi)^(3n) included
0123    * and the SHat of the process divided out to have everything expressed in the units of the ThePEG conventions, i.e.
0124    * without the Q^2 factor
0125    */
0126
0127   long double flatWeights(int n) const;
0128
0129 public:
0130
0131   /** @name Functions used by the persistent I/O system. */
0132   //@{
0133   /**
0134    * Function used to write out object persistently.
0135    * @param os the persistent output stream written to.
0136    */
0137   void persistentOutput(PersistentOStream & os) const;
0138
0139   /**
0140    * Function used to read in object persistently.
0141    * @param is the persistent input stream read from.
0142    * @param version the version number of the object when written.
0143    */
0144   void persistentInput(PersistentIStream & is, int version);
0145   //@}
0146
0147   /**
0148    * The standard Init function used to initialize the interfaces.
0149    * Called exactly once for each class by the class description system
0150    * before the main function starts or
0151    * when this class is dynamically loaded.
0152    */
0153   static void Init();
0154
0155 protected:
0156
0157   /** @name Clone Methods. */
0158   //@{
0159   /**
0160    * Make a simple clone of this object.
0161    * @return a pointer to the new object.
0162    */
0163   virtual IBPtr clone() const;
0164
0165   /** Make a clone of this object, possibly modifying the cloned object
0166    * to make it sane.
0167    * @return a pointer to the new object.
0168    */
0169   virtual IBPtr fullclone() const;
0170   //@}
0171
0172
0173 // If needed, insert declarations of virtual function defined in the
0174 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0175
0176
0177 private:
0178
0179   /**
0180    * The assignment operator is private and must never be called.
0181    * In fact, it should not even be implemented.
0182    */
0183   FlatInvertiblePhasespace & operator=(const FlatInvertiblePhasespace &) = delete;
0184
0185 };
0186
0187 }
0188
0189 #endif /* Herwig_FlatInvertiblePhasespace_H */