herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
GeneralHardME.h
1 // -*- C++ -*-
2 //
3 // GeneralHardME.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_GeneralHardME_H
10 #define HERWIG_GeneralHardME_H
11 //
12 // This is the declaration of the GeneralHardME class.
13 //
14 
15 #include "Herwig/MatrixElement/HwMEBase.h"
16 #include "ThePEG/Utilities/Exception.h"
17 #include "ThePEG/Persistency/PersistentOStream.h"
18 #include "ThePEG/Persistency/PersistentIStream.h"
19 #include "Herwig/Models/General/HPDiagram.h"
20 #include "Herwig/MatrixElement/ProductionMatrixElement.h"
21 #include "Herwig/MatrixElement/HardVertex.h"
22 #include "ThePEG/EventRecord/SpinInfo.h"
23 #include "ThePEG/PDF/PolarizedBeamParticleData.h"
24 #include "GeneralHardME.fh"
25 
26 namespace Herwig {
27 using namespace ThePEG;
28 using Helicity::VertexBasePtr;
29 
40 class GeneralHardME: public HwMEBase {
41 
42 public:
43 
47  typedef vector<HPDiagram>::size_type HPCount;
48 
52  enum ColourStructure {UNDEFINED,
53  Colour11to11,Colour11to33bar,Colour11to88,
54  Colour33to33,Colour33barto11,Colour33barto33bar,
55  Colour33barto66bar, Colour33barto6bar6,
56  Colour33to61, Colour3bar3barto6bar1,
57  Colour33to16, Colour3bar3barto16bar,
58  Colour38to3bar6, Colour38to63bar,
59  Colour33barto18,Colour33barto81,Colour33barto88,
60  Colour38to13,Colour38to31,
61  Colour38to83,Colour38to38,
62  Colour3bar3barto3bar3bar,
63  Colour3bar8to13bar,Colour3bar8to3bar1,
64  Colour3bar8to83bar,Colour3bar8to3bar8,
65  Colour88to11,Colour88to33bar,
66  Colour88to66bar,Colour88to88,
67  Colour88to18,Colour88to81,
68  Colour33to13bar,Colour33to3bar1,
69  Colour33to83bar,Colour33to3bar8,
70  Colour3bar3barto13,Colour3bar3barto31,
71  Colour3bar3barto83,Colour3bar3barto38,
72  Colour38to3bar3bar,Colour3bar8to33};
73 
74 public:
75 
79  GeneralHardME();
80 
81 public:
82 
89  virtual unsigned int orderInAlphaS() const;
90 
95  virtual unsigned int orderInAlphaEW() const;
96 
104  virtual double me2() const = 0;
105 
109  virtual Energy2 scale() const {
110  if(scaleChoice_==0) {
111  return scaleFactor_*sHat();
112  }
113  else if(scaleChoice_==1) {
114  Energy2 mbar = 0.5*(meMomenta()[2].mass2()+meMomenta()[3].mass2());
115  Energy2 t = 0.5*(tHat()-mbar);
116  Energy2 u = 0.5*(uHat()-mbar);
117  Energy2 s = 0.5*sHat();
118  return scaleFactor_*4.*s*t*u/(s*s+t*t+u*u);
119  }
120  else if(scaleChoice_ ==2) {
121  Energy2 scale1 = meMomenta()[2].mass2()+meMomenta()[2].perp2();
122  Energy2 scale2 = meMomenta()[3].mass2()+meMomenta()[3].perp2();
123  return scaleFactor_*max(scale1,scale2);
124  }
125  else {
126  assert(false);
127  return ZERO;
128  }
129  }
130 
134  virtual void getDiagrams() const;
135 
144  virtual Selector<DiagramIndex>
145  diagrams(const DiagramVector & dv) const;
146 
155  colourGeometries(tcDiagPtr diag) const;
157 
169  void setProcessInfo(const vector<HPDiagram> & process,
170  ColourStructure colour, bool debug,
171  unsigned int scaleOption,
172  double scaleFactor);
173 
174 public:
175 
182  void persistentOutput(PersistentOStream & os) const;
183 
189  void persistentInput(PersistentIStream & is, int version);
191 
198  static void Init();
199 
200 protected:
201 
208  virtual void doinitrun();
210 
211 protected:
212 
217  virtual void debug(double ) const {}
218 
219 protected:
220 
225  const vector<HPDiagram> & getProcessInfo() const {
226  return diagrams_;
227  }
228 
233  pair<long, long> getIncoming() const {
234  return incoming_;
235  }
236 
241  pair<long, long> getOutgoing() const {
242  return outgoing_;
243  }
244 
248  const vector<DVector> & getColourFactors() const {
249  return colour_;
250  }
251 
255  HPCount numberOfDiags() const {
256  return numberOfDiagrams_;
257  }
258 
262  size_t numberOfFlows() const {
263  return numberOfFlows_;
264  }
265 
269  bool debugME() const {
270  return debug_;
271  }
272 
280  unsigned int colourFlow() const {return flow_;}
281 
285  void colourFlow(unsigned int flow) const {flow_=flow;}
286 
290  unsigned int diagram() const {return diagram_;}
291 
295  void diagram(unsigned int diag) const {diagram_=diag;}
297 
301  double selectColourFlow(vector<double> & flow,
302  vector<double> & me,double average) const;
303 
307  vector<ProductionMatrixElement> & flowME() const {
308  return flowME_;
309  }
310 
314  vector<ProductionMatrixElement> & diagramME() const {
315  return diagramME_;
316  }
317 
321  ColourStructure colour() const {return colourStructure_;}
322 
327  ParticleVector output(4);
328  output[0] = subp->incoming().first;
329  output[1] = subp->incoming().second;
330  output[2] = subp->outgoing()[0];
331  output[3] = subp->outgoing()[1];
332  //ensure particle ordering is the same as it was when
333  //the diagrams were created
334  if( output[0]->id() != getIncoming().first )
335  swap(output[0], output[1]);
336  if( output[2]->id() != getOutgoing().first )
337  swap(output[2], output[3]);
338  // return answer
339  return output;
340  }
341 
345  void setRescaledMomenta(const ParticleVector & external) {
346  cPDVector data(4);
347  vector<Lorentz5Momentum> momenta(4);
348  for( size_t i = 0; i < 4; ++i ) {
349  data[i] = external[i]->dataPtr();
350  momenta[i] = external[i]->momentum();
351  }
352  rescaleMomenta(momenta, data);
353  }
354 
359  ParticleVector & external) {
360  HardVertexPtr hardvertex = new_ptr(HardVertex());
361  hardvertex->ME(me);
362  for(ParticleVector::size_type i = 0; i < 4; ++i) {
363  tSpinPtr spin = external[i]->spinInfo();
364  if(i<2) {
365  tcPolarizedBeamPDPtr beam =
366  dynamic_ptr_cast<tcPolarizedBeamPDPtr>(external[i]->dataPtr());
367  if(beam) spin->rhoMatrix() = beam->rhoMatrix();
368  }
369  spin->productionVertex(hardvertex);
370  }
371  }
372 
377  PDT::Spin out1, PDT::Spin out2) {
378  flowME().resize(numberOfFlows(),
379  ProductionMatrixElement(in1,in2,out1,out2));
380  diagramME().resize(numberOfDiags(),
381  ProductionMatrixElement(in1,in2,out1,out2));
382  }
383 
384 private:
385 
390  GeneralHardME & operator=(const GeneralHardME &) = delete;
391 
392 private:
393 
401  pair<long, long> incoming_;
402 
406  pair<long, long> outgoing_;
408 
416  vector<HPDiagram> diagrams_;
417 
423 
432 
436  vector<DVector> colour_;
437 
441  unsigned int numberOfFlows_;
443 
447  bool debug_;
448 
452  unsigned int scaleChoice_;
453 
457  double scaleFactor_;
458 
466  mutable unsigned int flow_;
467 
471  mutable unsigned int diagram_;
473 
481  mutable vector<ProductionMatrixElement> flowME_;
482 
486  mutable vector<ProductionMatrixElement> diagramME_;
488 
489 };
490 
493 class MEException : public Exception {};
494 
495 }
496 
497 #endif /* HERWIG_GeneralHardME_H */
unsigned int colourFlow() const
Set/Get Info on the selected diagram and colour flow.
ColourStructure
Enum for the possible colour structures.
Definition: GeneralHardME.h:52
This defines the GeneralHardME class that is designed to serve as a base class for matrix elements of...
Definition: GeneralHardME.h:40
ColourStructure colour() const
Access to the colour structure.
unsigned int flow_
Info on the selected diagram and colour flow.
void diagram(unsigned int diag) const
Diagram.
ThePEG::Ptr< SpinInfo >::transient_pointer tSpinPtr
void initializeMatrixElements(PDT::Spin in1, PDT::Spin in2, PDT::Spin out1, PDT::Spin out2)
Initialize the storage of the helicity matrix elements.
virtual Energy2 scale() const
Return the scale associated with the last set phase space point.
bool debug_
Whether to test the value of me2 against the analytical function.
vector< HPDiagram > diagrams_
Diagrams.
bool debugME() const
Whether to print the debug information.
size_t numberOfFlows() const
Access number of colour flows.
pair< long, long > outgoing_
Store the outgoing particles.
vector< ProductionMatrixElement > & flowME() const
Access to the colour flow matrix element.
vector< DVector > colour_
Store colour factors for ME calc.
vector< ProductionMatrixElement > diagramME_
Matrix elements for the different Feynman diagrams.
HPCount numberOfDiags() const
Get the number of diagrams in this process.
const vector< HPDiagram > & getProcessInfo() const
Access the HPDiagrams that store the required information to create the diagrams. ...
ParticleVector hardParticles(tSubProPtr subp)
Extract the paricles from the subprocess.
ColourStructure colourStructure_
Colour information.
The storage of the helicity amplitude expression for the matrix element of a hard process...
pair< long, long > getOutgoing() const
Return the outgoing pair.
double scaleFactor_
The scale factor.
pair< long, long > incoming_
External particles.
vector< ProductionMatrixElement > flowME_
Storage of the matrix elements.
pair< long, long > getIncoming() const
Return the incoming pair.
void createVertex(ProductionMatrixElement &me, ParticleVector &external)
Create the vertes.
ThePEG::Ptr< SubProcess >::transient_pointer tSubProPtr
void setRescaledMomenta(const ParticleVector &external)
Set the rescaled momenta.
const vector< DVector > & getColourFactors() const
Return the matrix of colour factors.
HPCount numberOfDiagrams_
Store the number of diagrams for fast retrieval.
Exception class to indicate a problem has occurred with setting up to matrix element.
vector< cPDPtr > cPDVector
unsigned int scaleChoice_
The scale chocie.
vector< DiagPtr > DiagramVector
unsigned int diagram() const
Diagram.
-*- C++ -*-
unsigned int diagram_
Diagram.
void colourFlow(unsigned int flow) const
Colour flow.
vector< HPDiagram >::size_type HPCount
Convenient typedef for size_type of HPDiagram vector.
Definition: GeneralHardME.h:47
The HwMEBase class serves a number of purposes.
Definition: HwMEBase.h:35
vector< PPtr > ParticleVector
vector< ProductionMatrixElement > & diagramME() const
Access to the diagram matrix element.
virtual void debug(double) const
A debugging function to test the value of me2 against an analytic function.
constexpr ZeroUnit ZERO
unsigned int numberOfFlows_
The number of colourflows.