Page MenuHomeHEPForge

No OneTemporary

diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
@@ -1,554 +1,560 @@
// -*- C++ -*-
//
// MatchboxAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 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_MatchboxAmplitude_H
#define HERWIG_MatchboxAmplitude_H
//
// This is the declaration of the MatchboxAmplitude class.
//
#include "ThePEG/MatrixElement/Amplitude.h"
#include "ThePEG/Handlers/LastXCombInfo.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/ColourBasis.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Process information with coupling order
*/
struct Process {
PDVector legs;
unsigned int orderInAlphaS;
unsigned int orderInAlphaEW;
Process()
: orderInAlphaS(0), orderInAlphaEW(0) {}
Process(const PDVector& p,
unsigned int oas,
unsigned int oae)
: legs(p), orderInAlphaS(oas), orderInAlphaEW(oae) {}
bool operator==(const Process& other) const {
return
legs == other.legs &&
orderInAlphaS == other.orderInAlphaS &&
orderInAlphaEW == other.orderInAlphaEW;
}
bool operator<(const Process& other) const {
if ( orderInAlphaS != other.orderInAlphaS )
return orderInAlphaS < other.orderInAlphaS;
if ( orderInAlphaEW != other.orderInAlphaEW )
return orderInAlphaEW < other.orderInAlphaEW;
return legs < other.legs;
}
void persistentOutput(PersistentOStream & os) const {
os << legs << orderInAlphaS << orderInAlphaEW;
}
void persistentInput(PersistentIStream & is) {
is >> legs >> orderInAlphaS >> orderInAlphaEW;
}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Enumerate the type of calculation required
*/
namespace ProcessType {
enum Types {
treeME2 = 0,
colourCorrelatedME2,
spinColourCorrelatedME2,
oneLoopInterference
};
}
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxAmplitude is the base class for amplitude
* implementations inside Matchbox.
*
* @see \ref MatchboxAmplitudeInterfaces "The interfaces"
* defined for MatchboxAmplitude.
*/
class MatchboxAmplitude:
public Amplitude,
public LastXCombInfo<StandardXComb>,
public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxAmplitude();
/**
* The destructor.
*/
virtual ~MatchboxAmplitude();
//@}
public:
/**
* Return the amplitude. Needs to be implemented from
* ThePEG::Amplitude but is actually ill-defined, as colours of the
* external particles are not specified. To this extent, this
* implementation just asserts.
*/
virtual Complex value(const tcPDVector & particles,
const vector<Lorentz5Momentum> & momenta,
const vector<int> & helicities);
/** @name Subprocess information */
//@{
/**
* Return true, if this amplitude can handle the given process.
*/
virtual bool canHandle(const PDVector& p,
Ptr<MatchboxFactory>::tptr) const { return canHandle(p); }
/**
* Return true, if this amplitude can handle the given process.
*/
virtual bool canHandle(const PDVector&) const { return false; }
/**
* Return the number of random numbers required to evaluate this
* amplitude at a fixed phase space point.
*/
virtual int nDimAdditional() const { return 0; }
/**
* Return a ME instance appropriate for this amplitude and the given
* subprocesses
*/
virtual Ptr<MatchboxMEBase>::ptr makeME(const PDVector&) const;
/**
* Set the (tree-level) order in \f$g_S\f$ in which this matrix
* element should be evaluated.
*/
virtual void orderInGs(unsigned int) {}
/**
* Return the (tree-level) order in \f$g_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGs() const = 0;
/**
* Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element should be evaluated.
*/
virtual void orderInGem(unsigned int) {}
/**
* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGem() const = 0;
/**
* Return the Herwig++ StandardModel object
*/
Ptr<StandardModel>::tcptr standardModel() {
if ( !hwStandardModel() )
hwStandardModel(dynamic_ptr_cast<Ptr<StandardModel>::tcptr>(HandlerBase::standardModel()));
return hwStandardModel();
}
/**
* Tell whether the outgoing partons should be sorted when determining
* allowed subprocesses. Otherwise, all permutations are counted as
* separate subprocesses.
*/
virtual bool sortOutgoing() { return true; }
/**
* Return true, if this amplitude already includes averaging over
* incoming parton's quantum numbers.
*/
virtual bool hasInitialAverage() const { return false; }
/**
* Return true, if this amplitude already includes symmetry factors
* for identical outgoing particles.
*/
virtual bool hasFinalStateSymmetry() const { return false; }
/**
* Return true, if this amplitude is handled by a BLHA one-loop provider
*/
virtual bool isOLPTree() const { return false; }
/**
* Return true, if this amplitude is handled by a BLHA one-loop provider
*/
virtual bool isOLPLoop() const { return false; }
/**
* Return true, if colour and spin correlated matrix elements should
* be ordered from the OLP
*/
virtual bool needsOLPCorrelators() const { return true; }
/**
* Write the order file header
*/
virtual void olpOrderFileHeader(ostream&) const;
/**
* Write the order file process list
*/
virtual void olpOrderFileProcesses(ostream&,
const map<pair<Process,int>,int>& procs) const;
/**
* Start the one loop provider, if appropriate, giving order and
* contract files
*/
virtual void signOLP(const string&, const string&) { }
/**
* Start the one loop provider, if appropriate
*/
virtual void startOLP(const string&, int& status) { status = -1; }
/**
* Start the one loop provider, if appropriate. This default
* implementation writes an BLHA 2.0 order file and starts the OLP
*/
virtual bool startOLP(const map<pair<Process,int>,int>& procs);
//@}
/** @name Colour basis. */
//@{
/**
* Return the colour basis.
*/
virtual Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }
/**
* Return true, if this amplitude will not require colour correlations.
*/
virtual bool noCorrelations() const { return !haveOneLoop(); }
/**
* Return true, if the colour basis is capable of assigning colour
* flows.
*/
virtual bool haveColourFlows() const {
return colourBasis() ? colourBasis()->haveColourFlows() : false;
}
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
*/
virtual Selector<const ColourLines *> colourGeometries(tcDiagPtr diag) const;
/**
* Return an ordering identifier for the current subprocess and
* colour absis tensor index.
*/
const string& colourOrderingString(size_t id) const;
/**
* Return an ordering identifier for the current subprocess and
* colour absis tensor index.
*/
const vector<vector<size_t> >& colourOrdering(size_t id) const;
//@}
/** @name Phasespace point, crossing and helicities */
//@{
/**
* Set the xcomb object.
*/
virtual void setXComb(tStdXCombPtr xc);
/**
* Return the momentum as crossed appropriate for this amplitude.
*/
Lorentz5Momentum amplitudeMomentum(int) const;
/**
* Perform a normal ordering of external legs and fill the
* crossing information as. This default implementation sorts
* lexicographically in (abs(colour)/spin/abs(charge)), putting pairs
* of particles/anti-particles where possible.
*/
virtual void fillCrossingMap(size_t shift = 0);
/**
* Generate the helicity combinations.
*/
virtual set<vector<int> > generateHelicities() const;
//@}
/** @name Tree-level amplitudes */
//@{
/**
* Calculate the tree level amplitudes for the phasespace point
* stored in lastXComb.
*/
virtual void prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr);
/**
* Return the matrix element squared.
*/
virtual double me2() const;
/**
* Return the colour correlated matrix element.
*/
virtual double colourCorrelatedME2(pair<int,int> ij) const;
/**
* Return the large-N colour correlated matrix element.
*/
virtual double largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr largeNBasis) const;
/**
* Return a positive helicity polarization vector for a gluon of
* momentum p (with reference vector n) to be used when evaluating
* spin correlations.
*/
virtual LorentzVector<Complex> plusPolarization(const Lorentz5Momentum& p,
const Lorentz5Momentum& n,
int id = -1) const;
/**
* Return the colour and spin correlated matrix element.
*/
virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& c) const;
/**
* Return true, if tree-level contributions will be evaluated at amplitude level.
*/
virtual bool treeAmplitudes() const { return true; }
/**
* Evaluate the amplitude for the given colour tensor id and
* helicity assignment
*/
virtual Complex evaluate(size_t, const vector<int>&, Complex&) { return 0.; }
//@}
/** @name One-loop amplitudes */
//@{
/**
* Return true, if this amplitude is capable of calculating one-loop
* (QCD) corrections.
*/
virtual bool haveOneLoop() const { return false; }
/**
* Return true, if this amplitude only provides
* one-loop (QCD) corrections.
*/
virtual bool onlyOneLoop() const { return false; }
/**
* Return true, if one-loop contributions will be evaluated at amplitude level.
*/
virtual bool oneLoopAmplitudes() const { return true; }
/**
* Return true, if one loop corrections have been calculated in
* dimensional reduction. Otherwise conventional dimensional
* regularization is assumed. Note that renormalization is always
* assumed to be MSbar.
*/
virtual bool isDR() const { return false; }
/**
+ * Return true, if the amplitude is DRbar renormalized, otherwise
+ * MSbar is assumed.
+ */
+ virtual bool isDRbar() const { return false; }
+
+ /**
* Return true, if one loop corrections are given in the conventions
* of the integrated dipoles.
*/
virtual bool isCS() const { return false; }
/**
* Return true, if one loop corrections are given in the conventions
* of BDK.
*/
virtual bool isBDK() const { return false; }
/**
* Return true, if one loop corrections are given in the conventions
* of everything expanded.
*/
virtual bool isExpanded() const { return false; }
/**
* Return the value of the dimensional regularization
* parameter. Note that renormalization scale dependence is fully
* restored in DipoleIOperator.
*/
virtual Energy2 mu2() const { return 0.*GeV2; }
/**
* If defined, return the coefficient of the pole in epsilon^2
*/
virtual double oneLoopDoublePole() const { return 0.; }
/**
* If defined, return the coefficient of the pole in epsilon
*/
virtual double oneLoopSinglePole() const { return 0.; }
/**
* Calculate the one-loop amplitudes for the phasespace point
* stored in lastXComb, if provided.
*/
virtual void prepareOneLoopAmplitudes(Ptr<MatchboxMEBase>::tcptr);
/**
* Return the one-loop/tree interference.
*/
virtual double oneLoopInterference() const;
/**
* Evaluate the amplitude for the given colour tensor id and
* helicity assignment
*/
virtual Complex evaluateOneLoop(size_t, const vector<int>&) { return 0.; }
//@}
/** @name Caching and helpers to setup amplitude objects. */
//@{
/**
* Flush all cashes.
*/
virtual void flushCaches() {}
/**
* Clone this amplitude.
*/
Ptr<MatchboxAmplitude>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
virtual void cloneDependencies(const std::string& prefix = "");
//@}
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();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* Recursively generate helicities
*/
void doGenerateHelicities(set<vector<int> >& res,
vector<int>& current,
size_t pos) const;
/**
* The colour basis implementation to be used.
*/
Ptr<ColourBasis>::ptr theColourBasis;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxAmplitude & operator=(const MatchboxAmplitude &);
};
inline PersistentOStream& operator<<(PersistentOStream& os,
const Process& h) {
h.persistentOutput(os);
return os;
}
inline PersistentIStream& operator>>(PersistentIStream& is,
Process& h) {
h.persistentInput(is);
return is;
}
}
#endif /* HERWIG_MatchboxAmplitude_H */
diff --git a/MatrixElement/Matchbox/InsertionOperators/DipoleIOperator.cc b/MatrixElement/Matchbox/InsertionOperators/DipoleIOperator.cc
--- a/MatrixElement/Matchbox/InsertionOperators/DipoleIOperator.cc
+++ b/MatrixElement/Matchbox/InsertionOperators/DipoleIOperator.cc
@@ -1,270 +1,270 @@
// -*- C++ -*-
//
// DipoleIOperator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 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 DipoleIOperator class.
//
#include "DipoleIOperator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
using namespace Herwig;
using Constants::pi;
DipoleIOperator::DipoleIOperator()
: MatchboxInsertionOperator(),
CA(-1.0), CF(-1.0),
gammaQuark(-1.0), gammaGluon(-1.0),
betaZero(-1.),
KQuark(-1.0), KGluon(-1.0) {}
DipoleIOperator::~DipoleIOperator() {}
IBPtr DipoleIOperator::clone() const {
return new_ptr(*this);
}
IBPtr DipoleIOperator::fullclone() const {
return new_ptr(*this);
}
void DipoleIOperator::setXComb(tStdXCombPtr xc) {
MatchboxInsertionOperator::setXComb(xc);
if ( CA < 0. ) {
CA = SM().Nc();
CF = (SM().Nc()*SM().Nc()-1.0)/(2.*SM().Nc());
gammaQuark = (3./2.)*CF;
gammaGluon = (11./6.)*CA - (1./3.)*lastBorn()->nLight();
betaZero = gammaGluon;
KQuark = (7./2.-sqr(pi)/6.)*CF;
KGluon = (67./18.-sqr(pi)/6.)*CA-(5./9.)*lastBorn()->nLight();
if ( isDR() ) {
gammaQuark -= CF/2.;
gammaGluon -= CA/6.;
}
}
}
bool DipoleIOperator::apply(const cPDVector& pd) const {
bool first = false;
bool second = false;
for ( cPDVector::const_iterator p = pd.begin();
p != pd.end(); ++p ) {
if ( !first ) {
if ( apply(*p) )
first = true;
} else {
if ( apply(*p) )
second = true;
}
}
return first && second;
}
bool DipoleIOperator::apply(tcPDPtr pd) const {
return
pd->mass() == ZERO &&
(abs(pd->id()) < 6 || pd->id() == ParticleID::g);
}
double DipoleIOperator::me2() const {
double res = 0.;
int idi = 0; int idj = 0;
Energy2 mu2 = lastBorn()->mu2();
for ( cPDVector::const_iterator i = mePartonData().begin();
i != mePartonData().end(); ++i, ++idi ) {
if ( !apply(*i) )
continue;
idj = 0;
for ( cPDVector::const_iterator j = mePartonData().begin();
j != mePartonData().end(); ++j, ++idj ) {
if ( !apply(*j) )
continue;
if ( i == j || lastBorn()->noDipole(idi,idj) )
continue;
double delta = 0.;
double xgammaGluon = gammaGluon;
double xgammaQuark = gammaQuark;
if ( isDR() ) {
xgammaGluon += CA/6.;
xgammaQuark += CF/2.;
}
if ( isBDK() ) {
assert(!isCS() && !isExpanded());
delta =
((**i).id() == ParticleID::g ? xgammaGluon : xgammaQuark) *
log(mu2/(2.*(meMomenta()[idi]*meMomenta()[idj])));
if ( (idi < 2 && idj < 2) ||
(idi > 1 && idj > 1) )
delta +=
((**i).id() == ParticleID::g ? CA : CF) * sqr(pi) / 2.;
}
if ( isExpanded() ) {
assert(!isCS() && !isBDK());
double theLog = log(mu2/(2.*(meMomenta()[idi]*meMomenta()[idj])));
delta =
((**i).id() == ParticleID::g ? CA : CF) * 0.5 * sqr(theLog) +
((**i).id() == ParticleID::g ? xgammaGluon : xgammaQuark) * theLog;
}
res +=
( ((**i).id() == ParticleID::g ? CA : CF) * (-sqr(pi)/3.) +
((**i).id() == ParticleID::g ? gammaGluon : gammaQuark) +
((**i).id() == ParticleID::g ? KGluon : KQuark) +
delta ) *
lastBorn()->colourCorrelatedME2(make_pair(idi,idj));
}
}
Energy2 muR2 =
lastBorn()->renormalizationScale()*
sqr(lastBorn()->renormalizationScaleFactor());
if ( muR2 != mu2 ) {
res -=
betaZero *
lastBorn()->orderInAlphaS() * log(muR2/mu2) *
lastBorn()->me2();
}
// include the finite renormalization for DR here; ATTENTION this
// has to be mentioned in the manual! see hep-ph/9305239 for
// details; this guarantees an expansion in alpha_s^\bar{MS} when
// using dimensional reduction
- if ( isDR() )
+ if ( isDR() && isDRbar() )
res -= (CA/6.)*lastBorn()->orderInAlphaS()*lastBorn()->me2();
res *= ( - lastBorn()->lastAlphaS() / (2.*pi) );
return res;
}
double DipoleIOperator::oneLoopDoublePole() const {
if ( !isExpanded() )
return 0.;
double res = 0.;
for ( cPDVector::const_iterator i = mePartonData().begin();
i != mePartonData().end(); ++i ) {
if ( !apply(*i) )
continue;
res += (**i).id() == ParticleID::g ? CA : CF;
}
res *= ( - lastBorn()->lastAlphaS() / (2.*pi) ) * ( - lastBorn()->me2() );
return res;
}
double DipoleIOperator::oneLoopSinglePole() const {
if ( !isExpanded() )
return 0.;
double res = 0.;
int idi = 0; int idj = 0;
Energy2 mu2 = lastBorn()->mu2();
for ( cPDVector::const_iterator i = mePartonData().begin();
i != mePartonData().end(); ++i, ++idi ) {
if ( !apply(*i) )
continue;
idj = 0;
for ( cPDVector::const_iterator j = mePartonData().begin();
j != mePartonData().end(); ++j, ++idj ) {
if ( !apply(*j) )
continue;
if ( i == j || lastBorn()->noDipole(idi,idj) )
continue;
double xgammaGluon = gammaGluon;
double xgammaQuark = gammaQuark;
if ( isDR() ) {
xgammaGluon += CA/6.;
xgammaQuark += CF/2.;
}
double theLog = log(mu2/(2.*(meMomenta()[idi]*meMomenta()[idj])));
double delta =
((**i).id() == ParticleID::g ? CA : CF) * theLog +
((**i).id() == ParticleID::g ? xgammaGluon : xgammaQuark);
res +=
delta * lastBorn()->colourCorrelatedME2(make_pair(idi,idj));
}
}
res *= ( - lastBorn()->lastAlphaS() / (2.*pi) );
return res;
}
void DipoleIOperator::persistentOutput(PersistentOStream & os) const {
os << CA << CF << gammaQuark << gammaGluon << betaZero
<< KQuark << KGluon;
}
void DipoleIOperator::persistentInput(PersistentIStream & is, int) {
is >> CA >> CF >> gammaQuark >> gammaGluon >> betaZero
>> KQuark >> KGluon;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<DipoleIOperator,MatchboxInsertionOperator>
describeHerwigDipoleIOperator("Herwig::DipoleIOperator", "HwMatchbox.so");
void DipoleIOperator::Init() {
static ClassDocumentation<DipoleIOperator> documentation
("DipoleIOperator");
DipoleRepository::registerInsertionOperator<0,DipoleIOperator>("LightIOperator");
}

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 10:04 AM (1 d, 11 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111038
Default Alt Text
(22 KB)

Event Timeline