herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
MatchboxMEBase.h
1 // -*- C++ -*-
2 //
3 // MatchboxMEBase.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4 // Copyright (C) 2002-2019 The Herwig Collaboration
5 //
6 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef HERWIG_MatchboxMEBase_H
10 #define HERWIG_MatchboxMEBase_H
11 //
12 // This is the declaration of the MatchboxMEBase class.
13 //
14 
15 #include "ThePEG/MatrixElement/MEBase.h"
16 #include "Herwig/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
17 #include "Herwig/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
18 #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
19 #include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.h"
20 #include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
21 #include "Herwig/MatrixElement/Matchbox/Base/MatchboxReweightBase.h"
22 #include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
23 #include "Herwig/MatrixElement/Matchbox/Base/MergerBase.h"
24 #include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
25 #include "Herwig/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h"
26 #include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
27 #include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
28 #include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
29 
30 namespace Herwig {
31 
32 using namespace ThePEG;
33 
45  public MEBase, public LastMatchboxXCombInfo {
46 
47 public:
48 
55 
59  virtual ~MatchboxMEBase();
61 
62 public:
63 
67  Ptr<MatchboxFactory>::tptr factory() const;
68 
71 
75  const Process& subProcess() const { return theSubprocess; }
76 
80  Process& subProcess() { return theSubprocess; }
81 
85  Ptr<Tree2toNGenerator>::tptr diagramGenerator() const;
86 
90  Ptr<ProcessData>::tptr processData() const;
91 
99  virtual bool noMirror () const { return true; }
100 
104  virtual void getDiagrams() const;
105  using MEBase::getDiagrams;
106 
113  virtual Selector<DiagramIndex> diagrams(const DiagramVector &) const;
114  using MEBase::diagrams;
115 
121  colourGeometries(tcDiagPtr diag) const;
122 
127  virtual bool canFillRhoMatrix() const {
128  if ( matchboxAmplitude() )
129  return matchboxAmplitude()->canFillRhoMatrix();
130  return false;
131  }
132 
136  virtual void constructVertex(tSubProPtr) {}
137 
141  virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);
142 
147  virtual unsigned int orderInAlphaS() const;
148  using MEBase::orderInAlphaS;
149 
154  virtual unsigned int orderInAlphaEW() const;
155  using MEBase::orderInAlphaEW;
156 
161  virtual bool hasInitialAverage() const {
162  return matchboxAmplitude() ? matchboxAmplitude()->hasInitialAverage() : false;
163  }
164 
169  virtual bool hasFinalStateSymmetry() const {
170  return matchboxAmplitude() ? matchboxAmplitude()->hasFinalStateSymmetry() : false;
171  }
172 
173 
178  virtual unsigned int getNLight() const;
179 
185  virtual vector<long> getNLightJetVec() const;
186 
192  virtual vector<long> getNHeavyJetVec() const;
193 
199  virtual vector<long> getNLightProtonVec() const;
200 
204  virtual bool isOLPTree() const {
205  return matchboxAmplitude() ? matchboxAmplitude()->isOLPTree() : false;
206  }
207 
211  virtual bool isOLPLoop() const {
212  return matchboxAmplitude() ? matchboxAmplitude()->isOLPLoop() : false;
213  }
214 
219  virtual bool needsOLPCorrelators() const {
220  return matchboxAmplitude() ? matchboxAmplitude()->needsOLPCorrelators() : true;
221  }
222 
226  const vector<int>& olpProcess() const { return theOLPProcess; }
227 
231  void olpProcess(int pType, int id) {
232  if ( theOLPProcess.empty() )
233  theOLPProcess.resize(5,0);
234  theOLPProcess[pType] = id;
235  }
236 
241  bool noCorrelations() const {
242  return theNoCorrelations;
243  }
244 
250  theNoCorrelations = true;
251  }
252 
258  theNoCorrelations = false;
259  }
260 
262 
265 
269  Ptr<MatchboxPhasespace>::tptr phasespace() const { return thePhasespace; }
270 
274  void phasespace(Ptr<MatchboxPhasespace>::ptr ps) { thePhasespace = ps; }
275 
280  virtual void setXComb(tStdXCombPtr xc);
281 
287  virtual bool keepRandomNumbers() const { return true; }
288 
297  double generateIncomingPartons(const double* r1, const double* r2);
298 
309  virtual bool generateKinematics(const double * r);
310 
318  virtual void setKinematics();
319 
324  virtual void clearKinematics();
325 
330  virtual int nDim() const;
331 
336  virtual int nDimBorn() const;
337 
346  virtual bool haveX1X2() const {
347  return
348  (phasespace() ? phasespace()->haveX1X2() : false) ||
349  diagrams().front()->partons().size() == 3;
350  }
351 
356  virtual bool wantCMS() const {
357  return
358  (phasespace() ? phasespace()->wantCMS() : true) &&
359  diagrams().front()->partons().size() != 3; }
360 
365  const vector<Lorentz5Momentum>& lastMEMomenta() const { return meMomenta(); }
366 
370  vector<Lorentz5Momentum>& lastMEMomenta() { return meMomenta(); }
371 
372 
377  int legsize() const {return int(meMomenta().size());}
378 
380 
383 
387  void scaleChoice(Ptr<MatchboxScaleChoice>::ptr sc) { theScaleChoice = sc; }
388 
392  Ptr<MatchboxScaleChoice>::tptr scaleChoice() const { return theScaleChoice; }
393 
397  void setScale(Energy2 ren=ZERO,Energy2 fac=ZERO) const;
398 
402  virtual bool hasRunningAlphaS() const {
403  if ( matchboxAmplitude() )
404  return matchboxAmplitude()->hasRunningAlphaS();
405  return false;
406  }
407 
411  virtual bool hasRunningAlphaEW() const {
412  if ( matchboxAmplitude() )
413  return matchboxAmplitude()->hasRunningAlphaEW();
414  return false;
415  }
416 
421  virtual Energy2 scale() const { return lastScale(); }
422 
426  virtual Energy2 factorizationScale() const;
427 
431  virtual double factorizationScaleFactor() const;
432 
433 
437  virtual double facFac() const{return factorizationScaleFactor();}
438 
442  virtual Energy2 renormalizationScale() const;
443 
447  virtual double renormalizationScaleFactor() const;
448 
449 
453  virtual double renFac() const{return renormalizationScaleFactor();}
454 
458  virtual Energy2 renormalizationScaleQED() const;
459 
463  virtual Energy2 showerScale() const;
464 
470  virtual void setVetoScales(tSubProPtr) const;
471 
475  bool fixedCouplings() const;
476 
480  bool fixedQEDCouplings() const;
481 
487  virtual double alphaS() const { return lastAlphaS(); }
488 
494  virtual double alphaEM() const { return lastAlphaEM(); }
495 
500  virtual bool havePDFWeight1() const;
501 
506  virtual bool havePDFWeight2() const;
507 
511  void getPDFWeight(Energy2 factorizationScale = ZERO) const;
512 
516  double pdf1(Energy2 factorizationScale = ZERO,
517  double xEx = 1., double xFactor = 1.) const;
518 
522  double pdf2(Energy2 factorizationScale = ZERO,
523  double xEx = 1., double xFactor = 1.) const;
524 
526 
529 
533  Ptr<MatchboxAmplitude>::tptr matchboxAmplitude() const { return theAmplitude; }
534 
538  void matchboxAmplitude(Ptr<MatchboxAmplitude>::ptr amp) { theAmplitude = amp; }
539 
545  virtual double me2() const;
546 
552  virtual double largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const;
553 
557  virtual double finalStateSymmetry() const;
558 
563  double me2Norm(unsigned int addAlphaS = 0) const;
564 
569  virtual CrossSection dSigHatDR() const;
570 
574  CrossSection prefactor()const;
575 
579  CrossSection dSigHatDRB() const ;
580 
584  CrossSection dSigHatDRV() const ;
585 
589  CrossSection dSigHatDRI() const ;
590 
596  CrossSection dSigHatDRAlphaDiff(double alpha) const ;
597 
598 
600 
603 
607  virtual double oneLoopInterference() const;
608 
613  virtual bool haveOneLoop() const;
614 
619  virtual bool onlyOneLoop() const;
620 
625  virtual bool isDRbar() const;
626 
633  virtual bool isDR() const;
634 
639  virtual bool isCS() const;
640 
645  virtual bool isBDK() const;
646 
651  virtual bool isExpanded() const;
652 
658  virtual Energy2 mu2() const;
659 
663  virtual double oneLoopDoublePole() const;
664 
668  virtual double oneLoopSinglePole() const;
669 
673  bool checkPoles() const;
674 
679 
683  double lower;
684 
688  double upper;
689 
693  map<double,double> bins;
694 
698  unsigned long sameSign;
699 
703  unsigned long oppositeSign;
704 
708  unsigned long nans;
709 
713  unsigned long overflow;
714 
718  unsigned long underflow;
719 
723  AccuracyHistogram(double low = -40.,
724  double up = 0.,
725  unsigned int nbins = 80);
726 
730  void book(double a, double b);
731 
735  void dump(const std::string& folder, const std::string& prefix,
736  const cPDVector& proc) const;
737 
741  void persistentOutput(PersistentOStream&) const;
742 
746  void persistentInput(PersistentIStream&);
747 
748  };
749 
753  void logPoles() const;
754 
758  const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const {
759  return theVirtuals;
760  }
761 
765  vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() {
766  return theVirtuals;
767  }
768 
772  void doOneLoop() { theOneLoop = true; }
773 
778  void noOneLoop() { theOneLoop = false; }
779 
783  bool oneLoop() const { return theOneLoop; }
784 
789  void doOneLoopNoBorn() { theOneLoop = true; theOneLoopNoBorn = true; }
790 
791  void noOneLoopNoBorn() { theOneLoop = false; theOneLoopNoBorn = false; }
792 
797  bool oneLoopNoBorn() const { return theOneLoopNoBorn || onlyOneLoop(); }
798 
803  void doOneLoopNoLoops() { theOneLoop = true; theOneLoopNoLoops = true; }
804 
809  bool oneLoopNoLoops() const { return theOneLoopNoLoops; }
810 
812 
815 
823  vector<Ptr<SubtractionDipole>::ptr>
824  getDipoles(const vector<Ptr<SubtractionDipole>::ptr>&,
825  const vector<Ptr<MatchboxMEBase>::ptr>&,bool slim=false) const;
826 
827 
835  virtual bool noDipole(int,int,int) const { return false; }
836 
844  virtual bool noDipole(int,int) const { return false; }
845 
851  virtual double colourCorrelatedME2(pair<int,int>) const;
852 
859  virtual double largeNColourCorrelatedME2(pair<int,int> ij,
860  Ptr<ColourBasis>::tptr largeNBasis) const;
861 
867  virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
868  const SpinCorrelationTensor& c) const;
869 
875  virtual double spinCorrelatedME2(pair<int,int> emitterSpectator,
876  const SpinCorrelationTensor& c) const;
877 
879 
882 
888  virtual void flushCaches();
889 
893  bool verbose() const;
894 
898  bool initVerbose() const;
899 
903  void print(ostream&) const;
904 
908  virtual void printLastEvent(ostream&) const;
909 
914  void logGenerateKinematics(const double * r) const;
915 
920  void logSetScale() const;
921 
926  void logPDFWeight() const;
927 
932  void logME2() const;
933 
938  void logDSigHatDR() const;
939 
941 
944 
948  void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }
949 
953  const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }
954 
958  vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }
959 
963  const MergerBasePtr merger() const;
964 
968  MergerBasePtr merger() ;
969 
973  void merger(MergerBasePtr v);
974 
976 
979 
986  virtual bool preInitialize() const { return true; }
987 
991  Ptr<MatchboxMEBase>::ptr cloneMe() const {
992  return dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(clone());
993  }
994 
998  void cloneDependencies(const std::string& prefix = "",bool slim = false );
999 
1003  void prepareXComb(MatchboxXCombData&) const;
1004 
1009  virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
1010  tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
1011  tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
1012  const PBPair & newPartonBins, tCutsPtr newCuts,
1013  const DiagramVector & newDiagrams, bool mir,
1014  const PartonPairVec& allPBins,
1015  tStdXCombPtr newHead = tStdXCombPtr(),
1016  tMEPtr newME = tMEPtr());
1017 
1022  virtual StdXCombPtr makeXComb(tStdXCombPtr newHead,
1023  const PBPair & newPartonBins,
1024  const DiagramVector & newDiagrams,
1025  tMEPtr newME = tMEPtr());
1026 
1028 
1029 public:
1030 
1037  void persistentOutput(PersistentOStream & os) const;
1038 
1044  void persistentInput(PersistentIStream & is, int version);
1046 
1053  static void Init();
1054 
1055 protected:
1056 
1063  virtual IBPtr clone() const;
1064 
1069  virtual IBPtr fullclone() const;
1071 
1072 protected:
1073 
1081  virtual void doinit();
1082 
1087  virtual void doinitrun();
1088 
1093  virtual void dofinish();
1095 
1096 private:
1097 
1101  Ptr<MatchboxPhasespace>::ptr thePhasespace;
1102 
1106  Ptr<MatchboxAmplitude>::ptr theAmplitude;
1107 
1111  Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
1112 
1116  vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
1117 
1122  vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
1123 
1124 private:
1125 
1130 
1135 
1141 
1147 
1151  vector<int> theOLPProcess;
1152 
1156  mutable map<cPDVector,AccuracyHistogram> epsilonSquarePoleHistograms;
1157 
1161  mutable map<cPDVector,AccuracyHistogram> epsilonPoleHistograms;
1162 
1168 
1172  mutable pair<bool,bool> theHavePDFs;
1173 
1177  mutable bool checkedPDFs;
1178 
1184  MergerBasePtr theMerger;
1185 
1186 
1187 private:
1188 
1193  MatchboxMEBase & operator=(const MatchboxMEBase &) = delete;
1194 
1195 };
1196 
1199  h.persistentOutput(os);
1200  return os;
1201 }
1202 
1205  h.persistentInput(is);
1206  return is;
1207 }
1208 
1209 }
1210 
1211 #endif /* HERWIG_MatchboxMEBase_H */
void persistentInput(PersistentIStream &)
Read from persistent istream.
vector< Ptr< MatchboxReweightBase >::ptr > theReweights
A vector of reweight objects the sum of which should be applied to reweight this matrix element...
void needsCorrelations()
Indicate that this is a virtual matrix element which does require colour correlators.
const DiagramVector & diagrams() const
virtual double facFac() const
Get the factorization scale factor.
unsigned long nans
The number of points being nan or inf.
bool theOneLoopNoLoops
True, if this matrix element includes one-loop corrections but no actual loop contributions (e...
bool theNoCorrelations
True, if this is a real emission matrix element which does not require colour correlators.
void addReweight(Ptr< MatchboxReweightBase >::ptr rw)
Insert a reweight object.
MatchboxMEBase is the base class for matrix elements in the context of the matchbox NLO interface...
pair< bool, bool > theHavePDFs
Flag which pdfs should be included.
virtual bool haveX1X2() const
Return true, if this matrix element will generate momenta for the incoming partons itself...
map< cPDVector, AccuracyHistogram > epsilonPoleHistograms
Histograms of epsilon pole cancellation.
Process theSubprocess
The subprocess to be considered.
virtual bool noDipole(int, int, int) const
If this matrix element is considered a real emission matrix element, but actually neglecting a subcla...
bool oneLoop() const
Return true, if this matrix element includes one-loop corrections.
Process & subProcess()
Access the subprocess.
virtual double alphaEM() const
Return the value of associated with the phase space point provided by the last call to setKinematics...
ThePEG::Ptr< InterfacedBase >::pointer IBPtr
MergerBasePtr theMerger
The merging helper to be used.
bool theOneLoopNoBorn
True, if this matrix element includes one-loop corrections but no Born contributions.
virtual unsigned int orderInAlphaEW() const=0
virtual bool noDipole(int, int) const
If this matrix element is considered an underlying Born matrix element in the context of a subtracted...
void noOneLoop()
Instruct this matrix element not to include one-loop corrections.
unsigned long sameSign
The number of points of same sign.
Ptr< MatchboxPhasespace >::ptr thePhasespace
The phase space generator to be used.
Ptr< MatchboxAmplitude >::tptr matchboxAmplitude() const
Return the amplitude.
Ptr< MatchboxMEBase >::ptr cloneMe() const
Clone this matrix element.
vector< Lorentz5Momentum > & lastMEMomenta()
Access the meMomenta.
map< cPDVector, AccuracyHistogram > epsilonSquarePoleHistograms
Histograms of epsilon^2 pole cancellation.
void phasespace(Ptr< MatchboxPhasespace >::ptr ps)
Set the phase space generator to be used.
Ptr< MatchboxAmplitude >::ptr theAmplitude
The amplitude to be used.
void doOneLoopNoLoops()
Instruct this matrix element to include one-loop corrections but no actual loop contributions.
void persistentOutput(PersistentOStream &) const
Write to persistent ostream.
bool oneLoopNoBorn() const
Return true, if this matrix element includes one-loop corrections but no Born contributions.
void olpProcess(int pType, int id)
Set the process index, if this is an OLP handled matrix element.
vector< int > theOLPProcess
The process index, if this is an OLP handled matrix element.
ostream & operator<<(ostream &, const DecayIntegrator &)
Output information on the DecayIntegrator for debugging purposes.
virtual bool preInitialize() const
Return true if this object needs to be initialized before all other objects (except those for which t...
size_type size() const
Ptr< MatchboxPhasespace >::tptr phasespace() const
Return the phase space generator to be used.
ThePEG::Ptr< EventHandler >::transient_pointer tEHPtr
unsigned long oppositeSign
The number of points of opposite sign.
Process information with coupling order.
ThePEG::Ptr< CascadeHandler >::transient_pointer tCascHdlPtr
Provide easy access to MatchboxXComb XComb extensions.
virtual double renFac() const
Get the renormalization scale factor.
PersistentIStream & operator>>(PersistentIStream &is, map< ShowerInteraction, T, Cmp, A > &m)
Read a map with ShowerInteraction as the key.
bool noCorrelations() const
Return true, if this is a real emission matrix element which does not require colour correlators...
virtual unsigned int orderInAlphaS() const=0
virtual void constructVertex(tSubProPtr)
construct the spin information for the interaction
virtual bool wantCMS() const
Return true, if this matrix element expects the incoming partons in their center-of-mass system...
map< double, double > bins
The bins, indexed by upper bound.
virtual Energy2 scale() const
Return the scale associated with the phase space point provided by the last call to setKinematics()...
virtual bool hasFinalStateSymmetry() const
Return true, if this amplitude already includes symmetry factors for identical outgoing particles...
virtual bool needsOLPCorrelators() const
Return true, if colour and spin correlated matrix elements should be ordered from the OLP...
void doOneLoop()
Instruct this matrix element to include one-loop corrections.
virtual bool keepRandomNumbers() const
Return true, if the XComb steering this matrix element should keep track of the random numbers used t...
ThePEG::Ptr< SubProcess >::transient_pointer tSubProPtr
virtual bool canFillRhoMatrix() const
Return true, if this amplitude is capable of consistently filling the rho matrices for the spin corre...
bool checkedPDFs
True, if already checked for which PDFs to include.
bool theOneLoop
True, if this matrix element includes one-loop corrections.
Matchbox extensions to StandardXComb.
vector< cPDPtr > cPDVector
virtual double alphaS() const
Return the value of associated with the phase space point provided by the last call to setKinematics...
const vector< int > & olpProcess() const
Return the process index, if this is an OLP handled matrix element.
virtual void getDiagrams() const=0
int legsize() const
leg size
virtual bool hasInitialAverage() const
Return true, if this amplitude already includes averaging over incoming parton&#39;s quantum numbers...
const vector< Lorentz5Momentum > & lastMEMomenta() const
Return the meMomenta as generated at the last phase space point.
vector< Ptr< MatchboxReweightBase >::ptr > & reweights()
Access the reweights.
Ptr< MatchboxScaleChoice >::tptr scaleChoice() const
Return the scale choice object.
unsigned long underflow
The underflow.
vector< DiagPtr > DiagramVector
const vector< Ptr< MatchboxInsertionOperator >::ptr > & virtuals() const
Return the virtual corrections.
void needsNoCorrelations()
Indicate that this is a real emission matrix element which does not require colour correlators...
-*- C++ -*-
virtual bool hasRunningAlphaS() const
Indicate that this matrix element is running alphas by itself.
const vector< Ptr< MatchboxReweightBase >::ptr > & reweights() const
Return the reweights.
unsigned long overflow
The overflow.
Simple histogram for accuracy checks.
void doOneLoopNoBorn()
Instruct this matrix element to include one-loop corrections but no Born contributions.
pair< cPDPtr, cPDPtr > cPDPair
bool oneLoopNoLoops() const
Return true, if this matrix element includes one-loop corrections but no actual loop contributions...
vector< Ptr< MatchboxInsertionOperator >::ptr > & virtuals()
Return the virtual corrections.
virtual bool noMirror() const
Return true, if this matrix element does not want to make use of mirroring processes; in this case al...
virtual bool isOLPTree() const
Return true, if this matrix element is handled by a BLHA one-loop provider.
Ptr< MatchboxScaleChoice >::ptr theScaleChoice
The scale choice object.
vector< Ptr< MatchboxInsertionOperator >::ptr > theVirtuals
The virtual corrections.
void matchboxAmplitude(Ptr< MatchboxAmplitude >::ptr amp)
Set the amplitude.
virtual bool isOLPLoop() const
Return true, if this matrix element is handled by a BLHA one-loop provider.
SpinCorrelationTensor represents a spin correlation tensor of the form.
void scaleChoice(Ptr< MatchboxScaleChoice >::ptr sc)
Set the scale choice object.
constexpr ZeroUnit ZERO
virtual bool hasRunningAlphaEW() const
Indicate that this matrix element is running alphaew by itself.
const Process & subProcess() const
Return the subprocess.
ThePEG::Ptr< PartonExtractor >::transient_pointer tPExtrPtr