Page MenuHomeHEPForge
Files F9501529

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/General/MEvv2ff.h b/MatrixElement/General/MEvv2ff.h
--- a/MatrixElement/General/MEvv2ff.h
+++ b/MatrixElement/General/MEvv2ff.h
@@ -1,193 +1,191 @@
// -*- C++ -*-
//
// MEvv2ff.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MEvv2ff_H
#define HERWIG_MEvv2ff_H
//
// This is the declaration of the MEvv2ff class.
//
#include "GeneralHardME.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVVertex.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
using ThePEG::Helicity::SpinorWaveFunction;
using ThePEG::Helicity::SpinorBarWaveFunction;
using ThePEG::Helicity::VectorWaveFunction;
/**
* This class is designed to implement the matrix element for the
* \f$2 \rightarrow 2\f$ process vector-vector to fermion-antifermion pair. It
* inherits from GeneralHardME and implements the me2() virtual function.
*
- * @see \ref MEvv2ffInterfaces "The Interfaces"
- * defined for MEvv2ff.
* @see GeneralHardME
*
*/
class MEvv2ff: public GeneralHardME {
public:
/** A Vector of VectorWaveFunction objects. */
typedef vector<VectorWaveFunction> VBVector;
/** A vector of SpinorBarWaveFunction objects. */
typedef vector<SpinorWaveFunction> SpinorVector;
/** A vector of SpinorBarWaveFunction objects. */
typedef vector<SpinorBarWaveFunction> SpinorBarVector;
public:
/** @name Virtual functions required by the MEBase 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:
/**
* Calculate the value of the matrix element
*/
ProductionMatrixElement vv2ffME(const VBVector & v1, const VBVector & v2,
const SpinorBarVector & sbar,
const SpinorVector & sp,
double & me2, bool first) const;
protected:
/**
* A debugging function to test the value of me2 against an
* analytic function.
* @param me2 The value of the \f$ |\bar{\mathcal{M}}|^2 \f$
*/
virtual void debug(double me2) 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 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();
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** 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 {return new_ptr(*this);}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEvv2ff & operator=(const MEvv2ff &);
private:
/** @name Dynamically casted vertices. */
//@{
/**
* Intermediate scalar
*/
vector<pair<AbstractVVSVertexPtr, AbstractFFSVertexPtr > > scalar_;
/**
* Intermediate fermion
*/
vector<pair<AbstractFFVVertexPtr, AbstractFFVVertexPtr> > fermion_;
/**
* Intermediate vector
*/
vector<pair<AbstractVVVVertexPtr, AbstractFFVVertexPtr> > vector_;
/**
* Intermediate tensor
*/
vector<pair<AbstractVVTVertexPtr, AbstractFFTVertexPtr> > tensor_;
/**
* Four point vertices
*/
vector<AbstractFFVVVertexPtr> four_;
//@}
};
}
#endif /* HERWIG_MEvv2ff_H */
diff --git a/MatrixElement/General/MEvv2rf.cc b/MatrixElement/General/MEvv2rf.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/General/MEvv2rf.cc
@@ -0,0 +1,377 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the MEvv2rf class.
+//
+
+#include "MEvv2rf.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+using ThePEG::Helicity::incoming;
+using ThePEG::Helicity::outgoing;
+
+IBPtr MEvv2rf::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr MEvv2rf::fullclone() const {
+ return new_ptr(*this);
+}
+
+void MEvv2rf::doinit() {
+ GeneralHardME::doinit();
+ scalar_ .resize(numberOfDiags());
+ fermion_.resize(numberOfDiags());
+ vector_ .resize(numberOfDiags());
+ four_ .resize(numberOfDiags());
+ initializeMatrixElements(PDT::Spin1 , PDT::Spin1,
+ PDT::Spin3Half, PDT::Spin1Half);
+
+ for( size_t i = 0; i < numberOfDiags(); ++i ) {
+ HPDiagram dg = getProcessInfo()[i];
+ if( dg.channelType == HPDiagram::tChannel ) {
+ AbstractRFVVertexPtr rfv;
+ AbstractFFVVertexPtr ffv;
+ if(dg.ordered.second) {
+ rfv = dynamic_ptr_cast<AbstractRFVVertexPtr>(dg.vertices.first);
+ ffv = dynamic_ptr_cast<AbstractFFVVertexPtr>(dg.vertices.second);
+ }
+ else {
+ rfv = dynamic_ptr_cast<AbstractRFVVertexPtr>(dg.vertices.second);
+ ffv = dynamic_ptr_cast<AbstractFFVVertexPtr>(dg.vertices.first );
+ }
+ fermion_[i] = make_pair(rfv, ffv);
+ }
+ else if( dg.channelType == HPDiagram::sChannel ) {
+ if( dg.intermediate->iSpin() == PDT::Spin0 ) {
+ AbstractVVSVertexPtr vvs =
+ dynamic_ptr_cast<AbstractVVSVertexPtr>(dg.vertices.first );
+ AbstractRFSVertexPtr rfs =
+ dynamic_ptr_cast<AbstractRFSVertexPtr>(dg.vertices.second);
+ scalar_[i] = make_pair(vvs,rfs);
+ }
+ else if( dg.intermediate->iSpin() == PDT::Spin1) {
+ AbstractVVVVertexPtr vvv =
+ dynamic_ptr_cast<AbstractVVVVertexPtr>(dg.vertices.first);
+ AbstractRFVVertexPtr rfv =
+ dynamic_ptr_cast<AbstractRFVVertexPtr>(dg.vertices.second);
+ vector_[i] = make_pair(vvv,rfv);
+ }
+ }
+ else if ( dg.channelType == HPDiagram::fourPoint) {
+ four_[i] = dynamic_ptr_cast<AbstractRFVVVertexPtr>(dg.vertices.first);
+ }
+ }
+}
+
+void MEvv2rf::persistentOutput(PersistentOStream & os) const {
+ os << scalar_ << fermion_ << vector_ << four_;
+}
+
+void MEvv2rf::persistentInput(PersistentIStream & is, int) {
+ is >> scalar_ >> fermion_ >> vector_ >> four_;
+ initializeMatrixElements(PDT::Spin1 , PDT::Spin1,
+ PDT::Spin3Half, PDT::Spin1Half);
+}
+
+//The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<MEvv2rf,GeneralHardME>
+describeHerwigMEvv2rf("Herwig::MEvv2rf", "Herwig.so");
+
+void MEvv2rf::Init() {
+
+ static ClassDocumentation<MEvv2rf> documentation
+ ("The MEvv2rf class handes the ME calculation for vv -> rf");
+
+}
+
+double MEvv2rf::me2() const {
+ // set up the vector wavefunctions
+ VBVector v1(2), v2(2);
+ for( size_t i = 0; i < 2; ++i ) {
+ v1[i] = VectorWaveFunction(rescaledMomenta()[0],mePartonData()[0], 2*i,
+ incoming);
+ v2[i] = VectorWaveFunction(rescaledMomenta()[1],mePartonData()[1], 2*i,
+ incoming);
+ }
+ // setup spinor wavefunctions and decide which case to use
+ bool massless = mePartonData()[2]->mass()==ZERO;
+ double full_me(0.);
+ if(mePartonData()[2]->id()<0) {
+ RSSpinorVector sp(4); SpinorBarVector sbar(2);
+ for( size_t i = 0; i < 4; ++i ) {
+ if(massless && (i==2||i==3)) continue;
+ sp[i] = RSSpinorWaveFunction(rescaledMomenta()[2], mePartonData()[2], i,
+ outgoing);
+ }
+ for( size_t i = 0; i < 2; ++i ) {
+ sbar[i] = SpinorBarWaveFunction(rescaledMomenta()[3], mePartonData()[3], i,
+ outgoing);
+ }
+ vv2frME(v1, v2, sbar, sp, full_me,true);
+ }
+ else {
+ SpinorVector sp(2); RSSpinorBarVector sbar(4);
+ for( size_t i = 0; i < 4; ++i ) {
+ if(massless && (i==2||i==3)) continue;
+ sbar[i] = RSSpinorBarWaveFunction(rescaledMomenta()[2], mePartonData()[2], i,
+ outgoing);
+ }
+ for( size_t i = 0; i < 2; ++i ) {
+ sp[i] = SpinorWaveFunction(rescaledMomenta()[3], mePartonData()[3], i,
+ outgoing);
+ }
+ vv2rfME(v1, v2, sbar, sp, full_me,true);
+ }
+ return full_me;
+}
+
+ProductionMatrixElement
+MEvv2rf::vv2rfME(const VBVector & v1, const VBVector & v2,
+ const RSSpinorBarVector & sbar,const SpinorVector & sp,
+ double & me2, bool first) const {
+ // scale
+ const Energy2 q2 = scale();
+ // whether or not rs fermion is massless
+ bool massless = mePartonData()[2]->mass()==ZERO;
+ // 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.);
+ //sum over vector helicities
+ for(unsigned int iv1 = 0; iv1 < 2; ++iv1) {
+ for(unsigned int iv2 = 0; iv2 < 2; ++iv2) {
+ //sum over fermion helicities
+ for(unsigned int of1 = 0; of1 < 4; ++of1) {
+ if(massless && (of1==1 || of1==2) ) continue;
+ for(unsigned int of2 = 0; of2 < 2; ++of2) {
+ vector<Complex> flows(numberOfFlows(),0.);
+ for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
+ Complex diag(0.);
+ const HPDiagram & current = getProcessInfo()[ix];
+ tPDPtr offshell = current.intermediate;
+ if(current.channelType == HPDiagram::tChannel &&
+ offshell->iSpin() == PDT::Spin1Half) {
+ if(current.ordered.second) {
+ SpinorBarWaveFunction interF = fermion_[ix].first->
+ evaluate(q2, 3, offshell, sbar[of1], v1[iv1]);
+ diag = fermion_[ix].second->
+ evaluate(q2, sp[of2], interF, v2[iv2]);
+ }
+ else {
+ SpinorWaveFunction interF = fermion_[ix].second->
+ evaluate(q2, 3, offshell, sp[of2], v1[iv1]);
+ diag = fermion_[ix].first->
+ evaluate(q2, interF, sbar[of1], v2[iv2]);
+ }
+ }
+ else if(current.channelType == HPDiagram::sChannel) {
+ if(offshell->iSpin() == PDT::Spin0) {
+ ScalarWaveFunction interS = scalar_[ix].first->
+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
+ diag = scalar_[ix].second->
+ evaluate(q2, sp[of2], sbar[of1], interS);
+ }
+ else if(offshell->iSpin() == PDT::Spin1) {
+ VectorWaveFunction interV = vector_[ix].first->
+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
+ diag = vector_[ix].second->
+ evaluate(q2, sp[of2], sbar[of1], interV);
+ }
+ }
+ else if(current.channelType == HPDiagram::fourPoint) {
+ diag = four_[ix]->evaluate(q2, sp[of2], sbar[of1], v1[iv1], v2[iv2]);
+ }
+ else
+ assert(false);
+ me[ix] += norm(diag);
+ diagramME()[ix](2*iv1, 2*iv2, of1, of2) = 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)
+ flowME()[iy](2*iv1, 2*iv2, of1, of2) = 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()];
+}
+
+ProductionMatrixElement
+MEvv2rf::vv2frME(const VBVector & v1, const VBVector & v2,
+ const SpinorBarVector & sbar,const RSSpinorVector & sp,
+ double & me2, bool first) const {
+ // scale
+ const Energy2 q2 = scale();
+ // whether or not rs fermion is massless
+ bool massless = mePartonData()[2]->mass()==ZERO;
+ // 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.);
+ //sum over vector helicities
+ for(unsigned int iv1 = 0; iv1 < 2; ++iv1) {
+ for(unsigned int iv2 = 0; iv2 < 2; ++iv2) {
+ //sum over fermion helicities
+ for(unsigned int of1 = 0; of1 < 4; ++of1) {
+ for(unsigned int of2 = 0; of2 < 2; ++of2) {
+ if(massless && (of2==1 || of2==2) ) continue;
+ vector<Complex> flows(numberOfFlows(),0.);
+ for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
+ Complex diag(0.);
+ const HPDiagram & current = getProcessInfo()[ix];
+ tPDPtr offshell = current.intermediate;
+ if(current.channelType == HPDiagram::tChannel &&
+ offshell->iSpin() == PDT::Spin1Half) {
+ if(current.ordered.second) {
+ SpinorWaveFunction interF = fermion_[ix].first->
+ evaluate(q2, 3, offshell, sp[of2], v1[iv1]);
+ diag = fermion_[ix].second->
+ evaluate(q2, interF, sbar[of1],v2[iv2]);
+ }
+ else {
+ SpinorBarWaveFunction interF = fermion_[ix].second->
+ evaluate(q2, 3, offshell, sbar[of1], v1[iv1]);
+ diag = fermion_[ix].first->
+ evaluate(q2, sp[of2], interF, v2[iv2]);
+ }
+ }
+ else if(current.channelType == HPDiagram::sChannel) {
+ if(offshell->iSpin() == PDT::Spin0) {
+ ScalarWaveFunction interS = scalar_[ix].first->
+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
+ diag = scalar_[ix].second->
+ evaluate(q2, sp[of2], sbar[of1], interS);
+ }
+ else if(offshell->iSpin() == PDT::Spin1) {
+ VectorWaveFunction interV = vector_[ix].first->
+ evaluate(q2, 1, offshell, v1[iv1], v2[iv2]);
+ diag = vector_[ix].second->
+ evaluate(q2, sp[of2], sbar[of1], interV);
+ }
+ }
+ else if(current.channelType == HPDiagram::fourPoint) {
+ diag = four_[ix]->evaluate(q2, sp[of2], sbar[of1], v1[iv1], v2[iv2]);
+ }
+ else
+ assert(false);
+ me[ix] += norm(diag);
+ diagramME()[ix](2*iv1, 2*iv2, of2, of1) = 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)
+ flowME()[iy](2*iv1, 2*iv2, of1, of2) = 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 MEvv2rf::constructVertex(tSubProPtr sub) {
+ ParticleVector ext = hardParticles(sub);
+ // vector wavefunctions are common do them first
+ // wavefunction with real momenta
+ VBVector v1, v2;
+ VectorWaveFunction(v1, ext[0], incoming, false, true);
+ VectorWaveFunction(v2, ext[1], incoming, false, true);
+ // rescale momenta
+ setRescaledMomenta(ext);
+ // wavefunctions with rescaled momenta
+ VectorWaveFunction v1r(rescaledMomenta()[0],
+ ext[0]->dataPtr(), incoming);
+ VectorWaveFunction v2r(rescaledMomenta()[1],
+ ext[1]->dataPtr(), incoming);
+ ProductionMatrixElement pme;
+ for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
+ v1r.reset(2*ihel);
+ v1[ihel] = v1r;
+ v2r.reset(2*ihel);
+ v2[ihel] = v2r;
+ }
+ // setup spinor wavefunctions and decide which case to use
+ bool massless = mePartonData()[2]->mass()==ZERO;
+ double dummy(0.);
+ if(ext[2]->id()<0) {
+ RSSpinorVector sp;
+ RSSpinorWaveFunction(sp, ext[2], outgoing, true);
+ SpinorBarVector sbar;
+ SpinorBarWaveFunction(sbar, ext[3], outgoing, true);
+ // wavefuncions with rescaled momenta
+ SpinorBarWaveFunction sbr(rescaledMomenta()[3],
+ ext[3]->dataPtr(), outgoing);
+ RSSpinorWaveFunction spr(rescaledMomenta()[2],
+ ext[2]->dataPtr(), outgoing);
+ for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
+ sbr.reset(ihel);
+ sbar[ihel] = sbr;
+ }
+ for( unsigned int ihel = 0; ihel < 4; ++ihel ) {
+ if(massless && (ihel==1 || ihel==2)) continue;
+ spr.reset(ihel);
+ sp[ihel] = spr;
+ }
+ pme = vv2frME(v1, v2, sbar, sp, dummy,false);
+ }
+ else {
+ SpinorVector sp;
+ SpinorWaveFunction(sp, ext[3], outgoing, true);
+ RSSpinorBarVector sbar;
+ RSSpinorBarWaveFunction(sbar, ext[2], outgoing, true);
+ // wavefuncions with rescaled momenta
+ RSSpinorBarWaveFunction sbr(rescaledMomenta()[2],
+ ext[2]->dataPtr(), outgoing);
+ SpinorWaveFunction spr(rescaledMomenta()[3],
+ ext[3]->dataPtr(), outgoing);
+ for( unsigned int ihel = 0; ihel < 4; ++ihel ) {
+ if(massless && (ihel==1 || ihel==2)) continue;
+ sbr.reset(ihel);
+ sbar[ihel] = sbr;
+ }
+ for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
+ spr.reset(ihel);
+ sp[ihel] = spr;
+ }
+ pme = vv2rfME(v1, v2, sbar, sp, dummy,false);
+ }
+ createVertex(pme,ext);
+}
diff --git a/MatrixElement/General/MEvv2rf.h b/MatrixElement/General/MEvv2rf.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/General/MEvv2rf.h
@@ -0,0 +1,182 @@
+// -*- C++ -*-
+#ifndef Herwig_MEvv2rf_H
+#define Herwig_MEvv2rf_H
+//
+// This is the declaration of the MEvv2rf class.
+//
+
+#include "GeneralHardME.h"
+#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractRFVVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractRFVVVertex.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/RSSpinorBarWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
+#include "Herwig/MatrixElement/ProductionMatrixElement.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * This class is designed to implement the matrix element for the
+ * \f$2 \rightarrow 2\f$ process vector-vector to fermion- RS fermion. It
+ * inherits from GeneralHardME and implements the me2() virtual function.
+ *
+ * @see GeneralHardME
+ *
+ */
+class MEvv2rf: public GeneralHardME {
+
+public:
+
+ /** A Vector of VectorWaveFunction objects. */
+ typedef vector<VectorWaveFunction> VBVector;
+
+ /** A vector of SpinorBarWaveFunction objects. */
+ typedef vector<SpinorWaveFunction> SpinorVector;
+
+ /** A vector of SpinorBarWaveFunction objects. */
+ typedef vector<SpinorBarWaveFunction> SpinorBarVector;
+
+ /** A vector of SpinorBarWaveFunction objects. */
+ typedef vector<RSSpinorWaveFunction> RSSpinorVector;
+
+ /** A vector of SpinorBarWaveFunction objects. */
+ typedef vector<RSSpinorBarWaveFunction> RSSpinorBarVector;
+
+public:
+
+ /** @name Virtual functions required by the MEBase 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);
+
+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:
+
+ /**
+ * Calculate the value of the matrix element
+ */
+ ProductionMatrixElement vv2rfME(const VBVector & v1, const VBVector & v2,
+ const RSSpinorBarVector & sbar,
+ const SpinorVector & sp,
+ double & me2, bool first) const;
+
+ /**
+ * Calculate the value of the matrix element
+ */
+ ProductionMatrixElement vv2frME(const VBVector & v1, const VBVector & v2,
+ const SpinorBarVector & sbar,
+ const RSSpinorVector & sp,
+ double & me2, bool first) const;
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ MEvv2rf & operator=(const MEvv2rf &);
+private:
+
+ /** @name Dynamically casted vertices. */
+ //@{
+ /**
+ * Intermediate scalar
+ */
+ vector<pair<AbstractVVSVertexPtr, AbstractRFSVertexPtr > > scalar_;
+
+ /**
+ * Intermediate fermion
+ */
+ vector<pair<AbstractRFVVertexPtr, AbstractFFVVertexPtr> > fermion_;
+
+ /**
+ * Intermediate vector
+ */
+ vector<pair<AbstractVVVVertexPtr, AbstractRFVVertexPtr> > vector_;
+
+ /**
+ * Four point vertices
+ */
+ vector<AbstractRFVVVertexPtr> four_;
+ //@}
+
+};
+
+}
+
+#endif /* Herwig_MEvv2rf_H */
diff --git a/MatrixElement/General/Makefile.am b/MatrixElement/General/Makefile.am
--- a/MatrixElement/General/Makefile.am
+++ b/MatrixElement/General/Makefile.am
@@ -1,20 +1,21 @@
noinst_LTLIBRARIES = libHwGeneralME.la
libHwGeneralME_la_SOURCES = \
GeneralHardME.cc GeneralHardME.h GeneralHardME.fh \
MEvv2ff.cc MEvv2ff.h \
+MEvv2rf.cc MEvv2rf.h \
MEvv2ss.cc MEvv2ss.h \
MEfv2fs.cc MEfv2fs.h \
MEff2ss.cc MEff2ss.h \
MEff2ff.cc MEff2ff.h \
MEff2vv.cc MEff2vv.h \
MEfv2vf.cc MEfv2vf.h \
MEff2vs.cc MEff2vs.h \
MEff2sv.cc MEff2sv.h \
MEvv2vv.cc MEvv2vv.h \
MEvv2vs.cc MEvv2vs.h \
MEff2tv.cc MEff2tv.h \
MEvv2tv.cc MEvv2tv.h \
MEfv2tf.cc MEfv2tf.h \
GeneralfftoVH.cc GeneralfftoVH.h GeneralfftoVH.fh\
GeneralfftoffH.cc GeneralfftoffH.h GeneralfftoffH.fh\
GeneralQQHiggs.cc GeneralQQHiggs.h GeneralQQHiggs.fh
diff --git a/Models/General/TwoToTwoProcessConstructor.cc b/Models/General/TwoToTwoProcessConstructor.cc
--- a/Models/General/TwoToTwoProcessConstructor.cc
+++ b/Models/General/TwoToTwoProcessConstructor.cc
@@ -1,643 +1,644 @@
// -*- C++ -*-
//
// TwoToTwoProcessConstructor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TwoToTwoProcessConstructor class.
//
#include "TwoToTwoProcessConstructor.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include <sstream>
using std::stringstream;
using namespace Herwig;
TwoToTwoProcessConstructor::TwoToTwoProcessConstructor() :
Nout_(0), nv_(0), allDiagrams_(true),
processOption_(0), scaleChoice_(0), scaleFactor_(1.)
{}
IBPtr TwoToTwoProcessConstructor::clone() const {
return new_ptr(*this);
}
IBPtr TwoToTwoProcessConstructor::fullclone() const {
return new_ptr(*this);
}
void TwoToTwoProcessConstructor::doinit() {
HardProcessConstructor::doinit();
if((processOption_==2 || processOption_==4) &&
outgoing_.size()!=2)
throw InitException()
<< "Exclusive processes require exactly"
<< " two outgoing particles but " << outgoing_.size()
<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
<< Exception::runerror;
if(processOption_==4 && incoming_.size()!=2)
throw InitException()
<< "Exclusive processes require exactly"
<< " two incoming particles but " << incoming_.size()
<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
<< Exception::runerror;
Nout_ = outgoing_.size();
PDVector::size_type ninc = incoming_.size();
// exit if nothing to do
if(Nout_==0||ninc==0) return;
//create vector of initial-state pairs
for(PDVector::size_type i = 0; i < ninc; ++i) {
for(PDVector::size_type j = 0; j < ninc; ++j) {
tPDPair inc = make_pair(incoming_[i], incoming_[j]);
if( (inc.first->iSpin() > inc.second->iSpin()) ||
(inc.first->iSpin() == inc.second->iSpin() &&
inc.first->id() < inc.second->id()) )
swap(inc.first, inc.second);
if( !HPC_helper::duplicateIncoming(inc,incPairs_) ) {
incPairs_.push_back(inc);
}
}
}
// excluded vertices
excludedVertexSet_ =
set<VertexBasePtr>(excludedVertexVector_.begin(),
excludedVertexVector_.end());
}
void TwoToTwoProcessConstructor::persistentOutput(PersistentOStream & os) const {
os << vertices_ << incoming_ << outgoing_
<< allDiagrams_ << processOption_
<< scaleChoice_ << scaleFactor_ << excluded_ << excludedExternal_
<< excludedVertexVector_ << excludedVertexSet_;
}
void TwoToTwoProcessConstructor::persistentInput(PersistentIStream & is, int) {
is >> vertices_ >> incoming_ >> outgoing_
>> allDiagrams_ >> processOption_
>> scaleChoice_ >> scaleFactor_ >> excluded_ >> excludedExternal_
>> excludedVertexVector_ >> excludedVertexSet_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<TwoToTwoProcessConstructor,HardProcessConstructor>
describeHerwigTwoToTwoProcessConstructor("Herwig::TwoToTwoProcessConstructor", "Herwig.so");
void TwoToTwoProcessConstructor::Init() {
static ClassDocumentation<TwoToTwoProcessConstructor> documentation
("TwoToTwoProcessConstructor constructs the possible diagrams for "
"a process given the external particles");
static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceIn
("Incoming",
"Pointers to incoming particles",
&TwoToTwoProcessConstructor::incoming_, -1, false, false, true, false);
static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceOut
("Outgoing",
"Pointers to incoming particles",
&TwoToTwoProcessConstructor::outgoing_, -1, false, false, true, false);
static Switch<TwoToTwoProcessConstructor,bool> interfaceIncludeAllDiagrams
("IncludeEW",
"Switch to decide which diagrams to include in ME calc.",
&TwoToTwoProcessConstructor::allDiagrams_, true, false, false);
static SwitchOption interfaceIncludeAllDiagramsNo
(interfaceIncludeAllDiagrams,
"No",
"Only include QCD diagrams",
false);
static SwitchOption interfaceIncludeAllDiagramsYes
(interfaceIncludeAllDiagrams,
"Yes",
"Include EW+QCD.",
true);
static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceProcesses
("Processes",
"Whether to generate inclusive or exclusive processes",
&TwoToTwoProcessConstructor::processOption_, 0, false, false);
static SwitchOption interfaceProcessesSingleParticleInclusive
(interfaceProcesses,
"SingleParticleInclusive",
"Require at least one particle from the list of outgoing particles"
" in the hard process",
0);
static SwitchOption interfaceProcessesTwoParticleInclusive
(interfaceProcesses,
"TwoParticleInclusive",
"Require that both the particles in the hard processes are in the"
" list of outgoing particles",
1);
static SwitchOption interfaceProcessesExclusive
(interfaceProcesses,
"Exclusive",
"Require that both the particles in the hard processes are in the"
" list of outgoing particles in every hard process",
2);
static SwitchOption interfaceProcessesVeryExclusive
(interfaceProcesses,
"VeryExclusive",
"Require that both the incoming and outgoing particles in the hard processes are in the"
" list of outgoing particles in every hard process",
4);
static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceScaleChoice
("ScaleChoice",
"&TwoToTwoProcessConstructor::scaleChoice_",
&TwoToTwoProcessConstructor::scaleChoice_, 0, false, false);
static SwitchOption interfaceScaleChoiceDefault
(interfaceScaleChoice,
"Default",
"Use if sHat if intermediates all colour neutral, otherwise the transverse mass",
0);
static SwitchOption interfaceScaleChoicesHat
(interfaceScaleChoice,
"sHat",
"Always use sHat",
1);
static SwitchOption interfaceScaleChoiceTransverseMass
(interfaceScaleChoice,
"TransverseMass",
"Always use the transverse mass",
2);
static SwitchOption interfaceScaleChoiceGeometicMean
(interfaceScaleChoice,
"MaxMT",
"Use the maximum of m^2+p_T^2 for the two particles",
3);
static Parameter<TwoToTwoProcessConstructor,double> interfaceScaleFactor
("ScaleFactor",
"The prefactor used in the scale calculation. The scale used is"
" that defined by scaleChoice multiplied by this prefactor",
&TwoToTwoProcessConstructor::scaleFactor_, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceExcluded
("Excluded",
"Particles which are not allowed as intermediates",
&TwoToTwoProcessConstructor::excluded_, -1, false, false, true, false, false);
static RefVector<TwoToTwoProcessConstructor,ParticleData> interfaceExcludedExternal
("ExcludedExternal",
"Particles which are not allowed as outgoing particles",
&TwoToTwoProcessConstructor::excludedExternal_, -1,
false, false, true, false, false);
static RefVector<TwoToTwoProcessConstructor,VertexBase> interfaceExcludedVertices
("ExcludedVertices",
"Vertices which are not included in the 2 -> 2 scatterings",
&TwoToTwoProcessConstructor::excludedVertexVector_, -1, false, false, true, true, false);
}
void TwoToTwoProcessConstructor::constructDiagrams() {
if(incPairs_.empty() || outgoing_.empty() || !subProcess() ) return;
nv_ = model()->numberOfVertices();
//make sure vertices are initialised
for(unsigned int ix = 0; ix < nv_; ++ix ) {
VertexBasePtr vertex = model()->vertex(ix);
if(excludedVertexSet_.find(vertex) !=
excludedVertexSet_.end()) continue;
vertices_.push_back(vertex);
}
nv_ = vertices_.size();
//Create necessary diagrams
vector<tcPDPair>::size_type is;
PDVector::size_type os;
for(is = 0; is < incPairs_.size(); ++is) {
tPDPair ppi = incPairs_[is];
for(os = 0; os < Nout_; ++os) {
long fs = outgoing_[os]->id();
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexA = vertices_[iv];
//This skips an effective vertex and the EW ones if
// we only want the strong diagrams
if( !allDiagrams_ && vertexA->orderInGs() == 0 )
continue;
if(vertexA->getNpoint() == 3) {
//scattering diagrams
createTChannels(ppi, fs, vertexA);
//resonance diagrams
if( vertexA->isIncoming(ppi.first) &&
vertexA->isIncoming(ppi.second) )
createSChannels(ppi, fs, vertexA);
}
else {
makeFourPointDiagrams(ppi.first->id(), ppi.second->id(),
fs, vertexA);
}
}
}
}
// need to find all of the diagrams that relate to the same process
// first insert them into a map which uses the '<' operator
// to sort the diagrams
multiset<HPDiagram> grouped;
HPDVector::iterator dit = processes_.begin();
HPDVector::iterator dend = processes_.end();
for( ; dit != dend; ++dit) {
grouped.insert(*dit);
}
assert( processes_.size() == grouped.size() );
processes_.clear();
typedef multiset<HPDiagram>::const_iterator set_iter;
set_iter it = grouped.begin(), iend = grouped.end();
while( it != iend ) {
pair<set_iter,set_iter> range = grouped.equal_range(*it);
set_iter itb = range.first;
HPDVector process;
for( ; itb != range.second; ++itb ) {
process.push_back(*itb);
}
// if inclusive enforce the exclusivity
if(processOption_==2 || processOption_==4) {
if(!((process[0].outgoing. first==outgoing_[0]->id()&&
process[0].outgoing.second==outgoing_[1]->id())||
(process[0].outgoing. first==outgoing_[1]->id()&&
process[0].outgoing.second==outgoing_[0]->id()))) {
process.clear();
it = range.second;
continue;
}
if(processOption_==4) {
if(!((process[0].incoming. first==incoming_[0]->id()&&
process[0].incoming.second==incoming_[1]->id())||
(process[0].incoming. first==incoming_[1]->id()&&
process[0].incoming.second==incoming_[0]->id()))) {
process.clear();
it = range.second;
continue;
}
}
}
// check no zero width s-channel intermediates
for( dit=process.begin(); dit != process.end(); ++dit) {
tPDPtr out1 = getParticleData(dit->outgoing.first );
tPDPtr out2 = getParticleData(dit->outgoing.second);
if(dit->channelType == HPDiagram::sChannel &&
dit->intermediate->width()==ZERO &&
dit->intermediate->mass() > out1->mass()+ out2->mass()) {
tPDPtr in1 = getParticleData(dit->incoming.first );
tPDPtr in2 = getParticleData(dit->incoming.second);
throw Exception() << "Process with zero width resonant intermediates\n"
<< dit->intermediate->PDGName()
<< " can be on-shell in the process "
<< in1 ->PDGName() << " " << in2->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< " but has zero width.\nEither set the width, enable "
<< "calculation of its decays, and hence the width,\n"
<< "or disable it as a potential intermediate using\n"
<< "insert " << fullName() << ":Excluded 0 "
<< dit->intermediate->fullName() << "\n---\n"
<< Exception::runerror;
}
}
if(find(excludedExternal_.begin(),excludedExternal_.end(),
getParticleData(process[0].outgoing. first))!=excludedExternal_.end()) {
process.clear();
it = range.second;
continue;
}
if(find(excludedExternal_.begin(),excludedExternal_.end(),
getParticleData(process[0].outgoing.second))!=excludedExternal_.end()) {
process.clear();
it = range.second;
continue;
}
// finally if the process is allow assign the colour flows
for(unsigned int ix=0;ix<process.size();++ix) assignToCF(process[ix]);
// create the matrix element
createMatrixElement(process);
process.clear();
it = range.second;
}
}
void TwoToTwoProcessConstructor::
createSChannels(tcPDPair inpp, long fs, tVertexBasePtr vertex) {
//Have 2 incoming particle and a vertex, find the possible offshell
//particles
pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
tPDSet offshells = search(vertex, inpp.first->id(), incoming,
inpp.second->id(), incoming, outgoing);
tPDSet::const_iterator it;
for(it = offshells.begin(); it != offshells.end(); ++it) {
if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexB = vertices_[iv];
if( vertexB->getNpoint() != 3) continue;
if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
tPDSet final;
if( vertexB->isOutgoing(getParticleData(fs)) &&
vertexB->isIncoming(*it) )
final = search(vertexB, (*it)->id(), incoming, fs,
outgoing, outgoing);
//Now make diagrams
if(!final.empty())
makeDiagrams(inc, fs, final, *it, HPDiagram::sChannel,
make_pair(vertex, vertexB), make_pair(true,true));
}
}
}
void TwoToTwoProcessConstructor::
createTChannels(tPDPair inpp, long fs, tVertexBasePtr vertex) {
pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
//first try a with c
tPDSet offshells = search(vertex, inpp.first->id(), incoming, fs,
outgoing, outgoing);
tPDSet::const_iterator it;
for(it = offshells.begin(); it != offshells.end(); ++it) {
if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexB = vertices_[iv];
if( vertexB->getNpoint() != 3 ) continue;
if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
tPDSet final;
if( vertexB->isIncoming(inpp.second) )
final = search(vertexB, inc.second, incoming, (*it)->id(),
incoming, outgoing);
if( !final.empty() )
makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
make_pair(vertex, vertexB), make_pair(true,true));
}
}
//now try b with c
offshells = search(vertex, inpp.second->id(), incoming, fs,
outgoing, incoming);
for(it = offshells.begin(); it != offshells.end(); ++it) {
if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexB = vertices_[iv];
if( vertexB->getNpoint() != 3 ) continue;
if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
tPDSet final;
if( vertexB->isIncoming(inpp.first) )
final = search(vertexB, inc.first, incoming, (*it)->id(),
outgoing, outgoing);
if( !final.empty() )
makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
make_pair(vertexB, vertex), make_pair(true, false));
}
}
}
void TwoToTwoProcessConstructor::makeFourPointDiagrams(long parta, long partb,
long partc, VBPtr vert) {
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
getParticleData(partc));
if(loc==outgoing_.end()) return;
}
tPDSet ext = search(vert, parta, incoming, partb,incoming, partc, outgoing);
if( ext.empty() ) return;
IDPair in(parta, partb);
for(tPDSet::const_iterator iter=ext.begin(); iter!=ext.end();
++iter) {
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
*iter);
if(loc==outgoing_.end()) continue;
}
HPDiagram nhp(in,make_pair(partc, (*iter)->id()));
nhp.vertices = make_pair(vert, vert);
nhp.channelType = HPDiagram::fourPoint;
fixFSOrder(nhp);
if(!checkOrder(nhp)) continue;
if( !duplicate(nhp, processes_) ) processes_.push_back(nhp);
}
}
void
TwoToTwoProcessConstructor::makeDiagrams(IDPair in, long out1, const tPDSet & out2,
PDPtr inter, HPDiagram::Channel chan,
VBPair vertexpair, BPair cross) {
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
getParticleData(out1));
if(loc==outgoing_.end()) return;
}
for(tPDSet::const_iterator it = out2.begin(); it != out2.end(); ++it) {
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
*it);
if(loc==outgoing_.end()) continue;
}
HPDiagram nhp( in,make_pair(out1, (*it)->id()) );
nhp.intermediate = inter;
nhp.vertices = vertexpair;
nhp.channelType = chan;
nhp.ordered = cross;
fixFSOrder(nhp);
if(!checkOrder(nhp)) continue;
if( !duplicate(nhp, processes_) ) processes_.push_back(nhp);
}
}
set<tPDPtr>
TwoToTwoProcessConstructor::search(VBPtr vertex, long part1, direction d1,
long part2, direction d2, direction d3) {
if(vertex->getNpoint() != 3) return tPDSet();
if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
vector<long> ext;
tPDSet third;
for(unsigned int ix = 0;ix < 3; ++ix) {
vector<long> pdlist = vertex->search(ix, part1);
ext.insert(ext.end(), pdlist.begin(), pdlist.end());
}
for(unsigned int ix = 0; ix < ext.size(); ix += 3) {
long id0 = ext.at(ix);
long id1 = ext.at(ix+1);
long id2 = ext.at(ix+2);
int pos;
if((id0 == part1 && id1 == part2) ||
(id0 == part2 && id1 == part1))
pos = ix + 2;
else if((id0 == part1 && id2 == part2) ||
(id0 == part2 && id2 == part1))
pos = ix + 1;
else if((id1 == part1 && id2 == part2) ||
(id1 == part2 && id2 == part1))
pos = ix;
else
pos = -1;
if(pos >= 0) {
tPDPtr p = getParticleData(ext[pos]);
if(d3 == incoming && p->CC()) p = p->CC();
third.insert(p);
}
}
return third;
}
set<tPDPtr>
TwoToTwoProcessConstructor::search(VBPtr vertex,
long part1, direction d1,
long part2, direction d2,
long part3, direction d3,
direction d4) {
if(vertex->getNpoint() != 4) return tPDSet();
if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
if(d3 == incoming && getParticleData(part3)->CC()) part3 = -part3;
vector<long> ext;
tPDSet fourth;
for(unsigned int ix = 0;ix < 4; ++ix) {
vector<long> pdlist = vertex->search(ix, part1);
ext.insert(ext.end(), pdlist.begin(), pdlist.end());
}
for(unsigned int ix = 0;ix < ext.size(); ix += 4) {
long id0 = ext.at(ix); long id1 = ext.at(ix + 1);
long id2 = ext.at(ix + 2); long id3 = ext.at(ix + 3);
int pos;
if((id0 == part1 && id1 == part2 && id2 == part3) ||
(id0 == part1 && id1 == part3 && id2 == part2) ||
(id0 == part2 && id1 == part1 && id2 == part3) ||
(id0 == part2 && id1 == part3 && id2 == part1) ||
(id0 == part3 && id1 == part1 && id2 == part2) ||
(id0 == part3 && id1 == part2 && id2 == part1))
pos = ix + 3;
else if((id0 == part1 && id1 == part2 && id3 == part3) ||
(id0 == part1 && id1 == part3 && id3 == part2) ||
(id0 == part2 && id1 == part1 && id3 == part3) ||
(id0 == part2 && id1 == part3 && id3 == part1) ||
(id0 == part3 && id1 == part1 && id3 == part2) ||
(id0 == part3 && id1 == part2 && id3 == part1))
pos = ix + 2;
else if((id0 == part1 && id2 == part2 && id3 == part3) ||
(id0 == part1 && id2 == part3 && id3 == part2) ||
(id0 == part2 && id2 == part1 && id3 == part3) ||
(id0 == part2 && id2 == part3 && id3 == part1) ||
(id0 == part3 && id2 == part1 && id3 == part2) ||
(id0 == part3 && id2 == part2 && id3 == part1))
pos = ix + 1;
else if((id1 == part1 && id2 == part2 && id3 == part3) ||
(id1 == part1 && id2 == part3 && id3 == part2) ||
(id1 == part2 && id2 == part1 && id3 == part3) ||
(id1 == part2 && id2 == part3 && id3 == part1) ||
(id1 == part3 && id2 == part1 && id3 == part2) ||
(id1 == part3 && id2 == part2 && id3 == part1))
pos = ix;
else
pos = -1;
if(pos >= 0) {
tPDPtr p = getParticleData(ext[pos]);
if(d4 == incoming && p->CC())
p = p->CC();
fourth.insert(p);
}
}
return fourth;
}
void
TwoToTwoProcessConstructor::createMatrixElement(const HPDVector & process) const {
if ( process.empty() ) return;
// external particles
tcPDVector extpart(4);
extpart[0] = getParticleData(process[0].incoming.first);
extpart[1] = getParticleData(process[0].incoming.second);
extpart[2] = getParticleData(process[0].outgoing.first);
extpart[3] = getParticleData(process[0].outgoing.second);
// create the object
string objectname ("/Herwig/MatrixElements/");
string classname = MEClassname(extpart, objectname);
GeneralHardMEPtr matrixElement = dynamic_ptr_cast<GeneralHardMEPtr>
(generator()->preinitCreate(classname, objectname));
if( !matrixElement ) {
std::stringstream message;
message << "TwoToTwoProcessConstructor::createMatrixElement "
<< "- No matrix element object could be created for "
<< "the process "
<< extpart[0]->PDGName() << " " << extpart[0]->iSpin() << ","
<< extpart[1]->PDGName() << " " << extpart[1]->iSpin() << "->"
<< extpart[2]->PDGName() << " " << extpart[2]->iSpin() << ","
<< extpart[3]->PDGName() << " " << extpart[3]->iSpin()
<< ". Constructed class name: \"" << classname << "\"";
generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
return;
}
// choice for the scale
unsigned int scale;
if(scaleChoice_==0) {
// check coloured initial and final state
bool inColour = ( extpart[0]->coloured() ||
extpart[1]->coloured());
bool outColour = ( extpart[2]->coloured() ||
extpart[3]->coloured());
if(inColour&&outColour) {
bool coloured = false;
for(unsigned int ix=0;ix<process.size();++ix) {
if(process[ix].intermediate&&
process[ix].intermediate->coloured()) {
coloured = true;
break;
}
}
scale = coloured ? 1 : 0;
}
else {
scale = 0;
}
}
else {
scale = scaleChoice_-1;
}
// set the information
matrixElement->setProcessInfo(process, colourFlow(extpart),
debug(), scale, scaleFactor_);
// insert it
generator()->preinitInterface(subProcess(), "MatrixElements",
subProcess()->MEs().size(),
"insert", matrixElement->fullName());
}
string TwoToTwoProcessConstructor::MEClassname(const vector<tcPDPtr> & extpart,
string & objname) const {
string classname("Herwig::ME");
for(vector<tcPDPtr>::size_type ix = 0; ix < extpart.size(); ++ix) {
if(ix == 2) classname += "2";
if(extpart[ix]->iSpin() == PDT::Spin0) classname += "s";
else if(extpart[ix]->iSpin() == PDT::Spin1) classname += "v";
else if(extpart[ix]->iSpin() == PDT::Spin1Half) classname += "f";
+ else if(extpart[ix]->iSpin() == PDT::Spin3Half) classname += "r";
else if(extpart[ix]->iSpin() == PDT::Spin2) classname += "t";
else {
std::stringstream message;
message << "MEClassname() : Encountered an unknown spin for "
<< extpart[ix]->PDGName() << " while constructing MatrixElement "
<< "classname " << extpart[ix]->iSpin();
generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
}
}
objname += "ME" + extpart[0]->PDGName() + extpart[1]->PDGName() + "2"
+ extpart[2]->PDGName() + extpart[3]->PDGName();
return classname;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:36 PM (15 h, 47 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486658
Default Alt Text
(49 KB)

Event Timeline