Matchbox/Utility/ColourBasis.h

0001 // -*- C++ -*-
0002 //
0003 // ColourBasis.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_ColourBasis_H
0010 #define HERWIG_ColourBasis_H
0011 //
0012 // This is the declaration of the ColourBasis class.
0013 //
0014
0015 #include "ThePEG/Handlers/HandlerBase.h"
0016
0017 #include "ThePEG/MatrixElement/Tree2toNDiagram.h"
0018 #include "ThePEG/MatrixElement/MEBase.h"
0019
0020 #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
0021 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
0022
0023 #include <iterator>
0024 #include <tuple>
0025
0026 namespace Herwig {
0027
0028 using std::iterator_traits;
0029 using std::distance;
0030
0031 using namespace ThePEG;
0032
0033 using boost::numeric::ublas::matrix;
0034 using boost::numeric::ublas::symmetric_matrix;
0035 using boost::numeric::ublas::compressed_matrix;
0036 using boost::numeric::ublas::upper;
0037
0038 /**
0039  * \ingroup Matchbox
0040  * \author Simon Platzer
0041  *
0042  * \brief ColourBasis is an interface to a colour basis
0043  * implementation.
0044  *
0045  */
0046 class ColourBasis: public HandlerBase {
0047
0048 public:
0049
0050   /** @name Standard constructors and destructors. */
0051   //@{
0052   /**
0053    * The default constructor.
0054    */
0055   ColourBasis();
0056
0057   /**
0058    * The destructor.
0059    */
0060   virtual ~ColourBasis();
0061   //@}
0062
0063 public:
0064
0065   /**
0066    * Return the factory which produced this matrix element
0067    */
0068   Ptr<MatchboxFactory>::tptr factory() const;
0069
0070   /**
0071    * Clone this colour basis.
0072    */
0073   Ptr<ColourBasis>::ptr cloneMe() const {
0074     return dynamic_ptr_cast<Ptr<ColourBasis>::ptr>(clone());
0075   }
0076
0077   /**
0078    * Clear this colour basis
0079    */
0080   virtual void clear();
0081
0082   /**
0083    * Prepare for the given sub process and return the basis
0084    * dimensionality.
0085    */
0086   size_t prepare(const cPDVector&, bool);
0087
0088   /**
0089    * Prepare for the given diagrams.
0090    */
0091   size_t prepare(const MEBase::DiagramVector&, bool);
0092
0093   /**
0094    * Return the index map.
0095    */
0096   const map<cPDVector,map<size_t,size_t> >& indexMap() const { return theIndexMap; }
0097
0098   /**
0099    * Return a map of basis tensor indices to vectors identifying a
0100    * certain ordering corresponding to the given colour structure. May
0101    * not be supported by all colour basis implementations.
0102    */
0103   virtual map<size_t,vector<vector<size_t> > > basisList(const vector<PDT::Colour>&) const {
0104     return map<size_t,vector<vector<size_t> > >();
0105   }
0106
0107   /**
0108    * Given a physical subprocess, a colour to amplitude label map and
0109    * a basis tensor index, return an identifier of the ordering
0110    * coresponding to the given colour structure. This will only return
0111    * sensible results for colour bases which implement the basisList
0112    * query.
0113    */
0114   const string& orderingString(const cPDVector& sub,
0115                    const map<size_t,size_t>& colourToAmplitude,
0116                    size_t tensorId);
0117
0118   /**
0119    * Given a physical subprocess, a colour to amplitude label map and
0120    * a basis tensor index, return an identifier of the ordering
0121    * coresponding to the given colour structure. This will only return
0122    * sensible results for colour bases which implement the basisList
0123    * query.
0124    */
0125   const set<vector<size_t> >& ordering(const cPDVector& sub,
0126                        const map<size_t,size_t>& colourToAmplitude,
0127                        size_t tensorId, size_t shift = 0);
0128
0129   /**
0130    * For the given subprocess and amplitude vectors
0131    * calculate the amplitude squared.
0132    */
0133   double me2(const cPDVector&, const map<vector<int>,CVector>&) const;
0134
0135   /**
0136    * For the given subprocess and amplitude vectors
0137    * calculate the interference.
0138    */
0139   double interference(const cPDVector&,
0140               const map<vector<int>,CVector>&,
0141               const map<vector<int>,CVector>&) const;
0142
0143   /**
0144    * For the given subprocess and amplitude vector
0145    * calculate the colour correlated amplitude.
0146    */
0147   double colourCorrelatedME2(const pair<size_t,size_t>&,
0148                  const cPDVector&,
0149                  const map<vector<int>,CVector>&) const;
0150
0151   /**
0152    * For the given subprocess and amplitude vector
0153    * calculate the amplitude squared.
0154    */
0155   Complex interference(const cPDVector&,
0156                const CVector&, const CVector&) const;
0157
0158   /**
0159    * For the given subprocess and amplitude vector
0160    * calculate the colour correlated amplitude.
0161    */
0162   Complex colourCorrelatedInterference(const pair<size_t,size_t>&,
0163                        const cPDVector&,
0164                        const CVector&, const CVector&) const;
0165
0166   /**
0167    * For the given subprocess and amplitude given as amp amp^\dagger
0168    * calculate the amplitude squared.
0169    */
0170   double me2(const cPDVector&, const matrix<Complex>&) const;
0171
0172   /**
0173    * For the given subprocess and amplitude given as amp amp^\dagger
0174    * calculate the colour correlated amplitude.
0175    */
0176   double colourCorrelatedME2(const pair<size_t,size_t>&,
0177                  const cPDVector&,
0178                  const matrix<Complex>&) const;
0179
0180   /**
0181    * Return the scalar product matrix for the given process.
0182    */
0183   const symmetric_matrix<double,upper>& scalarProducts(const cPDVector&) const;
0184
0185   /**
0186    * Return the correlator matrix for the given process.
0187    */
0188   const symmetric_matrix<double,upper>& correlator(const cPDVector&,
0189                            const pair<size_t,size_t>&) const;
0190
0191   /**
0192    * Return true, if the colour basis is capable of assigning colour
0193    * flows.
0194    */
0195   virtual bool haveColourFlows() const { return false; }
0196
0197   /**
0198    * Return a Selector with possible colour geometries for the selected
0199    * diagram weighted by their relative probabilities.
0200    */
0201   Selector<const ColourLines *> colourGeometries(tcDiagPtr diag,
0202                          const map<vector<int>,CVector>& amps);
0203
0204   /**
0205    * Return the colour tensor used for the selected colour flow
0206    */
0207   size_t tensorIdFromFlow(tcDiagPtr diag, const ColourLines * cl);
0208
0209   /**
0210    * Match colour representation.
0211    */
0212   struct matchRep {
0213     PDT::Colour m;
0214     matchRep(PDT::Colour n)
0215       : m(n) {}
0216     bool operator()(PDT::Colour c) const {
0217       return c == m;
0218     }
0219   };
0220
0221   /**
0222    * Return true, if this basis is running in large-N mode
0223    */
0224   virtual bool largeN() const { return theLargeN; }
0225
0226   /**
0227    * Switch to large n
0228    */
0229   void doLargeN(bool yes = true) { theLargeN = yes; }
0230
0231   /**
0232    * Convert particle data to colour information
0233    */
0234   vector<PDT::Colour> projectColour(const cPDVector&) const;
0235
0236   /**
0237    * Perform a normal ordering of the external legs. This default
0238    * implementation assumes normal ordered legs as 3 3bar ... 3 3bar 8 ... 8
0239    * while removing all non-coloured particles.
0240    */
0241   virtual vector<PDT::Colour> normalOrder(const vector<PDT::Colour>&) const;
0242
0243   /**
0244    * Determine the mapping of process to colour indices and return the
0245    * normal ordered vector of colour indices
0246    */
0247   vector<PDT::Colour> normalOrderMap(const cPDVector& sub);
0248
0249   /**
0250    * Get the normal ordered legs
0251    */
0252   const vector<PDT::Colour>& normalOrderedLegs(const cPDVector& sub) const;
0253
0254   /**
0255    * Generate the emission/splitting map for basis labels from sub-process
0256    * information; we consider the splitting of leg ij from subprocess subFrom
0257    * with the legs not participating in the splitting relabelled according to
0258    * emissionMap.
0259    */
0260   const std::tuple<vector<PDT::Colour>,vector<PDT::Colour>,
0261            size_t,size_t,size_t,map<size_t,size_t> >&
0262   normalOrderEmissionMap(const cPDVector& subFrom,
0263              const cPDVector& subTo,
0264              size_t ij, size_t i, size_t j,
0265              const map<size_t,size_t>& emissionMap);
0266
0267   /**
0268    * Return the colour charge matrix representation for the given splitting ij
0269    * -> i,j with other legs relabbeled as indicated by the dictionary map.
0270    */
0271   const pair<compressed_matrix<double>,vector<pair<size_t,size_t> > >&
0272   charge(const cPDVector& subFrom,
0273      const cPDVector& subTo,
0274      size_t ij, size_t i, size_t j,
0275      const map<size_t,size_t>& emissionMap);
0276
0277   /**
0278    * Convert the legs to a string.
0279    */
0280   string file(const vector<PDT::Colour>&) const;
0281
0282   /**
0283    * Calculate T_i^\dagger X T_j
0284    */
0285   void chargeProduct(const compressed_matrix<double>& ti,
0286              const vector<pair<size_t,size_t> >& tiNonZero,
0287              const symmetric_matrix<double,upper>& X,
0288              const compressed_matrix<double>& tj,
0289              const vector<pair<size_t,size_t> >& tjNonZero,
0290              symmetric_matrix<double,upper>& result) const;
0291
0292   /**
0293    * Calculate T_i X T_j^\dagger
0294    */
0295   void chargeProductAdd(const compressed_matrix<double>& ti,
0296             const vector<pair<size_t,size_t> >& tiNonZero,
0297             const matrix<Complex>& X,
0298             const compressed_matrix<double>& tj,
0299             const vector<pair<size_t,size_t> >& tjNonZero,
0300             matrix<Complex>& result,
0301             double factor = 1.) const;
0302
0303 public:
0304
0305   /**
0306    * Find a coloured path from a to b within the given diagram.
0307    */
0308   static list<pair<int,bool> > colouredPath(pair<int,bool> a, pair<int,bool> b,
0309                         Ptr<Tree2toNDiagram>::tcptr);
0310
0311   /**
0312    * Get all colour flows for the given diagram.
0313    */
0314   static list<list<list<pair<int,bool> > > > colourFlows(Ptr<Tree2toNDiagram>::tcptr);
0315
0316   /**
0317    * Convert a flow to a string representation appropriate for
0318    * ColourLines
0319    */
0320   static string cfstring(const list<list<pair<int,bool> > >&);
0321
0322   /**
0323    * Returns a map of how the order of the vectors of a colour basis are
0324    * changed when the indices are changed.
0325    */
0326   virtual map<size_t,size_t> indexChange(const vector<PDT::Colour>&,
0327                      const size_t,
0328                      const map<size_t,size_t>&) const {
0329     map<size_t,size_t> aMap;
0330     return aMap;
0331   }
0332
0333
0334 protected:
0335
0336   /**
0337    * Prepare the basis for the normal ordered legs and return the
0338    * dimensionality of the basis.
0339    */
0340   virtual size_t prepareBasis(const vector<PDT::Colour>&) = 0;
0341
0342   /**
0343    * Return the scalar product of basis tensors labelled a and b in
0344    * the basis used for the given normal ordered legs.
0345    */
0346   virtual double scalarProduct(size_t a, size_t b,
0347                    const vector<PDT::Colour>& abBasis) const = 0;
0348
0349   /**
0350    * Return the matrix element of a colour charge or quark splitting
0351    * <c_{n+1,a}|T_i|c_{n,b}> between basis tensors a and b, with respect to
0352    * aBasis and bBasis; k and l index the splitting product's labels, and dict
0353    * encodes how legs not participating in the splitting are relabeled to the
0354    * larger basis.
0355    */
0356   virtual double tMatrixElement(size_t i, size_t a, size_t b,
0357                 const vector<PDT::Colour>& aBasis,
0358                 const vector<PDT::Colour>& bBasis,
0359                 size_t k, size_t l,
0360                 const map<size_t,size_t>& dict) const = 0;
0361
0362   /**
0363    * Return true, if a large-N colour connection exists for the
0364    * given external legs and basis tensor.
0365    */
0366   virtual bool colourConnected(const cPDVector&,
0367                    const vector<PDT::Colour>&,
0368                    const pair<int,bool>&,
0369                    const pair<int,bool>&,
0370                    size_t) const;
0371
0372   /**
0373    * Return true, if a large-N colour connection exists for the
0374    * given external legs and basis tensor.
0375    */
0376   virtual bool colourConnected(const vector<PDT::Colour>&,
0377                    int, int, size_t) const {
0378     return false;
0379   }
0380
0381   /**
0382    * Match up colour flows for given diagram to basis tensors.
0383    */
0384   vector<string> makeFlows(Ptr<Tree2toNDiagram>::tcptr, size_t) const;
0385
0386   /**
0387    * Return the colour line map.
0388    */
0389   map<Ptr<Tree2toNDiagram>::tcptr,vector<ColourLines*> >&
0390   colourLineMap();
0391
0392   /**
0393    * Update the colour line map for a given diagram.
0394    */
0395   void updateColourLines(Ptr<Tree2toNDiagram>::tcptr);
0396
0397
0398 public:
0399
0400   /** @name Functions used by the persistent I/O system. */
0401   //@{
0402   /**
0403    * Function used to write out object persistently.
0404    * @param os the persistent output stream written to.
0405    */
0406   void persistentOutput(PersistentOStream & os) const;
0407
0408   /**
0409    * Function used to read in object persistently.
0410    * @param is the persistent input stream read from.
0411    * @param version the version number of the object when written.
0412    */
0413   void persistentInput(PersistentIStream & is, int version);
0414   //@}
0415
0416   /**
0417    * The standard Init function used to initialize the interfaces.
0418    * Called exactly once for each class by the class description system
0419    * before the main function starts or
0420    * when this class is dynamically loaded.
0421    */
0422   static void Init();
0423
0424
0425 // If needed, insert declarations of virtual function defined in the
0426 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0427
0428
0429 protected:
0430
0431   /** @name Standard Interfaced functions. */
0432   //@{
0433   /**
0434    * Initialize this object after the setup phase before saving an
0435    * EventGenerator to disk.
0436    * @throws InitException if object could not be initialized properly.
0437    */
0438   virtual void doinit();
0439
0440   /**
0441    * Initialize this object. Called in the run phase just before
0442    * a run begins.
0443    */
0444   virtual void doinitrun();
0445
0446   /**
0447    * Finalize this object. Called in the run phase just after a
0448    * run has ended. Used eg. to write out statistics.
0449    */
0450   virtual void dofinish();
0451   //@}
0452
0453 private:
0454
0455   typedef map<vector<PDT::Colour>,symmetric_matrix<double,upper> >
0456   ScalarProductMap;
0457
0458   typedef map<vector<PDT::Colour>,map<pair<size_t,size_t>,symmetric_matrix<double,upper> > >
0459   CorrelatorMap;
0460
0461   typedef map<std::tuple<vector<PDT::Colour>,vector<PDT::Colour>,
0462              size_t,size_t,size_t,map<size_t,size_t> >,
0463           pair<compressed_matrix<double>,vector<pair<size_t,size_t> > > > TSMap;
0464
0465   /**
0466    * True, if this basis is running in large-N mode
0467    */
0468   bool theLargeN;
0469
0470   /**
0471    * Map external legs to normal ordered versions
0472    */
0473   map<cPDVector,vector<PDT::Colour> > theNormalOrderedLegs;
0474
0475   /**
0476    * Index mappings to normal order from given leg assignments,
0477    * indexed by the original leg assignment.
0478    */
0479   map<cPDVector,map<size_t,size_t> > theIndexMap;
0480
0481   /**
0482    * Translations of emission maps
0483    */
0484   map<std::tuple<cPDVector,cPDVector,
0485          size_t,size_t,size_t,map<size_t,size_t> >,
0486       std::tuple<vector<PDT::Colour>,vector<PDT::Colour>,
0487          size_t,size_t,size_t,map<size_t,size_t> > >
0488   theEmissionMaps;
0489
0490   /**
0491    * The scalar product matrix S_n = <c_{n,a}|c_{n,b}> , indexed
0492    * by normal ordered leg assignments.
0493    */
0494   ScalarProductMap theScalarProducts;
0495
0496   /**
0497    * The correlator matrices T_i\cdot T_j -> T_i^\dagger S_{n+1} T_j
0498    * with T_i = <c_{n+1,a}|T_i|c_{n,b}> indexed by the `n' basis
0499    * normal ordered legs and indices i,j
0500    */
0501   CorrelatorMap theCorrelators;
0502
0503   /**
0504    * Colour charge or quark splitting matrix representations
0505    */
0506   TSMap theCharges;
0507
0508   /**
0509    * Map diagrams to colour flows indexed by basis tensor.
0510    */
0511   map<Ptr<Tree2toNDiagram>::tcptr,vector<string> > theFlowMap;
0512
0513   /**
0514    * Map diagrams to colour line objects.
0515    */
0516   map<Ptr<Tree2toNDiagram>::tcptr,vector<ColourLines*> > theColourLineMap;
0517
0518   /**
0519    * Store ordering identifiers
0520    */
0521   map<cPDVector,map<size_t,string> > theOrderingStringIdentifiers;
0522
0523   /**
0524    * Store ordering identifiers
0525    */
0526   map<cPDVector,map<size_t,set<vector<size_t> > > > theOrderingIdentifiers;
0527
0528   /**
0529    * Write out yet unknown basis computations.
0530    */
0531   void writeBasis(const string& prefix = "") const;
0532
0533   /**
0534    * Read in the basis computation which are supposed to be known.
0535    */
0536   void readBasis();
0537
0538   /**
0539    * Read in the basis computation which are supposed to be known.
0540    */
0541   bool readBasis(const vector<PDT::Colour>&);
0542
0543   /**
0544    * Gather any implementation dependend details when reading a basis
0545    */
0546   virtual void readBasisDetails(const vector<PDT::Colour>&) {}
0547
0548   /**
0549    * Write out symmetric matrices.
0550    */
0551   void write(const symmetric_matrix<double,upper>&, ostream&) const;
0552
0553   /**
0554    * Read in symmetric matrices.
0555    */
0556   void read(symmetric_matrix<double,upper>&, istream&);
0557
0558   /**
0559    * Write out compressed matrices.
0560    */
0561   void write(const compressed_matrix<double>&, ostream&,
0562          const vector<pair<size_t,size_t> >&) const;
0563
0564   /**
0565    * Read in compressed matrices.
0566    */
0567   void read(compressed_matrix<double>&, istream&,
0568         vector<pair<size_t,size_t> >&);
0569
0570   /**
0571    * True, if an attempt to read in basis information has been
0572    * completed.
0573    */
0574   bool didRead;
0575
0576   /**
0577    * True, if an attempt to write out basis information has been
0578    * completed.
0579    */
0580   mutable bool didWrite;
0581
0582   /**
0583    * Temporary storage.
0584    */
0585   matrix<double> tmp;
0586
0587   /**
0588    * The search path
0589    */
0590   string theSearchPath;
0591
0592   /**
0593    * The assignment operator is private and must never be called.
0594    * In fact, it should not even be implemented.
0595    */
0596   ColourBasis & operator=(const ColourBasis &) = delete;
0597
0598 };
0599
0600 }
0601
0602 #endif /* HERWIG_ColourBasis_H */