Page MenuHomeHEPForge
Files F8723834

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
@@ -1,1209 +1,1227 @@
// -*- 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) 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 {
+double MatchboxMEBase::pdf1(Energy2 fscale, double xEx) const {
assert(lastXCombPtr()->partonBins().first->pdf());
+ if ( xEx < 1. && lastX1() >= xEx ) {
+ return
+ ( ( 1. - lastX1() ) / ( 1. - xEx ) ) *
+ lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
+ lastPartons().first->dataPtr(),
+ fscale == ZERO ? lastScale() : fscale,
+ xEx)/xEx;
+ }
+
return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
lastPartons().first->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX1())/lastX1();
}
-double MatchboxMEBase::pdf2(Energy2 fscale) const {
+double MatchboxMEBase::pdf2(Energy2 fscale, double xEx) const {
assert(lastXCombPtr()->partonBins().second->pdf());
+ if ( xEx < 1. && lastX2() >= xEx ) {
+ return
+ ( ( 1. - lastX2() ) / ( 1. - xEx ) ) *
+ lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
+ lastPartons().second->dataPtr(),
+ fscale == ZERO ? lastScale() : fscale,
+ xEx)/xEx;
+ }
+
return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
lastPartons().second->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX2())/lastX2();
}
double MatchboxMEBase::me2() const {
if ( matchboxAmplitude() ) {
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()* lastMEPDFWeight() * xme2;
if ( oneLoop() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * vme2;
if ( !onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).setXComb(lastXCombPtr());
res += (**v).dSigHatDR();
}
}
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
lastMECrossSection(res);
return lastMECrossSection();
}
double MatchboxMEBase::oneLoopInterference() const {
if ( matchboxAmplitude() ) {
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(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,949 +1,951 @@
// -*- 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 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;
+ double pdf1(Energy2 factorizationScale = ZERO,
+ double xEx = 1.) const;
/**
* Supply the PDF weight for the second incoming parton.
*/
- double pdf2(Energy2 factorizationScale = ZERO) const;
+ double pdf2(Energy2 factorizationScale = ZERO,
+ double xEx = 1.) const;
//@}
/** @name Amplitude information and matrix element evaluation */
//@{
/**
* Return the amplitude.
*/
Ptr<MatchboxAmplitude>::tptr matchboxAmplitude() const { return theAmplitude; }
/**
* Set the amplitude.
*/
void matchboxAmplitude(Ptr<MatchboxAmplitude>::ptr amp) { theAmplitude = amp; }
/**
* Return the matrix element for the kinematical configuation
* previously provided by the last call to setKinematics(), suitably
* scaled by sHat() to give a dimension-less number.
*/
virtual double me2() const;
/**
* Return the symmetry factor for identical final state particles.
*/
virtual double finalStateSymmetry() const;
/**
* Return the normalizing factor for the matrix element averaged
* over quantum numbers and including running couplings.
*/
double me2Norm(unsigned int addAlphaS = 0) const;
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR() const;
//@}
/** @name One-loop corrections */
//@{
/**
* Return the one-loop/tree interference.
*/
virtual double oneLoopInterference() const;
/**
* Return true, if this matrix element is capable of calculating
* one-loop (QCD) corrections.
*/
virtual bool haveOneLoop() const;
/**
* Return true, if this matrix element only provides
* one-loop (QCD) corrections.
*/
virtual bool onlyOneLoop() const;
/**
* Return true, if one loop corrections have been calculated in
* dimensional reduction. Otherwise conventional dimensional
* regularization is assumed. Note that renormalization is always
* assumed to be MSbar.
*/
virtual bool isDR() const;
/**
* Return true, if one loop corrections are given in the conventions
* of the integrated dipoles.
*/
virtual bool isCS() const;
/**
* Return true, if one loop corrections are given in the conventions
* of BDK.
*/
virtual bool isBDK() const;
/**
* Return true, if one loop corrections are given in the conventions
* of everything expanded.
*/
virtual bool isExpanded() const;
/**
* Return the value of the dimensional regularization
* parameter. Note that renormalization scale dependence is fully
* restored in DipoleIOperator.
*/
virtual Energy2 mu2() const;
/**
* If defined, return the coefficient of the pole in epsilon^2
*/
virtual double oneLoopDoublePole() const;
/**
* If defined, return the coefficient of the pole in epsilon
*/
virtual double oneLoopSinglePole() const;
/**
* Return true, if cancellationn of epsilon poles should be checked.
*/
bool checkPoles() const { 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/SubtractedME.cc b/MatrixElement/Matchbox/Base/SubtractedME.cc
--- a/MatrixElement/Matchbox/Base/SubtractedME.cc
+++ b/MatrixElement/Matchbox/Base/SubtractedME.cc
@@ -1,669 +1,692 @@
// -*- 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) 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(const vector<tStdXCombPtr>& dep) {
+
+ if ( inclusive() && !lastXComb().lastProjector() )
+ return 1.;
+
+ if ( virtualShowerSubtraction() && !lastXComb().lastProjector() ) {
+ assert(!showerApproximation()->belowCutoff());
+ return 0.;
+ }
+
if ( realShowerSubtraction() ) {
assert(showerApproximation());
bool below = showerApproximation()->belowCutoff();
if ( below )
return 0.;
return 1.;
}
- if ( virtualShowerSubtraction() ) {
- assert(showerApproximation());
- bool above = !showerApproximation()->belowCutoff();
- if ( above )
- return 0.;
+
+ if ( virtualShowerSubtraction() || inclusive() ) {
+ if ( virtualShowerSubtraction() ) {
+ assert(showerApproximation());
+ bool above = !showerApproximation()->belowCutoff();
+ if ( above )
+ return 0.;
+ }
+ double sum = 0.;
+ size_t n = 0;
+ for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
+ if ( (**d).lastCrossSection() != ZERO ) {
+ sum += (**d).lastME2();
+ ++n;
+ }
+ }
+ return
+ n * lastXComb().lastProjector()->lastME2() / sum;
}
- double sum = 0.;
- for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d )
- sum += (**d).lastME2();
- return
- dep.size() * lastXComb().lastProjector()->lastME2() / sum;
+
+ return 1.;
+
}
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 ( inclusive() && !lastXComb().lastProjector() )
+ return 0.;
+
+ if ( virtualShowerSubtraction() && !lastXComb().lastProjector() ) {
+ assert(!showerApproximation()->belowCutoff());
+ return 0.;
}
- if ( xc != lastXComb().lastProjector() )
- return 0.;
- 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;
+
+ if ( virtualShowerSubtraction() || inclusive() ) {
+ if ( xc != lastXComb().lastProjector() )
+ return 0.;
+ size_t n = 0;
+ for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
+ if ( (**d).lastCrossSection() != ZERO ) {
+ ++n;
+ }
+ }
+ return n;
}
- return w * dep.size();
+
+ return 1.;
+
}
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/Dipoles/SubtractionDipole.cc b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
@@ -1,1041 +1,1043 @@
// -*- 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());
lastXCombPtr()->didGenerateKinematics();
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();
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;
}
- lastME2(me2());
lastMECrossSection(-showerApproximation()->dSigHatDR());
return lastMECrossSection();
}
double xme2 = 0.0;
- if ( !showerKernel() )
- xme2 = me2();
- else
- xme2 = me2Avg(-underlyingBornME()->me2());
+ if ( lastME2() == 0.0 ) {
+ if ( !showerKernel() )
+ xme2 = me2();
+ else
+ xme2 = me2Avg(-underlyingBornME()->me2());
+ } else {
+ xme2 = lastME2();
+ }
if ( xme2 == 0.0 ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
lastME2(xme2);
CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
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 ( 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);
+ res -= shower;
}
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/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc
--- a/MatrixElement/Matchbox/MatchboxFactory.cc
+++ b/MatrixElement/Matchbox/MatchboxFactory.cc
@@ -1,1164 +1,1177 @@
// -*- C++ -*-
//
// MatchboxFactory.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxFactory class.
//
#include "MatchboxFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SU2Helper.h"
#include <boost/progress.hpp>
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
using std::ostream_iterator;
MatchboxFactory::MatchboxFactory()
: SubProcessHandler(), theNLight(0),
theOrderInAlphaS(0), theOrderInAlphaEW(0),
theBornContributions(true), theVirtualContributions(true),
theRealContributions(true), theIndependentVirtuals(false),
theSubProcessGroups(false), theInclusive(false),
theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0),
theFixedCouplings(false), theFixedQEDCouplings(false), theVetoScales(false),
theVerbose(false), theInitVerbose(false), theSubtractionData(""), theCheckPoles(false) {}
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() ) {
+ generator()->log() << "preparing Born "
+ << (virtualContributions() ? "and virtual" : "")
+ << " matrix elements." << flush;
+ }
+
if ( (bornContributions() && !virtualContributions()) || independentVirtuals() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
if ( (**born).onlyOneLoop() )
continue;
Ptr<MatchboxMEBase>::ptr bornme = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( independentVirtuals() )
pname += ".Born";
if ( ! (generator()->preinitRegister(bornme,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
bornme->cloneDependencies();
MEs().push_back(bornme);
}
}
if ( virtualContributions() ) {
- 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<MatchboxMEBase>::ptr nlo = (**born).cloneMe();
- string pname = fullName() + "/" + (**born).name();\
+ string pname = fullName() + "/" + (**born).name();
if ( !independentVirtuals() )
pname += ".BornVirtual";
else
pname += ".Virtual";
if ( ! (generator()->preinitRegister(nlo,pname) ) )
throw InitException() << "NLO ME " << pname << " already existing.";
nlo->virtuals().clear();
if ( !nlo->onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
if ( nlo->virtuals().empty() )
throw InitException() << "No insertion operators have been found for "
<< (**born).name() << "\n";
if ( checkPoles() ) {
if ( !virtualsAreExpanded ) {
throw InitException() << "Cannot check epsilon poles if virtuals are not in `expanded' convention.\n";
}
nlo->doCheckPoles();
}
}
if ( !bornContributions() || independentVirtuals() ) {
nlo->doOneLoopNoBorn();
} else {
nlo->doOneLoop();
}
nlo->cloneDependencies();
bornVirtualMEs().push_back(nlo);
MEs().push_back(nlo);
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
theSplittingDipoles.clear();
set<cPDVector> bornProcs;
if ( showerApproximation() )
if ( showerApproximation()->needsSplittingGenerator() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born )
for ( MEBase::DiagramVector::const_iterator d = (**born).diagrams().begin();
d != (**born).diagrams().end(); ++d )
bornProcs.insert((**d).partons());
}
if ( realContributions() ) {
generator()->log() << "preparing real emission matrix elements." << 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() + ".Real";
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);
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;
+ if ( showerApproximation()->needsTildeXCombs() ) {
+ channel.tildeXCombs.clear();
+ assert(!channel.dipole->partnerDipoles().empty());
+ for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator p =
+ channel.dipole->partnerDipoles().begin();
+ p != channel.dipole->partnerDipoles().end(); ++p ) {
+ StdXCombPtr txc = channel.dipole->makeBornXComb(channel.realXComb);
+ if ( txc )
+ channel.tildeXCombs.push_back(txc);
+ }
+ }
res.push_back(channel);
}
}
return res;
}
void MatchboxFactory::print(ostream& os) const {
os << "--- MatchboxFactory setup -----------------------------------------------------------\n";
if ( !amplitudes().empty() ) {
os << " generated Born matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = bornMEs().begin();
m != bornMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
os << " generated real emission matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = realEmissionMEs().begin();
m != realEmissionMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
}
os << " generated Born+virtual matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator bv
= bornVirtualMEs().begin(); bv != bornVirtualMEs().end(); ++bv ) {
(**bv).print(os);
}
os << " generated subtracted matrix elements:\n";
for ( vector<Ptr<SubtractedME>::ptr>::const_iterator sub
= subtractedMEs().begin(); sub != subtractedMEs().end(); ++sub ) {
os << " '" << (**sub).name() << "'\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxFactory::doinit() {
setup();
if ( initVerbose() )
print(Repository::clog());
SubProcessHandler::doinit();
}
void MatchboxFactory::persistentOutput(PersistentOStream & os) const {
os << theDiagramGenerator << theProcessData
<< theNLight << theOrderInAlphaS << theOrderInAlphaEW
<< theBornContributions << theVirtualContributions
<< theRealContributions << theIndependentVirtuals << theSubProcessGroups << theInclusive
<< thePhasespace << theScaleChoice
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theFixedCouplings << theFixedQEDCouplings << theVetoScales
<< theAmplitudes << theCache
<< theBornMEs << theVirtuals << theRealEmissionMEs
<< theBornVirtualMEs << theSubtractedMEs << theFiniteRealMEs
<< theVerbose << theInitVerbose << theSubtractionData << theCheckPoles
<< theParticleGroups << process << realEmissionProcess
<< theShowerApproximation << theSplittingDipoles;
}
void MatchboxFactory::persistentInput(PersistentIStream & is, int) {
is >> theDiagramGenerator >> theProcessData
>> theNLight >> theOrderInAlphaS >> theOrderInAlphaEW
>> theBornContributions >> theVirtualContributions
>> theRealContributions >> theIndependentVirtuals >> theSubProcessGroups >> theInclusive
>> thePhasespace >> theScaleChoice
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theFixedCouplings >> theFixedQEDCouplings >> theVetoScales
>> theAmplitudes >> theCache
>> theBornMEs >> theVirtuals >> theRealEmissionMEs
>> theBornVirtualMEs >> theSubtractedMEs >> theFiniteRealMEs
>> 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> interfaceIndependentVirtuals
("IndependentVirtuals",
"Switch on or off virtual contributions as separate subprocesses.",
&MatchboxFactory::theIndependentVirtuals, true, false, false);
static SwitchOption interfaceIndependentVirtualsOn
(interfaceIndependentVirtuals,
"On",
"Switch on virtual contributions as separate subprocesses.",
true);
static SwitchOption interfaceIndependentVirtualsOff
(interfaceIndependentVirtuals,
"Off",
"Switch off virtual contributions as separate subprocesses.",
false);
static Switch<MatchboxFactory,bool> interfaceSubProcessGroups
("SubProcessGroups",
"Switch on or off production of sub-process groups.",
&MatchboxFactory::theSubProcessGroups, false, false, false);
static SwitchOption interfaceSubProcessGroupsOn
(interfaceSubProcessGroups,
"On",
"On",
true);
static SwitchOption interfaceSubProcessGroupsOff
(interfaceSubProcessGroups,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceInclusive
("Inclusive",
"Switch on or off production of inclusive cross section.",
&MatchboxFactory::theInclusive, false, false, false);
static SwitchOption interfaceInclusiveOn
(interfaceInclusive,
"On",
"On",
true);
static SwitchOption interfaceInclusiveOff
(interfaceInclusive,
"Off",
"Off",
false);
static Reference<MatchboxFactory,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"Set the phasespace generator.",
&MatchboxFactory::thePhasespace, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxScaleChoice> interfaceScaleChoice
("ScaleChoice",
"Set the scale choice object.",
&MatchboxFactory::theScaleChoice, false, false, true, true, false);
static Parameter<MatchboxFactory,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&MatchboxFactory::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&MatchboxFactory::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceFixedCouplings
("FixedCouplings",
"Switch on or off fixed couplings.",
&MatchboxFactory::theFixedCouplings, true, false, false);
static SwitchOption interfaceFixedCouplingsOn
(interfaceFixedCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedCouplingsOff
(interfaceFixedCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceFixedQEDCouplings
("FixedQEDCouplings",
"Switch on or off fixed QED couplings.",
&MatchboxFactory::theFixedQEDCouplings, true, false, false);
static SwitchOption interfaceFixedQEDCouplingsOn
(interfaceFixedQEDCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedQEDCouplingsOff
(interfaceFixedQEDCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceVetoScales
("VetoScales",
"Switch on or setting veto scales.",
&MatchboxFactory::theVetoScales, false, false, false);
static SwitchOption interfaceVetoScalesOn
(interfaceVetoScales,
"On",
"On",
true);
static SwitchOption interfaceVetoScalesOff
(interfaceVetoScales,
"Off",
"Off",
false);
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceAmplitudes
("Amplitudes",
"The amplitude objects.",
&MatchboxFactory::theAmplitudes, -1, false, false, true, true, false);
static 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,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,776 +1,781 @@
// -*- 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/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 virtual contributions should be treated as independent subprocesses
*/
bool independentVirtuals() const { return theIndependentVirtuals; }
/**
* Switch on/off virtual contributions should be treated as independent subprocesses
*/
void setIndependentVirtuals(bool on = true) { theIndependentVirtuals = on; }
/**
* Return true, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool subProcessGroups() const { return theSubProcessGroups; }
/**
* Switch on or off producing subprocess groups.
*/
void setSubProcessGroups(bool on = true) { theSubProcessGroups = on; }
/**
* Return true, if 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<MatchboxMEBase>::ptr>& bornVirtualMEs() const { return theBornVirtualMEs; }
/**
* Access the produced NLO matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr>& bornVirtualMEs() { return theBornVirtualMEs; }
/**
* Return the real emission matrix elements to be considered
*/
const vector<Ptr<MatchboxMEBase>::ptr>& realEmissionMEs() const { return theRealEmissionMEs; }
/**
* Access the real emission matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr>& realEmissionMEs() { return theRealEmissionMEs; }
/**
* Return the produced subtracted matrix elements
*/
const vector<Ptr<SubtractedME>::ptr>& subtractedMEs() const { return theSubtractedMEs; }
/**
* Access the produced subtracted matrix elements
*/
vector<Ptr<SubtractedME>::ptr>& subtractedMEs() { return theSubtractedMEs; }
/**
* Return the produced finite real emission matrix elements
*/
const vector<Ptr<MatchboxMEBase>::ptr>& finiteRealMEs() const { return theFiniteRealMEs; }
/**
* Access the produced finite real emission elements
*/
vector<Ptr<MatchboxMEBase>::ptr>& finiteRealMEs() { return theFiniteRealMEs; }
/**
* Return the map of Born processes to splitting dipoles
*/
const map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >& splittingDipoles() const {
return theSplittingDipoles;
}
/**
* Identify a splitting channel
*/
struct SplittingChannel {
/**
* The Born XComb
*/
StdXCombPtr bornXComb;
/**
* The real XComb
*/
StdXCombPtr realXComb;
/**
+ * The set of tilde XCombs to consider for the real xcomb
+ */
+ vector<StdXCombPtr> tildeXCombs;
+
+ /**
* The dipole in charge of the splitting
*/
Ptr<SubtractionDipole>::ptr dipole;
};
/**
* 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 virtual contributions should be treated as independent subprocesses
*/
bool theIndependentVirtuals;
/**
* True, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool theSubProcessGroups;
/**
* True, if the integral over the unresolved emission should be
* calculated.
*/
bool theInclusive;
/**
* The phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The scale choice object
*/
Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
/**
* Use non-running couplings.
*/
bool theFixedCouplings;
/**
* Use non-running couplings.
*/
bool theFixedQEDCouplings;
/**
* True, if veto scales should be set
* for the real emission
*/
bool theVetoScales;
/**
* The amplitudes to be considered
*/
vector<Ptr<MatchboxAmplitude>::ptr> theAmplitudes;
/**
* The 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<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/DipoleMatching.cc b/MatrixElement/Matchbox/Matching/DipoleMatching.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Matching/DipoleMatching.cc
@@ -0,0 +1,88 @@
+// -*- C++ -*-
+//
+// DipoleMatching.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 DipoleMatching class.
+//
+
+#include "DipoleMatching.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
+
+using namespace Herwig;
+
+DipoleMatching::DipoleMatching() {}
+
+DipoleMatching::~DipoleMatching() {}
+
+IBPtr DipoleMatching::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr DipoleMatching::fullclone() const {
+ return new_ptr(*this);
+}
+
+CrossSection DipoleMatching::dSigHatDR() const {
+ double pdfFactor = 1.;
+ if ( showerScalesInSubtraction() ) {
+ double bornPDF = bornPDFWeight(showerScalesInSubtraction());
+ double bornPDFHard = bornPDFWeight(false);
+ pdfFactor = bornPDFHard / bornPDF;
+ }
+ return
+ sqr(hbarc) *
+ realXComb()->jacobian() *
+ realPDFWeight(showerScalesInSubtraction()) * pdfFactor *
+ couplingWeight(showerScalesInSubtraction()) *
+ dipole()->me2() /
+ (2. * realXComb()->lastSHat());
+}
+
+double DipoleMatching::me2() const {
+ throw Exception() << "Not intented to use. Disable the ShowerApproximationGenerator."
+ << Exception::abortnow;
+ return 0.;
+}
+
+// If needed, insert default implementations of virtual function defined
+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
+
+
+void DipoleMatching::persistentOutput(PersistentOStream &) const {}
+
+void DipoleMatching::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<DipoleMatching,Herwig::ShowerApproximation>
+ describeHerwigDipoleMatching("Herwig::DipoleMatching", "HwMatchbox.so");
+
+void DipoleMatching::Init() {
+
+ static ClassDocumentation<DipoleMatching> documentation
+ ("DipoleMatching implements naive NLO matching with the dipole shower.");
+
+}
+
diff --git a/MatrixElement/Matchbox/Matching/DipoleMatching.h b/MatrixElement/Matchbox/Matching/DipoleMatching.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Matching/DipoleMatching.h
@@ -0,0 +1,122 @@
+// -*- C++ -*-
+//
+// DipoleMatching.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_DipoleMatching_H
+#define Herwig_DipoleMatching_H
+//
+// This is the declaration of the DipoleMatching class.
+//
+
+#include "Herwig++/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * \ingroup Matchbox
+ * \author Simon Platzer
+ *
+ * \brief DipoleMatching implements naive NLO matching with the dipole shower.
+ *
+ */
+class DipoleMatching: public Herwig::ShowerApproximation {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ DipoleMatching();
+
+ /**
+ * The destructor.
+ */
+ virtual ~DipoleMatching();
+ //@}
+
+public:
+
+ /**
+ * Return the shower approximation to the real emission cross
+ * section for the given pair of Born and real emission
+ * configurations.
+ */
+ virtual CrossSection dSigHatDR() const;
+
+ /**
+ * Return the shower approximation splitting kernel for the given
+ * pair of Born and real emission configurations in units of the
+ * Born center of mass energy squared, and including a weight to
+ * project onto the splitting given by the dipole used.
+ */
+ virtual double me2() const;
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name 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.
+ */
+ DipoleMatching & operator=(const DipoleMatching &);
+
+};
+
+}
+
+#endif /* Herwig_DipoleMatching_H */
diff --git a/MatrixElement/Matchbox/Matching/MEMatching.cc b/MatrixElement/Matchbox/Matching/MEMatching.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Matching/MEMatching.cc
@@ -0,0 +1,173 @@
+// -*- C++ -*-
+//
+// MEMatching.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 MEMatching class.
+//
+
+#include "MEMatching.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+#include "ThePEG/PDT/EnumParticles.h"
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
+
+using namespace Herwig;
+
+MEMatching::MEMatching()
+ : theBornScreening(true),
+ theScreeningPower(2.0) {}
+
+MEMatching::~MEMatching() {}
+
+IBPtr MEMatching::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr MEMatching::fullclone() const {
+ return new_ptr(*this);
+}
+
+double MEMatching::channelWeight(int emitter, int emission, int spectator) const {
+ // do the most simple thing for the time being; needs fixing later
+ if ( realCXComb()->mePartonData()[emission]->id() == ParticleID::g ) {
+ Energy2 pipk =
+ realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[spectator];
+ Energy2 pipj =
+ realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[emission];
+ Energy2 pjpk =
+ realCXComb()->meMomenta()[emission] * realCXComb()->meMomenta()[spectator];
+ return GeV2 * pipk / ( pipj * ( pipj + pjpk ) );
+ }
+ return
+ GeV2 / (realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[emission]);
+}
+
+double MEMatching::channelWeight() const {
+ double currentChannel = channelWeight(dipole()->realEmitter(),
+ dipole()->realEmission(),
+ dipole()->realSpectator());
+ if ( currentChannel == 0. )
+ return 0.;
+ double sum = 0.;
+ for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator dip =
+ dipole()->partnerDipoles().begin();
+ dip != dipole()->partnerDipoles().end(); ++dip )
+ sum += channelWeight((**dip).realEmitter(),
+ (**dip).realEmission(),
+ (**dip).realSpectator());
+ assert(sum > 0.0);
+ return currentChannel / sum;
+}
+
+double MEMatching::screeningME2() const {
+ return
+ pow(sqr(dipole()->lastPt())/bornXComb()->lastSHat(),screeningPower()) *
+ dipole()->underlyingBornME()->me2Norm();
+}
+
+CrossSection MEMatching::dSigHatDR() const {
+ double pdfFactor = 1.;
+ double bornPDF = bornPDFWeight(showerScalesInSubtraction());
+ double bornPDFHard = bornPDF;
+ if ( showerScalesInSubtraction() )
+ bornPDFHard = bornPDFWeight(false);
+ if ( bornScreening() ) {
+ double bornME2 = dipole()->underlyingBornME()->me2();
+ double screenME2 = screeningME2();
+ pdfFactor = bornME2 * bornPDFHard / ( bornME2 * bornPDF + screenME2 );
+ } else {
+ pdfFactor = bornPDFHard / bornPDF;
+ }
+ assert(realXComb()->lastME2() > 0.0);
+ return
+ sqr(hbarc) *
+ realXComb()->jacobian() *
+ realPDFWeight(showerScalesInSubtraction()) *
+ couplingWeight(showerScalesInSubtraction()) *
+ pdfFactor *
+ channelWeight() * realXComb()->lastME2() /
+ (2. * realXComb()->lastSHat());
+}
+
+double MEMatching::me2() const {
+ double bornPDF = bornPDFWeight(showerScalesInSplitting());
+ double realPDF = realPDFWeight(showerScalesInSplitting());
+ assert(bornXComb()->lastME2() > 0.0);
+ double den =
+ bornXComb()->lastME2() * bornPDF;
+ if ( bornScreening() )
+ den += screeningME2();
+ double num =
+ dipole()->realEmissionME()->me2() * realPDF;
+ num *= pow(bornXComb()->lastSHat()/realXComb()->lastSHat(),2.*(realCXComb()->mePartonData().size())-8.);
+ return
+ (num/den) *
+ (bornXComb()->lastSHat()/realXComb()->lastSHat()) *
+ couplingWeight(showerScalesInSplitting());
+}
+
+// If needed, insert default implementations of virtual function defined
+// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
+
+
+void MEMatching::persistentOutput(PersistentOStream & os) const {
+ os << theBornScreening << theScreeningPower;
+}
+
+void MEMatching::persistentInput(PersistentIStream & is, int) {
+ is >> theBornScreening >> theScreeningPower;
+}
+
+
+// *** 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<MEMatching,Herwig::ShowerApproximation>
+ describeHerwigMEMatching("Herwig::MEMatching", "HwMatchbox.so");
+
+void MEMatching::Init() {
+
+ static ClassDocumentation<MEMatching> documentation
+ ("MEMatching implements NLO matching with matrix element correction (aka Powheg).");
+
+ static Switch<MEMatching,bool> interfaceBornScreening
+ ("BornScreening",
+ "Switch on or off Born screening",
+ &MEMatching::theBornScreening, true, false, false);
+ static SwitchOption interfaceBornScreeningOn
+ (interfaceBornScreening,
+ "On",
+ "Perform Born screening",
+ true);
+ static SwitchOption interfaceBornScreeningOff
+ (interfaceBornScreening,
+ "Off",
+ "Do not perform Born screening",
+ false);
+
+ static Parameter<MEMatching,double> interfaceScreeningPower
+ ("ScreeningPower",
+ "Set the power of pt used in the screening term",
+ &MEMatching::theScreeningPower, 2.0, 1.0, 0,
+ false, false, Interface::lowerlim);
+
+}
+
diff --git a/MatrixElement/Matchbox/Matching/MEMatching.h b/MatrixElement/Matchbox/Matching/MEMatching.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/Matchbox/Matching/MEMatching.h
@@ -0,0 +1,165 @@
+// -*- C++ -*-
+//
+// MEMatching.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_MEMatching_H
+#define Herwig_MEMatching_H
+//
+// This is the declaration of the MEMatching class.
+//
+
+#include "Herwig++/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * \ingroup Matchbox
+ * \author Simon Platzer
+ *
+ * \brief MEMatching implements NLO matching with matrix element correction (aka Powheg).
+ *
+ */
+class MEMatching: public Herwig::ShowerApproximation {
+
+public:
+
+ /** @name Standard constructors and destructors. */
+ //@{
+ /**
+ * The default constructor.
+ */
+ MEMatching();
+
+ /**
+ * The destructor.
+ */
+ virtual ~MEMatching();
+ //@}
+
+public:
+
+ /**
+ * Return true, if this shower approximation will require a
+ * splitting generator
+ */
+ virtual bool needsSplittingGenerator() const { return true; }
+
+public:
+
+ /**
+ * Return the shower approximation to the real emission cross
+ * section for the given pair of Born and real emission
+ * configurations.
+ */
+ virtual CrossSection dSigHatDR() const;
+
+ /**
+ * Return the shower approximation splitting kernel for the given
+ * pair of Born and real emission configurations in units of the
+ * Born center of mass energy squared, and including a weight to
+ * project onto the splitting given by the dipole used.
+ */
+ virtual double me2() const;
+
+ /**
+ * Generate a weight for the given dipole channel
+ */
+ double channelWeight(int emitter, int emission, int spectator) const;
+
+ /**
+ * Generate a normalized weight taking into account all channels
+ */
+ double channelWeight() const;
+
+ /**
+ * Return true, if 'Born screening' is taken into account
+ */
+ bool bornScreening() const { return theBornScreening; }
+
+ /**
+ * Return the power of pt used in the screening term
+ */
+ double screeningPower() const { return theScreeningPower; }
+
+ /**
+ * Return the screening `matrix element squared'
+ */
+ double screeningME2() 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:
+
+ /**
+ * True, if 'Born screening' should be done
+ */
+ bool theBornScreening;
+
+ /**
+ * The power of pt used in the screening term
+ */
+ double theScreeningPower;
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ MEMatching & operator=(const MEMatching &);
+
+};
+
+}
+
+#endif /* Herwig_MEMatching_H */
diff --git a/MatrixElement/Matchbox/Matching/Makefile.am b/MatrixElement/Matchbox/Matching/Makefile.am
--- a/MatrixElement/Matchbox/Matching/Makefile.am
+++ b/MatrixElement/Matchbox/Matching/Makefile.am
@@ -1,11 +1,15 @@
if WANT_DIPOLE
noinst_LTLIBRARIES = libHwMatchboxMatching.la
endif
libHwMatchboxMatching_la_SOURCES = \
ShowerApproximation.h \
ShowerApproximation.cc \
ShowerApproximationKernel.h \
ShowerApproximationKernel.cc \
ShowerApproximationGenerator.h \
-ShowerApproximationGenerator.cc
+ShowerApproximationGenerator.cc \
+DipoleMatching.h \
+DipoleMatching.cc \
+MEMatching.h \
+MEMatching.cc
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximation.cc b/MatrixElement/Matchbox/Matching/ShowerApproximation.cc
--- a/MatrixElement/Matchbox/Matching/ShowerApproximation.cc
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximation.cc
@@ -1,66 +1,252 @@
// -*- C++ -*-
//
// ShowerApproximation.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 ShowerApproximation class.
//
#include "ShowerApproximation.h"
#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
using namespace Herwig;
ShowerApproximation::ShowerApproximation()
- : HandlerBase(), theBelowCutoff(false) {}
+ : HandlerBase(), theBelowCutoff(false),
+ theFFPtCut(1.0*GeV), theFIPtCut(1.0*GeV), theIIPtCut(1.0*GeV),
+ theShowerScalesInSubtraction(false),
+ theShowerScalesInSplitting(true),
+ theRestrictPhasespace(true), theHardScaleFactor(1.0),
+ theExtrapolationX(0.65) {}
ShowerApproximation::~ShowerApproximation() {}
void ShowerApproximation::setDipole(Ptr<SubtractionDipole>::tcptr dip) { theDipole = dip; }
Ptr<SubtractionDipole>::tcptr ShowerApproximation::dipole() const { return theDipole; }
+bool ShowerApproximation::isAboveCutoff() const {
+
+ if ( dipole()->bornEmitter() > 1 &&
+ dipole()->bornSpectator() > 1 ) {
+ return dipole()->lastPt() > ffPtCut();
+ } else if ( ( dipole()->bornEmitter() > 1 &&
+ dipole()->bornSpectator() < 2 ) ||
+ ( dipole()->bornEmitter() < 2 &&
+ dipole()->bornSpectator() > 1 ) ) {
+ return dipole()->lastPt() > fiPtCut();
+ } else {
+ assert(dipole()->bornEmitter() < 2 &&
+ dipole()->bornSpectator() < 2);
+ return dipole()->lastPt() > iiPtCut();
+ }
+
+ return true;
+
+}
+
+bool ShowerApproximation::isInShowerPhasespace() const {
+
+ if ( !isAboveCutoff() )
+ return false;
+ if ( !restrictPhasespace() )
+ return true;
+
+ Energy maxPt = generator()->maximumCMEnergy();
+ vector<Lorentz5Momentum>::const_iterator p =
+ bornCXComb()->meMomenta().begin() + 2;
+ cPDVector::const_iterator pp =
+ bornCXComb()->mePartonData().begin() + 2;
+ for ( ; p != bornCXComb()->meMomenta().end(); ++p, ++pp )
+ if ( (**pp).coloured() )
+ maxPt = min(maxPt,p->perp());
+ if ( maxPt == generator()->maximumCMEnergy() )
+ maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m();
+ maxPt *= sqrt(hardScaleFactor());
+
+ return dipole()->lastPt() <= maxPt;
+
+}
+
+double ShowerApproximation::couplingWeight(bool showerscales) const {
+ if ( !showerscales )
+ return 1.;
+ double hardAlpha = dipole()->realEmissionME()->lastAlphaS();
+ Energy2 mur = sqr(dipole()->lastPt());
+ mur *= dipole()->realEmissionME()->renormalizationScaleFactor();
+ double runAlpha = SM().alphaS(mur);
+ return runAlpha/hardAlpha;
+}
+
+double ShowerApproximation::bornPDFWeight(bool showerscales) const {
+ if ( !bornCXComb()->mePartonData()[0]->coloured() &&
+ !bornCXComb()->mePartonData()[1]->coloured() )
+ return 1.;
+ Energy2 muf;
+ if ( showerscales ) {
+ muf = sqr(dipole()->lastPt());
+ } else {
+ muf = dipole()->underlyingBornME()->factorizationScale();
+ }
+ muf *= dipole()->underlyingBornME()->factorizationScaleFactor();
+ double pdfweight = 1.;
+ if ( bornCXComb()->mePartonData()[0]->coloured() &&
+ dipole()->underlyingBornME()->havePDFWeight1() )
+ pdfweight *= dipole()->underlyingBornME()->pdf1(muf,theExtrapolationX);
+ if ( bornCXComb()->mePartonData()[1]->coloured() &&
+ dipole()->underlyingBornME()->havePDFWeight2() )
+ pdfweight *= dipole()->underlyingBornME()->pdf2(muf,theExtrapolationX);
+ return pdfweight;
+}
+
+double ShowerApproximation::realPDFWeight(bool showerscales) const {
+ if ( !realCXComb()->mePartonData()[0]->coloured() &&
+ !realCXComb()->mePartonData()[1]->coloured() )
+ return 1.;
+ Energy2 muf;
+ if ( showerscales ) {
+ muf = sqr(dipole()->lastPt());
+ } else {
+ muf = dipole()->realEmissionME()->factorizationScale();
+ }
+ muf *= dipole()->realEmissionME()->factorizationScaleFactor();
+ double pdfweight = 1.;
+ if ( realCXComb()->mePartonData()[0]->coloured() &&
+ dipole()->realEmissionME()->havePDFWeight1() )
+ pdfweight *= dipole()->realEmissionME()->pdf1(muf,theExtrapolationX);
+ if ( realCXComb()->mePartonData()[1]->coloured() &&
+ dipole()->realEmissionME()->havePDFWeight2() )
+ pdfweight *= dipole()->realEmissionME()->pdf2(muf,theExtrapolationX);
+ return pdfweight;
+}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void ShowerApproximation::persistentOutput(PersistentOStream & os) const {
- os << theBornXComb << theRealXComb << theDipole << theBelowCutoff;
+ os << theBornXComb << theRealXComb << theTildeXCombs << theDipole << theBelowCutoff
+ << ounit(theFFPtCut,GeV) << ounit(theFIPtCut,GeV)
+ << ounit(theIIPtCut,GeV) << theShowerScalesInSubtraction
+ << theShowerScalesInSplitting
+ << theRestrictPhasespace << theHardScaleFactor
+ << theExtrapolationX;
}
void ShowerApproximation::persistentInput(PersistentIStream & is, int) {
- is >> theBornXComb >> theRealXComb >> theDipole >> theBelowCutoff;
+ is >> theBornXComb >> theRealXComb >> theTildeXCombs >> theDipole >> theBelowCutoff
+ >> iunit(theFFPtCut,GeV) >> iunit(theFIPtCut,GeV)
+ >> iunit(theIIPtCut,GeV) >> theShowerScalesInSubtraction
+ >> theShowerScalesInSplitting
+ >> theRestrictPhasespace >> theHardScaleFactor
+ >> theExtrapolationX;
}
// *** 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<ShowerApproximation,HandlerBase>
describeHerwigShowerApproximation("Herwig::ShowerApproximation", "HwMatchbox.so");
void ShowerApproximation::Init() {
static ClassDocumentation<ShowerApproximation> documentation
("ShowerApproximation describes the shower emission to be used "
"in NLO matching.");
+ static Parameter<ShowerApproximation,Energy> interfaceFFPtCut
+ ("FFPtCut",
+ "Set the pt infrared cutoff",
+ &ShowerApproximation::theFFPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
+ false, false, Interface::lowerlim);
+
+ static Parameter<ShowerApproximation,Energy> interfaceFIPtCut
+ ("FIPtCut",
+ "Set the pt infrared cutoff",
+ &ShowerApproximation::theFIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
+ false, false, Interface::lowerlim);
+
+ static Parameter<ShowerApproximation,Energy> interfaceIIPtCut
+ ("IIPtCut",
+ "Set the pt infrared cutoff",
+ &ShowerApproximation::theIIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
+ false, false, Interface::lowerlim);
+
+ static Switch<ShowerApproximation,bool> interfaceShowerScalesInSubtraction
+ ("ShowerScalesInSubtraction",
+ "Switch on or off shower scales in the matching subtraction",
+ &ShowerApproximation::theShowerScalesInSubtraction, true, false, false);
+ static SwitchOption interfaceShowerScalesInSubtractionOn
+ (interfaceShowerScalesInSubtraction,
+ "On",
+ "Use shower scales in the matching subtraction",
+ true);
+ static SwitchOption interfaceShowerScalesInSubtractionOff
+ (interfaceShowerScalesInSubtraction,
+ "Off",
+ "Use hard process scales in the matching subtraction",
+ false);
+
+ static Switch<ShowerApproximation,bool> interfaceShowerScalesInSplitting
+ ("ShowerScalesInSplitting",
+ "Switch on or off shower scales in the splitting generation",
+ &ShowerApproximation::theShowerScalesInSplitting, true, false, false);
+ static SwitchOption interfaceShowerScalesInSplittingOn
+ (interfaceShowerScalesInSplitting,
+ "On",
+ "Use shower scales in the matching splitting generation",
+ true);
+ static SwitchOption interfaceShowerScalesInSplittingOff
+ (interfaceShowerScalesInSplitting,
+ "Off",
+ "Use hard process scales in the matching splitting generation",
+ false);
+
+ static Switch<ShowerApproximation,bool> interfaceRestrictPhasespace
+ ("RestrictPhasespace",
+ "Switch on or off phasespace restrictions",
+ &ShowerApproximation::theRestrictPhasespace, true, false, false);
+ static SwitchOption interfaceRestrictPhasespaceOn
+ (interfaceRestrictPhasespace,
+ "On",
+ "Perform phasespace restrictions",
+ true);
+ static SwitchOption interfaceRestrictPhasespaceOff
+ (interfaceRestrictPhasespace,
+ "Off",
+ "Do not perform phasespace restrictions",
+ false);
+
+ static Parameter<ShowerApproximation,double> interfaceHardScaleFactor
+ ("HardScaleFactor",
+ "The hard scale factor.",
+ &ShowerApproximation::theHardScaleFactor, 1.0, 0.0, 0,
+ false, false, Interface::lowerlim);
+
+ static Parameter<ShowerApproximation,double> interfaceExtrapolationX
+ ("ExtrapolationX",
+ "The x from which on extrapolation should be performed.",
+ &ShowerApproximation::theExtrapolationX, 0.65, 0.0, 1.0,
+ false, false, Interface::limited);
+
}
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,202 +1,325 @@
// -*- 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; }
+ /**
+ * Return true, if this shower approximation will require tilde
+ * XCombs for the real phase space point generated
+ */
+ virtual bool needsTildeXCombs() 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; }
/**
+ * Return the XComb object describing the Born process
+ */
+ tcStdXCombPtr bornCXComb() 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; }
/**
+ * Return the XComb object describing the real emission process
+ */
+ tcStdXCombPtr realCXComb() const { return theRealXComb; }
+
+ /**
+ * Set the tilde xcomb objects associated to the real xcomb
+ */
+ void setTildeXCombs(const vector<StdXCombPtr>& xc) { theTildeXCombs = xc; }
+
+ /**
+ * Return the tilde xcomb objects associated to the real xcomb
+ */
+ const vector<StdXCombPtr>& tildeXCombs() const { return theTildeXCombs; }
+
+ /**
* 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; }
+ /**
+ * Return the pt cut to be applied for final-final dipoles.
+ */
+ Energy ffPtCut() const { return theFFPtCut; }
+
+ /**
+ * Return the pt cut to be applied for final-initial dipoles.
+ */
+ Energy fiPtCut() const { return theFIPtCut; }
+
+ /**
+ * Return the pt cut to be applied for initial-initial dipoles.
+ */
+ Energy iiPtCut() const { return theIIPtCut; }
+
public:
/**
+ * Return true, if the phase space restrictions of the dipole shower should
+ * be applied.
+ */
+ bool restrictPhasespace() const { return theRestrictPhasespace; }
+
+ /**
+ * Return the scale factor for the hard scale
+ */
+ double hardScaleFactor() const { return theHardScaleFactor; }
+
+ /**
* 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;
+ virtual bool isInShowerPhasespace() const;
/**
* 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;
+ virtual bool isAboveCutoff() const;
/**
* Return the shower approximation to the real emission cross
* section for the given pair of Born and real emission
* configurations.
*/
virtual CrossSection dSigHatDR() const = 0;
/**
* Return the shower approximation splitting kernel for the given
* pair of Born and real emission configurations in units of the
* Born center of mass energy squared, and including a weight to
* project onto the splitting given by the dipole used.
*/
virtual double me2() const = 0;
+ /**
+ * Return true, if the shower scales should be used in the subtraction
+ */
+ bool showerScalesInSubtraction() const { return theShowerScalesInSubtraction; }
+
+ /**
+ * Return true, if the shower scales should be used in splitting generation
+ */
+ bool showerScalesInSplitting() const { return theShowerScalesInSplitting; }
+
+ /**
+ * Return the running coupling weight
+ */
+ double couplingWeight(bool showerscales) const;
+
+ /**
+ * Return the Born PDF weight
+ */
+ double bornPDFWeight(bool showerscales) const;
+
+ /**
+ * Return the real emission PDF weight
+ */
+ double realPDFWeight(bool showerscales) 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).
private:
/**
* The XComb object describing the Born process
*/
tStdXCombPtr theBornXComb;
/**
* The XComb object describing the real emission process
*/
tStdXCombPtr theRealXComb;
/**
+ * The tilde xcomb objects associated to the real xcomb
+ */
+ vector<StdXCombPtr> theTildeXCombs;
+
+ /**
* 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;
+ /**
+ * The pt cut to be applied for final-final dipoles.
+ */
+ Energy theFFPtCut;
+
+ /**
+ * The pt cut to be applied for final-initial dipoles.
+ */
+ Energy theFIPtCut;
+
+ /**
+ * The pt cut to be applied for initial-initial dipoles.
+ */
+ Energy theIIPtCut;
+
+ /**
+ * True, if the shower scales should be used in the subtraction
+ */
+ bool theShowerScalesInSubtraction;
+
+ /**
+ * True, if the shower scales should be used in splitting generation
+ */
+ bool theShowerScalesInSplitting;
+
+ /**
+ * True, if the phase space restrictions of the dipole shower should
+ * be applied.
+ */
+ bool theRestrictPhasespace;
+
+ /**
+ * The scale factor for the hard scale
+ */
+ double theHardScaleFactor;
+
+ /**
+ * The x value from which on we extrapolate PDFs for numerically stable ratios.
+ */
+ double theExtrapolationX;
+
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/Matching/ShowerApproximationGenerator.cc b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
--- a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
@@ -1,418 +1,397 @@
// -*- C++ -*-
//
// ShowerApproximationGenerator.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 ShowerApproximationGenerator class.
//
#include "ShowerApproximationGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
ShowerApproximationGenerator::ShowerApproximationGenerator()
: thePresamplingPoints(10000), theMaxTry(100000) {}
ShowerApproximationGenerator::~ShowerApproximationGenerator() {}
double ShowerApproximationGenerator::generateFraction(tcPDPtr pd, double r, double xmin) const {
if ( pd->coloured() || pd->id() == ParticleID::gamma ) {
return pow(xmin,r);
}
double x0 = 1.e-5;
return 1. + x0 - x0*pow((1.+x0)/x0,r);
}
double ShowerApproximationGenerator::invertFraction(tcPDPtr pd, double x, double xmin) const {
if ( pd->coloured() || pd->id() == ParticleID::gamma ) {
return log(x)/log(xmin);
}
double x0 = 1.e-5;
return log((1.-x+x0)/x0)/log((1.+x0)/x0);
}
bool ShowerApproximationGenerator::prepare() {
tSubProPtr oldSub = lastIncomingXComb->subProcess();
tcStdXCombPtr cIncomingXC = lastIncomingXComb;
theLastMomenta.resize(cIncomingXC->mePartonData().size());
theLastRandomNumbers.resize(thePhasespace->nDim(theLastMomenta.size()-2) + 2);
double x1 =
oldSub->incoming().first->momentum().plus() /
lastIncomingXComb->lastParticles().first->momentum().plus();
theLastRandomNumbers[0] = invertFraction(oldSub->incoming().first->dataPtr(),x1,
lastIncomingXComb->cuts()->x1Min());
double x2 =
oldSub->incoming().second->momentum().minus() /
lastIncomingXComb->lastParticles().second->momentum().minus();
theLastRandomNumbers[1] = invertFraction(oldSub->incoming().second->dataPtr(),x2,
lastIncomingXComb->cuts()->x2Min());
Boost toCMS =
(oldSub->incoming().first->momentum() +
oldSub->incoming().second->momentum()).findBoostToCM();
theLastMomenta[0] = oldSub->incoming().first->momentum();
theLastMomenta[0].boost(toCMS);
theLastMomenta[1] = oldSub->incoming().second->momentum();
theLastMomenta[1].boost(toCMS);
ParticleVector::const_iterator out = oldSub->outgoing().begin();
vector<Lorentz5Momentum>::iterator p = theLastMomenta.begin() + 2;
for ( ; out != oldSub->outgoing().end(); ++out, ++p ) {
*p = (**out).momentum();
p->boost(toCMS);
}
theLastPartons.first =
oldSub->incoming().first->data().produceParticle(oldSub->incoming().first->momentum());
theLastPartons.second =
oldSub->incoming().second->data().produceParticle(oldSub->incoming().second->momentum());
thePhasespace->invertKinematics(theLastMomenta,&theLastRandomNumbers[2]);
theLastBornXComb->clean();
theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
theLastMomenta,theLastRandomNumbers);
if ( !theLastBornXComb->cuts()->initSubProcess(theLastBornXComb->lastSHat(),
theLastBornXComb->lastY(),
theLastBornXComb->mirror()) )
return false;
theLastBornME->setXComb(theLastBornXComb);
if ( !theLastBornME->generateKinematics(&theLastRandomNumbers[2]) )
return false;
CrossSection bornXS = theLastBornME->dSigHatDR();
if ( bornXS == ZERO )
return false;
return true;
}
bool ShowerApproximationGenerator::generate(const vector<double>& r) {
theLastBornXComb->clean();
double x = generateFraction(theLastPartons.first->dataPtr(),r[0],
lastIncomingXComb->cuts()->x1Min());
Energy Q = lastIncomingXComb->lastParticles().first->momentum().plus();
Energy mass = theLastPartons.first->dataPtr()->mass();
double xi = (4.*sqr(x*Q) - sqr(mass))/(4.*sqr(Q)*x);
Lorentz5Momentum p1(ZERO,ZERO,xi*Q);
p1.setMass(mass); p1.rescaleEnergy();
theLastPartons.first->set5Momentum(p1);
x = generateFraction(theLastPartons.second->dataPtr(),r[1],
lastIncomingXComb->cuts()->x1Min());
Q = lastIncomingXComb->lastParticles().second->momentum().minus();
mass = theLastPartons.second->dataPtr()->mass();
xi = (4.*sqr(x*Q) - sqr(mass))/(4.*sqr(Q)*x);
Lorentz5Momentum p2(ZERO,ZERO,-xi*Q);
p2.setMass(mass); p2.rescaleEnergy();
theLastPartons.second->set5Momentum(p2);
theLastPresamplingMomenta.resize(theLastMomenta.size());
Boost toCMS =
(theLastPartons.first->momentum() +
theLastPartons.second->momentum()).findBoostToCM();
theLastPresamplingMomenta[0] = theLastPartons.first->momentum();
theLastPresamplingMomenta[0].boost(toCMS);
theLastPresamplingMomenta[1] = theLastPartons.second->momentum();
theLastPresamplingMomenta[1].boost(toCMS);
theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
theLastPresamplingMomenta,r);
if ( !theLastBornXComb->cuts()->initSubProcess(theLastBornXComb->lastSHat(),
theLastBornXComb->lastY(),
theLastBornXComb->mirror()) )
return false;
theLastBornME->setXComb(theLastBornXComb);
if ( !theLastBornME->generateKinematics(&r[2]) )
return false;
CrossSection bornXS = theLastBornME->dSigHatDR();
if ( bornXS == ZERO )
return false;
return true;
}
void ShowerApproximationGenerator::restore() {
tSubProPtr oldSub = lastIncomingXComb->subProcess();
theLastPartons.first->set5Momentum(oldSub->incoming().first->momentum());
theLastPartons.second->set5Momentum(oldSub->incoming().second->momentum());
theLastBornXComb->clean();
theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
theLastMomenta,theLastRandomNumbers);
theLastBornME->setXComb(theLastBornXComb);
theLastBornME->generateKinematics(&theLastRandomNumbers[2]);
theLastBornME->dSigHatDR();
}
void ShowerApproximationGenerator::
handle(EventHandler & eh, const tPVector &,
const Hint &) {
lastIncomingXComb = dynamic_ptr_cast<tStdXCombPtr>(eh.lastXCombPtr());
if ( !lastIncomingXComb )
throw Exception() << "expecting a standard event handler"
<< Exception::abortnow;
if ( lastIncomingXComb->lastProjector() )
lastIncomingXComb = lastIncomingXComb->lastProjector();
const StandardXComb& xc = *lastIncomingXComb;
map<cPDVector,set<Ptr<ShowerApproximationKernel>::ptr> >::const_iterator
kernelit = theKernelMap.find(xc.mePartonData());
if ( kernelit == theKernelMap.end() ) {
list<MatchboxFactory::SplittingChannel> channels =
theFactory->getSplittingChannels(lastIncomingXComb);
set<Ptr<ShowerApproximationKernel>::ptr> newKernels;
for ( list<MatchboxFactory::SplittingChannel>::const_iterator c =
channels.begin(); c != channels.end(); ++c ) {
Ptr<ShowerApproximationKernel>::ptr kernel =
new_ptr(ShowerApproximationKernel());
kernel->setBornXComb(c->bornXComb);
kernel->setRealXComb(c->realXComb);
+ kernel->setTildeXCombs(c->tildeXCombs);
kernel->setDipole(c->dipole);
kernel->showerApproximation(theShowerApproximation);
kernel->presamplingPoints(thePresamplingPoints);
kernel->maxtry(theMaxTry);
if ( kernel->dipole()->bornEmitter() > 1 &&
kernel->dipole()->bornSpectator() > 1 ) {
- kernel->ptCut(theFFPtCut);
+ kernel->ptCut(ffPtCut());
} else if ( ( kernel->dipole()->bornEmitter() > 1 &&
kernel->dipole()->bornSpectator() < 2 ) ||
( kernel->dipole()->bornEmitter() < 2 &&
kernel->dipole()->bornSpectator() > 1 ) ) {
- kernel->ptCut(theFIPtCut);
+ kernel->ptCut(fiPtCut());
} else {
assert(kernel->dipole()->bornEmitter() < 2 &&
kernel->dipole()->bornSpectator() < 2);
- kernel->ptCut(theIIPtCut);
+ kernel->ptCut(iiPtCut());
}
newKernels.insert(kernel);
}
theKernelMap[xc.mePartonData()] = newKernels;
kernelit = theKernelMap.find(xc.mePartonData());
}
if ( kernelit->second.empty() )
return;
const set<Ptr<ShowerApproximationKernel>::ptr>& kernels = kernelit->second;
theLastBornME = (**kernels.begin()).dipole()->underlyingBornME();
theLastBornXComb = (**kernels.begin()).bornXComb();
if ( !prepare() )
return;
Energy winnerPt = ZERO;
Ptr<ShowerApproximationKernel>::ptr winnerKernel;
for ( set<Ptr<ShowerApproximationKernel>::ptr>::const_iterator k =
kernels.begin(); k != kernels.end(); ++k ) {
if ( (**k).generate() != 0. ) {
winnerKernel = *k;
winnerPt = max(winnerPt,winnerKernel->dipole()->lastPt());
}
}
if ( !winnerKernel || winnerPt == ZERO )
return;
SubProPtr oldSub = lastIncomingXComb->subProcess();
SubProPtr newSub;
try {
newSub = winnerKernel->realXComb()->construct();
} catch(Veto&) {
return;
}
tParticleSet firstS = oldSub->incoming().first->siblings();
assert(firstS.empty() || firstS.size() == 1);
if ( !firstS.empty() ) {
eh.currentStep()->removeParticle(*firstS.begin());
}
tParticleSet secondS = oldSub->incoming().second->siblings();
assert(secondS.empty() || secondS.size() == 1);
if ( !secondS.empty() ) {
eh.currentStep()->removeParticle(*secondS.begin());
}
// prevent the colliding particles from disappearing
// in the initial state and appearing
// in the final state when we've cut off all their
// (physical) children; only applies to the case
// where we have a parton extractor not build from
// noPDF, so check wether the incoming particle
// doesnt equal the incoming parton -- this needs fixing in ThePEG
PPtr dummy = new_ptr(Particle(getParticleData(ParticleID::gamma)));
bool usedDummy = false;
if ( eh.currentStep()->incoming().first != oldSub->incoming().first ) {
eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().first,dummy);
usedDummy = true;
}
if ( eh.currentStep()->incoming().second != oldSub->incoming().second ) {
eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().second,dummy);
usedDummy = true;
}
eh.currentStep()->removeSubProcess(oldSub);
eh.currentStep()->addSubProcess(newSub);
// get rid of the dummy
if ( usedDummy ) {
eh.currentStep()->removeParticle(dummy);
}
eh.select(winnerKernel->realXComb());
winnerKernel->realXComb()->
recreatePartonBinInstances(winnerKernel->realXComb()->lastScale());
eh.lastExtractor()->constructRemnants(winnerKernel->realXComb()->partonBinInstances(),
newSub, eh.currentStep());
}
IBPtr ShowerApproximationGenerator::clone() const {
return new_ptr(*this);
}
IBPtr ShowerApproximationGenerator::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 ShowerApproximationGenerator::persistentOutput(PersistentOStream & os) const {
os << theShowerApproximation << thePhasespace
<< theFactory << theKernelMap
- << ounit(theFFPtCut,GeV) << ounit(theFIPtCut,GeV)
- << ounit(theIIPtCut,GeV)
<< thePresamplingPoints << theMaxTry;
}
void ShowerApproximationGenerator::persistentInput(PersistentIStream & is, int) {
is >> theShowerApproximation >> thePhasespace
>> theFactory >> theKernelMap
- >> iunit(theFFPtCut,GeV) >> iunit(theFIPtCut,GeV)
- >> iunit(theIIPtCut,GeV)
>> thePresamplingPoints >> theMaxTry;
}
// *** 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<ShowerApproximationGenerator,StepHandler>
describeHerwigShowerApproximationGenerator("Herwig::ShowerApproximationGenerator", "HwMatchbox.so");
void ShowerApproximationGenerator::Init() {
static ClassDocumentation<ShowerApproximationGenerator> documentation
("ShowerApproximationGenerator generates emissions according to a "
"shower approximation entering a NLO matching.");
static Reference<ShowerApproximationGenerator,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to sample.",
&ShowerApproximationGenerator::theShowerApproximation, false, false, true, false, false);
static Reference<ShowerApproximationGenerator,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"The phase space generator to use.",
&ShowerApproximationGenerator::thePhasespace, false, false, true, false, false);
static Reference<ShowerApproximationGenerator,MatchboxFactory> interfaceFactory
("Factory",
"The factory object to use.",
&ShowerApproximationGenerator::theFactory, false, false, true, false, false);
- static Parameter<ShowerApproximationGenerator,Energy> interfaceFFPtCut
- ("FFPtCut",
- "Set the pt infrared cutoff",
- &ShowerApproximationGenerator::theFFPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
- false, false, Interface::lowerlim);
-
- static Parameter<ShowerApproximationGenerator,Energy> interfaceFIPtCut
- ("FIPtCut",
- "Set the pt infrared cutoff",
- &ShowerApproximationGenerator::theFIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
- false, false, Interface::lowerlim);
-
- static Parameter<ShowerApproximationGenerator,Energy> interfaceIIPtCut
- ("IIPtCut",
- "Set the pt infrared cutoff",
- &ShowerApproximationGenerator::theIIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
- false, false, Interface::lowerlim);
-
static Parameter<ShowerApproximationGenerator,unsigned long> interfacePresamplingPoints
("PresamplingPoints",
"Set the number of presampling points.",
&ShowerApproximationGenerator::thePresamplingPoints, 10000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerApproximationGenerator,unsigned long> interfaceMaxTry
("MaxTry",
"Set the number of maximum attempts.",
&ShowerApproximationGenerator::theMaxTry, 100000, 1, 0,
false, false, Interface::lowerlim);
}
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.h b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.h
--- a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.h
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.h
@@ -1,245 +1,245 @@
// -*- C++ -*-
//
// ShowerApproximationGenerator.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_ShowerApproximationGenerator_H
#define Herwig_ShowerApproximationGenerator_H
//
// This is the declaration of the ShowerApproximationGenerator class.
//
#include "ThePEG/Handlers/StepHandler.h"
#include "Herwig++/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief ShowerApproximationGenerator generates emissions according to a
* shower approximation entering a NLO matching.
*
*/
class ShowerApproximationGenerator: public StepHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ShowerApproximationGenerator();
/**
* The destructor.
*/
virtual ~ShowerApproximationGenerator();
//@}
public:
/** @name Virtual functions required by the StepHandler class. */
//@{
/**
* The main function called by the EventHandler class to
* perform a step. Given the current state of an Event, this function
* performs the event generation step and includes the result in a new
* Step object int the Event record.
* @param eh the EventHandler in charge of the Event generation.
* @param tagged if not empty these are the only particles which should
* be considered by the StepHandler.
* @param hint a Hint object with possible information from previously
* performed steps.
* @throws Veto if the StepHandler requires the current step to be discarded.
* @throws Stop if the generation of the current Event should be stopped
* after this call.
* @throws Exception if something goes wrong.
*/
virtual void handle(EventHandler & eh, const tPVector & tagged,
const Hint & hint);
//@}
public:
/**
* Fill information on the Born process
*/
bool prepare();
/**
* Generate a Born phase space point while kernels are being
* presampled
*/
bool generate(const vector<double>&);
/**
* Restore information on the Born process
*/
void restore();
protected:
/**
* Generate a momentum fraction for the given parton species
*/
double generateFraction(tcPDPtr, double, double) const;
/**
* Invert the momentum fraction for the given parton species
*/
double invertFraction(tcPDPtr, double, double) const;
+ /**
+ * Return the pt cut to be applied for final-final dipoles.
+ */
+ Energy ffPtCut() const { return theShowerApproximation->ffPtCut();; }
+
+ /**
+ * Return the pt cut to be applied for final-initial dipoles.
+ */
+ Energy fiPtCut() const { return theShowerApproximation->fiPtCut(); }
+
+ /**
+ * Return the pt cut to be applied for initial-initial dipoles.
+ */
+ Energy iiPtCut() const { return theShowerApproximation->iiPtCut(); }
+
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 shower approximation to consider
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* The (invertible) phase space generator to use
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The factory object to fetch splitting channels from
*/
Ptr<MatchboxFactory>::ptr theFactory;
/**
* Map hard processes to the respective kernels.
*/
map<cPDVector,set<Ptr<ShowerApproximationKernel>::ptr> > theKernelMap;
/**
- * The pt cut to be applied for final-final dipoles.
- */
- Energy theFFPtCut;
-
- /**
- * The pt cut to be applied for final-initial dipoles.
- */
- Energy theFIPtCut;
-
- /**
- * The pt cut to be applied for initial-initial dipoles.
- */
- Energy theIIPtCut;
-
- /**
* The number of points to presample this splitting generator.
*/
unsigned long thePresamplingPoints;
/**
* The maximum number of trials to generate a splitting.
*/
unsigned long theMaxTry;
/**
* The last external Born XComb dealt with
*/
tStdXCombPtr lastIncomingXComb;
// the next three are filled from the incoming xcomb in the prepare method
/**
* The last internal Born matrix element dealt with
*/
Ptr<MatchboxMEBase>::ptr theLastBornME;
/**
* The last Born phase space point
*/
vector<Lorentz5Momentum> theLastMomenta;
/**
* The last Born phase space point used while presampling
*/
vector<Lorentz5Momentum> theLastPresamplingMomenta;
/**
* The random numbers which have produced the last Born phase space
* point.
*/
vector<double> theLastRandomNumbers;
/**
* The last internal Born XComb dealt with
*/
tStdXCombPtr theLastBornXComb;
/**
* The last internal incoming partons dealt with
*/
PPair theLastPartons;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerApproximationGenerator & operator=(const ShowerApproximationGenerator &);
};
}
#endif /* Herwig_ShowerApproximationGenerator_H */
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
--- a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
@@ -1,183 +1,190 @@
// -*- C++ -*-
//
// ShowerApproximationKernel.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 ShowerApproximationKernel class.
//
#include "ShowerApproximationKernel.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/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ShowerApproximationGenerator.h"
using namespace Herwig;
ShowerApproximationKernel::ShowerApproximationKernel()
: thePresampling(false), thePresamplingPoints(10000),
theMaxTry(100000), sampler(0) {}
ShowerApproximationKernel::~ShowerApproximationKernel() {}
IBPtr ShowerApproximationKernel::clone() const {
return new_ptr(*this);
}
IBPtr ShowerApproximationKernel::fullclone() const {
return new_ptr(*this);
}
void ShowerApproximationKernel::showerApproximationGenerator(Ptr<ShowerApproximationGenerator>::tptr gen) {
theShowerApproximationGenerator = gen;
}
Ptr<ShowerApproximationGenerator>::tptr ShowerApproximationKernel::showerApproximationGenerator() const {
return theShowerApproximationGenerator;
}
const vector<bool>& ShowerApproximationKernel::sampleFlags() {
if ( !theFlags.empty() )
return theFlags;
theFlags.resize(nDim(),false);
for ( int k = nDimBorn();
k < nDimBorn() + dipole()->nDimRadiation(); ++k )
theFlags[k] = true;
return theFlags;
}
const pair<vector<double>,vector<double> >& ShowerApproximationKernel::support() {
if ( !theSupport.first.empty() )
return theSupport;
vector<double> l(nDim(),0.0);
vector<double> u(nDim(),1.0);
theSupport.first = l;
theSupport.second = u;
return theSupport;
}
const vector<double>& ShowerApproximationKernel::parameterPoint() {
theLastParameterPoint.resize(nDim());
copy(bornCXComb()->lastRandomNumbers().begin(),
bornCXComb()->lastRandomNumbers().end(),
theLastParameterPoint.begin());
theLastParameterPoint[evolutionVariable()] = 1.;
return theLastParameterPoint;
}
void ShowerApproximationKernel::startPresampling() {
thePresampling = true;
}
void ShowerApproximationKernel::stopPresampling() {
showerApproximationGenerator()->restore();
thePresampling = false;
}
double ShowerApproximationKernel::evaluate(const vector<double>& r) {
if ( presampling() ) {
theLastBornPoint.resize(nDimBorn());
copy(r.begin(),r.begin()+nDimBorn(),theLastBornPoint.begin());
if ( !showerApproximationGenerator()->generate(theLastBornPoint) )
return 0.;
}
+ assert(dipole()->splitting());
realXComb()->clean();
dipole()->setXComb(realXComb());
+ for ( vector<StdXCombPtr>::const_iterator t = tildeXCombs().begin();
+ t != tildeXCombs().end(); ++t ) {
+ (**t).clean();
+ (**t).matrixElement()->setXComb(*t);
+ }
if ( !dipole()->generateKinematics(&r[nDimBorn()]) )
return 0.;
double jac = dipole()->invertedTildeKinematics()->jacobian();
showerApproximation()->setBornXComb(bornXComb());
showerApproximation()->setRealXComb(realXComb());
+ showerApproximation()->setTildeXCombs(tildeXCombs());
showerApproximation()->setDipole(dipole());
if ( !showerApproximation()->isInShowerPhasespace() )
return 0.;
return showerApproximation()->me2() * jac;
}
double ShowerApproximationKernel::generate() {
if ( !sampler ) {
sampler = new ExponentialGenerator();
sampler->sampling_parameters().maxtry = maxtry();
sampler->sampling_parameters().presampling_points = presamplingPoints();
sampler->function(this);
sampler->initialize();
}
double res = 0.;
while (true) {
try {
res = sampler->generate();
} catch (exsample::exponential_regenerate&) {
continue;
} catch (exsample::hit_and_miss_maxtry&) {
throw Veto();
} catch (exsample::selection_maxtry&) {
throw Veto();
}
break;
}
return res;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void ShowerApproximationKernel::persistentOutput(PersistentOStream & os) const {
os << theDipole << theShowerApproximation
- << theBornXComb << theRealXComb
+ << theBornXComb << theRealXComb << theTildeXCombs
<< thePresamplingPoints << theMaxTry << theFlags << theSupport
<< theLastParameterPoint << theShowerApproximationGenerator;
}
void ShowerApproximationKernel::persistentInput(PersistentIStream & is, int) {
is >> theDipole >> theShowerApproximation
- >> theBornXComb >> theRealXComb
+ >> theBornXComb >> theRealXComb >> theTildeXCombs
>> thePresamplingPoints >> theMaxTry >> theFlags >> theSupport
>> theLastParameterPoint >> theShowerApproximationGenerator;
}
// *** 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<ShowerApproximationKernel,HandlerBase>
describeHerwigShowerApproximationKernel("Herwig::ShowerApproximationKernel", "HwMatchbox.so");
void ShowerApproximationKernel::Init() {
static ClassDocumentation<ShowerApproximationKernel> documentation
("ShowerApproximationKernel generates emissions according to a "
"shower approximation entering a NLO matching.");
}
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
--- a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
@@ -1,407 +1,422 @@
// -*- C++ -*-
//
// ShowerApproximationKernel.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_ShowerApproximationKernel_H
#define Herwig_ShowerApproximationKernel_H
//
// This is the declaration of the ShowerApproximationKernel class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig++/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/Exsample2/exsample/exponential_generator.h"
namespace Herwig {
using namespace ThePEG;
class ShowerApproximationGenerator;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief ShowerApproximationKernel generates emissions according to a
* shower approximation entering a NLO matching.
*
*/
class ShowerApproximationKernel: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ShowerApproximationKernel();
/**
* The destructor.
*/
virtual ~ShowerApproximationKernel();
//@}
public:
/**
* Set the XComb object describing the Born process
*/
void setBornXComb(tStdXCombPtr xc) { theBornXComb = xc; }
/**
* Return the XComb object describing the Born process
*/
tcStdXCombPtr bornCXComb() const { return theBornXComb; }
/**
* 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
*/
tcStdXCombPtr realCXComb() const { return theRealXComb; }
/**
* Return the XComb object describing the real emission process
*/
tStdXCombPtr realXComb() const { return theRealXComb; }
/**
+ * Set the tilde xcomb objects associated to the real xcomb
+ */
+ void setTildeXCombs(const vector<StdXCombPtr>& xc) { theTildeXCombs = xc; }
+
+ /**
+ * Return the tilde xcomb objects associated to the real xcomb
+ */
+ const vector<StdXCombPtr>& tildeXCombs() const { return theTildeXCombs; }
+
+ /**
* Set the dipole in charge for the emission
*/
void setDipole(Ptr<SubtractionDipole>::tptr dip) { theDipole = dip; }
/**
* Return the dipole in charge for the emission
*/
Ptr<SubtractionDipole>::tptr dipole() const { return theDipole; }
/**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr app) { theShowerApproximation = app; }
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
/**
* Set the shower approximation generator.
*/
void showerApproximationGenerator(Ptr<ShowerApproximationGenerator>::tptr);
/**
* Return the shower approximation generator.
*/
Ptr<ShowerApproximationGenerator>::tptr showerApproximationGenerator() const;
/**
* Generate the next emission
*/
double generate();
public:
/**
* Set a pt cut on the dipole to generate the radiation
*/
void ptCut(Energy pt) { dipole()->ptCut(pt); }
/**
* Return the number of random numbers
* needed to sample this kernel.
*/
int nDim() const {
return
nDimBorn() +
dipole()->nDimRadiation();
}
/**
* Return the number of random numbers
* needed to sample the Born process.
*/
int nDimBorn() const {
return bornCXComb()->lastRandomNumbers().size();
}
/**
* Flag, which variables are free variables.
*/
const vector<bool>& sampleFlags();
/**
* Return the support of the splitting kernel.
* The lower bound on the first variable is
* assumed to correspond to the cutoff on the
* evolution variable.
*/
const pair<vector<double>,vector<double> >& support();
/**
* Return the parameter point associated to the splitting
* previously supplied through fixParameters.
*/
const vector<double>& parameterPoint();
/**
* Indicate that presampling of this kernel
* will be performed in the next calls to
* evaluate until stopPresampling() is called.
*/
void startPresampling();
/**
* Indicate that presampling of this kernel
* is done until startPresampling() is called.
*/
void stopPresampling();
/**
* Return true, if currently being presampled
*/
bool presampling() const { return thePresampling; }
/**
* Return the number of points to presample this
* splitting generator.
*/
unsigned long presamplingPoints() const { return thePresamplingPoints; }
/**
* Return the maximum number of trials
* to generate a splitting.
*/
unsigned long maxtry() const { return theMaxTry; }
/**
* Set the number of points to presample this
* splitting generator.
*/
void presamplingPoints(unsigned long p) { thePresamplingPoints = p; }
/**
* Set the maximum number of trials
* to generate a splitting.
*/
void maxtry(unsigned long p) { theMaxTry = p; }
/**
* Evalute the splitting kernel.
*/
double evaluate(const vector<double>&);
/**
* Return the index of the random number corresponding
* to the evolution variable.
*/
int evolutionVariable() const {
return
nDimBorn() +
dipole()->invertedTildeKinematics()->evolutionVariable();
}
/**
* Return the cutoff on the evolution
* random number corresponding to the pt cut.
*/
double evolutionCutoff() const {
return dipole()->invertedTildeKinematics()->evolutionCutoff();
}
public:
/**@name Wrap to the exsample2 interface until this is finally cleaned up. */
//@{
inline const vector<bool>& variable_flags () {
return sampleFlags();
}
inline size_t evolution_variable () const { return evolutionVariable(); }
inline double evolution_cutoff () const {
return evolutionCutoff();
}
inline const vector<double>& parameter_point () {
return parameterPoint();
}
inline void start_presampling () {
startPresampling();
}
inline void stop_presampling () {
stopPresampling();
}
inline size_t dimension () const {
return nDim();
}
inline unsigned long presampling_points () const {
return presamplingPoints();
}
//@}
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 dipole in charge of the emission
*/
Ptr<SubtractionDipole>::ptr theDipole;
/**
* The shower approximation to consider
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* The XComb off which radiation will be generated
*/
StdXCombPtr theBornXComb;
/**
* The XComb describing the process after radiation
*/
StdXCombPtr theRealXComb;
/**
+ * The tilde xcomb objects associated to the real xcomb
+ */
+ vector<StdXCombPtr> theTildeXCombs;
+
+ /**
* True, if currently being presampled
*/
bool thePresampling;
/**
* The number of points to presample this
* splitting generator.
*/
unsigned long thePresamplingPoints;
/**
* The maximum number of trials
* to generate a splitting.
*/
unsigned long theMaxTry;
/**
* The sampling flags
*/
vector<bool> theFlags;
/**
* The support.
*/
pair<vector<double>,vector<double> > theSupport;
/**
* The shower approximation generator.
*/
Ptr<ShowerApproximationGenerator>::tptr theShowerApproximationGenerator;
/**
* The last parameter point
*/
vector<double> theLastParameterPoint;
/**
* The last random numbers used for Born sampling
*/
vector<double> theLastBornPoint;
/**
* Define the Sudakov sampler
*/
typedef
exsample::exponential_generator<ShowerApproximationKernel,UseRandom>
ExponentialGenerator;
/**
* Define a pointer to the Sudakov sampler
*/
typedef
exsample::exponential_generator<ShowerApproximationKernel,UseRandom>*
ExponentialGeneratorPtr;
/**
* The Sudakov sampler
*/
ExponentialGeneratorPtr sampler;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerApproximationKernel & operator=(const ShowerApproximationKernel &);
};
}
#endif /* Herwig_ShowerApproximationKernel_H */

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:43 PM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242504
Default Alt Text
(251 KB)

Event Timeline