Page MenuHomeHEPForge
Files F8310524

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,325 +1,346 @@
// -*- C++ -*-
//
// MEff2ss.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEff2ss class.
//
#include "MEff2ss.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using ThePEG::Helicity::VectorWaveFunction;
using ThePEG::Helicity::TensorWaveFunction;
using ThePEG::Helicity::incoming;
using ThePEG::Helicity::outgoing;
void MEff2ss::doinit() {
GeneralHardME::doinit();
fermion_.resize(numberOfDiags());
+ RSfermion_.resize(numberOfDiags());
scalar_ .resize(numberOfDiags());
vector_ .resize(numberOfDiags());
tensor_ .resize(numberOfDiags());
fourPoint_.resize(numberOfDiags());
initializeMatrixElements(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 if(current.intermediate->iSpin() == PDT::Spin3Half)
+ RSfermion_[i] =
+ make_pair(dynamic_ptr_cast<AbstractRFSVertexPtr>(current.vertices.first),
+ dynamic_ptr_cast<AbstractRFSVertexPtr>(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::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 if(current.channelType == HPDiagram::fourPoint) {
fourPoint_[i] = dynamic_ptr_cast<AbstractFFSSVertexPtr>(current.vertices.first);
}
else
throw InitException() << "MEFF2ss:doinit() - Cannot find correct "
<< "channel from diagram. Vertex not cast! "
<< Exception::runerror;
}
}
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.channelType == HPDiagram::tChannel) {
if(internal->CC()) internal=internal->CC();
- unsigned int iopt = ( abs(sbar[if2].particle()->id()) == abs(internal->id()) ||
- abs(sp[if1] .particle()->id()) == abs(internal->id())) ? 5 : 3;
- SpinorBarWaveFunction interFB;
- if(current.ordered.second) {
- if(iopt==3) {
- interFB = fermion_[ix].second->
- evaluate(q2, iopt, internal, sbar[if2], sca2);
+ if(internal->iSpin() == PDT::Spin1Half) {
+ unsigned int iopt = ( abs(sbar[if2].particle()->id()) == abs(internal->id()) ||
+ abs(sp[if1] .particle()->id()) == abs(internal->id())) ? 5 : 3;
+ SpinorBarWaveFunction interFB;
+ if(current.ordered.second) {
+ if(iopt==3) {
+ interFB = fermion_[ix].second->
+ evaluate(q2, iopt, internal, sbar[if2], sca2);
+ }
+ else {
+ interFB = fermion_[ix].second->
+ evaluate(q2, iopt, internal, sbar[if2], sca2, 0.*GeV, 0.*GeV);
+ }
+ diag = fermion_[ix].first->evaluate(q2, sp[if1], interFB, sca1);
}
else {
- interFB = fermion_[ix].second->
- evaluate(q2, iopt, internal, sbar[if2], sca2, 0.*GeV, 0.*GeV);
+ if(iopt==3) {
+ interFB = fermion_[ix].second->
+ evaluate(q2, iopt, internal, sbar[if2], sca1);
+ }
+ else {
+ interFB = fermion_[ix].second->
+ evaluate(q2, iopt, internal, sbar[if2], sca1, 0.*GeV, 0.*GeV);
+ }
+ diag = fermion_[ix].first->evaluate(q2, sp[if1], interFB, sca2);
}
- diag = fermion_[ix].first->evaluate(q2, sp[if1], interFB, sca1);
}
- else {
- if(iopt==3) {
- interFB = fermion_[ix].second->
- evaluate(q2, iopt, internal, sbar[if2], sca1);
+ else if (internal->iSpin() == PDT::Spin3Half) {
+ RSSpinorBarWaveFunction interFB;
+ if(current.ordered.second) {
+ interFB = RSfermion_[ix].second->
+ evaluate(q2, 3, internal, sbar[if2], sca2);
+ diag = RSfermion_[ix].first->evaluate(q2, sp[if1], interFB, sca1);
}
- else {
- interFB = fermion_[ix].second->
- evaluate(q2, iopt, internal, sbar[if2], sca1, 0.*GeV, 0.*GeV);
+ else {
+ interFB = RSfermion_[ix].second->
+ evaluate(q2, 3, internal, sbar[if2], sca1);
+ diag = RSfermion_[ix].first->evaluate(q2, sp[if1], interFB, sca2);
}
- diag = fermion_[ix].first->evaluate(q2, sp[if1], interFB, sca2);
}
+ else
+ assert(false);
}
else if(current.channelType == HPDiagram::sChannel) {
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);
}
}
else if(current.channelType == HPDiagram::fourPoint) {
diag = fourPoint_[ix]->evaluate(q2,sp[if1], sbar[if2],sca1,sca2);
}
// 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) {
assert(current.colourFlow[iy].first<flows.size());
flows[current.colourFlow[iy].first] +=
current.colourFlow[iy].second * diag;
}
}
// MEs for the different colour flows
for(unsigned int iy = 0; iy < numberOfFlows(); ++iy)
flowME()[iy](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_ << scalar_ << vector_ << tensor_ << fourPoint_;
+ os << fermion_ << scalar_ << vector_ << tensor_ << fourPoint_ << RSfermion_;
}
void MEff2ss::persistentInput(PersistentIStream & is, int) {
- is >> fermion_ >> scalar_ >> vector_ >> tensor_ >> fourPoint_;
+ is >> fermion_ >> scalar_ >> vector_ >> tensor_ >> fourPoint_ >> RSfermion_;
initializeMatrixElements(PDT::Spin1Half, PDT::Spin1Half,
PDT::Spin0 , PDT::Spin0 );
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEff2ss,GeneralHardME>
describeHerwigMEff2ss("Herwig::MEff2ss", "Herwig.so");
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,201 +1,201 @@
// -*- C++ -*-
//
// MEff2ss.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_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/AbstractRFSVertex.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/Vertex/AbstractFFSSVertex.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), fourPoint_(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();
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
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
+ * vector
+ */
+ vector<pair<AbstractFFSVertexPtr, AbstractSSSVertexPtr> > scalar_;
+
+ /**
+ * 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_;
+ vector<pair<AbstractFFVVertexPtr, AbstractVSSVertexPtr> > vector_;
/**
* Storage for dynamically cast vertices for a diagram with intermediate
- * vector
+ * fermion
*/
- vector<pair<AbstractFFVVertexPtr, AbstractVSSVertexPtr> > vector_;
+ vector<pair<AbstractRFSVertexPtr, AbstractRFSVertexPtr> > RSfermion_;
/**
* Storage for dynamically cast vertices for a diagram with intermediate
* tensor
*/
vector<pair<AbstractFFTVertexPtr, AbstractSSTVertexPtr> > tensor_;
/**
* Storage for dynamically cast 4 point vertices
*/
vector<AbstractFFSSVertexPtr> fourPoint_;
};
}
#endif /* HERWIG_MEff2ss_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:44 PM (5 h, 37 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023857
Default Alt Text
(20 KB)

Event Timeline