Matchbox/Utility/DensityOperator.h

0001 // -*- C++ -*-
0002 //
0003 // DensityOperator.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 2 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef Herwig_DensityOperator_H
0010 #define Herwig_DensityOperator_H
0011 //
0012 // This is the declaration of the DensityOperator class.
0013 //
0014
0015 #include "ThePEG/Handlers/HandlerBase.h"
0016
0017 #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h"
0018
0019 #include <tuple>
0020 #include <boost/numeric/ublas/matrix.hpp>
0021 #include <boost/numeric/ublas/vector.hpp>
0022
0023 namespace Herwig {
0024
0025 using namespace ThePEG;
0026
0027 typedef boost::numeric::ublas::vector<Complex> CVector;
0028
0029
0030
0031 /**
0032  * Here is the documentation of the DensityOperator class.
0033  *
0034  * @see \ref DensityOperatorInterfaces "The interfaces"
0035  * defined for DensityOperator.
0036  */
0037 class DensityOperator: public HandlerBase {
0038
0039 public:
0040
0041   /**
0042    * The default constructor.
0043    */
0044   DensityOperator();
0045
0046 public:
0047
0048   /**
0049    * Clears theDensityOperatorMap.
0050    */
0051   void clear();
0052
0053   /**
0054    * Prepare for the given sub process.
0055    */
0056   void prepare(const cPDVector&);
0057
0058   /**
0059    * Fill the density operator for the given hard subprocess, summing over all
0060    * helicity configurations.
0061    */
0062   void fill(const Ptr<MatchboxXComb>::ptr,
0063         const cPDVector&, const vector<Lorentz5Momentum>& momenta);
0064
0065   /**
0066    * Evolve the density operator, by
0067    * M_{n+1} = -\sum_{i,k}{-4*pi*alpha_s/Ti2*V_{ij,k} T_{i,n}M_nT_{k,n}^\dag},
0068    * see arXiv:1206.0180 eq. (5), note that the pi*pj factor is assumed to be
0069    * included in V_{ij,k}.
0070    */
0071   void evolve(const map<pair<size_t,size_t>,Complex>& Vijk,
0072           const cPDVector& before,
0073           const cPDVector& after,
0074           const map<std::tuple<size_t,size_t,size_t>,map<size_t,size_t> >& emissionsMap,
0075           const bool splitAGluon,
0076           const bool initialGluonSplitting);
0077
0078   /**
0079    * Calculate the colour matrix element correction.
0080    * -(1+delta(i,gluon))/Ti^2 Tr(Sn+1 Ti Mn Tk^dagger)/Tr(Sn Mn)
0081    * where the bracket in front compensates for the gluon symmetry factor,
0082    * Ti^2 is C_f or C_a, Sn+1 is the matrix of scalar products, and
0083    * Ti is the radiation matrix.
0084    * The first arg contains (emitter index, spectator index, emission pid)
0085    *
0086    */
0087   double colourMatrixElementCorrection(const std::tuple<size_t,size_t,long>& ikemission,
0088                        const cPDVector& particles);
0089
0090   /**
0091    * Checking colour conservation for the colour matrix element corrections.
0092    */
0093   void colourConservation(const cPDVector& particles);
0094
0095   /**
0096    * Get the colour basis.
0097    */
0098   Ptr<ColourBasis>::tptr colourBasis() { return theColourBasis; }
0099
0100   /**
0101    * Get the colour basis.
0102    */
0103   const Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }
0104
0105   /**
0106    * Set the colour basis.
0107    */
0108   void colourBasis(Ptr<ColourBasis>::ptr ptr) { theColourBasis = ptr; }
0109
0110   /**
0111    * Get the correlator map.
0112    */
0113   const map<pair<vector<PDT::Colour>,pair<size_t,size_t> >,double>& correlatorMap() const {
0114     return theCorrelatorMap;
0115   }
0116
0117   /** @name Functions used by the persistent I/O system. */
0118   //@{
0119   /**
0120    * Function used to write out object persistently.
0121    * @param os the persistent output stream written to.
0122    */
0123   void persistentOutput(PersistentOStream & os) const;
0124
0125   /**
0126    * Function used to read in object persistently.
0127    * @param is the persistent input stream read from.
0128    * @param version the version number of the object when written.
0129    */
0130   void persistentInput(PersistentIStream & is, int version);
0131   //@}
0132
0133   /**
0134    * The standard Init function used to initialize the interfaces.
0135    * Called exactly once for each class by the class description system
0136    * before the main function starts or
0137    * when this class is dynamically loaded.
0138    */
0139   static void Init();
0140
0141 protected:
0142
0143   /** @name Clone Methods. */
0144   //@{
0145   /**
0146    * Make a simple clone of this object.
0147    * @return a pointer to the new object.
0148    */
0149   virtual IBPtr clone() const;
0150
0151   /** Make a clone of this object, possibly modifying the cloned object
0152    * to make it sane.
0153    * @return a pointer to the new object.
0154    */
0155   virtual IBPtr fullclone() const;
0156   //@}
0157
0158
0159 // If needed, insert declarations of virtual function defined in the
0160 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0161
0162
0163 private:
0164
0165   /**
0166   * Number of colours used in colourNorm.
0167   */
0168   double Nc;
0169
0170   /**
0171    * QCD vertex normalization.
0172    */
0173   double TR;
0174
0175   /**
0176    * Normalization of colour charges \mathbf{T}_{ij}^2.
0177    */
0178   double colourNorm(const cPDPtr particle);
0179
0180   /**
0181    * Fast evaluation of Tij*Mn, where a Tij is the matrix from ColourBasis::charge,
0182    * which is a sparse matrix, and Mn is the density operator, a dense matrix.
0183    *
0184    */
0185   matrix<Complex> prodSparseDense(const compressed_matrix<double>&,
0186                   const matrix<Complex>&);
0187   /**
0188    * Fast evaluation of TijMn*Tkdagger, where a TijMn is the result from the method
0189    * prodSparseDense, a dense matrix, and Tkdagger is the transponse conjugate of
0190    * the matrix from ColourBasis::charge, a sparse matrix.
0191    *
0192    */
0193   matrix<Complex> prodDenseSparse(const matrix<Complex>&,
0194                   const compressed_matrix<double>&);
0195
0196   /**
0197    * Boosts a vector of momenta to the rest frame of the initial pair
0198    * of particles (the first 2 elements of the argument vector). Returns
0199    * the boosted vectors
0200    */
0201   vector<Lorentz5Momentum> boostToRestFrame(const vector<Lorentz5Momentum>& momenta);
0202
0203   /**
0204    * Boosts a vector of momenta to the rest frame of the initial pair
0205    */
0206   bool compareMomentum(const Lorentz5Momentum& p, const Lorentz5Momentum& q);
0207
0208   /**
0209    * Mapping of colour structures to density operator matrices.
0210    *
0211    */
0212   map<vector<PDT::Colour>,matrix<Complex> > theDensityOperatorMap;
0213
0214   /**
0215    * Mapping of colour structures and legs to colour correlators.
0216    */
0217   map<pair<vector<PDT::Colour>,pair<size_t,size_t> >,double> theCorrelatorMap;
0218
0219   /**
0220    * A map from the hard subprocess particles to a map of amplitude colour
0221    * basis order to the normal ordered colour basis.
0222    */
0223   map<cPDVector, map<size_t,size_t> > theColourBasisToColourBasisMap;
0224
0225   /**
0226    * Colour basis used.
0227    */
0228   Ptr<ColourBasis>::ptr theColourBasis;
0229
0230   /**
0231    * The assignment operator is private and must never be called.
0232    * In fact, it should not even be implemented.
0233    */
0234   DensityOperator & operator=(const DensityOperator &) = delete;
0235
0236 };
0237
0238 }
0239
0240 #endif /* Herwig_DensityOperator_H */