Page MenuHomeHEPForge
Files F7877931

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
@@ -1,1201 +1,1281 @@
// -*- C++ -*-
//
// MatchboxMEBase.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 MatchboxMEBase class.
//
#include "MatchboxMEBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PDF.h"
#include "ThePEG/PDT/PDT.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
MatchboxMEBase::MatchboxMEBase()
: MEBase(),
theOneLoop(false),
theOneLoopNoBorn(false) {}
MatchboxMEBase::~MatchboxMEBase() {}
Ptr<MatchboxFactory>::tcptr MatchboxMEBase::factory() const { return theFactory; }
void MatchboxMEBase::factory(Ptr<MatchboxFactory>::tcptr f) { theFactory = f; }
Ptr<Tree2toNGenerator>::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); }
Ptr<ProcessData>::tptr MatchboxMEBase::processData() const { return factory()->processData(); }
unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); }
double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); }
double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); }
bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); }
bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); }
bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); }
bool MatchboxMEBase::verbose() const { return factory()->verbose(); }
bool MatchboxMEBase::initVerbose() const { return factory()->initVerbose(); }
void MatchboxMEBase::getDiagrams() const {
if ( diagramGenerator() && processData() ) {
vector<Ptr<Tree2toNDiagram>::ptr> diags;
for ( vector<PDVector>::const_iterator p = subProcesses().begin();
p != subProcesses().end(); ++p ) {
vector<Ptr<Tree2toNDiagram>::ptr>& res =
processData()->diagramMap()[*p];
if ( res.empty() ) {
res = diagramGenerator()->generate(*p,orderInAlphaS(),orderInAlphaEW());
}
copy(res.begin(),res.end(),back_inserter(diags));
processData()->fillMassGenerators(*p);
}
if ( diags.empty() )
return;
for ( vector<Ptr<Tree2toNDiagram>::ptr>::iterator d = diags.begin();
d != diags.end(); ++d ) {
add(*d);
}
return;
}
throw Exception()
<< "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
<< "Please check your setup." << Exception::abortnow;
}
Selector<MEBase::DiagramIndex>
MatchboxMEBase::diagrams(const DiagramVector & diags) const {
if ( phasespace() ) {
return phasespace()->selectDiagrams(diags);
}
throw Exception()
<< "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return Selector<MEBase::DiagramIndex>();
}
Selector<const ColourLines *>
MatchboxMEBase::colourGeometries(tcDiagPtr diag) const {
if ( matchboxAmplitude() ) {
if ( !matchboxAmplitude()->haveColourFlows() )
throw Exception() << "A colour flow implementation is not present."
<< Exception::abortnow;
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
return matchboxAmplitude()->colourGeometries(diag);
}
throw Exception()
<< "MatchboxMEBase::colourGeometries() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return Selector<const ColourLines *>();
}
unsigned int MatchboxMEBase::orderInAlphaS() const {
if ( matchboxAmplitude() ) {
return matchboxAmplitude()->orderInGs();
}
throw Exception()
<< "MatchboxMEBase::orderInAlphaS() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0;
}
unsigned int MatchboxMEBase::orderInAlphaEW() const {
if ( matchboxAmplitude() ) {
return matchboxAmplitude()->orderInGem();
}
throw Exception()
<< "MatchboxMEBase::orderInAlphaEW() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0;
}
void MatchboxMEBase::setXComb(tStdXCombPtr xc) {
MEBase::setXComb(xc);
lastMatchboxXComb(xc);
if ( phasespace() )
phasespace()->setXComb(xc);
if ( scaleChoice() )
scaleChoice()->setXComb(xc);
if ( matchboxAmplitude() )
matchboxAmplitude()->setXComb(xc);
}
double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) {
// shamelessly stolen from PartonExtractor.cc
Energy2 shmax = lastCuts().sHatMax();
Energy2 shmin = lastCuts().sHatMin();
Energy2 sh = shmin*pow(shmax/shmin, *r1);
double ymax = lastCuts().yHatMax();
double ymin = lastCuts().yHatMin();
double km = log(shmax/shmin);
ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));
double y = ymin + (*r2)*(ymax - ymin);
double x1 = exp(-0.5*log(lastS()/sh) + y);
double x2 = exp(-0.5*log(lastS()/sh) - y);
Lorentz5Momentum P1 = lastParticles().first->momentum();
LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy());
p1.rotateY(P1.theta());
p1.rotateZ(P1.phi());
meMomenta()[0] = p1;
Lorentz5Momentum P2 = lastParticles().second->momentum();
LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy());
p2.rotateY(P2.theta());
p2.rotateZ(P2.phi());
meMomenta()[1] = p2;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
return km*(ymax - ymin);
}
bool MatchboxMEBase::generateKinematics(const double * r) {
if ( phasespace() ) {
jacobian(phasespace()->generateKinematics(r,meMomenta()));
if ( jacobian() == 0.0 )
return false;
setScale();
logGenerateKinematics(r);
assert(lastMatchboxXComb());
if ( nDimAmplitude() > 0 ) {
amplitudeRandomNumbers().resize(nDimAmplitude());
copy(r + nDimPhasespace(),
r + nDimPhasespace() + nDimAmplitude(),
amplitudeRandomNumbers().begin());
}
if ( nDimInsertions() > 0 ) {
insertionRandomNumbers().resize(nDimInsertions());
copy(r + nDimPhasespace() + nDimAmplitude(),
r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(),
insertionRandomNumbers().begin());
}
return true;
}
throw Exception()
<< "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return false;
}
int MatchboxMEBase::nDim() const {
if ( lastMatchboxXComb() )
return nDimPhasespace() + nDimAmplitude() + nDimInsertions();
int ampAdd = 0;
if ( matchboxAmplitude() ) {
ampAdd = matchboxAmplitude()->nDimAdditional();
}
int insertionAdd = 0;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
insertionAdd = max(insertionAdd,(**v).nDimAdditional());
}
return nDimBorn() + ampAdd + insertionAdd;
}
int MatchboxMEBase::nDimBorn() const {
if ( lastMatchboxXComb() )
return nDimPhasespace();
if ( phasespace() ) {
size_t nout = diagrams().front()->partons().size()-2;
int n = phasespace()->nDim(nout);
if ( phasespace()->useMassGenerators() ) {
for ( cPDVector::const_iterator pd =
diagrams().front()->partons().begin();
pd != diagrams().front()->partons().end(); ++pd ) {
if ( processData()->massGenerator(*pd) ||
(**pd).width() != ZERO ) {
++n;
}
}
}
return n;
}
throw Exception()
<< "MatchboxMEBase::nDim() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0;
}
void MatchboxMEBase::setScale() const {
if ( haveX1X2() ) {
lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
}
Energy2 fscale = factorizationScale()*sqr(factorizationScaleFactor());
Energy2 rscale = renormalizationScale()*sqr(renormalizationScaleFactor());
Energy2 ewrscale = renormalizationScaleQED();
lastXCombPtr()->lastScale(fscale);
if ( !fixedCouplings() ) {
if ( rscale > lastCuts().scaleMin() )
lastXCombPtr()->lastAlphaS(SM().alphaS(rscale));
else
lastXCombPtr()->lastAlphaS(SM().alphaS(lastCuts().scaleMin()));
} else {
lastXCombPtr()->lastAlphaS(SM().alphaS());
}
if ( !fixedQEDCouplings() ) {
lastXCombPtr()->lastAlphaEM(SM().alphaEM(ewrscale));
} else {
lastXCombPtr()->lastAlphaEM(SM().alphaEMMZ());
}
logSetScale();
}
Energy2 MatchboxMEBase::factorizationScale() const {
if ( scaleChoice() ) {
return scaleChoice()->factorizationScale();
}
throw Exception()
<< "MatchboxMEBase::factorizationScale() expects a MatchboxScaleChoice object.\n"
<< "Please check your setup." << Exception::abortnow;
return ZERO;
}
Energy2 MatchboxMEBase::renormalizationScale() const {
if ( scaleChoice() ) {
return scaleChoice()->renormalizationScale();
}
throw Exception()
<< "MatchboxMEBase::renormalizationScale() expects a MatchboxScaleChoice object.\n"
<< "Please check your setup." << Exception::abortnow;
return ZERO;
}
Energy2 MatchboxMEBase::renormalizationScaleQED() const {
if ( scaleChoice() ) {
return scaleChoice()->renormalizationScaleQED();
}
return renormalizationScale();
}
void MatchboxMEBase::setVetoScales(tSubProPtr) const {}
void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const {
if ( !mePartonData()[0]->coloured() &&
!mePartonData()[1]->coloured() ) {
lastMEPDFWeight(1.0);
logPDFWeight();
return;
}
double w = 1.;
if ( mePartonData()[0]->coloured() && havePDFWeight1() )
w *= pdf1(factorizationScale);
if ( mePartonData()[1]->coloured() && havePDFWeight2() )
w *= pdf2(factorizationScale);
lastMEPDFWeight(w);
logPDFWeight();
}
double MatchboxMEBase::pdf1(Energy2 fscale, double xEx) const {
assert(lastXCombPtr()->partonBins().first->pdf());
if ( xEx < 1. && lastX1() >= xEx ) {
return
( ( 1. - lastX1() ) / ( 1. - xEx ) ) *
lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
lastPartons().first->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
xEx)/xEx;
}
return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
lastPartons().first->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX1())/lastX1();
}
double MatchboxMEBase::pdf2(Energy2 fscale, double xEx) const {
assert(lastXCombPtr()->partonBins().second->pdf());
if ( xEx < 1. && lastX2() >= xEx ) {
return
( ( 1. - lastX2() ) / ( 1. - xEx ) ) *
lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
lastPartons().second->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
xEx)/xEx;
}
return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
lastPartons().second->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX2())/lastX2();
}
double MatchboxMEBase::me2() const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->me2()*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::finalStateSymmetry() const {
if ( symmetryFactor() > 0.0 )
return symmetryFactor();
double sFactor = 1.;
map<long,int> counts;
cPDVector checkData;
copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData));
cPDVector::iterator p = checkData.begin();
while ( !checkData.empty() ) {
if ( counts.find((**p).id()) != counts.end() ) {
counts[(**p).id()] += 1;
} else {
counts[(**p).id()] = 1;
}
checkData.erase(p);
p = checkData.begin();
continue;
}
for ( map<long,int>::const_iterator c = counts.begin();
c != counts.end(); ++c ) {
if ( c->second == 1 )
continue;
if ( c->second == 2 )
sFactor /= 2.;
else if ( c->second == 3 )
sFactor /= 6.;
else if ( c->second == 4 )
sFactor /= 24.;
}
symmetryFactor(sFactor);
return symmetryFactor();
}
double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const {
// assume that we always have incoming
// spin-1/2 or massless spin-1 particles
double fac = 1./4.;
if ( orderInAlphaS() > 0 || addAlphaS != 0 )
fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS));
if ( orderInAlphaEW() > 0 )
fac *= pow(lastAlphaEM()/SM().alphaEM(),double(orderInAlphaEW()));
if ( mePartonData()[0]->iColour() == PDT::Colour3 ||
mePartonData()[0]->iColour() == PDT::Colour3bar )
fac /= SM().Nc();
else if ( mePartonData()[0]->iColour() == PDT::Colour8 )
fac /= (SM().Nc()*SM().Nc()-1.);
if ( mePartonData()[1]->iColour() == PDT::Colour3 ||
mePartonData()[1]->iColour() == PDT::Colour3bar )
fac /= SM().Nc();
else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
fac /= (SM().Nc()*SM().Nc()-1.);
return finalStateSymmetry()*fac;
}
CrossSection MatchboxMEBase::dSigHatDR() const {
getPDFWeight();
if ( !lastXCombPtr()->willPassCuts() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double xme2 = me2();
lastME2(xme2);
if ( xme2 == 0. && !oneLoopNoBorn() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double vme2 = 0.;
if ( oneLoop() )
vme2 = oneLoopInterference();
CrossSection res = ZERO;
if ( !oneLoopNoBorn() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * xme2;
if ( oneLoop() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * vme2;
if ( !onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).setXComb(lastXCombPtr());
res += (**v).dSigHatDR();
}
+ if ( checkPoles() )
+ logPoles();
}
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
lastMECrossSection(res);
return lastMECrossSection();
}
double MatchboxMEBase::oneLoopInterference() const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->oneLoopAmplitudes() )
matchboxAmplitude()->prepareOneLoopAmplitudes(this);
lastME2(matchboxAmplitude()->oneLoopInterference()*
crossingSign()*
me2Norm(1));
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
+MatchboxMEBase::AccuracyHistogram::AccuracyHistogram(double low,
+ double up,
+ unsigned int nbins)
+ : lower(low), sameSign(0), oppositeSign(0), nans(0) {
+
+ double step = (up-low)/nbins;
+
+ for ( unsigned int k = 1; k <= nbins; ++k )
+ bins[lower + k*step] = 0.0;
+
+}
+
+void MatchboxMEBase::AccuracyHistogram::book(double a, double b) {
+ if ( isnan(a) || isnan(b) ||
+ isinf(a) || isinf(b) ) {
+ ++nans;
+ return;
+ }
+ if ( a*b >= 0. )
+ ++sameSign;
+ if ( a*b < 0. )
+ ++oppositeSign;
+ double r = 1.;
+ if ( abs(a) != 0.0 )
+ r = 1.-abs(b/a);
+ else if ( abs(b) != 0.0 )
+ r = abs(b);
+ map<double,double>::iterator bin =
+ bins.upper_bound(log(r));
+ if ( bin == bins.end() )
+ return;
+ bin->second += 1.;
+}
+
+void MatchboxMEBase::AccuracyHistogram::dump(const std::string& prefix,
+ const cPDVector& proc) const {
+ ostringstream fname("");
+ for ( cPDVector::const_iterator p = proc.begin();
+ p != proc.end(); ++p )
+ fname << (**p).PDGName();
+ ofstream out((prefix+fname.str()+".dat").c_str());
+ out << "# same sign : " << sameSign << " opposite sign : "
+ << oppositeSign << " nans : " << nans << "\n";
+ for ( map<double,double>::const_iterator b = bins.begin();
+ b != bins.end(); ++b ) {
+ map<double,double>::const_iterator bp = b; --bp;
+ if ( b->second != 0. ) {
+ if ( b != bins.begin() )
+ out << bp->first;
+ else
+ out << lower;
+ out << " " << b->first
+ << " " << b->second
+ << "\n" << flush;
+ }
+ }
+}
+
+void MatchboxMEBase::AccuracyHistogram::persistentOutput(PersistentOStream& os) const {
+ os << lower << bins;
+}
+
+void MatchboxMEBase::AccuracyHistogram::persistentInput(PersistentIStream& is) {
+ is >> lower >> bins;
+}
+
void MatchboxMEBase::logPoles() const {
double res2me = oneLoopDoublePole();
double res1me = oneLoopSinglePole();
double res2i = 0.;
double res1i = 0.;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
res2i += (**v).oneLoopDoublePole();
res1i += (**v).oneLoopSinglePole();
}
- double diff2 = abs(res2me) != 0. ? 1.-abs(res2i/res2me) : abs(res2i)-abs(res2me);
- double diff1 = abs(res1me) != 0. ? 1.-abs(res1i/res1me) : abs(res1i)-abs(res1me);
- generator()->log()
- << "check "
- << log10(abs(diff2)) << " " << log10(abs(diff1)) << "\n"
- << flush;
+ epsilonSquarePoleHistograms[mePartonData()].book(res2me,res2i);
+ epsilonPoleHistograms[mePartonData()].book(res1me,res1i);
}
bool MatchboxMEBase::haveOneLoop() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->haveOneLoop();
return false;
}
bool MatchboxMEBase::onlyOneLoop() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->onlyOneLoop();
return false;
}
bool MatchboxMEBase::isDR() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isDR();
return false;
}
bool MatchboxMEBase::isCS() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isCS();
return false;
}
bool MatchboxMEBase::isBDK() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isBDK();
return false;
}
bool MatchboxMEBase::isExpanded() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isExpanded();
return false;
}
Energy2 MatchboxMEBase::mu2() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->mu2();
return 0*GeV2;
}
double MatchboxMEBase::oneLoopDoublePole() const {
if ( matchboxAmplitude() ) {
return
matchboxAmplitude()->oneLoopDoublePole()*
crossingSign()*
me2Norm(1);
}
return 0.;
}
double MatchboxMEBase::oneLoopSinglePole() const {
if ( matchboxAmplitude() ) {
return
matchboxAmplitude()->oneLoopSinglePole()*
crossingSign()*
me2Norm(1);
}
return 0.;
}
vector<Ptr<SubtractionDipole>::ptr>
MatchboxMEBase::getDipoles(const vector<Ptr<SubtractionDipole>::ptr>& dipoles,
const vector<Ptr<MatchboxMEBase>::ptr>& borns) const {
vector<Ptr<SubtractionDipole>::ptr> res;
// keep track of the dipoles we already did set up
set<pair<pair<pair<int,int>,int>,pair<Ptr<MatchboxMEBase>::tptr,Ptr<SubtractionDipole>::tptr> > > done;
cPDVector rep = diagrams().front()->partons();
int nreal = rep.size();
// now loop over configs
for ( int emitter = 0; emitter < nreal; ++emitter ) {
for ( int spectator = 0; spectator < nreal; ++spectator ) {
if ( emitter == spectator )
continue;
for ( int emission = 2; emission < nreal; ++emission ) {
if ( emission == emitter || emission == spectator )
continue;
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator b =
borns.begin(); b != borns.end(); ++b ) {
if ( (**b).onlyOneLoop() )
continue;
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
dipoles.begin(); d != dipoles.end(); ++d ) {
if ( !rep[emitter]->coloured() ||
!rep[emission]->coloured() ||
!rep[spectator]->coloured() ) {
continue;
}
if ( noDipole(emitter,emission,spectator) ) {
continue;
}
if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(*b,*d)))
!= done.end() ) {
continue;
}
if ( !(**d).canHandle(rep,emitter,emission,spectator) ) {
continue;
}
// now get to work
Ptr<SubtractionDipole>::ptr nDipole = (**d).cloneMe();
nDipole->realEmitter(emitter);
nDipole->realEmission(emission);
nDipole->realSpectator(spectator);
nDipole->realEmissionME(const_cast<MatchboxMEBase*>(this));
nDipole->underlyingBornME(*b);
nDipole->setupBookkeeping();
if ( !(nDipole->empty()) ) {
res.push_back(nDipole);
done.insert(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(*b,*d)));
if ( nDipole->isSymmetric() )
done.insert(make_pair(make_pair(make_pair(emission,emitter),spectator),make_pair(*b,*d)));
ostringstream dname;
dname << fullName() << "." << (**b).name() << "."
<< (**d).name() << ".[("
<< emitter << "," << emission << ")," << spectator << "]";
if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) )
throw InitException() << "Dipole " << dname.str() << " already existing.";
nDipole->cloneDependencies(dname.str());
}
}
}
}
}
}
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = res.begin();
d != res.end(); ++d )
(**d).partnerDipoles(res);
return res;
}
double MatchboxMEBase::colourCorrelatedME2(pair<int,int> ij) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->colourCorrelatedME2(ij)*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr largeNBasis) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
largeNBasis->prepare(mePartonData(),false);
lastME2(matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
void MatchboxMEBase::flushCaches() {
MEBase::flushCaches();
if ( matchboxAmplitude() )
matchboxAmplitude()->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).flushCaches();
}
}
void MatchboxMEBase::print(ostream& os) const {
os << "--- MatchboxMEBase setup -------------------------------------------------------\n";
os << " '" << name() << "' for subprocesses:\n";
for ( vector<PDVector>::const_iterator p = subProcesses().begin();
p != subProcesses().end(); ++p ) {
os << " ";
for ( PDVector::const_iterator pp = p->begin();
pp != p->end(); ++pp ) {
os << (**pp).PDGName() << " ";
if ( pp == p->begin() + 1 )
os << "-> ";
}
os << "\n";
}
os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections";
if ( oneLoopNoBorn() )
os << " without Born contributions";
os << "\n";
if ( oneLoop() && !onlyOneLoop() ) {
os << " using insertion operators\n";
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
os << " '" << (**v).name() << "' with "
<< ((**v).isDR() ? "" : "C") << "DR/";
if ( (**v).isCS() )
os << "CS";
if ( (**v).isBDK() )
os << "BDK";
if ( (**v).isExpanded() )
os << "expanded";
os << " conventions\n";
}
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxMEBase::printLastEvent(ostream& os) const {
os << "--- MatchboxMEBase last event information --------------------------------------\n";
os << " for matrix element '" << name() << "'\n";
os << " process considered:\n ";
int in = 0;
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
if ( ++in == 2 )
os << " -> ";
}
os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n"
<< " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2)
<< " x1 = " << lastX1() << " x2 = " << lastX2()
<< " alphaS = " << lastAlphaS() << "\n";
os << " momenta/GeV generated from random numbers\n ";
copy(lastXComb().lastRandomNumbers().begin(),
lastXComb().lastRandomNumbers().end(),ostream_iterator<double>(os," "));
os << ":\n ";
for ( vector<Lorentz5Momentum>::const_iterator p = meMomenta().begin();
p != meMomenta().end(); ++p ) {
os << (*p/GeV) << "\n ";
}
os << "last cross section/nb calculated was:\n "
<< (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n";
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxMEBase::logGenerateKinematics(const double * r) const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' generated kinematics\nfrom "
<< nDim() << " random numbers:\n";
copy(r,r+nDim(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "storing phase space information in XComb "
<< lastXCombPtr() << "\n";
generator()->log() << "generated phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
cPDVector::const_iterator dit = mePartonData().begin();
for ( ; pit != meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
<< "and Jacobian = " << jacobian() << " sHat/GeV2 = "
<< (lastSHat()/GeV2) << "\n" << flush;
}
void MatchboxMEBase::logSetScale() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n"
<< "scale/GeV2 = " << (scale()/GeV2) << " xi_R = "
<< renormalizationScaleFactor() << " xi_F = "
<< factorizationScaleFactor() << "\n"
<< "alpha_s = " << lastAlphaS() << "\n" << flush;
}
void MatchboxMEBase::logPDFWeight() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' calculated pdf weight = "
<< lastMEPDFWeight() << " from XComb "
<< lastXCombPtr() << "\n"
<< "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) "
<< "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n"
<< flush;
}
void MatchboxMEBase::logME2() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' evaluated me2 using XComb "
<< lastXCombPtr() << "\n"
<< "and phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
cPDVector::const_iterator dit = mePartonData().begin();
for ( ; pit != meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
<< "sHat/GeV2 = " << (lastSHat()/GeV2)
<< " me2 = " << lastME2() << "\n" << flush;
}
void MatchboxMEBase::logDSigHatDR() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' evaluated cross section using XComb "
<< lastXCombPtr() << "\n"
<< "Jacobian = " << jacobian() << " sHat/GeV2 = "
<< (lastSHat()/GeV2) << " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void MatchboxMEBase::cloneDependencies(const std::string& prefix) {
if ( phasespace() ) {
Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
throw InitException() << "Phasespace generator " << pname.str() << " already existing.";
myPhasespace->cloneDependencies(pname.str());
phasespace(myPhasespace);
}
theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
if ( matchboxAmplitude() ) {
Ptr<MatchboxAmplitude>::ptr myAmplitude = matchboxAmplitude()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myAmplitude->name();
if ( ! (generator()->preinitRegister(myAmplitude,pname.str()) ) )
throw InitException() << "Amplitude " << pname.str() << " already existing.";
myAmplitude->cloneDependencies(pname.str());
matchboxAmplitude(myAmplitude);
amplitude(myAmplitude);
matchboxAmplitude()->orderInGs(orderInAlphaS());
matchboxAmplitude()->orderInGem(orderInAlphaEW());
}
if ( scaleChoice() ) {
Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
throw InitException() << "Scale choice " << pname.str() << " already existing.";
scaleChoice(myScaleChoice);
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw InitException() << "Reweight " << pname.str() << " already existing.";
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
Ptr<MatchboxInsertionOperator>::ptr myIOP = (**v).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**v).name();
if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
throw InitException() << "Insertion operator " << pname.str() << " already existing.";
*v = myIOP;
}
}
void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const {
if ( phasespace() ) {
size_t nout = diagrams().front()->partons().size()-2;
xc.nDimPhasespace(phasespace()->nDim(nout));
}
if ( matchboxAmplitude() ) {
xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional());
if ( matchboxAmplitude()->colourBasis() ) {
size_t cdim =
matchboxAmplitude()->colourBasis()->prepare(diagrams(),matchboxAmplitude()->noCorrelations());
xc.colourBasisDim(cdim);
}
}
int insertionAdd = 0;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
insertionAdd = max(insertionAdd,(**v).nDimAdditional());
}
xc.nDimInsertions(insertionAdd);
xc.nLight(getNLight());
if ( initVerbose() ) {
string fname = name() + ".diagrams";
ifstream test(fname.c_str());
if ( !test ) {
test.close();
ofstream out(fname.c_str());
for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
d != diagrams().end(); ++d ) {
DiagramDrawer::drawDiag(out,dynamic_cast<const Tree2toNDiagram&>(**d));
out << "\n";
}
}
}
}
StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec&,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
Ptr<MatchboxXComb>::ptr xc =
new_ptr(MatchboxXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts, newME,
newDiagrams, mir,
newHead));
prepareXComb(*xc);
return xc;
}
StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
Ptr<MatchboxXComb>::ptr xc =
new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams));
prepareXComb(*xc);
return xc;
}
void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theFactory << thePhasespace
<< theAmplitude << theScaleChoice << theVirtuals
<< theReweights << theSubprocesses << theOneLoop
- << theOneLoopNoBorn;
+ << theOneLoopNoBorn
+ << epsilonSquarePoleHistograms << epsilonPoleHistograms;
}
void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theFactory >> thePhasespace
>> theAmplitude >> theScaleChoice >> theVirtuals
>> theReweights >> theSubprocesses >> theOneLoop
- >> theOneLoopNoBorn;
+ >> theOneLoopNoBorn
+ >> epsilonSquarePoleHistograms >> epsilonPoleHistograms;
lastMatchboxXComb(theLastXComb);
}
void MatchboxMEBase::Init() {
static ClassDocumentation<MatchboxMEBase> documentation
("MatchboxMEBase is the base class for matrix elements "
"in the context of the matchbox NLO interface.");
}
IBPtr MatchboxMEBase::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxMEBase::fullclone() const {
return new_ptr(*this);
}
void MatchboxMEBase::doinit() {
MEBase::doinit();
if ( !theAmplitude )
theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
}
+void MatchboxMEBase::dofinish() {
+ MEBase::dofinish();
+ for ( map<cPDVector,AccuracyHistogram>::const_iterator
+ b = epsilonSquarePoleHistograms.begin();
+ b != epsilonSquarePoleHistograms.end(); ++b ) {
+ b->second.dump(factory()->poleData() + "epsilonSquarePoles-",b->first);
+ }
+ for ( map<cPDVector,AccuracyHistogram>::const_iterator
+ b = epsilonPoleHistograms.begin();
+ b != epsilonPoleHistograms.end(); ++b ) {
+ b->second.dump(factory()->poleData() + "epsilonPoles-",b->first);
+ }
+}
+
// *** 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<MatchboxMEBase,MEBase>
describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.h b/MatrixElement/Matchbox/Base/MatchboxMEBase.h
--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.h
+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.h
@@ -1,822 +1,910 @@
// -*- C++ -*-
//
// MatchboxMEBase.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_MatchboxMEBase_H
#define HERWIG_MatchboxMEBase_H
//
// This is the declaration of the MatchboxMEBase class.
//
#include "ThePEG/MatrixElement/MEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/ProcessData.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxReweightBase.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig++/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig++/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxMEBase is the base class for matrix elements
* in the context of the matchbox NLO interface.
*
* @see \ref MatchboxMEBaseInterfaces "The interfaces"
* defined for MatchboxMEBase.
*/
class MatchboxMEBase:
public MEBase, public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxMEBase();
/**
* The destructor.
*/
virtual ~MatchboxMEBase();
//@}
public:
/**
* Return the factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tcptr factory() const;
/**
* Set the factory which produced this matrix element
*/
void factory(Ptr<MatchboxFactory>::tcptr f);
/** @name Subprocess and diagram information. */
//@{
/**
* Return the subprocesses.
*/
const vector<PDVector>& subProcesses() const { return theSubprocesses; }
/**
* Access the subprocesses.
*/
vector<PDVector>& subProcesses() { return theSubprocesses; }
/**
* Return the diagram generator.
*/
Ptr<Tree2toNGenerator>::tptr diagramGenerator() const;
/**
* Return the process data.
*/
Ptr<ProcessData>::tptr processData() const;
/**
* Return true, if this matrix element does not want to
* make use of mirroring processes; in this case all
* possible partonic subprocesses with a fixed assignment
* of incoming particles need to be provided through the diagrams
* added with the add(...) method.
*/
virtual bool noMirror () const { return true; }
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
using MEBase::getDiagrams;
/**
* With the information previously supplied with the
* setKinematics(...) method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector &) const;
using MEBase::diagrams;
/**
* 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 the order in \f$\alpha_S\f$ in which this matrix element
* is given.
*/
virtual unsigned int orderInAlphaS() const;
using MEBase::orderInAlphaS;
/**
* Return the order in \f$\alpha_{EM}\f$ in which this matrix
* element is given. Returns 0.
*/
virtual unsigned int orderInAlphaEW() const;
using MEBase::orderInAlphaEW;
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
virtual unsigned int getNLight() const;
//@}
/** @name Phasespace generation */
//@{
/**
* Return the phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::tptr phasespace() const { return thePhasespace; }
/**
* Set the phase space generator to be used.
*/
void phasespace(Ptr<MatchboxPhasespace>::ptr ps) { thePhasespace = ps; }
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr xc);
/**
* Return true, if the XComb steering this matrix element
* should keep track of the random numbers used to generate
* the last phase space point
*/
virtual bool keepRandomNumbers() const { return true; }
/**
* Generate incoming parton momenta. This default
* implementation performs the standard mapping
* from x1,x2 -> tau,y making 1/tau flat; incoming
* parton momenta are stored in meMomenta()[0,1],
* only massless partons are supported so far;
* return the Jacobian of the mapping
*/
double generateIncomingPartons(const double* r1, const double* r2);
/**
* Generate internal degrees of freedom given nDim() uniform random
* numbers in the interval ]0,1[. To help the phase space generator,
* the 'dSigHatDR' should be a smooth function of these numbers,
* although this is not strictly necessary. The return value should
* be true of the generation succeeded. If so the generated momenta
* should be stored in the meMomenta() vector. Derived classes
* must call this method once internal degrees of freedom are setup
* and finally return the result of this method.
*/
virtual bool generateKinematics(const double * r);
/**
* The number of internal degreed of freedom used in the matrix
* element.
*/
virtual int nDim() const;
/**
* The number of internal degrees of freedom used in the matrix
* element for generating a Born phase space point
*/
virtual int nDimBorn() const;
/**
* Return true, if this matrix element will generate momenta for the
* incoming partons itself. The matrix element is required to store
* the incoming parton momenta in meMomenta()[0,1]. No mapping in
* tau and y is performed by the PartonExtractor object, if a
* derived class returns true here. The phase space jacobian is to
* include a factor 1/(x1 x2).
*/
virtual bool haveX1X2() const {
return
(phasespace() ? phasespace()->haveX1X2() : false) ||
diagrams().front()->partons().size() == 3;
}
/**
* Return true, if this matrix element expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS() const {
return
(phasespace() ? phasespace()->wantCMS() : true) &&
diagrams().front()->partons().size() != 3; }
/**
* Return the meMomenta as generated at the last
* phase space point.
*/
const vector<Lorentz5Momentum>& lastMEMomenta() const { return meMomenta(); }
/**
* Access the meMomenta.
*/
vector<Lorentz5Momentum>& lastMEMomenta() { return meMomenta(); }
//@}
/** @name Scale choices, couplings and PDFs */
//@{
/**
* Set the scale choice object
*/
void scaleChoice(Ptr<MatchboxScaleChoice>::ptr sc) { theScaleChoice = sc; }
/**
* Return the scale choice object
*/
Ptr<MatchboxScaleChoice>::tptr scaleChoice() const { return theScaleChoice; }
/**
* Set scales and alphaS
*/
void setScale() const;
/**
* Return the scale associated with the phase space point provided
* by the last call to setKinematics().
*/
virtual Energy2 scale() const { return lastScale(); }
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 factorizationScale() const;
/**
* Get the factorization scale factor
*/
virtual double factorizationScaleFactor() const;
/**
* Return the (QCD) renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScale() const;
/**
* Get the renormalization scale factor
*/
virtual double renormalizationScaleFactor() const;
/**
* Return the QED renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScaleQED() const;
/**
* Set veto scales on the particles at the given
* SubProcess which has been generated using this
* matrix element.
*/
virtual void setVetoScales(tSubProPtr) const;
/**
* Return true, if fixed couplings are used.
*/
bool fixedCouplings() const;
/**
* Return true, if fixed couplings are used.
*/
bool fixedQEDCouplings() const;
/**
* Return the value of \f$\alpha_S\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaS(scale()).
*/
virtual double alphaS() const { return lastAlphaS(); }
/**
* Return the value of \f$\alpha_EM\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaEM(scale()).
*/
virtual double alphaEM() const { return lastAlphaEM(); }
/**
* Return true, if this matrix element provides the PDF
* weight for the first incoming parton itself.
*/
virtual bool havePDFWeight1() const {
return diagrams().front()->partons()[0]->coloured();
}
/**
* Return true, if this matrix element provides the PDF
* weight for the second incoming parton itself.
*/
virtual bool havePDFWeight2() const {
return diagrams().front()->partons()[1]->coloured();
}
/**
* Set the PDF weight.
*/
void getPDFWeight(Energy2 factorizationScale = ZERO) const;
/**
* Supply the PDF weight for the first incoming parton.
*/
double pdf1(Energy2 factorizationScale = ZERO,
double xEx = 1.) const;
/**
* Supply the PDF weight for the second incoming parton.
*/
double pdf2(Energy2 factorizationScale = ZERO,
double xEx = 1.) const;
//@}
/** @name Amplitude information and matrix element evaluation */
//@{
/**
* Return the amplitude.
*/
Ptr<MatchboxAmplitude>::tptr matchboxAmplitude() const { return theAmplitude; }
/**
* Set the amplitude.
*/
void matchboxAmplitude(Ptr<MatchboxAmplitude>::ptr amp) { theAmplitude = amp; }
/**
* Return the matrix element for the kinematical configuation
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
*/
virtual double me2() const;
/**
* Return the symmetry factor for identical final state particles.
*/
virtual double finalStateSymmetry() const;
/**
* Return the normalizing factor for the matrix element averaged
* over quantum numbers and including running couplings.
*/
double me2Norm(unsigned int addAlphaS = 0) const;
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR() const;
//@}
/** @name One-loop corrections */
//@{
/**
* Return the one-loop/tree interference.
*/
virtual double oneLoopInterference() const;
/**
* Return true, if this matrix element is capable of calculating
* one-loop (QCD) corrections.
*/
virtual bool haveOneLoop() const;
/**
* Return true, if this matrix element only provides
* one-loop (QCD) corrections.
*/
virtual bool onlyOneLoop() const;
/**
* 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 true, if one loop corrections are given in the conventions
* of the integrated dipoles.
*/
virtual bool isCS() const;
/**
* Return true, if one loop corrections are given in the conventions
* of BDK.
*/
virtual bool isBDK() const;
/**
* Return true, if one loop corrections are given in the conventions
* of everything expanded.
*/
virtual bool isExpanded() const;
/**
* Return the value of the dimensional regularization
* parameter. Note that renormalization scale dependence is fully
* restored in DipoleIOperator.
*/
virtual Energy2 mu2() const;
/**
* If defined, return the coefficient of the pole in epsilon^2
*/
virtual double oneLoopDoublePole() const;
/**
* If defined, return the coefficient of the pole in epsilon
*/
virtual double oneLoopSinglePole() const;
/**
* Return true, if cancellationn of epsilon poles should be checked.
*/
bool checkPoles() const;
/**
- * Perform the check of epsilon pole cancellation. Results will be
- * written to the log file, one per phasespace point.
+ * Simple histogram for accuracy checks
+ */
+ struct AccuracyHistogram {
+
+ /**
+ * The lower bound
+ */
+ double lower;
+
+ /**
+ * The bins, indexed by upper bound.
+ */
+ map<double,double> bins;
+
+ /**
+ * The number of points of same sign
+ */
+ unsigned long sameSign;
+
+ /**
+ * The number of points of opposite sign
+ */
+ unsigned long oppositeSign;
+
+ /**
+ * The number of points being nan or inf
+ */
+ unsigned long nans;
+
+ /**
+ * Constructor
+ */
+ AccuracyHistogram(double low = -40.,
+ double up = 0.,
+ unsigned int nbins = 80);
+
+ /**
+ * Book two values to be checked for numerical compatibility
+ */
+ void book(double a, double b);
+
+ /**
+ * Write to file.
+ */
+ void dump(const std::string& prefix,
+ const cPDVector& proc) const;
+
+ /**
+ * Write to persistent ostream
+ */
+ void persistentOutput(PersistentOStream&) const;
+
+ /**
+ * Read from persistent istream
+ */
+ void persistentInput(PersistentIStream&);
+
+ };
+
+ /**
+ * Perform the check of epsilon pole cancellation.
*/
void logPoles() const;
/**
* Return the virtual corrections
*/
const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const {
return theVirtuals;
}
/**
* Return the virtual corrections
*/
vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() {
return theVirtuals;
}
/**
* Instruct this matrix element to include one-loop corrections
*/
void doOneLoop() { theOneLoop = true; }
/**
* Return true, if this matrix element includes one-loop corrections
*/
bool oneLoop() const { return theOneLoop; }
/**
* Instruct this matrix element to include one-loop corrections but
* no Born contributions
*/
void doOneLoopNoBorn() { theOneLoop = true; theOneLoopNoBorn = true; }
/**
* Return true, if this matrix element includes one-loop corrections
* but no Born contributions
*/
bool oneLoopNoBorn() const { return theOneLoopNoBorn || onlyOneLoop(); }
//@}
/** @name Dipole subtraction */
//@{
/**
* If this matrix element is considered a real
* emission matrix element, return all subtraction
* dipoles needed given a set of subtraction terms
* and underlying Born matrix elements to choose
* from.
*/
vector<Ptr<SubtractionDipole>::ptr>
getDipoles(const vector<Ptr<SubtractionDipole>::ptr>&,
const vector<Ptr<MatchboxMEBase>::ptr>&) const;
/**
* If this matrix element is considered a real emission matrix
* element, but actually neglecting a subclass of the contributing
* diagrams, return true if the given emitter-emission-spectator
* configuration should not be considered when setting up
* subtraction dipoles.
*/
virtual bool noDipole(int,int,int) const { return false; }
/**
* If this matrix element is considered an underlying Born matrix
* element in the context of a subtracted real emission, but
* actually neglecting a subclass of the contributing diagrams,
* return true if the given emitter-spectator configuration
* should not be considered when setting up subtraction dipoles.
*/
virtual bool noDipole(int,int) const { return false; }
/**
* Return the colour correlated matrix element squared with
* respect to the given two partons as appearing in mePartonData(),
* suitably scaled by sHat() to give a dimension-less number.
*/
virtual double colourCorrelatedME2(pair<int,int>) const;
/**
* Return the colour correlated matrix element squared in the
* large-N approximation with respect to the given two partons as
* appearing in mePartonData(), suitably scaled by sHat() to give a
* dimension-less number.
*/
virtual double largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr largeNBasis) const;
/**
* Return the colour and spin correlated matrix element squared for
* the gluon indexed by the first argument using the given
* correlation tensor.
*/
virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& c) const;
//@}
/** @name Caching and diagnostic information */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches();
/**
* Return true, if verbose
*/
bool verbose() const;
/**
* Return true, if verbose
*/
bool initVerbose() const;
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Print debug information on the last event
*/
virtual void printLastEvent(ostream&) const;
/**
* Write out diagnostic information for
* generateKinematics
*/
void logGenerateKinematics(const double * r) const;
/**
* Write out diagnostic information for
* setting scales
*/
void logSetScale() const;
/**
* Write out diagnostic information for
* pdf evaluation
*/
void logPDFWeight() const;
/**
* Write out diagnostic information for
* me2 evaluation
*/
void logME2() const;
/**
* Write out diagnostic information
* for dsigdr evaluation
*/
void logDSigHatDR() const;
//@}
/** @name Reweight objects */
//@{
/**
* Insert a reweight object
*/
void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }
/**
* Return the reweights
*/
const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }
/**
* Access the reweights
*/
vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }
//@}
/** @name Methods used to setup MatchboxMEBase objects */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return true; }
/**
* Clone this matrix element.
*/
Ptr<MatchboxMEBase>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
void cloneDependencies(const std::string& prefix = "");
/**
* Prepare an xcomb
*/
void prepareXComb(MatchboxXCombData&) const;
/**
* For the given event generation setup return a xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allPBins,
tStdXCombPtr newHead = tStdXCombPtr(),
tMEPtr newME = tMEPtr());
/**
* For the given event generation setup return a dependent xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME = tMEPtr());
//@}
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();
+
+ /**
+ * Finalize this object. Called in the run phase just after a
+ * run has ended. Used eg. to write out statistics.
+ */
+ virtual void dofinish();
//@}
private:
/**
* The factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tcptr theFactory;
/**
* The phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The amplitude to be used
*/
Ptr<MatchboxAmplitude>::ptr theAmplitude;
/**
* The scale choice object
*/
Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
/**
* The virtual corrections.
*/
vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
/**
* A vector of reweight objects the sum of which
* should be applied to reweight this matrix element
*/
vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
private:
/**
* The subprocesses to be considered.
*/
vector<PDVector> theSubprocesses;
/**
* True, if this matrix element includes one-loop corrections
*/
bool theOneLoop;
/**
* True, if this matrix element includes one-loop corrections
* but no Born contributions
*/
bool theOneLoopNoBorn;
+ /**
+ * Histograms of epsilon^2 pole cancellation
+ */
+ mutable map<cPDVector,AccuracyHistogram> epsilonSquarePoleHistograms;
+
+ /**
+ * Histograms of epsilon pole cancellation
+ */
+ mutable map<cPDVector,AccuracyHistogram> epsilonPoleHistograms;
+
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxMEBase & operator=(const MatchboxMEBase &);
};
+inline PersistentOStream& operator<<(PersistentOStream& os,
+ const MatchboxMEBase::AccuracyHistogram& h) {
+ h.persistentOutput(os);
+ return os;
+}
+
+inline PersistentIStream& operator>>(PersistentIStream& is,
+ MatchboxMEBase::AccuracyHistogram& h) {
+ h.persistentInput(is);
+ return is;
+}
+
+
}
#endif /* HERWIG_MatchboxMEBase_H */
diff --git a/MatrixElement/Matchbox/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc
--- a/MatrixElement/Matchbox/MatchboxFactory.cc
+++ b/MatrixElement/Matchbox/MatchboxFactory.cc
@@ -1,1265 +1,1259 @@
// -*- C++ -*-
//
// MatchboxFactory.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 MatchboxFactory class.
//
#include "MatchboxFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SU2Helper.h"
#include <boost/progress.hpp>
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
using std::ostream_iterator;
MatchboxFactory::MatchboxFactory()
: SubProcessHandler(), theNLight(0),
theOrderInAlphaS(0), theOrderInAlphaEW(0),
theBornContributions(true), theVirtualContributions(true),
theRealContributions(true), theIndependentVirtuals(false),
theSubProcessGroups(false), theInclusive(false),
theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0),
theFixedCouplings(false), theFixedQEDCouplings(false), theVetoScales(false),
- theVerbose(false), theInitVerbose(false), theSubtractionData(""),
- theCheckPoles(false), theRealEmissionScales(false),
- theAllProcesses(false) {}
+ theVerbose(false), theInitVerbose(false), theSubtractionData(""), thePoleData(""),
+ theRealEmissionScales(false), theAllProcesses(false) {}
MatchboxFactory::~MatchboxFactory() {}
IBPtr MatchboxFactory::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxFactory::fullclone() const {
return new_ptr(*this);
}
void MatchboxFactory::prepareME(Ptr<MatchboxMEBase>::ptr me) const {
Ptr<MatchboxAmplitude>::ptr amp =
dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>((*me).amplitude());
me->matchboxAmplitude(amp);
me->factory(this);
if ( phasespace() && !me->phasespace() )
me->phasespace(phasespace());
if ( scaleChoice() && !me->scaleChoice() )
me->scaleChoice(scaleChoice());
}
struct LegIndex {
int spin;
int charge;
int colour;
int isSameAs;
int isSameFamilyAs;
inline bool operator==(const LegIndex& other) const {
return
spin == other.spin &&
charge == other.charge &&
colour == other.colour &&
isSameAs == other.isSameAs &&
isSameFamilyAs == other.isSameFamilyAs;
}
inline bool operator<(const LegIndex& other) const {
if ( spin != other.spin )
return spin < other.spin;
if ( charge != other.charge )
return charge < other.charge;
if ( colour != other.colour )
return colour < other.colour;
if ( isSameAs != other.isSameAs )
return isSameAs < other.isSameAs;
if ( isSameFamilyAs != other.isSameFamilyAs )
return isSameFamilyAs < other.isSameFamilyAs;
return false;
}
};
vector<LegIndex> makeIndex(const PDVector& proc) {
map<long,int> idMap;
map<int,int> familyIdMap;
int lastId = 0;
int lastFamilyId = 0;
vector<LegIndex> res;
for ( PDVector::const_iterator p = proc.begin();
p != proc.end(); ++p ) {
int id;
if ( idMap.find((**p).id()) != idMap.end() ) {
id = idMap[(**p).id()];
} else {
id = lastId;
idMap[(**p).id()] = lastId;
++lastId;
}
int familyId;
if ( familyIdMap.find(SU2Helper::family(*p)) != familyIdMap.end() ) {
familyId = familyIdMap[SU2Helper::family(*p)];
} else {
familyId = lastFamilyId;
familyIdMap[SU2Helper::family(*p)] = lastFamilyId;
++lastFamilyId;
}
LegIndex idx;
idx.spin = (**p).iSpin();
idx.charge = (**p).iCharge();
idx.colour = (**p).iColour();
idx.isSameAs = id;
idx.isSameFamilyAs = familyId;
res.push_back(idx);
}
return res;
}
string pid(const vector<LegIndex>& key) {
ostringstream res;
for ( vector<LegIndex>::const_iterator k =
key.begin(); k != key.end(); ++k )
res << k->spin << k->charge
<< k->colour << k->isSameAs
<< k->isSameFamilyAs;
return res.str();
}
vector<Ptr<MatchboxMEBase>::ptr> MatchboxFactory::
makeMEs(const vector<string>& proc, unsigned int orderas) const {
typedef vector<LegIndex> QNKey;
generator()->log() << "determining subprocesses for ";
copy(proc.begin(),proc.end(),ostream_iterator<string>(generator()->log()," "));
generator()->log() << flush;
map<Ptr<MatchboxAmplitude>::ptr,map<QNKey,vector<PDVector> > > ampProcs;
set<PDVector> processes = makeSubProcesses(proc);
set<PDVector> unsortedProcesses = makeUnsortedSubProcesses(proc);
vector<Ptr<MatchboxAmplitude>::ptr> matchAmplitudes;
unsigned int lowestAsOrder =
allProcesses() ? 0 : orderas;
unsigned int highestAsOrder = orderas;
unsigned int lowestAeOrder =
allProcesses() ? 0 : orderInAlphaEW();
unsigned int highestAeOrder = orderInAlphaEW();
for ( unsigned int oas = lowestAsOrder; oas <= highestAsOrder; ++oas ) {
for ( unsigned int oae = lowestAeOrder; oae <= highestAeOrder; ++oae ) {
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
(**amp).orderInGs(oas);
(**amp).orderInGem(oae);
if ( (**amp).orderInGs() != oas ||
(**amp).orderInGem() != oae ) {
continue;
}
matchAmplitudes.push_back(*amp);
}
}
}
size_t combinations = processes.size()*matchAmplitudes.size()
+ unsortedProcesses.size()*matchAmplitudes.size();
size_t procCount = 0;
boost::progress_display * progressBar =
new boost::progress_display(combinations,generator()->log());
for ( unsigned int oas = lowestAsOrder; oas <= highestAsOrder; ++oas ) {
for ( unsigned int oae = lowestAeOrder; oae <= highestAeOrder; ++oae ) {
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= matchAmplitudes.begin(); amp != matchAmplitudes.end(); ++amp ) {
(**amp).orderInGs(oas);
(**amp).orderInGem(oae);
for ( set<PDVector>::const_iterator p = processes.begin();
p != processes.end(); ++p ) {
++(*progressBar);
if ( !(**amp).canHandle(*p) || !(**amp).sortOutgoing() )
continue;
QNKey key = makeIndex(*p);
++procCount;
ampProcs[*amp][key].push_back(*p);
}
for ( set<PDVector>::const_iterator p = unsortedProcesses.begin();
p != unsortedProcesses.end(); ++p ) {
++(*progressBar);
if ( !(**amp).canHandle(*p) || (**amp).sortOutgoing() )
continue;
QNKey key = makeIndex(*p);
++procCount;
ampProcs[*amp][key].push_back(*p);
}
}
}
}
delete progressBar;
generator()->log() << flush;
vector<Ptr<MatchboxMEBase>::ptr> res;
for ( map<Ptr<MatchboxAmplitude>::ptr,map<QNKey,vector<PDVector> > >::const_iterator
ap = ampProcs.begin(); ap != ampProcs.end(); ++ap ) {
for ( map<QNKey,vector<PDVector> >::const_iterator m = ap->second.begin();
m != ap->second.end(); ++m ) {
Ptr<MatchboxMEBase>::ptr me = ap->first->makeME(m->second);
me->subProcesses() = m->second;
me->amplitude(ap->first);
prepareME(me);
string pname = "ME" + ap->first->name() + pid(m->first);
if ( ! (generator()->preinitRegister(me,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
if ( me->diagrams().empty() )continue;
res.push_back(me);
}
}
generator()->log() << "created " << res.size()
<< " matrix element objects for "
<< procCount << " subprocesses.\n";
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
return res;
}
void MatchboxFactory::setup() {
if ( bornMEs().empty() ) {
if ( particleGroups().find("j") == particleGroups().end() )
throw InitException() << "Could not find a jet particle group named 'j'";
// rebind the particle data objects
for ( map<string,PDVector>::iterator g = particleGroups().begin();
g != particleGroups().end(); ++g )
for ( PDVector::iterator p = g->second.begin();
p != g->second.end(); ++p ) {
#ifndef NDEBUG
long checkid = (**p).id();
#endif
*p = getParticleData((**p).id());
assert((**p).id() == checkid);
}
const PDVector& partons = particleGroups()["j"];
unsigned int nl = 0;
for ( PDVector::const_iterator p = partons.begin();
p != partons.end(); ++p )
if ( abs((**p).id()) < 6 )
++nl;
nLight(nl/2);
vector<Ptr<MatchboxMEBase>::ptr> ames = makeMEs(process,orderInAlphaS());
copy(ames.begin(),ames.end(),back_inserter(bornMEs()));
if ( realContributions() && realEmissionMEs().empty() ) {
vector<string> rproc = process;
if ( realEmissionProcess.empty() ) {
rproc.push_back("j");
} else {
rproc = realEmissionProcess;
}
ames = makeMEs(rproc,orderInAlphaS()+1);
copy(ames.begin(),ames.end(),back_inserter(realEmissionMEs()));
}
}
// check if we have virtual contributions
bool haveVirtuals = true;
// check DR conventions of virtual contributions
bool virtualsAreDR = false;
bool virtualsAreCDR = false;
// check finite term conventions of virtual contributions
bool virtualsAreCS = false;
bool virtualsAreBDK = false;
bool virtualsAreExpanded = false;
// check and prepare the Born and virtual matrix elements
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
prepareME(*born);
haveVirtuals &= (**born).haveOneLoop();
if ( (**born).haveOneLoop() ) {
virtualsAreDR |= (**born).isDR();
virtualsAreCDR |= !(**born).isDR();
virtualsAreCS |= (**born).isCS();
virtualsAreBDK |= (**born).isBDK();
virtualsAreExpanded |= (**born).isExpanded();
}
}
// check the additional insertion operators
if ( !virtuals().empty() )
haveVirtuals = true;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
virtualsAreDR |= (**virt).isDR();
virtualsAreCDR |= !(**virt).isDR();
virtualsAreCS |= (**virt).isCS();
virtualsAreBDK |= (**virt).isBDK();
virtualsAreExpanded |= (**virt).isExpanded();
}
// check for consistent conventions on virtuals, if we are to include them
if ( virtualContributions() ) {
if ( virtualsAreDR && virtualsAreCDR ) {
throw InitException() << "Virtual corrections use inconsistent regularization schemes.\n";
}
if ( (virtualsAreCS && virtualsAreBDK) ||
(virtualsAreCS && virtualsAreExpanded) ||
(virtualsAreBDK && virtualsAreExpanded) ||
(!virtualsAreCS && !virtualsAreBDK && !virtualsAreExpanded) ) {
throw InitException() << "Virtual corrections use inconsistent conventions on finite terms.\n";
}
if ( !haveVirtuals ) {
throw InitException() << "Could not find amplitudes for all virtual contributions needed.\n";
}
}
// prepare dipole insertion operators
if ( virtualContributions() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( virtualsAreDR )
(**virt).useDR();
else
(**virt).useCDR();
if ( virtualsAreCS )
(**virt).useCS();
if ( virtualsAreBDK )
(**virt).useBDK();
if ( virtualsAreExpanded )
(**virt).useExpanded();
}
}
// prepare the real emission matrix elements
if ( realContributions() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
prepareME(*real);
}
}
// start creating matrix elements
MEs().clear();
// setup born and virtual contributions
if ( bornContributions() || virtualContributions() ) {
generator()->log() << "preparing Born"
<< (virtualContributions() ? " and virtual" : "")
<< " matrix elements.\n" << flush;
}
if ( (bornContributions() && !virtualContributions()) ||
(bornContributions() && virtualContributions() && independentVirtuals()) ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
if ( (**born).onlyOneLoop() )
continue;
Ptr<MatchboxMEBase>::ptr bornme = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( independentVirtuals() )
pname += ".Born";
if ( ! (generator()->preinitRegister(bornme,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
bornme->cloneDependencies();
MEs().push_back(bornme);
}
}
if ( virtualContributions() ) {
bornVirtualMEs().clear();
boost::progress_display * progressBar =
new boost::progress_display(bornMEs().size(),generator()->log());
+ if ( thePoleData != "" )
+ if ( thePoleData[thePoleData.size()-1] != '/' )
+ thePoleData += "/";
+
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
Ptr<MatchboxMEBase>::ptr nlo = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( !independentVirtuals() )
pname += ".BornVirtual";
else
pname += ".Virtual";
if ( ! (generator()->preinitRegister(nlo,pname) ) )
throw InitException() << "NLO ME " << pname << " already existing.";
nlo->virtuals().clear();
if ( !nlo->onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
if ( nlo->virtuals().empty() )
throw InitException() << "No insertion operators have been found for "
<< (**born).name() << "\n";
if ( checkPoles() ) {
if ( !virtualsAreExpanded ) {
throw InitException() << "Cannot check epsilon poles if virtuals are not in `expanded' convention.\n";
}
}
}
if ( !bornContributions() || independentVirtuals() ) {
nlo->doOneLoopNoBorn();
} else {
nlo->doOneLoop();
}
nlo->cloneDependencies();
bornVirtualMEs().push_back(nlo);
MEs().push_back(nlo);
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
theSplittingDipoles.clear();
set<cPDVector> bornProcs;
if ( showerApproximation() )
if ( showerApproximation()->needsSplittingGenerator() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born )
for ( MEBase::DiagramVector::const_iterator d = (**born).diagrams().begin();
d != (**born).diagrams().end(); ++d )
bornProcs.insert((**d).partons());
}
if ( realContributions() ) {
generator()->log() << "preparing real emission matrix elements.\n" << flush;
if ( theSubtractionData != "" )
if ( theSubtractionData[theSubtractionData.size()-1] != '/' )
theSubtractionData += "/";
subtractedMEs().clear();
boost::progress_display * progressBar =
new boost::progress_display(realEmissionMEs().size(),generator()->log());
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
Ptr<SubtractedME>::ptr sub = new_ptr(SubtractedME());
string pname = fullName() + "/" + (**real).name() + ".SubtractedReal";
if ( ! (generator()->preinitRegister(sub,pname) ) )
throw InitException() << "Subtracted ME " << pname << " already existing.";
sub->factory(this);
sub->head(*real);
sub->dependent().clear();
sub->getDipoles();
if ( sub->dependent().empty() ) {
// finite real contribution
Ptr<MatchboxMEBase>::ptr fme =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(sub->head())->cloneMe();
string qname = fullName() + "/" + (**real).name() + ".FiniteReal";
if ( ! (generator()->preinitRegister(fme,qname) ) )
throw InitException() << "ME " << qname << " already existing.";
MEs().push_back(fme);
finiteRealMEs().push_back(fme);
sub->head(tMEPtr());
continue;
}
if ( realEmissionScales() )
sub->doRealEmissionScales();
subtractedMEs().push_back(sub);
MEs().push_back(sub);
if ( showerApproximation() ) {
if ( virtualContributions() ) {
Ptr<SubtractedME>::ptr subv = new_ptr(*sub);
string vname = sub->fullName() + ".SubtractionIntegral";
if ( ! (generator()->preinitRegister(subv,vname) ) )
throw InitException() << "Subtracted ME " << vname << " already existing.";
subv->cloneDependencies(vname);
subv->doVirtualShowerSubtraction();
subtractedMEs().push_back(subv);
MEs().push_back(subv);
}
sub->doRealShowerSubtraction();
if ( showerApproximation()->needsSplittingGenerator() )
for ( set<cPDVector>::const_iterator p = bornProcs.begin();
p != bornProcs.end(); ++p ) {
vector<Ptr<SubtractionDipole>::ptr> sdip = sub->splitDipoles(*p);
set<Ptr<SubtractionDipole>::ptr>& dips = theSplittingDipoles[*p];
copy(sdip.begin(),sdip.end(),inserter(dips,dips.begin()));
}
}
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
if ( !theSplittingDipoles.empty() ) {
map<Ptr<SubtractionDipole>::ptr,Ptr<SubtractionDipole>::ptr> cloneMap;
for ( map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::const_iterator sd = theSplittingDipoles.begin();
sd != theSplittingDipoles.end(); ++sd ) {
for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator d = sd->second.begin();
d != sd->second.end(); ++d ) {
cloneMap[*d] = Ptr<SubtractionDipole>::ptr();
}
}
for ( map<Ptr<SubtractionDipole>::ptr,Ptr<SubtractionDipole>::ptr>::iterator cd =
cloneMap.begin(); cd != cloneMap.end(); ++cd ) {
Ptr<SubtractionDipole>::ptr cloned = cd->first->cloneMe();
string dname = cd->first->fullName() + ".splitting";
if ( ! (generator()->preinitRegister(cloned,dname)) )
throw InitException() << "Dipole '" << dname << "' already existing.";
cloned->cloneDependencies();
cloned->showerApproximation(Ptr<ShowerApproximation>::tptr());
cloned->doSplitting();
cd->second = cloned;
}
for ( map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::iterator sd = theSplittingDipoles.begin();
sd != theSplittingDipoles.end(); ++sd ) {
set<Ptr<SubtractionDipole>::ptr> cloned;
for ( set<Ptr<SubtractionDipole>::ptr>::iterator d = sd->second.begin();
d != sd->second.end(); ++d ) {
cloned.insert(cloneMap[*d]);
}
sd->second = cloned;
}
}
generator()->log() << "Process setup finished.\n" << flush;
}
void MatchboxFactory::SplittingChannel::print(ostream& os) const {
os << "--- SplittingChannel setup -----------------------------------------------------\n";
os << " Born process ";
const StandardXComb& bxc = *bornXComb;
os << bxc.mePartonData()[0]->PDGName() << " "
<< bxc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = bxc.mePartonData().begin() + 2;
p != bxc.mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "\n";
os << " to real emission process ";
const StandardXComb& rxc = *realXComb;
os << rxc.mePartonData()[0]->PDGName() << " "
<< rxc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = rxc.mePartonData().begin() + 2;
p != rxc.mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "\n";
os << " with dipole:\n";
dipole->print(os);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
list<MatchboxFactory::SplittingChannel>
MatchboxFactory::getSplittingChannels(tStdXCombPtr xcptr) const {
if ( xcptr->lastProjector() )
xcptr = xcptr->lastProjector();
const StandardXComb& xc = *xcptr;
cPDVector proc = xc.mePartonData();
map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::const_iterator splitEntries
= splittingDipoles().find(proc);
list<SplittingChannel> res;
if ( splitEntries == splittingDipoles().end() )
return res;
const set<Ptr<SubtractionDipole>::ptr>& splitDipoles = splitEntries->second;
SplittingChannel channel;
if ( !splitDipoles.empty() ) {
Ptr<MatchboxMEBase>::tptr bornME =
const_ptr_cast<Ptr<MatchboxMEBase>::tptr>((**splitDipoles.begin()).underlyingBornME());
channel.bornXComb =
bornME->makeXComb(xc.maxEnergy(),xc.particles(),xc.eventHandlerPtr(),
const_ptr_cast<tSubHdlPtr>(xc.subProcessHandler()),
xc.pExtractor(),xc.CKKWHandler(),
xc.partonBins(),xc.cuts(),xc.diagrams(),xc.mirror(),
PartonPairVec());
}
for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator sd =
splitDipoles.begin(); sd != splitDipoles.end(); ++sd ) {
channel.dipole = *sd;
vector<StdXCombPtr> realXCombs = (**sd).makeRealXCombs(channel.bornXComb);
for ( vector<StdXCombPtr>::const_iterator rxc = realXCombs.begin();
rxc != realXCombs.end(); ++rxc ) {
channel.realXComb = *rxc;
if ( showerApproximation()->needsTildeXCombs() ) {
channel.tildeXCombs.clear();
assert(!channel.dipole->partnerDipoles().empty());
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator p =
channel.dipole->partnerDipoles().begin();
p != channel.dipole->partnerDipoles().end(); ++p ) {
StdXCombPtr txc = channel.dipole->makeBornXComb(channel.realXComb);
if ( txc )
channel.tildeXCombs.push_back(txc);
}
}
res.push_back(channel);
}
}
if ( initVerbose() ) {
generator()->log()
<< "--- MatchboxFactory splitting channels ----------------------------------------------\n";
const StandardXComb& bxc = *xcptr;
generator()->log() << " hard process handled is: ";
generator()->log() << bxc.mePartonData()[0]->PDGName() << " "
<< bxc.mePartonData()[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = bxc.mePartonData().begin() + 2;
p != bxc.mePartonData().end(); ++p ) {
generator()->log() << (**p).PDGName() << " ";
}
generator()->log() << "\n";
for ( list<MatchboxFactory::SplittingChannel>::const_iterator sp =
res.begin(); sp != res.end(); ++sp ) {
sp->print(generator()->log());
}
generator()->log()
<< "-------------------------------------------------------------------------------------\n"
<< flush;
}
return res;
}
void MatchboxFactory::print(ostream& os) const {
os << "--- MatchboxFactory setup -----------------------------------------------------------\n";
if ( !amplitudes().empty() ) {
os << " generated Born matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = bornMEs().begin();
m != bornMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
os << " generated real emission matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = realEmissionMEs().begin();
m != realEmissionMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
}
os << " generated Born+virtual matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator bv
= bornVirtualMEs().begin(); bv != bornVirtualMEs().end(); ++bv ) {
(**bv).print(os);
}
os << " generated subtracted matrix elements:\n";
for ( vector<Ptr<SubtractedME>::ptr>::const_iterator sub
= subtractedMEs().begin(); sub != subtractedMEs().end(); ++sub ) {
os << " '" << (**sub).name() << "'\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxFactory::doinit() {
setup();
if ( initVerbose() )
print(Repository::clog());
SubProcessHandler::doinit();
}
void MatchboxFactory::persistentOutput(PersistentOStream & os) const {
os << theDiagramGenerator << theProcessData
<< theNLight << theOrderInAlphaS << theOrderInAlphaEW
<< theBornContributions << theVirtualContributions
<< theRealContributions << theIndependentVirtuals << theSubProcessGroups << theInclusive
<< thePhasespace << theScaleChoice
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theFixedCouplings << theFixedQEDCouplings << theVetoScales
<< theAmplitudes
<< theBornMEs << theVirtuals << theRealEmissionMEs
<< theBornVirtualMEs << theSubtractedMEs << theFiniteRealMEs
- << theVerbose << theInitVerbose << theSubtractionData << theCheckPoles
+ << theVerbose << theInitVerbose << theSubtractionData << thePoleData
<< theParticleGroups << process << realEmissionProcess
<< theShowerApproximation << theSplittingDipoles
<< theRealEmissionScales << theAllProcesses;
}
void MatchboxFactory::persistentInput(PersistentIStream & is, int) {
is >> theDiagramGenerator >> theProcessData
>> theNLight >> theOrderInAlphaS >> theOrderInAlphaEW
>> theBornContributions >> theVirtualContributions
>> theRealContributions >> theIndependentVirtuals >> theSubProcessGroups >> theInclusive
>> thePhasespace >> theScaleChoice
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theFixedCouplings >> theFixedQEDCouplings >> theVetoScales
>> theAmplitudes
>> theBornMEs >> theVirtuals >> theRealEmissionMEs
>> theBornVirtualMEs >> theSubtractedMEs >> theFiniteRealMEs
- >> theVerbose >> theInitVerbose >> theSubtractionData >> theCheckPoles
+ >> theVerbose >> theInitVerbose >> theSubtractionData >> thePoleData
>> theParticleGroups >> process >> realEmissionProcess
>> theShowerApproximation >> theSplittingDipoles
>> theRealEmissionScales >> theAllProcesses;
}
string MatchboxFactory::startParticleGroup(string name) {
particleGroupName = StringUtils::stripws(name);
particleGroup.clear();
return "";
}
string MatchboxFactory::endParticleGroup(string) {
if ( particleGroup.empty() )
throw InitException() << "Empty particle group.";
particleGroups()[particleGroupName] = particleGroup;
particleGroup.clear();
return "";
}
string MatchboxFactory::doProcess(string in) {
process = StringUtils::split(in);
if ( process.size() < 3 )
throw InitException() << "Invalid process.";
for ( vector<string>::iterator p = process.begin();
p != process.end(); ++p ) {
*p = StringUtils::stripws(*p);
}
return "";
}
string MatchboxFactory::doSingleRealProcess(string in) {
realEmissionProcess = StringUtils::split(in);
if ( realEmissionProcess.size() < 3 )
throw InitException() << "Invalid process.";
for ( vector<string>::iterator p = realEmissionProcess.begin();
p != realEmissionProcess.end(); ++p ) {
*p = StringUtils::stripws(*p);
}
return "";
}
struct SortPID {
inline bool operator()(PDPtr a, PDPtr b) const {
return a->id() < b->id();
}
};
set<PDVector> MatchboxFactory::
makeSubProcesses(const vector<string>& proc, bool sorted) const {
if ( proc.empty() )
throw InitException() << "No process specified.";
vector<PDVector> allProcs(1);
size_t pos = 0;
typedef map<string,PDVector>::const_iterator GroupIterator;
while ( pos < proc.size() ) {
GroupIterator git =
particleGroups().find(proc[pos]);
if ( git == particleGroups().end() ) {
throw InitException() << "particle group '"
<< proc[pos] << "' not defined.";
}
vector<PDVector> mine;
for ( vector<PDVector>::const_iterator i = allProcs.begin();
i != allProcs.end(); ++i ) {
for ( PDVector::const_iterator p = git->second.begin();
p != git->second.end(); ++p ) {
PDVector v = *i;
v.push_back(*p);
mine.push_back(v);
}
}
allProcs = mine;
++pos;
}
set<PDVector> allCheckedProcs;
for ( vector<PDVector>::const_iterator p = allProcs.begin();
p != allProcs.end(); ++p ) {
int charge = -(*p)[0]->iCharge() -(*p)[1]->iCharge();
for ( size_t k = 2; k < (*p).size(); ++k )
charge += (*p)[k]->iCharge();
if ( charge != 0 )
continue;
PDVector pr = *p;
if ( sorted )
sort(pr.begin()+2,pr.end(),SortPID());
allCheckedProcs.insert(pr);
}
return allCheckedProcs;
}
void MatchboxFactory::Init() {
static ClassDocumentation<MatchboxFactory> documentation
("MatchboxFactory",
"NLO QCD corrections have been calculated "
"using Matchbox \\cite{Platzer:2011bc}",
"%\\cite{Platzer:2011bc}\n"
"\\bibitem{Platzer:2011bc}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Dipole Showers and Automated NLO Matching in Herwig++,''\n"
"arXiv:1109.6256 [hep-ph].\n"
"%%CITATION = ARXIV:1109.6256;%%");
static Reference<MatchboxFactory,Tree2toNGenerator> interfaceDiagramGenerator
("DiagramGenerator",
"Set the diagram generator.",
&MatchboxFactory::theDiagramGenerator, false, false, true, true, false);
static Reference<MatchboxFactory,ProcessData> interfaceProcessData
("ProcessData",
"Set the process data object to be used.",
&MatchboxFactory::theProcessData, false, false, true, true, false);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaS
("OrderInAlphaS",
"The order in alpha_s to consider.",
&MatchboxFactory::theOrderInAlphaS, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaEW
("OrderInAlphaEW",
"The order in alpha_EW",
&MatchboxFactory::theOrderInAlphaEW, 2, 0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceBornContributions
("BornContributions",
"Switch on or off the Born contributions.",
&MatchboxFactory::theBornContributions, true, false, false);
static SwitchOption interfaceBornContributionsOn
(interfaceBornContributions,
"On",
"Switch on Born contributions.",
true);
static SwitchOption interfaceBornContributionsOff
(interfaceBornContributions,
"Off",
"Switch off Born contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceVirtualContributions
("VirtualContributions",
"Switch on or off the virtual contributions.",
&MatchboxFactory::theVirtualContributions, true, false, false);
static SwitchOption interfaceVirtualContributionsOn
(interfaceVirtualContributions,
"On",
"Switch on virtual contributions.",
true);
static SwitchOption interfaceVirtualContributionsOff
(interfaceVirtualContributions,
"Off",
"Switch off virtual contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceRealContributions
("RealContributions",
"Switch on or off the real contributions.",
&MatchboxFactory::theRealContributions, true, false, false);
static SwitchOption interfaceRealContributionsOn
(interfaceRealContributions,
"On",
"Switch on real contributions.",
true);
static SwitchOption interfaceRealContributionsOff
(interfaceRealContributions,
"Off",
"Switch off real contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceIndependentVirtuals
("IndependentVirtuals",
"Switch on or off virtual contributions as separate subprocesses.",
&MatchboxFactory::theIndependentVirtuals, true, false, false);
static SwitchOption interfaceIndependentVirtualsOn
(interfaceIndependentVirtuals,
"On",
"Switch on virtual contributions as separate subprocesses.",
true);
static SwitchOption interfaceIndependentVirtualsOff
(interfaceIndependentVirtuals,
"Off",
"Switch off virtual contributions as separate subprocesses.",
false);
static Switch<MatchboxFactory,bool> interfaceSubProcessGroups
("SubProcessGroups",
"Switch on or off production of sub-process groups.",
&MatchboxFactory::theSubProcessGroups, false, false, false);
static SwitchOption interfaceSubProcessGroupsOn
(interfaceSubProcessGroups,
"On",
"On",
true);
static SwitchOption interfaceSubProcessGroupsOff
(interfaceSubProcessGroups,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceInclusive
("Inclusive",
"Switch on or off production of inclusive cross section.",
&MatchboxFactory::theInclusive, false, false, false);
static SwitchOption interfaceInclusiveOn
(interfaceInclusive,
"On",
"On",
true);
static SwitchOption interfaceInclusiveOff
(interfaceInclusive,
"Off",
"Off",
false);
static Reference<MatchboxFactory,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"Set the phasespace generator.",
&MatchboxFactory::thePhasespace, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxScaleChoice> interfaceScaleChoice
("ScaleChoice",
"Set the scale choice object.",
&MatchboxFactory::theScaleChoice, false, false, true, true, false);
static Parameter<MatchboxFactory,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&MatchboxFactory::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&MatchboxFactory::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceFixedCouplings
("FixedCouplings",
"Switch on or off fixed couplings.",
&MatchboxFactory::theFixedCouplings, true, false, false);
static SwitchOption interfaceFixedCouplingsOn
(interfaceFixedCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedCouplingsOff
(interfaceFixedCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceFixedQEDCouplings
("FixedQEDCouplings",
"Switch on or off fixed QED couplings.",
&MatchboxFactory::theFixedQEDCouplings, true, false, false);
static SwitchOption interfaceFixedQEDCouplingsOn
(interfaceFixedQEDCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedQEDCouplingsOff
(interfaceFixedQEDCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceVetoScales
("VetoScales",
"Switch on or setting veto scales.",
&MatchboxFactory::theVetoScales, false, false, false);
static SwitchOption interfaceVetoScalesOn
(interfaceVetoScales,
"On",
"On",
true);
static SwitchOption interfaceVetoScalesOff
(interfaceVetoScales,
"Off",
"Off",
false);
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceAmplitudes
("Amplitudes",
"The amplitude objects.",
&MatchboxFactory::theAmplitudes, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceBornMEs
("BornMEs",
"The Born matrix elements to be used",
&MatchboxFactory::theBornMEs, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxInsertionOperator> interfaceVirtuals
("Virtuals",
"The virtual corrections to include",
&MatchboxFactory::theVirtuals, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceRealEmissionMEs
("RealEmissionMEs",
"The RealEmission matrix elements to be used",
&MatchboxFactory::theRealEmissionMEs, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceBornVirtuals
("BornVirtualMEs",
"The generated Born/virtual contributions",
&MatchboxFactory::theBornVirtualMEs, -1, false, true, true, true, false);
static RefVector<MatchboxFactory,SubtractedME> interfaceSubtractedMEs
("SubtractedMEs",
"The generated subtracted real emission contributions",
&MatchboxFactory::theSubtractedMEs, -1, false, true, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceFiniteRealMEs
("FiniteRealMEs",
"The generated finite real contributions",
&MatchboxFactory::theFiniteRealMEs, -1, false, true, true, true, false);
static Switch<MatchboxFactory,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&MatchboxFactory::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"On",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceInitVerbose
("InitVerbose",
"Print setup information.",
&MatchboxFactory::theInitVerbose, false, false, false);
static SwitchOption interfaceInitVerboseOn
(interfaceInitVerbose,
"On",
"On",
true);
static SwitchOption interfaceInitVerboseOff
(interfaceInitVerbose,
"Off",
"Off",
false);
static Parameter<MatchboxFactory,string> interfaceSubtractionData
("SubtractionData",
"Prefix for subtraction check data.",
&MatchboxFactory::theSubtractionData, "",
false, false);
- static Switch<MatchboxFactory,bool> interfaceCheckPoles
- ("CheckPoles",
- "Switch on or off checks of epsilon poles.",
- &MatchboxFactory::theCheckPoles, true, false, false);
- static SwitchOption interfaceCheckPolesOn
- (interfaceCheckPoles,
- "On",
- "On",
- true);
- static SwitchOption interfaceCheckPolesOff
- (interfaceCheckPoles,
- "Off",
- "Off",
- false);
+ static Parameter<MatchboxFactory,string> interfacePoleData
+ ("PoleData",
+ "Prefix for subtraction check data.",
+ &MatchboxFactory::thePoleData, "",
+ false, false);
static RefVector<MatchboxFactory,ParticleData> interfaceParticleGroup
("ParticleGroup",
"The particle group just started.",
&MatchboxFactory::particleGroup, -1, false, false, true, false, false);
static Command<MatchboxFactory> interfaceStartParticleGroup
("StartParticleGroup",
"Start a particle group.",
&MatchboxFactory::startParticleGroup, false);
static Command<MatchboxFactory> interfaceEndParticleGroup
("EndParticleGroup",
"End a particle group.",
&MatchboxFactory::endParticleGroup, false);
static Command<MatchboxFactory> interfaceProcess
("Process",
"Set the process to consider.",
&MatchboxFactory::doProcess, false);
static Command<MatchboxFactory> interfaceSingleRealProcess
("SingleRealProcess",
"Set the process to consider.",
&MatchboxFactory::doSingleRealProcess, false);
static Reference<MatchboxFactory,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&MatchboxFactory::theShowerApproximation, false, false, true, true, false);
static Switch<MatchboxFactory,bool> interfaceRealEmissionScales
("RealEmissionScales",
"Switch on or off calculation of subtraction scales from real emission kinematics.",
&MatchboxFactory::theRealEmissionScales, false, false, false);
static SwitchOption interfaceRealEmissionScalesOn
(interfaceRealEmissionScales,
"On",
"On",
true);
static SwitchOption interfaceRealEmissionScalesOff
(interfaceRealEmissionScales,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceAllProcesses
("AllProcesses",
"Consider all processes up to a maximum coupling order specified by the coupling order interfaces.",
&MatchboxFactory::theAllProcesses, false, false, false);
static SwitchOption interfaceAllProcessesYes
(interfaceAllProcesses,
"Yes",
"Include all processes.",
true);
static SwitchOption interfaceAllProcessesNo
(interfaceAllProcesses,
"No",
"Only consider processes matching the exact order in the couplings.",
false);
}
// *** 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<MatchboxFactory,SubProcessHandler>
describeHerwigMatchboxFactory("Herwig::MatchboxFactory", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/MatchboxFactory.h b/MatrixElement/Matchbox/MatchboxFactory.h
--- a/MatrixElement/Matchbox/MatchboxFactory.h
+++ b/MatrixElement/Matchbox/MatchboxFactory.h
@@ -1,809 +1,814 @@
// -*- C++ -*-
//
// MatchboxFactory.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_MatchboxFactory_H
#define HERWIG_MatchboxFactory_H
//
// This is the declaration of the MatchboxFactory class.
//
#include "ThePEG/Handlers/SubProcessHandler.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/ProcessData.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Base/SubtractedME.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxFactory automatically sets up a NLO
* QCD calculation carried out in dipole subtraction.
*
* @see \ref MatchboxFactoryInterfaces "The interfaces"
* defined for MatchboxFactory.
*/
class MatchboxFactory: public SubProcessHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxFactory();
/**
* The destructor.
*/
virtual ~MatchboxFactory();
//@}
public:
/** @name Process and diagram information */
//@{
/**
* Return the diagram generator.
*/
Ptr<Tree2toNGenerator>::tptr diagramGenerator() const { return theDiagramGenerator; }
/**
* Set the diagram generator.
*/
void diagramGenerator(Ptr<Tree2toNGenerator>::ptr dg) { theDiagramGenerator = dg; }
/**
* Return the process data.
*/
Ptr<ProcessData>::tptr processData() const { return theProcessData; }
/**
* Set the process data.
*/
void processData(Ptr<ProcessData>::ptr pd) { theProcessData = pd; }
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
unsigned int nLight() const { return theNLight; }
/**
* Set the number of light flavours, this matrix
* element is calculated for.
*/
void nLight(unsigned int n) { theNLight = n; }
/**
* Return the order in \f$\alpha_S\f$.
*/
unsigned int orderInAlphaS() const { return theOrderInAlphaS; }
/**
* Set the order in \f$\alpha_S\f$.
*/
void orderInAlphaS(unsigned int o) { theOrderInAlphaS = o; }
/**
* Return the order in \f$\alpha_{EM}\f$.
*/
unsigned int orderInAlphaEW() const { return theOrderInAlphaEW; }
/**
* Set the order in \f$\alpha_{EM}\f$.
*/
void orderInAlphaEW(unsigned int o) { theOrderInAlphaEW = o; }
/**
* Return true, if all processes up to a maximum order are considered
*/
bool allProcesses() const { return theAllProcesses; }
/**
* Switch on/off inclusino off all processes up to a maximum order
*/
void setAllProcesses(bool on = true) { theAllProcesses = on; }
/**
* Return true, if Born contributions should be included.
*/
bool bornContributions() const { return theBornContributions; }
/**
* Switch on or off Born contributions
*/
void setBornContributions(bool on = true) { theBornContributions = on; }
/**
* Return true, if virtual contributions should be included.
*/
bool virtualContributions() const { return theVirtualContributions; }
/**
* Switch on or off virtual contributions
*/
void setVirtualContributions(bool on = true) { theVirtualContributions = on; }
/**
* Return true, if subtracted real emission contributions should be included.
*/
bool realContributions() const { return theRealContributions; }
/**
* Switch on or off subtracted real emission contributions
*/
void setRealContributions(bool on = true) { theRealContributions = on; }
/**
* Return true, if virtual contributions should be treated as independent subprocesses
*/
bool independentVirtuals() const { return theIndependentVirtuals; }
/**
* Switch on/off virtual contributions should be treated as independent subprocesses
*/
void setIndependentVirtuals(bool on = true) { theIndependentVirtuals = on; }
/**
* Return true, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool subProcessGroups() const { return theSubProcessGroups; }
/**
* Switch on or off producing subprocess groups.
*/
void setSubProcessGroups(bool on = true) { theSubProcessGroups = on; }
/**
* Return true, if subtraction scales should be caluclated from real emission kinematics
*/
bool realEmissionScales() const { return theRealEmissionScales; }
/**
* Switch on/off that subtraction scales should be caluclated from real emission kinematics
*/
void setRealEmissionScales(bool on = true) { theRealEmissionScales = on; }
/**
* Return true, if the integral over the unresolved emission should be
* calculated.
*/
bool inclusive() const { return theInclusive; }
/**
* Switch on or off inclusive mode.
*/
void setInclusive(bool on = true) { theInclusive = on; }
/**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr app) { theShowerApproximation = app; }
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
//@}
/** @name Phasespace generation and scale choice */
//@{
/**
* Return the phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::tptr phasespace() const { return thePhasespace; }
/**
* Set the phase space generator to be used.
*/
void phasespace(Ptr<MatchboxPhasespace>::ptr ps) { thePhasespace = ps; }
/**
* Set the scale choice object
*/
void scaleChoice(Ptr<MatchboxScaleChoice>::ptr sc) { theScaleChoice = sc; }
/**
* Return the scale choice object
*/
Ptr<MatchboxScaleChoice>::tptr scaleChoice() const { return theScaleChoice; }
/**
* Get the factorization scale factor
*/
double factorizationScaleFactor() const { return theFactorizationScaleFactor; }
/**
* Set the factorization scale factor
*/
void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }
/**
* Get the renormalization scale factor
*/
double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }
/**
* Set the renormalization scale factor
*/
void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }
/**
* Return true, if fixed couplings are used.
*/
bool fixedCouplings() const { return theFixedCouplings; }
/**
* Switch on fixed couplings.
*/
void setFixedCouplings(bool on = true) { theFixedCouplings = on; }
/**
* Return true, if fixed couplings are used.
*/
bool fixedQEDCouplings() const { return theFixedQEDCouplings; }
/**
* Switch on fixed couplings.
*/
void setFixedQEDCouplings(bool on = true) { theFixedQEDCouplings = on; }
/**
* Return true, if veto scales should be set
* for the real emission
*/
bool vetoScales() const { return theVetoScales; }
/**
* Switch on setting veto scales
*/
void doVetoScales() { theVetoScales = true; }
/**
* Switch off setting veto scales
*/
void noVetoScales() { theVetoScales = true; }
//@}
/** @name Amplitudes and caching */
//@{
/**
* Return the amplitudes to be considered
*/
const vector<Ptr<MatchboxAmplitude>::ptr>& amplitudes() const { return theAmplitudes; }
/**
* Access the amplitudes to be considered
*/
vector<Ptr<MatchboxAmplitude>::ptr>& amplitudes() { return theAmplitudes; }
//@}
/** @name Matrix element objects. */
//@{
/**
* Return the Born matrix elements to be considered
*/
const vector<Ptr<MatchboxMEBase>::ptr>& bornMEs() const { return theBornMEs; }
/**
* Access the Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr>& bornMEs() { return theBornMEs; }
/**
* Return the virtual corrections to be considered
*/
const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const { return theVirtuals; }
/**
* Access the virtual corrections to be considered
*/
vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() { return theVirtuals; }
/**
* Return the produced NLO matrix elements
*/
const vector<Ptr<MatchboxMEBase>::ptr>& bornVirtualMEs() const { return theBornVirtualMEs; }
/**
* Access the produced NLO matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr>& bornVirtualMEs() { return theBornVirtualMEs; }
/**
* Return the real emission matrix elements to be considered
*/
const vector<Ptr<MatchboxMEBase>::ptr>& realEmissionMEs() const { return theRealEmissionMEs; }
/**
* Access the real emission matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr>& realEmissionMEs() { return theRealEmissionMEs; }
/**
* Return the produced subtracted matrix elements
*/
const vector<Ptr<SubtractedME>::ptr>& subtractedMEs() const { return theSubtractedMEs; }
/**
* Access the produced subtracted matrix elements
*/
vector<Ptr<SubtractedME>::ptr>& subtractedMEs() { return theSubtractedMEs; }
/**
* Return the produced finite real emission matrix elements
*/
const vector<Ptr<MatchboxMEBase>::ptr>& finiteRealMEs() const { return theFiniteRealMEs; }
/**
* Access the produced finite real emission elements
*/
vector<Ptr<MatchboxMEBase>::ptr>& finiteRealMEs() { return theFiniteRealMEs; }
/**
* Return the map of Born processes to splitting dipoles
*/
const map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >& splittingDipoles() const {
return theSplittingDipoles;
}
/**
* Identify a splitting channel
*/
struct SplittingChannel {
/**
* The Born XComb
*/
StdXCombPtr bornXComb;
/**
* The real XComb
*/
StdXCombPtr realXComb;
/**
* The set of tilde XCombs to consider for the real xcomb
*/
vector<StdXCombPtr> tildeXCombs;
/**
* The dipole in charge of the splitting
*/
Ptr<SubtractionDipole>::ptr dipole;
/**
* Dump the setup
*/
void print(ostream&) const;
};
/**
* Generate all splitting channels for the Born process handled by
* the given XComb
*/
list<SplittingChannel> getSplittingChannels(tStdXCombPtr xc) const;
//@}
/** @name Setup the matrix elements */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return true; }
/**
* Prepare a matrix element.
*/
void prepareME(Ptr<MatchboxMEBase>::ptr) const;
/**
* Setup everything
*/
virtual void setup();
//@}
/** @name Diagnostic information */
//@{
/**
* Return true, if verbose
*/
bool verbose() const { return theVerbose; }
/**
* Switch on diagnostic information.
*/
void setVerbose(bool on = true) { theVerbose = on; }
/**
* Return true, if verbose while initializing
*/
bool initVerbose() const { return theInitVerbose || verbose(); }
/**
* Switch on diagnostic information while initializing
*/
void setInitVerbose(bool on = true) { theInitVerbose = on; }
/**
* Dump the setup
*/
void print(ostream&) const;
/**
* Return the subtraction data prefix.
*/
const string& subtractionData() const { return theSubtractionData; }
/**
* Set the subtraction data prefix.
*/
void subtractionData(const string& s) { theSubtractionData = s; }
/**
+ * Return the pole data prefix.
+ */
+ const string& poleData() const { return thePoleData; }
+
+ /**
+ * Set the pole data prefix.
+ */
+ void poleData(const string& s) { thePoleData = s; }
+
+ /**
* Return true, if cancellationn of epsilon poles should be checked.
*/
- bool checkPoles() const { return theCheckPoles; }
-
- /**
- * Switch on checking of epsilon pole cancellation.
- */
- void doCheckPoles() { theCheckPoles = true; }
+ bool checkPoles() const { return poleData() != ""; }
//@}
/** @name Process generation */
//@{
/**
* Return the particle groups.
*/
const map<string,PDVector>& particleGroups() const { return theParticleGroups; }
/**
* Access the particle groups.
*/
map<string,PDVector>& particleGroups() { return theParticleGroups; }
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The diagram generator.
*/
Ptr<Tree2toNGenerator>::ptr theDiagramGenerator;
/**
* The process data object to be used
*/
Ptr<ProcessData>::ptr theProcessData;
/**
* The number of light flavours, this matrix
* element is calculated for.
*/
unsigned int theNLight;
/**
* The order in \f$\alpha_S\f$.
*/
unsigned int theOrderInAlphaS;
/**
* The order in \f$\alpha_{EM}\f$.
*/
unsigned int theOrderInAlphaEW;
/**
* Switch on or off Born contributions
*/
bool theBornContributions;
/**
* Switch on or off virtual contributions
*/
bool theVirtualContributions;
/**
* Switch on or off subtracted real emission contributions should be included.
*/
bool theRealContributions;
/**
* True if virtual contributions should be treated as independent subprocesses
*/
bool theIndependentVirtuals;
/**
* True, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool theSubProcessGroups;
/**
* True, if the integral over the unresolved emission should be
* calculated.
*/
bool theInclusive;
/**
* The phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The scale choice object
*/
Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
/**
* Use non-running couplings.
*/
bool theFixedCouplings;
/**
* Use non-running couplings.
*/
bool theFixedQEDCouplings;
/**
* True, if veto scales should be set
* for the real emission
*/
bool theVetoScales;
/**
* The amplitudes to be considered
*/
vector<Ptr<MatchboxAmplitude>::ptr> theAmplitudes;
/**
* The Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theBornMEs;
/**
* The virtual corrections to be considered
*/
vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
/**
* The real emission matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theRealEmissionMEs;
/**
* The produced NLO matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr> theBornVirtualMEs;
/**
* The produced subtracted matrix elements
*/
vector<Ptr<SubtractedME>::ptr> theSubtractedMEs;
/**
* The produced finite real emission matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr> theFiniteRealMEs;
/**
* Switch on or off verbosity
*/
bool theVerbose;
/**
* True, if verbose while initializing
*/
bool theInitVerbose;
/**
* Prefix for subtraction data
*/
string theSubtractionData;
/**
+ * Prefix for pole data.
+ */
+ string thePoleData;
+
+ /**
* Command to limit the real emission process to be considered.
*/
string doSingleRealProcess(string);
/**
* The real emission process to be included; if empty, all possible
* ones will be considered.
*/
vector<string> realEmissionProcess;
/**
- * True, if cancellationn of epsilon poles should be checked.
- */
- bool theCheckPoles;
-
- /**
* Particle groups.
*/
map<string,PDVector> theParticleGroups;
/**
* Command to start a particle group.
*/
string startParticleGroup(string);
/**
* The name of the particle group currently edited.
*/
string particleGroupName;
/**
* The particle group currently edited.
*/
PDVector particleGroup;
/**
* Command to end a particle group.
*/
string endParticleGroup(string);
/**
* Command to set the process.
*/
string doProcess(string);
/**
* The process to consider in terms of particle groups.
*/
vector<string> process;
/**
* Generate subprocesses.
*/
set<PDVector> makeSubProcesses(const vector<string>&, bool sorted = true) const;
/**
* Generate subprocesses with all permutations of outgoing partons.
*/
set<PDVector> makeUnsortedSubProcesses(const vector<string>& data) const {
return makeSubProcesses(data, false);
}
/**
* Generate matrix element objects for the given process.
*/
vector<Ptr<MatchboxMEBase>::ptr> makeMEs(const vector<string>&,
unsigned int orderas) const;
/**
* The shower approximation.
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* The map of Born processes to splitting dipoles
*/
map<cPDVector,set<Ptr<SubtractionDipole>::ptr> > theSplittingDipoles;
/**
* True, if subtraction scales should be caluclated from real emission kinematics
*/
bool theRealEmissionScales;
/**
* Consider all processes with order in couplings specifying the
* maximum order.
*/
bool theAllProcesses;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxFactory & operator=(const MatchboxFactory &);
};
}
#endif /* HERWIG_MatchboxFactory_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:45 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805229
Default Alt Text
(124 KB)

Event Timeline