herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
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
26namespace Herwig {
27using namespace ThePEG;
28using Helicity::VertexBasePtr;
29
40class GeneralHardME: public HwMEBase {
41
42public:
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
74public:
75
80
81public:
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
145 diagrams(const DiagramVector & dv) const;
146
155 colourGeometries(tcDiagPtr diag) const;
157
169 void setProcessInfo(const vector<HPDiagram> & process,
171 unsigned int scaleOption,
172 double scaleFactor);
173
174public:
175
183
189 void persistentInput(PersistentIStream & is, int version);
191
198 static void Init();
199
200protected:
201
208 virtual void doinitrun();
210
211protected:
212
217 virtual void debug(double ) const {}
218
219protected:
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
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
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
384private:
385
391
392private:
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
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
493class MEException : public Exception {};
494
495}
496
497#endif /* HERWIG_GeneralHardME_H */
This defines the GeneralHardME class that is designed to serve as a base class for matrix elements of...
Definition: GeneralHardME.h:40
virtual Selector< DiagramIndex > diagrams(const DiagramVector &dv) const
Get diagram selector.
pair< long, long > incoming_
External particles.
virtual Selector< const ColourLines * > colourGeometries(tcDiagPtr diag) const
Return a Selector with possible colour geometries for the selected diagram weighted by their relative...
void diagram(unsigned int diag) const
Diagram.
void initializeMatrixElements(PDT::Spin in1, PDT::Spin in2, PDT::Spin out1, PDT::Spin out2)
Initialize the storage of the helicity matrix elements.
unsigned int diagram() const
Diagram.
vector< ProductionMatrixElement > & flowME() const
Access to the colour flow matrix element.
virtual void getDiagrams() const
Add all possible diagrams with the add() function.
virtual unsigned int orderInAlphaS() const
Return the order in in which this matrix element is given.
vector< ProductionMatrixElement > & diagramME() const
Access to the diagram matrix element.
virtual Energy2 scale() const
Return the scale associated with the last set phase space point.
GeneralHardME()
The default constructor.
void setProcessInfo(const vector< HPDiagram > &process, ColourStructure colour, bool debug, unsigned int scaleOption, double scaleFactor)
Set the diagrams and matrix of colour factors.
ColourStructure
Enum for the possible colour structures.
Definition: GeneralHardME.h:52
virtual void doinitrun()
Initialize this object.
GeneralHardME & operator=(const GeneralHardME &)=delete
The assignment operator is private and must never be called.
ParticleVector hardParticles(tSubProPtr subp)
Extract the paricles from the subprocess.
pair< long, long > getIncoming() const
Return the incoming pair.
static void Init()
The standard Init function used to initialize the interfaces.
double selectColourFlow(vector< double > &flow, vector< double > &me, double average) const
Calculate weight and select colour flow.
bool debug_
Whether to test the value of me2 against the analytical function.
vector< DVector > colour_
Store colour factors for ME calc.
vector< ProductionMatrixElement > diagramME_
Matrix elements for the different Feynman diagrams.
virtual double me2() const =0
The matrix element for the kinematical configuration previously provided by the last call to setKinem...
unsigned int colourFlow() const
Set/Get Info on the selected diagram and colour flow.
const vector< HPDiagram > & getProcessInfo() const
Access the HPDiagrams that store the required information to create the diagrams.
unsigned int numberOfFlows_
The number of colourflows.
HPCount numberOfDiags() const
Get the number of diagrams in this process.
ColourStructure colour() const
Access to the colour structure.
vector< ProductionMatrixElement > flowME_
Storage of the matrix elements.
const vector< DVector > & getColourFactors() const
Return the matrix of colour factors.
HPCount numberOfDiagrams_
Store the number of diagrams for fast retrieval.
size_t numberOfFlows() const
Access number of colour flows.
unsigned int diagram_
Diagram.
void colourFlow(unsigned int flow) const
Colour flow.
unsigned int scaleChoice_
The scale chocie.
virtual unsigned int orderInAlphaEW() const
Return the order in in which this matrix element is given.
double scaleFactor_
The scale factor.
vector< HPDiagram >::size_type HPCount
Convenient typedef for size_type of HPDiagram vector.
Definition: GeneralHardME.h:47
void persistentInput(PersistentIStream &is, int version)
Function used to read in object persistently.
pair< long, long > outgoing_
Store the outgoing particles.
void persistentOutput(PersistentOStream &os) const
Function used to write out object persistently.
bool debugME() const
Whether to print the debug information.
unsigned int flow_
Info on the selected diagram and colour flow.
void setRescaledMomenta(const ParticleVector &external)
Set the rescaled momenta.
virtual void debug(double) const
A debugging function to test the value of me2 against an analytic function.
void createVertex(ProductionMatrixElement &me, ParticleVector &external)
Create the vertes.
vector< HPDiagram > diagrams_
Diagrams.
pair< long, long > getOutgoing() const
Return the outgoing pair.
ColourStructure colourStructure_
Colour information.
The HwMEBase class serves a number of purposes.
Definition: HwMEBase.h:35
Energy2 uHat() const
Return the of the last set phase space point.
Definition: HwMEBase.h:176
Energy2 tHat() const
Return the of the last set phase space point.
Definition: HwMEBase.h:171
bool rescaleMomenta(const vector< Lorentz5Momentum > &, const cPDVector &)
rescale the momenta for the computation of the helicity matrix element
Exception class to indicate a problem has occurred with setting up to matrix element.
The storage of the helicity amplitude expression for the matrix element of a hard process.
vector< Lorentz5Momentum > & meMomenta()
vector< DiagPtr > DiagramVector
Energy2 sHat() const
-*- C++ -*-
ThePEG::Ptr< SubProcess >::transient_pointer tSubProPtr
vector< PPtr > ParticleVector
constexpr ZeroUnit ZERO
vector< cPDPtr > cPDVector
ThePEG::Ptr< SpinInfo >::transient_pointer tSpinPtr
void debug() const