Matchbox/Base/MatchboxAmplitude.h

0001 // -*- C++ -*-
0002 //
0003 // MatchboxAmplitude.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_MatchboxAmplitude_H
0010 #define HERWIG_MatchboxAmplitude_H
0011 //
0012 // This is the declaration of the MatchboxAmplitude class.
0013 //
0014
0015 #include "ThePEG/MatrixElement/Amplitude.h"
0016 #include "ThePEG/Handlers/LastXCombInfo.h"
0017 #include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h"
0018 #include "Herwig/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
0019 #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
0020 #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
0021 #include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
0022 #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
0023 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
0024 #include "ThePEG/Persistency/PersistentOStream.h"
0025 #include "ThePEG/Persistency/PersistentIStream.h"
0026
0027 namespace Herwig {
0028
0029 using namespace ThePEG;
0030
0031 /**
0032  * \ingroup Matchbox
0033  * \author Simon Platzer
0034  *
0035  * \brief Process information with coupling order
0036  */
0037 struct Process {
0038
0039   PDVector legs;
0040   unsigned int orderInAlphaS;
0041   unsigned int orderInAlphaEW;
0042
0043   Process()
0044     : orderInAlphaS(0), orderInAlphaEW(0) {}
0045
0046   Process(const PDVector& p,
0047       unsigned int oas,
0048       unsigned int oae)
0049     : legs(p), orderInAlphaS(oas), orderInAlphaEW(oae) {}
0050
0051   bool operator==(const Process& other) const {
0052     return
0053       legs == other.legs &&
0054       orderInAlphaS == other.orderInAlphaS &&
0055       orderInAlphaEW == other.orderInAlphaEW;
0056   }
0057
0058   bool operator<(const Process& other) const {
0059     if ( orderInAlphaS != other.orderInAlphaS )
0060       return orderInAlphaS < other.orderInAlphaS;
0061     if ( orderInAlphaEW != other.orderInAlphaEW )
0062       return orderInAlphaEW < other.orderInAlphaEW;
0063     return legs < other.legs;
0064   }
0065
0066   void persistentOutput(PersistentOStream & os) const {
0067     os << legs << orderInAlphaS << orderInAlphaEW;
0068   }
0069
0070   void persistentInput(PersistentIStream & is) {
0071     is >> legs >> orderInAlphaS >> orderInAlphaEW;
0072   }
0073
0074 };
0075
0076 /**
0077  * \ingroup Matchbox
0078  * \author Simon Platzer
0079  *
0080  * \brief Enumerate the type of calculation required
0081  */
0082 namespace ProcessType {
0083
0084   enum Types {
0085
0086     treeME2 = 0,
0087     colourCorrelatedME2,
0088     spinColourCorrelatedME2,
0089     oneLoopInterference,
0090     loopInducedME2,
0091     spinCorrelatedME2
0092
0093   };
0094
0095 }
0096
0097
0098 /**
0099  * \ingroup Matchbox
0100  * \author Simon Platzer
0101  *
0102  * \brief MatchboxAmplitude is the base class for amplitude
0103  * implementations inside Matchbox.
0104  *
0105  * @see \ref MatchboxAmplitudeInterfaces "The interfaces"
0106  * defined for MatchboxAmplitude.
0107  */
0108 class MatchboxAmplitude:
0109     public Amplitude,
0110     public LastXCombInfo<StandardXComb>,
0111     public LastMatchboxXCombInfo {
0112
0113 public:
0114
0115   /**
0116    * The default constructor.
0117    */
0118   MatchboxAmplitude();
0119
0120 public:
0121
0122   /**
0123    * Return the amplitude. Needs to be implemented from
0124    * ThePEG::Amplitude but is actually ill-defined, as colours of the
0125    * external particles are not specified. To this extent, this
0126    * implementation just asserts.
0127    */
0128   virtual Complex value(const tcPDVector & particles,
0129             const vector<Lorentz5Momentum> & momenta,
0130             const vector<int> & helicities);
0131
0132   /**
0133    * Return the factory which produced this matrix element
0134    */
0135   Ptr<MatchboxFactory>::tptr factory() const;
0136
0137   /** @name Subprocess information */
0138   //@{
0139
0140   /**
0141    * Return true, if this amplitude can handle the given process.
0142    */
0143   virtual bool canHandle(const PDVector& p,
0144              Ptr<MatchboxFactory>::tptr,
0145              bool) const { return canHandle(p); }
0146
0147   /**
0148    * Return true, if this amplitude can handle the given process.
0149    */
0150   virtual bool canHandle(const PDVector&) const { return false; }
0151
0152   /**
0153    * Return the number of random numbers required to evaluate this
0154    * amplitude at a fixed phase space point.
0155    */
0156   virtual int nDimAdditional() const { return 0; }
0157
0158   /**
0159    * Return a ME instance appropriate for this amplitude and the given
0160    * subprocesses
0161    */
0162   virtual Ptr<MatchboxMEBase>::ptr makeME(const PDVector&) const;
0163
0164   /**
0165    * Set the (tree-level) order in \f$g_S\f$ in which this matrix
0166    * element should be evaluated.
0167    */
0168   virtual void orderInGs(unsigned int) {}
0169
0170   /**
0171    * Return the (tree-level) order in \f$g_S\f$ in which this matrix
0172    * element is given.
0173    */
0174   virtual unsigned int orderInGs() const = 0;
0175
0176   /**
0177    * Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
0178    * element should be evaluated.
0179    */
0180   virtual void orderInGem(unsigned int) {}
0181
0182   /**
0183    * Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
0184    * element is given.
0185    */
0186   virtual unsigned int orderInGem() const = 0;
0187
0188   /**
0189    * Return the Herwig StandardModel object
0190    */
0191   Ptr<StandardModel>::tcptr standardModel() {
0192     if ( !hwStandardModel() )
0193       hwStandardModel(dynamic_ptr_cast<Ptr<StandardModel>::tcptr>(HandlerBase::standardModel()));
0194     return hwStandardModel();
0195   }
0196
0197   /**
0198    * Return true, if this amplitude already includes averaging over
0199    * incoming parton's quantum numbers.
0200    */
0201   virtual bool hasInitialAverage() const { return false; }
0202
0203   /**
0204    * Return true, if this amplitude already includes symmetry factors
0205    * for identical outgoing particles.
0206    */
0207   virtual bool hasFinalStateSymmetry() const { return false; }
0208
0209   /**
0210    * Return true, if this amplitude is handled by a BLHA one-loop provider
0211    */
0212   virtual bool isOLPTree() const { return false; }
0213
0214   /**
0215    * Return true, if this amplitude is handled by a BLHA one-loop provider
0216    */
0217   virtual bool isOLPLoop() const { return false; }
0218
0219   /**
0220    * Return true, if colour and spin correlated matrix elements should
0221    * be ordered from the OLP
0222    */
0223   virtual bool needsOLPCorrelators() const { return true; }
0224
0225   /**
0226    * Write the order file header
0227    */
0228   virtual void olpOrderFileHeader(ostream&) const;
0229
0230   /**
0231    * Write the order file process list
0232    */
0233   virtual void olpOrderFileProcesses(ostream&,
0234                      const map<pair<Process,int>,int>& procs) const;
0235
0236   /**
0237    * Start the one loop provider, if appropriate, giving order and
0238    * contract files
0239    */
0240   virtual void signOLP(const string&, const string&) { }
0241
0242   /**
0243    * Start the one loop provider, if appropriate
0244    */
0245   virtual void startOLP(const string&, int& status) { status = -1; }
0246
0247   /**
0248    * Start the one loop provider, if appropriate. This default
0249    * implementation writes an BLHA 2.0 order file and starts the OLP
0250    */
0251   virtual bool startOLP(const map<pair<Process,int>,int>& procs);
0252
0253   /**
0254    * Return true, if this amplitude needs to initialize an external
0255    * code.
0256    */
0257   virtual bool isExternal() const { return false; }
0258
0259   /**
0260    * Initialize this amplitude
0261    */
0262   virtual bool initializeExternal() { return false; }
0263
0264   /**
0265    * Return a generic process id for the given process
0266    */
0267   virtual int externalId(const cPDVector&) { return 0; }
0268
0269   /**
0270    * Return the map with masses to be used for amplitude evaluation
0271    */
0272   const map<long,Energy>& reshuffleMasses() const { return theReshuffleMasses; }
0273
0274   /**
0275    * Check if reshuffling is needed at all
0276    */
0277   void checkReshuffling(Ptr<MatchboxPhasespace>::tptr);
0278
0279   /**
0280    * Return true, if this amplitude makes use of amplitudeMomenta
0281    */
0282   virtual bool hasAmplitudeMomenta() const { return false; }
0283
0284   //@}
0285
0286   /** @name Colour basis. */
0287   //@{
0288
0289   /**
0290    * Return the colour basis.
0291    */
0292   virtual Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }
0293
0294   /**
0295    * Return true, if the colour basis is capable of assigning colour
0296    * flows.
0297    */
0298   virtual bool haveColourFlows() const {
0299     return colourBasis() ? colourBasis()->haveColourFlows() : false;
0300   }
0301
0302   /**
0303    * Return a Selector with possible colour geometries for the selected
0304    * diagram weighted by their relative probabilities.
0305    */
0306   virtual Selector<const ColourLines *> colourGeometries(tcDiagPtr diag) const;
0307
0308   /**
0309    * Return an ordering identifier for the current subprocess and
0310    * colour absis tensor index.
0311    */
0312   const string& colourOrderingString(size_t id) const;
0313
0314   /**
0315    * Return an ordering identifier for the current subprocess and
0316    * colour absis tensor index.
0317    */
0318   const set<vector<size_t> >& colourOrdering(size_t id) const;
0319
0320   //@}
0321
0322   /** @name Phasespace point, crossing and helicities */
0323   //@{
0324
0325   /**
0326    * Set the xcomb object.
0327    */
0328   virtual void setXComb(tStdXCombPtr xc);
0329
0330   /**
0331    * Return the momentum as crossed appropriate for this amplitude.
0332    */
0333   Lorentz5Momentum amplitudeMomentum(int) const;
0334
0335   /**
0336    * Perform a normal ordering of external legs and fill the
0337    * crossing information as. This default implementation sorts
0338    * lexicographically in (abs(colour)/spin/abs(charge)), putting pairs
0339    * of particles/anti-particles where possible.
0340    */
0341   virtual void fillCrossingMap(size_t shift = 0);
0342
0343   /**
0344    * Generate the helicity combinations.
0345    */
0346   virtual set<vector<int> > generateHelicities() const;
0347
0348   /**
0349    * Return the helicity combination of the physical process in the
0350    * conventions used by the spin correlation algorithm.
0351    */
0352   virtual vector<unsigned int> physicalHelicities(const vector<int>&) const;
0353
0354   //@}
0355
0356   /** @name Tree-level amplitudes */
0357   //@{
0358
0359   /**
0360    * Calculate the tree level amplitudes for the phasespace point
0361    * stored in lastXComb.
0362    */
0363   virtual void prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr);
0364
0365   /**
0366    * Return the matrix element squared.
0367    */
0368   virtual double me2() const;
0369
0370   /**
0371    * Return the colour charge of a given leg
0372    */
0373   double colourCharge(tcPDPtr) const;
0374
0375   /**
0376    * Return the large-N charge of a given leg
0377    */
0378   double largeNColourCharge(tcPDPtr) const;
0379
0380   /**
0381    * Return the largeN matrix element squared.
0382    */
0383   virtual double largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const;
0384
0385   /**
0386    * Return the colour correlated matrix element.
0387    */
0388   virtual double colourCorrelatedME2(pair<int,int> ij) const;
0389
0390   /**
0391    * Return the large-N colour correlated matrix element.
0392    */
0393   virtual double largeNColourCorrelatedME2(pair<int,int> ij,
0394                        Ptr<ColourBasis>::tptr largeNBasis) const;
0395
0396   /**
0397    * Return true if trivial colour configuration.
0398    */
0399   bool trivialColourLegs() const { return theTrivialColourLegs; }
0400
0401   /**
0402    * Return true, if this amplitude is capable of consistently filling
0403    * the rho matrices for the spin correllations
0404    */
0405   virtual bool canFillRhoMatrix() const { return false; }
0406
0407   /**
0408    * Return a positive helicity polarization vector for a gluon of
0409    * momentum p (with reference vector n) to be used when evaluating
0410    * spin correlations.
0411    */
0412   virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
0413                           const Lorentz5Momentum& n,
0414                           int id = -1) const;
0415
0416   /**
0417    * Return the colour and spin correlated matrix element.
0418    */
0419   virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
0420                      const SpinCorrelationTensor& c) const;
0421
0422   /**
0423    * Return the spin correlated matrix element.
0424    */
0425   virtual double spinCorrelatedME2(pair<int,int> emitterSpectator,
0426                    const SpinCorrelationTensor& c) const;
0427
0428
0429   /**
0430    * Return true, if tree-level contributions will be evaluated at amplitude level.
0431    */
0432   virtual bool treeAmplitudes() const { return true; }
0433
0434   /**
0435    * Evaluate the amplitude for the given colour tensor id and
0436    * helicity assignment
0437    */
0438   virtual Complex evaluate(size_t, const vector<int>&, Complex&) { return 0.; }
0439
0440   //@}
0441
0442   /** @name One-loop amplitudes */
0443   //@{
0444
0445   /**
0446    * Return the one-loop amplitude, if applicable.
0447    */
0448   virtual Ptr<MatchboxAmplitude>::tptr oneLoopAmplitude() const {
0449     return Ptr<MatchboxAmplitude>::tptr();
0450   }
0451
0452   /**
0453    * Diasble one-loop functionality if not needed.
0454    */
0455   virtual void disableOneLoop() {}
0456
0457   /**
0458    * Return true, if this amplitude is capable of calculating one-loop
0459    * (QCD) corrections.
0460    */
0461   virtual bool haveOneLoop() const { return false; }
0462
0463   /**
0464    * Return true, if this amplitude only provides
0465    * one-loop (QCD) corrections.
0466    */
0467   virtual bool onlyOneLoop() const { return false; }
0468
0469   /**
0470    * Return true, if one-loop contributions will be evaluated at amplitude level.
0471    */
0472   virtual bool oneLoopAmplitudes() const { return true; }
0473
0474   /**
0475    * Return true, if one loop corrections have been calculated in
0476    * dimensional reduction. Otherwise conventional dimensional
0477    * regularization is assumed. Note that renormalization is always
0478    * assumed to be MSbar.
0479    */
0480   virtual bool isDR() const { return false; }
0481
0482   /**
0483    * Return true, if the amplitude is DRbar renormalized, otherwise
0484    * MSbar is assumed.
0485    */
0486   virtual bool isDRbar() const { return true; }
0487
0488   /**
0489    * Return true, if one loop corrections are given in the conventions
0490    * of the integrated dipoles.
0491    */
0492   virtual bool isCS() const { return false; }
0493
0494   /**
0495    * Return true, if one loop corrections are given in the conventions
0496    * of BDK.
0497    */
0498   virtual bool isBDK() const { return false; }
0499
0500   /**
0501    * Return true, if one loop corrections are given in the conventions
0502    * of everything expanded.
0503    */
0504   virtual bool isExpanded() const { return false; }
0505
0506   /**
0507    * Return the value of the dimensional regularization
0508    * parameter. Note that renormalization scale dependence is fully
0509    * restored in DipoleIOperator.
0510    */
0511   virtual Energy2 mu2() const { return 0.*GeV2; }
0512
0513   /**
0514    * Indicate that this amplitude is running alphas by itself.
0515    */
0516   virtual bool hasRunningAlphaS() const { return false; }
0517
0518   /**
0519    * Indicate that this amplitude is running alphaew by itself.
0520    */
0521   virtual bool hasRunningAlphaEW() const { return false; }
0522
0523   /**
0524    * If defined, return the coefficient of the pole in epsilon^2
0525    */
0526   virtual double oneLoopDoublePole() const { return 0.; }
0527
0528   /**
0529    * If defined, return the coefficient of the pole in epsilon
0530    */
0531   virtual double oneLoopSinglePole() const { return 0.; }
0532
0533   /**
0534    * Calculate the one-loop amplitudes for the phasespace point
0535    * stored in lastXComb, if provided.
0536    */
0537   virtual void prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr);
0538
0539   /**
0540    * Return the one-loop/tree interference.
0541    */
0542   virtual double oneLoopInterference() const;
0543
0544   /**
0545    * Evaluate the amplitude for the given colour tensor id and
0546    * helicity assignment
0547    */
0548   virtual Complex evaluateOneLoop(size_t, const vector<int>&) { return 0.; }
0549
0550   //@}
0551
0552   /** @name Caching and helpers to setup amplitude objects. */
0553   //@{
0554
0555   /**
0556    * Flush all cashes.
0557    */
0558   virtual void flushCaches() {}
0559
0560   /**
0561    * Clone this amplitude.
0562    */
0563   Ptr<MatchboxAmplitude>::ptr cloneMe() const {
0564     return dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(clone());
0565   }
0566
0567   /**
0568    * Clone the dependencies, using a given prefix.
0569    */
0570   virtual void cloneDependencies(const std::string& prefix="" , bool slim=false);
0571
0572   //@}
0573
0574 public:
0575
0576   /** @name Functions used by the persistent I/O system. */
0577   //@{
0578   /**
0579    * Function used to write out object persistently.
0580    * @param os the persistent output stream written to.
0581    */
0582   void persistentOutput(PersistentOStream & os) const;
0583
0584   /**
0585    * Function used to read in object persistently.
0586    * @param is the persistent input stream read from.
0587    * @param version the version number of the object when written.
0588    */
0589   void persistentInput(PersistentIStream & is, int version);
0590   //@}
0591
0592   /**
0593    * The standard Init function used to initialize the interfaces.
0594    * Called exactly once for each class by the class description system
0595    * before the main function starts or
0596    * when this class is dynamically loaded.
0597    */
0598   static void Init();
0599
0600 // If needed, insert declarations of virtual function defined in the
0601 // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
0602
0603 protected:
0604
0605   /** @name Standard Interfaced functions. */
0606   //@{
0607
0608   /**
0609    * Initialize this object after the setup phase before saving an
0610    * EventGenerator to disk.
0611    * @throws InitException if object could not be initialized properly.
0612    */
0613   virtual void doinit();
0614
0615   /**
0616    * Initialize this object. Called in the run phase just before
0617    * a run begins.
0618    */
0619   virtual void doinitrun();
0620   //@}
0621
0622 private:
0623
0624   /**
0625    * Recursively generate helicities
0626    */
0627   void doGenerateHelicities(set<vector<int> >& res,
0628                 vector<int>& current,
0629                 size_t pos) const;
0630
0631   /**
0632    * The colour basis implementation to be used.
0633    */
0634   Ptr<ColourBasis>::ptr theColourBasis;
0635
0636   /**
0637    * The number of points after which helicity combinations wil be
0638    * cleaned up
0639    */
0640   int theCleanupAfter;
0641
0642   /**
0643    * The number of points that are calculated before a certain
0644    * helicity is excluded.  Needed in pp->V
0645    */
0646   int treeLevelHelicityPoints;
0647
0648   /**
0649    * The number of points that are calculated before a certain
0650    * helicity is excluded.  Needed in pp->V
0651    */
0652   int oneLoopHelicityPoints;
0653
0654   /**
0655    * The map with masses to be used for amplitude evaluation
0656    */
0657   map<long,Energy> theReshuffleMasses;
0658
0659   /**
0660    * True if trivial colour configuration.
0661    */
0662   bool theTrivialColourLegs;
0663
0664   /**
0665    * A command to fill the reshuffle mass map
0666    */
0667   string doReshuffle(string);
0668
0669   /**
0670    * A command to fill the reshuffle mass map
0671    */
0672   string doMassless(string);
0673
0674   /**
0675    * A command to fill the reshuffle mass map
0676    */
0677   string doOnShell(string);
0678
0679   /**
0680    * Clear the reshuffling map
0681    */
0682   string doClearReshuffling(string);
0683
0684   /**
0685    * The assignment operator is private and must never be called.
0686    * In fact, it should not even be implemented.
0687    */
0688   MatchboxAmplitude & operator=(const MatchboxAmplitude &) = delete;
0689
0690 };
0691
0692 inline PersistentOStream& operator<<(PersistentOStream& os,
0693                      const Process& h) {
0694   h.persistentOutput(os);
0695   return os;
0696 }
0697
0698 inline PersistentIStream& operator>>(PersistentIStream& is,
0699                      Process& h) {
0700   h.persistentInput(is);
0701   return is;
0702 }
0703
0704 }
0705
0706 #endif /* HERWIG_MatchboxAmplitude_H */