Page MenuHomeHEPForge
Files F9501438

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,1692 +1,1692 @@
// -*- 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 "ThePEG/EventRecord/SubProcess.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig++/Utilities/RunDirectories.h"
#include "Herwig++/MatrixElement/ProductionMatrixElement.h"
#include "Herwig++/MatrixElement/HardVertex.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
MatchboxMEBase::MatchboxMEBase()
: MEBase(),
theOneLoop(false),
theOneLoopNoBorn(false),
theOneLoopNoLoops(false),
theNoCorrelations(false),
theHavePDFs(false,false), checkedPDFs(false),
theDiagramWeightVerboseDown(10000000000000.),
theDiagramWeightVerboseUp(0.) {}
MatchboxMEBase::~MatchboxMEBase() {}
Ptr<MatchboxFactory>::tptr MatchboxMEBase::factory() const { return theFactory; }
void MatchboxMEBase::factory(Ptr<MatchboxFactory>::tptr f) { theFactory = f; }
Ptr<Tree2toNGenerator>::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); }
Ptr<ProcessData>::tptr MatchboxMEBase::processData() const { return factory()->processData(); }
unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); }
vector<int> MatchboxMEBase::getNLightJetVec() const { return factory()->nLightJetVec(); }
vector<int> MatchboxMEBase::getNHeavyJetVec() const { return factory()->nHeavyJetVec(); }
vector<int> MatchboxMEBase::getNLightProtonVec() const { return factory()->nLightProtonVec(); }
double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); }
double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); }
bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); }
bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); }
bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); }
bool MatchboxMEBase::verbose() const { return factory()->verbose(); }
bool MatchboxMEBase::initVerbose() const { return factory()->initVerbose(); }
void MatchboxMEBase::getDiagrams() const {
if ( diagramGenerator() && processData() ) {
vector<Ptr<Tree2toNDiagram>::ptr> diags;
vector<Ptr<Tree2toNDiagram>::ptr>& res =
processData()->diagramMap()[subProcess().legs];
if ( res.empty() ) {
res = diagramGenerator()->generate(subProcess().legs,orderInAlphaS(),orderInAlphaEW());
}
copy(res.begin(),res.end(),back_inserter(diags));
processData()->fillMassGenerators(subProcess().legs);
if ( diags.empty() )
return;
for ( vector<Ptr<Tree2toNDiagram>::ptr>::iterator d = diags.begin();
d != diags.end(); ++d ) {
add(*d);
}
return;
}
throw Exception()
<< "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
<< "Please check your setup." << Exception::abortnow;
}
Selector<MEBase::DiagramIndex>
MatchboxMEBase::diagrams(const DiagramVector & diags) const {
if ( phasespace() ) {
return phasespace()->selectDiagrams(diags);
}
throw Exception()
<< "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return Selector<MEBase::DiagramIndex>();
}
Selector<const ColourLines *>
MatchboxMEBase::colourGeometries(tcDiagPtr diag) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->haveColourFlows() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
return matchboxAmplitude()->colourGeometries(diag);
}
}
Ptr<Tree2toNDiagram>::tcptr tdiag =
dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>(diag);
assert(diag && processData());
vector<ColourLines*>& flows = processData()->colourFlowMap()[tdiag];
if ( flows.empty() ) {
list<list<list<pair<int,bool> > > > cflows =
ColourBasis::colourFlows(tdiag);
for ( list<list<list<pair<int,bool> > > >::const_iterator fit =
cflows.begin(); fit != cflows.end(); ++fit ) {
flows.push_back(new ColourLines(ColourBasis::cfstring(*fit)));
}
}
Selector<const ColourLines *> res;
for ( vector<ColourLines*>::const_iterator f = flows.begin();
f != flows.end(); ++f )
res.insert(1.0,*f);
return res;
}
void MatchboxMEBase::constructVertex(tSubProPtr sub, const ColourLines* cl) {
if ( !canFillRhoMatrix() || !factory()->spinCorrelations() )
return;
assert(matchboxAmplitude());
assert(matchboxAmplitude()->colourBasis());
// get the colour structure for the selected colour flow
size_t cStructure =
matchboxAmplitude()->colourBasis()->tensorIdFromFlow(lastXComb().lastDiagram(),cl);
// hard process for processing the spin info
tPVector hard;
hard.push_back(sub->incoming().first);
hard.push_back(sub->incoming().second);
vector<PDT::Spin> out;
for ( size_t k = 0; k < sub->outgoing().size(); ++k ) {
out.push_back(sub->outgoing()[k]->data().iSpin());
hard.push_back(sub->outgoing()[k]);
}
// calculate dummy wave functions to fill the spin info
static vector<VectorWaveFunction> dummyPolarizations;
static vector<SpinorWaveFunction> dummySpinors;
static vector<SpinorBarWaveFunction> dummyBarSpinors;
for ( size_t k = 0; k < hard.size(); ++k ) {
if ( hard[k]->data().iSpin() == PDT::Spin1Half ) {
if ( hard[k]->id() > 0 && k > 1 ) {
SpinorBarWaveFunction(dummyBarSpinors,hard[k],
outgoing, true);
} else if ( hard[k]->id() < 0 && k > 1 ) {
SpinorWaveFunction(dummySpinors,hard[k],
outgoing, true);
} else if ( hard[k]->id() > 0 && k < 2 ) {
SpinorWaveFunction(dummySpinors,hard[k],
incoming, false);
} else if ( hard[k]->id() < 0 && k < 2 ) {
SpinorBarWaveFunction(dummyBarSpinors,hard[k],
incoming, false);
}
} else if ( hard[k]->data().iSpin() == PDT::Spin1 ) {
VectorWaveFunction(dummyPolarizations,hard[k],
k > 1 ? outgoing : incoming,
k > 1 ? true : false,
hard[k]->data().hardProcessMass() == ZERO);
} else assert(false);
}
// fill the production matrix element
ProductionMatrixElement pMe(mePartonData()[0]->iSpin(),
mePartonData()[1]->iSpin(),
out);
for ( map<vector<int>,CVector>::const_iterator lamp = lastLargeNAmplitudes().begin();
lamp != lastLargeNAmplitudes().end(); ++lamp ) {
vector<unsigned int> pMeHelicities
= matchboxAmplitude()->physicalHelicities(lamp->first);
pMe(pMeHelicities) = lamp->second[cStructure];
}
// set the spin information
HardVertexPtr hardvertex = new_ptr(HardVertex());
hardvertex->ME(pMe);
if ( sub->incoming().first->spinInfo() )
sub->incoming().first->spinInfo()->productionVertex(hardvertex);
if ( sub->incoming().second->spinInfo() )
sub->incoming().second->spinInfo()->productionVertex(hardvertex);
for ( ParticleVector::const_iterator p = sub->outgoing().begin();
p != sub->outgoing().end(); ++p ) {
if ( (**p).spinInfo() )
(**p).spinInfo()->productionVertex(hardvertex);
}
}
unsigned int MatchboxMEBase::orderInAlphaS() const {
return subProcess().orderInAlphaS;
}
unsigned int MatchboxMEBase::orderInAlphaEW() const {
return subProcess().orderInAlphaEW;
}
void MatchboxMEBase::setXComb(tStdXCombPtr xc) {
MEBase::setXComb(xc);
lastMatchboxXComb(xc);
if ( phasespace() )
phasespace()->setXComb(xc);
if ( scaleChoice() )
scaleChoice()->setXComb(xc);
if ( matchboxAmplitude() )
matchboxAmplitude()->setXComb(xc);
}
double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) {
// shamelessly stolen from PartonExtractor.cc
Energy2 shmax = lastCuts().sHatMax();
Energy2 shmin = lastCuts().sHatMin();
Energy2 sh = shmin*pow(shmax/shmin, *r1);
double ymax = lastCuts().yHatMax();
double ymin = lastCuts().yHatMin();
double km = log(shmax/shmin);
ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));
double y = ymin + (*r2)*(ymax - ymin);
double x1 = exp(-0.5*log(lastS()/sh) + y);
double x2 = exp(-0.5*log(lastS()/sh) - y);
Lorentz5Momentum P1 = lastParticles().first->momentum();
LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy());
p1.rotateY(P1.theta());
p1.rotateZ(P1.phi());
meMomenta()[0] = p1;
Lorentz5Momentum P2 = lastParticles().second->momentum();
LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy());
p2.rotateY(P2.theta());
p2.rotateZ(P2.phi());
meMomenta()[1] = p2;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
return km*(ymax - ymin);
}
bool MatchboxMEBase::generateKinematics(const double * r) {
if ( phasespace() ) {
jacobian(phasespace()->generateKinematics(r,meMomenta()));
if ( jacobian() == 0.0 )
return false;
setScale();
logGenerateKinematics(r);
assert(lastMatchboxXComb());
if ( nDimAmplitude() > 0 ) {
amplitudeRandomNumbers().resize(nDimAmplitude());
copy(r + nDimPhasespace(),
r + nDimPhasespace() + nDimAmplitude(),
amplitudeRandomNumbers().begin());
}
if ( nDimInsertions() > 0 ) {
insertionRandomNumbers().resize(nDimInsertions());
copy(r + nDimPhasespace() + nDimAmplitude(),
r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(),
insertionRandomNumbers().begin());
}
return true;
}
throw Exception()
<< "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return false;
}
int MatchboxMEBase::nDim() const {
if ( lastMatchboxXComb() )
return nDimPhasespace() + nDimAmplitude() + nDimInsertions();
int ampAdd = 0;
if ( matchboxAmplitude() ) {
ampAdd = matchboxAmplitude()->nDimAdditional();
}
int insertionAdd = 0;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
insertionAdd = max(insertionAdd,(**v).nDimAdditional());
}
return nDimBorn() + ampAdd + insertionAdd;
}
int MatchboxMEBase::nDimBorn() const {
if ( lastMatchboxXComb() )
return nDimPhasespace();
if ( phasespace() )
return phasespace()->nDim(diagrams().front()->partons());
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 fcscale = factorizationScale();
Energy2 fscale = fcscale*sqr(factorizationScaleFactor());
Energy2 rscale = renormalizationScale()*sqr(renormalizationScaleFactor());
Energy2 ewrscale = renormalizationScaleQED();
lastXCombPtr()->lastScale(fscale);
lastXCombPtr()->lastCentralScale(fcscale);
lastMatchboxXComb()->lastRenormalizationScale(rscale);
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().alphaEMME(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 {}
bool MatchboxMEBase::havePDFWeight1() const {
if ( checkedPDFs )
return theHavePDFs.first;
theHavePDFs.first =
factory()->isIncoming(mePartonData()[0]) &&
lastXCombPtr()->partonBins().first->pdf();
theHavePDFs.second =
factory()->isIncoming(mePartonData()[1]) &&
lastXCombPtr()->partonBins().second->pdf();
checkedPDFs = true;
return theHavePDFs.first;
}
bool MatchboxMEBase::havePDFWeight2() const {
if ( checkedPDFs )
return theHavePDFs.second;
theHavePDFs.first =
factory()->isIncoming(mePartonData()[0]) &&
lastXCombPtr()->partonBins().first->pdf();
theHavePDFs.second =
factory()->isIncoming(mePartonData()[1]) &&
lastXCombPtr()->partonBins().second->pdf();
checkedPDFs = true;
return theHavePDFs.second;
}
void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const {
if ( !havePDFWeight1() && !havePDFWeight2() ) {
lastMEPDFWeight(1.0);
logPDFWeight();
return;
}
double w = 1.;
if ( havePDFWeight1() )
w *= pdf1(factorizationScale);
if ( havePDFWeight2() )
w *= pdf2(factorizationScale);
lastMEPDFWeight(w);
logPDFWeight();
}
-double MatchboxMEBase::pdf1(Energy2 fscale, double xEx) const {
+double MatchboxMEBase::pdf1(Energy2 fscale, double xEx, double xFactor) const {
assert(lastXCombPtr()->partonBins().first->pdf());
- if ( xEx < 1. && lastX1() >= xEx ) {
+ if ( xEx < 1. && lastX1()*xFactor >= xEx ) {
return
- ( ( 1. - lastX1() ) / ( 1. - xEx ) ) *
+ ( ( 1. - lastX1()*xFactor ) / ( 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();
+ lastX1()*xFactor)/lastX1()/xFactor;
}
-double MatchboxMEBase::pdf2(Energy2 fscale, double xEx) const {
+double MatchboxMEBase::pdf2(Energy2 fscale, double xEx, double xFactor) const {
assert(lastXCombPtr()->partonBins().second->pdf());
- if ( xEx < 1. && lastX2() >= xEx ) {
+ if ( xEx < 1. && lastX2()*xFactor >= xEx ) {
return
- ( ( 1. - lastX2() ) / ( 1. - xEx ) ) *
+ ( ( 1. - lastX2()*xFactor ) / ( 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();
+ lastX2()*xFactor)/lastX2()/xFactor;
}
double MatchboxMEBase::me2() const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
double res =
matchboxAmplitude()->me2()*
me2Norm();
return res;
}
throw Exception()
<< "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::largeNME2(Ptr<ColourBasis>::tptr largeNBasis) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() ) {
largeNBasis->prepare(mePartonData(),false);
matchboxAmplitude()->prepareAmplitudes(this);
}
double res =
matchboxAmplitude()->largeNME2(largeNBasis)*
me2Norm();
return res;
}
throw Exception()
<< "MatchboxMEBase::largeNME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::finalStateSymmetry() const {
if ( symmetryFactor() > 0.0 )
return symmetryFactor();
double sFactor = 1.;
map<long,int> counts;
cPDVector checkData;
copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData));
cPDVector::iterator p = checkData.begin();
while ( !checkData.empty() ) {
if ( counts.find((**p).id()) != counts.end() ) {
counts[(**p).id()] += 1;
} else {
counts[(**p).id()] = 1;
}
checkData.erase(p);
p = checkData.begin();
continue;
}
for ( map<long,int>::const_iterator c = counts.begin();
c != counts.end(); ++c ) {
if ( c->second == 1 )
continue;
if ( c->second == 2 )
sFactor /= 2.;
else if ( c->second == 3 )
sFactor /= 6.;
else if ( c->second == 4 )
sFactor /= 24.;
}
symmetryFactor(sFactor);
return symmetryFactor();
}
double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const {
// assume that we always have incoming
// spin-1/2 or massless spin-1 particles
double fac = 1./4.;
if ( hasInitialAverage() )
fac = 1.;
double couplings = 1.0;
if ( (orderInAlphaS() > 0 || addAlphaS != 0) && !hasRunningAlphaS() ) {
fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS));
couplings *= pow(lastAlphaS(),double(orderInAlphaS()+addAlphaS));
}
if ( orderInAlphaEW() > 0 && !hasRunningAlphaEW() ) {
fac *= pow(lastAlphaEM()/SM().alphaEMMZ(),double(orderInAlphaEW()));
couplings *= pow(lastAlphaEM(),double(orderInAlphaEW()));
}
lastMECouplings(couplings);
if ( !hasInitialAverage() ) {
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 !hasFinalStateSymmetry() ? finalStateSymmetry()*fac : fac;
}
CrossSection MatchboxMEBase::dSigHatDR() const {
getPDFWeight();
if ( !lastXCombPtr()->willPassCuts() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double xme2 = me2();
lastME2(xme2);
if (factory()->verboseDia()){
double diagweightsum = 0.0;
for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
d != diagrams().end(); ++d ) {
diagweightsum += phasespace()->diagramWeight(dynamic_cast<const Tree2toNDiagram&>(**d));
}
double piWeight = pow(2.*Constants::pi,(int)(3*(meMomenta().size()-2)-4));
double units = pow(lastSHat() / GeV2, mePartonData().size() - 4.);
bookMEoverDiaWeight(log(xme2/(diagweightsum*piWeight*units)));//
}
if ( xme2 == 0. && !oneLoopNoBorn() ) {
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double vme2 = 0.;
if ( oneLoop() && !oneLoopNoLoops() )
vme2 = oneLoopInterference();
CrossSection res = ZERO;
if ( !oneLoopNoBorn() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * xme2;
if ( oneLoop() && !oneLoopNoLoops() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * vme2;
if ( !onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).setXComb(lastXCombPtr());
res += (**v).dSigHatDR();
}
if ( checkPoles() )
logPoles();
}
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
lastMECrossSection(res);
return lastMECrossSection();
}
double MatchboxMEBase::oneLoopInterference() const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->oneLoopAmplitudes() )
matchboxAmplitude()->prepareOneLoopAmplitudes(this);
double res =
matchboxAmplitude()->oneLoopInterference()*
me2Norm(1);
return res;
}
throw Exception()
<< "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
MatchboxMEBase::AccuracyHistogram::AccuracyHistogram(double low,
double up,
unsigned int nbins)
: lower(low), upper(up),
sameSign(0), oppositeSign(0), nans(0),
overflow(0), underflow(0) {
double step = (up-low)/nbins;
for ( unsigned int k = 1; k <= nbins; ++k )
bins[lower + k*step] = 0.0;
}
void MatchboxMEBase::AccuracyHistogram::book(double a, double b) {
if ( isnan(a) || isnan(b) ||
isinf(a) || isinf(b) ) {
++nans;
return;
}
if ( a*b >= 0. )
++sameSign;
if ( a*b < 0. )
++oppositeSign;
double r = 1.;
if ( abs(a) != 0.0 )
r = abs(1.-abs(b/a));
else if ( abs(b) != 0.0 )
r = abs(b);
if ( log10(r) < lower || r == 0.0 ) {
++underflow;
return;
}
if ( log10(r) > upper ) {
++overflow;
return;
}
map<double,double>::iterator bin =
bins.upper_bound(log10(r));
if ( bin == bins.end() )
return;
bin->second += 1.;
}
void MatchboxMEBase::AccuracyHistogram::dump(const std::string& folder, const std::string& prefix,
const cPDVector& proc) const {
ostringstream fname("");
for ( cPDVector::const_iterator p = proc.begin();
p != proc.end(); ++p )
fname << (**p).PDGName();
ofstream out((folder+"/"+prefix+fname.str()+".dat").c_str());
out << "# same sign : " << sameSign << " opposite sign : "
<< oppositeSign << " nans : " << nans
<< " overflow : " << overflow
<< " underflow : " << underflow << "\n";
for ( map<double,double>::const_iterator b = bins.begin();
b != bins.end(); ++b ) {
map<double,double>::const_iterator bp = b; --bp;
if ( b->second != 0. ) {
if ( b != bins.begin() )
out << bp->first;
else
out << lower;
out << " " << b->first
<< " " << b->second
<< "\n" << flush;
}
}
ofstream gpout((folder+"/"+prefix+fname.str()+".gp").c_str());
gpout << "set terminal png\n"
<< "set xlabel 'accuracy of pole cancellation [decimal places]'\n"
<< "set ylabel 'counts\n"
<< "set xrange [-20:0]\n"
<< "set output '" << prefix << fname.str() << ".png'\n"
<< "plot '" << prefix << fname.str() << ".dat' using (0.5*($1+$2)):3 with linespoints pt 7 ps 1 not";
}
void MatchboxMEBase::AccuracyHistogram::persistentOutput(PersistentOStream& os) const {
os << lower << upper << bins
<< sameSign << oppositeSign << nans
<< overflow << underflow;
}
void MatchboxMEBase::AccuracyHistogram::persistentInput(PersistentIStream& is) {
is >> lower >> upper >> bins
>> sameSign >> oppositeSign >> nans
>> overflow >> underflow;
}
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();
}
if (res2me != 0.0 || res2i != 0.0) epsilonSquarePoleHistograms[mePartonData()].book(res2me,res2i);
if (res1me != 0.0 || res1i != 0.0) epsilonPoleHistograms[mePartonData()].book(res1me,res1i);
}
bool MatchboxMEBase::haveOneLoop() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->haveOneLoop();
return false;
}
bool MatchboxMEBase::onlyOneLoop() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->onlyOneLoop();
return false;
}
bool MatchboxMEBase::isDRbar() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isDRbar();
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()*
me2Norm(1);
}
return 0.;
}
double MatchboxMEBase::oneLoopSinglePole() const {
if ( matchboxAmplitude() ) {
return
matchboxAmplitude()->oneLoopSinglePole()*
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 ) {
list<Ptr<SubtractionDipole>::ptr> matchDipoles;
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
dipoles.begin(); d != dipoles.end(); ++d ) {
if ( !(**d).canHandleEmitter(rep,emitter) )
continue;
matchDipoles.push_back(*d);
}
if ( matchDipoles.empty() )
continue;
for ( int emission = 2; emission < nreal; ++emission ) {
if ( emission == emitter )
continue;
list<Ptr<SubtractionDipole>::ptr> matchDipoles2;
for ( list<Ptr<SubtractionDipole>::ptr>::const_iterator d =
matchDipoles.begin(); d != matchDipoles.end(); ++d ) {
if ( !(**d).canHandleSplitting(rep,emitter,emission) )
continue;
matchDipoles2.push_back(*d);
}
if ( matchDipoles2.empty() )
continue;
map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo> mergeInfo;
for ( DiagramVector::const_iterator d = diagrams().begin(); d != diagrams().end(); ++d ) {
Ptr<Tree2toNDiagram>::ptr check =
new_ptr(Tree2toNDiagram(*dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d)));
map<int,int> theMergeLegs;
for ( unsigned int i = 0; i < check->external().size(); ++i )
theMergeLegs[i] = -1;
int theEmitter = check->mergeEmission(emitter,emission,theMergeLegs);
// no underlying Born
if ( theEmitter == -1 )
continue;
SubtractionDipole::MergeInfo info;
info.diagram = check;
info.emitter = theEmitter;
info.mergeLegs = theMergeLegs;
mergeInfo[*d] = info;
}
if ( mergeInfo.empty() )
continue;
for ( int spectator = 0; spectator < nreal; ++spectator ) {
if ( spectator == emitter || spectator == emission )
continue;
list<Ptr<SubtractionDipole>::ptr> matchDipoles3;
for ( list<Ptr<SubtractionDipole>::ptr>::const_iterator d =
matchDipoles2.begin(); d != matchDipoles2.end(); ++d ) {
if ( !(**d).canHandleSpectator(rep,spectator) )
continue;
matchDipoles3.push_back(*d);
}
if ( matchDipoles3.empty() )
continue;
if ( noDipole(emitter,emission,spectator) )
continue;
for ( list<Ptr<SubtractionDipole>::ptr>::const_iterator d =
matchDipoles3.begin(); d != matchDipoles3.end(); ++d ) {
if ( !(**d).canHandle(rep,emitter,emission,spectator) )
continue;
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator b =
borns.begin(); b != borns.end(); ++b ) {
if ( (**b).onlyOneLoop() )
continue;
if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(*b,*d)))
!= done.end() )
continue;
// now get to work
(**d).clearBookkeeping();
(**d).factory(factory());
(**d).realEmitter(emitter);
(**d).realEmission(emission);
(**d).realSpectator(spectator);
(**d).realEmissionME(const_cast<MatchboxMEBase*>(this));
(**d).underlyingBornME(*b);
(**d).setupBookkeeping(mergeInfo);
if ( !((**d).empty()) ) {
res.push_back((**d).cloneMe());
Ptr<SubtractionDipole>::tptr nDipole = res.back();
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() << "MatchboxMEBase::getDipoles(): Dipole " << dname.str() << " already existing.";
if ( !factory()->reweighters().empty() ) {
for ( vector<ReweightPtr>::const_iterator rw = factory()->reweighters().begin();
rw != factory()->reweighters().end(); ++rw )
nDipole->addReweighter(*rw);
}
if ( !factory()->preweighters().empty() ) {
for ( vector<ReweightPtr>::const_iterator rw = factory()->preweighters().begin();
rw != factory()->preweighters().end(); ++rw )
nDipole->addPreweighter(*rw);
}
nDipole->cloneDependencies(dname.str());
}
}
}
}
}
}
vector<Ptr<SubtractionDipole>::tptr> partners;
copy(res.begin(),res.end(),back_inserter(partners));
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = res.begin();
d != res.end(); ++d )
(**d).partnerDipoles(partners);
return res;
}
double MatchboxMEBase::colourCorrelatedME2(pair<int,int> ij) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
double res =
matchboxAmplitude()->colourCorrelatedME2(ij)*
me2Norm();
return res;
}
throw Exception()
<< "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr largeNBasis) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() ) {
largeNBasis->prepare(mePartonData(),false);
matchboxAmplitude()->prepareAmplitudes(this);
}
double res =
matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)*
me2Norm();
return res;
}
throw Exception()
<< "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
double res =
matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
me2Norm();
return res;
}
throw Exception()
<< "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::spinCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
double res =
matchboxAmplitude()->spinCorrelatedME2(ij,c)*
me2Norm();
return res;
}
throw Exception()
<< "MatchboxMEBase::spinCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
void MatchboxMEBase::flushCaches() {
MEBase::flushCaches();
if ( matchboxAmplitude() )
matchboxAmplitude()->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).flushCaches();
}
}
void MatchboxMEBase::print(ostream& os) const {
os << "--- MatchboxMEBase setup -------------------------------------------------------\n";
os << " '" << name() << "' for subprocess:\n";
os << " ";
for ( PDVector::const_iterator pp = subProcess().legs.begin();
pp != subProcess().legs.end(); ++pp ) {
os << (**pp).PDGName() << " ";
if ( pp == subProcess().legs.begin() + 1 )
os << "-> ";
}
os << "\n";
os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections";
if ( oneLoopNoBorn() )
os << " without Born contributions";
if ( oneLoopNoLoops() )
os << " without loop contributions";
os << "\n";
if ( oneLoop() && !onlyOneLoop() ) {
os << " using insertion operators\n";
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
os << " '" << (**v).name() << "' with "
<< ((**v).isDR() ? "" : "C") << "DR/";
if ( (**v).isCS() )
os << "CS";
if ( (**v).isBDK() )
os << "BDK";
if ( (**v).isExpanded() )
os << "expanded";
os << " conventions\n";
}
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxMEBase::printLastEvent(ostream& os) const {
os << "--- MatchboxMEBase last event information --------------------------------------\n";
os << " for matrix element '" << name() << "'\n";
os << " process considered:\n ";
int in = 0;
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
if ( ++in == 2 )
os << " -> ";
}
os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n"
<< " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2)
<< " x1 = " << lastX1() << " x2 = " << lastX2()
<< " alphaS = " << lastAlphaS() << "\n";
os << " momenta/GeV generated from random numbers\n ";
copy(lastXComb().lastRandomNumbers().begin(),
lastXComb().lastRandomNumbers().end(),ostream_iterator<double>(os," "));
os << ":\n ";
for ( vector<Lorentz5Momentum>::const_iterator p = meMomenta().begin();
p != meMomenta().end(); ++p ) {
os << (*p/GeV) << "\n ";
}
os << "last cross section/nb calculated was:\n "
<< (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n";
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxMEBase::logGenerateKinematics(const double * r) const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' generated kinematics\nfrom "
<< nDim() << " random numbers:\n";
copy(r,r+nDim(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "storing phase space information in XComb "
<< lastXCombPtr() << "\n";
generator()->log() << "generated phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
cPDVector::const_iterator dit = mePartonData().begin();
for ( ; pit != meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
<< "and Jacobian = " << jacobian() << " sHat/GeV2 = "
<< (lastSHat()/GeV2) << "\n" << flush;
}
void MatchboxMEBase::logSetScale() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n"
<< "scale/GeV2 = " << (scale()/GeV2) << " xi_R = "
<< renormalizationScaleFactor() << " xi_F = "
<< factorizationScaleFactor() << "\n"
<< "alpha_s = " << lastAlphaS() << "\n" << flush;
}
void MatchboxMEBase::logPDFWeight() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' calculated pdf weight = "
<< lastMEPDFWeight() << " from XComb "
<< lastXCombPtr() << "\n"
<< "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) "
<< "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n"
<< flush;
}
void MatchboxMEBase::logME2() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' evaluated me2 using XComb "
<< lastXCombPtr() << "\n"
<< "and phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
cPDVector::const_iterator dit = mePartonData().begin();
for ( ; pit != meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
<< "sHat/GeV2 = " << (lastSHat()/GeV2)
<< " me2 = " << lastME2() << "\n" << flush;
}
void MatchboxMEBase::logDSigHatDR() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' evaluated cross section using XComb "
<< lastXCombPtr() << "\n"
<< "Jacobian = " << jacobian() << " sHat/GeV2 = "
<< (lastSHat()/GeV2) << " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void MatchboxMEBase::cloneDependencies(const std::string& prefix) {
if ( phasespace() ) {
Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
throw InitException() << "MatchboxMEBase::cloneDependencies(): 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() << "MatchboxMEBase::cloneDependencies(): Amplitude " << pname.str() << " already existing.";
myAmplitude->cloneDependencies(pname.str());
matchboxAmplitude(myAmplitude);
amplitude(myAmplitude);
matchboxAmplitude()->orderInGs(orderInAlphaS());
matchboxAmplitude()->orderInGem(orderInAlphaEW());
}
if ( scaleChoice() ) {
Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
throw InitException() << "MatchboxMEBase::cloneDependencies(): 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() << "MatchboxMEBase::cloneDependencies(): 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() << "MatchboxMEBase::cloneDependencies(): Insertion operator " << pname.str() << " already existing.";
*v = myIOP;
}
}
void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const {
// fixme We need to pass on the partons from the xcmob here, not
// assuming one subprocess per matrix element
if ( phasespace() )
xc.nDimPhasespace(phasespace()->nDim(diagrams().front()->partons()));
if ( matchboxAmplitude() ) {
xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional());
if ( matchboxAmplitude()->colourBasis() ) {
size_t cdim =
matchboxAmplitude()->colourBasis()->prepare(diagrams(),noCorrelations());
xc.colourBasisDim(cdim);
}
if ( matchboxAmplitude()->isExternal() ) {
xc.externalId(matchboxAmplitude()->externalId(diagrams().front()->partons()));
}
}
int insertionAdd = 0;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
insertionAdd = max(insertionAdd,(**v).nDimAdditional());
}
xc.nDimInsertions(insertionAdd);
xc.nLight(getNLight());
for (size_t inlv=0; inlv<getNLightJetVec().size(); ++inlv)
xc.nLightJetVec(getNLightJetVec()[inlv]);
for (size_t inhv=0; inhv<getNHeavyJetVec().size(); ++inhv)
xc.nHeavyJetVec(getNHeavyJetVec()[inhv]);
for (size_t inlpv=0; inlpv<getNLightProtonVec().size(); ++inlpv)
xc.nLightProtonVec(getNLightProtonVec()[inlpv]);
xc.olpId(olpProcess());
if ( initVerbose() ) {
string fname = name() + ".diagrams";
ifstream test(fname.c_str());
if ( !test ) {
test.close();
ofstream out(fname.c_str());
for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
d != diagrams().end(); ++d ) {
DiagramDrawer::drawDiag(out,dynamic_cast<const Tree2toNDiagram&>(**d));
out << "\n";
}
}
}
}
StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec&,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
Ptr<MatchboxXComb>::ptr xc =
new_ptr(MatchboxXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts, newME,
newDiagrams, mir,
newHead));
prepareXComb(*xc);
return xc;
}
StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
Ptr<MatchboxXComb>::ptr xc =
new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams));
prepareXComb(*xc);
return xc;
}
void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theFactory << thePhasespace
<< theAmplitude << theScaleChoice << theVirtuals
<< theReweights << theSubprocess << theOneLoop
<< theOneLoopNoBorn << theOneLoopNoLoops
<< epsilonSquarePoleHistograms << epsilonPoleHistograms
<< theOLPProcess << theNoCorrelations
<< theHavePDFs << checkedPDFs<<theDiagramWeightVerboseDown<<theDiagramWeightVerboseUp;
}
void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theFactory >> thePhasespace
>> theAmplitude >> theScaleChoice >> theVirtuals
>> theReweights >> theSubprocess >> theOneLoop
>> theOneLoopNoBorn >> theOneLoopNoLoops
>> epsilonSquarePoleHistograms >> epsilonPoleHistograms
>> theOLPProcess >> theNoCorrelations
>> theHavePDFs >> checkedPDFs>>theDiagramWeightVerboseDown>>theDiagramWeightVerboseUp;
lastMatchboxXComb(theLastXComb);
}
void MatchboxMEBase::Init() {
static ClassDocumentation<MatchboxMEBase> documentation
("MatchboxMEBase is the base class for matrix elements "
"in the context of the matchbox NLO interface.");
}
IBPtr MatchboxMEBase::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxMEBase::fullclone() const {
return new_ptr(*this);
}
void MatchboxMEBase::doinit() {
MEBase::doinit();
if ( !theAmplitude )
theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
if ( matchboxAmplitude() )
matchboxAmplitude()->init();
if ( phasespace() ) {
phasespace()->init();
matchboxAmplitude()->checkReshuffling(phasespace());
}
if ( scaleChoice() ) {
scaleChoice()->init();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).init();
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).init();
}
}
void MatchboxMEBase::bookMEoverDiaWeight(double x) const {
if (MEoverDiaWeight.size()==0){
theDiagramWeightVerboseDown=min(theDiagramWeightVerboseDown,x*0.9);
theDiagramWeightVerboseUp=max(theDiagramWeightVerboseUp,x*1.1);
}
map<double,double>::iterator bx =MEoverDiaWeight.upper_bound(x);
if ( bx == MEoverDiaWeight.end() ) {
return;
}
bx->second += 1.;
Nevents++;
if (int(Nevents)%1000==0){
ofstream out((RunDirectories::runStorage()+"/"+name()+"-MeoDiaW.dat").c_str());
int i=0;
double m=0.;
for ( map<double,double>::const_iterator bx = MEoverDiaWeight.begin();bx != MEoverDiaWeight.end(); ++bx,i++ ) {
out << " " << bx->first<<" "<<( bx->second/double(Nevents))<<"\n ";
m=max(m,bx->second/double(Nevents));
}
out.close();
ofstream gpout((RunDirectories::runStorage()+"/"+name()+"-MeoDiaW.gp").c_str());
gpout << "set terminal epslatex color solid\n"
<< "set output '" << name()<<"-MeoDiaW"<< "-plot.tex'\n"
<< "#set logscale x\n"
<< "set xrange [" << theDiagramWeightVerboseDown << ":" << theDiagramWeightVerboseUp << "]\n"
<< "set yrange [0.:"<<(m*0.95)<<"]\n"
<< "set xlabel '$log(ME/\\sum DiaW)$'\n"
<< "set size 0.7,0.7\n"
<< "plot 1 w lines lc rgbcolor \"#DDDDDD\" notitle, '" << name()<<"-MeoDiaW"
<< ".dat' with histeps lc rgbcolor \"#00AACC\" t '$"<<name()<<"$'";
gpout.close();
}
}
void MatchboxMEBase::doinitrun() {
MEBase::doinitrun();
if ( matchboxAmplitude() )
matchboxAmplitude()->initrun();
if ( phasespace() ) {
phasespace()->initrun();
}
if ( scaleChoice() ) {
scaleChoice()->initrun();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).initrun();
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).initrun();
}
if ( factory()->verboseDia() ) {
for ( int k = 0; k < factory()->diagramWeightVerboseNBins() ; ++k ) {
MEoverDiaWeight[theDiagramWeightVerboseDown+
double(k)*(theDiagramWeightVerboseUp-
theDiagramWeightVerboseDown)
/double(factory()->diagramWeightVerboseNBins()) ] = 0.;
}
Nevents=0.;
ofstream out("DiagramWeights.sh");
out<<"P=$(pwd)"
<<"\ncd "<<RunDirectories::runStorage()
<<"\nrm -f DiagramWeights.tex"
<<"\n echo \"\\documentclass{article}\" >> DiagramWeights.tex"
<<"\n echo \"\\usepackage{amsmath,amsfonts,amssymb,graphicx,color}\" >> DiagramWeights.tex"
<<"\n echo \"\\usepackage[left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}\" >> DiagramWeights.tex"
<<"\n echo \"\\begin{document}\" >> DiagramWeights.tex"
<<"\n echo \"\\setlength{\\parindent}{0cm}\" >> DiagramWeights.tex"
<<"\n\n for i in $(ls *.gp | sed s/'\\.gp'//g) ; "
<<"\n do"
<<"\n echo \"\\input{\"\"$i\"-plot\"}\" >> DiagramWeights.tex"
<<"\n done"
<<"\n echo \"\\end{document}\" >> DiagramWeights.tex "
<<"\n for i in *.gp ; do "
<<"\n gnuplot $i "
<<"\n done "
<<"\n pdflatex DiagramWeights.tex \ncp DiagramWeights.pdf $P";
out.close();
}
}
void MatchboxMEBase::dofinish() {
MEBase::dofinish();
for ( map<cPDVector,AccuracyHistogram>::const_iterator
b = epsilonSquarePoleHistograms.begin();
b != epsilonSquarePoleHistograms.end(); ++b ) {
b->second.dump(factory()->poleData(),"epsilonSquarePoles-",b->first);
}
for ( map<cPDVector,AccuracyHistogram>::const_iterator
b = epsilonPoleHistograms.begin();
b != epsilonPoleHistograms.end(); ++b ) {
b->second.dump(factory()->poleData(),"epsilonPoles-",b->first);
}
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MatchboxMEBase,MEBase>
describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "Herwig.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,1137 +1,1137 @@
// -*- C++ -*-
//
// MatchboxMEBase.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MatchboxMEBase_H
#define HERWIG_MatchboxMEBase_H
//
// This is the declaration of the MatchboxMEBase class.
//
#include "ThePEG/MatrixElement/MEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/ProcessData.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxReweightBase.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig++/MatrixElement/Matchbox/InsertionOperators/MatchboxInsertionOperator.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig++/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXComb.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxMEBase is the base class for matrix elements
* in the context of the matchbox NLO interface.
*
* @see \ref MatchboxMEBaseInterfaces "The interfaces"
* defined for MatchboxMEBase.
*/
class MatchboxMEBase:
public MEBase, public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxMEBase();
/**
* The destructor.
*/
virtual ~MatchboxMEBase();
//@}
public:
/**
* Return the factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tptr factory() const;
/**
* Set the factory which produced this matrix element
*/
void factory(Ptr<MatchboxFactory>::tptr f);
/** @name Subprocess and diagram information. */
//@{
/**
* Return the subprocess.
*/
const Process& subProcess() const { return theSubprocess; }
/**
* Access the subprocess.
*/
Process& subProcess() { return theSubprocess; }
/**
* Return the diagram generator.
*/
Ptr<Tree2toNGenerator>::tptr diagramGenerator() const;
/**
* Return the process data.
*/
Ptr<ProcessData>::tptr processData() const;
/**
* Return true, if this matrix element does not want to
* make use of mirroring processes; in this case all
* possible partonic subprocesses with a fixed assignment
* of incoming particles need to be provided through the diagrams
* added with the add(...) method.
*/
virtual bool noMirror () const { return true; }
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
using MEBase::getDiagrams;
/**
* With the information previously supplied with the
* setKinematics(...) method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector &) const;
using MEBase::diagrams;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Return true, if this amplitude is capable of consistently filling
* the rho matrices for the spin correllations
*/
virtual bool canFillRhoMatrix() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->canFillRhoMatrix();
return false;
}
/**
* construct the spin information for the interaction
*/
virtual void constructVertex(tSubProPtr) {}
/**
* construct the spin information for the interaction
*/
virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);
/**
* 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 true, if this amplitude already includes averaging over
* incoming parton's quantum numbers.
*/
virtual bool hasInitialAverage() const {
return matchboxAmplitude() ? matchboxAmplitude()->hasInitialAverage() : false;
}
/**
* Return true, if this amplitude already includes symmetry factors
* for identical outgoing particles.
*/
virtual bool hasFinalStateSymmetry() const {
return matchboxAmplitude() ? matchboxAmplitude()->hasFinalStateSymmetry() : false;
}
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
virtual unsigned int getNLight() const;
/**
* Return the vector that contains the PDG ids of
* the light flavours, which are contained in the
* jet particle group.
*/
virtual vector<int> getNLightJetVec() const;
/**
* Return the vector that contains the PDG ids of
* the heavy flavours, which are contained in the
* jet particle group.
*/
virtual vector<int> getNHeavyJetVec() const;
/**
* Return the vector that contains the PDG ids of
* the light flavours, which are contained in the
* proton particle group.
*/
virtual vector<int> getNLightProtonVec() const;
/**
* Return true, if this matrix element is handled by a BLHA one-loop provider
*/
virtual bool isOLPTree() const {
return matchboxAmplitude() ? matchboxAmplitude()->isOLPTree() : false;
}
/**
* Return true, if this matrix element is handled by a BLHA one-loop provider
*/
virtual bool isOLPLoop() const {
return matchboxAmplitude() ? matchboxAmplitude()->isOLPLoop() : false;
}
/**
* Return true, if colour and spin correlated matrix elements should
* be ordered from the OLP
*/
virtual bool needsOLPCorrelators() const {
return matchboxAmplitude() ? matchboxAmplitude()->needsOLPCorrelators() : true;
}
/**
* Return the process index, if this is an OLP handled matrix element
*/
const vector<int>& olpProcess() const { return theOLPProcess; }
/**
* Set the process index, if this is an OLP handled matrix element
*/
void olpProcess(int pType, int id) {
if ( theOLPProcess.empty() )
theOLPProcess.resize(5,0);
theOLPProcess[pType] = id;
}
/**
* Return true, if this is a real emission matrix element which does
* not require colour correlators.
*/
bool noCorrelations() const {
return theNoCorrelations;
}
/**
* Indicate that this is a real emission matrix element which does
* not require colour correlators.
*/
void needsNoCorrelations() {
theNoCorrelations = true;
}
/**
* Indicate that this is a virtual matrix element which does
* require colour correlators.
*/
void needsCorrelations() {
theNoCorrelations = false;
}
//@}
/** @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;
/**
* Indicate that this matrix element is running alphas by itself.
*/
virtual bool hasRunningAlphaS() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->hasRunningAlphaS();
return false;
}
/**
* Indicate that this matrix element is running alphaew by itself.
*/
virtual bool hasRunningAlphaEW() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->hasRunningAlphaEW();
return false;
}
/**
* Return the scale associated with the phase space point provided
* by the last call to setKinematics().
*/
virtual Energy2 scale() const { return lastScale(); }
/**
* Return the renormalization scale for the last generated phasespace point.
*/
virtual Energy2 factorizationScale() const;
/**
* Get the factorization scale factor
*/
virtual double factorizationScaleFactor() const;
/**
* Return the (QCD) renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScale() const;
/**
* Get the renormalization scale factor
*/
virtual double renormalizationScaleFactor() const;
/**
* Return the QED renormalization scale for the last generated phasespace point.
*/
virtual Energy2 renormalizationScaleQED() const;
/**
* Set veto scales on the particles at the given
* SubProcess which has been generated using this
* matrix element.
*/
virtual void setVetoScales(tSubProPtr) const;
/**
* Return true, if fixed couplings are used.
*/
bool fixedCouplings() const;
/**
* Return true, if fixed couplings are used.
*/
bool fixedQEDCouplings() const;
/**
* Return the value of \f$\alpha_S\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaS(scale()).
*/
virtual double alphaS() const { return lastAlphaS(); }
/**
* Return the value of \f$\alpha_EM\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaEM(scale()).
*/
virtual double alphaEM() const { return lastAlphaEM(); }
/**
* Return true, if this matrix element provides the PDF
* weight for the first incoming parton itself.
*/
virtual bool havePDFWeight1() const;
/**
* Return true, if this matrix element provides the PDF
* weight for the second incoming parton itself.
*/
virtual bool havePDFWeight2() const;
/**
* Set the PDF weight.
*/
void getPDFWeight(Energy2 factorizationScale = ZERO) const;
/**
* Supply the PDF weight for the first incoming parton.
*/
double pdf1(Energy2 factorizationScale = ZERO,
- double xEx = 1.) const;
+ double xEx = 1., double xFactor = 1.) const;
/**
* Supply the PDF weight for the second incoming parton.
*/
double pdf2(Energy2 factorizationScale = ZERO,
- double xEx = 1.) const;
+ double xEx = 1., double xFactor = 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 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 largeNME2(Ptr<ColourBasis>::tptr largeNBasis) 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 the amplitude is DRbar renormalized, otherwise
* MSbar is assumed.
*/
virtual bool isDRbar() 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;
/**
* Simple histogram for accuracy checks
*/
struct AccuracyHistogram {
/**
* The lower bound
*/
double lower;
/**
* The upper bound
*/
double upper;
/**
* The bins, indexed by upper bound.
*/
map<double,double> bins;
/**
* The number of points of same sign
*/
unsigned long sameSign;
/**
* The number of points of opposite sign
*/
unsigned long oppositeSign;
/**
* The number of points being nan or inf
*/
unsigned long nans;
/**
* The overflow
*/
unsigned long overflow;
/**
* The underflow
*/
unsigned long underflow;
/**
* Constructor
*/
AccuracyHistogram(double low = -40.,
double up = 0.,
unsigned int nbins = 80);
/**
* Book two values to be checked for numerical compatibility
*/
void book(double a, double b);
/**
* Write to file.
*/
void dump(const std::string& folder, const std::string& prefix,
const cPDVector& proc) const;
/**
* Write to persistent ostream
*/
void persistentOutput(PersistentOStream&) const;
/**
* Read from persistent istream
*/
void persistentInput(PersistentIStream&);
};
/**
* Perform the check of epsilon pole cancellation.
*/
void logPoles() const;
/**
* Return the virtual corrections
*/
const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const {
return theVirtuals;
}
/**
* Return the virtual corrections
*/
vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() {
return theVirtuals;
}
/**
* Instruct this matrix element to include one-loop corrections
*/
void doOneLoop() { theOneLoop = true; }
/**
* Return true, if this matrix element includes one-loop corrections
*/
bool oneLoop() const { return theOneLoop; }
/**
* Instruct this matrix element to include one-loop corrections but
* no Born contributions
*/
void doOneLoopNoBorn() { theOneLoop = true; theOneLoopNoBorn = true; }
/**
* Return true, if this matrix element includes one-loop corrections
* but no Born contributions
*/
bool oneLoopNoBorn() const { return theOneLoopNoBorn || onlyOneLoop(); }
/**
* Instruct this matrix element to include one-loop corrections but
* no actual loop contributions
*/
void doOneLoopNoLoops() { theOneLoop = true; theOneLoopNoLoops = true; }
/**
* Return true, if this matrix element includes one-loop corrections
* but no actual loop contributions
*/
bool oneLoopNoLoops() const { return theOneLoopNoLoops; }
//@}
/** @name Dipole subtraction */
//@{
/**
* If this matrix element is considered a real
* emission matrix element, return all subtraction
* dipoles needed given a set of subtraction terms
* and underlying Born matrix elements to choose
* from.
*/
vector<Ptr<SubtractionDipole>::ptr>
getDipoles(const vector<Ptr<SubtractionDipole>::ptr>&,
const vector<Ptr<MatchboxMEBase>::ptr>&) const;
/**
* If this matrix element is considered a real emission matrix
* element, but actually neglecting a subclass of the contributing
* diagrams, return true if the given emitter-emission-spectator
* configuration should not be considered when setting up
* subtraction dipoles.
*/
virtual bool noDipole(int,int,int) const { return false; }
/**
* If this matrix element is considered an underlying Born matrix
* element in the context of a subtracted real emission, but
* actually neglecting a subclass of the contributing diagrams,
* return true if the given emitter-spectator configuration
* should not be considered when setting up subtraction dipoles.
*/
virtual bool noDipole(int,int) const { return false; }
/**
* Return the colour correlated matrix element squared with
* respect to the given two partons as appearing in mePartonData(),
* suitably scaled by sHat() to give a dimension-less number.
*/
virtual double colourCorrelatedME2(pair<int,int>) const;
/**
* Return the colour correlated matrix element squared in the
* large-N approximation with respect to the given two partons as
* appearing in mePartonData(), suitably scaled by sHat() to give a
* dimension-less number.
*/
virtual double largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr largeNBasis) const;
/**
* Return the colour and spin correlated matrix element squared for
* the gluon indexed by the first argument using the given
* correlation tensor.
*/
virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& c) const;
/**
* Return the spin correlated matrix element squared for
* the vector boson indexed by the first argument using the given
* correlation tensor.
*/
virtual double spinCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& c) const;
//@}
/** @name Caching and diagnostic information */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches();
/**
* Return true, if verbose
*/
bool verbose() const;
/**
* Return true, if verbose
*/
bool initVerbose() const;
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Print debug information on the last event
*/
virtual void printLastEvent(ostream&) const;
/**
* Write out diagnostic information for
* generateKinematics
*/
void logGenerateKinematics(const double * r) const;
/**
* Write out diagnostic information for
* setting scales
*/
void logSetScale() const;
/**
* Write out diagnostic information for
* pdf evaluation
*/
void logPDFWeight() const;
/**
* Write out diagnostic information for
* me2 evaluation
*/
void logME2() const;
/**
* Write out diagnostic information
* for dsigdr evaluation
*/
void logDSigHatDR() const;
//@}
/** @name Reweight objects */
//@{
/**
* Insert a reweight object
*/
void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }
/**
* Return the reweights
*/
const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }
/**
* Access the reweights
*/
vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }
//@}
/** @name Methods used to setup MatchboxMEBase objects */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return true; }
/**
* Clone this matrix element.
*/
Ptr<MatchboxMEBase>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
void cloneDependencies(const std::string& prefix = "");
/**
* Prepare an xcomb
*/
void prepareXComb(MatchboxXCombData&) const;
/**
* For the given event generation setup return a xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allPBins,
tStdXCombPtr newHead = tStdXCombPtr(),
tMEPtr newME = tMEPtr());
/**
* For the given event generation setup return a dependent xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME = tMEPtr());
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tptr theFactory;
/**
* The phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The amplitude to be used
*/
Ptr<MatchboxAmplitude>::ptr theAmplitude;
/**
* The scale choice object
*/
Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
/**
* The virtual corrections.
*/
vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
/**
* A vector of reweight objects the sum of which
* should be applied to reweight this matrix element
*/
vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
private:
/**
* The subprocess to be considered.
*/
Process theSubprocess;
/**
* 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;
/**
* True, if this matrix element includes one-loop corrections
* but no actual loop contributions (e.g. finite collinear terms)
*/
bool theOneLoopNoLoops;
/**
* The process index, if this is an OLP handled matrix element
*/
vector<int> theOLPProcess;
/**
* Histograms of epsilon^2 pole cancellation
*/
mutable map<cPDVector,AccuracyHistogram> epsilonSquarePoleHistograms;
/**
* Histograms of epsilon pole cancellation
*/
mutable map<cPDVector,AccuracyHistogram> epsilonPoleHistograms;
/**
* True, if this is a real emission matrix element which does
* not require colour correlators.
*/
bool theNoCorrelations;
/**
* Flag which pdfs should be included.
*/
mutable pair<bool,bool> theHavePDFs;
/**
* True, if already checked for which PDFs to include.
*/
mutable bool checkedPDFs;
/**
* Diagnostic Diagram for TreePhaseSpace
*/
void bookMEoverDiaWeight(double x) const;
mutable map<double,double > MEoverDiaWeight;
mutable int Nevents;
/**
* Range of diagram weight verbosity
*/
mutable double theDiagramWeightVerboseDown, theDiagramWeightVerboseUp;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxMEBase & operator=(const MatchboxMEBase &);
};
inline PersistentOStream& operator<<(PersistentOStream& os,
const MatchboxMEBase::AccuracyHistogram& h) {
h.persistentOutput(os);
return os;
}
inline PersistentIStream& operator>>(PersistentIStream& is,
MatchboxMEBase::AccuracyHistogram& h) {
h.persistentInput(is);
return is;
}
}
#endif /* HERWIG_MatchboxMEBase_H */
diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.cc b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
--- a/MatrixElement/Matchbox/Matching/QTildeMatching.cc
+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.cc
@@ -1,445 +1,492 @@
// -*- C++ -*-
//
// QTildeMatching.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 QTildeMatching class.
//
#include "QTildeMatching.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.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"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
using namespace Herwig;
-QTildeMatching::QTildeMatching() {}
+QTildeMatching::QTildeMatching()
+ : theCorrectForXZMismatch(false) {}
QTildeMatching::~QTildeMatching() {}
IBPtr QTildeMatching::clone() const {
return new_ptr(*this);
}
IBPtr QTildeMatching::fullclone() const {
return new_ptr(*this);
}
void QTildeMatching::checkCutoff() {
if ( showerTildeKinematics() ) {
showerTildeKinematics()->
prepare(realCXComb(),bornCXComb());
showerTildeKinematics()->dipole(dipole());
showerTildeKinematics()->getShowerVariables();
}
}
void QTildeMatching::getShowerVariables() {
// already filled from checkCutoff in this case
if ( showerTildeKinematics() )
return;
// get the shower variables
calculateShowerVariables();
// check for the cutoff
dipole()->isAboveCutoff(isAboveCutoff());
// get the hard scale
dipole()->showerHardScale(hardScale());
// check for phase space
dipole()->isInShowerPhasespace(isInShowerPhasespace());
}
bool QTildeMatching::isInShowerPhasespace() const {
assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) &&
"implementation only provided for default and pt cutoff");
Energy qtildeHard = ZERO;
Energy qtilde = dipole()->showerScale();
assert(!dipole()->showerParameters().empty());
double z = dipole()->showerParameters()[0];
// FF
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
qtildeHard =
theQTildeFinder->
calculateFinalFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
bornCXComb()->meMomenta()[dipole()->bornSpectator()],
bornCXComb()->mePartonData()[dipole()->bornEmitter()]->iColour() == PDT::Colour3).first;
}
// FI
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
qtildeHard =
theQTildeFinder->
calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornSpectator()],
bornCXComb()->meMomenta()[dipole()->bornEmitter()],false).second;
}
// IF
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
qtildeHard =
theQTildeFinder->
calculateInitialFinalScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
bornCXComb()->meMomenta()[dipole()->bornSpectator()],false).first;
if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
return false;
}
// II
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
qtildeHard =
theQTildeFinder->
calculateInitialInitialScales(bornCXComb()->meMomenta()[dipole()->bornEmitter()],
bornCXComb()->meMomenta()[dipole()->bornSpectator()]).first;
if ( z < (dipole()->bornEmitter() == 0 ? bornCXComb()->lastX1() : bornCXComb()->lastX2()) )
return false;
}
Energy Qg = theQTildeSudakov->kinScale();
Energy2 pt2 = ZERO;
if ( dipole()->bornEmitter() > 1 ) {
Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu);
else
pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg);
}
if ( dipole()->bornEmitter() < 2 ) {
pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg);
}
if ( pt2 < theQTildeSudakov->pT2min() )
return false;
return qtilde <= qtildeHard && sqrt(pt2) < dipole()->showerHardScale();
}
bool QTildeMatching::isAboveCutoff() const {
assert((theQTildeSudakov->cutOffOption() == 0 || theQTildeSudakov->cutOffOption() == 2) &&
"implementation only provided for default and pt cutoff");
Energy qtilde = dipole()->showerScale();
assert(!dipole()->showerParameters().empty());
double z = dipole()->showerParameters()[0];
Energy Qg = theQTildeSudakov->kinScale();
if ( dipole()->bornEmitter() > 1 ) {
Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
return sqr(z*(1.-z)*qtilde) - sqr(mu) >= theQTildeSudakov->pT2min();
else
return sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg) >= theQTildeSudakov->pT2min();
}
if ( dipole()->bornEmitter() < 2 ) {
return
sqr((1.-z)*qtilde) - z*sqr(Qg) >= theQTildeSudakov->pT2min();
}
return false;
}
CrossSection QTildeMatching::dSigHatDR() const {
assert(!dipole()->showerParameters().empty());
pair<Energy2,double> vars =
make_pair(sqr(dipole()->showerScale()),
dipole()->showerParameters()[0]);
pair<int,int> ij(dipole()->bornEmitter(),
dipole()->bornSpectator());
double ccme2 =
dipole()->underlyingBornME()->largeNColourCorrelatedME2(ij,theLargeNBasis);
ccme2 *=
dipole()->underlyingBornME()->me2() /
dipole()->underlyingBornME()->largeNME2(theLargeNBasis);
Energy2 prop = ZERO;
if ( dipole()->bornEmitter() > 1 ) {
prop =
(realCXComb()->meMomenta()[dipole()->realEmitter()] +
realCXComb()->meMomenta()[dipole()->realEmission()]).m2()
- bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2();
} else {
prop =
2.*vars.second*(realCXComb()->meMomenta()[dipole()->realEmitter()]*
realCXComb()->meMomenta()[dipole()->realEmission()]);
}
// note alphas included downstream from subtractionScaleWeight()
double xme2 = -8.*Constants::pi*ccme2*splitFn(vars)*realXComb()->lastSHat()/prop;
xme2 *=
pow(realCXComb()->lastSHat() / bornCXComb()->lastSHat(),
bornCXComb()->mePartonData().size()-4.);
double bornPDF = bornPDFWeight(dipole()->underlyingBornME()->lastScale());
if ( bornPDF == 0.0 )
return ZERO;
xme2 *= bornPDF;
+ // take care of mismatch between z and x as we are approaching the
+ // hard phase space boundary
+ // TODO get rid of this useless scale option business and simplify PDF handling in here
+ if ( dipole()->bornEmitter() < 2 && theCorrectForXZMismatch ) {
+ Energy2 emissionScale = ZERO;
+ if ( emissionScaleInSubtraction() == showerScale ) {
+ emissionScale = showerFactorizationScale();
+ } else if ( emissionScaleInSubtraction() == realScale ) {
+ emissionScale = dipole()->realEmissionME()->lastScale();
+ } else if ( emissionScaleInSubtraction() == bornScale ) {
+ emissionScale = dipole()->underlyingBornME()->lastScale();
+ }
+ double xzMismatch =
+ dipole()->subtractionParameters()[0] / dipole()->showerParameters()[0];
+ double realCorrectedPDF =
+ dipole()->bornEmitter() == 0 ?
+ dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,
+ xzMismatch) :
+ dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,
+ xzMismatch);
+ double realPDF =
+ dipole()->bornEmitter() == 0 ?
+ dipole()->realEmissionME()->pdf1(emissionScale,theExtrapolationX,1.0) :
+ dipole()->realEmissionME()->pdf2(emissionScale,theExtrapolationX,1.0);
+ if ( realPDF == 0.0 || realCorrectedPDF == 0.0 )
+ return ZERO;
+ xme2 *= realCorrectedPDF / realPDF;
+ }
+
Energy qtilde = sqrt(vars.first);
double z = vars.second;
Energy2 pt2 = ZERO;
Energy Qg = theQTildeSudakov->kinScale();
if ( dipole()->bornEmitter() > 1 ) {
Energy mu = max(Qg,realCXComb()->meMomenta()[dipole()->realEmitter()].mass());
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g )
pt2 = sqr(z*(1.-z)*qtilde) - sqr(mu);
else
pt2 = sqr(z*(1.-z)*qtilde) - sqr((1.-z)*mu) - z*sqr(Qg);
}
if ( dipole()->bornEmitter() < 2 ) {
pt2 = sqr((1.-z)*qtilde) - z*sqr(Qg);
}
assert(pt2 >= ZERO);
xme2 *= hardScaleProfile(dipole()->showerHardScale(),sqrt(pt2));
CrossSection res =
sqr(hbarc) *
realXComb()->jacobian() *
subtractionScaleWeight() *
xme2 /
(2. * realXComb()->lastSHat());
return res;
}
double QTildeMatching::me2() const {
throw Exception() << "QTildeMatching::me2(): Not intented to use. Disable the ShowerApproximationGenerator."
<< Exception::abortnow;
return 0.;
}
void QTildeMatching::calculateShowerVariables() const {
Lorentz5Momentum n;
Energy2 Q2 = ZERO;
const Lorentz5Momentum& pb = bornCXComb()->meMomenta()[dipole()->bornEmitter()];
const Lorentz5Momentum& pc = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
if ( dipole()->bornEmitter() > 1 ) {
Q2 = (pb+pc).m2();
} else {
Q2 = -(pb-pc).m2();
}
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() > 1 ) {
double b = sqr(bornCXComb()->meMomenta()[dipole()->bornEmitter()].m())/Q2;
double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
double lambda = sqrt(1.+sqr(b)+sqr(c)-2.*b-2.*c-2.*b*c);
n = (1.-0.5*(1.-b+c-lambda))*pc - 0.5*(1.-b+c-lambda)*pb;
}
if ( dipole()->bornEmitter() > 1 && dipole()->bornSpectator() < 2 ) {
n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
}
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() > 1 ) {
double c = sqr(bornCXComb()->meMomenta()[dipole()->bornSpectator()].m())/Q2;
n = (1.+c)*pc - c*pb;
}
if ( dipole()->bornEmitter() < 2 && dipole()->bornSpectator() < 2 ) {
n = bornCXComb()->meMomenta()[dipole()->bornSpectator()];
}
// the light-cone condition is numerically not very stable, so we
// explicitly push it on the light-cone here
n.setMass(ZERO);
n.rescaleEnergy();
double z = 0.0;
if ( dipole()->bornEmitter() > 1 ) {
z = 1. -
(n*realCXComb()->meMomenta()[dipole()->realEmission()])/
(n*bornCXComb()->meMomenta()[dipole()->bornEmitter()]);
} else {
z = 1. -
(n*realCXComb()->meMomenta()[dipole()->realEmission()])/
(n*realCXComb()->meMomenta()[dipole()->realEmitter()]);
}
Energy2 qtilde2 = ZERO;
Energy2 q2 = ZERO;
if ( dipole()->bornEmitter() > 1 ) {
q2 =
(realCXComb()->meMomenta()[dipole()->realEmitter()] + realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
qtilde2 = (q2 - bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(z*(1.-z));
} else {
q2 =
-(realCXComb()->meMomenta()[dipole()->realEmitter()] - realCXComb()->meMomenta()[dipole()->realEmission()]).m2();
qtilde2 = (q2 + bornCXComb()->meMomenta()[dipole()->bornEmitter()].m2())/(1.-z);
}
assert(qtilde2 >= ZERO && z >= 0.0 && z <= 1.0);
dipole()->showerScale(sqrt(qtilde2));
dipole()->showerParameters().resize(1);
dipole()->showerParameters()[0] = z;
}
double QTildeMatching::splitFn(const pair<Energy2,double>& vars) const {
const Energy2& qtilde2 = vars.first;
const double& z = vars.second;
double Nc = SM().Nc();
// final state branching
if ( dipole()->bornEmitter() > 1 ) {
// final state quark quark branching
if ( abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 7 ) {
Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
return
((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z);
}
// final state gluon branching
if ( bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == ParticleID::g ) {
if ( realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
// ATTENTION the factor 2 here is intentional as it cancels to the 1/2
// stemming from the large-N colour correlator
return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z));
}
if ( abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
Energy m = realCXComb()->mePartonData()[dipole()->realEmission()]->hardProcessMass();
return (1./2.)*(1.-2.*z*(1.-z)+2.*sqr(m)/(z*(1.-z)*qtilde2));
}
}
// final state squark branching
if ((abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 1000000 &&
abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 1000007) ||
(abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) > 2000000 &&
abs(bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id()) < 2000007)){
Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
return ((sqr(Nc)-1.)/Nc)*(z-sqr(m)/(z*qtilde2))/(1.-z);
}
// final state gluino branching
if (bornCXComb()->mePartonData()[dipole()->bornEmitter()]->id() == 1000021){
Energy m = bornCXComb()->mePartonData()[dipole()->bornEmitter()]->hardProcessMass();
return Nc*(1.+sqr(z)-2.*sqr(m)/(z*qtilde2))/(1.-z);
}
}
// initial state branching
if ( dipole()->bornEmitter() < 2 ) {
// g/g
if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
// see above for factor of 2
return 2.*Nc*(z/(1.-z)+(1.-z)/z+z*(1.-z));
}
// q/q
if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
realCXComb()->mePartonData()[dipole()->realEmission()]->id() == ParticleID::g ) {
return
((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(z))/(1.-z);
}
// g/q
if ( realCXComb()->mePartonData()[dipole()->realEmitter()]->id() == ParticleID::g &&
abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
return (1./2.)*(1.-2.*z*(1.-z));
}
// q/g
if ( abs(realCXComb()->mePartonData()[dipole()->realEmitter()]->id()) < 7 &&
abs(realCXComb()->mePartonData()[dipole()->realEmission()]->id()) < 7 ) {
return
((sqr(Nc)-1.)/(2.*Nc))*(1+sqr(1.-z))/z;
}
}
return 0.0;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void QTildeMatching::doinit() {
assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
theShowerHandler->init();
theQTildeFinder->init();
theQTildeSudakov->init();
if ( theShowerHandler->scaleFactorOption() < 2 ) {
hardScaleFactor(theShowerHandler->hardScaleFactor());
factorizationScaleFactor(theShowerHandler->factorizationScaleFactor());
renormalizationScaleFactor(theShowerHandler->renormalizationScaleFactor());
}
ShowerApproximation::doinit();
}
void QTildeMatching::doinitrun() {
assert(theShowerHandler && theQTildeFinder && theQTildeSudakov);
theShowerHandler->initrun();
theQTildeFinder->initrun();
theQTildeSudakov->initrun();
ShowerApproximation::doinitrun();
}
void QTildeMatching::persistentOutput(PersistentOStream & os) const {
- os << theQTildeFinder << theQTildeSudakov << theShowerHandler;
+ os << theQTildeFinder << theQTildeSudakov
+ << theShowerHandler << theCorrectForXZMismatch;
}
void QTildeMatching::persistentInput(PersistentIStream & is, int) {
- is >> theQTildeFinder >> theQTildeSudakov >> theShowerHandler;
+ is >> theQTildeFinder >> theQTildeSudakov
+ >> theShowerHandler >> theCorrectForXZMismatch;
}
// *** 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<QTildeMatching,Herwig::ShowerApproximation>
describeHerwigQTildeMatching("Herwig::QTildeMatching", "HwShower.so HwQTildeMatching.so");
void QTildeMatching::Init() {
static ClassDocumentation<QTildeMatching> documentation
("QTildeMatching implements NLO matching with the default shower.");
static Reference<QTildeMatching,QTildeFinder> interfaceQTildeFinder
("QTildeFinder",
"Set the partner finder to calculate hard scales.",
&QTildeMatching::theQTildeFinder, false, false, true, false, false);
static Reference<QTildeMatching,QTildeSudakov> interfaceQTildeSudakov
("QTildeSudakov",
"Set the partner finder to calculate hard scales.",
&QTildeMatching::theQTildeSudakov, false, false, true, false, false);
static Reference<QTildeMatching,ShowerHandler> interfaceShowerHandler
("ShowerHandler",
"",
&QTildeMatching::theShowerHandler, false, false, true, true, false);
+ static Switch<QTildeMatching,bool> interfaceCorrectForXZMismatch
+ ("CorrectForXZMismatch",
+ "Correct for x/z mismatch near hard phase space boundary.",
+ &QTildeMatching::theCorrectForXZMismatch, false, false, false);
+ static SwitchOption interfaceCorrectForXZMismatchYes
+ (interfaceCorrectForXZMismatch,
+ "Yes",
+ "Include the correction factor.",
+ true);
+ static SwitchOption interfaceCorrectForXZMismatchNo
+ (interfaceCorrectForXZMismatch,
+ "No",
+ "Do not include the correction factor.",
+ false);
+
}
diff --git a/MatrixElement/Matchbox/Matching/QTildeMatching.h b/MatrixElement/Matchbox/Matching/QTildeMatching.h
--- a/MatrixElement/Matchbox/Matching/QTildeMatching.h
+++ b/MatrixElement/Matchbox/Matching/QTildeMatching.h
@@ -1,195 +1,201 @@
// -*- C++ -*-
//
// QTildeMatching.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_QTildeMatching_H
#define Herwig_QTildeMatching_H
//
// This is the declaration of the QTildeMatching class.
//
#include "Herwig++/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
#include "Herwig++/Shower/ShowerHandler.h"
#include "Herwig++/Shower/Default/QTildeFinder.h"
#include "Herwig++/Shower/Default/QTildeSudakov.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief QTildeMatching implements NLO matching with the default shower.
*
*/
class QTildeMatching: public Herwig::ShowerApproximation {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
QTildeMatching();
/**
* The destructor.
*/
virtual ~QTildeMatching();
//@}
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;
/**
* Determine if the configuration is below or above the cutoff.
*/
virtual void checkCutoff();
/**
* Determine all kinematic variables which are not provided by the
* dipole kinematics; store all shower variables in the respective
* dipole object for later use.
*/
virtual void getShowerVariables();
protected:
/**
* 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;
/**
* 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;
/**
* Calculate qtilde^2 and z for the splitting considered
*/
void calculateShowerVariables() const;
/**
* Return the splitting function as a function of the kinematic
* variables
*/
double splitFn(const pair<Energy2,double>&) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
QTildeMatching & operator=(const QTildeMatching &);
/**
* The shower handler to be used
*/
Ptr<ShowerHandler>::ptr theShowerHandler;
/**
* The qtilde partner finder for calculating the hard scales
*/
Ptr<QTildeFinder>::ptr theQTildeFinder;
/**
* The qtilde Sudakov to access the cutoff
*/
Ptr<QTildeSudakov>::ptr theQTildeSudakov;
+ /**
+ * True, if PDF weight should be corrected for z/x mismatch at the
+ * hard phase space boundary
+ */
+ bool theCorrectForXZMismatch;
+
};
}
#endif /* Herwig_QTildeMatching_H */
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,713 +1,712 @@
// -*- 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/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/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
using namespace Herwig;
ShowerApproximation::ShowerApproximation()
: HandlerBase(),
- theBelowCutoff(false),
+ theExtrapolationX(1.0), theBelowCutoff(false),
theFFPtCut(1.0*GeV), theFFScreeningScale(ZERO),
theFIPtCut(1.0*GeV), theFIScreeningScale(ZERO),
theIIPtCut(1.0*GeV), theIIScreeningScale(ZERO),
theRestrictPhasespace(true), theHardScaleFactor(1.0),
theRenormalizationScaleFactor(1.0), theFactorizationScaleFactor(1.0),
- theExtrapolationX(1.0),
theRealEmissionScaleInSubtraction(showerScale),
theBornScaleInSubtraction(showerScale),
theEmissionScaleInSubtraction(showerScale),
theRealEmissionScaleInSplitting(showerScale),
theBornScaleInSplitting(showerScale),
theEmissionScaleInSplitting(showerScale),
theRenormalizationScaleFreeze(1.*GeV),
theFactorizationScaleFreeze(1.*GeV),
theProfileScales(true),
theProfileRho(0.3), maxPtIsMuF(false) {}
ShowerApproximation::~ShowerApproximation() {}
void ShowerApproximation::setLargeNBasis() {
assert(dipole()->realEmissionME()->matchboxAmplitude());
if ( !dipole()->realEmissionME()->matchboxAmplitude()->treeAmplitudes() )
return;
if ( !theLargeNBasis ) {
if ( !dipole()->realEmissionME()->matchboxAmplitude()->colourBasis() )
throw Exception() << "ShowerApproximation::setLargeNBasis(): Expecting a colour basis object."
<< Exception::abortnow;
theLargeNBasis =
dipole()->realEmissionME()->matchboxAmplitude()->colourBasis()->cloneMe();
theLargeNBasis->clear();
theLargeNBasis->doLargeN();
}
}
void ShowerApproximation::setDipole(Ptr<SubtractionDipole>::tptr dip) {
theDipole = dip;
setLargeNBasis();
}
Ptr<SubtractionDipole>::tptr ShowerApproximation::dipole() const { return theDipole; }
Ptr<TildeKinematics>::tptr
ShowerApproximation::showerTildeKinematics() const {
return Ptr<TildeKinematics>::tptr();
}
Ptr<InvertedTildeKinematics>::tptr
ShowerApproximation::showerInvertedTildeKinematics() const {
return Ptr<InvertedTildeKinematics>::tptr();
}
void ShowerApproximation::checkCutoff() {
assert(!showerTildeKinematics());
}
void ShowerApproximation::getShowerVariables() {
// check for the cutoff
dipole()->isAboveCutoff(isAboveCutoff());
// get the hard scale
dipole()->showerHardScale(hardScale());
// set the shower scale and variables for completeness
dipole()->showerScale(dipole()->lastPt());
dipole()->showerParameters().resize(1);
dipole()->showerParameters()[0] = dipole()->lastZ();
// check for phase space
dipole()->isInShowerPhasespace(isInShowerPhasespace());
}
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;
}
Energy ShowerApproximation::hardScale() const {
if ( !maxPtIsMuF ) {
if ( !bornCXComb()->mePartonData()[0]->coloured() &&
!bornCXComb()->mePartonData()[1]->coloured() ) {
Energy maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m();
maxPt *= hardScaleFactor();
return maxPt;
}
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->mt());
if ( maxPt == generator()->maximumCMEnergy() )
maxPt = (bornCXComb()->meMomenta()[0] + bornCXComb()->meMomenta()[1]).m();
maxPt *= hardScaleFactor();
return maxPt;
} else {
return hardScaleFactor()*sqrt(bornCXComb()->lastCentralScale());
}
}
double ShowerApproximation::hardScaleProfile(Energy hard, Energy soft) const {
double x = soft/hard;
if ( theProfileScales ) {
if ( x > 1. ) {
return 0.;
} else if ( x <= 1. && x > 1. - theProfileRho ) {
return sqr(1.-x)/(2.*sqr(theProfileRho));
} else if ( x <= 1. - theProfileRho &&
x > 1. - 2.*theProfileRho ) {
return 1. - sqr(1.-2.*theProfileRho-x)/(2.*sqr(theProfileRho));
} else {
return 1.;
}
}
if ( x <= 1. )
return 1.;
return 0.;
}
bool ShowerApproximation::isInShowerPhasespace() const {
if ( !dipole()->isAboveCutoff() )
return false;
if ( !restrictPhasespace() )
return true;
InvertedTildeKinematics& kinematics =
const_cast<InvertedTildeKinematics&>(*dipole()->invertedTildeKinematics());
tcStdXCombPtr tmpreal = kinematics.realXComb();
tcStdXCombPtr tmpborn = kinematics.bornXComb();
Ptr<SubtractionDipole>::tptr tmpdip = kinematics.dipole();
Energy hard = dipole()->showerHardScale();
Energy pt = dipole()->lastPt();
double z = dipole()->lastZ();
pair<double,double> zbounds(0.,1.);
kinematics.dipole(const_ptr_cast<Ptr<SubtractionDipole>::tptr>(theDipole));
kinematics.prepare(realCXComb(),bornCXComb());
if ( pt > hard ) {
kinematics.dipole(tmpdip);
kinematics.prepare(tmpreal,tmpborn);
return false;
}
try {
zbounds = kinematics.zBounds(pt,hard);
} catch(...) {
kinematics.dipole(tmpdip);
kinematics.prepare(tmpreal,tmpborn);
throw;
}
kinematics.dipole(tmpdip);
kinematics.prepare(tmpreal,tmpborn);
return z > zbounds.first && z < zbounds.second;
}
Energy2 ShowerApproximation::showerEmissionScale() const {
Energy2 mur = sqr(dipole()->lastPt());
if ( dipole()->bornEmitter() > 1 &&
dipole()->bornSpectator() > 1 ) {
return mur + sqr(ffScreeningScale());
} else if ( ( dipole()->bornEmitter() > 1 &&
dipole()->bornSpectator() < 2 ) ||
( dipole()->bornEmitter() < 2 &&
dipole()->bornSpectator() > 1 ) ) {
return mur + sqr(fiScreeningScale());
} else {
assert(dipole()->bornEmitter() < 2 &&
dipole()->bornSpectator() < 2);
return mur + sqr(iiScreeningScale());
}
return mur;
}
Energy2 ShowerApproximation::bornRenormalizationScale() const {
return
sqr(dipole()->underlyingBornME()->renormalizationScaleFactor()) *
dipole()->underlyingBornME()->renormalizationScale();
}
Energy2 ShowerApproximation::bornFactorizationScale() const {
return
sqr(dipole()->underlyingBornME()->factorizationScaleFactor()) *
dipole()->underlyingBornME()->factorizationScale();
}
Energy2 ShowerApproximation::realRenormalizationScale() const {
return
sqr(dipole()->realEmissionME()->renormalizationScaleFactor()) *
dipole()->realEmissionME()->renormalizationScale();
}
Energy2 ShowerApproximation::realFactorizationScale() const {
return
sqr(dipole()->realEmissionME()->factorizationScaleFactor()) *
dipole()->realEmissionME()->factorizationScale();
}
double ShowerApproximation::bornPDFWeight(Energy2 muf) const {
if ( !bornCXComb()->mePartonData()[0]->coloured() &&
!bornCXComb()->mePartonData()[1]->coloured() )
return 1.;
if ( muf < sqr(theFactorizationScaleFreeze) )
muf = sqr(theFactorizationScaleFreeze);
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(Energy2 muf) const {
if ( !realCXComb()->mePartonData()[0]->coloured() &&
!realCXComb()->mePartonData()[1]->coloured() )
return 1.;
if ( muf < sqr(theFactorizationScaleFreeze) )
muf = sqr(theFactorizationScaleFreeze);
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;
}
double ShowerApproximation::scaleWeight(int rScale, int bScale, int eScale) const {
double emissionAlpha = 1.;
Energy2 emissionScale = ZERO;
Energy2 showerscale = ZERO;
if ( eScale == showerScale || bScale == showerScale || eScale == showerScale ) {
showerscale = showerRenormalizationScale();
if ( showerscale < sqr(theRenormalizationScaleFreeze) )
showerscale = sqr(theFactorizationScaleFreeze);
}
if ( eScale == showerScale ) {
emissionAlpha = SM().alphaS(showerscale);
emissionScale = showerFactorizationScale();
} else if ( eScale == realScale ) {
emissionAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS();
emissionScale = dipole()->realEmissionME()->lastScale();
} else if ( eScale == bornScale ) {
emissionAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS();
emissionScale = dipole()->underlyingBornME()->lastScale();
}
double emissionPDF = realPDFWeight(emissionScale);
double couplingFactor = 1.;
if ( bScale != rScale ) {
double bornAlpha = 1.;
if ( bScale == showerScale ) {
bornAlpha = SM().alphaS(showerscale);
} else if ( bScale == realScale ) {
bornAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS();
} else if ( bScale == bornScale ) {
bornAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS();
}
double realAlpha = 1.;
if ( rScale == showerScale ) {
realAlpha = SM().alphaS(showerscale);
} else if ( rScale == realScale ) {
realAlpha = dipole()->realEmissionME()->lastXComb().lastAlphaS();
} else if ( rScale == bornScale ) {
realAlpha = dipole()->underlyingBornME()->lastXComb().lastAlphaS();
}
couplingFactor *=
pow(realAlpha/bornAlpha,(double)(dipole()->underlyingBornME()->orderInAlphaS()));
}
Energy2 hardScale = ZERO;
if ( bScale == showerScale ) {
hardScale = showerFactorizationScale();
} else if ( bScale == realScale ) {
hardScale = dipole()->realEmissionME()->lastScale();
} else if ( bScale == bornScale ) {
hardScale = dipole()->underlyingBornME()->lastScale();
}
double bornPDF = bornPDFWeight(hardScale);
if ( bornPDF < 1e-8 )
bornPDF = 0.0;
if ( emissionPDF < 1e-8 )
emissionPDF = 0.0;
if ( emissionPDF == 0.0 || bornPDF == 0.0 )
return 0.0;
return
emissionAlpha * emissionPDF *
couplingFactor / bornPDF;
}
double ShowerApproximation::channelWeight(int emitter, int emission,
int spectator, int) const {
double cfac = 1.;
double Nc = generator()->standardModel()->Nc();
if (realCXComb()->mePartonData()[emitter]->iColour() == PDT::Colour8){
if (realCXComb()->mePartonData()[emission]->iColour() == PDT::Colour8)
cfac = Nc;
else if ( realCXComb()->mePartonData()[emission]->iColour() == PDT::Colour3 ||
realCXComb()->mePartonData()[emission]->iColour() == PDT::Colour3bar)
cfac = 0.5;
else assert(false);
}
else if ((realCXComb()->mePartonData()[emitter] ->iColour() == PDT::Colour3 ||
realCXComb()->mePartonData()[emitter] ->iColour() == PDT::Colour3bar))
cfac = (sqr(Nc)-1.)/(2.*Nc);
else assert(false);
// 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 cfac *GeV2 * pipk / ( pipj * ( pipj + pjpk ) );
}
return
cfac * GeV2 / (realCXComb()->meMomenta()[emitter] * realCXComb()->meMomenta()[emission]);
}
double ShowerApproximation::channelWeight() const {
double currentChannel = channelWeight(dipole()->realEmitter(),
dipole()->realEmission(),
dipole()->realSpectator(),
dipole()->bornEmitter());
if ( currentChannel == 0. )
return 0.;
double sum = 0.;
for ( vector<Ptr<SubtractionDipole>::tptr>::const_iterator dip =
dipole()->partnerDipoles().begin();
dip != dipole()->partnerDipoles().end(); ++dip )
sum += channelWeight((**dip).realEmitter(),
(**dip).realEmission(),
(**dip).realSpectator(),
(**dip).bornEmitter());
assert(sum > 0.0);
return currentChannel / sum;
}
// 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 << theLargeNBasis
<< theBornXComb << theRealXComb << theTildeXCombs << theDipole << theBelowCutoff
<< ounit(theFFPtCut,GeV) << ounit(theFFScreeningScale,GeV)
<< ounit(theFIPtCut,GeV) << ounit(theFIScreeningScale,GeV)
<< ounit(theIIPtCut,GeV) << ounit(theIIScreeningScale,GeV)
<< theRestrictPhasespace << theHardScaleFactor
<< theRenormalizationScaleFactor << theFactorizationScaleFactor
<< theExtrapolationX
<< theRealEmissionScaleInSubtraction << theBornScaleInSubtraction
<< theEmissionScaleInSubtraction << theRealEmissionScaleInSplitting
<< theBornScaleInSplitting << theEmissionScaleInSplitting
<< ounit(theRenormalizationScaleFreeze,GeV)
<< ounit(theFactorizationScaleFreeze,GeV)
<< theProfileScales << theProfileRho << maxPtIsMuF;
}
void ShowerApproximation::persistentInput(PersistentIStream & is, int) {
is >> theLargeNBasis
>> theBornXComb >> theRealXComb >> theTildeXCombs >> theDipole >> theBelowCutoff
>> iunit(theFFPtCut,GeV) >> iunit(theFFScreeningScale,GeV)
>> iunit(theFIPtCut,GeV) >> iunit(theFIScreeningScale,GeV)
>> iunit(theIIPtCut,GeV) >> iunit(theIIScreeningScale,GeV)
>> theRestrictPhasespace >> theHardScaleFactor
>> theRenormalizationScaleFactor >> theFactorizationScaleFactor
>> theExtrapolationX
>> theRealEmissionScaleInSubtraction >> theBornScaleInSubtraction
>> theEmissionScaleInSubtraction >> theRealEmissionScaleInSplitting
>> theBornScaleInSplitting >> theEmissionScaleInSplitting
>> iunit(theRenormalizationScaleFreeze,GeV)
>> iunit(theFactorizationScaleFreeze,GeV)
>> theProfileScales >> theProfileRho >> maxPtIsMuF;
}
// *** 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", "Herwig.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 Parameter<ShowerApproximation,Energy> interfaceFFScreeningScale
("FFScreeningScale",
"Set the screening scale",
&ShowerApproximation::theFFScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<ShowerApproximation,Energy> interfaceFIScreeningScale
("FIScreeningScale",
"Set the screening scale",
&ShowerApproximation::theFIScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<ShowerApproximation,Energy> interfaceIIScreeningScale
("IIScreeningScale",
"Set the screening scale",
&ShowerApproximation::theIIScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
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> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The hard scale factor.",
&ShowerApproximation::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerApproximation,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The hard scale factor.",
&ShowerApproximation::theFactorizationScaleFactor, 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, 1.0, 0.0, 1.0,
false, false, Interface::limited);
static Switch<ShowerApproximation,int> interfaceRealEmissionScaleInSubtraction
("RealEmissionScaleInSubtraction",
"Set the scale choice for the real emission cross section in the matching subtraction.",
&ShowerApproximation::theRealEmissionScaleInSubtraction, showerScale, false, false);
static SwitchOption interfaceRealEmissionScaleInSubtractionRealScale
(interfaceRealEmissionScaleInSubtraction,
"RealScale",
"Use the real emission scale.",
realScale);
static SwitchOption interfaceRealEmissionScaleInSubtractionBornScale
(interfaceRealEmissionScaleInSubtraction,
"BornScale",
"Use the Born scale.",
bornScale);
static SwitchOption interfaceRealEmissionScaleInSubtractionShowerScale
(interfaceRealEmissionScaleInSubtraction,
"ShowerScale",
"Use the shower scale",
showerScale);
static Switch<ShowerApproximation,int> interfaceBornScaleInSubtraction
("BornScaleInSubtraction",
"Set the scale choice for the Born cross section in the matching subtraction.",
&ShowerApproximation::theBornScaleInSubtraction, showerScale, false, false);
static SwitchOption interfaceBornScaleInSubtractionRealScale
(interfaceBornScaleInSubtraction,
"RealScale",
"Use the real emission scale.",
realScale);
static SwitchOption interfaceBornScaleInSubtractionBornScale
(interfaceBornScaleInSubtraction,
"BornScale",
"Use the Born scale.",
bornScale);
static SwitchOption interfaceBornScaleInSubtractionShowerScale
(interfaceBornScaleInSubtraction,
"ShowerScale",
"Use the shower scale",
showerScale);
static Switch<ShowerApproximation,int> interfaceEmissionScaleInSubtraction
("EmissionScaleInSubtraction",
"Set the scale choice for the emission in the matching subtraction.",
&ShowerApproximation::theEmissionScaleInSubtraction, showerScale, false, false);
static SwitchOption interfaceEmissionScaleInSubtractionRealScale
(interfaceEmissionScaleInSubtraction,
"RealScale",
"Use the real emission scale.",
realScale);
static SwitchOption interfaceEmissionScaleInSubtractionEmissionScale
(interfaceEmissionScaleInSubtraction,
"BornScale",
"Use the Born scale.",
bornScale);
static SwitchOption interfaceEmissionScaleInSubtractionShowerScale
(interfaceEmissionScaleInSubtraction,
"ShowerScale",
"Use the shower scale",
showerScale);
static Switch<ShowerApproximation,int> interfaceRealEmissionScaleInSplitting
("RealEmissionScaleInSplitting",
"Set the scale choice for the real emission cross section in the splitting.",
&ShowerApproximation::theRealEmissionScaleInSplitting, showerScale, false, false);
static SwitchOption interfaceRealEmissionScaleInSplittingRealScale
(interfaceRealEmissionScaleInSplitting,
"RealScale",
"Use the real emission scale.",
realScale);
static SwitchOption interfaceRealEmissionScaleInSplittingBornScale
(interfaceRealEmissionScaleInSplitting,
"BornScale",
"Use the Born scale.",
bornScale);
static SwitchOption interfaceRealEmissionScaleInSplittingShowerScale
(interfaceRealEmissionScaleInSplitting,
"ShowerScale",
"Use the shower scale",
showerScale);
static Switch<ShowerApproximation,int> interfaceBornScaleInSplitting
("BornScaleInSplitting",
"Set the scale choice for the Born cross section in the splitting.",
&ShowerApproximation::theBornScaleInSplitting, showerScale, false, false);
static SwitchOption interfaceBornScaleInSplittingRealScale
(interfaceBornScaleInSplitting,
"RealScale",
"Use the real emission scale.",
realScale);
static SwitchOption interfaceBornScaleInSplittingBornScale
(interfaceBornScaleInSplitting,
"BornScale",
"Use the Born scale.",
bornScale);
static SwitchOption interfaceBornScaleInSplittingShowerScale
(interfaceBornScaleInSplitting,
"ShowerScale",
"Use the shower scale",
showerScale);
static Switch<ShowerApproximation,int> interfaceEmissionScaleInSplitting
("EmissionScaleInSplitting",
"Set the scale choice for the emission in the splitting.",
&ShowerApproximation::theEmissionScaleInSplitting, showerScale, false, false);
static SwitchOption interfaceEmissionScaleInSplittingRealScale
(interfaceEmissionScaleInSplitting,
"RealScale",
"Use the real emission scale.",
realScale);
static SwitchOption interfaceEmissionScaleInSplittingEmissionScale
(interfaceEmissionScaleInSplitting,
"BornScale",
"Use the Born scale.",
bornScale);
static SwitchOption interfaceEmissionScaleInSplittingShowerScale
(interfaceEmissionScaleInSplitting,
"ShowerScale",
"Use the shower scale",
showerScale);
static Parameter<ShowerApproximation,Energy> interfaceRenormalizationScaleFreeze
("RenormalizationScaleFreeze",
"The freezing scale for the renormalization scale.",
&ShowerApproximation::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<ShowerApproximation,Energy> interfaceFactorizationScaleFreeze
("FactorizationScaleFreeze",
"The freezing scale for the factorization scale.",
&ShowerApproximation::theFactorizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<ShowerApproximation,bool> interfaceProfileScales
("ProfileScales",
"Switch on or off use of profile scales.",
&ShowerApproximation::theProfileScales, true, false, false);
static SwitchOption interfaceProfileScalesYes
(interfaceProfileScales,
"Yes",
"Use profile scales.",
true);
static SwitchOption interfaceProfileScalesNo
(interfaceProfileScales,
"No",
"Use a hard cutoff.",
false);
static Parameter<ShowerApproximation,double> interfaceProfileRho
("ProfileRho",
"The rho parameter of the profile scales.",
&ShowerApproximation::theProfileRho, 0.3, 0.0, 1.0,
false, false, Interface::limited);
static Reference<ShowerApproximation,ColourBasis> interfaceLargeNBasis
("LargeNBasis",
"Set the large-N colour basis implementation.",
&ShowerApproximation::theLargeNBasis, false, false, true, true, false);
static Switch<ShowerApproximation,bool> interfaceMaxPtIsMuF
("MaxPtIsMuF",
"",
&ShowerApproximation::maxPtIsMuF, false, false, false);
static SwitchOption interfaceMaxPtIsMuFYes
(interfaceMaxPtIsMuF,
"Yes",
"",
true);
static SwitchOption interfaceMaxPtIsMuFNo
(interfaceMaxPtIsMuF,
"No",
"",
false);
}
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,640 +1,640 @@
// -*- 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"
#include "Herwig++/MatrixElement/Matchbox/Utility/ColourBasis.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.fh"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.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
* H events
*/
virtual bool hasHEvents() const { return true; }
/**
* Return true, if this shower approximation will require tilde
* XCombs for the real phase space point generated
*/
virtual bool needsTildeXCombs() const { return false; }
/**
* Return true, if this shower approximation will require
* a truncated parton shower
*/
virtual bool needsTruncatedShower() const { return false; }
/**
* Return the tilde kinematics object returning the shower
* kinematics parametrization if different from the nominal dipole
* mappings.
*/
virtual Ptr<TildeKinematics>::tptr showerTildeKinematics() const;
/**
* Return the tilde kinematics object returning the shower
* kinematics parametrization if different from the nominal dipole
* mappings.
*/
virtual Ptr<InvertedTildeKinematics>::tptr showerInvertedTildeKinematics() const;
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>::tptr);
/**
* Return the dipole in charge for the emission
*/
Ptr<SubtractionDipole>::tptr dipole() const;
/**
* Return true, if this matching is capable of spin correlations.
*/
virtual bool hasSpinCorrelations() const { return false; }
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; }
/**
* Return the screening scale to be applied for final-final dipoles.
*/
Energy ffScreeningScale() const { return theFFScreeningScale; }
/**
* Return the screening scale to be applied for final-initial dipoles.
*/
Energy fiScreeningScale() const { return theFIScreeningScale; }
/**
* Return the screening scale to be applied for initial-initial dipoles.
*/
Energy iiScreeningScale() const { return theIIScreeningScale; }
/**
* Return the shower renormalization scale
*/
virtual Energy2 showerEmissionScale() const;
/**
* Return the shower renormalization scale
*/
Energy2 showerRenormalizationScale() const {
return sqr(renormalizationScaleFactor())*showerEmissionScale();
}
/**
* Return the shower factorization scale
*/
Energy2 showerFactorizationScale() const {
return sqr(factorizationScaleFactor())*showerEmissionScale();
}
/**
* Return the Born renormalization scale
*/
Energy2 bornRenormalizationScale() const;
/**
* Return the Born factorization scale
*/
Energy2 bornFactorizationScale() const;
/**
* Return the real emission renormalization scale
*/
Energy2 realRenormalizationScale() const;
/**
* Return the real emission factorization scale
*/
Energy2 realFactorizationScale() const;
/**
* Enumerate possible scale choices
*/
enum ScaleChoices {
bornScale = 0,
/** Use the born scales */
realScale = 1,
/** Use the real scales */
showerScale = 2
/** Use the shower scales */
};
/**
* Return the scale choice in the real emission cross section to be
* used in the matching subtraction.
*/
int realEmissionScaleInSubtraction() const { return theRealEmissionScaleInSubtraction; }
/**
* Return the scale choice in the born cross section to be
* used in the matching subtraction.
*/
int bornScaleInSubtraction() const { return theBornScaleInSubtraction; }
/**
* Return the scale choice in the emission contribution to be
* used in the matching subtraction.
*/
int emissionScaleInSubtraction() const { return theEmissionScaleInSubtraction; }
/**
* Return the scale choice in the real emission cross section to be
* used in the splitting.
*/
int realEmissionScaleInSplitting() const { return theRealEmissionScaleInSplitting; }
/**
* Return the scale choice in the born cross section to be
* used in the splitting.
*/
int bornScaleInSplitting() const { return theBornScaleInSplitting; }
/**
* Return the scale choice in the emission contribution to be
* used in the splitting.
*/
int emissionScaleInSplitting() const { return theEmissionScaleInSplitting; }
/**
* Return the scale weight
*/
double scaleWeight(int rScale, int bScale, int eScale) const;
/**
* Return the scale weight for the matching subtraction
*/
double subtractionScaleWeight() const {
return scaleWeight(realEmissionScaleInSubtraction(),
bornScaleInSubtraction(),
emissionScaleInSubtraction());
}
/**
* Return the scale weight for the splitting
*/
double splittingScaleWeight() const {
return scaleWeight(realEmissionScaleInSplitting(),
bornScaleInSplitting(),
emissionScaleInSplitting());
}
public:
/**
* Return true, if the phase space restrictions of the dipole shower should
* be applied.
*/
bool restrictPhasespace() const { return theRestrictPhasespace; }
/**
* Return true if we are to use profile scales
*/
bool profileScales() const { return theProfileScales; }
/**
* Return the scale factor for the hard scale
*/
double hardScaleFactor() const { return theHardScaleFactor; }
/**
* Set the scale factor for the hard scale
*/
void hardScaleFactor(double f) { theHardScaleFactor = f; }
/**
* Return a scale profile towards the hard scale
*/
virtual double hardScaleProfile(Energy hard, Energy soft) const;
/**
* Get the factorization scale factor
*/
double factorizationScaleFactor() const { return theFactorizationScaleFactor; }
/**
* Get the renormalization scale factor
*/
double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }
/**
* Set the factorization scale factor
*/
void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }
/**
* Set the renormalization scale factor
*/
void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }
/**
* Determine if the configuration is below or above the cutoff.
*/
virtual void checkCutoff();
/**
* Determine all kinematic variables which are not provided by the
* dipole kinematics; store all shower variables in the respective
* dipole object for later use.
*/
virtual void getShowerVariables();
/**
* 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 the Born PDF weight
*/
double bornPDFWeight(Energy2 muF) const;
/**
* Return the real emission PDF weight
*/
double realPDFWeight(Energy2 muF) const;
protected:
/**
* 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;
/**
* 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;
/**
* Return the relevant hard scale
*/
virtual Energy hardScale() const;
public:
/**
* Generate a weight for the given dipole channel
*/
virtual double channelWeight(int emitter, int emission,
int spectator, int bemitter) const;
/**
* Generate a normalized weight taking into account all channels
*/
virtual double channelWeight() 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).
public:
/**
* A large-N colour basis to be used when reproducing the shower
* kernels.
*/
Ptr<ColourBasis>::tptr largeNBasis() const { return theLargeNBasis; }
protected:
/**
* A large-N colour basis to be used when reproducing the shower
* kernels.
*/
Ptr<ColourBasis>::ptr theLargeNBasis;
/**
* Set the large-N basis
*/
void setLargeNBasis();
+ /**
+ * The x value from which on we extrapolate PDFs for numerically stable ratios.
+ */
+ double theExtrapolationX;
+
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>::tptr 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;
/**
* An optional screening scale for final-final dipoles; see
* DipoleSplittingKernel
*/
Energy theFFScreeningScale;
/**
* The pt cut to be applied for final-initial dipoles.
*/
Energy theFIPtCut;
/**
* An optional screening scale for final-initial dipoles; see
* DipoleSplittingKernel
*/
Energy theFIScreeningScale;
/**
* The pt cut to be applied for initial-initial dipoles.
*/
Energy theIIPtCut;
/**
* An optional screening scale for initial-initial dipoles; see
* DipoleSplittingKernel
*/
Energy theIIScreeningScale;
/**
* 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 scale factor for the renormalization scale
*/
double theRenormalizationScaleFactor;
/**
* The scale factor for the factorization scale
*/
double theFactorizationScaleFactor;
/**
- * The x value from which on we extrapolate PDFs for numerically stable ratios.
- */
- double theExtrapolationX;
-
- /**
* The scale choice in the real emission cross section to be
* used in the matching subtraction.
*/
int theRealEmissionScaleInSubtraction;
/**
* The scale choice in the born cross section to be
* used in the matching subtraction.
*/
int theBornScaleInSubtraction;
/**
* The scale choice in the emission contribution to be
* used in the matching subtraction.
*/
int theEmissionScaleInSubtraction;
/**
* The scale choice in the real emission cross section to be
* used in the splitting.
*/
int theRealEmissionScaleInSplitting;
/**
* The scale choice in the born cross section to be
* used in the splitting.
*/
int theBornScaleInSplitting;
/**
* The scale choice in the emission contribution to be
* used in the splitting.
*/
int theEmissionScaleInSplitting;
/**
* A freezing value for the renormalization scale
*/
Energy theRenormalizationScaleFreeze;
/**
* A freezing value for the factorization scale
*/
Energy theFactorizationScaleFreeze;
/**
* True if we are to use profile scales
*/
bool theProfileScales;
/**
* The profile scale parameter
*/
double theProfileRho;
/**
* True if maximum pt should be deduced from the factorization scale
*/
bool maxPtIsMuF;
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 */

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:27 PM (4 h, 33 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486597
Default Alt Text
(148 KB)

Event Timeline