Page MenuHomeHEPForge

No OneTemporary

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/MatrixElement/Matchbox/Base/Makefile.am b/MatrixElement/Matchbox/Base/Makefile.am
--- a/MatrixElement/Matchbox/Base/Makefile.am
+++ b/MatrixElement/Matchbox/Base/Makefile.am
@@ -1,18 +1,16 @@
if WANT_DIPOLE
noinst_LTLIBRARIES = libHwMatchboxBase.la
endif
libHwMatchboxBase_la_SOURCES = \
DipoleRepository.cc \
DipoleRepository.h \
MatchboxAmplitude.cc \
MatchboxAmplitude.h \
MatchboxMEBase.cc \
MatchboxMEBase.fh \
MatchboxMEBase.h \
-MatchboxNLOME.cc \
-MatchboxNLOME.h \
MatchboxReweightBase.cc \
MatchboxReweightBase.h \
SubtractedME.cc \
SubtractedME.h
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,1172 +1,1209 @@
// -*- 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 <iterator>
using std::ostream_iterator;
using namespace Herwig;
MatchboxMEBase::MatchboxMEBase()
: MEBase(),
theFactorizationScaleFactor(1.0),
theRenormalizationScaleFactor(1.0),
theVerbose(false),
theFixedCouplings(false), theFixedQEDCouplings(false),
theNLight(0), theGetColourCorrelatedMEs(false),
theCheckPoles(false), theOneLoop(false),
theOneLoopNoBorn(false) {}
MatchboxMEBase::~MatchboxMEBase() {}
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 =
theProcessData->diagramMap()[*p];
if ( res.empty() ) {
res = diagramGenerator()->generate(*p,orderInAlphaS(),orderInAlphaEW());
}
copy(res.begin(),res.end(),back_inserter(diags));
}
if ( diags.empty() )
return;
for ( vector<Ptr<Tree2toNDiagram>::ptr>::iterator d = diags.begin();
d != diags.end(); ++d ) {
add(*d);
}
if ( theVerbose ) {
string fname = name() + ".diagrams";
ifstream test(fname.c_str());
if ( !test ) {
test.close();
ofstream out(fname.c_str());
for ( vector<Ptr<Tree2toNDiagram>::ptr>::const_iterator d = diags.begin();
d != diags.end(); ++d ) {
DiagramDrawer::drawDiag(out,**d);
out << "\n";
}
}
}
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);
if ( phasespace() && !xc->head() )
phasespace()->prepare(xc,theVerbose);
if ( scaleChoice() )
scaleChoice()->setXComb(xc);
if ( cache() )
cache()->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);
int ampShift = 0;
int nBorn = nDimBorn();
if ( matchboxAmplitude() ) {
ampShift = matchboxAmplitude()->nDimAdditional();
if ( ampShift )
matchboxAmplitude()->additionalKinematics(r + nBorn);
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).additionalKinematics(r + nBorn + ampShift);
}
return true;
}
throw Exception()
<< "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return false;
}
int MatchboxMEBase::nDim() const {
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 ( phasespace() ) {
size_t nout = diagrams().front()->partons().size()-2;
int n = phasespace()->nDim(nout);
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()*factorizationScaleFactor();
Energy2 rscale = renormalizationScale()*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 subpro) const {
- for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
- theReweights.begin(); rw != theReweights.end(); ++rw ) {
- (**rw).setXComb(lastXCombPtr());
- if ( !(**rw).apply() )
- continue;
- (**rw).setVetoScales(subpro);
- }
-}
+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) const {
assert(lastXCombPtr()->partonBins().first->pdf());
return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
lastPartons().first->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX1())/lastX1();
}
double MatchboxMEBase::pdf2(Energy2 fscale) const {
assert(lastXCombPtr()->partonBins().second->pdf());
return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
lastPartons().second->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX2())/lastX2();
}
double MatchboxMEBase::me2() const {
if ( matchboxAmplitude() ) {
double res;
if ( !calculateME2(res) )
return res;
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->me2()*
matchboxAmplitude()->lastCrossingSign()*
me2Norm());
cacheME2(lastME2());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::finalStateSymmetry() const {
map<tStdXCombPtr,double>::const_iterator s =
symmetryFactors.find(lastXCombPtr());
if ( s != symmetryFactors.end() )
return s->second;
double symmetryFactor = 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 )
symmetryFactor /= 2.;
else if ( c->second == 3 )
symmetryFactor /= 6.;
else if ( c->second == 4 )
symmetryFactor /= 24.;
}
symmetryFactors[lastXCombPtr()] = symmetryFactor;
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;
}
void MatchboxMEBase::storeColourCorrelatedMEs(double xme2) const {
map<pair<int,int>,double>& ccs =
colourCorrelatedMEs[lastXCombPtr()];
int n = meMomenta().size();
for ( int i = 0; i < n; ++i ) {
if ( !mePartonData()[i]->coloured() )
continue;
for ( int j = i+1; j < n; ++j ) {
if ( !mePartonData()[j]->coloured() )
continue;
if ( noDipole(i,j) )
continue;
ccs[make_pair(i,j)] = colourCorrelatedME2(make_pair(i,j));
}
}
lastXCombPtr()->meta(MatchboxMetaKeys::ColourCorrelatedMEs,ccs);
double& bme2 = bornMEs[lastXCombPtr()];
bme2 = xme2 >= 0. ? xme2 : me2();
lastXCombPtr()->meta(MatchboxMetaKeys::BornME,bme2);
}
CrossSection MatchboxMEBase::dSigHatDR() const {
getPDFWeight();
if ( !lastXComb().willPassCuts() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double xme2 = me2();
lastME2(xme2);
if ( xme2 == 0. && !oneLoopNoBorn() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
if ( getColourCorrelatedMEs() )
storeColourCorrelatedMEs(xme2);
double vme2 = 0.;
if ( oneLoop() )
vme2 = oneLoopInterference();
CrossSection res = ZERO;
if ( !oneLoopNoBorn() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian() * xme2;
if ( oneLoop() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian() * vme2;
if ( !onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).setXComb(lastXCombPtr());
res += (**v).dSigHatDR();
}
}
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 * lastMEPDFWeight());
return lastMECrossSection();
}
double MatchboxMEBase::oneLoopInterference() const {
if ( matchboxAmplitude() ) {
double res;
if ( !calculateME2(res,make_pair<int,int>(-1,-1)) )
return res;
if ( matchboxAmplitude()->oneLoopAmplitudes() )
matchboxAmplitude()->prepareOneLoopAmplitudes(this);
lastME2(matchboxAmplitude()->oneLoopInterference()*
matchboxAmplitude()->lastCrossingSign()*
me2Norm(1));
cacheME2(lastME2(),make_pair<int,int>(-1,-1));
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
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;
}
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()*
matchboxAmplitude()->lastCrossingSign()*
me2Norm(1);
}
return 0.;
}
double MatchboxMEBase::oneLoopSinglePole() const {
if ( matchboxAmplitude() ) {
return
matchboxAmplitude()->oneLoopSinglePole()*
matchboxAmplitude()->lastCrossingSign()*
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() ) {
double res;
if ( !calculateME2(res,ij) )
return res;
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->colourCorrelatedME2(ij)*
matchboxAmplitude()->lastCrossingSign()*
me2Norm());
cacheME2(lastME2(),ij);
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::colourCorrelatedME2() 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() ) {
double res;
if ( !calculateME2(res,ij) )
return res;
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
matchboxAmplitude()->lastCrossingSign()*
me2Norm());
cacheME2(lastME2(),ij);
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 ( cache() )
cache()->flush();
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(meInfo().begin(),meInfo().end(),ostream_iterator<double>(os," "));
+ 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 ( !theVerbose )
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 ( !theVerbose )
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 ( !theVerbose )
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 ( !theVerbose )
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 ( !theVerbose )
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());
matchboxAmplitude()->processData(processData());
}
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;
(**v).setBorn(this);
}
}
void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
os << theReweights << thePhasespace << theAmplitude << theScaleChoice
<< theDiagramGenerator << theProcessData << theSubprocesses
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theVerbose << theCache << theVirtuals << theFixedCouplings << theFixedQEDCouplings
<< theNLight << theGetColourCorrelatedMEs << theCheckPoles
<< theOneLoop << theOneLoopNoBorn << symmetryFactors;
}
void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
is >> theReweights >> thePhasespace >> theAmplitude >> theScaleChoice
>> theDiagramGenerator >> theProcessData >> theSubprocesses
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theVerbose >> theCache >> theVirtuals >> theFixedCouplings >> theFixedQEDCouplings
>> theNLight >> theGetColourCorrelatedMEs >> theCheckPoles
>> theOneLoop >> theOneLoopNoBorn >> symmetryFactors;
}
void MatchboxMEBase::Init() {
static ClassDocumentation<MatchboxMEBase> documentation
("MatchboxMEBase is the base class for matrix elements "
"in the context of the matchbox NLO interface.");
static RefVector<MatchboxMEBase,MatchboxReweightBase> interfaceReweights
("Reweights",
"Reweight objects to be applied to this matrix element.",
&MatchboxMEBase::theReweights, -1, false, false, true, true, false);
static Reference<MatchboxMEBase,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"Set the phasespace generator to be used.",
&MatchboxMEBase::thePhasespace, false, false, true, true, false);
static Reference<MatchboxMEBase,Tree2toNGenerator> interfaceDiagramGenerator
("DiagramGenerator",
"Set the diagram generator to be used.",
&MatchboxMEBase::theDiagramGenerator, false, false, true, true, false);
static Reference<MatchboxMEBase,ProcessData> interfaceProcessData
("ProcessData",
"Set the process data object to be used.",
&MatchboxMEBase::theProcessData, false, false, true, true, false);
static Reference<MatchboxMEBase,MatchboxScaleChoice> interfaceScaleChoice
("ScaleChoice",
"Set the scale choice to be used.",
&MatchboxMEBase::theScaleChoice, false, false, true, true, false);
static Reference<MatchboxMEBase,MatchboxMECache> interfaceMECache
("Cache",
"Set the cache object to be used.",
&MatchboxMEBase::theCache, false, false, true, true, false);
static RefVector<MatchboxMEBase,MatchboxInsertionOperator> interfaceVirtuals
("Virtuals",
"The virtual corrections to be added.",
&MatchboxMEBase::theVirtuals, -1, false, false, true, true, false);
static Parameter<MatchboxMEBase,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&MatchboxMEBase::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxMEBase,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&MatchboxMEBase::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxMEBase,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&MatchboxMEBase::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"On",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"Off",
false);
static Switch<MatchboxMEBase,bool> interfaceFixedCouplings
("FixedCouplings",
"Indicate that no running couplings should be used.",
&MatchboxMEBase::theFixedCouplings, false, false, false);
static SwitchOption interfaceFixedCouplingsOn
(interfaceFixedCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedCouplingsOff
(interfaceFixedCouplings,
"Off",
"Off",
false);
static Switch<MatchboxMEBase,bool> interfaceFixedQEDCouplings
("FixedQEDCouplings",
"Indicate that no running QED couplings should be used.",
&MatchboxMEBase::theFixedQEDCouplings, false, false, false);
static SwitchOption interfaceFixedQEDCouplingsOn
(interfaceFixedQEDCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedQEDCouplingsOff
(interfaceFixedQEDCouplings,
"Off",
"Off",
false);
}
IBPtr MatchboxMEBase::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxMEBase::fullclone() const {
return new_ptr(*this);
}
void MatchboxMEBase::doinit() {
MEBase::doinit();
theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->colourBasis() ) {
size_t dim =
matchboxAmplitude()->colourBasis()->prepare(diagrams(),matchboxAmplitude()->noCorrelations());
matchboxAmplitude()->colourBasisDim(dim);
}
matchboxAmplitude()->nLight(nLight());
}
}
// *** 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,944 +1,949 @@
// -*- 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/MatchboxMECache.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"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Keys for XComb meta information
*/
struct MatchboxMetaKeys {
enum Keys {
BornME,
ColourCorrelatedMEs
};
};
/**
* \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:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxMEBase();
/**
* The destructor.
*/
virtual ~MatchboxMEBase();
//@}
public:
/** @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 theDiagramGenerator; }
/**
* Set the diagram generator.
*/
void diagramGenerator(Ptr<Tree2toNGenerator>::ptr gen) { theDiagramGenerator = gen; }
/**
* Return the process data.
*/
Ptr<ProcessData>::tptr processData() const { return theProcessData; }
/**
* Set the process data.
*/
void processData(Ptr<ProcessData>::ptr pd) { theProcessData = pd; }
/**
* 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 nLight() const { return theNLight; }
/**
* Set the number of light flavours, this matrix
* element is calculated for.
*/
void nLight(unsigned int n) { theNLight = n; }
//@}
/** @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 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; }
-
- /**
* 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 theFactorizationScaleFactor; }
/**
* Set the factorization scale factor
*/
void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }
/**
* 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 theRenormalizationScaleFactor; }
/**
* Set the renormalization scale factor
*/
void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }
/**
* 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 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 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) const;
/**
* Supply the PDF weight for the second incoming parton.
*/
double pdf2(Energy2 factorizationScale = ZERO) 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 { return theCheckPoles; }
/**
* Switch on checking of epsilon pole cancellation.
*/
void doCheckPoles() { theCheckPoles = true; }
/**
* Perform the check of epsilon pole cancellation. Results will be
* written to the log file, one per phasespace point.
*/
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 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;
/**
* Return true, if colour correlated matrix elements should be calculated
* for later use
*/
bool getColourCorrelatedMEs() const { return theGetColourCorrelatedMEs; }
/**
* Switch on calculation of colour correlated matrix elements for
* later use
*/
void doColourCorrelatedMEs() { theGetColourCorrelatedMEs = true; }
/**
* Calculate colour correlated matrix elements for later use
*/
void storeColourCorrelatedMEs(double xme2 = -1.) const;
//@}
/** @name Caching and diagnostic information */
//@{
/**
* Set the ME cache object
*/
void cache(Ptr<MatchboxMECache>::ptr c) { theCache = c; }
/**
* Get the ME cache object
*/
Ptr<MatchboxMECache>::tptr cache() const { return theCache; }
/**
* 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 the matrix element needs to be
* recalculated for the given phase space point.
* If not, return the cached value in the given reference.
*/
bool calculateME2(double& xme2,
const pair<int,int>& corr = make_pair(0,0)) const {
if ( !cache() ) {
xme2 = 0.;
return true;
}
cache()->setXComb(lastXCombPtr());
return cache()->calculateME2(xme2,corr);
}
/**
* Cache a calculated matrix element
* for the last phase space point.
*/
void cacheME2(double xme2,
const pair<int,int>& corr = make_pair(0,0)) const {
if ( !cache() )
return;
cache()->cacheME2(xme2,corr);
}
/**
* Return true, if verbose
*/
bool verbose() const { return theVerbose; }
/**
* Switch on diagnostic information.
*/
void setVerbose(bool on = true) { theVerbose = on; }
/**
+ * 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 = "");
//@}
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();
//@}
protected:
/**
* The final state symmetry factors.
*/
mutable map<tStdXCombPtr,double> symmetryFactors;
private:
/**
* A vector of reweight objects the sum of which
* should be applied to reweight this matrix element
*/
vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
/**
* The phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The amplitude to be used
*/
Ptr<MatchboxAmplitude>::ptr theAmplitude;
/**
* The diagram generator to be used.
*/
Ptr<Tree2toNGenerator>::ptr theDiagramGenerator;
/**
* The process data object to be used
*/
Ptr<ProcessData>::ptr theProcessData;
/**
* The scale choice object
*/
Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
/**
* The ME cache object
*/
Ptr<MatchboxMECache>::ptr theCache;
/**
* The virtual corrections.
*/
vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
/**
* The subprocesses to be considered, if a diagram generator is
* present.
*/
vector<PDVector> theSubprocesses;
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
/**
* Wether or not diagnostic information
* should be written to the generator log
*/
bool theVerbose;
/**
* Use non-running couplings.
*/
bool theFixedCouplings;
/**
* Use non-running couplings.
*/
bool theFixedQEDCouplings;
/**
* The number of light flavours, this matrix
* element is calculated for.
*/
unsigned int theNLight;
/**
* True, if colour correlated matrix elements should be calculated
* for later use
*/
bool theGetColourCorrelatedMEs;
/**
* True, if cancellationn of epsilon poles should be checked.
*/
bool theCheckPoles;
/**
* 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;
/**
* Map xcomb's to storage of Born matrix elements squared.
*/
mutable map<tStdXCombPtr,double> bornMEs;
/**
* Map xcomb's to storage of colour correlated matrix elements.
*/
mutable map<tStdXCombPtr,map<pair<int,int>,double> > colourCorrelatedMEs;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxMEBase & operator=(const MatchboxMEBase &);
};
}
#endif /* HERWIG_MatchboxMEBase_H */
diff --git a/MatrixElement/Matchbox/Base/MatchboxNLOME.cc b/MatrixElement/Matchbox/Base/MatchboxNLOME.cc
deleted file mode 100644
--- a/MatrixElement/Matchbox/Base/MatchboxNLOME.cc
+++ /dev/null
@@ -1,270 +0,0 @@
-// -*- C++ -*-
-//
-// MatchboxNLOME.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.
-//
-//
-// This is the implementation of the non-inlined, non-templated member
-// functions of the MatchboxNLOME class.
-//
-
-#include "MatchboxNLOME.h"
-#include "ThePEG/Interface/ClassDocumentation.h"
-#include "ThePEG/Interface/Reference.h"
-#include "ThePEG/Interface/RefVector.h"
-#include "ThePEG/Utilities/DescribeClass.h"
-#include "ThePEG/Repository/Repository.h"
-#include "ThePEG/Persistency/PersistentOStream.h"
-#include "ThePEG/Persistency/PersistentIStream.h"
-#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
-
-#include <iterator>
-using std::ostream_iterator;
-
-using namespace Herwig;
-
-MatchboxNLOME::MatchboxNLOME()
- : MEBase(), theNDim(-1), theCheckPoles(false) {}
-
-MatchboxNLOME::~MatchboxNLOME() {}
-
-
-IBPtr MatchboxNLOME::clone() const {
- return new_ptr(*this);
-}
-
-IBPtr MatchboxNLOME::fullclone() const {
- return new_ptr(*this);
-}
-
-void MatchboxNLOME::cloneDependencies(const std::string& prefix) {
-
- Ptr<MatchboxMEBase>::ptr myMatrixElement = theMatrixElement->cloneMe();
- ostringstream pname;
- pname << (prefix == "" ? fullName() : prefix) << "/" << myMatrixElement->name();
- if ( ! (generator()->preinitRegister(myMatrixElement,pname.str()) ) )
- throw InitException() << "Matrix element " << pname.str() << " already existing.";
- myMatrixElement->cloneDependencies(pname.str());
- theMatrixElement = myMatrixElement;
-
- 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;
- (**v).setBorn(theMatrixElement);
- }
-
-}
-
-bool MatchboxNLOME::generateKinematics(const double * r) {
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
- virtuals().begin(); v != virtuals().end(); ++v )
- if ( (**v).nDimAdditional() )
- (**v).additionalKinematics(r+theMatrixElement->nDim());
- bool ret = matrixElement()->generateKinematics(r);
- jacobian(matrixElement()->lastXComb().jacobian());
- return ret;
-}
-
-void MatchboxNLOME::logPoles() const {
- double res2me = matrixElement()->oneLoopDoublePole();
- double res1me = matrixElement()->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;
-}
-
-double MatchboxNLOME::me2() const {
- if ( !matrixElement()->onlyOneLoop() && checkPoles() )
- logPoles();
- double res = !matrixElement()->onlyOneLoop() ? matrixElement()->me2() : 0.;
- if ( matrixElement()->haveOneLoop() ) {
- res +=
- matrixElement()->oneLoopInterference();
- }
- if ( !matrixElement()->onlyOneLoop() )
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
- virtuals().begin(); v != virtuals().end(); ++v )
- res += (**v).me2();
- return res;
-}
-
-CrossSection MatchboxNLOME::dSigHatDR() const {
- CrossSection res =
- !matrixElement()->onlyOneLoop() ? matrixElement()->dSigHatDR() : ZERO;
- if ( res == ZERO && !matrixElement()->onlyOneLoop() )
- return res;
- if ( !matrixElement()->onlyOneLoop() && checkPoles() )
- logPoles();
- if ( matrixElement()->haveOneLoop() ) {
- if ( matrixElement()->onlyOneLoop() )
- matrixElement()->getPDFWeight();
- double vme =
- matrixElement()->oneLoopInterference();
- res +=
- sqr(hbarc) * vme *
- jacobian() * lastMEPDFWeight() /
- (2.*lastSHat());
- }
- if ( !matrixElement()->onlyOneLoop() )
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
- virtuals().begin(); v != virtuals().end(); ++v )
- res += (**v).dSigHatDR();
- lastMECrossSection(res);
- return res;
-}
-
-void MatchboxNLOME::flushCaches() {
- theMatrixElement->flushCaches();
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
- virtuals().begin(); v != virtuals().end(); ++v )
- (**v).flushCaches();
-}
-
-void MatchboxNLOME::setXComb(tStdXCombPtr xc) {
-
- MEBase::setXComb(xc);
- theMatrixElement->setXComb(xc);
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
- virtuals().begin(); v != virtuals().end(); ++v )
- (**v).setXComb(xc);
-}
-
-void MatchboxNLOME::getNDim() const {
- int maxadd = 0;
- for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
- virtuals().begin(); v != virtuals().end(); ++v ) {
- if ( (**v).nDimAdditional() > 1 ) {
- throw InitException() << "at most one additional random number supported for "
- << "virtual corrections currently";
- }
- maxadd = max(maxadd,(**v).nDimAdditional());
- }
- theNDim = theMatrixElement->nDim() + maxadd;
-}
-
-void MatchboxNLOME::print(ostream& os) const {
-
- os << "--- MatchboxNLOME setup --------------------------------------------------------\n";
-
- os << " '" << name() << "' using\n"
- << " matrix element '" << matrixElement()->name() << "' and 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 MatchboxNLOME::printLastEvent(ostream& os) const {
-
- os << "--- MatchboxNLOME 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(meInfo().begin(),meInfo().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 MatchboxNLOME::doinit() {
- MEBase::doinit();
-}
-
-void MatchboxNLOME::doinitrun() {
- MEBase::doinitrun();
-}
-
-void MatchboxNLOME::persistentOutput(PersistentOStream & os) const {
- os << theMatrixElement << theVirtuals << theNDim << theCheckPoles;
-}
-
-void MatchboxNLOME::persistentInput(PersistentIStream & is, int) {
- is >> theMatrixElement >> theVirtuals >> theNDim >> theCheckPoles;
-}
-
-void MatchboxNLOME::Init() {
-
- static ClassDocumentation<MatchboxNLOME> documentation
- ("MatchboxNLOME");
-
-
- static Reference<MatchboxNLOME,MatchboxMEBase> interfaceBornME
- ("BornME",
- "The Born matrix element",
- &MatchboxNLOME::theMatrixElement, false, false, true, false, false);
-
-
- static RefVector<MatchboxNLOME,MatchboxInsertionOperator> interfaceVirtuals
- ("Virtuals",
- "The virtual corrections to be added.",
- &MatchboxNLOME::theVirtuals, -1, false, false, true, true, 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<MatchboxNLOME,MEBase>
-describeMatchboxNLOME("Herwig::MatchboxNLOME", "HwMatchbox.so");
-
diff --git a/MatrixElement/Matchbox/Base/MatchboxNLOME.h b/MatrixElement/Matchbox/Base/MatchboxNLOME.h
deleted file mode 100644
--- a/MatrixElement/Matchbox/Base/MatchboxNLOME.h
+++ /dev/null
@@ -1,444 +0,0 @@
-// -*- C++ -*-
-//
-// MatchboxNLOME.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_MatchboxNLOME_H
-#define HERWIG_MatchboxNLOME_H
-//
-// This is the declaration of the MatchboxNLOME class.
-//
-
-#include "ThePEG/MatrixElement/MEBase.h"
-#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
-#include "Herwig++/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h"
-
-namespace Herwig {
-
-using namespace ThePEG;
-
-/**
- * \ingroup Matchbox
- * \author Simon Platzer
- *
- * \brief MatchboxNLOME assembles a MatchboxMEBase
- * and several MatchboxInsertionOperator objects to calculate
- * the Born + finite virtual cross section in a NLO
- * QCD calculation.
- *
- * @see \ref MatchboxNLOMEInterfaces "The interfaces"
- * defined for MatchboxNLOME.
- */
-class MatchboxNLOME: public MEBase {
-
-public:
-
- /** @name Standard constructors and destructors. */
- //@{
- /**
- * The default constructor.
- */
- MatchboxNLOME();
-
- /**
- * The destructor.
- */
- virtual ~MatchboxNLOME();
- //@}
-
-public:
-
- /** @name Subprocess and diagram information. */
- //@{
-
- /**
- * Add all possible diagrams with the add() function.
- */
- virtual void getDiagrams() const {
- theMatrixElement->diagrams();
- useDiagrams(theMatrixElement);
- }
-
- /**
- * 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; }
-
- /**
- * 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 & dv) const {
- return matrixElement()->diagrams(dv);
- }
-
- /**
- * Select a diagram. Default version uses diagrams(const
- * DiagramVector &) to select a diagram according to the
- * weights. This is the only method used that should be outside of
- * MEBase.
- */
- virtual DiagramIndex diagram(const DiagramVector & dv) const {
- return matrixElement()->diagram(dv);
- }
-
- /**
- * 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 matrixElement()->colourGeometries(diag); }
-
- /**
- * Select a ColpurLines geometry. The default version returns a
- * colour geometry selected among the ones returned from
- * colourGeometries(tcDiagPtr).
- */
- virtual const ColourLines &
- selectColourGeometry(tcDiagPtr diag) const { return matrixElement()->selectColourGeometry(diag); }
-
- /**
- * Return the order in \f$\alpha_S\f$ in which this matrix element
- * is given.
- */
- virtual unsigned int orderInAlphaS() const {
- return
- virtuals().empty() ? theMatrixElement->orderInAlphaS() :
- theMatrixElement->orderInAlphaS() + 1;
- }
-
- /**
- * Return the order in \f$\alpha_{EM}\f$ in which this matrix
- * element is given. Returns 0.
- */
- virtual unsigned int orderInAlphaEW() const { return theMatrixElement->orderInAlphaEW(); }
-
- //@}
-
- /** @name Phasespace generation */
- //@{
-
- /**
- * Set the XComb object to be used in the next call to
- * generateKinematics() and dSigHatDR().
- */
- virtual void setXComb(tStdXCombPtr);
-
- /**
- * Set the typed and momenta of the incoming and outgoing partons to
- * be used in subsequent calls to me() and colourGeometries()
- * according to the associated XComb object. If the function is
- * overridden in a sub class the new function must call the base
- * class one first.
- */
- virtual void setKinematics() { matrixElement()->setKinematics(); }
-
- /**
- * 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.
- */
- virtual bool generateKinematics(const double * r);
-
- /**
- * The number of internal degreed of freedom used in the matrix
- * element. This default version returns 0;
- */
- virtual int nDim() const { if ( theNDim < 0 ) getNDim(); return theNDim; }
-
- /**
- * Get the ranom numbers needed
- */
- void getNDim() 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 theMatrixElement->haveX1X2(); }
-
- /**
- * Return true, if this matrix element expects
- * the incoming partons in their center-of-mass system
- */
- virtual bool wantCMS () const { return theMatrixElement->wantCMS(); }
-
- /**
- * 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; }
-
- /**
- * Clear the information previously provided by a call to
- * setKinematics(...).
- */
- virtual void clearKinematics() { matrixElement()->clearKinematics(); }
-
- //@}
-
- /** @name Scale choices, couplings and PDFs */
- //@{
-
- /**
- * Return the scale associated with the phase space point provided
- * by the last call to setKinematics().
- */
- virtual Energy2 scale() const { return matrixElement()->scale(); }
-
- /**
- * 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 matrixElement()->alphaS(); }
-
- /**
- * 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 matrixElement()->alphaEM(); }
-
- /**
- * Return true, if this matrix element provides the PDF
- * weight for the first incoming parton itself.
- */
- virtual bool havePDFWeight1() const { return theMatrixElement->havePDFWeight1(); }
-
- /**
- * Return true, if this matrix element provides the PDF
- * weight for the second incoming parton itself.
- */
- virtual bool havePDFWeight2() const { return theMatrixElement->havePDFWeight2(); }
-
- //@}
-
- /** @name Matrix elements and evaluation */
- //@{
-
- /**
- * Return the Born matrix element
- */
- Ptr<MatchboxMEBase>::tcptr matrixElement() const { return theMatrixElement; }
-
- /**
- * Return the Born matrix element
- */
- Ptr<MatchboxMEBase>::tptr matrixElement() { return theMatrixElement; }
-
- /**
- * Set the Born matrix element
- */
- void matrixElement(Ptr<MatchboxMEBase>::ptr b) { theMatrixElement = b; }
-
- /**
- * Return the virtual corrections
- */
- const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const {
- return theVirtuals;
- }
-
- /**
- * Return the virtual corrections
- */
- vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() {
- return theVirtuals;
- }
-
- /**
- * 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; }
-
- /**
- * Perform the check of epsilon pole cancellation. Results will be
- * written to the log file, one per phasespace point.
- */
- void logPoles() const;
-
- /**
- * 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 matrix element squared differential in the variables
- * given by the last call to generateKinematics().
- */
- virtual CrossSection dSigHatDR() 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();
-
- /**
- * Dump the setup to an ostream
- */
- void print(ostream&) const;
-
- /**
- * Print information on last event generated.
- */
- void printLastEvent(ostream& os) const;
-
- //@}
-
- /** @name Methods used to setup NLO ME objects */
- //@{
-
- /**
- * Clone the dependencies, using a given prefix.
- */
- void cloneDependencies(const std::string& prefix = "");
-
- //@}
-
- /** @name Methods required to setup the event record */
- //@{
-
- /**
- * construct the spin information for the interaction
- */
- virtual void constructVertex(tSubProPtr sub) { matrixElement()->constructVertex(sub); }
-
- /**
- * Comlete a SubProcess object using the internal degrees of freedom
- * generated in the last generateKinematics() (and possible other
- * degrees of freedom which was intergated over in dSigHatDR(). This
- * default version does nothing. Will be made purely virtual in the
- * future.
- */
- virtual void generateSubCollision(SubProcess & sub) { matrixElement()->generateSubCollision(sub); }
-
- //@}
-
-public:
-
- /** @name Functions used by the persistent I/O system. */
- //@{
- /**
- * Function used to write out object persistently.
- * @param os the persistent output stream written to.
- */
- void persistentOutput(PersistentOStream & os) const;
-
- /**
- * Function used to read in object persistently.
- * @param is the persistent input stream read from.
- * @param version the version number of the object when written.
- */
- void persistentInput(PersistentIStream & is, int version);
- //@}
-
- /**
- * The standard Init function used to initialize the interfaces.
- * Called exactly once for each class by the class description system
- * before the main function starts or
- * when this class is dynamically loaded.
- */
- static void Init();
-
-protected:
-
- /** @name Clone Methods. */
- //@{
- /**
- * Make a simple clone of this object.
- * @return a pointer to the new object.
- */
- virtual IBPtr clone() const;
-
- /** Make a clone of this object, possibly modifying the cloned object
- * to make it sane.
- * @return a pointer to the new object.
- */
- virtual IBPtr fullclone() const;
- //@}
-
-
-protected:
-
- /** @name Standard Interfaced functions. */
- //@{
- /**
- * Initialize this object after the setup phase before saving an
- * EventGenerator to disk.
- * @throws InitException if object could not be initialized properly.
- */
- virtual void doinit();
-
- /**
- * Initialize this object. Called in the run phase just before
- * a run begins.
- */
- virtual void doinitrun();
- //@}
-
-private:
-
- /**
- * The Born matrix element
- */
- Ptr<MatchboxMEBase>::ptr theMatrixElement;
-
- /**
- * The virtual corrections.
- */
- vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
-
- /**
- * The number of random variables needed
- */
- mutable int theNDim;
-
- /**
- * True, if cancellationn of epsilon poles should be checked.
- */
- bool theCheckPoles;
-
-private:
-
- /**
- * The assignment operator is private and must never be called.
- * In fact, it should not even be implemented.
- */
- MatchboxNLOME & operator=(const MatchboxNLOME &);
-
-};
-
-}
-
-#endif /* HERWIG_MatchboxNLOME_H */
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.cc b/MatrixElement/Matchbox/Base/SubtractedME.cc
--- a/MatrixElement/Matchbox/Base/SubtractedME.cc
+++ b/MatrixElement/Matchbox/Base/SubtractedME.cc
@@ -1,704 +1,669 @@
// -*- C++ -*-
//
// SubtractedME.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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SubtractedME class.
//
#include "SubtractedME.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
using namespace Herwig;
SubtractedME::SubtractedME()
: MEGroup(),
theSubtractionData(""),
theVerbose(false),
theSubProcessGroups(false), theInclusive(false),
theVetoScales(false),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false) {}
SubtractedME::~SubtractedME() {}
IBPtr SubtractedME::clone() const {
return new_ptr(*this);
}
IBPtr SubtractedME::fullclone() const {
return new_ptr(*this);
}
void SubtractedME::setXComb(tStdXCombPtr xc) {
MEGroup::setXComb(xc);
}
MEBase::DiagramVector SubtractedME::dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const {
Ptr<SubtractionDipole>::tptr dipole =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(depME);
if ( !dipole ) {
Throw<InitException>() << "A dependent matrix element of SubtractedME "
<< "has not been derived from SubtractionDipole. "
<< "Please check the corresponding input file.";
}
return dipole->underlyingBornDiagrams(proc);
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::dipoles() {
if ( allDipoles().empty() ) {
allDipoles() = DipoleRepository::dipoles();
}
if ( dependent().empty() )
getDipoles();
vector<Ptr<SubtractionDipole>::ptr> res;
for ( MEVector::const_iterator k = dependent().begin();
k != dependent().end(); ++k )
res.push_back(dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(*k));
return res;
}
void SubtractedME::getDipoles() {
if ( !dependent().empty() )
return;
if ( allDipoles().empty() ) {
allDipoles() = DipoleRepository::dipoles();
}
Ptr<MatchboxMEBase>::tptr real =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( allDipoles().empty() || theBorns.empty() || !real )
Throw<InitException>() << "The SubtractedME '"
<< name() << "' could not generate "
<< "subtraction terms for the real emission "
<< "matrix element '" << real->name() << "'. "
<< "Please check the corresponding input file.";
Ptr<MatchboxMEBase>::ptr myRealEmissionME = real->cloneMe();
ostringstream pname;
pname << fullName() << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myRealEmissionME->cloneDependencies(pname.str());
head(myRealEmissionME);
real = myRealEmissionME;
MEVector dipMEs;
vector<Ptr<SubtractionDipole>::ptr> genDipoles
= real->getDipoles(allDipoles(),theBorns);
if ( theSubtractionData != "" ) {
theSubProcessGroups = true;
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
genDipoles.begin(); d != genDipoles.end(); ++d )
(**d).doTestSubtraction();
}
if ( genDipoles.empty() ) {
// probably finite real contribution, but warn
generator()->log() << "\nWarning: No subtraction dipoles could be found for the processes:\n";
for ( vector<PDVector>::const_iterator s = real->subProcesses().begin();
s != real->subProcesses().end(); ++s ) {
generator()->log() << (*s)[0]->PDGName() << " " << (*s)[1]->PDGName()
<< " -> ";
for ( PDVector::const_iterator p = s->begin() + 2; p != s->end(); ++p )
generator()->log() << (**p).PDGName() << " ";
generator()->log() << "\n" << flush;
}
generator()->log() << "Assuming finite tree-level O(alphaS) correction.\n";
}
dipMEs.resize(genDipoles.size());
copy(genDipoles.begin(),genDipoles.end(),dipMEs.begin());
dependent() = dipMEs;
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::splitDipoles(const cPDVector& born) {
vector<Ptr<SubtractionDipole>::ptr> dips = dipoles();
vector<Ptr<SubtractionDipole>::ptr> res;
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = dips.begin();
d != dips.end(); ++d ) {
for ( DiagramVector::const_iterator p = (**d).underlyingBornME()->diagrams().begin();
p != (**d).underlyingBornME()->diagrams().end(); ++p )
if ( born == (**p).partons() ) {
res.push_back(*d);
break;
}
}
return res;
}
void SubtractedME::showerApproximation(Ptr<ShowerApproximation>::tptr app) {
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->showerApproximation(app);
}
}
void SubtractedME::doRealShowerSubtraction() {
theRealShowerSubtraction = true;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->doRealShowerSubtraction();
}
}
void SubtractedME::doVirtualShowerSubtraction() {
theVirtualShowerSubtraction = true;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->doVirtualShowerSubtraction();
}
}
-void SubtractedME::setVetoScales(tSubProPtr subpro) const {
-
- if ( !vetoScales() )
- return;
-
- Ptr<SubtractionDipole>::tptr dipole;
- Energy pt;
-
- assert(head()->noMirror());
-
- for ( MEVector::const_iterator d = dependent().begin();
- d != dependent().end(); ++d ) {
-
- dipole = dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(*d);
- assert(dipole);
- pt = dipole->lastPt();
-
- if ( dipole->realEmitter() == 0 ||
- dipole->realSpectator() == 0 ) {
- if ( subpro->incoming().first->vetoScale() < 0.0*GeV2 ||
- subpro->incoming().first->vetoScale() > sqr(pt) )
- subpro->incoming().first->vetoScale(sqr(pt));
- }
-
- if ( dipole->realEmitter() == 1 ||
- dipole->realSpectator() == 1 ) {
- if ( subpro->incoming().second->vetoScale() < 0.0*GeV2 ||
- subpro->incoming().second->vetoScale() > sqr(pt) )
- subpro->incoming().second->vetoScale(sqr(pt));
- }
-
- if ( dipole->realEmitter() > 1 ) {
- if ( subpro->outgoing()[dipole->realEmitter()-2]->vetoScale() < 0.0*GeV2 ||
- subpro->outgoing()[dipole->realEmitter()-2]->vetoScale() > sqr(pt) )
- subpro->outgoing()[dipole->realEmitter()-2]->vetoScale(sqr(pt));
- }
-
- if ( dipole->realSpectator() > 1 ) {
- if ( subpro->outgoing()[dipole->realSpectator()-2]->vetoScale() < 0.0*GeV2 ||
- subpro->outgoing()[dipole->realSpectator()-2]->vetoScale() > sqr(pt) )
- subpro->outgoing()[dipole->realSpectator()-2]->vetoScale(sqr(pt));
- }
-
- if ( subpro->outgoing()[dipole->realEmission()-2]->vetoScale() < 0.0*GeV2 ||
- subpro->outgoing()[dipole->realEmission()-2]->vetoScale() > sqr(pt) )
- subpro->outgoing()[dipole->realEmission()-2]->vetoScale(sqr(pt));
-
- }
-
-}
+void SubtractedME::setVetoScales(tSubProPtr) const {}
void SubtractedME::fillProjectors() {
if ( !inclusive() && !virtualShowerSubtraction() )
return;
Ptr<StdXCombGroup>::tptr group =
dynamic_ptr_cast<Ptr<StdXCombGroup>::tptr>(lastXCombPtr());
for ( vector<StdXCombPtr>::const_iterator d = group->dependent().begin();
d != group->dependent().end(); ++d ) {
if ( (**d).lastCrossSection() != ZERO )
lastXCombPtr()->projectors().insert(1.,*d);
}
}
-double SubtractedME::reweightHead() {
- double sum = 0.;
+double SubtractedME::reweightHead(const vector<tStdXCombPtr>& dep) {
if ( realShowerSubtraction() ) {
assert(showerApproximation());
bool below = showerApproximation()->belowCutoff();
- showerApproximation()->resetBelowCutoff();
if ( below )
return 0.;
- return lastXComb().projectors().size();
+ return 1.;
}
if ( virtualShowerSubtraction() ) {
assert(showerApproximation());
bool above = !showerApproximation()->belowCutoff();
- showerApproximation()->resetBelowCutoff();
if ( above )
return 0.;
}
- for ( Selector<tStdXCombPtr>::const_iterator d =
- lastXComb().projectors().begin(); d != lastXComb().projectors().end(); ++d )
- sum += d->second->lastME2();
+ double sum = 0.;
+ for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d )
+ sum += (**d).lastME2();
return
- lastXComb().projectors().size() * lastXComb().lastProjector()->lastME2() / sum;
+ dep.size() * lastXComb().lastProjector()->lastME2() / sum;
}
-double SubtractedME::reweightDependent(tStdXCombPtr xc) {
+double SubtractedME::reweightDependent(tStdXCombPtr xc, const vector<tStdXCombPtr>& dep) {
+ if ( realShowerSubtraction() ) {
+ assert(showerApproximation());
+ double sum = 0.;
+ for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d )
+ sum += (**d).lastME2();
+ return xc->lastME2()/sum;
+ }
if ( xc != lastXComb().lastProjector() )
return 0.;
- return lastXComb().projectors().size();
+ double w = 1.;
+ if ( virtualShowerSubtraction() ) {
+ assert(showerApproximation());
+ double sum = 0.;
+ for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d )
+ sum += (**d).lastME2();
+ assert(xc->meInfo().size() == 1);
+ double r = xc->meInfo()[0];
+ w -= xc->lastME2()*r/sum;
+ }
+ return w * dep.size();
}
void SubtractedME::doinit() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinit();
return;
}
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theReal ) {
getDipoles();
}
if ( theVerbose )
print(Repository::clog());
MEGroup::doinit();
}
void SubtractedME::doinitrun() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinitrun();
return;
}
MEGroup::doinitrun();
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theSubtractionData != "" ) {
set<cPDVector> procs;
for ( DiagramVector::const_iterator d = head()->diagrams().begin();
d != head()->diagrams().end(); ++d ) {
if ( procs.find((**d).partons()) == procs.end() )
procs.insert((**d).partons());
}
for ( set<cPDVector>::const_iterator p = procs.begin();
p != procs.end(); ++p ) {
for ( size_t i = 0; i < (*p).size(); ++i ) {
if ( !(*p)[i]->coloured() )
continue;
if ( i > 1 &&
(*p)[i]->id() == ParticleID::g ) {
softHistograms[SoftSubtractionIndex(*p,i)] = SubtractionHistogram(0.001,10.);
if ( theReal->phasespace() )
theReal->phasespace()->singularLimit(i,i);
}
for ( size_t j = i+1; j < (*p).size(); ++j ) {
if ( !(*p)[j]->coloured() )
continue;
long iid = (*p)[i]->id();
long jid = (*p)[j]->id();
if ( i < 2 && j < 2 )
continue;
if ( i < 2 && j > 1 ) {
if ( abs(iid) < 7 && abs(jid) < 7 && iid != jid )
continue;
}
if ( i > 1 && j > 1 ) {
if ( abs(iid) < 7 && abs(jid) < 7 && iid + jid != 0 )
continue;
}
bool haveDipole = false;
for ( MEVector::const_iterator k = dependent().begin();
k != dependent().end(); ++k ) {
const SubtractionDipole& dip = dynamic_cast<const SubtractionDipole&>(**k);
if ( ( (size_t)(dip.realEmitter()) == i && (size_t)(dip.realEmission()) == j ) ||
( (size_t)(dip.realEmitter()) == j && (size_t)(dip.realEmission()) == i ) ) {
haveDipole = true;
break;
}
}
if ( !haveDipole )
continue;
collinearHistograms[CollinearSubtractionIndex(*p,make_pair(i,j))] = SubtractionHistogram(0.001,20.);
if ( theReal->phasespace() )
theReal->phasespace()->singularLimit(i,j);
}
}
}
}
}
void SubtractedME::dofinish() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::dofinish();
return;
}
MEGroup::dofinish();
for ( map<CollinearSubtractionIndex,SubtractionHistogram>::
const_iterator b = collinearHistograms.begin();
b != collinearHistograms.end(); ++b ) {
b->second.dump(theSubtractionData,
b->first.first,
b->first.second.first,
b->first.second.second);
}
for ( map<SoftSubtractionIndex,SubtractionHistogram>::
const_iterator b = softHistograms.begin();
b != softHistograms.end(); ++b ) {
b->second.dump(theSubtractionData,
b->first.first,
b->first.second,
b->first.second);
}
}
void SubtractedME::rebind(const TranslationMap & trans) {
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
*d = trans.translate(*d);
MEGroup::rebind(trans);
}
IVector SubtractedME::getReferences() {
IVector ret = MEGroup::getReferences();
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
ret.push_back(*d);
return ret;
}
void SubtractedME::print(ostream& os) const {
os << "--- SubtractedME setup ---------------------------------------------------------\n";
os << " '" << name() << "' subtracting real emission\n '"
<< head()->name() << "' using the dipoles:\n";
for ( MEVector::const_iterator d = dependent().begin();
d != dependent().end(); ++d )
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*d)->print(os);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractedME::printLastEvent(ostream& os) const {
os << "--- SubtractedME last event information ----------------------------------------\n";
os << " for subtracted matrix element '" << name() << "'\n";
os << " real emission event information:\n";
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head())->printLastEvent(os);
os << " dipoles event information:\n";
for ( MEVector::const_iterator d = dependent().begin();
d != dependent().end(); ++d )
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*d)->printLastEvent(os);
os << "--- end SubtractedME last event information ------------------------------------\n\n\n";
os << flush;
}
void SubtractedME::lastEventStatistics() {
MEGroup::lastEventStatistics();
if ( !generator() )
return;
/*
if ( theVerbose )
printLastEvent(generator()->log());
*/
if ( !collinearHistograms.empty() )
lastEventSubtraction();
}
SubtractedME::SubtractionHistogram::
SubtractionHistogram(double low,
double up,
unsigned int nbins)
: lower(low) {
nbins = nbins + 1;
double c = log10(up/low) / (nbins-1.);
for ( unsigned int k = 1; k < nbins; ++k ) {
bins[low*pow(10.0,k*c)] = make_pair(Constants::MaxDouble,0.);
}
}
void SubtractedME::SubtractionHistogram::
dump(const std::string& prefix,
const cPDVector& proc,
int i, int j) const {
ostringstream fname("");
for ( cPDVector::const_iterator p = proc.begin();
p != proc.end(); ++p )
fname << (**p).PDGName();
fname << "-" << i << "-" << j;
ofstream out((prefix+fname.str()+".dat").c_str());
for ( map<double,pair<double,double> >::const_iterator b = bins.begin();
b != bins.end(); ++b ) {
map<double,pair<double,double> >::const_iterator bp = b; --bp;
if ( b->second.first != Constants::MaxDouble &&
b->second.second != 0. ) {
if ( b != bins.begin() )
out << bp->first;
else
out << lower;
out << " " << b->first
<< " " << b->second.first
<< " " << b->second.second
<< "\n" << flush;
}
}
ofstream gpout((prefix+fname.str()+".gp").c_str());
gpout << "set terminal epslatex color solid;\n"
<< "set output '" << fname.str() << "-plot.tex';\n"
<< "set log x;\n"
<< "set size 0.5,0.6;\n"
<< "set yrange [0:2];\n"
<< "set xrange [0.001:10];\n";
if ( i != j ) {
gpout << "set xlabel '$\\sqrt{s_{" << i << j << "}}/{\\rm GeV}$'\n";
} else {
gpout << "set xlabel '$E_{" << i << "}/{\\rm GeV}$'\n";
}
gpout << "plot 1 w lines lc rgbcolor \"#DDDDDD\" notitle, '" << fname.str()
<< ".dat' u (($1+$2)/2.):3:($4 < 4. ? $4 : 4.) w filledcurves lc rgbcolor \"#00AACC\" t "
<< "'$";
for ( size_t k = 0; k < proc.size(); k++ ) {
if ( k == 2 )
gpout << "\\to ";
gpout << (proc[k]->id() < 0 ? "\\bar{" : "")
<< (proc[k]->id() < 0 ? proc[k]->CC()->PDGName() : proc[k]->PDGName())
<< (proc[k]->id() < 0 ? "}" : "") << " ";
}
gpout << "$';\n";
gpout << "reset;\n";
}
void SubtractedME::lastEventSubtraction() {
tStdXCombGroupPtr xc = dynamic_ptr_cast<tStdXCombGroupPtr>(lastXCombPtr());
CrossSection xcme2 = xc->lastHeadCrossSection();
CrossSection xcdip = ZERO;
if ( xcme2 == ZERO )
return;
for ( vector<StdXCombPtr>::const_iterator d = xc->dependent().begin();
d != xc->dependent().end(); ++d ) {
if ( !(*d) )
continue;
if ( !(**d).matrixElement()->apply() )
continue;
xcdip += (**d).lastCrossSection();
}
if ( theReal->phasespace() ) {
size_t i = theReal->phasespace()->lastSingularLimit().first;
size_t j = theReal->phasespace()->lastSingularLimit().second;
if ( i == j &&
softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,abs(xcdip)/abs(xcme2));
}
if ( i != j &&
collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j)))
!= collinearHistograms.end() ) {
double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV;
collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))].
book(s,abs(xcdip)/abs(xcme2));
}
return;
}
for ( size_t i = 0; i < meMomenta().size(); ++i ) {
if ( i > 1 ) {
if ( softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,abs(xcdip)/abs(xcme2));
}
}
for ( size_t j = i+1; j < meMomenta().size(); ++j ) {
if ( collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j)))
== collinearHistograms.end() )
continue;
double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV;
collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))].
book(s,abs(xcdip)/abs(xcme2));
}
}
}
void SubtractedME::persistentOutput(PersistentOStream & os) const {
os << theDipoles << theBorns << theSubtractionData
<< theVerbose << theInclusive << theSubProcessGroups << theVetoScales
<< theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction;
}
void SubtractedME::persistentInput(PersistentIStream & is, int) {
is >> theDipoles >> theBorns >> theSubtractionData
>> theVerbose >> theInclusive >> theSubProcessGroups >> theVetoScales
>> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction;
}
void SubtractedME::Init() {
static ClassDocumentation<SubtractedME> documentation
("SubtractedME represents a subtracted real emission matrix element.");
static RefVector<SubtractedME,MatchboxMEBase> interfaceBorns
("Borns",
"The underlying Born matrix elements to be considered",
&SubtractedME::theBorns, -1, false, false, true, true, false);
static Parameter<SubtractedME,string> interfaceSubtractionData
("SubtractionData",
"File to dump subtraction check to.",
&SubtractedME::theSubtractionData, "",
false, false);
static Switch<SubtractedME,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&SubtractedME::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"On",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"Off",
false);
static Switch<SubtractedME,bool> interfaceSubProcessGroups
("SubProcessGroups",
"Switch on or off production of sub-process groups.",
&SubtractedME::theSubProcessGroups, false, false, false);
static SwitchOption interfaceSubProcessGroupsOn
(interfaceSubProcessGroups,
"On",
"On",
true);
static SwitchOption interfaceSubProcessGroupsOff
(interfaceSubProcessGroups,
"Off",
"Off",
false);
static Switch<SubtractedME,bool> interfaceInclusive
("Inclusive",
"Switch on or off production of inclusive cross section.",
&SubtractedME::theInclusive, false, false, false);
static SwitchOption interfaceInclusiveOn
(interfaceInclusive,
"On",
"On",
true);
static SwitchOption interfaceInclusiveOff
(interfaceInclusive,
"Off",
"Off",
false);
static Switch<SubtractedME,bool> interfaceVetoScales
("VetoScales",
"Switch on or off production of sub-process groups.",
&SubtractedME::theVetoScales, false, false, false);
static SwitchOption interfaceVetoScalesOn
(interfaceVetoScales,
"On",
"On",
true);
static SwitchOption interfaceVetoScalesOff
(interfaceVetoScales,
"Off",
"Off",
false);
static Reference<SubtractedME,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&SubtractedME::theShowerApproximation, false, false, true, true, 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<SubtractedME,MEGroup>
describeHerwigSubtractedME("Herwig::SubtractedME", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.h b/MatrixElement/Matchbox/Base/SubtractedME.h
--- a/MatrixElement/Matchbox/Base/SubtractedME.h
+++ b/MatrixElement/Matchbox/Base/SubtractedME.h
@@ -1,516 +1,534 @@
// -*- C++ -*-
//
// SubtractedME.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_SubtractedME_H
#define HERWIG_SubtractedME_H
//
// This is the declaration of the SubtractedME class.
//
#include "ThePEG/MatrixElement/MEGroup.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief SubtractedME represents a subtracted real emission matrix element.
*
* @see \ref SubtractedMEInterfaces "The interfaces"
* defined for SubtractedME.
*/
class SubtractedME: public MEGroup {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SubtractedME();
/**
* The destructor.
*/
virtual ~SubtractedME();
//@}
public:
/** @name Phasespace and subprocess information */
//@{
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr);
/**
* Return true, if the same additional random numbers
* should be presented to any of the dependent
* matrix elements.
*/
virtual bool uniformAdditional() const { return true; }
/**
+ * 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; }
+
+ /**
* Given a process from the head matrix element,
* return a list of diagrams which should be considered for
* the given dependent matrix element.
*/
virtual MEBase::DiagramVector dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const;
/**
* 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.
*/
virtual bool subProcessGroups() const { return theSubProcessGroups; }
/**
* Switch on or off producing subprocess groups.
*/
void setSubProcessGroups(bool on = true) { theSubProcessGroups = on; }
/**
* Return true, if one of the dependent subprocesses should be
* constructed in place of the one driven by the head matrix element
* or a full subprocess group.
*/
virtual bool selectDependentSubProcess() const { return theInclusive; }
/**
* 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; }
/**
* Fill the projectors object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
virtual void fillProjectors();
/**
* Return true, if this MEGroup will reweight the contributing cross
* sections.
*/
virtual bool groupReweighted() const { return inclusive() || showerApproximation(); }
/**
* Reweight the head cross section
*/
- virtual double reweightHead();
+ virtual double reweightHead(const vector<tStdXCombPtr>&);
/**
* Reweight the dependent cross section
*/
- virtual double reweightDependent(tStdXCombPtr);
+ virtual double reweightDependent(tStdXCombPtr, const vector<tStdXCombPtr>&);
//@}
/** @name Methods relevant to matching */
//@{
/**
+ * Inform this matrix element that a new phase space
+ * point is about to be generated, so all caches should
+ * be flushed.
+ */
+ virtual void flushCaches() {
+ MEGroup::flushCaches();
+ if ( showerApproximation() )
+ showerApproximation()->resetBelowCutoff();
+ }
+
+ /**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr);
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
/**
* Indicate that the shower real emission contribution should be subtracted.
*/
void doRealShowerSubtraction();
/**
* Return true, if the shower real emission contribution should be subtracted.
*/
bool realShowerSubtraction() const { return theRealShowerSubtraction; }
/**
* Indicate that the shower virtual contribution should be subtracted.
*/
void doVirtualShowerSubtraction();
/**
* Return true, if the shower virtual contribution should be subtracted.
*/
bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
//@}
/** @name Matrix element and dipole information */
//@{
/**
* Return the subtraction dipoles.
*/
vector<Ptr<SubtractionDipole>::ptr> dipoles();
/**
* Return the underlying born matrix elements.
*/
const vector<Ptr<MatchboxMEBase>::ptr>& borns() const { return theBorns; }
/**
* Access the underlying born matrix elements.
*/
vector<Ptr<MatchboxMEBase>::ptr>& borns() { return theBorns; }
/**
* Build up dipoles needed.
*/
void getDipoles();
/**
* Return all dipoles matching the given Born process
*/
vector<Ptr<SubtractionDipole>::ptr> splitDipoles(const cPDVector&);
/**
* Return the subtraction dipoles.
*/
const vector<Ptr<SubtractionDipole>::ptr>& allDipoles() const { return theDipoles; }
/**
* Access the subtraction dipoles.
*/
vector<Ptr<SubtractionDipole>::ptr>& allDipoles() { return theDipoles; }
//@}
/** @name Veto scale settings */
//@{
/**
* 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 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 Diagnostic information */
//@{
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Collect information on the last evaluated phasespace
* point for verification or debugging purposes. This
* only called, if the StdXCombGroup did accumulate
* a non-zero cross section from this ME group.
*/
virtual void lastEventStatistics();
/**
* Print debug information on the last event
*/
void printLastEvent(ostream&) const;
/**
* Check the subtraction for the last event
*/
void lastEventSubtraction();
/**
* Set a file to print subtraction check to
*/
void subtractionData(string n) { theSubtractionData = n; }
/**
* Return true, if verbose
*/
bool verbose() const { return theVerbose; }
/**
* Switch on or off verbosity for this subtracted ME
*/
void setVerbose(bool on = true) { theVerbose = on; }
//@}
/** @name Setup of Subtracted ME 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; }
//@}
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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans);
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
//@}
private:
/**
* Pointer to the head real emission ME casted to a MatchboxMEBase
* object.
*/
Ptr<MatchboxMEBase>::ptr theReal;
/**
* The dipoles to be considered; the dipoles generated
* can be accessed throught the dependent() matrxi element
* vector, provided the head() is a MatchboxMEBase object.
*/
vector<Ptr<SubtractionDipole>::ptr> theDipoles;
/**
* The underlying Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theBorns;
/**
* File name to dump subtraction check to
*/
string theSubtractionData;
/**
* Simple envelope histogram to keep track of subtraction
*/
struct SubtractionHistogram {
/**
* The lower bound
*/
double lower;
/**
* The bins, indexed by upper bound.
*/
map<double,pair<double,double> > bins;
/**
* Constructor
*/
SubtractionHistogram(double low = 0.001,
double up = 10.,
unsigned int nbins = 100);
/**
* Book an event.
*/
void book(double inv, double ratio) {
map<double,pair<double,double> >::iterator b =
bins.upper_bound(inv);
if ( b == bins.end() ) return;
b->second.first = min(b->second.first,abs(ratio));
b->second.second = max(b->second.second,abs(ratio));
}
/**
* Write to file given name and invariant.
*/
void dump(const std::string& prefix,
const cPDVector& proc,
int i, int j) const;
};
/**
* Define the key for the collinear subtraction data.
*/
typedef pair<cPDVector,pair<size_t, size_t> > CollinearSubtractionIndex;
/**
* subtraction data for collinear limits.
*/
map<CollinearSubtractionIndex,SubtractionHistogram> collinearHistograms;
/**
* Define the key for the soft subtraction data.
*/
typedef pair<cPDVector,size_t> SoftSubtractionIndex;
/**
* subtraction data for soft limits.
*/
map<SoftSubtractionIndex,SubtractionHistogram> softHistograms;
/**
* Switch to print full information on the
* last phase space point.
*/
bool theVerbose;
/**
* 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;
/**
* True, if veto scales should be set
* for the real emission
*/
bool theVetoScales;
/**
* The shower approximation.
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* True, if the shower real emission contribution should be subtracted.
*/
bool theRealShowerSubtraction;
/**
* True, if the shower virtual contribution should be subtracted.
*/
bool theVirtualShowerSubtraction;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SubtractedME & operator=(const SubtractedME &);
};
}
#endif /* HERWIG_SubtractedME_H */
diff --git a/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEPP2llbarJet.h b/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEPP2llbarJet.h
--- a/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEPP2llbarJet.h
+++ b/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEPP2llbarJet.h
@@ -1,237 +1,237 @@
// -*- C++ -*-
//
// MatchboxMEPP2llbarJet.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_MatchboxMEPP2llbarJet_H
#define HERWIG_MatchboxMEPP2llbarJet_H
//
// This is the declaration of the MatchboxMEPP2llbarJet class.
//
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbarqqbarg.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxMEPP2llbarJet implements charged lepton
* pair plus jet production in hadron-hadron collisions.
*
* @see \ref MatchboxMEPP2llbarJetInterfaces "The interfaces"
* defined for MatchboxMEPP2llbarJet.
*/
class MatchboxMEPP2llbarJet: public MatchboxMEBase, public MatchboxMEllbarqqbarg {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxMEPP2llbarJet();
/**
* The destructor.
*/
virtual ~MatchboxMEPP2llbarJet();
//@}
public:
/**
* The number of internal degreed of freedom used in the matrix
* element.
*/
- virtual int nDim() const { return phasespace() ? phasespace()->nDim(3) : 5; }
+ virtual int nDimBorn() const { return phasespace() ? phasespace()->nDim(3) : 5; }
/**
* Return the order in \f$\alpha_S\f$ in which this matrix element
* is given.
*/
virtual unsigned int orderInAlphaS() const { return 1; }
/**
* Return the order in \f$\alpha_{EM}\f$ in which this matrix
* element is given. Returns 0.
*/
virtual unsigned int orderInAlphaEW() const { return 2; }
/**
* 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);
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 factorizationScale() const;
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScale() const;
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
virtual unsigned int nLight() const { return theQuarkFlavours.size(); }
/**
* 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 and spin correlated matrix element squared for
* the gluon indexed by the first argument with and the
* given correlation tensor build as p1^\mu p2^\nu
*/
virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& p2) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* The lepton flavours to be considered.
*/
PDVector theLeptonFlavours;
/**
* The quark flavours to be considered.
*/
PDVector theQuarkFlavours;
/**
* The Z mass squared
*/
Energy2 theZMass2;
/**
* The Z width squared
*/
Energy2 theZWidth2;
private:
/**
* A user defined scale; if zero, the center of mass
* energy is used.
*/
Energy theUserScale;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class with persistent data.
*/
static AbstractClassDescription<MatchboxMEPP2llbarJet> initMatchboxMEPP2llbarJet;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxMEPP2llbarJet & operator=(const MatchboxMEPP2llbarJet &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MatchboxMEPP2llbarJet. */
template <>
struct BaseClassTrait<Herwig::MatchboxMEPP2llbarJet,1> {
/** Typedef of the first base class of MatchboxMEPP2llbarJet. */
typedef Herwig::MatchboxMEBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MatchboxMEPP2llbarJet class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MatchboxMEPP2llbarJet>
: public ClassTraitsBase<Herwig::MatchboxMEPP2llbarJet> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MatchboxMEPP2llbarJet"; }
/**
* The name of a file containing the dynamic library where the class
* MatchboxMEPP2llbarJet is implemented. It may also include several, space-separated,
* libraries if the class MatchboxMEPP2llbarJet depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwMatchbox.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MatchboxMEPP2llbarJet_H */
diff --git a/MatrixElement/Matchbox/Builtin/Processes/MatchboxMElP2lJetJet.h b/MatrixElement/Matchbox/Builtin/Processes/MatchboxMElP2lJetJet.h
--- a/MatrixElement/Matchbox/Builtin/Processes/MatchboxMElP2lJetJet.h
+++ b/MatrixElement/Matchbox/Builtin/Processes/MatchboxMElP2lJetJet.h
@@ -1,243 +1,243 @@
// -*- C++ -*-
//
// MatchboxMElP2lJetJet.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MatchboxMElP2lJetJet_H
#define HERWIG_MatchboxMElP2lJetJet_H
//
// This is the declaration of the MatchboxMElP2lJetJet class.
//
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbarqqbarg.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer, Luca D'Errico
*
* \brief MatchboxMElP2lJetJet implements the matrix element
* for neutral current, charged lepton 1+1 jet DIS with light quarks.
*
* @see \ref MatchboxMElP2lJetJetInterfaces "The interfaces"
* defined for MatchboxMElP2lJetJet.
*/
class MatchboxMElP2lJetJet: public MatchboxMEBase, public MatchboxMEllbarqqbarg {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxMElP2lJetJet();
/**
* The destructor.
*/
virtual ~MatchboxMElP2lJetJet();
//@}
public:
/**
* The number of internal degreed of freedom used in the matrix
* element.
*/
- virtual int nDim() const { return phasespace() ? phasespace()->nDim(3) : 5; }
+ virtual int nDimBorn() const { return phasespace() ? phasespace()->nDim(3) : 5; }
/**
* Return the order in \f$\alpha_S\f$ in which this matrix element
* is given.
*/
virtual unsigned int orderInAlphaS() const { return 1; }
/**
* Return the order in \f$\alpha_{EM}\f$ in which this matrix
* element is given. Returns 0.
*/
virtual unsigned int orderInAlphaEW() const { return 2; }
/**
* 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);
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 factorizationScale() const;
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScale() const;
/**
* Return true, if this matrix element provides the PDF
* weight for the first incoming parton itself.
*/
virtual bool havePDFWeight1() const { return false; }
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
virtual unsigned int nLight() const { return theQuarkFlavours.size(); }
/**
* 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 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;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* The lepton flavours to be considered.
*/
PDVector theLeptonFlavours;
/**
* The quark flavours to be considered.
*/
PDVector theQuarkFlavours;
/**
* The Z mass squared
*/
Energy2 theZMass2;
/**
* The Z width squared
*/
Energy2 theZWidth2;
private:
/**
* A user defined scale; if zero, the center of mass
* energy is used.
*/
Energy theUserScale;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class with persistent data.
*/
static AbstractClassDescription<MatchboxMElP2lJetJet> initMatchboxMElP2lJetJet;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxMElP2lJetJet & operator=(const MatchboxMElP2lJetJet &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MatchboxMElP2lJetJet. */
template <>
struct BaseClassTrait<Herwig::MatchboxMElP2lJetJet,1> {
/** Typedef of the first base class of MatchboxMElP2lJetJet. */
typedef Herwig::MatchboxMEBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MatchboxMElP2lJetJet class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MatchboxMElP2lJetJet>
: public ClassTraitsBase<Herwig::MatchboxMElP2lJetJet> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MatchboxMElP2lJetJet"; }
/**
* The name of a file containing the dynamic library where the class
* MatchboxMElP2lJetJet is implemented. It may also include several, space-separated,
* libraries if the class MatchboxMElP2lJetJet depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwMatchbox.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MatchboxMElP2lJetJet_H */
diff --git a/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbar2qqbarg.h b/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbar2qqbarg.h
--- a/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbar2qqbarg.h
+++ b/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbar2qqbarg.h
@@ -1,282 +1,282 @@
// -*- C++ -*-
//
// MatchboxMEllbar2qqbarg.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_MatchboxMEllbar2qqbarg_H
#define HERWIG_MatchboxMEllbar2qqbarg_H
//
// This is the declaration of the MatchboxMEllbar2qqbarg class.
//
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Builtin/Processes/MatchboxMEllbarqqbarg.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer, Martin Stoll
*
* \brief MatchboxMEllbar2qqbarg implements the
* matrix element for charged lepton annihilation
* into a light q qbar pair and a gluon.
*
* @see \ref MatchboxMEllbar2qqbargInterfaces "The interfaces"
* defined for MatchboxMEllbar2qqbarg.
*/
class MatchboxMEllbar2qqbarg: public MatchboxMEBase, public MatchboxMEllbarqqbarg {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxMEllbar2qqbarg();
/**
* The destructor.
*/
virtual ~MatchboxMEllbar2qqbarg();
//@}
public:
/**
* The number of internal degreed of freedom used in the matrix
* element.
*/
- virtual int nDim() const { return phasespace() ? phasespace()->nDim(3) : 5; }
+ virtual int nDimBorn() const { return phasespace() ? phasespace()->nDim(3) : 5; }
/**
* Return the order in \f$\alpha_S\f$ in which this matrix element
* is given.
*/
virtual unsigned int orderInAlphaS() const { return 1; }
/**
* Return the order in \f$\alpha_{EM}\f$ in which this matrix
* element is given. Returns 0.
*/
virtual unsigned int orderInAlphaEW() const { return 2; }
/**
* 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);
/**
* 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 renormalization scale for the last generated phasespace point.
*/
virtual Energy2 factorizationScale() const;
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScale() const;
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
virtual unsigned int nLight() const { return theQuarkFlavours.size(); }
/**
* Return true, if this matrix element provides the PDF
* weight for the first incoming parton itself.
*/
virtual bool havePDFWeight1() const { return false; }
/**
* Return true, if this matrix element provides the PDF
* weight for the second incoming parton itself.
*/
virtual bool havePDFWeight2() 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 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;
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* 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;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
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 lepton flavours to be considered.
*/
PDVector theLeptonFlavours;
/**
* The quark flavours to be considered.
*/
PDVector theQuarkFlavours;
/**
* A user defined scale; if zero, the center of mass
* energy is used.
*/
Energy theUserScale;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MatchboxMEllbar2qqbarg> initMatchboxMEllbar2qqbarg;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxMEllbar2qqbarg & operator=(const MatchboxMEllbar2qqbarg &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MatchboxMEllbar2qqbarg. */
template <>
struct BaseClassTrait<Herwig::MatchboxMEllbar2qqbarg,1> {
/** Typedef of the first base class of MatchboxMEllbar2qqbarg. */
typedef Herwig::MatchboxMEBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MatchboxMEllbar2qqbarg class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MatchboxMEllbar2qqbarg>
: public ClassTraitsBase<Herwig::MatchboxMEllbar2qqbarg> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MatchboxMEllbar2qqbarg"; }
/**
* The name of a file containing the dynamic library where the class
* MatchboxMEllbar2qqbarg is implemented. It may also include several, space-separated,
* libraries if the class MatchboxMEllbar2qqbarg depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwMatchbox.so"; }
};
/** @endcond */
}
#endif /* HERWIG_MatchboxMEllbar2qqbarg_H */
diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
@@ -1,1040 +1,1040 @@
// -*- C++ -*-
//
// SubtractionDipole.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 SubtractionDipole class.
//
#include "SubtractionDipole.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
SubtractionDipole::SubtractionDipole()
: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
theIgnoreCuts(false), theShowerKernel(false),
theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
theBornEmitter(-1), theBornSpectator(-1),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false) {}
SubtractionDipole::~SubtractionDipole() {}
void SubtractionDipole::setupBookkeeping() {
theMergingMap.clear();
theSplittingMap.clear();
int xemitter = -1;
int xspectator = -1;
map<int,int> mergeLegs;
map<int,int> remapLegs;
map<int,int> realBornMap;
map<int,int> bornRealMap;
for ( DiagramVector::const_iterator d =
theRealEmissionME->diagrams().begin();
d != theRealEmissionME->diagrams().end(); ++d ) {
Ptr<Tree2toNDiagram>::ptr check =
new_ptr(*dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d));
map<int,int> theMergeLegs;
map<int,int> theRemapLegs;
for ( unsigned int i = 0; i < check->external().size(); ++i )
theMergeLegs[i] = -1;
int theEmitter = check->mergeEmission(realEmitter(),realEmission(),theMergeLegs);
// no underlying Born
if ( theEmitter == -1 )
continue;
DiagramVector::const_iterator bd =
theUnderlyingBornME->diagrams().end();
for ( DiagramVector::const_iterator b =
theUnderlyingBornME->diagrams().begin();
b != theUnderlyingBornME->diagrams().end(); ++b )
if ( check->isSame(*b,theRemapLegs) ) {
bd = b; break;
}
// no underlying Born
if ( bd == theUnderlyingBornME->diagrams().end() )
continue;
if ( xemitter == -1 ) {
xemitter = theEmitter;
mergeLegs = theMergeLegs;
remapLegs = theRemapLegs;
assert(remapLegs.find(xemitter) != remapLegs.end());
xemitter = remapLegs[xemitter];
// work out the leg remapping real -> born
for ( map<int,int>::const_iterator k = mergeLegs.begin();
k != mergeLegs.end(); ++k ) {
assert(remapLegs.find(k->second) != remapLegs.end());
realBornMap[k->first] = remapLegs[k->second];
}
// work out the leg remapping born -> real
for ( map<int,int>::const_iterator k = realBornMap.begin();
k != realBornMap.end(); ++k ) {
bornRealMap[k->second] = k->first;
}
// work out the spectator
assert(mergeLegs.find(realSpectator()) != mergeLegs.end());
assert(remapLegs.find(mergeLegs[realSpectator()]) != remapLegs.end());
xspectator = realBornMap[realSpectator()];
}
RealEmissionKey realKey = realEmissionKey((**d).partons(),realEmitter(),realEmission(),realSpectator());
UnderlyingBornKey bornKey = underlyingBornKey((**bd).partons(),xemitter,xspectator);
if ( theMergingMap.find(realKey) == theMergingMap.end() )
theMergingMap.insert(make_pair(realKey,make_pair(bornKey,realBornMap)));
theSplittingMap.insert(make_pair(bornKey,make_pair(realKey,bornRealMap)));
}
if ( theSplittingMap.empty() )
return;
theIndexMap.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator s =
theSplittingMap.begin(); s != theSplittingMap.end(); ++s ) {
theIndexMap[process(s->first)] = make_pair(emitter(s->first),spectator(s->first));
}
theUnderlyingBornDiagrams.clear();
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator r =
theMergingMap.begin(); r != theMergingMap.end(); ++r ) {
DiagramVector borns;
for ( DiagramVector::const_iterator d = theUnderlyingBornME->diagrams().begin();
d != theUnderlyingBornME->diagrams().end(); ++d ) {
if ( (**d).partons() == process(r->second.first) )
borns.push_back(*d);
}
theUnderlyingBornDiagrams[process(r->first)] = borns;
}
theRealEmissionDiagrams.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator r =
theSplittingMap.begin(); r != theSplittingMap.end(); ++r ) {
DiagramVector reals;
for ( DiagramVector::const_iterator d = theRealEmissionME->diagrams().begin();
d != theRealEmissionME->diagrams().end(); ++d ) {
if ( (**d).partons() == process(r->second.first) )
reals.push_back(*d);
}
theRealEmissionDiagrams[process(r->first)] = reals;
}
}
void SubtractionDipole::subtractionBookkeeping() {
if ( theMergingMap.empty() )
setupBookkeeping();
lastRealEmissionKey =
realEmissionKey(lastHeadXComb().mePartonData(),realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() ) {
theApply = false;
return;
}
theApply = true;
lastUnderlyingBornKey = k->second.first;
bornEmitter(emitter(lastUnderlyingBornKey));
bornSpectator(spectator(lastUnderlyingBornKey));
}
void SubtractionDipole::splittingBookkeeping() {
if ( theMergingMap.empty() )
setupBookkeeping();
map<cPDVector,pair<int,int> >::const_iterator esit =
theIndexMap.find(lastHeadXComb().mePartonData());
if ( esit == theIndexMap.end() ) {
theApply = false;
return;
}
theApply = true;
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(lastHeadXComb().mePartonData(),bornEmitter(),bornSpectator());
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
assert(kr.first != kr.second);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
assert(lastRealEmissionInfo != kr.second);
lastRealEmissionKey = lastRealEmissionInfo->second.first;
realEmitter(emitter(lastRealEmissionKey));
realEmission(emission(lastRealEmissionKey));
realSpectator(spectator(lastRealEmissionKey));
}
StdXCombPtr SubtractionDipole::makeBornXComb(tStdXCombPtr realXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();
if ( theMergingMap.find(realEmissionKey(proc,realEmitter(),
realEmission(),realSpectator())) ==
theMergingMap.end() )
return StdXCombPtr();
PartonPairVec pbs = realXC->pExtractor()->getPartons(realXC->maxEnergy(),
realXC->particles(),
*(realXC->cuts()));
DiagramVector bornDiags = underlyingBornDiagrams(proc);
assert(!bornDiags.empty());
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == bornDiags.front()->partons()[0] &&
ppit->second->parton() == bornDiags.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
StdXCombPtr res =
new_ptr(StandardXComb(realXC,*ppit,this,bornDiags));
return res;
}
vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();
if ( theSplittingMap.find(underlyingBornKey(proc,bornEmitter(),bornSpectator())) ==
theSplittingMap.end() )
return vector<StdXCombPtr>();
PartonPairVec pbs = bornXC->pExtractor()->getPartons(bornXC->maxEnergy(),
bornXC->particles(),
*(bornXC->cuts()));
DiagramVector realDiags = realEmissionDiagrams(proc);
assert(!realDiags.empty());
vector<StdXCombPtr> res;
map<cPDVector,DiagramVector> realProcs;
for ( MEBase::DiagramVector::const_iterator d = realDiags.begin();
d != realDiags.end(); ++d ) {
realProcs[(**d).partons()].push_back(*d);
}
for ( map<cPDVector,DiagramVector>::const_iterator pr =
realProcs.begin(); pr != realProcs.end(); ++pr ) {
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == pr->second.front()->partons()[0] &&
ppit->second->parton() == pr->second.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
StdXCombPtr rxc =
new_ptr(StandardXComb(bornXC,*ppit,this,pr->second));
res.push_back(rxc);
}
return res;
}
const MEBase::DiagramVector& SubtractionDipole::underlyingBornDiagrams(const cPDVector& real) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theUnderlyingBornDiagrams.find(real);
if (k == theUnderlyingBornDiagrams.end() )
return empty;
return k->second;
}
const MEBase::DiagramVector& SubtractionDipole::realEmissionDiagrams(const cPDVector& born) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theRealEmissionDiagrams.find(born);
if ( k == theRealEmissionDiagrams.end() )
return empty;
return k->second;
}
void SubtractionDipole::getDiagrams() const {
if ( splitting() ) {
realEmissionME()->diagrams();
useDiagrams(realEmissionME());
} else {
underlyingBornME()->diagrams();
useDiagrams(underlyingBornME());
}
}
Selector<MEBase::DiagramIndex> SubtractionDipole::diagrams(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->prepare(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagrams(dv);
}
MEBase::DiagramIndex SubtractionDipole::diagram(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->prepare(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagram(dv);
}
Selector<const ColourLines *>
SubtractionDipole::colourGeometries(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->colourGeometries(diag) :
underlyingBornME()->colourGeometries(diag);
}
const ColourLines &
SubtractionDipole::selectColourGeometry(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->selectColourGeometry(diag) :
underlyingBornME()->selectColourGeometry(diag);
}
void SubtractionDipole::flushCaches() {
theUnderlyingBornME->flushCaches();
theRealEmissionME->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
}
void SubtractionDipole::setXComb(tStdXCombPtr xc) {
if ( !xc ) {
theApply = false;
return;
} else {
theApply = true;
}
MEBase::setXComb(xc);
if ( splitting() ) {
realEmissionME()->setXComb(xc);
underlyingBornME()->setXComb(xc->head());
splittingBookkeeping();
} else {
realEmissionME()->setXComb(xc->head());
underlyingBornME()->setXComb(xc);
subtractionBookkeeping();
}
if ( !apply() )
return;
}
void SubtractionDipole::setKinematics() {
MEBase::setKinematics();
if ( splitting() )
realEmissionME()->setKinematics();
else
underlyingBornME()->setKinematics();
}
bool SubtractionDipole::generateKinematics(const double * r) {
+ if ( lastXCombPtr()->kinematicsGenerated() )
+ return true;
if ( splitting() ) {
if ( !generateRadiationKinematics(r) )
return false;
realEmissionME()->lastXCombPtr()->setIncomingPartons();
realEmissionME()->setScale();
+ double jac = jacobian();
+ jac *= pow(underlyingBornME()->lastXComb().lastSHat() / realEmissionME()->lastXComb().lastSHat(),
+ realEmissionME()->lastXComb().mePartonData().size()-4.);
+ jacobian(jac);
assert(lastXCombPtr() == realEmissionME()->lastXCombPtr());
return true;
}
if ( !generateTildeKinematics() )
return false;
underlyingBornME()->setScale();
assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
underlyingBornME()->lastXCombPtr()->setIncomingPartons();
+ lastXCombPtr()->didGenerateKinematics();
return true;
}
int SubtractionDipole::nDim() const {
if ( !splitting() )
return underlyingBornME()->nDim();
- int alldim =
- underlyingBornME()->nDim() + nDimRadiation();
- // don't mess with random number for collinear remainders
- // see MEGroup.cc:55 and MatchboxNLOME
- if ( underlyingBornME()->diagrams().front()->partons()[0]->coloured() ||
- underlyingBornME()->diagrams().front()->partons()[1]->coloured() )
- ++alldim;
- return alldim;
+ return underlyingBornME()->nDim() + nDimRadiation();
}
void SubtractionDipole::clearKinematics() {
MEBase::clearKinematics();
if ( splitting() )
realEmissionME()->clearKinematics();
else
underlyingBornME()->clearKinematics();
}
void SubtractionDipole::tildeKinematics(Ptr<TildeKinematics>::tptr tk) {
theTildeKinematics = tk;
}
bool SubtractionDipole::generateTildeKinematics() {
assert(!splitting());
if ( !tildeKinematics() ) {
jacobian(0.0);
return false;
}
theTildeKinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());
if ( !theTildeKinematics->doMap() ) {
jacobian(0.0);
return false;
}
theLastSubtractionScale = theTildeKinematics->lastScale();
theLastSubtractionPt = theTildeKinematics->lastPt();
meMomenta().resize(lastHeadXComb().meMomenta().size() - 1);
assert(mergingMap().find(lastRealEmissionKey) != mergingMap().end());
map<int,int>& trans = theMergingMap[lastRealEmissionKey].second;
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == realEmitter() || k == realEmission() || k == realSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( theTildeKinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = theTildeKinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[bornEmitter()] = tildeKinematics()->bornEmitterMomentum();
meMomenta()[bornSpectator()] = tildeKinematics()->bornSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).mass());
p->rescaleRho();
}
jacobian(realEmissionME()->lastXComb().jacobian());
logGenerateTildeKinematics();
return true;
}
void SubtractionDipole::invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr itk) {
theInvertedTildeKinematics = itk;
}
int SubtractionDipole::nDimRadiation() const {
return invertedTildeKinematics() ?
invertedTildeKinematics()->nDimRadiation() :
0;
}
bool SubtractionDipole::generateRadiationKinematics(const double * r) {
assert(splitting());
if ( !invertedTildeKinematics() ) {
jacobian(0.0);
return false;
}
theInvertedTildeKinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());
if ( !theInvertedTildeKinematics->doMap(r) ) {
jacobian(0.0);
return false;
}
theLastSplittingScale = theInvertedTildeKinematics->lastScale();
theLastSplittingPt = theInvertedTildeKinematics->lastPt();
meMomenta().resize(lastHeadXComb().meMomenta().size() + 1);
assert(splittingMap().find(lastUnderlyingBornKey) != splittingMap().end());
map<int,int>& trans = const_cast<map<int,int>&>(lastRealEmissionInfo->second.second);
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == bornEmitter() || k == bornSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( invertedTildeKinematics()->doesTransform() && k > 1 )
meMomenta()[trans[k]] = invertedTildeKinematics()->transform(meMomenta()[trans[k]]);
}
meMomenta()[realEmitter()] = invertedTildeKinematics()->realEmitterMomentum();
meMomenta()[realEmission()] = invertedTildeKinematics()->realEmissionMomentum();
meMomenta()[realSpectator()] = invertedTildeKinematics()->realSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).mass());
p->rescaleRho();
}
jacobian(underlyingBornME()->lastXComb().jacobian() *
invertedTildeKinematics()->jacobian());
logGenerateRadiationKinematics(r);
return true;
}
void SubtractionDipole::ptCut(Energy cut) {
theInvertedTildeKinematics->ptCut(cut);
}
CrossSection SubtractionDipole::dSigHatDR(Energy2 factorizationScale) const {
+ double pdfweight = 1.;
+
+ double jac = jacobian();
+
+ if ( splitting() && jac == 0.0 ) {
+ lastMECrossSection(ZERO);
+ lastME2(0.0);
+ return ZERO;
+ }
+
+ if ( havePDFWeight1() || havePDFWeight2() ) {
+ if ( splitting() )
+ realEmissionME()->getPDFWeight(factorizationScale);
+ pdfweight = realEmissionME()->lastXComb().lastMEPDFWeight();
+ lastMEPDFWeight(pdfweight);
+ }
+
if ( showerApproximation() && realShowerSubtraction() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
if ( !showerApproximation()->isAboveCutoff() )
showerApproximation()->wasBelowCutoff();
if ( !showerApproximation()->isInShowerPhasespace() ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
- return -showerApproximation()->dSigHatDR();
- }
-
- double pdfweight = 1.;
-
- double jac = jacobian();
-
- if ( splitting() && jac == 0.0 )
- return ZERO;
-
- if ( splitting() )
- jac *= pow(underlyingBornME()->lastXComb().lastSHat() / realEmissionME()->lastXComb().lastSHat(),
- realEmissionME()->lastXComb().mePartonData().size()-4.);
-
- if ( havePDFWeight1() || havePDFWeight2() ) {
- if ( splitting() )
- realEmissionME()->getPDFWeight(factorizationScale);
- pdfweight = realEmissionME()->lastXComb().lastMEPDFWeight();
- lastMEPDFWeight(pdfweight);
- }
-
- CrossSection res = ZERO;
-
- if ( showerApproximation() && virtualShowerSubtraction() ) {
- if ( !splitting() ) {
- showerApproximation()->setBornXComb(lastXCombPtr());
- showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
- } else {
- showerApproximation()->setBornXComb(underlyingBornME()->lastXCombPtr());
- showerApproximation()->setRealXComb(lastXCombPtr());
- }
- showerApproximation()->setDipole(this);
- if ( !showerApproximation()->isAboveCutoff() )
- showerApproximation()->wasBelowCutoff();
- if ( showerApproximation()->isInShowerPhasespace() )
- res = showerApproximation()->dSigHatDR();
+ lastME2(me2());
+ lastMECrossSection(-showerApproximation()->dSigHatDR());
+ return lastMECrossSection();
}
double xme2 = 0.0;
if ( !showerKernel() )
xme2 = me2();
else
xme2 = me2Avg(-underlyingBornME()->me2());
if ( xme2 == 0.0 ) {
- lastMECrossSection(res);
+ lastMECrossSection(ZERO);
lastME2(0.0);
- return res;
+ return ZERO;
}
lastME2(xme2);
- res -=
+ CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
- double weight = 0.0;
- bool applied = false;
- for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
- theReweights.begin(); rw != theReweights.end(); ++rw ) {
- (**rw).setXComb(theRealEmissionME->lastXCombPtr());
- if ( !(**rw).apply() )
- continue;
- weight += (**rw).evaluate();
- applied = true;
+ if ( !showerApproximation() ) {
+ double weight = 0.0;
+ bool applied = false;
+ for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
+ theReweights.begin(); rw != theReweights.end(); ++rw ) {
+ (**rw).setXComb(theRealEmissionME->lastXCombPtr());
+ if ( !(**rw).apply() )
+ continue;
+ weight += (**rw).evaluate();
+ applied = true;
+ }
+ if ( applied )
+ res *= weight;
}
- if ( applied )
- res *= weight;
- lastMECrossSection(res);
+ if ( showerApproximation() && virtualShowerSubtraction() ) {
+ assert(!splitting());
+ showerApproximation()->setBornXComb(lastXCombPtr());
+ showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
+ showerApproximation()->setDipole(this);
+ if ( !showerApproximation()->isAboveCutoff() )
+ showerApproximation()->wasBelowCutoff();
+ CrossSection shower = ZERO;
+ if ( showerApproximation()->isInShowerPhasespace() )
+ shower = showerApproximation()->dSigHatDR();
+ vector<double> ratio(1,shower/res);
+ meInfo(ratio);
+ }
+
+ lastMECrossSection(-res);
logDSigHatDR(jac);
return lastMECrossSection();
}
void SubtractionDipole::print(ostream& os) const {
os << "--- SubtractionDipole setup ----------------------------------------------------\n";
os << " subtraction '" << name() << "'\n for real emission '"
<< theRealEmissionME->name() << "'\n using underlying Born '"
<< theUnderlyingBornME->name() << "'\n";
os << " tilde kinematics are '"
<< (theTildeKinematics ? theTildeKinematics->name() : "")
<< " '\n inverted tilde kinematics are '"
<< (theInvertedTildeKinematics ? theInvertedTildeKinematics->name() : "") << "'\n";
os << " the following subtraction mappings have been found:\n";
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator m =
theMergingMap.begin(); m != theMergingMap.end(); ++m ) {
os << " " << process(m->second.first)[0]->PDGName() << " "
<< process(m->second.first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->second.first).begin() + 2;
p != process(m->second.first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[" << emitter(m->second.first) << "," << spectator(m->second.first) << "] <=> ";
os << process(m->first)[0]->PDGName() << " "
<< process(m->first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->first).begin() + 2;
p != process(m->first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[(" << emitter(m->first) << "," << emission(m->first) << ")," << spectator(m->first) << "]\n"
<< " non-dipole momenta ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->second << " ";
}
os << ") <=> ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->first << " ";
}
os << ")\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractionDipole::printLastEvent(ostream& os) const {
os << "--- SubtractionDipole last event information -----------------------------------\n";
os << " for dipole '" << name() << "' applying ["
<< bornEmitter() << "," << bornSpectator() << "] <=> [("
<< realEmitter() << "," << realEmission() << ")," << realSpectator() << "]\n"
<< " evaluated the cross section/nb " << (lastMECrossSection()/nanobarn) << "\n"
<< " with subtraction parameters x[0] = " << subtractionParameters()[0]
<< " x[1] = " << subtractionParameters()[1] << "\n";
os << " the last real emission event was:\n";
realEmissionME()->printLastEvent(os);
os << " the last underlying Born event was:\n";
underlyingBornME()->printLastEvent(os);
os << "--- end SubtractionDipole last event information -------------------------------\n";
os << flush;
}
void SubtractionDipole::logME2() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated me2 using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "Born phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = bornxc->meMomenta().begin();
cPDVector::const_iterator dit = bornxc->mePartonData().begin();
for ( ; pit != bornxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << bornxc->lastX1() << " x2 = " << bornxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (bornxc->lastSHat()/GeV2) << "\n";
generator()->log() << "Real emission phase space point (in GeV):\n";
pit = realxc->meMomenta().begin();
dit = realxc->mePartonData().begin();
for ( ; pit != realxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << realxc->lastX1() << " x2 = " << realxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (realxc->lastSHat()/GeV2) << "\n";
generator()->log() << "me2 = " << lastME2() << "\n" << flush;
}
void SubtractionDipole::logDSigHatDR(double effectiveJac) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated cross section using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n"
<< "Jacobian = " << jacobian()
<< " effective Jacobian = " << effectiveJac << "\n"
<< "Born sHat/GeV2 = " << (bornxc->lastSHat()/GeV2)
<< " real sHat/GeV2 = " << (realxc->lastSHat()/GeV2)
<< " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void SubtractionDipole::logGenerateTildeKinematics() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating tilde kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with real xcomb " << lastHeadXCombPtr() << " born xcomb "
<< lastXCombPtr() << "\n"
<< "from real emission phase space point:\n";
Lorentz5Momentum rSum;
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
size_t count = 0;
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
rSum -= *pr;
} else {
rSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (rSum/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n"
<< "with scale/GeV = " << (theLastSubtractionScale/GeV)
<< "and pt/GeV = " << (theLastSubtractionPt/GeV) << "\n";
generator()->log() << "generated tilde kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
count = 0;
Lorentz5Momentum bSum;
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
bSum -= *pr;
} else {
bSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (bSum/GeV) << "\n";
generator()->log() << "Jacobian = " << jacobian() << "\n" << flush;
}
void SubtractionDipole::logGenerateRadiationKinematics(const double * r) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating radiation kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with born xcomb " << lastHeadXCombPtr() << " real xcomb "
<< lastXCombPtr() << "\n"
<< "from random numbers:\n";
copy(r,r+nDimRadiation(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "and born phase space point:\n";
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n" << flush;
generator()->log() << "scales: scale/GeV = " << (theLastSplittingScale/GeV)
<< " pt/GeV = " << (theLastSplittingPt/GeV) << "\n" << flush;
generator()->log() << "generated real emission kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "Jacobian = "
<< jacobian() << " = "
<< underlyingBornME()->lastXComb().jacobian()
<< "|Born * "
<< invertedTildeKinematics()->jacobian()
<< "|Radiation\n" << flush;
}
void SubtractionDipole::cloneDependencies(const std::string& prefix) {
if ( underlyingBornME() ) {
Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myUnderlyingBornME->cloneDependencies(pname.str());
underlyingBornME(myUnderlyingBornME);
}
if ( realEmissionME() ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = realEmissionME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myRealEmissionME->cloneDependencies(pname.str());
realEmissionME(myRealEmissionME);
}
if ( tildeKinematics() ) {
Ptr<TildeKinematics>::ptr myTildeKinematics = tildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myTildeKinematics->name();
if ( ! (generator()->preinitRegister(myTildeKinematics,pname.str()) ) )
throw InitException() << "Tilde kinematics " << pname.str() << " already existing.";
myTildeKinematics->dipole(this);
tildeKinematics(myTildeKinematics);
}
if ( invertedTildeKinematics() ) {
Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
throw InitException() << "Inverted tilde kinematics " << pname.str() << " already existing.";
myInvertedTildeKinematics->dipole(this);
invertedTildeKinematics(myInvertedTildeKinematics);
}
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;
}
}
void SubtractionDipole::constructVertex(tSubProPtr sub) {
if ( splitting() )
realEmissionME()->constructVertex(sub);
else
underlyingBornME()->constructVertex(sub);
}
void SubtractionDipole::generateSubCollision(SubProcess & sub) {
if ( splitting() )
realEmissionME()->generateSubCollision(sub);
else
underlyingBornME()->generateSubCollision(sub);
}
void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
os << theSplitting << theApply << theSubtractionTest << theIgnoreCuts
<< theShowerKernel << theRealEmissionME << theUnderlyingBornME
<< theTildeKinematics << theInvertedTildeKinematics << theReweights
<< theRealEmitter << theRealEmission << theRealSpectator
<< theSubtractionParameters
<< theMergingMap << theSplittingMap << theIndexMap
<< theUnderlyingBornDiagrams << theRealEmissionDiagrams
<< lastRealEmissionKey << lastUnderlyingBornKey
<< theBornEmitter << theBornSpectator
<< theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction
<< thePartners;
}
void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
is >> theSplitting >> theApply >> theSubtractionTest >> theIgnoreCuts
>> theShowerKernel >> theRealEmissionME >> theUnderlyingBornME
>> theTildeKinematics >> theInvertedTildeKinematics >> theReweights
>> theRealEmitter >> theRealEmission >> theRealSpectator
>> theSubtractionParameters
>> theMergingMap >> theSplittingMap >> theIndexMap
>> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
>> lastRealEmissionKey >> lastUnderlyingBornKey
>> theBornEmitter >> theBornSpectator
>> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction
>> thePartners;
}
void SubtractionDipole::Init() {
static ClassDocumentation<SubtractionDipole> documentation
("SubtractionDipole represents a dipole subtraction "
"term in the formalism of Catani and Seymour.");
static Reference<SubtractionDipole,MatchboxMEBase> interfaceUnderlyingBornME
("UnderlyingBornME",
"The underlying Born matrix element.",
&SubtractionDipole::theUnderlyingBornME, false, false, true, true, false);
static Reference<SubtractionDipole,MatchboxMEBase> interfaceRealEmissionME
("RealEmissionME",
"The real emission matrix element.",
&SubtractionDipole::theRealEmissionME, false, false, true, true, false);
static RefVector<SubtractionDipole,SubtractionDipole> interfacePartners
("Partners",
"The partner dipoles found along with this dipole.",
&SubtractionDipole::thePartners, -1, false, false, true, false, false);
static Reference<SubtractionDipole,TildeKinematics> interfaceTildeKinematics
("TildeKinematics",
"Set the TildeKinematics object to be used.",
&SubtractionDipole::theTildeKinematics, false, false, true, false, false);
static Reference<SubtractionDipole,InvertedTildeKinematics> interfaceInvertedTildeKinematics
("InvertedTildeKinematics",
"Set the InvertedTildeKinematics object to be used.",
&SubtractionDipole::theInvertedTildeKinematics, false, false, true, true, false);
static RefVector<SubtractionDipole,MatchboxReweightBase> interfaceReweights
("Reweights",
"Reweight objects to be applied to this matrix element.",
&SubtractionDipole::theReweights, -1, false, false, true, true, false);
static Reference<SubtractionDipole,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&SubtractionDipole::theShowerApproximation, false, false, true, true, 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).
DescribeAbstractClass<SubtractionDipole,MEBase>
describeSubtractionDipole("Herwig::SubtractionDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Makefile.am b/MatrixElement/Matchbox/Makefile.am
--- a/MatrixElement/Matchbox/Makefile.am
+++ b/MatrixElement/Matchbox/Makefile.am
@@ -1,23 +1,23 @@
SUBDIRS = \
Base Utility Phasespace \
-Dipoles InsertionOperators Powheg Matching \
-Builtin
+Dipoles InsertionOperators Matching \
+Builtin Tests
if WANT_DIPOLE
pkglib_LTLIBRARIES = HwMatchbox.la
endif
HwMatchbox_la_LDFLAGS = -module -version-info 2:0:0
HwMatchbox_la_LIBADD = \
Base/libHwMatchboxBase.la \
Utility/libHwMatchboxUtility.la \
Phasespace/libHwMatchboxPhasespace.la \
Dipoles/libHwMatchboxDipoles.la \
InsertionOperators/libHwMatchboxInsertionOperators.la \
-Powheg/libHwMatchboxPowheg.la \
Matching/libHwMatchboxMatching.la \
Builtin/Processes/libHwMatchboxBuiltinProcesses.la \
-Builtin/Amplitudes/libHwMatchboxBuiltinAmplitudes.la
+Builtin/Amplitudes/libHwMatchboxBuiltinAmplitudes.la \
+Tests/libHwMatchboxTests.la
HwMatchbox_la_SOURCES = \
MatchboxFactory.h MatchboxFactory.cc
diff --git a/MatrixElement/Matchbox/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc
--- a/MatrixElement/Matchbox/MatchboxFactory.cc
+++ b/MatrixElement/Matchbox/MatchboxFactory.cc
@@ -1,1097 +1,1142 @@
// -*- 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), theSubProcessGroups(false), theInclusive(false),
theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0),
theFixedCouplings(false), theFixedQEDCouplings(false), theVetoScales(false),
- theVerbose(false), theSubtractionData(""), theCheckPoles(false) {}
+ theVerbose(false), theInitVerbose(false), theSubtractionData(""), theCheckPoles(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);
if ( diagramGenerator() && !me->diagramGenerator() )
me->diagramGenerator(diagramGenerator());
if ( processData() && !me->processData() )
me->processData(processData());
if ( me->nLight() == 0 )
me->nLight(nLight());
if ( phasespace() && !me->phasespace() )
me->phasespace(phasespace());
if ( scaleChoice() && !me->scaleChoice() )
me->scaleChoice(scaleChoice());
if ( me->factorizationScaleFactor() == 1.0 )
me->factorizationScaleFactor(factorizationScaleFactor());
if ( me->renormalizationScaleFactor() == 1.0 )
me->renormalizationScaleFactor(renormalizationScaleFactor());
if ( fixedCouplings() )
me->setFixedCouplings();
if ( fixedQEDCouplings() )
me->setFixedQEDCouplings();
if ( cache() && !me->cache() )
me->cache(cache());
if ( verbose() )
me->setVerbose();
}
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);
vector<Ptr<MatchboxAmplitude>::ptr> matchAmplitudes;
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
(**amp).orderInGs(orderas);
(**amp).orderInGem(orderInAlphaEW());
if ( (**amp).orderInGs() != orderas ||
(**amp).orderInGem() != orderInAlphaEW() ) {
continue;
}
matchAmplitudes.push_back(*amp);
}
size_t combinations = processes.size()*matchAmplitudes.size();
size_t procCount = 0;
boost::progress_display * progressBar =
new boost::progress_display(combinations,generator()->log());
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= matchAmplitudes.begin(); amp != matchAmplitudes.end(); ++amp ) {
(**amp).orderInGs(orderas);
(**amp).orderInGem(orderInAlphaEW());
for ( set<PDVector>::const_iterator p = processes.begin();
p != processes.end(); ++p ) {
++(*progressBar);
if ( !(**amp).canHandle(*p) )
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);
string pname = "ME" + ap->first->name() + pid(m->first);
if ( ! (generator()->preinitRegister(me,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
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 ( !amplitudes().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() ) {
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() ) {
- throw InitException() << "Virtual corrections without Born contributions not yet supported.\n";
- }
-
if ( bornContributions() && !virtualContributions() ) {
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 ( ! (generator()->preinitRegister(bornme,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
bornme->cloneDependencies();
MEs().push_back(bornme);
}
}
- if ( bornContributions() && virtualContributions() ) {
+ if ( virtualContributions() ) {
generator()->log() << "preparing Born "
<< (virtualContributions() ? "and virtual" : "")
<< " matrix elements." << flush;
bornVirtualMEs().clear();
boost::progress_display * progressBar =
new boost::progress_display(bornMEs().size(),generator()->log());
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
- Ptr<MatchboxNLOME>::ptr nlo = new_ptr(MatchboxNLOME());
+ Ptr<MatchboxMEBase>::ptr nlo = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( ! (generator()->preinitRegister(nlo,pname) ) )
throw InitException() << "NLO ME " << pname << " already existing.";
- nlo->matrixElement(*born);
nlo->virtuals().clear();
- if ( !nlo->matrixElement()->onlyOneLoop() ) {
+ 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";
}
nlo->doCheckPoles();
}
}
+ if ( !bornContributions() ) {
+ 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() ) {
- 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 ( 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." << 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();
if ( ! (generator()->preinitRegister(sub,pname) ) )
throw InitException() << "Subtracted ME " << pname << " already existing.";
sub->borns() = bornMEs();
sub->head(*real);
sub->allDipoles().clear();
sub->dependent().clear();
if ( subtractionData() != "" )
sub->subtractionData(subtractionData());
sub->getDipoles();
if ( sub->dependent().empty() ) {
// finite real contribution
MEs().push_back(sub->head());
Ptr<MatchboxMEBase>::ptr fme = dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(sub->head());
finiteRealMEs().push_back(fme);
sub->head(tMEPtr());
continue;
}
if ( verbose() )
sub->setVerbose();
if ( subProcessGroups() )
sub->setSubProcessGroups();
if ( inclusive() )
sub->setInclusive();
if ( vetoScales() )
sub->doVetoScales();
subtractedMEs().push_back(sub);
MEs().push_back(sub);
if ( showerApproximation() ) {
sub->showerApproximation(showerApproximation());
Ptr<SubtractedME>::ptr subv = new_ptr(*sub);
string vname = sub->fullName() + ".vsub";
if ( ! (generator()->preinitRegister(subv,pname) ) )
throw InitException() << "Subtracted ME " << vname << " already existing.";
sub->doRealShowerSubtraction();
subv->doVirtualShowerSubtraction();
subtractedMEs().push_back(subv);
MEs().push_back(subv);
- 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()));
- }
+ 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;
}
+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;
+
+
+ for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator sd =
+ splitDipoles.begin(); sd != splitDipoles.end(); ++sd ) {
+
+ Ptr<MatchboxMEBase>::tptr bornME =
+ const_ptr_cast<Ptr<MatchboxMEBase>::tptr>((**sd).underlyingBornME());
+
+ SplittingChannel channel;
+ channel.dipole = *sd;
+ 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());
+
+ vector<StdXCombPtr> realXCombs = (**sd).makeRealXCombs(channel.bornXComb);
+
+ for ( vector<StdXCombPtr>::const_iterator rxc = realXCombs.begin();
+ rxc != realXCombs.end(); ++rxc ) {
+ channel.realXComb = *rxc;
+ res.push_back(channel);
+ }
+
+ }
+
+ 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 ) {
- os << " '" << (**m).name() << "' for subprocesses:\n";
- for ( vector<PDVector>::const_iterator p = (**m).subProcesses().begin();
- p != (**m).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";
- }
+ (**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 ) {
- os << " '" << (**m).name() << "' for subprocesses:\n";
- for ( vector<PDVector>::const_iterator p = (**m).subProcesses().begin();
- p != (**m).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";
- }
+ (**m).print(os);
}
os << flush;
}
os << " generated Born+virtual matrix elements:\n";
- for ( vector<Ptr<MatchboxNLOME>::ptr>::const_iterator bv
+ 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 ( theVerbose )
+ if ( initVerbose() )
print(Repository::clog());
SubProcessHandler::doinit();
}
void MatchboxFactory::persistentOutput(PersistentOStream & os) const {
os << theDiagramGenerator << theProcessData
<< theNLight << theOrderInAlphaS << theOrderInAlphaEW
<< theBornContributions << theVirtualContributions
<< theRealContributions << theSubProcessGroups << theInclusive
<< thePhasespace << theScaleChoice
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theFixedCouplings << theFixedQEDCouplings << theVetoScales
<< theAmplitudes << theCache
<< theBornMEs << theVirtuals << theRealEmissionMEs
<< theBornVirtualMEs << theSubtractedMEs << theFiniteRealMEs
- << theVerbose << theSubtractionData << theCheckPoles
+ << theVerbose << theInitVerbose << theSubtractionData << theCheckPoles
<< theParticleGroups << process << realEmissionProcess
<< theShowerApproximation << theSplittingDipoles;
}
void MatchboxFactory::persistentInput(PersistentIStream & is, int) {
is >> theDiagramGenerator >> theProcessData
>> theNLight >> theOrderInAlphaS >> theOrderInAlphaEW
>> theBornContributions >> theVirtualContributions
>> theRealContributions >> theSubProcessGroups >> theInclusive
>> thePhasespace >> theScaleChoice
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theFixedCouplings >> theFixedQEDCouplings >> theVetoScales
>> theAmplitudes >> theCache
>> theBornMEs >> theVirtuals >> theRealEmissionMEs
>> theBornVirtualMEs >> theSubtractedMEs >> theFiniteRealMEs
- >> theVerbose >> theSubtractionData >> theCheckPoles
+ >> theVerbose >> theInitVerbose >> theSubtractionData >> theCheckPoles
>> theParticleGroups >> process >> realEmissionProcess
>> theShowerApproximation >> theSplittingDipoles;
}
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) 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;
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> 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 Reference<MatchboxFactory,MatchboxMECache> interfaceCache
("Cache",
"Set the matrix element cache object.",
&MatchboxFactory::theCache, 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,MatchboxNLOME> interfaceBornVirtuals
+ 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 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);
}
// *** 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,719 +1,761 @@
// -*- 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/Utility/MatchboxMECache.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
-#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxNLOME.h"
#include "Herwig++/MatrixElement/Matchbox/Base/SubtractedME.h"
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 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 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 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; }
/**
* Set the ME cache object
*/
void cache(Ptr<MatchboxMECache>::ptr c) { theCache = c; }
/**
* Get the ME cache object
*/
Ptr<MatchboxMECache>::tptr cache() const { return theCache; }
//@}
/** @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<MatchboxNLOME>::ptr>& bornVirtualMEs() const { return theBornVirtualMEs; }
+ const vector<Ptr<MatchboxMEBase>::ptr>& bornVirtualMEs() const { return theBornVirtualMEs; }
/**
* Access the produced NLO matrix elements
*/
- vector<Ptr<MatchboxNLOME>::ptr>& bornVirtualMEs() { return theBornVirtualMEs; }
+ 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 dipole in charge of the splitting
+ */
+ Ptr<SubtractionDipole>::ptr dipole;
+
+ };
+
+ /**
+ * 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
*/
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 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; }
//@}
/** @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 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 ME cache object
*/
Ptr<MatchboxMECache>::ptr theCache;
/**
* 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<MatchboxNLOME>::ptr> theBornVirtualMEs;
+ 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;
/**
* 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>&) const;
/**
* 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;
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 */
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximation.h b/MatrixElement/Matchbox/Matching/ShowerApproximation.h
--- a/MatrixElement/Matchbox/Matching/ShowerApproximation.h
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximation.h
@@ -1,194 +1,202 @@
// -*- C++ -*-
//
// ShowerApproximation.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_ShowerApproximation_H
#define Herwig_ShowerApproximation_H
//
// This is the declaration of the ShowerApproximation class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief ShowerApproximation describes the shower emission to be used
* in NLO matching.
*
*/
class ShowerApproximation: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ShowerApproximation();
/**
* The destructor.
*/
virtual ~ShowerApproximation();
//@}
public:
/**
+ * Return true, if this shower approximation will require a
+ * splitting generator
+ */
+ virtual bool needsSplittingGenerator() const { return false; }
+
+public:
+
+ /**
* Set the XComb object describing the Born process
*/
void setBornXComb(tStdXCombPtr xc) { theBornXComb = xc; }
/**
* Return the XComb object describing the Born process
*/
tStdXCombPtr bornXComb() const { return theBornXComb; }
/**
* Set the XComb object describing the real emission process
*/
void setRealXComb(tStdXCombPtr xc) { theRealXComb = xc; }
/**
* Return the XComb object describing the real emission process
*/
tStdXCombPtr realXComb() const { return theRealXComb; }
/**
* Set the dipole in charge for the emission
*/
void setDipole(Ptr<SubtractionDipole>::tcptr);
/**
* Return the dipole in charge for the emission
*/
Ptr<SubtractionDipole>::tcptr dipole() const;
public:
/**
* Return true if one of the recently encountered configutations was
* below the infrared cutoff.
*/
bool belowCutoff() const { return theBelowCutoff; }
/**
* Indicate that one of the recently encountered configutations was
* below the infrared cutoff.
*/
void wasBelowCutoff() { theBelowCutoff = true; }
/**
* Reset the below cutoff flag.
*/
void resetBelowCutoff() { theBelowCutoff = false; }
public:
/**
* Return true, if the shower was able to generate an emission
* leading from the given Born to the given real emission process.
*/
virtual bool isInShowerPhasespace() const = 0;
/**
* Return true, if the shower emission leading from the given Born
* to the given real emission process would have been generated
* above the shower's infrared cutoff.
*/
virtual bool isAboveCutoff() const = 0;
/**
* Return the shower approximation to the real emission cross
- * section for the given pair of Bron and real emission
+ * section for the given pair of Born and real emission
* configurations.
*/
virtual CrossSection dSigHatDR() const = 0;
/**
* Return the shower approximated real emission cross section
* divided by the Born cross section and the radiation phasespace
- * jacobian for the given pair of Bron and real emission
+ * jacobian for the given pair of Born and real emission
* configurations.
*/
virtual double me2() const = 0;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The XComb object describing the Born process
*/
tStdXCombPtr theBornXComb;
/**
* The XComb object describing the real emission process
*/
tStdXCombPtr theRealXComb;
/**
* The dipole in charge for the emission
*/
Ptr<SubtractionDipole>::tcptr theDipole;
/**
* True if one of the recently encountered configutations was below
* the infrared cutoff.
*/
bool theBelowCutoff;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerApproximation & operator=(const ShowerApproximation &);
};
}
#endif /* Herwig_ShowerApproximation_H */
diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
@@ -0,0 +1,422 @@
+// -*- C++ -*-
+//
+// FlatInvertiblePhasespace.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 FlatInvertiblePhasespace class.
+//
+
+#include "FlatInvertiblePhasespace.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "Herwig++/Utilities/GSLBisection.h"
+#include "ThePEG/Cuts/Cuts.h"
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+FlatInvertiblePhasespace::FlatInvertiblePhasespace() {}
+
+FlatInvertiblePhasespace::~FlatInvertiblePhasespace() {}
+
+IBPtr FlatInvertiblePhasespace::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr FlatInvertiblePhasespace::fullclone() const {
+ return new_ptr(*this);
+}
+
+void FlatInvertiblePhasespace::prepare(tStdXCombPtr xc, bool) {
+ theLastXComb = xc;
+}
+
+double FlatInvertiblePhasespace::bisect(double v, double n,
+ double target, double maxLevel) const {
+
+ if ( v != 0.0 && v != 1.0 ) {
+
+ double level = 0;
+ double left = 0;
+ double right = 1;
+
+ double checkV = -1.;
+ double u = -1;
+
+ while ( level < maxLevel ) {
+
+ u = (left+right)*pow(0.5,level+1.);
+ checkV =
+ pow(u,n+1.)*(n+2.-(n+1.)*u);
+
+ if ( log10(abs(1.-checkV/v)) <= target )
+ break;
+
+ left *= 2.;
+ right *= 2.;
+
+ if ( v <= checkV ) {
+ right -= 1.;
+ ++level;
+ }
+
+ if ( v > checkV ) {
+ left += 1.;
+ ++level;
+ }
+
+ }
+
+ return u;
+
+ }
+
+ return v;
+
+}
+
+static double flatWeights[7] = {
+
+ -1.,-1.,
+ 0.039788735772973833942,
+ 0.00012598255637968550463,
+ 1.3296564302788840628E-7,
+ 7.0167897579949011130E-11,
+ 2.2217170114046130768E-14
+
+};
+
+double FlatInvertiblePhasespace::generateIntermediates(vector<Energy>& K,
+ const double* r) const {
+
+ size_t n = K.size() + 1;
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ double u = bisect(r[i-2],n-1-i);
+ K[i-1] = sqrt(u*sqr(K[i-2]));
+ }
+
+ return flatWeights[n];
+
+}
+
+double FlatInvertiblePhasespace::invertIntermediates(const vector<Energy>& K,
+ double* r) const {
+
+ size_t n = K.size() + 1;
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ double u = sqr(K[i-1]/K[i-2]);
+ r[i-2] = (n+1-i)*pow(u,(double)(n-i)) - (n-i)*pow(u,(double)(n+1-i));
+ }
+
+ return flatWeights[n];
+
+}
+
+double FlatInvertiblePhasespace::generateIntermediates(vector<Energy>& M,
+ const vector<Energy>& m,
+ const double* r) const {
+
+ size_t n = M.size() + 1;
+
+ vector<Energy> K = M;
+ for ( size_t i = 1; i <= n; ++i )
+ K[0] -= m[i-1];
+
+ double w0 = generateIntermediates(K,r);
+
+ M = K;
+ for ( size_t i = 1; i <= n-1; ++i ) {
+ for ( size_t k = i; k <= n; ++k )
+ M[i-1] += m[k-1];
+ }
+
+ double weight = 8.*w0*rho(M[n-2],m[n-1],m[n-2]);
+
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ weight *=
+ (rho(M[i-2],M[i-1],m[i-2])/rho(K[i-2],K[i-1],ZERO)) * (M[i-1]/K[i-1]);
+ }
+
+ weight *= pow(K[0]/M[0],2.*n-4.);
+
+ return weight;
+
+}
+
+double FlatInvertiblePhasespace::invertIntermediates(const vector<Energy>& M,
+ const vector<Energy>& m,
+ double* r) const {
+
+ size_t n = M.size() + 1;
+
+ vector<Energy> K = M;
+ for ( size_t i = 1; i <= n-1; ++i ) {
+ for ( size_t k = i; k <= n; ++k )
+ K[i-1] -= m[k-1];
+ }
+
+ double w0 = invertIntermediates(K,r);
+
+ double weight = 8.*w0*rho(M[n-2],m[n-1],m[n-2]);
+
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ weight *=
+ (rho(M[i-2],M[i-1],m[i-2])/rho(K[i-2],K[i-1],ZERO)) * (M[i-1]/K[i-1]);
+ }
+
+ weight *= pow(K[0]/M[0],2.*n-4.);
+
+ return weight;
+
+}
+
+
+double FlatInvertiblePhasespace::generateKinematics(vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ const double* r) const {
+
+ vector<Energy> m;
+ for ( vector<Lorentz5Momentum>::const_iterator p =
+ P.begin(); p != P.end(); ++p )
+ m.push_back(p->mass());
+
+ size_t n = P.size();
+ vector<Energy> M(n-1);
+ M[0] = Ecm;
+
+ double weight = generateIntermediates(M,m,r);
+
+ M.push_back(m.back());
+
+ Lorentz5Momentum Q(M[0]);
+ Lorentz5Momentum nextQ;
+
+ for ( size_t i = 2; i <= n; ++i ) {
+
+ Energy q = 4.*M[i-2]*rho(M[i-2],M[i-1],m[i-2]);
+
+ double c = 2.*r[n-6+2*i]-1.;
+ double s = sqrt(1.-sqr(c));
+ double phi = 2.*Constants::pi*r[n-5+2*i];
+ double cphi = cos(phi);
+ double sphi = sqrt(1.-sqr(cphi));
+ if ( phi > Constants::pi )
+ sphi = -sphi;
+
+ P[i-2].setX(q*cphi*s);
+ P[i-2].setY(q*sphi*s);
+ P[i-2].setZ(q*c);
+ P[i-2].rescaleEnergy();
+ P[i-2].boost(Q.boostVector());
+ P[i-2].rescaleEnergy();
+
+ nextQ = Q - P[i-2];
+ nextQ.setMass(M[i-1]);
+ nextQ.rescaleEnergy();
+
+ Q = nextQ;
+
+ }
+
+ P.back() = Q;
+
+ return weight;
+
+}
+
+double FlatInvertiblePhasespace::invertKinematics(const vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ double* r) const {
+
+ vector<Energy> m;
+ for ( vector<Lorentz5Momentum>::const_iterator p =
+ P.begin(); p != P.end(); ++p )
+ m.push_back(p->mass());
+
+ size_t n = P.size();
+ vector<Energy> M(n-1);
+ M[0] = Ecm;
+
+ vector<Lorentz5Momentum> Q(n-1);
+ Q[0] = Lorentz5Momentum(M[0]);
+
+ for ( size_t i = 2; i <= n-1; ++i ) {
+ for ( size_t k = i; k <= n; ++k )
+ Q[i-1] += P[k-1];
+ M[i-1] = Q[i-1].m();
+ }
+
+ double weight = invertIntermediates(M,m,r);
+
+ for ( size_t i = 2; i <= n; ++i ) {
+ Lorentz5Momentum p = P[i-2];
+ p.boost(-Q[i-2].boostVector());
+ r[n-6+2*i] = (p.cosTheta()+1.)/2.;
+ double phi = p.phi();
+ if ( phi < 0. )
+ phi = 2.*Constants::pi + phi;
+ r[n-5+2*i] = phi/(2.*Constants::pi);
+ }
+
+ return weight;
+
+}
+
+
+double FlatInvertiblePhasespace::generateKinematics(const double* r,
+ vector<Lorentz5Momentum>& momenta) {
+
+ cPDVector::const_iterator pd = mePartonData().begin();
+ vector<Lorentz5Momentum>::iterator p = momenta.begin();
+ for ( ; pd != mePartonData().end(); ++pd, ++p )
+ p->setMass((**pd).mass());
+
+ Energy ECMMin = ZERO;
+ pd = mePartonData().begin() + 2;
+ for ( ; pd != mePartonData().end(); ++pd )
+ ECMMin += (**pd).mass();
+
+ double weight = 1.;
+
+ double tauMin = sqr(ECMMin)/lastXCombPtr()->lastS();
+ double tau = 1.;
+
+ if ( momenta.size() > 3 ) {
+ tau = tauMin + r[1]*(1.-tauMin);
+ } else {
+ tau = tauMin;
+ }
+
+ double ltau = log(tau)/2.;
+ weight *= -2.*ltau;
+
+ if ( momenta.size() > 3 ) {
+ weight *= 1.-tauMin;
+ } else {
+ weight *= 2.*Constants::pi;
+ }
+
+ double y = ltau - 2.*r[0]*ltau;
+ double x1 = sqrt(tau)*exp(y);
+ double x2 = sqrt(tau)*exp(-y);
+
+ ThreeVector<Energy> p1 =
+ x1*(lastXCombPtr()->lastParticles().first->momentum().vect());
+
+ ThreeVector<Energy> p2 =
+ x2*(lastXCombPtr()->lastParticles().second->momentum().vect());
+
+ momenta[0] = Lorentz5Momentum(momenta[0].mass(),p1);
+ momenta[1] = Lorentz5Momentum(momenta[1].mass(),p2);
+
+ lastXCombPtr()->lastX1X2(make_pair(x1,x2));
+ lastXCombPtr()->lastSHat((momenta[0]+momenta[1]).m2());
+
+ if ( momenta.size() == 3 ) {
+ ThreeVector<Energy> q = p1 + p2;
+ momenta[2] = Lorentz5Momentum(momenta[2].mass(),q);
+ return weight;
+ }
+
+ vector<Lorentz5Momentum> cmMomenta;
+ copy(momenta.begin()+2,momenta.end(),back_inserter(cmMomenta));
+
+ weight *= generateKinematics(cmMomenta,sqrt(lastXCombPtr()->lastSHat()),r+2);
+
+ Boost toLab = (momenta[0]+momenta[1]).boostVector();
+
+ vector<Lorentz5Momentum>::iterator pcm = cmMomenta.begin();
+ vector<Lorentz5Momentum>::iterator plab = momenta.begin() + 2;
+
+ for ( ; pcm != cmMomenta.end(); ++pcm, ++plab )
+ *plab = pcm->boost(toLab);
+
+ return weight/(x1*x2);
+
+}
+
+double FlatInvertiblePhasespace::invertKinematics(const vector<Lorentz5Momentum>& momenta,
+ double* r) const {
+
+ Energy ECMMin = ZERO;
+ cPDVector::const_iterator pd = mePartonData().begin() + 2;
+ for ( ; pd != mePartonData().end(); ++pd )
+ ECMMin += (**pd).mass();
+
+ double weight = 1.;
+
+ double tauMin = sqr(ECMMin)/lastXCombPtr()->lastS();
+ double tau = lastXCombPtr()->lastTau();
+
+ if ( momenta.size() > 3 ) {
+ r[1] = (tau - tauMin)/(1.-tauMin);
+ }
+
+ double ltau = log(tau)/2.;
+ weight *= -2.*ltau;
+
+ if ( momenta.size() > 3 ) {
+ weight *= 1.-tauMin;
+ } else {
+ weight *= 2.*Constants::pi;
+ }
+
+ double y = lastXCombPtr()->lastY();
+ r[0] = (ltau - y)/(2.*ltau);
+
+ if ( momenta.size() == 3 )
+ return weight;
+
+ vector<Lorentz5Momentum> cmMomenta;
+ copy(momenta.begin()+2,momenta.end(),back_inserter(cmMomenta));
+
+ Boost toCM = -(momenta[0]+momenta[1]).boostVector();
+
+ for ( vector<Lorentz5Momentum>::iterator p = cmMomenta.begin();
+ p != cmMomenta.end(); ++p )
+ p->boost(toCM);
+
+ weight *= invertKinematics(cmMomenta,sqrt(lastXCombPtr()->lastSHat()),r+2);
+
+ double x1 = lastXCombPtr()->lastX1();
+ double x2 = lastXCombPtr()->lastX2();
+
+ return weight/(x1*x2);
+
+}
+
+// If needed, insert default implementations of virtual function defined
+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
+
+
+void FlatInvertiblePhasespace::persistentOutput(PersistentOStream &) const {}
+
+void FlatInvertiblePhasespace::persistentInput(PersistentIStream &, int) {}
+
+
+// *** 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<FlatInvertiblePhasespace,MatchboxPhasespace>
+ describeHerwigFlatInvertiblePhasespace("Herwig::FlatInvertiblePhasespace", "HwMatchbox.so");
+
+void FlatInvertiblePhasespace::Init() {
+
+ static ClassDocumentation<FlatInvertiblePhasespace> documentation
+ ("FlatInvertiblePhasespace implements flat, invertible phase space generation.");
+
+}
+
diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
@@ -0,0 +1,209 @@
+// -*- C++ -*-
+//
+// FlatInvertiblePhasespace.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_FlatInvertiblePhasespace_H
+#define Herwig_FlatInvertiblePhasespace_H
+//
+// This is the declaration of the FlatInvertiblePhasespace class.
+//
+
+#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * \ingroup Matchbox
+ * \author Simon Platzer
+ *
+ * \brief FlatInvertiblePhasespace implements flat, invertible phase space generation.
+ *
+ */
+class FlatInvertiblePhasespace: public MatchboxPhasespace {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ FlatInvertiblePhasespace();
+
+ /**
+ * The destructor.
+ */
+ virtual ~FlatInvertiblePhasespace();
+ //@}
+
+public:
+
+ /**
+ * Prepare a phase space generator for the given xcomb object.
+ */
+ virtual void prepare(tStdXCombPtr, bool verbose = false);
+
+ /**
+ * Generate a phase space point and return its weight.
+ */
+ virtual double generateKinematics(const double*,
+ vector<Lorentz5Momentum>& momenta);
+
+ /**
+ * Return the number of random numbers required to produce a given
+ * multiplicity final state.
+ */
+ virtual int nDim(int nFinal) const {
+ if ( nFinal == 1 )
+ return 1;
+ return 3*nFinal - 2;
+ }
+
+ /**
+ * Return true, if this phasespace generator will generate incoming
+ * partons itself.
+ */
+ virtual bool haveX1X2() const { return true; }
+
+ /**
+ * Return true, if this phase space generator expects
+ * the incoming partons in their center-of-mass system
+ */
+ virtual bool wantCMS() const { return false; }
+
+public:
+
+ /**
+ * Return true, if this phase space generator is invertible
+ */
+ virtual bool isInvertible() const { return true; }
+
+ /**
+ * Invert the given phase space point to the random numbers which
+ * would have generated it.
+ */
+ virtual double invertKinematics(const vector<Lorentz5Momentum>&,
+ double*) const;
+
+private:
+
+ /**
+ * Solve v = (n+2) * u^(n+1) - (n+1) * u^(n+2) for u
+ */
+ double bisect(double v, double n,
+ double target = -16., double maxLevel = 80.) const;
+
+ /**
+ * Return rho
+ */
+ double rho(Energy M, Energy N, Energy m) const {
+ return sqrt((sqr(M)-sqr(N+m))*(sqr(M)-sqr(N-m)))/(8.*sqr(M));
+ }
+
+ /**
+ * Generate intermediate masses for a massless final state
+ */
+ double generateIntermediates(vector<Energy>& K,
+ const double* r) const;
+
+ /**
+ * Invert intermediate masses for a massless final state
+ */
+ double invertIntermediates(const vector<Energy>& K,
+ double* r) const;
+
+ /**
+ * Generate intermediate masses for a massive final state
+ */
+ double generateIntermediates(vector<Energy>& M,
+ const vector<Energy>& m,
+ const double* r) const;
+
+ /**
+ * Invert intermediate masses for a massive final state
+ */
+ double invertIntermediates(const vector<Energy>& M,
+ const vector<Energy>& m,
+ double* r) const;
+
+ /**
+ * Generate momenta in the CMS
+ */
+ double generateKinematics(vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ const double* r) const;
+
+ /**
+ * Invert momenta in the CMS
+ */
+ double invertKinematics(const vector<Lorentz5Momentum>& P,
+ Energy Ecm,
+ double* r) const;
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+
+// If needed, insert declarations of virtual function defined in the
+// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
+
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ FlatInvertiblePhasespace & operator=(const FlatInvertiblePhasespace &);
+
+};
+
+}
+
+#endif /* Herwig_FlatInvertiblePhasespace_H */
diff --git a/MatrixElement/Matchbox/Phasespace/Makefile.am b/MatrixElement/Matchbox/Phasespace/Makefile.am
--- a/MatrixElement/Matchbox/Phasespace/Makefile.am
+++ b/MatrixElement/Matchbox/Phasespace/Makefile.am
@@ -1,50 +1,52 @@
if WANT_DIPOLE
noinst_LTLIBRARIES = libHwMatchboxPhasespace.la
endif
libHwMatchboxPhasespace_la_SOURCES = \
FFLightInvertedTildeKinematics.cc \
FFLightInvertedTildeKinematics.h \
FFLightTildeKinematics.cc \
FFLightTildeKinematics.h \
FFMassiveInvertedTildeKinematics.cc \
FFMassiveInvertedTildeKinematics.h \
FFMassiveTildeKinematics.cc \
FFMassiveTildeKinematics.h \
FILightInvertedTildeKinematics.cc \
FILightInvertedTildeKinematics.h \
FILightTildeKinematics.cc \
FILightTildeKinematics.h \
FIMassiveInvertedTildeKinematics.cc \
FIMassiveInvertedTildeKinematics.h \
FIMassiveTildeKinematics.cc \
FIMassiveTildeKinematics.h \
IFLightInvertedTildeKinematics.cc \
IFLightInvertedTildeKinematics.h \
IFLightTildeKinematics.cc \
IFLightTildeKinematics.h \
IFMassiveInvertedTildeKinematics.cc \
IFMassiveInvertedTildeKinematics.h \
IFMassiveTildeKinematics.cc \
IFMassiveTildeKinematics.h \
IILightInvertedTildeKinematics.cc \
IILightInvertedTildeKinematics.h \
IILightTildeKinematics.cc \
IILightTildeKinematics.h \
InvertedTildeKinematics.cc \
InvertedTildeKinematics.h \
MatchboxPhasespace.cc \
MatchboxPhasespace.h \
MatchboxRambo.cc \
MatchboxRambo.h \
PhasespaceHelpers.cc \
PhasespaceHelpers.h \
RandomHelpers.h \
TildeKinematics.cc \
TildeKinematics.h \
TreePhasespace.cc \
TreePhasespace.h \
TreePhasespaceChannels.cc \
TreePhasespaceChannels.h \
MatchboxReference.h \
-MatchboxReference.cc
+MatchboxReference.cc \
+FlatInvertiblePhasespace.h \
+FlatInvertiblePhasespace.cc
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
--- a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
@@ -1,260 +1,276 @@
// -*- C++ -*-
//
// MatchboxPhasespace.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_MatchboxPhasespace_H
#define HERWIG_MatchboxPhasespace_H
//
// This is the declaration of the MatchboxPhasespace class.
//
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Wrap around a vector of random numbers to behave as a stream
* of those.
*/
struct StreamingRnd {
/**
* The random numbers
*/
const double* numbers;
/**
* The number of random numbers available.
*/
size_t nRnd;
/**
* Default constructor.
*/
StreamingRnd()
: numbers(0), nRnd(0) {}
/**
* Construct from random numbers.
*/
explicit StreamingRnd(const double* newNumbers,
size_t n)
: numbers(newNumbers), nRnd(n) {}
/**
* Return next random number
*/
inline double operator()() {
assert(numbers && nRnd > 0);
const double ret = numbers[0];
++numbers; --nRnd;
return ret;
}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxPhasespace defines an abstract interface to a phase
* space generator.
*
*/
class MatchboxPhasespace:
public HandlerBase, public LastXCombInfo<StandardXComb> {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxPhasespace();
/**
* The destructor.
*/
virtual ~MatchboxPhasespace();
//@}
public:
/**
* Prepare a phase space generator for the given xcomb object.
*/
virtual void prepare(tStdXCombPtr, bool verbose = false) = 0;
/**
* Generate a phase space point and return its weight.
*/
virtual double generateKinematics(const double*,
vector<Lorentz5Momentum>& momenta) = 0;
/**
* Return the number of random numbers required to produce a given
* multiplicity final state.
*/
virtual int nDim(int nFinal) const = 0;
/**
* Return true, if this phasespace generator will generate incoming
* partons itself.
*/
virtual bool haveX1X2() const { return false; }
/**
* Return true, if this phase space generator expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS() const { return true; }
/**
* Fill a diagram selector for the last phase space point.
*/
virtual Selector<MEBase::DiagramIndex> selectDiagrams(const MEBase::DiagramVector&) const;
/**
* Return the momentum and weight appropriate to the given timelike
* branch of the diagram.
*/
pair<double,Lorentz5Momentum> timeLikeWeight(const Tree2toNDiagram& diag,
int branch, double flatCut) const;
/**
* Return the weight appropriate to the given spacelike branch of
* the diagram.
*/
double spaceLikeWeight(const Tree2toNDiagram& diag,
const Lorentz5Momentum& incoming,
int branch, double flatCut) const;
/**
* Return the weight appropriate to the given diagram.
*/
double diagramWeight(const Tree2toNDiagram& diag) const {
assert( !diagramWeights.empty() );
return diagramWeights.find(diag.id())->second;
}
/**
* Fill the diagram weights.
*/
void fillDiagramWeights(double flatCut = 0.0);
/**
* Clone this phase space generator.
*/
Ptr<MatchboxPhasespace>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxPhasespace>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
virtual void cloneDependencies(const std::string& prefix = "");
public:
/**
+ * Return true, if this phase space generator is invertible
+ */
+ virtual bool isInvertible() const { return false; }
+
+ /**
+ * Invert the given phase space point to the random numbers which
+ * would have generated it.
+ */
+ virtual double invertKinematics(const vector<Lorentz5Momentum>&,
+ double*) const {
+ return 0.;
+ }
+
+public:
+
+ /**
* Limit phasespace generation to a given collinear or soft limit.
*/
void singularLimit(size_t i, size_t j) {
if ( i > j )
swap(i,j);
singularLimits.insert(make_pair(i,j));
}
/**
* Return the last matched singular limit.
*/
const pair<size_t,size_t>& lastSingularLimit() const {
assert(theLastSingularLimit != singularLimits.end());
return *theLastSingularLimit;
}
/**
* Return true, if constraints on phasespace generation have been met.
*/
bool matchConstraints(const vector<Lorentz5Momentum>& momenta);
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);
//@}
public:
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The diagram weights indexed by diagram id.
*/
map<int,double> diagramWeights;
/**
* If not empty, the entries here serve to limit phasespace
* generation to the corresponding collinear limits, or soft limits
* if both pair entries are equal.
*/
set<pair<size_t,size_t> > singularLimits;
/**
* The last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator theLastSingularLimit;
/**
* A cutoff below which a region is considered singular.
*/
Energy singularCutoff;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxPhasespace & operator=(const MatchboxPhasespace &);
};
}
#endif /* HERWIG_MatchboxPhasespace_H */
diff --git a/MatrixElement/Matchbox/Tests/HardProcessAnalysis.cc b/MatrixElement/Matchbox/Tests/HardProcessAnalysis.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Tests/HardProcessAnalysis.cc
@@ -0,0 +1,211 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the HardProcessAnalysis class.
+//
+
+#include "HardProcessAnalysis.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+#include "ThePEG/EventRecord/SubProcess.h"
+#include "ThePEG/EventRecord/SubProcessGroup.h"
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+HardProcessAnalysis::HardProcessAnalysis() {}
+
+HardProcessAnalysis::~HardProcessAnalysis() {}
+
+
+
+#ifndef LWH_AIAnalysisFactory_H
+#ifndef LWH
+#define LWH ThePEGLWH
+#endif
+#include "ThePEG/Analysis/LWH/AnalysisFactory.h"
+#endif
+
+HardProcessAnalysis::Histograms::Histograms(Energy ECM) {
+
+ size_t nbins = 100;
+
+ vector<double> logBins(nbins+1);
+ double logLow = 0.1;
+ double logUp = ECM/GeV;
+ double cLog = log10(logUp/logLow)/nbins;
+ for ( size_t k = 0; k < nbins+1; ++k )
+ logBins[k] = logLow*pow(10.0,cLog*k);
+
+ energy = new_ptr(Histogram(logBins));
+
+ logUp = ECM/GeV/4.;
+ cLog = log10(logUp/logLow)/nbins;
+ for ( size_t k = 0; k < nbins+1; ++k )
+ logBins[k] = logLow*pow(10.0,cLog*k);
+
+ transverse = new_ptr(Histogram(logBins));
+
+ cosTheta = new_ptr(Histogram(-1.,1.,nbins));
+
+ rapidity = new_ptr(Histogram(-7.,7.,nbins));
+
+ phi = new_ptr(Histogram(-Constants::pi,Constants::pi,nbins));
+
+}
+
+void HardProcessAnalysis::Histograms::fill(const Lorentz5Momentum& p, double weight) {
+
+ energy->addWeighted(p.t()/GeV,weight);
+ transverse->addWeighted(p.perp()/GeV,weight);
+ cosTheta->addWeighted(p.cosTheta(),weight);
+ rapidity->addWeighted(p.rapidity(),weight);
+ phi->addWeighted(p.phi(),weight);
+
+}
+
+void HardProcessAnalysis::Histograms::finalize(const string& subpro,
+ size_t legid) {
+
+ energy->normaliseToCrossSection();
+ transverse->normaliseToCrossSection();
+ cosTheta->normaliseToCrossSection();
+ rapidity->normaliseToCrossSection();
+ phi->normaliseToCrossSection();
+
+ ostringstream prefix;
+ prefix << subpro << "_" << legid;
+
+ string energyName = prefix.str() + "_energy.dat";
+ ofstream energyOut(energyName.c_str());
+ energy->rivetOutput(energyOut,"HardProcessAnalysis",prefix.str() + "_energy",
+ "Energy of " + prefix.str(),"$E$/GeV","${\\rm d}\\sigma/{\\rm d}E$/(nb/GeV)");
+
+ string transverseName = prefix.str() + "_transverse.dat";
+ ofstream transverseOut(transverseName.c_str());
+ transverse->rivetOutput(transverseOut,"HardProcessAnalysis",prefix.str() + "_transverse",
+ "Transverse momentum of " + prefix.str(),"$p_\\perp$/GeV","${\\rm d}\\sigma/{\\rm d}p_\\perp$/(nb/GeV)");
+
+ string costhetaName = prefix.str() + "_costheta.dat";
+ ofstream costhetaOut(costhetaName.c_str());
+ cosTheta->rivetOutput(costhetaOut,"HardProcessAnalysis",prefix.str() + "_costheta",
+ "Polar angle of " + prefix.str(),"$\\cos\\theta$","${\\rm d}\\sigma/{\\rm d}\\cos\\theta$/nb");
+
+ string rapidityName = prefix.str() + "_rapidity.dat";
+ ofstream rapidityOut(rapidityName.c_str());
+ rapidity->rivetOutput(rapidityOut,"HardProcessAnalysis",prefix.str() + "_rapidity",
+ "Rapidity of " + prefix.str(),"$y$","${\\rm d}\\sigma/{\\rm d}y$/nb");
+
+ string phiName = prefix.str() + "_phi.dat";
+ ofstream phiOut(phiName.c_str());
+ phi->rivetOutput(phiOut,"HardProcessAnalysis",prefix.str() + "_phi",
+ "Azimuthal angle of " + prefix.str(),"$\\phi$","${\\rm d}\\sigma/{\\rm d}\\phi$/nb");
+
+}
+
+struct SortIdEnergy {
+ inline bool operator()(PPtr a, PPtr b) const {
+ if ( a->id() < b->id() )
+ return true;
+ if ( a->momentum().t() > b->momentum().t() )
+ return true;
+ return false;
+ }
+};
+
+struct GetName {
+ inline string operator()(PPtr p) const {
+ return p->PDGName();
+ }
+};
+
+void HardProcessAnalysis::fill(PPair in, ParticleVector out, double weight) {
+ sort(out.begin(),out.end(),SortIdEnergy());
+ vector<string> proc;
+ proc.push_back(in.first->PDGName());
+ proc.push_back(in.second->PDGName());
+ std::transform(out.begin(),out.end(),
+ back_inserter(proc),GetName());
+ vector<Histograms>& data = histogramData[proc];
+ if ( data.empty() )
+ data.resize(out.size(),Histograms(generator()->maximumCMEnergy()));
+ ParticleVector::const_iterator p = out.begin();
+ vector<Histograms>::iterator h = data.begin();
+ for ( ; p != out.end(); ++p, ++h )
+ h->fill((**p).momentum(),weight);
+}
+
+void HardProcessAnalysis::analyze(tEventPtr event, long ieve, int loop, int state) {
+ AnalysisHandler::analyze(event, ieve, loop, state);
+ tSubProPtr sub = event->primarySubProcess();
+ Ptr<SubProcessGroup>::tptr grp =
+ dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);
+ fill(sub->incoming(),sub->outgoing(),event->weight()*sub->groupWeight());
+ if ( grp ) {
+ for ( SubProcessVector::const_iterator s = grp->dependent().begin();
+ s != grp->dependent().end(); ++s ) {
+ fill((**s).incoming(),(**s).outgoing(),event->weight()*(**s).groupWeight());
+ }
+ }
+}
+
+void HardProcessAnalysis::dofinish() {
+ AnalysisHandler::dofinish();
+ for ( map<vector<string>,vector<Histograms> >::iterator h =
+ histogramData.begin(); h != histogramData.end(); ++h ) {
+ string subpro;
+ for ( vector<string>::const_iterator p = h->first.begin();
+ p != h->first.end(); ++p ) {
+ subpro += *p + (p != --(h->first.end()) ? "_" : "");
+ }
+ for ( size_t k = 0; k < h->second.size(); ++k )
+ h->second[k].finalize(subpro,k+2);
+ }
+}
+
+void HardProcessAnalysis::doinitrun() {
+ AnalysisHandler::doinitrun();
+ // *** ATTENTION *** histogramFactory().registerClient(this); // Initialize histograms.
+ // *** ATTENTION *** histogramFactory().mkdirs("/SomeDir"); // Put histograms in specal directory.
+}
+
+
+IBPtr HardProcessAnalysis::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr HardProcessAnalysis::fullclone() const {
+ return new_ptr(*this);
+}
+
+
+// If needed, insert default implementations of virtual function defined
+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
+
+
+void HardProcessAnalysis::persistentOutput(PersistentOStream &) const {}
+
+void HardProcessAnalysis::persistentInput(PersistentIStream &, int) {}
+
+
+// *** 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<HardProcessAnalysis,AnalysisHandler>
+ describeHerwigHardProcessAnalysis("Herwig::HardProcessAnalysis", "HardProcessAnalysis.so");
+
+void HardProcessAnalysis::Init() {
+
+ static ClassDocumentation<HardProcessAnalysis> documentation
+ ("There is no documentation for the HardProcessAnalysis class");
+
+}
+
diff --git a/MatrixElement/Matchbox/Tests/HardProcessAnalysis.h b/MatrixElement/Matchbox/Tests/HardProcessAnalysis.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Tests/HardProcessAnalysis.h
@@ -0,0 +1,201 @@
+// -*- C++ -*-
+#ifndef Herwig_HardProcessAnalysis_H
+#define Herwig_HardProcessAnalysis_H
+//
+// This is the declaration of the HardProcessAnalysis class.
+//
+
+#include "ThePEG/Handlers/AnalysisHandler.h"
+#include "Herwig++/Utilities/Histogram.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the HardProcessAnalysis class.
+ *
+ * @see \ref HardProcessAnalysisInterfaces "The interfaces"
+ * defined for HardProcessAnalysis.
+ */
+class HardProcessAnalysis: public AnalysisHandler {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ HardProcessAnalysis();
+
+ /**
+ * The destructor.
+ */
+ virtual ~HardProcessAnalysis();
+ //@}
+
+public:
+
+ /** @name Virtual functions required by the AnalysisHandler class. */
+ //@{
+ /**
+ * Analyze a given Event. Note that a fully generated event
+ * may be presented several times, if it has been manipulated in
+ * between. The default version of this function will call transform
+ * to make a lorentz transformation of the whole event, then extract
+ * all final state particles and call analyze(tPVector) of this
+ * analysis object and those of all associated analysis objects. The
+ * default version will not, however, do anything on events which
+ * have not been fully generated, or have been manipulated in any
+ * way.
+ * @param event pointer to the Event to be analyzed.
+ * @param ieve the event number.
+ * @param loop the number of times this event has been presented.
+ * If negative the event is now fully generated.
+ * @param state a number different from zero if the event has been
+ * manipulated in some way since it was last presented.
+ */
+ virtual void analyze(tEventPtr event, long ieve, int loop, int state);
+
+ //@}
+
+protected:
+
+ /**
+ * Initialize this object. Called in the run phase just before
+ * a run begins.
+ */
+ virtual void doinitrun();
+
+ /**
+ * Finalize this object. Called in the run phase just after a
+ * run has ended. Used eg. to write out statistics.
+ */
+ virtual void dofinish();
+
+
+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;
+ //@}
+
+
+// If needed, insert declarations of virtual function defined in the
+// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
+
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ HardProcessAnalysis & operator=(const HardProcessAnalysis &);
+
+ /**
+ * Differential information per outgoing parton
+ */
+ struct Histograms {
+
+ /**
+ * The constructor
+ */
+ Histograms() {}
+
+ /**
+ * The constructor
+ */
+ explicit Histograms(Energy ECM);
+
+ /**
+ * Analyse given momentum
+ */
+ void fill(const Lorentz5Momentum& p, double weight);
+
+ /**
+ * Finalize given process id and cross section.
+ */
+ void finalize(const string& subpro,
+ size_t legid);
+
+ /**
+ * Energy spectrum
+ */
+ HistogramPtr energy;
+
+ /**
+ * Pt spectrum
+ */
+ HistogramPtr transverse;
+
+ /**
+ * Polar angle distribution
+ */
+ HistogramPtr cosTheta;
+
+ /**
+ * Rapidity distribution
+ */
+ HistogramPtr rapidity;
+
+ /**
+ * Azimuthal angle distribution
+ */
+ HistogramPtr phi;
+
+ };
+
+ /**
+ * Histograms per subprocess
+ */
+ map<vector<string>,vector<Histograms> > histogramData;
+
+ /**
+ * Analyze a given final state
+ */
+ void fill(PPair, ParticleVector, double);
+
+};
+
+}
+
+#endif /* Herwig_HardProcessAnalysis_H */
diff --git a/MatrixElement/Matchbox/Tests/Makefile.am b/MatrixElement/Matchbox/Tests/Makefile.am
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Tests/Makefile.am
@@ -0,0 +1,6 @@
+if WANT_DIPOLE
+noinst_LTLIBRARIES = libHwMatchboxTests.la
+endif
+
+libHwMatchboxTests_la_SOURCES = \
+HardProcessAnalysis.cc HardProcessAnalysis.h
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1,190 +1,191 @@
dnl Process this file with autoconf to produce a configure script.
AC_PREREQ([2.59])
AC_INIT([Herwig++],[devel],[herwig@projects.hepforge.org],[Herwig++])
AC_CONFIG_SRCDIR([Utilities/HerwigStrategy.cc])
AC_CONFIG_AUX_DIR([Config])
AC_CONFIG_MACRO_DIR([m4])
AC_CONFIG_HEADERS([Config/config.h])
dnl AC_PRESERVE_HELP_ORDER
AC_CANONICAL_HOST
case "${host}" in
*-darwin[[0156]].*)
AC_MSG_ERROR([Herwig++ requires OS X 10.3 or later])
;;
*-darwin7.*)
if test "x$MACOSX_DEPLOYMENT_TARGET" != "x10.3"; then
AC_MSG_ERROR(
[Please add 'MACOSX_DEPLOYMENT_TARGET=10.3' to the configure line.])
fi
;;
esac
dnl === disable debug symbols by default =====
if test "x$CXXFLAGS" = "x"; then
CXXFLAGS=-O3
fi
if test "x$CFLAGS" = "x"; then
CFLAGS=-O3
fi
dnl Looptools manual requires optimization off
if test "x$FCFLAGS" = "x"; then
FCFLAGS=-O0
fi
dnl ==========================================
AC_LANG([C++])
AM_INIT_AUTOMAKE([1.9 gnu dist-bzip2 -Wall])
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
m4_ifdef([AM_PROG_AR], [AM_PROG_AR])
dnl Checks for programs.
AC_PROG_CXX([g++])
AC_PROG_INSTALL
AC_PROG_MAKE_SET
AC_PROG_LN_S
dnl modified search order
AC_PROG_FC([gfortran g95 g77])
dnl xlf95 f95 fort ifort ifc efc pgf95 lf95 ftn xlf90 f90 pgf90 pghpf epcf90 xlf f77 frt pgf77 cf77 fort77 fl32 af77])
AC_LANG_PUSH([Fortran])
AC_MSG_CHECKING([if the Fortran compiler ($FC) works])
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([A Fortran compiler is required to build Herwig++.])
]
)
AC_LANG_POP([Fortran])
LT_PREREQ([2.2])
LT_INIT([disable-static dlopen pic-only])
dnl ####################################
dnl ####################################
dnl for Doc/fixinterfaces.pl
AC_PATH_PROG(PERL, perl)
HERWIG_CHECK_GSL
HERWIG_CHECK_THEPEG
HERWIG_CHECK_BOOST
HERWIG_COMPILERFLAGS
HERWIG_LOOPTOOLS
HERWIG_PDF_PATH
FASTJET_CHECK_FASTJET
HERWIG_VERSIONSTRING
HERWIG_CHECK_ABS_BUG
HERWIG_ENABLE_MODELS
HERWIG_ENABLE_DIPOLE
SHARED_FLAG=-shared
AM_CONDITIONAL(NEED_APPLE_FIXES,
[test "xx${host/darwin/foundit}xx" != "xx${host}xx"])
if test "xx${host/darwin/foundit}xx" != "xx${host}xx"; then
APPLE_DSO_FLAGS=-Wl,-undefined,dynamic_lookup
SHARED_FLAG=-bundle
fi
AC_SUBST([APPLE_DSO_FLAGS])
AC_SUBST([SHARED_FLAG])
AC_CONFIG_FILES([UnderlyingEvent/Makefile
Models/Makefile
Models/StandardModel/Makefile
Models/RSModel/Makefile
Models/General/Makefile
Models/Susy/Makefile
Models/Susy/NMSSM/Makefile
Models/UED/Makefile
Models/LH/Makefile
Models/LHTP/Makefile
Models/Transplanckian/Makefile
Models/Leptoquarks/Makefile
Models/Zprime/Makefile
Models/TTbAsymm/Makefile
Models/ADD/Makefile
Models/Sextet/Makefile
Decay/Makefile
Decay/FormFactors/Makefile
Decay/Tau/Makefile
Decay/Baryon/Makefile
Decay/VectorMeson/Makefile
Decay/Perturbative/Makefile
Decay/ScalarMeson/Makefile
Decay/TensorMeson/Makefile
Decay/WeakCurrents/Makefile
Decay/Partonic/Makefile
Decay/General/Makefile
Decay/Radiation/Makefile
Doc/refman.conf
Doc/refman.h
PDT/Makefile
PDF/Makefile
MatrixElement/Makefile
MatrixElement/General/Makefile
MatrixElement/Lepton/Makefile
MatrixElement/Hadron/Makefile
MatrixElement/DIS/Makefile
MatrixElement/Powheg/Makefile
MatrixElement/Gamma/Makefile
MatrixElement/Matchbox/Makefile
MatrixElement/Matchbox/Base/Makefile
MatrixElement/Matchbox/Utility/Makefile
MatrixElement/Matchbox/Phasespace/Makefile
MatrixElement/Matchbox/Dipoles/Makefile
MatrixElement/Matchbox/InsertionOperators/Makefile
MatrixElement/Matchbox/Powheg/Makefile
MatrixElement/Matchbox/Matching/Makefile
MatrixElement/Matchbox/Builtin/Makefile
MatrixElement/Matchbox/Builtin/Processes/Makefile
MatrixElement/Matchbox/Builtin/Amplitudes/Makefile
+ MatrixElement/Matchbox/Tests/Makefile
Exsample2/Makefile
Shower/SplittingFunctions/Makefile
Shower/Default/Makefile
Shower/Base/Makefile
Shower/Makefile
DipoleShower/Makefile
DipoleShower/Base/Makefile
DipoleShower/Kernels/Makefile
DipoleShower/Kinematics/Makefile
DipoleShower/Utility/Makefile
DipoleShower/AlphaS/Makefile
Utilities/Makefile
Hadronization/Makefile
lib/Makefile
include/Makefile
src/Makefile
src/defaults/Makefile
src/herwig-config
Doc/Makefile
Doc/HerwigDefaults.in
Looptools/Makefile
Analysis/Makefile
src/Makefile-UserModules
src/defaults/Analysis.in
Contrib/Makefile
Contrib/make_makefiles.sh
Tests/Makefile
Makefile])
AC_CONFIG_FILES([Doc/fixinterfaces.pl],[chmod +x Doc/fixinterfaces.pl])
HERWIG_OVERVIEW
AC_CONFIG_COMMANDS([summary],[cat config.herwig])
AC_OUTPUT

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 2:21 PM (12 h, 5 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4013584
Default Alt Text
(287 KB)

Event Timeline