Page MenuHomeHEPForge
Files F10881516

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/General/MEff2ss.cc b/MatrixElement/General/MEff2ss.cc
--- a/MatrixElement/General/MEff2ss.cc
+++ b/MatrixElement/General/MEff2ss.cc
@@ -1,300 +1,310 @@
// -*- C++ -*-
//
// MEff2ss.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 MEff2ss class.
//
#include "MEff2ss.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using ThePEG::Helicity::VectorWaveFunction;
using ThePEG::Helicity::TensorWaveFunction;
using ThePEG::Helicity::incoming;
using ThePEG::Helicity::outgoing;
void MEff2ss::doinit() {
GeneralHardME::doinit();
fermion_.resize(numberOfDiags());
+ scalar_ .resize(numberOfDiags());
vector_ .resize(numberOfDiags());
tensor_ .resize(numberOfDiags());
flowME().resize(numberOfFlows(),
ProductionMatrixElement(PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin0 , PDT::Spin0 ));
diagramME().resize(numberOfDiags(),
ProductionMatrixElement(PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin0 , PDT::Spin0 ));
for(HPCount i = 0; i < numberOfDiags(); ++i) {
const HPDiagram & current = getProcessInfo()[i];
if(current.channelType == HPDiagram::tChannel) {
if(current.intermediate->iSpin() == PDT::Spin1Half)
fermion_[i] =
make_pair(dynamic_ptr_cast<AbstractFFSVertexPtr>(current.vertices.first),
dynamic_ptr_cast<AbstractFFSVertexPtr>(current.vertices.second));
else
throw InitException() << "MEFF2ss:doinit() - t-channel"
<< " intermediate must be a fermion "
<< Exception::runerror;
}
else if(current.channelType == HPDiagram::sChannel) {
- if(current.intermediate->iSpin() == PDT::Spin1)
+ if(current.intermediate->iSpin() == PDT::Spin0)
+ scalar_[i] =
+ make_pair(dynamic_ptr_cast<AbstractFFSVertexPtr>(current.vertices.first),
+ dynamic_ptr_cast<AbstractSSSVertexPtr>(current.vertices.second));
+ else if(current.intermediate->iSpin() == PDT::Spin1)
vector_[i] =
make_pair(dynamic_ptr_cast<AbstractFFVVertexPtr>(current.vertices.first),
dynamic_ptr_cast<AbstractVSSVertexPtr>(current.vertices.second));
else if(current.intermediate->iSpin() == PDT::Spin2)
tensor_[i] =
make_pair(dynamic_ptr_cast<AbstractFFTVertexPtr>(current.vertices.first),
dynamic_ptr_cast<AbstractSSTVertexPtr>(current.vertices.second));
else
throw InitException() << "MEFF2ss:doinit() - s-channel"
<< " intermediate must be a vector or tensor "
<< Exception::runerror;
}
else
throw InitException() << "MEFF2ss:doinit() - Cannot find correct "
<< "channel from diagram. Vertex not cast! "
<< Exception::runerror;
}
}
void MEff2ss::doinitrun() {
GeneralHardME::doinitrun();
flowME().resize(numberOfFlows(),
ProductionMatrixElement(PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin0 , PDT::Spin0 ));
diagramME().resize(numberOfDiags(),
ProductionMatrixElement(PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin0 , PDT::Spin0 ));
}
double MEff2ss::me2() const {
// first setup wavefunctions for external particles
SpinorVector sp(2);
SpinorBarVector sbar(2);
for( unsigned int i = 0; i < 2; ++i ) {
sp[i] = SpinorWaveFunction (rescaledMomenta()[0],
mePartonData()[0], i, incoming);
sbar[i] = SpinorBarWaveFunction(rescaledMomenta()[1],
mePartonData()[1], i, incoming);
}
ScalarWaveFunction sca1(rescaledMomenta()[2],
mePartonData()[2], 1., outgoing);
ScalarWaveFunction sca2(rescaledMomenta()[3],
mePartonData()[3], 1., outgoing);
// calculate the ME
double full_me(0.);
ff2ssME(sp, sbar, sca1, sca2, full_me,true);
// debugging tests if needed
#ifndef NDEBUG
if( debugME() ) debug(full_me);
#endif
// return the answer
return full_me;
}
ProductionMatrixElement
MEff2ss::ff2ssME(const SpinorVector & sp, const SpinorBarVector & sbar,
const ScalarWaveFunction & sca1,
const ScalarWaveFunction & sca2,
double & me2, bool first) const {
// scale
const Energy2 q2(scale());
// 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 if1 = 0; if1 < 2; ++if1) {
for(unsigned int if2 = 0; if2 < 2; ++if2) {
vector<Complex> flows(numberOfFlows(),0.);
for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
Complex diag(0.);
const HPDiagram & current = getProcessInfo()[ix];
tcPDPtr internal(current.intermediate);
if(current.channelType == HPDiagram::tChannel &&
internal->iSpin() == PDT::Spin1Half) {
if(current.ordered.second) {
SpinorBarWaveFunction interFB = fermion_[ix].second->
evaluate(q2, 3, internal, sbar[if2], sca2);
diag = fermion_[ix].first->evaluate(q2, sp[if1], interFB, sca1);
}
else {
SpinorBarWaveFunction interFB = fermion_[ix].second->
evaluate(q2, 3, internal, sbar[if2], sca1);
diag = fermion_[ix].first->evaluate(q2, sp[if1], interFB, sca2);
}
}
else if(current.channelType == HPDiagram::sChannel) {
- if(internal->iSpin() == PDT::Spin1) {
+ if(internal->iSpin() == PDT::Spin0) {
+ ScalarWaveFunction interS = scalar_[ix].first->
+ evaluate(q2, 1, internal, sp[if1], sbar[if2]);
+ diag = scalar_[ix].second->evaluate(q2, interS, sca2, sca1);
+ }
+ else if(internal->iSpin() == PDT::Spin1) {
VectorWaveFunction interV = vector_[ix].first->
evaluate(q2, 1, internal, sp[if1], sbar[if2]);
diag = vector_[ix].second->evaluate(q2, interV, sca2, sca1);
}
else if(internal->iSpin() == PDT::Spin2) {
TensorWaveFunction interT = tensor_[ix].first->
evaluate(q2, 1, internal, sp[if1], sbar[if2]);
diag = tensor_[ix].second ->evaluate(q2, sca2, sca1, interT);
}
}
// diagram
me[ix] += norm(diag);
diagramME()[ix](if1,if2,0,0) = diag;
// contributions to the different colour flows
for(unsigned int iy = 0; iy < current.colourFlow.size(); ++iy) {
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](if1,if2,0,0) = flows[iy];
// contribution to the squared matrix element
for(unsigned int ii = 0; ii < numberOfFlows(); ++ii)
for(unsigned int 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 MEff2ss::persistentOutput(PersistentOStream & os) const {
- os << fermion_ << vector_ << tensor_;
+ os << fermion_ << scalar_ << vector_ << tensor_;
}
void MEff2ss::persistentInput(PersistentIStream & is, int) {
- is >> fermion_ >> vector_ >> tensor_;
+ is >> fermion_ >> scalar_ >> vector_ >> tensor_;
}
ClassDescription<MEff2ss> MEff2ss::initMEff2ss;
// Definition of the static class description member.
void MEff2ss::Init() {
static ClassDocumentation<MEff2ss> documentation
("MEff2ss implements the ME calculation of the fermion-antifermion "
"to scalar-scalar hard process.");
}
void MEff2ss::constructVertex(tSubProPtr sub) {
//get particles
ParticleVector ext = hardParticles(sub);
//First calculate wave functions with off-shell momenta
//to calculate correct spin information
SpinorVector sp;
SpinorBarVector sbar;
SpinorWaveFunction (sp , ext[0], incoming, false);
SpinorBarWaveFunction (sbar, ext[1], incoming, false);
ScalarWaveFunction sca1( ext[2], outgoing, true);
ScalarWaveFunction sca2( ext[3], outgoing, true);
// Need to use rescale momenta to calculate matrix element
setRescaledMomenta(ext);
// wave functions with rescaled momenta
SpinorWaveFunction spr(rescaledMomenta()[0],
ext[0]->dataPtr(), incoming);
SpinorBarWaveFunction sbr(rescaledMomenta()[1],
ext[1]->dataPtr(), incoming);
sca1 = ScalarWaveFunction(rescaledMomenta()[2],
ext[2]->dataPtr(), outgoing);
sca2 = ScalarWaveFunction(rescaledMomenta()[3],
ext[3]->dataPtr(), outgoing);
for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
spr.reset(ihel);
sp[ihel] = spr;
sbr.reset(ihel);
sbar[ihel] = sbr;
}
double dummy(0.);
ProductionMatrixElement pme = ff2ssME(sp, sbar, sca1, sca2, dummy,false);
#ifndef NDEBUG
if( debugME() ) debug(dummy/36.);
#endif
createVertex(pme,ext);
}
void MEff2ss::debug(double me2) const {
if( !generator()->logfile().is_open() ) return;
long id1 = mePartonData()[0]->id();
long id2 = mePartonData()[1]->id();
long id3 = mePartonData()[2]->id();
long id4 = mePartonData()[3]->id();
if( (abs(id1) != 1 && abs(id1) != 2) || (abs(id2) != 1 && abs(id2) != 2) ||
( abs(id3) != 1000001 && abs(id3) != 1000002 &&
abs(id3) != 2000001 && abs(id3) != 2000002 ) ||
( abs(id4) != 1000001 && abs(id4) != 1000002 &&
abs(id4) != 2000001 && abs(id4) != 2000002 ) ) return;
tcSMPtr sm = generator()->standardModel();
double gs4 = sqr( 4.*Constants::pi*sm->alphaS(scale()) );
int Nc = sm->Nc();
double Cf = (sqr(Nc) - 1)/2./Nc;
Energy2 s(sHat());
Energy2 mgos = sqr( getParticleData(ParticleID::SUSY_g)->mass());
Energy2 m3s = sqr(mePartonData()[2]->mass());
Energy2 m4s = sqr(mePartonData()[3]->mass());
Energy4 spt2 = uHat()*tHat() - m3s*m4s;
Energy2 tgl(tHat() - mgos), ugl(uHat() - mgos);
unsigned int alpha = abs(id3)/1000000;
unsigned int beta = abs(id4)/1000000;
bool iflav = ( abs(id1) == abs(id2) );
unsigned int oflav = ( abs(id3) - abs(id1) ) % 10;
double analytic(0.);
if( alpha != beta ) {
if( ( id1 > 0 && id2 > 0) ||
( id1 < 0 && id2 < 0) ) {
analytic = spt2/sqr(tgl);
if( iflav ) analytic += spt2/sqr(ugl);
}
else {
analytic = s*mgos/sqr(tgl);
}
}
else {
if( oflav != 0 ) {
analytic = 2.*spt2/sqr(s);
}
else if( ( id1 > 0 && id2 > 0) ||
( id1 < 0 && id2 < 0) ) {
analytic = s*mgos/sqr(tgl);
if( iflav ) {
analytic += s*mgos/sqr(ugl) - 2.*s*mgos/Nc/tgl/ugl;
}
analytic /= ( iflav ? 2. : 1.);
}
else {
analytic = spt2/sqr(tgl);
if( iflav ) {
analytic += 2.*spt2/sqr(s) - 2.*spt2/Nc/s/tgl;
}
}
}
analytic *= gs4*Cf/2./Nc;
double diff = abs(analytic - me2);
if( diff > 1e-4 ) {
generator()->log()
<< mePartonData()[0]->PDGName() << ","
<< mePartonData()[1]->PDGName() << "->"
<< mePartonData()[2]->PDGName() << ","
<< mePartonData()[3]->PDGName() << " difference: "
<< setprecision(10) << diff << " ratio: " << analytic/me2
<< '\n';
}
}
diff --git a/MatrixElement/General/MEff2ss.h b/MatrixElement/General/MEff2ss.h
--- a/MatrixElement/General/MEff2ss.h
+++ b/MatrixElement/General/MEff2ss.h
@@ -1,230 +1,234 @@
// -*- C++ -*-
//
// MEff2ss.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MEff2ss_H
#define HERWIG_MEff2ss_H
//
// This is the declaration of the MEff2ss class.
//
#include "GeneralHardME.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.h"
-#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
-#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
-#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
using ThePEG::Helicity::SpinorWaveFunction;
using ThePEG::Helicity::SpinorBarWaveFunction;
using ThePEG::Helicity::ScalarWaveFunction;
/**
* The MEff2ss class is designed to implement the matrix element for a
* fermion-antifermion to scalar-scalar hard process. It inherits from
* GeneralHardME and implements the appropriate virtual functions for this
* specific spin combination.
*
* @see \ref MEff2ssInterfaces "The interfaces"
* defined for MEff2ss.
* @see GeneralHardME
*/
class MEff2ss: public GeneralHardME {
public:
/** Vector of SpinorWaveFunctions objects */
typedef vector<SpinorWaveFunction> SpinorVector;
/** Vector of SpinorBarWaveFunction objects. */
typedef vector<SpinorBarWaveFunction> SpinorBarVector;
public:
/**
* The default constructor.
*/
MEff2ss() : fermion_(0), vector_(0), tensor_(0) {}
/** @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);
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();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
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 static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEff2ss> initMEff2ss;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEff2ss & operator=(const MEff2ss &);
private:
/**
* Calculate the matrix element
* @param sp A vector of SpinorWaveFunction objects
* @param sbar A vector of SpinorBarWaveFunction objects
* @param sca1 A ScalarWaveFunction for an outgoing scalar
* @param sca2 A ScalarWaveFunction for the other outgoing scalar
* @param me2 The spin averaged matrix element
* @param first Whether or not first call to decide if colour decomposition etc
* should be calculated
*/
ProductionMatrixElement ff2ssME(const SpinorVector & sp,
const SpinorBarVector & sbar,
const ScalarWaveFunction & sca1,
const ScalarWaveFunction & sca2,
double & me2, bool first) const;
private:
/**
* Storage for dynamically cast vertices for a diagram with intermediate
* fermion
*/
vector<pair<AbstractFFSVertexPtr, AbstractFFSVertexPtr> > fermion_;
/**
* Storage for dynamically cast vertices for a diagram with intermediate
* vector
*/
+ vector<pair<AbstractFFSVertexPtr, AbstractSSSVertexPtr> > scalar_;
+
+ /**
+ * Storage for dynamically cast vertices for a diagram with intermediate
+ * vector
+ */
vector<pair<AbstractFFVVertexPtr, AbstractVSSVertexPtr> > vector_;
/**
* Storage for dynamically cast vertices for a diagram with intermediate
* tensor
*/
vector<pair<AbstractFFTVertexPtr, AbstractSSTVertexPtr> > tensor_;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEff2ss. */
template <>
struct BaseClassTrait<Herwig::MEff2ss,1> {
/** Typedef of the first base class of MEff2ss. */
typedef Herwig::GeneralHardME NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEff2ss class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEff2ss>
: public ClassTraitsBase<Herwig::MEff2ss> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEff2ss"; }
};
/** @endcond */
}
#endif /* HERWIG_MEff2ss_H */
diff --git a/MatrixElement/General/MEvv2ss.cc b/MatrixElement/General/MEvv2ss.cc
--- a/MatrixElement/General/MEvv2ss.cc
+++ b/MatrixElement/General/MEvv2ss.cc
@@ -1,262 +1,275 @@
// -*- C++ -*-
//
// MEvv2ss.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 MEvv2ss class.
//
#include "MEvv2ss.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using ThePEG::Helicity::TensorWaveFunction;
using ThePEG::Helicity::incoming;
using ThePEG::Helicity::outgoing;
using ThePEG::Helicity::SpinfoPtr;
void MEvv2ss::doinit() {
GeneralHardME::doinit();
- scalar_.resize(numberOfDiags());
- vector_.resize(numberOfDiags());
- tensor_.resize(numberOfDiags());
+ scalar1_.resize(numberOfDiags());
+ scalar2_.resize(numberOfDiags());
+ vector_ .resize(numberOfDiags());
+ tensor_ .resize(numberOfDiags());
flowME().resize(numberOfFlows(),
ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
PDT::Spin0, PDT::Spin0));
diagramME().resize(numberOfDiags(),
ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
PDT::Spin0, PDT::Spin0));
for(size_t i = 0; i < numberOfDiags(); ++i ) {
HPDiagram dg = getProcessInfo()[i];
if( !dg.intermediate ) {
contact_ = dynamic_ptr_cast<AbstractVVSSVertexPtr>(dg.vertices.first);
}
else if(dg.channelType == HPDiagram::tChannel) {
AbstractVSSVertexPtr vss1 =
dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.first);
AbstractVSSVertexPtr vss2 =
dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.second);
- scalar_[i] = make_pair(vss1, vss2);
+ scalar2_[i] = make_pair(vss1, vss2);
}
else {
- if( dg.intermediate->iSpin() == PDT::Spin1 ) {
+ if( dg.intermediate->iSpin() == PDT::Spin0 ) {
+ AbstractVVSVertexPtr vvs =
+ dynamic_ptr_cast<AbstractVVSVertexPtr>(dg.vertices.first);
+ AbstractSSSVertexPtr sss =
+ dynamic_ptr_cast<AbstractSSSVertexPtr>(dg.vertices.second);
+ scalar1_[i] = make_pair(vvs, sss);
+ }
+ else if( dg.intermediate->iSpin() == PDT::Spin1 ) {
AbstractVVVVertexPtr vvv =
dynamic_ptr_cast<AbstractVVVVertexPtr>(dg.vertices.first);
AbstractVSSVertexPtr vss =
dynamic_ptr_cast<AbstractVSSVertexPtr>(dg.vertices.second);
vector_[i] = make_pair(vvv, vss);
}
else if( dg.intermediate->iSpin() == PDT::Spin2 ) {
AbstractVVTVertexPtr vvt =
dynamic_ptr_cast<AbstractVVTVertexPtr>(dg.vertices.first);
AbstractSSTVertexPtr sst =
dynamic_ptr_cast<AbstractSSTVertexPtr>(dg.vertices.second);
tensor_[i] = make_pair(vvt, sst);
}
}
}
}
void MEvv2ss::doinitrun() {
GeneralHardME::doinitrun();
flowME().resize(numberOfFlows(),
ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
PDT::Spin0, PDT::Spin0));
diagramME().resize(numberOfDiags(),
ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
PDT::Spin0, PDT::Spin0));
}
double MEvv2ss::me2() const {
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);
}
ScalarWaveFunction sca1(rescaledMomenta()[2],mePartonData()[2],
Complex(1.,0.),outgoing);
ScalarWaveFunction sca2(rescaledMomenta()[3],mePartonData()[3],
Complex(1.,0.),outgoing);
double full_me(0.);
vv2ssME(v1, v2, sca1, sca2, full_me , true);
#ifndef NDEBUG
if( debugME() ) debug(full_me);
#endif
return full_me;
}
ProductionMatrixElement
MEvv2ss::vv2ssME(const VBVector & v1, const VBVector & v2,
const ScalarWaveFunction & sca1,
const ScalarWaveFunction & sca2,
double & me2, bool first) const {
const Energy2 m2(scale());
const Energy masst = sca1.mass(), massu = sca2.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.);
//loop over vector helicities
for(unsigned int iv1 = 0; iv1 < 2; ++iv1) {
for(unsigned int iv2 = 0; iv2 < 2; ++iv2) {
vector<Complex> flows(numberOfFlows(),0.);
// loop over diagrams
for(HPCount ix = 0; ix < numberOfDiags(); ++ix) {
Complex diag(0.);
const HPDiagram & current = getProcessInfo()[ix];
// do four-point diag first
if(current.channelType == HPDiagram::fourPoint) {
diag = contact_->evaluate(m2, v1[iv1], v2[iv2], sca1, sca2);
}
else {
tcPDPtr offshell = current.intermediate;
if(current.channelType == HPDiagram::tChannel) {
if(current.ordered.second) {
- ScalarWaveFunction interS = scalar_[ix].first->
+ ScalarWaveFunction interS = scalar2_[ix].first->
evaluate(m2, 3, offshell, v1[iv1], sca1, masst);
- diag = scalar_[ix].second->evaluate(m2, v2[iv2], interS, sca2);
+ diag = scalar2_[ix].second->evaluate(m2, v2[iv2], interS, sca2);
}
else {
- ScalarWaveFunction interS = scalar_[ix].first->
+ ScalarWaveFunction interS = scalar2_[ix].first->
evaluate(m2, 3, offshell, v1[iv1], sca2, massu);
- diag = scalar_[ix].second->evaluate(m2, v2[iv2], interS, sca1);
+ diag = scalar2_[ix].second->evaluate(m2, v2[iv2], interS, sca1);
}
}
else if(current.channelType == HPDiagram::sChannel) {
- if(offshell->iSpin() == PDT::Spin1) {
+ if(offshell->iSpin() == PDT::Spin0) {
+ ScalarWaveFunction interS = scalar1_[ix].first->
+ evaluate(m2, 1, offshell, v1[iv1], v2[iv2]);
+ diag = scalar1_[ix].second->evaluate(m2, interS, sca1, sca2);
+ }
+ else if(offshell->iSpin() == PDT::Spin1) {
VectorWaveFunction interV = vector_[ix].first->
evaluate(m2, 1, offshell, v1[iv1], v2[iv2]);
diag = vector_[ix].second->evaluate(m2, interV, sca1, sca2);
}
else if(offshell->iSpin() == PDT::Spin2) {
TensorWaveFunction interT = tensor_[ix].first->
evaluate(m2, 1, offshell, v1[iv1], v2[iv2]);
diag = tensor_[ix].second->evaluate(m2, sca1, sca2, interT);
}
}
else
diag = 0.;
}
me[ix] += norm(diag);
diagramME()[ix](2*iv1, 2*iv2, 0, 0) = diag;
//Compute flows
for(size_t iy = 0; iy < current.colourFlow.size(); ++iy) {
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, 0, 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 MEvv2ss::persistentOutput(PersistentOStream & os) const {
- os << scalar_ << vector_ << tensor_ << contact_;
+ os << scalar1_ << scalar2_ << vector_ << tensor_ << contact_;
}
void MEvv2ss::persistentInput(PersistentIStream & is, int) {
- is >> scalar_ >> vector_ >> tensor_ >> contact_;
+ is >> scalar1_ >> scalar2_ >> vector_ >> tensor_ >> contact_;
}
ClassDescription<MEvv2ss> MEvv2ss::initMEvv2ss;
// Definition of the static class description member.
void MEvv2ss::Init() {
static ClassDocumentation<MEvv2ss> documentation
("This class implements the ME for the vector-vector to scalar-scalar "
"hard-process");
}
void MEvv2ss::constructVertex(tSubProPtr sub) {
ParticleVector ext = hardParticles(sub);
VBVector v1, v2;
// set up the wavefunctions with real momenta
VectorWaveFunction(v1, ext[0], incoming, false, true);
VectorWaveFunction(v2, ext[1], incoming, false, true);
ScalarWaveFunction sca1(ext[2], outgoing, true);
ScalarWaveFunction sca2(ext[3], outgoing, true);
// calculate rescaled moment
setRescaledMomenta(ext);
// wavefunctions with rescaled momenta
VectorWaveFunction v1r (rescaledMomenta()[0],
ext[0]->dataPtr(), incoming);
VectorWaveFunction v2r (rescaledMomenta()[1],
ext[1]->dataPtr(), incoming);
sca1 = ScalarWaveFunction(rescaledMomenta()[2],
ext[2]->dataPtr(), outgoing);
sca2 = ScalarWaveFunction(rescaledMomenta()[3],
ext[3]->dataPtr(), outgoing);
for( unsigned int ihel = 0; ihel < 2; ++ihel ) {
v1r.reset(2*ihel);
v1[ihel] = v1r;
v2r.reset(2*ihel);
v2[ihel] = v2r;
}
double dummy(0.);
ProductionMatrixElement pme = vv2ssME(v1, v2, sca1, sca2, dummy , false);
#ifndef NDEBUG
if( debugME() ) debug(dummy);
#endif
createVertex(pme,ext);
}
void MEvv2ss::debug(double me2) const {
if( !generator()->logfile().is_open() ) return;
//SUSY gg>~q~q
long id3 = abs(mePartonData()[2]->id());
long id4 = abs(mePartonData()[3]->id());
if( mePartonData()[0]->id() != 21 || mePartonData()[1]->id() != 21 ||
(id3 < 1000001 && id3 > 1000006 ) || (id3 < 2000001 && id3 > 2000006 ) ||
(id4 < 1000001 && id4 > 1000006 ) || (id4 < 2000001 && id4 > 2000006 ) )
return;
tcSMPtr sm = generator()->standardModel();
double gs4 = sqr( 4.*Constants::pi*sm->alphaS(scale()) );
int Nc = sm->Nc();
Energy4 s2 = sqr(sHat());
Energy2 m3s = meMomenta()[2].m2();
Energy2 m4s = meMomenta()[3].m2();
Energy4 spt2 = uHat()*tHat() - m3s*m4s;
Energy4 t3s = sqr(tHat() - m3s);
Energy4 u4s = sqr(uHat() - m4s);
double analytic = gs4*Nc*( sqr(spt2) + s2*m3s*m4s ) *
( u4s + t3s - s2/sqr(Nc) )/2./(sqr(Nc) - 1.)/s2/t3s/u4s;
double diff = abs(analytic - me2);
if( diff > 1e-4 ) {
generator()->log()
<< mePartonData()[0]->PDGName() << ","
<< mePartonData()[1]->PDGName() << "->"
<< mePartonData()[2]->PDGName() << ","
<< mePartonData()[3]->PDGName() << " difference: "
<< setprecision(10) << diff << " ratio: " << analytic/me2 << '\n';
}
}
diff --git a/MatrixElement/General/MEvv2ss.h b/MatrixElement/General/MEvv2ss.h
--- a/MatrixElement/General/MEvv2ss.h
+++ b/MatrixElement/General/MEvv2ss.h
@@ -1,227 +1,233 @@
// -*- C++ -*-
//
// MEvv2ss.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MEvv2ss_H
#define HERWIG_MEvv2ss_H
//
// This is the declaration of the MEvv2ss class.
//
#include "GeneralHardME.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
using ThePEG::Helicity::VectorWaveFunction;
using ThePEG::Helicity::ScalarWaveFunction;
/**
* This is the implementation of the matrix element for the process
* vector-vector to scalar-scalar. It inherits from GeneralHardME and
* implements the required virtual functions.
*
* @see \ref MEff2ffInterfaces "The Interfaces"
* defined for MEff2ff.
* @see GeneralHardME
*/
class MEvv2ss: public GeneralHardME {
public:
/** A vector of VectorWaveFunction objects*/
typedef vector<VectorWaveFunction> VBVector;
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;
//@}
/**
* Set the Hardvertex for the spin correlations
* @param sub
*/
virtual void constructVertex(tSubProPtr sub);
private:
/**
* Calculate the matrix element.
* @param v1 A vector of VectorWaveFunction objects for the first boson
* @param v2 A vector of VectorWaveFunction objects for the second boson
* @param sca1 A ScalarWaveFunction for the first outgoing
* @param sca2 A ScalarWaveFunction for the second outgoing
* @param me2 The value of the spin-summed matrix element squared
* (to be calculated)
* @param first Whether or not first call to decide if colour decomposition etc
* should be calculated
*/
ProductionMatrixElement vv2ssME(const VBVector & v1, const VBVector & v2,
const ScalarWaveFunction & sca1,
const ScalarWaveFunction & sca2,
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();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
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 static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEvv2ss> initMEvv2ss;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEvv2ss & operator=(const MEvv2ss &);
private:
/** @name The dynamically casted vertices. */
//@{
/**
+ * Intermediate s-channel scalar
+ */
+ vector<pair<AbstractVVSVertexPtr, AbstractSSSVertexPtr> > scalar1_;
+
+ /**
* Intermediate t-channel scalar
*/
- vector<pair<AbstractVSSVertexPtr, AbstractVSSVertexPtr> > scalar_;
+ vector<pair<AbstractVSSVertexPtr, AbstractVSSVertexPtr> > scalar2_;
/**
* Intermediate s-channel vector
*/
vector<pair<AbstractVVVVertexPtr, AbstractVSSVertexPtr> > vector_;
/**
* Intermediate s-channel tensor
*/
vector<pair<AbstractVVTVertexPtr, AbstractSSTVertexPtr> > tensor_;
/**
* The contact vertex
*/
AbstractVVSSVertexPtr contact_;
//@}
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEvv2ss. */
template <>
struct BaseClassTrait<Herwig::MEvv2ss,1> {
/** Typedef of the first base class of MEvv2ss. */
typedef Herwig::GeneralHardME NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEvv2ss class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEvv2ss>
: public ClassTraitsBase<Herwig::MEvv2ss> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEvv2ss"; }
};
/** @endcond */
}
#endif /* HERWIG_MEvv2ss_H */
diff --git a/MatrixElement/General/MEvv2vs.cc b/MatrixElement/General/MEvv2vs.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/General/MEvv2vs.cc
@@ -0,0 +1,261 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the MEvv2vs class.
+//
+
+#include "MEvv2vs.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;
+using ThePEG::Helicity::SpinfoPtr;
+
+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_;
+}
+
+void MEvv2vs::persistentInput(PersistentIStream & is, int) {
+ is >> scalar_ >> vector_;
+}
+
+ClassDescription<MEvv2vs> MEvv2vs::initMEvv2vs;
+// Definition of the static class description member.
+
+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());
+ flowME().resize(numberOfFlows(),
+ ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
+ PDT::Spin1, PDT::Spin0));
+ diagramME().resize(numberOfDiags(),
+ ProductionMatrixElement(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) {
+ 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);
+ }
+ }
+}
+
+void MEvv2vs::doinitrun() {
+ GeneralHardME::doinitrun();
+ flowME().resize(numberOfFlows(),
+ ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
+ PDT::Spin1, PDT::Spin0));
+ diagramME().resize(numberOfDiags(),
+ ProductionMatrixElement(PDT::Spin1, PDT::Spin1,
+ PDT::Spin1, PDT::Spin0));
+}
+
+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);
+ vb[i] = VectorWaveFunction(rescaledMomenta()[1], mePartonData()[1], 2*i,
+ 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);
+ }
+ 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) {
+ 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, 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
+ 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)
+ 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*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
new file mode 100644
--- /dev/null
+++ b/MatrixElement/General/MEvv2vs.h
@@ -0,0 +1,204 @@
+// -*- 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 "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();
+
+ /**
+ * Initialize this object. Called in the run phase just before
+ * a run begins.
+ */
+ virtual void doinitrun();
+ //@}
+
+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_;
+
+};
+
+}
+
+#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 */
diff --git a/MatrixElement/General/Makefile.am b/MatrixElement/General/Makefile.am
--- a/MatrixElement/General/Makefile.am
+++ b/MatrixElement/General/Makefile.am
@@ -1,19 +1,20 @@
noinst_LTLIBRARIES = libHwGeneralME.la
libHwGeneralME_la_SOURCES = \
GeneralHardME.cc GeneralHardME.h GeneralHardME.fh \
MEvv2ff.cc MEvv2ff.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/ResonantProcessConstructor.cc b/Models/General/ResonantProcessConstructor.cc
--- a/Models/General/ResonantProcessConstructor.cc
+++ b/Models/General/ResonantProcessConstructor.cc
@@ -1,309 +1,366 @@
// -*- C++ -*-
//
// ResonantProcessConstructor.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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 ResonantProcessConstructor class.
//
#include "ResonantProcessConstructor.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
using namespace Herwig;
IBPtr ResonantProcessConstructor::clone() const {
return new_ptr(*this);
}
IBPtr ResonantProcessConstructor::fullclone() const {
return new_ptr(*this);
}
void ResonantProcessConstructor::persistentOutput(PersistentOStream & os) const {
- os << theIncoming << theIntermediates << theOutgoing;
+ os << incoming_ << intermediates_ << outgoing_ << processOption_;
}
void ResonantProcessConstructor::persistentInput(PersistentIStream & is, int) {
- is >> theIncoming >> theIntermediates >> theOutgoing;
+ is >> incoming_ >> intermediates_ >> outgoing_ >> processOption_;
}
ClassDescription<ResonantProcessConstructor>
ResonantProcessConstructor::initResonantProcessConstructor;
// Definition of the static class description member.
void ResonantProcessConstructor::Init() {
static ClassDocumentation<ResonantProcessConstructor> documentation
("This class is designed solely to contruct resonant processes using"
"a provided set of intermediate particles");
static RefVector<ResonantProcessConstructor, ParticleData> interfaceOffshell
("Intermediates",
"A vector of offshell particles for resonant diagrams",
- &ResonantProcessConstructor::theIntermediates, -1, false, false, true,
+ &ResonantProcessConstructor::intermediates_, -1, false, false, true,
false);
static RefVector<ResonantProcessConstructor, ParticleData> interfaceIncoming
("Incoming",
"A vector of incoming particles for resonant diagrams",
- &ResonantProcessConstructor::theIncoming, -1, false, false, true,
+ &ResonantProcessConstructor::incoming_, -1, false, false, true,
false);
static RefVector<ResonantProcessConstructor, ParticleData> interfaceOutgoing
("Outgoing",
"A vector of outgoin particles for resonant diagrams",
- &ResonantProcessConstructor::theOutgoing, -1, false, false, true,
+ &ResonantProcessConstructor::outgoing_, -1, false, false, true,
false);
+
+ static Switch<ResonantProcessConstructor,unsigned int> interfaceProcesses
+ ("Processes",
+ "Whether to generate inclusive or exclusive processes",
+ &ResonantProcessConstructor::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 interfaceProcessesInclusive
+ (interfaceProcesses,
+ "Inclusive",
+ "Generate all modes which are allowed for the on-shell intermediate particle",
+ 3);
+}
+
+void ResonantProcessConstructor::doinit() {
+ HardProcessConstructor::doinit();
+ if(processOption_==2&&outgoing_.size()!=2)
+ throw InitException()
+ << "Exclusive processes require exactly"
+ << " two outgoing particles but " << outgoing_.size()
+ << "have been inserted in ResonantProcessConstructor::doinit()."
+ << Exception::runerror;
}
namespace {
// Helper functor for find_if in duplicate function.
class SameIncomingAs {
public:
SameIncomingAs(tPDPair in) : a(in.first->id()), b(in.second->id()) {}
bool operator()(tPDPair ppair) const {
long id1(ppair.first->id()), id2(ppair.second->id());
return ( id1 == a && id2 == b ) || ( id1 == b && id2 == a );
}
private:
long a, b;
};
bool duplicateIncoming(tPDPair ppair,const vector<tPDPair> &incPairs) {
vector<tPDPair>::const_iterator it =
find_if( incPairs.begin(), incPairs.end(), SameIncomingAs(ppair) );
return it != incPairs.end();
}
}
void ResonantProcessConstructor::constructDiagrams() {
- size_t ninc = theIncoming.size() , ninter = theIntermediates.size();
+ size_t ninc = incoming_.size() , ninter = intermediates_.size();
if(ninc == 0 || ninter == 0 ) return;
// find the incoming particle pairs
vector<tPDPair> incPairs;
for(PDVector::size_type i = 0; i < ninc; ++i) {
for(PDVector::size_type j = 0; j < ninc; ++j) {
- tPDPair inc = make_pair(theIncoming[i], theIncoming[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( !duplicateIncoming(inc,incPairs) ) {
incPairs.push_back(inc);
}
}
}
//Remove matrix elements already present
EGPtr eg = generator();
int nme = subProcess()->MEs().size();
for(int ix = nme - 1; ix >= 0; --ix)
eg->preinitInterface(subProcess(), "MatrixElements", ix, "erase", "");
// Get a pointer to the model in use
model()->init();
size_t nvertices = model()->numberOfVertices();
//init the vertices
for(size_t iv = 0; iv < nvertices; ++iv)
model()->vertex(iv)->init();
//To construct resonant diagrams loop over the incoming particles, intermediates
//and vertices to find allowed diagrams. Need to exclude the diagrams that have
//the intermediate as an external particle as well
for(vector<tcPDPair>::size_type is = 0; is < incPairs.size(); ++is) {
tPDPair ppi = incPairs[is];
for(vector<PDPtr>::size_type ik = 0; ik < ninter ; ++ik) {
- long ipart = theIntermediates[ik]->id();
+ long ipart = intermediates_[ik]->id();
for(size_t iv = 0; iv < nvertices; ++iv) {
VBPtr vertex = model()->vertex(iv);
if(vertex->getNpoint() > 3) continue;
long part1 = ppi.first->CC() ? -ppi.first->id() : ppi.first->id();
long part2 = ppi.second->CC() ? -ppi.second->id() : ppi.second->id();
if(vertex->allowed(part1, part2, ipart) ||
vertex->allowed(part1, ipart, part2) ||
vertex->allowed(part2, part1, ipart) ||
vertex->allowed(part2, ipart, part1) ||
vertex->allowed(ipart, part1, part2) ||
vertex->allowed(ipart, part2, part1) ) {
constructVertex2(make_pair(ppi.first->id(), ppi.second->id()), vertex,
- theIntermediates[ik]);
+ intermediates_[ik]);
}
}
}
}
//Create matrix element for each process
- const HPDVector::size_type ndiags = theDiagrams.size();
+ const HPDVector::size_type ndiags = diagrams_.size();
for(HPDVector::size_type ix = 0; ix < ndiags; ++ix)
- createMatrixElement(theDiagrams[ix]);
+ createMatrixElement(diagrams_[ix]);
}
void ResonantProcessConstructor::
constructVertex2(IDPair in, VertexBasePtr vertex,
PDPtr partc) {
//We have the left hand part of the diagram, just need all the possibilities
//for the RHS
size_t nvertices = model()->numberOfVertices();
- if(!theOutgoing.empty()) {
- for(size_t io = 0; io < theOutgoing.size(); ++io) {
- tcPDPtr outa = theOutgoing[io];
+ if(processOption_!=3) {
+ for(size_t io = 0; io < outgoing_.size(); ++io) {
+ tcPDPtr outa = outgoing_[io];
for(size_t iv = 0; iv < nvertices; ++iv) {
VBPtr vertex2 = model()->vertex(iv);
if(vertex2->getNpoint() > 3) continue;
tPDSet outb = search(vertex2, partc->id(), incoming, outa->id(), outgoing,
outgoing);
for(tPDSet::const_iterator ita = outb.begin(); ita != outb.end(); ++ita)
makeResonantDiagram(in, partc, outa->id(),(**ita).id(),
make_pair(vertex, vertex2));
}
}
}
else {
for(size_t iv = 0; iv < nvertices; ++iv) {
VBPtr vertex2 = model()->vertex(iv);
if(vertex2->getNpoint() > 3) continue;
for(unsigned int ix = 0;ix < 3; ++ix) {
vector<long> pdlist = vertex2->search(ix, partc->id());
for(unsigned int iy=0;iy<pdlist.size();iy+=3) {
long out1 = ix==0 ? pdlist.at(iy+1) : pdlist.at(iy );
long out2 = ix==2 ? pdlist.at(iy+1) : pdlist.at(iy+2);
if(partc->mass() < getParticleData(out1)->mass() +
getParticleData(out2)->mass()) continue;
makeResonantDiagram(in, partc, out1, out2,
make_pair(vertex, vertex2));
}
}
}
}
}
void ResonantProcessConstructor::
makeResonantDiagram(IDPair in, PDPtr offshell, long outa, long outb,
VBPair vertpair) {
assert(vertpair.first && vertpair.second);
if( abs(outa) == abs(offshell->id()) ||
abs(outb) == abs(offshell->id())) return;
HPDiagram newdiag(in,make_pair(outa,outb));
newdiag.intermediate = offshell;
newdiag.vertices = vertpair;
newdiag.channelType = HPDiagram::sChannel;
fixFSOrder(newdiag);
assignToCF(newdiag);
- if(!duplicate(newdiag,theDiagrams))
- theDiagrams.push_back(newdiag);
+ if(duplicate(newdiag,diagrams_)) return;
+ // if inclusive enforce the exclusivity
+ if(processOption_==1) {
+ PDVector::const_iterator loc =
+ std::find(outgoing_.begin(),outgoing_.end(),
+ getParticleData(newdiag.outgoing. first));
+ if(loc==outgoing_.end()) return;
+ loc =
+ std::find(outgoing_.begin(),outgoing_.end(),
+ getParticleData(newdiag.outgoing.second));
+ if(loc==outgoing_.end()) return;
+ }
+ else if(processOption_==2) {
+ if(!((newdiag.outgoing. first==outgoing_[0]->id()&&
+ newdiag.outgoing.second==outgoing_[1]->id())||
+ (newdiag.outgoing. first==outgoing_[1]->id()&&
+ newdiag.outgoing.second==outgoing_[0]->id())))
+ return;
+ }
+ // add to the list
+ diagrams_.push_back(newdiag);
}
set<tPDPtr>
ResonantProcessConstructor::search(VBPtr vertex, long part1, direction d1,
long part2, direction d2, direction d3) {
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;
}
IDPair ResonantProcessConstructor::
find(long part, const vector<PDPtr> & out) const {
vector<PDPtr>::size_type iloc(0);
bool found(false);
do {
if(out[iloc]->id() == part) found = true;
else ++iloc;
}
while(found == false && iloc < out.size());
//found offshell
IDPair outids;
if(iloc == 0)
outids = make_pair(out[1]->id(), out[2]->id());
else if(iloc == 1)
outids = make_pair(out[0]->id(), out[2]->id());
else
outids = make_pair(out[0]->id(), out[1]->id());
return outids;
}
void ResonantProcessConstructor::
createMatrixElement(const HPDiagram & diag) const {
vector<tcPDPtr> extpart(4);
extpart[0] = getParticleData(diag.incoming.first);
extpart[1] = getParticleData(diag.incoming.second);
extpart[2] = getParticleData(diag.outgoing.first);
extpart[3] = getParticleData(diag.outgoing.second);
string objectname ("/Herwig/MatrixElements/");
string classname = MEClassname(extpart, diag.intermediate, objectname);
GeneralHardMEPtr matrixElement = dynamic_ptr_cast<GeneralHardMEPtr>
(generator()->preinitCreate(classname, objectname));
if( !matrixElement ) {
throw RPConstructorError()
<< "ResonantProcessConstructor::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 << "\"\n"
<< Exception::warning;
return;
}
matrixElement->setProcessInfo(HPDVector(1, diag),
colourFlow(extpart), debug(),0);
generator()->preinitInterface(subProcess(), "MatrixElements",
subProcess()->MEs().size(),
"insert", matrixElement->fullName());
}
string ResonantProcessConstructor::
MEClassname(const vector<tcPDPtr> & extpart, tcPDPtr inter,
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::Spin2) classname += "t";
else
throw RPConstructorError()
<< "MEClassname() : Encountered an unknown spin for "
<< extpart[ix]->PDGName() << " while constructing MatrixElement "
<< "classname " << extpart[ix]->iSpin() << Exception::warning;
}
objname += "ME" + extpart[0]->PDGName() + extpart[1]->PDGName() + "2"
+ inter->PDGName() + "2"
+ extpart[2]->PDGName() + extpart[3]->PDGName();
return classname;
}
diff --git a/Models/General/ResonantProcessConstructor.h b/Models/General/ResonantProcessConstructor.h
--- a/Models/General/ResonantProcessConstructor.h
+++ b/Models/General/ResonantProcessConstructor.h
@@ -1,215 +1,234 @@
// -*- C++ -*-
//
// ResonantProcessConstructor.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ResonantProcessConstructor_H
#define HERWIG_ResonantProcessConstructor_H
//
// This is the declaration of the ResonantProcessConstructor class.
//
#include "HardProcessConstructor.h"
#include "ThePEG/Utilities/Exception.h"
#include "ResonantProcessConstructor.fh"
namespace Herwig {
using namespace ThePEG;
/**
* This class is designed to construct the diagrams for resonant processes
* using a provdided set of particles as interemdiates.
*
* @see \ref ResonantProcessConstructorInterfaces "The interfaces"
* defined for ResonantProcessConstructor.
* @see HardProcessConstructor
*/
class ResonantProcessConstructor: public HardProcessConstructor {
public:
/** Set of ParticleData pointers */
typedef set<tPDPtr> tPDSet;
/** Nested vector of doubles. */
typedef vector<vector<double> > CFMatrix;
/** Enumeration for the direction */
enum direction {incoming, outgoing};
public:
/**
* The default constructor.
*/
ResonantProcessConstructor() :
- theIncoming(0), theIntermediates(0), theOutgoing(0), theDiagrams(0)
+ processOption_(0), incoming_(0), intermediates_(0),
+ outgoing_(0), diagrams_(0)
{}
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();
/**
* The main function to create the resonant diagrams
*/
void constructDiagrams() ;
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 HardProcessConstructor 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:
/**
* Utility function to help second vertex
*/
void constructVertex2(IDPair in, VertexBasePtr vertex,
PDPtr partc);
/**
* Function to create the appropriate diagrams
*/
void makeResonantDiagram(IDPair in, PDPtr offshell, long outa,
long outb, VBPair vertices);
/**
* Given a vertex and 2 particle id's find the possible states
* that can be the 3rd external particle
* @param vertex Pointer to the vertex
* @param part1 id of first particle
* @param d1 direction of particle one
* @param part2 id of other particle
* @param d2 direction of particle two
* @param d3 required direction of 3rd state (default = outgoing)
* @return container of third particles
*/
tPDSet search(VertexBasePtr vertex, long part1, direction d1,
long part2, direction d2, direction d3 = outgoing);
/**
* Return the pair of outgoing particles from the list
*/
IDPair find(long part, const PDVector & out) const;
/**
* Create a matrix element from the given resonant process diagram
*/
void createMatrixElement(const HPDiagram & diag) const;
/**
* Create the correct classname and objectname for a matrix element
*/
string MEClassname(const tcPDVector & extpart, tcPDPtr inter,
string & objname) const;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<ResonantProcessConstructor>
initResonantProcessConstructor;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ResonantProcessConstructor & operator=(const ResonantProcessConstructor &);
private:
+
+ /**
+ * Which types of processes to generate
+ */
+ unsigned int processOption_;
/**
* Storage for the required intermediate particles
*/
- vector<PDPtr> theIncoming;
+ vector<PDPtr> incoming_;
/**
* Storage for the required intermediate particles
*/
- vector<PDPtr> theIntermediates;
+ vector<PDPtr> intermediates_;
/**
* Storage for the required intermediate particles
*/
- vector<PDPtr> theOutgoing;
+ vector<PDPtr> outgoing_;
/**
* Storage for the diagrams
*/
- vector<HPDiagram> theDiagrams;
+ vector<HPDiagram> diagrams_;
+
};
/** Exception class indicating setup problem. */
class RPConstructorError : public Exception {};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ResonantProcessConstructor. */
template <>
struct BaseClassTrait<Herwig::ResonantProcessConstructor,1> {
/** Typedef of the first base class of ResonantProcessConstructor. */
typedef Herwig::HardProcessConstructor NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ResonantProcessConstructor class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ResonantProcessConstructor>
: public ClassTraitsBase<Herwig::ResonantProcessConstructor> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ResonantProcessConstructor"; }
};
/** @endcond */
}
#endif /* HERWIG_ResonantProcessConstructor_H */
diff --git a/Models/General/TwoToTwoProcessConstructor.cc b/Models/General/TwoToTwoProcessConstructor.cc
--- a/Models/General/TwoToTwoProcessConstructor.cc
+++ b/Models/General/TwoToTwoProcessConstructor.cc
@@ -1,619 +1,620 @@
// -*- C++ -*-
//
// TwoToTwoProcessConstructor.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 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/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/Switch.h"
using namespace Herwig;
namespace {
// Helper functor for find_if in duplicate function.
class SameIncomingAs {
public:
SameIncomingAs(tPDPair in) : a(in.first->id()), b(in.second->id()) {}
bool operator()(tPDPair ppair) const {
long id1(ppair.first->id()), id2(ppair.second->id());
return ( id1 == a && id2 == b ) || ( id1 == b && id2 == a );
}
private:
long a, b;
};
bool duplicateIncoming(tPDPair ppair, const vector<tPDPair> & incPairs ) {
vector<tPDPair>::const_iterator it =
find_if( incPairs.begin(), incPairs.end(), SameIncomingAs(ppair) );
return it != incPairs.end();
}
}
TwoToTwoProcessConstructor::TwoToTwoProcessConstructor() :
Nout_(0), nv_(0), allDiagrams_(true),
processOption_(0), scaleChoice_(0)
{}
IBPtr TwoToTwoProcessConstructor::clone() const {
return new_ptr(*this);
}
IBPtr TwoToTwoProcessConstructor::fullclone() const {
return new_ptr(*this);
}
void TwoToTwoProcessConstructor::doinit() {
HardProcessConstructor::doinit();
if(processOption_==2&&outgoing_.size()!=2)
throw InitException()
<< "Exclusive processes require exactly"
<< " two outgoing particles but " << outgoing_.size()
- << "have been inserted." << Exception::runerror;
+ << "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( !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_ << excluded_ << excludedExternal_
<< excludedVertexVector_ << excludedVertexSet_;
}
void TwoToTwoProcessConstructor::persistentInput(PersistentIStream & is, int) {
is >> vertices_ >> incoming_ >> outgoing_
>> allDiagrams_ >> processOption_
>> scaleChoice_ >> excluded_ >> excludedExternal_
>> excludedVertexVector_ >> excludedVertexSet_;
}
ClassDescription<TwoToTwoProcessConstructor>
TwoToTwoProcessConstructor::initTwoToTwoProcessConstructor;
// Definition of the static class description member.
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 interfaceIncludeAllDiagramsOff
(interfaceIncludeAllDiagrams,
"No",
"Only include QCD diagrams",
false);
static SwitchOption interfaceIncludeAllDiagramsOn
(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 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 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);
}
namespace {
// Helper functor for find_if below.
class SameProcessAs {
public:
SameProcessAs(const HPDiagram & diag) : a(diag) {}
bool operator()(const HPDiagram & b) const {
return a.sameProcess(b);
}
private:
HPDiagram a;
};
}
void TwoToTwoProcessConstructor::constructDiagrams() {
if(incPairs_.empty() || outgoing_.empty()) return;
// delete the matrix elements we already have
int nme = subProcess()->MEs().size();
for(int ix = nme - 1; ix >= 0; --ix)
generator()->preinitInterface(subProcess(), "MatrixElements",
ix, "erase", "");
model()->init();
nv_ = model()->numberOfVertices();
//make sure vertices are initialised
for(unsigned int ix = 0; ix < nv_; ++ix ) {
VertexBasePtr vertex = model()->vertex(ix);
vertex->init();
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
multimap<HPDiagram,HPDiagram > grouped;
HPDVector::iterator dit = processes_.begin();
HPDVector::iterator dend = processes_.end();
for( ; dit != dend; ++dit) {
grouped.insert(make_pair(*dit, *dit));
}
assert( processes_.size() == grouped.size() );
processes_.clear();
typedef multimap<HPDiagram, HPDiagram>::const_iterator map_iter;
map_iter it = grouped.begin();
map_iter iend = grouped.end();
while( it != iend ) {
pair< map_iter, map_iter> range = grouped.equal_range(it->first);
map_iter itb = range.first;
HPDVector process;
for( ; itb != range.second; ++itb ) {
process.push_back(itb->second);
}
// if inclusive enforce the exclusivity
if(processOption_==2) {
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(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;
}
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( !duplicate(nhp, processes_) ) {
assignToCF(nhp);
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( !duplicate(nhp, processes_) ) {
assignToCF(nhp);
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 ) {
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);
// 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::Spin2) classname += "t";
else {
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;
}
diff --git a/Models/StandardModel/StandardModel.h b/Models/StandardModel/StandardModel.h
--- a/Models/StandardModel/StandardModel.h
+++ b/Models/StandardModel/StandardModel.h
@@ -1,473 +1,473 @@
// -*- C++ -*-
//
// StandardModel.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_StandardModel_H
#define HERWIG_StandardModel_H
//
// This is the declaration of the StandardModel class.
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "Herwig++/Models/StandardModel/RunningMassBase.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractSSSSVertex.h"
#include "Herwig++/Models/General/ModelGenerator.fh"
#include "StandardModel.fh"
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
/** \ingroup Models
*
* This is the Herwig++ StandardModel class which inherits from ThePEG
* Standard Model class and implements additional Standard Model couplings,
* access to vertices for helicity amplitude calculations etc.
*
* @see StandardModelBase
*/
class StandardModel: public StandardModelBase {
/**
* Some typedefs for the pointers.
*/
//@{
/**
* Pointer to the RunningMassBase object
*/
typedef Ptr<Herwig::RunningMassBase>::pointer runPtr;
/**
* Transient pointer to the RunningMassBase object
*/
typedef Ptr<Herwig::RunningMassBase>::transient_pointer trunPtr;
//@}
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor
*/
StandardModel();
/**
* Copy-constructor.
*/
StandardModel(const StandardModel &);
/**
* Destructor
*/
virtual ~StandardModel();
//@}
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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
public:
/**
* The left and right couplings of the Z^0 including sin and cos theta_W.
*/
//@{
/**
* The left-handed coupling of a neutrino
*/
double lnu() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(vnu()+anu());
}
/**
* The left-handed coupling of a charged lepton.
*/
double le() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(ve()+ae());
}
/**
* The left-handed coupling of an up type quark.
*/
double lu() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(vu()+au());
}
/**
* The left-handed coupling of a down type quark.
*/
double ld() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(vd()+ad());
}
/**
* The right-handed coupling of a neutrino
*/
double rnu() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(vnu()-anu());
}
/**
* The right-handed coupling of a charged lepton.
*/
double re() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(ve()-ae());
}
/**
* The right-handed coupling of an up type quark.
*/
double ru() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(vu()-au());
}
/**
* The right-handed coupling of a down type quark.
*/
double rd() const {
return 0.25/sqrt(sin2ThetaW()*(1.-sin2ThetaW()))*(vd()-ad());
}
//@}
/**
* Pointers to the objects handling the vertices.
*/
//@{
/**
* Pointer to the fermion-fermion-Z vertex
*/
tAbstractFFVVertexPtr vertexFFZ() const {
return FFZVertex_;
}
/**
* Pointer to the fermion-fermion-photon vertex
*/
tAbstractFFVVertexPtr vertexFFP() const {
return FFPVertex_;
}
/**
* Pointer to the fermion-fermion-gluon vertex
*/
tAbstractFFVVertexPtr vertexFFG() const {
return FFGVertex_;
}
/**
* Pointer to the fermion-fermion-W vertex
*/
tAbstractFFVVertexPtr vertexFFW() const {
return FFWVertex_;
}
/**
* Pointer to the fermion-fermion-Higgs vertex
*/
virtual tAbstractFFSVertexPtr vertexFFH() const {
return FFHVertex_;
}
/**
* Pointer to the triple gluon vertex
*/
tAbstractVVVVertexPtr vertexGGG() const {
return GGGVertex_;
}
/**
* Pointer to the triple electroweak gauge boson vertex.
*/
tAbstractVVVVertexPtr vertexWWW() const {
return WWWVertex_;
}
/**
* Pointer to the two electroweak gauge boson Higgs vertex.
*/
virtual tAbstractVVSVertexPtr vertexWWH() const {
return WWHVertex_;
}
/**
* Pointer to the quartic electroweak gauge boson vertex.
*/
tAbstractVVVVVertexPtr vertexWWWW() const {
return WWWWVertex_;
}
/**
* Pointer to the quartic gluon vertex
*/
tAbstractVVVVVertexPtr vertexGGGG() const {
return GGGGVertex_;
}
/**
* Pointer to the quartic gluon vertex
*/
virtual tAbstractVVSVertexPtr vertexHGG() const {
return HGGVertex_;
}
/**
* Pointer to the quartic gluon vertex
*/
tAbstractVVSVertexPtr vertexHPP() const {
return HPPVertex_;
}
/**
* Pointer to the triple Higgs vertex
*/
tAbstractSSSVertexPtr vertexHHH() const {
return HHHVertex_;
}
/**
* Pointer to the WWHH vertex
*/
tAbstractVVSSVertexPtr vertexWWHH() const {
return WWHHVertex_;
}
/**
* Total number of vertices
*/
unsigned int numberOfVertices() const {
return vertexList_.size();
}
/**
* Access to a vertex from the list
*/
- tVertexBasePtr vertex(unsigned int ix) {
+ tVertexBasePtr vertex(unsigned int ix) const {
return vertexList_[ix];
}
//@}
/**
* Return the running mass for a given scale \f$q^2\f$ and particle type.
* @param scale The scale \f$q^2\f$.
* @param part The ParticleData object for the particle
*/
Energy mass(Energy2 scale,tcPDPtr part) const {
return runningMass_->value(scale,part);
}
/**
* Return a pointer to the object handling the running mass.
*/
trunPtr massPtr() const {
return runningMass_;
}
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:
/**
* Initialize this object after the setup phase before saving and
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* Add a vertex to the list
*/
void addVertex(VertexBasePtr in) {
vertexList_.push_back(in);
}
private:
/**
* Describe a concrete class with persistent data.
*/
static ClassDescription<StandardModel> initStandardModel;
/**
* Private and non-existent assignment operator.
*/
StandardModel & operator=(const StandardModel &);
private:
/**
* Pointers to the vertices for Standard Model helicity amplitude
* calculations.
*/
//@{
/**
* Pointer to the fermion-fermion-Z vertex
*/
AbstractFFVVertexPtr FFZVertex_;
/**
* Pointer to the fermion-fermion-photon vertex
*/
AbstractFFVVertexPtr FFPVertex_;
/**
* Pointer to the fermion-fermion-gluon vertex
*/
AbstractFFVVertexPtr FFGVertex_;
/**
* Pointer to the fermion-fermion-W vertex
*/
AbstractFFVVertexPtr FFWVertex_;
/**
* Pointer to the fermion-fermion-Higgs vertex
*/
AbstractFFSVertexPtr FFHVertex_;
/**
* Pointer to the two electroweak gauge boson Higgs vertex.
*/
AbstractVVSVertexPtr WWHVertex_;
/**
* Pointer to the triple gluon vertex
*/
AbstractVVVVertexPtr GGGVertex_;
/**
* Pointer to the triple electroweak gauge boson vertex.
*/
AbstractVVVVertexPtr WWWVertex_;
/**
* Pointer to the quartic gluon vertex
*/
AbstractVVVVVertexPtr GGGGVertex_;
/**
* Pointer to the quartic electroweak gauge boson vertex.
*/
AbstractVVVVVertexPtr WWWWVertex_;
/**
* Pointer to higgs-gluon-gluon vertex
*/
AbstractVVSVertexPtr HGGVertex_;
/**
* Pointer to higgs-gamma-gamma vertex
*/
AbstractVVSVertexPtr HPPVertex_;
/**
* Pointer to triple Higgs vertex
*/
AbstractSSSVertexPtr HHHVertex_;
/**
* Pointer to WWHH vertex
*/
AbstractVVSSVertexPtr WWHHVertex_;
/**
* Full list of vertices as a vector to allow searching
*/
vector<VertexBasePtr> vertexList_;
//@}
/**
* The running mass.
*/
runPtr runningMass_;
/**
* Pointer to ModelGenerator Class
*/
ModelGeneratorPtr modelGenerator_;
};
}
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* The following template specialization informs ThePEG about the
* base class of StandardModel.
*/
template <>
struct BaseClassTrait<Herwig::StandardModel,1> {
/** Typedef of the base class of StandardModel. */
typedef StandardModelBase NthBase;
};
/**
* The following template specialization informs ThePEG about the
* name of this class and the shared object where it is defined.
*/
template <>
struct ClassTraits<Herwig::StandardModel>
: public ClassTraitsBase<Herwig::StandardModel> {
/**
* Return the class name.
*/
static string className() { return "Herwig::StandardModel"; }
};
/** @endcond */
}
#endif /* HERWIG_StandardModel_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:28 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983045
Default Alt Text
(105 KB)

Event Timeline