Page MenuHomeHEPForge
Files F7879015

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/General/MEvv2vs.cc b/MatrixElement/General/MEvv2vs.cc
--- a/MatrixElement/General/MEvv2vs.cc
+++ b/MatrixElement/General/MEvv2vs.cc
@@ -1,250 +1,259 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEvv2vs class.
//
#include "MEvv2vs.h"
+#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using ThePEG::Helicity::ScalarWaveFunction;
using ThePEG::Helicity::VectorWaveFunction;
using ThePEG::Helicity::incoming;
using ThePEG::Helicity::outgoing;
IBPtr MEvv2vs::clone() const {
return new_ptr(*this);
}
IBPtr MEvv2vs::fullclone() const {
return new_ptr(*this);
}
void MEvv2vs::persistentOutput(PersistentOStream & os) const {
- os << scalar_ << vector_;
+ os << scalar_ << vector_ << fourPointVertex_;
}
void MEvv2vs::persistentInput(PersistentIStream & is, int) {
- is >> scalar_ >> vector_;
+ is >> scalar_ >> vector_ >> fourPointVertex_;
initializeMatrixElements(PDT::Spin1, PDT::Spin1,
PDT::Spin1, PDT::Spin0);
}
-ClassDescription<MEvv2vs> MEvv2vs::initMEvv2vs;
-// Definition of the static class description member.
+// The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<MEvv2vs,GeneralHardME>
+describeHerwigMEvv2vs("Herwig::MEvv2vs", "Herwig.so");
void MEvv2vs::Init() {
static ClassDocumentation<MEvv2vs> documentation
("The MEvv2vs class implements the general matrix elements"
" for vector vector -> vector scalar");
}
void MEvv2vs::doinit() {
GeneralHardME::doinit();
scalar_.resize(numberOfDiags());
vector_.resize(numberOfDiags());
initializeMatrixElements(PDT::Spin1, PDT::Spin1,
PDT::Spin1, PDT::Spin0);
for(size_t i = 0; i < numberOfDiags(); ++i) {
HPDiagram diag = getProcessInfo()[i];
tcPDPtr offshell = diag.intermediate;
- assert(offshell);
- if(offshell->iSpin() == PDT::Spin0) {
+ if(!offshell) {
+ fourPointVertex_ =
+ dynamic_ptr_cast<AbstractVVVSVertexPtr>(diag.vertices.first);
+ }
+ else if(offshell->iSpin() == PDT::Spin0) {
AbstractVVSVertexPtr vert1;
AbstractVSSVertexPtr vert2;
if(diag.channelType == HPDiagram::sChannel ||
(diag.channelType == HPDiagram::tChannel && diag.ordered.second)) {
vert1 = dynamic_ptr_cast<AbstractVVSVertexPtr>(diag.vertices.first );
vert2 = dynamic_ptr_cast<AbstractVSSVertexPtr>(diag.vertices.second);
}
else {
vert1 = dynamic_ptr_cast<AbstractVVSVertexPtr>(diag.vertices.second);
vert2 = dynamic_ptr_cast<AbstractVSSVertexPtr>(diag.vertices.first );
}
scalar_[i] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1) {
AbstractVVVVertexPtr vert1;
AbstractVVSVertexPtr vert2;
if(diag.channelType == HPDiagram::sChannel ||
(diag.channelType == HPDiagram::tChannel && diag.ordered.second)) {
vert1 = dynamic_ptr_cast<AbstractVVVVertexPtr>(diag.vertices.first );
vert2 = dynamic_ptr_cast<AbstractVVSVertexPtr>(diag.vertices.second);
}
else {
vert1 = dynamic_ptr_cast<AbstractVVVVertexPtr>(diag.vertices.second);
vert2 = dynamic_ptr_cast<AbstractVVSVertexPtr>(diag.vertices.first );
}
vector_[i] = make_pair(vert1, vert2);
}
}
}
double MEvv2vs::me2() const {
VBVector va(2), vb(2), vc(3);
for(unsigned int i = 0; i < 2; ++i) {
va[i] = VectorWaveFunction(rescaledMomenta()[0], mePartonData()[0], 2*i,
- incoming);
+ incoming);
vb[i] = VectorWaveFunction(rescaledMomenta()[1], mePartonData()[1], 2*i,
- incoming);
+ incoming);
}
//always 0 and 2 polarisations
for(unsigned int i = 0; i < 2; ++i) {
vc[2*i] = VectorWaveFunction(rescaledMomenta()[2], mePartonData()[2], 2*i,
- outgoing);
+ outgoing);
}
bool mc = !(mePartonData()[2]->mass() > ZERO);
//massive vector, also 1
if( !mc )
vc[1] = VectorWaveFunction(rescaledMomenta()[2], mePartonData()[2], 1,
outgoing);
ScalarWaveFunction sd(rescaledMomenta()[3], mePartonData()[3], outgoing);
double full_me(0.);
vv2vsHeME(va, vb, vc, mc, sd, full_me,true);
return full_me;
}
ProductionMatrixElement
MEvv2vs::vv2vsHeME(VBVector & vin1, VBVector & vin2,
VBVector & vout1, bool mc, ScalarWaveFunction & sd,
double & me2, bool first) const {
const Energy2 q2(scale());
const Energy mass = vout1[0].mass();
// weights for the selection of the diagram
vector<double> me(numberOfDiags(), 0.);
// weights for the selection of the colour flow
vector<double> flow(numberOfFlows(),0.);
// flow over the helicities and diagrams
for(unsigned int ihel1 = 0; ihel1 < 2; ++ihel1) {
for(unsigned int ihel2 = 0; ihel2 < 2; ++ihel2) {
for(unsigned int ohel1 = 0; ohel1 < 3; ++ohel1) {
if(mc && ohel1 == 1) ++ohel1;
vector<Complex> flows(numberOfFlows(),0.);
for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
Complex diag(0.);
const HPDiagram & current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
- if(!offshell) continue;
- if(current.channelType == HPDiagram::sChannel) {
- if(offshell->iSpin() == PDT::Spin0) {
+ if(!offshell) {
+ diag = fourPointVertex_->evaluate(q2, vin1[ihel1], vin2[ihel2],
+ vout1[ohel1], sd);
+ }
+ else if(current.channelType == HPDiagram::sChannel) {
+ if(offshell->iSpin() == PDT::Spin0) {
+ ScalarWaveFunction interS = scalar_[ix].first->
+ evaluate(q2, 1, offshell,vin1[ihel1], vin2[ihel2]);
+ diag = scalar_[ix].second->
+ evaluate(q2, vout1[ohel1], sd, interS);
+ }
+ else if(offshell->iSpin() == PDT::Spin1) {
+ VectorWaveFunction interV = vector_[ix].first->
+ evaluate(q2, 1, offshell, vin1[ihel1], vin2[ihel2]);
+ diag = vector_[ix].second->
+ evaluate(q2, vout1[ohel1], interV, sd);
+ }
+ else
+ assert(false);
+ }
+ else if(current.channelType == HPDiagram::tChannel) {
+ if(offshell->iSpin() == PDT::Spin0) {
+ if(current.ordered.second) {
+ ScalarWaveFunction interS = scalar_[ix].
+ first->evaluate(q2, 3, offshell, vin1[ihel1],vout1[ohel1], mass);
+ diag = scalar_[ix].second->
+ evaluate(q2, vin2[ihel2], sd, interS);
+ }
+ else {
ScalarWaveFunction interS = scalar_[ix].first->
- evaluate(q2, 1, offshell,vin1[ihel1], vin2[ihel2]);
+ evaluate(q2, 3, offshell, vin2[ihel2],vout1[ohel1], mass);
diag = scalar_[ix].second->
- evaluate(q2, vout1[ohel1], sd, interS);
+ evaluate(q2, vin1[ihel1], sd, interS);
}
- else if(offshell->iSpin() == PDT::Spin1) {
+ }
+ else if(offshell->iSpin() == PDT::Spin1) {
+ if(current.ordered.second) {
+ VectorWaveFunction interV = vector_[ix].
+ first->evaluate(q2, 3, offshell, vin1[ihel1],vout1[ohel1], mass);
+ diag = vector_[ix].second->
+ evaluate(q2, vin2[ihel2], interV, sd);
+ }
+ else {
VectorWaveFunction interV = vector_[ix].first->
- evaluate(q2, 1, offshell, vin1[ihel1], vin2[ihel2]);
+ evaluate(q2, 3, offshell, vin2[ihel2],vout1[ohel1], mass);
diag = vector_[ix].second->
- evaluate(q2, vout1[ohel1], interV, sd);
+ evaluate(q2, vin1[ihel1], interV, sd);
}
- else
- assert(false);
}
- else if(current.channelType == HPDiagram::tChannel) {
- if(offshell->iSpin() == PDT::Spin0) {
- if(current.ordered.second) {
- ScalarWaveFunction interS = scalar_[ix].
- first->evaluate(q2, 3, offshell, vin1[ihel1],vout1[ohel1], mass);
- diag = scalar_[ix].second->
- evaluate(q2, vin2[ihel2], sd, interS);
- }
- else {
- ScalarWaveFunction interS = scalar_[ix].first->
- evaluate(q2, 3, offshell, vin2[ihel2],vout1[ohel1], mass);
- diag = scalar_[ix].second->
- evaluate(q2, vin1[ihel1], sd, interS);
- }
- }
- else if(offshell->iSpin() == PDT::Spin1) {
- if(current.ordered.second) {
- VectorWaveFunction interV = vector_[ix].
- first->evaluate(q2, 3, offshell, vin1[ihel1],vout1[ohel1], mass);
- diag = vector_[ix].second->
- evaluate(q2, vin2[ihel2], interV, sd);
- }
- else {
- VectorWaveFunction interV = vector_[ix].first->
- evaluate(q2, 3, offshell, vin2[ihel2],vout1[ohel1], mass);
- diag = vector_[ix].second->
- evaluate(q2, vin1[ihel1], interV, sd);
- }
- }
- else
- assert(false);
- }
- else
+ else
assert(false);
- me[ix] += norm(diag);
- diagramME()[ix](2*ihel1, 2*ihel2, ohel1,0) = diag;
- //Compute flows
- for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
- assert(current.colourFlow[iy].first<flows.size());
- flows[current.colourFlow[iy].first] +=
- current.colourFlow[iy].second * diag;
- }
+ }
+ else
+ assert(false);
+ me[ix] += norm(diag);
+ diagramME()[ix](2*ihel1, 2*ihel2, ohel1,0) = diag;
+ //Compute flows
+ for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
+ assert(current.colourFlow[iy].first<flows.size());
+ flows[current.colourFlow[iy].first] +=
+ current.colourFlow[iy].second * diag;
+ }
}
// MEs for the different colour flows
- for(unsigned int iy = 0; iy < numberOfFlows(); ++iy)
+ for(unsigned int iy = 0; iy < numberOfFlows(); ++iy)
flowME()[iy](2*ihel1, 2*ihel2, ohel1, 0) = flows[iy];
//Now add flows to me2 with appropriate colour factors
for(size_t ii = 0; ii < numberOfFlows(); ++ii)
for(size_t ij = 0; ij < numberOfFlows(); ++ij)
me2 += getColourFactors()[ii][ij]*(flows[ii]*conj(flows[ij])).real();
// contribution to the colour flow
for(unsigned int ii = 0; ii < numberOfFlows(); ++ii) {
flow[ii] += getColourFactors()[ii][ii]*norm(flows[ii]);
}
}
}
}
// if not computing the cross section return the selected colour flow
if(!first) return flowME()[colourFlow()];
me2 = selectColourFlow(flow,me,me2);
return flowME()[colourFlow()];
}
void MEvv2vs::constructVertex(tSubProPtr sub) {
ParticleVector ext = hardParticles(sub);
// set wave functions with real momenta
VBVector v1, v2, v3;
VectorWaveFunction(v1, ext[0], incoming, false, true);
VectorWaveFunction(v2, ext[1], incoming, false, true);
//function to calculate me2 expects massless incoming vectors
// and this constructor sets the '1' polarisation at element [2]
//in the vector
bool mc = !(ext[2]->data().mass() > ZERO);
VectorWaveFunction(v3, ext[2], outgoing, true, mc);
ScalarWaveFunction sd(ext[3], outgoing, true);
// Need to use rescale momenta to calculate matrix element
setRescaledMomenta(ext);
// wave functions with rescaled momenta
VectorWaveFunction vr1(rescaledMomenta()[0],
ext[0]->dataPtr(), incoming);
VectorWaveFunction vr2(rescaledMomenta()[1],
ext[1]->dataPtr(), incoming);
VectorWaveFunction vr3(rescaledMomenta()[2],
ext[2]->dataPtr(), outgoing);
ScalarWaveFunction sr4(rescaledMomenta()[3],
ext[3]->dataPtr(), outgoing);
for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
vr1.reset(2*ihel);
v1[ihel] = vr1;
vr2.reset(2*ihel);
v2[ihel] = vr2;
vr3.reset(2*ihel);
v3[2*ihel] = vr3;
}
if( !mc ) {
vr3.reset(1);
v3[1] = vr3;
}
double dummy(0.);
ProductionMatrixElement pme = vv2vsHeME(v1, v2, v3, mc, sd, dummy,false);
createVertex(pme,ext);
}
diff --git a/MatrixElement/General/MEvv2vs.h b/MatrixElement/General/MEvv2vs.h
--- a/MatrixElement/General/MEvv2vs.h
+++ b/MatrixElement/General/MEvv2vs.h
@@ -1,198 +1,163 @@
// -*- C++ -*-
#ifndef HERWIG_MEvv2vs_H
#define HERWIG_MEvv2vs_H
//
// This is the declaration of the MEvv2vs class.
//
#include "GeneralHardME.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractVVVSVertex.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::VectorWaveFunction;
using Helicity::ScalarWaveFunction;
/**
* This is the implementation of the matrix element for
* \f$2\to 2\f$ massless vector-boson pair to a vector and scalar boson.
* It inherits from GeneralHardME and implements the appropriate virtual
* member functions.
*
* @see \ref MEvv2vsInterfaces "The interfaces"
* defined for MEvv2vs.
*/
class MEvv2vs: public GeneralHardME {
public:
/**
* Typedef for VectorWaveFunction
*/
typedef vector<VectorWaveFunction> VBVector;
public:
/** @name Virtual functions required by the GeneralHardME class. */
//@{
/**
* The matrix element for the kinematical configuration
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
* @return the matrix element scaled with sHat() to give a
* dimensionless number.
*/
virtual double me2() const;
//@}
/**
* Construct the vertex information for the spin correlations
* @param sub Pointer to the relevent SubProcess
*/
virtual void constructVertex(tSubProPtr sub);
private:
/**
* Compute the matrix element for \f$V\, V\to V\, V\f$
* @param vin1 VectorWaveFunctions for first incoming particle
* @param vin2 VectorWaveFunctions for second incoming particle
* @param vout1 VectorWaveFunctions for first outgoing particle
* @param mc Whether vout1 is massless or not
* @param sout2 ScalarWaveFunction for outgoing particle
* @param me2 colour averaged, spin summed ME
* @param first Whether or not first call to decide if colour decomposition etc
* should be calculated
* @return ProductionMatrixElement containing results of
* helicity calculations
*/
ProductionMatrixElement
vv2vsHeME(VBVector & vin1, VBVector & vin2,
VBVector & vout1, bool mc, ScalarWaveFunction & sout2,
double & me2, bool first ) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
- * The static object used to initialize the description of this class.
- * Indicates that this is a concrete class with persistent data.
- */
- static ClassDescription<MEvv2vs> initMEvv2vs;
-
- /**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEvv2vs & operator=(const MEvv2vs &);
private:
/**
* Store the dynamically casted VVSVertex and VSSVertex pointers
*/
vector<pair<AbstractVVSVertexPtr, AbstractVSSVertexPtr> > scalar_;
/**
* Store the dynamically casted VVVVertex and VVSVertex pointers
*/
vector<pair<AbstractVVVVertexPtr, AbstractVVSVertexPtr> > vector_;
+ /**
+ * Store the dynamically casted VVVSVertex pointer
+ */
+ AbstractVVVSVertexPtr fourPointVertex_;
+
};
}
-#include "ThePEG/Utilities/ClassTraits.h"
-
-namespace ThePEG {
-
-/** @cond TRAITSPECIALIZATIONS */
-
-/** This template specialization informs ThePEG about the
- * base classes of MEvv2vs. */
-template <>
-struct BaseClassTrait<Herwig::MEvv2vs,1> {
- /** Typedef of the first base class of MEvv2vs. */
- typedef Herwig::GeneralHardME NthBase;
-};
-
-/** This template specialization informs ThePEG about the name of
- * the MEvv2vs class and the shared object where it is defined. */
-template <>
-struct ClassTraits<Herwig::MEvv2vs>
- : public ClassTraitsBase<Herwig::MEvv2vs> {
- /** Return a platform-independent class name */
- static string className() { return "Herwig::MEvv2vs"; }
- /**
- * The name of a file containing the dynamic library where the class
- * MEvv2vs is implemented. It may also include several, space-separated,
- * libraries if the class MEvv2vs depends on other classes (base classes
- * excepted). In this case the listed libraries will be dynamically
- * linked in the order they are specified.
- */
- static string library() { return "MEvv2vs.so"; }
-};
-
-/** @endcond */
-
-}
-
#endif /* HERWIG_MEvv2vs_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:15 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3790171
Default Alt Text
(17 KB)

Event Timeline