Page MenuHomeHEPForge
Files F7878348

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,1653 +1,1652 @@
// -*- 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/MatrixElement/Matchbox/Base/MergerBase.h"
#include "Herwig/API/RunDirectories.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
#include "Herwig/MatrixElement/HardVertex.h"
#include <cctype>
#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) {}
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<long> MatchboxMEBase::getNLightJetVec() const { return factory()->nLightJetVec(); }
vector<long> MatchboxMEBase::getNHeavyJetVec() const { return factory()->nHeavyJetVec(); }
vector<long> 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 (auto const & d : diags )
add(d);
return;
}
throw Exception()
<< "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
<< "Please check your setup." << Exception::runerror;
}
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::runerror;
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 ( auto const & fit : cflows)
flows.push_back(new ColourLines(ColourBasis::cfstring(fit)));
}
Selector<const ColourLines *> res;
for ( auto const & f : flows ) 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 ( auto const & p : sub->outgoing() ) {
out.push_back(p->data().iSpin());
hard.push_back(p);
}
// 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 if (hard[k]->data().iSpin() == PDT::Spin0 ) {
ScalarWaveFunction(hard[k],k > 1 ? outgoing : incoming,
k > 1 ? true : false);
}
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 ( auto const & p : sub->outgoing() )
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);
if (theMerger){
theMerger->setME(this);
theMerger->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();
if (theMerger&&!theMerger->generateKinematics(r)){
return false;
}
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::runerror;
return false;
}
int MatchboxMEBase::nDim() const {
if ( lastMatchboxXComb() )
return nDimPhasespace() + nDimAmplitude() + nDimInsertions();
int ampAdd = 0;
if ( matchboxAmplitude() ) {
ampAdd = matchboxAmplitude()->nDimAdditional();
}
int insertionAdd = 0;
for ( auto const & v : virtuals() ) {
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::runerror;
return 0;
}
void MatchboxMEBase::setScale(Energy2 ren, Energy2 fac) const {
if ( haveX1X2() ) {
lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
}
Energy2 fcscale = (fac == ZERO) ? factorizationScale() : fac;
Energy2 fscale = fcscale*sqr(factorizationScaleFactor());
Energy2 rscale = (ren == ZERO ? renormalizationScale() : ren)*sqr(renormalizationScaleFactor());
Energy2 ewrscale = renormalizationScaleQED();
lastXCombPtr()->lastScale(fscale);
lastXCombPtr()->lastCentralScale(fcscale);
lastXCombPtr()->lastShowerScale(showerScale());
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::runerror;
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::runerror;
return ZERO;
}
Energy2 MatchboxMEBase::renormalizationScaleQED() const {
if ( scaleChoice() ) {
return scaleChoice()->renormalizationScaleQED();
}
return renormalizationScale();
}
Energy2 MatchboxMEBase::showerScale() const {
if ( scaleChoice() ) {
return scaleChoice()->showerScale();
}
throw Exception()
<< "MatchboxMEBase::showerScale() expects a MatchboxScaleChoice object.\n"
<< "Please check your setup." << Exception::runerror;
return ZERO;
}
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, double xFactor) const {
assert(lastXCombPtr()->partonBins().first->pdf());
if ( xEx < 1. && lastX1()*xFactor >= xEx ) {
return
( ( 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()*xFactor)/lastX1()/xFactor;
}
double MatchboxMEBase::pdf2(Energy2 fscale, double xEx, double xFactor) const {
assert(lastXCombPtr()->partonBins().second->pdf());
if ( xEx < 1. && lastX2()*xFactor >= xEx ) {
return
( ( 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()*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::runerror;
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::runerror;
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 ( auto const & c : counts) {
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::prefactor()const{
return (sqr(hbarc)/(2.*lastSHat())) *jacobian()* lastMEPDFWeight();
}
CrossSection MatchboxMEBase::dSigHatDRB() const {
getPDFWeight();
lastME2(me2());
return oneLoopNoBorn()?ZERO:prefactor() * lastME2();
}
CrossSection MatchboxMEBase::dSigHatDRV() const {
getPDFWeight();
lastME2(me2());
return ( oneLoop() && !oneLoopNoLoops() )?(prefactor() * oneLoopInterference()):ZERO;
}
CrossSection MatchboxMEBase::dSigHatDRI() const {
getPDFWeight();
lastME2(me2());
CrossSection res=ZERO;
if (oneLoop() &&!onlyOneLoop()) {
for ( auto const & v : virtuals()) {
v->setXComb(lastXCombPtr());
res += v->dSigHatDR();
}
if ( checkPoles() && oneLoop() )
logPoles();
}
return res;
}
CrossSection MatchboxMEBase::dSigHatDRAlphaDiff(double alpha) const {
getPDFWeight();
lastME2(me2());
CrossSection res=ZERO;
for ( auto const & v: virtuals() ) {
v->setXComb(lastXCombPtr());
res+=v->dSigHatDRAlphaDiff( alpha);
}
return res;
}
CrossSection MatchboxMEBase::dSigHatDR() const {
getPDFWeight();
if (theMerger){
lastMECrossSection(theMerger->MergingDSigDR());
return lastMECrossSection();
}else if (lastXCombPtr()->willPassCuts() ) {
lastMECrossSection(dSigHatDRB()+
dSigHatDRV()+
dSigHatDRI());
return lastMECrossSection();
}
lastME2(ZERO);
lastMECrossSection(ZERO);
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::runerror;
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 ( ! (isfinite(a) && isfinite(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 ( auto const & v : virtuals()) {
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<SubtractionDipolePtr>
MatchboxMEBase::getDipoles(const vector<SubtractionDipolePtr>& dipoles,
const vector<MatchboxMEBasePtr> & borns,bool slim) const {
vector<SubtractionDipolePtr> 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<SubtractionDipolePtr> matchDipoles;
for ( auto const & d : dipoles ) {
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<SubtractionDipolePtr> matchDipoles2;
for ( auto const & d : matchDipoles ) {
if ( !d->canHandleSplitting(rep,emitter,emission) )
continue;
matchDipoles2.push_back(d);
}
if ( matchDipoles2.empty() )
continue;
map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo> mergeInfo;
for ( auto const & d : diagrams() ) {
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<SubtractionDipolePtr> matchDipoles3;
for ( auto const & d : matchDipoles2 ) {
if ( ! d->canHandleSpectator(rep,spectator) )
continue;
matchDipoles3.push_back(d);
}
if ( matchDipoles3.empty() )
continue;
if ( noDipole(emitter,emission,spectator) )
continue;
for ( auto const & d : matchDipoles3 ) {
if ( !d->canHandle(rep,emitter,emission,spectator) )
continue;
for ( auto const & b : borns ) {
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,slim);
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;
if ( theMerger) {
dname << fullName();
if (theOneLoopNoBorn) dname << ".virtual" << "." ;
dname << b->name() << "."
<< d->name() << ".[("
<< emitter << "," << emission << ")," << spectator << "]";
} else {
dname << fullName() << "." << b->name() << "."
<< d->name() << ".[("
<< emitter << "," << emission << ")," << spectator << "]";
}
if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) )
throw Exception() << "MatchboxMEBase::getDipoles(): Dipole " << dname.str() << " already existing." << Exception::runerror;
if ( !factory()->reweighters().empty() ) {
for ( auto const & rw : factory()->reweighters())
nDipole->addReweighter(rw);
}
if ( !factory()->preweighters().empty() ) {
for ( auto const & rw : factory()->preweighters() )
nDipole->addPreweighter(rw);
}
nDipole->cloneDependencies(dname.str(),slim);
}
}
}
}
}
}
vector<Ptr<SubtractionDipole>::tptr> partners;
copy(res.begin(),res.end(),back_inserter(partners));
for ( auto const & d : res )
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::runerror;
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::runerror;
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::runerror;
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::runerror;
return 0.;
}
-void MatchboxMEBase::flushCaches() {
+void MatchboxMEBase::flushCaches() {
+ if ( theMerger )theMerger->flushCaches();
MEBase::flushCaches();
if ( matchboxAmplitude() )
matchboxAmplitude()->flushCaches();
for ( auto const & r : reweights() )
r->flushCaches();
for ( auto const & v : virtuals())
v->flushCaches();
- // The Merger cache is flushed in
- // generateKinematics
}
void MatchboxMEBase::setKinematics() {
MEBase::setKinematics();
if ( theMerger )
theMerger->setKinematics();
}
void MatchboxMEBase::clearKinematics() {
MEBase::clearKinematics();
if ( theMerger )
theMerger->clearKinematics();
}
const MergerBasePtr MatchboxMEBase::merger() const {
return theMerger;
}
MergerBasePtr MatchboxMEBase::merger() {
return theMerger;
}
void MatchboxMEBase::merger(MergerBasePtr v) {
theMerger = v;
}
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) << "\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,bool slim) {
if ( phasespace() && !slim ) {
Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
throw Exception() << "MatchboxMEBase::cloneDependencies(): Phasespace generator " << pname.str() << " already existing." << Exception::runerror;
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 Exception() << "MatchboxMEBase::cloneDependencies(): Amplitude " << pname.str() << " already existing." << Exception::runerror;
}
myAmplitude->cloneDependencies(pname.str(),slim);
matchboxAmplitude(myAmplitude);
amplitude(myAmplitude);
matchboxAmplitude()->orderInGs(orderInAlphaS());
matchboxAmplitude()->orderInGem(orderInAlphaEW());
}
if ( scaleChoice() &&!slim ) {
Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
throw Exception() << "MatchboxMEBase::cloneDependencies(): Scale choice " << pname.str() << " already existing." << Exception::runerror;
scaleChoice(myScaleChoice);
}
for ( auto & rw : theReweights ) {
Ptr<MatchboxReweightBase>::ptr myReweight = rw->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << rw->name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw Exception() << "MatchboxMEBase::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror;
myReweight->cloneDependencies(pname.str());
rw = myReweight;
}
for ( auto & v : virtuals()) {
Ptr<MatchboxInsertionOperator>::ptr myIOP = v->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << v->name();
if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
throw Exception() << "MatchboxMEBase::cloneDependencies(): Insertion operator " << pname.str() << " already existing." << Exception::runerror;
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 ( auto const & v : virtuals() )
insertionAdd = max(insertionAdd,v->nDimAdditional());
xc.nDimInsertions(insertionAdd);
xc.nLight(getNLight());
if(xc.nLightJetVec().empty())
for (auto const & id : getNLightJetVec())
xc.nLightJetVec( id );
if(xc.nHeavyJetVec().empty())
for (auto const & id :getNHeavyJetVec())
xc.nHeavyJetVec(id);
if(xc.nLightProtonVec().empty())
for (auto const & id : getNLightProtonVec())
xc.nLightProtonVec(id);
xc.olpId(olpProcess());
if ( initVerbose() ) {
ostringstream fname_strm;
// only allow alphanumeric, / and _ in filename
for (const char c : name()) {
switch (c) {
case '+' : fname_strm << "+"; break;
case '-' : fname_strm << "-"; break;
case '~' : fname_strm << "_tilde"; break;
case ']' : break;
case ',' : fname_strm << "__"; break;
default : fname_strm << (isalnum(c) ? c : '_'); break;
}
}
fname_strm << ".diagrams";
const string fname = fname_strm.str();
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
<< theMerger
<< theOLPProcess << theNoCorrelations
<< theHavePDFs << checkedPDFs;
}
void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theFactory >> thePhasespace
>> theAmplitude >> theScaleChoice >> theVirtuals
>> theReweights >> theSubprocess >> theOneLoop
>> theOneLoopNoBorn >> theOneLoopNoLoops
>> epsilonSquarePoleHistograms >> epsilonPoleHistograms
>> theMerger
>> theOLPProcess >> theNoCorrelations
>> theHavePDFs >> checkedPDFs;
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 (auto const & rw : theReweights)
rw->init();
for (auto const & v : virtuals() )
v->init();
}
void MatchboxMEBase::doinitrun() {
MEBase::doinitrun();
if ( matchboxAmplitude() )
matchboxAmplitude()->initrun();
if ( phasespace() )
phasespace()->initrun();
if ( scaleChoice() )
scaleChoice()->initrun();
for (auto const & rw : theReweights)
rw->initrun();
for (auto const & v : virtuals() )
v->initrun();
}
void MatchboxMEBase::dofinish() {
MEBase::dofinish();
for (auto const & b : epsilonSquarePoleHistograms ) {
b.second.dump(factory()->poleData(),"epsilonSquarePoles-",b.first);
}
for (auto const & b : epsilonPoleHistograms ) {
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/Utility/MatchboxXCombData.cc b/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc
--- a/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc
+++ b/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc
@@ -1,335 +1,347 @@
// -*- C++ -*-
//
// MatchboxXCombData.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 MatchboxXCombData class.
//
#include "MatchboxXCombData.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
MatchboxXCombData::MatchboxXCombData()
: theCrossingSign(1.0), theCalculateTreeAmplitudes(true),
theCalculateOneLoopAmplitudes(true), theCalculateTreeME2(true),
theLastTreeME2(0.0), theCalculateLargeNME2(true),
theLastLargeNME2(0.0), theCalculateOneLoopInterference(true),
theLastOneLoopInterference(0.0), theCalculateOneLoopPoles(true),
theLastOneLoopPoles(0.0,0.0),
theColourBasisDim(0), theNDimPhasespace(0),
theNDimAmplitude(0), theNDimInsertions(0),
theSymmetryFactor(0.0), theOLPMomenta(0),
filledOLPMomenta(false), theExternalId(0),
theInitialized(false), filledExternalMomenta(false) {
flushCaches();
}
unsigned int MatchboxXCombData::theNLight(0);
vector<long> MatchboxXCombData::theNLightJetVec=vector<long> ();
vector<long> MatchboxXCombData::theNHeavyJetVec=vector<long>() ;
vector<long> MatchboxXCombData::theNLightProtonVec=vector<long>() ;
MatchboxXCombData::~MatchboxXCombData() {
if ( theOLPMomenta ) {
delete[] theOLPMomenta;
theOLPMomenta = 0;
}
for ( vector<double*>::iterator k = theExternalMomenta.begin();
k != theExternalMomenta.end(); ++k ) {
delete[] *k;
*k = 0;
}
theExternalMomenta.clear();
}
MatchboxXCombData::MatchboxXCombData(tMEPtr newME)
: theCrossingSign(1.0), theCalculateTreeAmplitudes(true),
theCalculateOneLoopAmplitudes(true), theCalculateTreeME2(true),
theLastTreeME2(0.0), theCalculateLargeNME2(true),
theLastLargeNME2(0.0), theCalculateOneLoopInterference(true),
theLastOneLoopInterference(0.0), theCalculateOneLoopPoles(true),
theLastOneLoopPoles(0.0,0.0),
theColourBasisDim(0), theNDimPhasespace(0),
theNDimAmplitude(0), theNDimInsertions(0),
theSymmetryFactor(0.0), theOLPMomenta(0),
filledOLPMomenta(false), theExternalId(0),
theInitialized(false), filledExternalMomenta(false) {
flushCaches();
theMatchboxME = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(newME);
theSubtractionDipole = dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(newME);
Ptr<SubtractedME>::tptr subme =
dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(newME);
if ( !theMatchboxME && !theSubtractionDipole && subme )
theMatchboxME = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(subme->head());
assert(theMatchboxME || theSubtractionDipole);
if ( theMatchboxME )
theFactory = theMatchboxME->factory();
else if ( theSubtractionDipole )
theFactory = theSubtractionDipole->realEmissionME()->factory();
assert(theFactory);
}
double MatchboxXCombData::ReshuffleEquation::operator() (double xi) const {
double res = -q/GeV;
cPDVector::const_iterator dit = dBegin;
vector<Lorentz5Momentum>::const_iterator mit = mBegin;
for ( ; dit != dEnd; ++dit, ++mit ) {
map<long,Energy>::const_iterator rm =
reshuffleMap->find((**dit).id());
Energy2 tmass2 =
rm == reshuffleMap->end() ?
mit->mass2() : sqr(rm->second);
res += sqrt(sqr(xi)*(mit->vect().mag2()) + tmass2)/GeV;
}
return res;
}
void MatchboxXCombData::reshuffle(vector<Lorentz5Momentum>& momenta,
const cPDVector& mePartonData,
const map<long,Energy>& reshuffleMap) const {
if ( momenta.size() == 3 ) // nothing to do; don't throw an exception
return;
bool needDoSomething = false;
for ( cPDVector::const_iterator d = mePartonData.begin() + 2;
d != mePartonData.end(); ++d )
needDoSomething |= reshuffleMap.find((**d).id()) != reshuffleMap.end();
if ( !needDoSomething )
return;
Lorentz5Momentum Q(ZERO,ZERO,ZERO,ZERO);
for ( vector<Lorentz5Momentum>::const_iterator p = momenta.begin()+2;
p != momenta.end(); ++p )
Q += *p;
Boost beta = Q.findBoostToCM();
bool needBoost = (beta.mag2() > Constants::epsilon);
if ( needBoost ) {
for ( vector<Lorentz5Momentum>::iterator p = momenta.begin()+2;
p != momenta.end(); ++p )
p->boost(beta);
}
ReshuffleEquation solve(Q.m(),
mePartonData.begin() + 2,
mePartonData.end(),
momenta.begin() + 2, &reshuffleMap);
double xi = 1.;
GSLBisection solver(1e-10,1e-8,10000);
try {
xi = solver.value(solve,0.0,1.1);
} catch (GSLBisection::GSLerror) {
throw Veto();
} catch (GSLBisection::IntervalError) {
throw Veto();
}
cPDVector::const_iterator d = mePartonData.begin() + 2;
vector<Lorentz5Momentum>::iterator p = momenta.begin() + 2;
for ( ; d != mePartonData.end(); ++d, ++p ) {
p->setVect(xi*p->vect());
map<long,Energy>::const_iterator mit =
reshuffleMap.find((**d).id());
Energy2 newMass2 = mit == reshuffleMap.end() ? p->mass2() : sqr(mit->second);
p->setE(sqrt(p->vect().mag2() + newMass2));
p->rescaleMass();
}
if ( needBoost ) {
for ( vector<Lorentz5Momentum>::iterator p = momenta.begin()+2;
p != momenta.end(); ++p )
p->boost(-beta);
}
}
void MatchboxXCombData::fillOLPMomenta(const vector<Lorentz5Momentum>& memomenta,
const cPDVector& mePartonData,
const map<long,Energy>& reshuffleMap) {
if ( filledOLPMomenta )
return;
if ( !reshuffleMap.empty() && memomenta.size() > 3 ) {
vector<Lorentz5Momentum> reshuffled = memomenta;
reshuffle(reshuffled,mePartonData,reshuffleMap);
fillOLPMomenta(reshuffled);
return;
}
if ( !theOLPMomenta ) {
theOLPMomenta = new double[5*memomenta.size()];
}
for ( size_t p = 0; p < memomenta.size(); ++p ) {
theOLPMomenta[5*p] = memomenta[p].t()/GeV;
theOLPMomenta[5*p+1] = memomenta[p].x()/GeV;
theOLPMomenta[5*p+2] = memomenta[p].y()/GeV;
theOLPMomenta[5*p+3] = memomenta[p].z()/GeV;
theOLPMomenta[5*p+4] = memomenta[p].mass()/GeV;
}
filledOLPMomenta = true;
}
void MatchboxXCombData::fillExternalMomenta(const vector<Lorentz5Momentum>& memomenta) {
if ( filledExternalMomenta )
return;
if ( theExternalMomenta.empty() ) {
theExternalMomenta.resize(memomenta.size());
for ( size_t k = 0; k < memomenta.size(); ++k )
theExternalMomenta[k] = new double[4];
}
for ( size_t p = 0; p < memomenta.size(); ++p ) {
theExternalMomenta[p][0] = memomenta[p].t()/GeV;
theExternalMomenta[p][1] = memomenta[p].x()/GeV;
theExternalMomenta[p][2] = memomenta[p].y()/GeV;
theExternalMomenta[p][3] = memomenta[p].z()/GeV;
}
filledExternalMomenta = true;
}
Ptr<MatchboxFactory>::tcptr MatchboxXCombData::factory() const {
return theFactory;
}
Ptr<MatchboxMEBase>::tptr MatchboxXCombData::matchboxME() const {
return theMatchboxME;
}
Ptr<SubtractionDipole>::tptr MatchboxXCombData::subtractionDipole() const {
return theSubtractionDipole;
}
void MatchboxXCombData::flushCaches() {
theCalculateTreeAmplitudes = true;
theCalculateOneLoopAmplitudes = true;
theCalculateTreeME2 = true;
theCalculateLargeNME2 = true;
theCalculateOneLoopInterference = true;
theCalculateOneLoopPoles = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateColourCorrelators.begin();
f != theCalculateColourCorrelators.end(); ++f )
f->second = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateLargeNColourCorrelators.begin();
f != theCalculateLargeNColourCorrelators.end(); ++f )
f->second = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateColourSpinCorrelators.begin();
f != theCalculateColourSpinCorrelators.end(); ++f )
f->second = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateSpinCorrelators.begin();
f != theCalculateSpinCorrelators.end(); ++f )
f->second = true;
filledOLPMomenta = false;
filledExternalMomenta = false;
+ theLastAmplitudes.clear();
+ theLastLargeNAmplitudes.clear();
+ theLastOneLoopAmplitudes.clear();
+ theColourCorrelators.clear();
+ theLargeNColourCorrelators.clear();
+ theColourSpinCorrelators.clear();
+ theSpinCorrelators.clear();
+ theAmplitudeRandomNumbers.clear();
+ theInsertionRandomNumbers.clear();
+ theDiagramWeights.clear();
+ theHelJamp.clear();
+ theLNHelJamp.clear();
}
void MatchboxXCombData::putCVector(PersistentOStream& os, const CVector& v) {
size_t n = v.size();
os << n;
for ( size_t k = 0; k < n; k++ )
os << v(k);
}
void MatchboxXCombData::getCVector(PersistentIStream& is, CVector& v) {
size_t n; is >> n;
v.resize(n);
Complex value;
for ( size_t k = 0; k < n; k++ ) {
is >> value; v(k) = value;
}
}
void MatchboxXCombData::putAmplitudeMap(PersistentOStream& os, const map<vector<int>,CVector>& amps) {
os << amps.size();
for ( map<vector<int>,CVector>::const_iterator a = amps.begin();
a != amps.end(); ++a ) {
os << a->first;
putCVector(os,a->second);
}
}
void MatchboxXCombData::getAmplitudeMap(PersistentIStream& is, map<vector<int>,CVector>& amps) {
size_t n; is >> n;
for ( size_t k = 0; k < n; ++k ) {
vector<int> hel; is >> hel;
getCVector(is,amps[hel]);
}
}
void MatchboxXCombData::persistentOutput(PersistentOStream & os) const {
os << theFactory << theMatchboxME << theSubtractionDipole << theCrossingMap
<< theAmplitudeToColourMap << theColourToAmplitudeMap
<< theCrossingSign << theAmplitudePartonData << ounit(theAmplitudeMomenta,GeV)
<< ounit(theLastRenormalizationScale,GeV2)
<< theCalculateTreeAmplitudes /* << theLastAmplitudes << theLastLargeNAmplitudes */
<< theCalculateOneLoopAmplitudes /* << theLastOneLoopAmplitudes */
<< theCalculateTreeME2 << theLastTreeME2
<< theCalculateLargeNME2 << theLastLargeNME2
<< theCalculateOneLoopInterference
<< theLastOneLoopInterference << theCalculateOneLoopPoles
<< theLastOneLoopPoles << theCalculateColourCorrelators
<< theColourCorrelators << theCalculateLargeNColourCorrelators
<< theLargeNColourCorrelators << theCalculateColourSpinCorrelators
<< theColourSpinCorrelators << theCalculateSpinCorrelators
<< theSpinCorrelators << theNLight
<< theNLightJetVec << theNHeavyJetVec << theNLightProtonVec
<< theColourBasisDim
<< theNDimPhasespace << theNDimAmplitude << theNDimInsertions
<< theAmplitudeRandomNumbers << theInsertionRandomNumbers
<< theDiagramWeights << theSingularLimits// << theLastSingularLimit
<< theStandardModel << theSymmetryFactor
<< theOLPId << theExternalId;
putAmplitudeMap(os,theLastAmplitudes);
putAmplitudeMap(os,theLastLargeNAmplitudes);
putAmplitudeMap(os,theLastOneLoopAmplitudes);
}
void MatchboxXCombData::persistentInput(PersistentIStream & is, int) {
is >> theFactory >> theMatchboxME >> theSubtractionDipole >> theCrossingMap
>> theAmplitudeToColourMap >> theColourToAmplitudeMap
>> theCrossingSign >> theAmplitudePartonData >> iunit(theAmplitudeMomenta,GeV)
>> iunit(theLastRenormalizationScale,GeV2)
>> theCalculateTreeAmplitudes /* >> theLastAmplitudes >> theLastLargeNAmplitudes */
>> theCalculateOneLoopAmplitudes /* >> theLastOneLoopAmplitudes */
>> theCalculateTreeME2 >> theLastTreeME2
>> theCalculateLargeNME2 >> theLastLargeNME2
>> theCalculateOneLoopInterference
>> theLastOneLoopInterference >> theCalculateOneLoopPoles
>> theLastOneLoopPoles >> theCalculateColourCorrelators
>> theColourCorrelators >> theCalculateLargeNColourCorrelators
>> theLargeNColourCorrelators >> theCalculateColourSpinCorrelators
>> theColourSpinCorrelators >> theCalculateSpinCorrelators
>> theSpinCorrelators >> theNLight
>> theNLightJetVec >> theNHeavyJetVec >> theNLightProtonVec
>> theColourBasisDim
>> theNDimPhasespace >> theNDimAmplitude >> theNDimInsertions
>> theAmplitudeRandomNumbers >> theInsertionRandomNumbers
>> theDiagramWeights >> theSingularLimits// >> theLastSingularLimit
>> theStandardModel >> theSymmetryFactor
>> theOLPId >> theExternalId;
getAmplitudeMap(is,theLastAmplitudes);
getAmplitudeMap(is,theLastLargeNAmplitudes);
getAmplitudeMap(is,theLastOneLoopAmplitudes);
}
diff --git a/Shower/Dipole/Merging/Merger.cc b/Shower/Dipole/Merging/Merger.cc
--- a/Shower/Dipole/Merging/Merger.cc
+++ b/Shower/Dipole/Merging/Merger.cc
@@ -1,1492 +1,1492 @@
// -*- C++ -*-
//
// Merger.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright ( C ) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL , see COPYING for details.
// Please respect the MCnet academic guidelines , see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined , non-templated member
// functions of the Merger class.
//
#include "Merger.h"
#include "Node.h"
#include "MergingFactory.h"
// other includes when needed below.
using namespace Herwig;
IBPtr Merger::clone() const {
return new_ptr( *this );
}
IBPtr Merger::fullclone() const {
return new_ptr( *this );
}
namespace {
double decideClustering(const NodePtr sub,const NodePtr head,bool& pro){
if( sub != head ){// at least one history step -> unitarisation
if ( UseRandom::rndbool() ){ pro = true; return -2.; }
else{ pro = false; return 2.; }
} // no ordered history -> no projection
else{ pro = false; return 1.; }
}
}
CrossSection Merger::MergingDSigDRBornStandard( ){
// get the history for the process
const NodePtr productionNode =
currentNode()-> getHistory( true, DSH()->hardScaleFactor() );
// decide if to cluster
weight = decideClustering(productionNode, currentNode(), projected);
// check if we only want to calculate the current multiplicity.
if(notOnlyMulti()) return ZERO;
// Check for cuts on the production proces.
if ( !productionNode->xcomb()->willPassCuts() ) return ZERO;
// calculate the staring scale for the production node
Energy startscale = CKKW_StartScale( productionNode );
// fill history with caluclation of sudakov supression
fillHistory( startscale , productionNode , currentNode() );
// fill the projector -> return the scale of the last splitting
currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
// the weight has three components to get shower history weight
weight *= history.back().weight* // Sudakov suppression
alphaReweight()* // alpha_s reweight
pdfReweight(); // pdf reweight
// If weight is zero return.
if( weight == ZERO ) return ZERO;
//calculate the cross section
return weight*TreedSigDR( startscale , 1. );
}
CrossSection Merger::MergingDSigDRVirtualStandard( ){
// get the history for the process
const NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() );
// decide if to cluster
weight = decideClustering(productionNode,currentNode(),projected);
// Check for cuts on the production proces.
if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
// calculate the staring scale
Energy startscale = CKKW_StartScale( productionNode );
// fill history with caluclation of sudakov supression
fillHistory( startscale , productionNode , currentNode() );
// fill the projector -> return the scale of the last splitting
currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
// the weight has three components to get shower history weight
double ww1 = history.back().weight;
double ww2 = alphaReweight();
double ww3 = pdfReweight();
weight *= ww1*ww2*ww3;
// If weight is zero return.
if( weight == 0. )return ZERO;
// calculate the cross section for virtual contribution
// and various insertion operators.
CrossSection matrixElement = LoopdSigDR( startscale );
// Now calculate the expansion of the shower history.
// first: the born contibution:
CrossSection bornWeight = currentME()->dSigHatDRB();
// second: expansion of pdf ,alpha_s-ratio and sudakov suppression.
double w1 = -sumPdfReweightExpansion();
double w2 = -sumAlphaSReweightExpansion();
double w3 = -sumFillHistoryExpansion();
// put together the expansion weights.
CrossSection expansionweight = ( w1+w2+w3 )*(theShowerExpansionWeights?1.:0.)
*bornWeight
*SM().alphaS()/( 2.*ThePEG::Constants::pi );
// [ DEBUG ]
if ( currentNode()->legsize() == 5 && Debug::level > 2 )
debugVirt(weight,w1,w2,w3,matrixElement,ww1,ww2,ww3,productionNode,bornWeight);
// return with correction that ME was calculated with fixed alpha_s
return weight* as( startscale*DSH()->renFac() )/
SM().alphaS()* ( matrixElement+expansionweight );
}
CrossSection Merger::MergingDSigDRRealStandard(){
if ( currentNode()->children().empty() ) {
throw Exception()
<< "Real emission contribution without underlying born."
<< "These are finite contibutions already handled in LO merging."
<< Exception::abortnow;
}
// check for IR Safe Cutoff
if( !currentNode()->allAbove( theIRSafePT ) )return ZERO;
auto inOutPair = currentNode()->getInOut();
NodePtr randomChild = currentNode()->randomChild();
bool meRegion =matrixElementRegion(
inOutPair.first ,
inOutPair.second ,
randomChild->pT() ,
theMergePt );
if ( meRegion )return MergingDSigDRRealAllAbove( );
if ( UseRandom::rndbool() )
return 2.*MergingDSigDRRealBelowSubReal( );
return 2.*MergingDSigDRRealBelowSubInt( );
}
CrossSection Merger::MergingDSigDRRealAllAbove( ){
//If all dipoles pts are above , we cluster to this dipole.
NodePtr CLNode = currentNode()->randomChild();
// Check if phase space poing is in ME region--> else rm PSP
if ( !CLNode->children().empty() ) {
auto inOutPair = CLNode->getInOut();
NodePtr randomChild = CLNode->randomChild();
if( !matrixElementRegion( inOutPair.first ,
inOutPair.second ,
randomChild->pT() ,
theMergePt ) )return ZERO;
}
// first find the history for the acctual Node
NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() );
// check if the randomly choosen CLNode is part of the history.
// If CLNode is not part of the history , dont calculate the Real contribution
// else multiply the real contribution with N ( number of children ).
// this returns the sudakov suppression according to
// the clustering of the born parts.
bool inhist = CLNode->isInHistoryOf( productionNode );
// get the history for the clustered Node.
productionNode = CLNode-> getHistory( false , DSH()->hardScaleFactor() );
// decide if to cluster
weight = decideClustering(productionNode,CLNode,projected);
// Check for cuts on the production process.
if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
// calculate the staring scale
Energy startscale = CKKW_StartScale( productionNode );
// fill history with caluclation of sudakov supression
fillHistory( startscale , productionNode , CLNode );
// fill the projector -> return the scale of the last splitting
currentNode()->runningPt( fillProjector( projected ? 2 : 1 ) );
// the weight has three components to get shower history weight
weight *= history.back().weight*alphaReweight()*pdfReweight();
if( weight == 0. )return ZERO;
// The inhist flag produces the correct cluster sensity.
CrossSection me = ( inhist?TreedSigDR( startscale ):ZERO );
// calculate the dipole
CrossSection dip = CLNode->calcDip( startscale* currentME()->renFac() );
CrossSection res = weight*as( startscale*DSH()->renFac() )/SM().alphaS()*
currentNode()->children().size()*( me -dip );
// [ DEBUG ]
if ( currentNode()->legsize() == 6&&Debug::level > 2 )
debugReal("RAA",weight,me,dip);
return res;
}
CrossSection Merger::MergingDSigDRRealBelowSubReal( ){
NodePtr HistNode = currentNode()->randomChild();
if ( !HistNode->children().empty() ) {
auto inOutPair = HistNode->getInOut();
NodePtr randomChild = HistNode->randomChild(); // Here we make sure that clustering and splitting are invertible
if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO;
}
const NodePtr productionNode = HistNode-> getHistory( false , DSH()->hardScaleFactor() );
weight = decideClustering(productionNode,HistNode,projected);
if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
Energy startscale = CKKW_StartScale( productionNode );
fillHistory( startscale , productionNode , HistNode );
currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
weight *= history.back().weight*alphaReweight()*pdfReweight();
if( weight == 0. )return ZERO;
CrossSection sumPS = ZERO;
for( auto const & child : currentNode()->children() ){
if ( child->allAbove( mergePt() ) ){
if( ( child )->pT()>mergePt()/theRealSubtractionRatio )
sumPS += child->calcPs( startscale* currentME()->renFac() );
else
sumPS += child->calcDip( startscale* currentME()->renFac() );
}else{
assert( child->pT()>mergePt() );
}
}
CrossSection me = TreedSigDR( startscale );
// [ DEBUG ]
if ( currentNode()->legsize() == 6&&Debug::level > 2 )
debugReal("RBSR",weight,me,sumPS);
//Here we subtract the PS ( and below the dynamical cutoff the Dip )
return weight*as( startscale*DSH()->renFac() )/SM().alphaS()*
( me-sumPS );
}
CrossSection Merger::MergingDSigDRRealBelowSubInt( ){
if( currentNode()->children().empty() )return ZERO;
NodePtr CLNode = currentNode()->randomChild();
if( CLNode->pT()<mergePt()/theRealSubtractionRatio )return ZERO;
if ( !CLNode->children().empty() ) {
auto inOutPair = CLNode->getInOut( );
NodePtr randomChild = CLNode->randomChild(); // Here we make sure that clustering and splitting are invertible
if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO;
}
const NodePtr productionNode = CLNode-> getHistory( false , DSH()->hardScaleFactor() );
weight = decideClustering(productionNode,CLNode,projected);
if ( !CLNode->allAbove( mergePt() ) )return ZERO;
if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
Energy startscale = CKKW_StartScale( productionNode );
fillHistory( startscale , productionNode , CLNode );
currentNode()->runningPt( fillProjector( projected ? 2 : 1 ) );
weight *= history.back().weight*alphaReweight()*pdfReweight();
if( weight == 0. )return ZERO;
pair<CrossSection , CrossSection> DipAndPs =
CLNode->calcDipandPS( startscale* currentME()->renFac() );
// [ DEBUG ]
if ( currentNode()->legsize() == 6&&Debug::level > 2 )
debugReal("RBSI",weight,DipAndPs.second,DipAndPs.first);
//Here we add the PS and subtrac the Dip ( only above the dynamical cutoff )
return weight*as( startscale*DSH()->renFac() )/SM().alphaS()*
currentNode()->children().size()*( DipAndPs.second-DipAndPs.first );
}
CrossSection Merger::MergingDSigDRBornGamma( ){
double weightCL = 0.;
weight = 1.;
if ( !currentNode()->children().empty() ) {
auto const inOutPair = currentNode()->getInOut();
// Here we make sure that clustering and splitting are invertible.
NodePtr randomChild = currentNode()->randomChild();
// Check if point is part of the ME region.
if( !matrixElementRegion( inOutPair.first ,
inOutPair.second ,
randomChild->pT() ,
theMergePt ) )weight *= 0.;
}
const NodePtr productionNode = currentNode()->getHistory( true , DSH()->hardScaleFactor() );
NodePtr CLNode;
NodePtr BornCL;
if( !currentNode()->children().empty() ){
if ( UseRandom::rndbool() ){
CLNode = currentNode()->randomChild();
bool inhist = CLNode->isInHistoryOf( productionNode );
weight *= inhist?( -2.*int( currentNode()->children().size() ) ):0.;
projected = true;
weightCL = 2.*int( currentNode()->children().size() );
BornCL = CLNode-> getHistory( false , DSH()->hardScaleFactor() );
}else{
weight = 2.;
projected = false;
}
}else{
weight = 1.;
projected = false;
}
if ( treefactory()->onlymulti() != -1&&
treefactory()->onlymulti() !=
int( currentNode()->legsize()-(projected ? 1 : 0) ) )
return ZERO;
if( !currentNode()->allAbove( mergePt()*(1.-1e-6) ) )weight = 0.;
if( CLNode&&!CLNode->allAbove( mergePt()*(1.-1e-6) ) )weightCL = 0.;
if ( !productionNode->xcomb()->willPassCuts() ){
return ZERO;
}
CrossSection res = ZERO;
bool maxMulti = currentNode()->legsize() == int( maxLegsLO() );
if( weight != 0. ){
Energy startscale = CKKW_StartScale( productionNode );
fillHistory( startscale , productionNode , currentNode() );
currentNode()->runningPt( fillProjector( (projected ? 1 : 0) ) );
weight *= history.back().weight*alphaReweight()*pdfReweight();
if( weight == 0.&&weightCL == 0. )return ZERO;
res += weight*TreedSigDR( startscale , ( !maxMulti&&!projected )?theGamma:1. );
}
if( CLNode&&theGamma != 1. ){
Energy startscale = CKKW_StartScale( BornCL );
fillHistory( startscale , BornCL , CLNode );
currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
weightCL *= history.back().weight*alphaReweight()*pdfReweight();
CrossSection diff = ZERO;
currentME()->factory()->setAlphaParameter( 1. );
diff -= weightCL*CLNode->dipole()->dSigHatDR( sqr( startscale* currentME()->renFac() ) );
currentME()->factory()->setAlphaParameter( theGamma );
string alp = ( CLNode->dipole()->aboveAlpha()?"Above":"Below" );
diff += weightCL*CLNode->dipole()->dSigHatDR( sqr( startscale* currentME()->renFac() ) );
currentME()->factory()->setAlphaParameter( 1. );
res += diff;
}
return res;
}
/// calculate the history weighted born cross section
#include "ThePEG/Handlers/SamplerBase.h"
CrossSection Merger::MergingDSigDRBornCheapME( ){
// Check if phase space poing is in ME region
if ( !currentNode()->children().empty() ) {
auto const & inOutPair = currentNode()->getInOut();
NodePtr randomChild = currentNode()->randomChild(); // Here we make sure that clustering and splitting are invertible
if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO;
}
CrossSection meweight=(SamplerBase::runLevel() == SamplerBase::RunMode)?TreedSigDR( 60.*GeV, 1. ):1.*nanobarn;
if (SamplerBase::runLevel() == SamplerBase::RunMode && theHighMeWeightMap.count(currentNode()) ) {
if (abs(meweight)>abs(theHighMeWeightMap[currentNode()])) {
theHighMeWeightMap[currentNode()]=abs(meweight);
}else{
if (meweight/theHighMeWeightMap[currentNode()]<UseRandom::rnd()) {
return ZERO;
}
}
}else{
theHighMeWeightMap.insert({currentNode(),abs(meweight)});
}
// get the history for the process
const NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() );
// decide if to cluster
weight = decideClustering(productionNode,currentNode(),projected);
if ( treefactory()->onlymulti() != -1&&
treefactory()->onlymulti() != int( currentNode()->legsize()-( projected ? 1 : 0 ) ) )
return ZERO;
// Check for cuts on the production proces.
if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
// calculate the staring scale
Energy startscale = CKKW_StartScale( productionNode );
// fill history with caluclation of sudakov supression
fillHistory( startscale , productionNode , currentNode() );
// fill the projector -> return the scale of the last splitting
currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
// the weight has three components to get shower history weight
weight *= history.back().weight*alphaReweight()*pdfReweight();
if( weight == 0. )return ZERO;
//calculate the cross section
return weight*TreedSigDR( startscale , 1. )*theHighMeWeightMap[currentNode()]/meweight;
}
CrossSection Merger::TreedSigDR( Energy startscale , double diffAlpha ){
currentME()->setScale( sqr( startscale ) , sqr( startscale ) );
CrossSection res = currentME()->dSigHatDRB();
if ( diffAlpha != 1. ) {
res += currentME()->dSigHatDRAlphaDiff( diffAlpha );
}
if( std::isnan( double( res/nanobarn ) ) ){
generator()->logWarning(Exception()
<< "Merger: TreedSigDR is nan"
<< Exception::warning);
res = ZERO;};
return res;
}
CrossSection Merger::LoopdSigDR( Energy startscale ){
currentME()->setScale( sqr( startscale ) , sqr( startscale ) );
currentME()->doOneLoopNoBorn();
CrossSection res = currentME()->dSigHatDRV()+
currentME()->dSigHatDRI();
currentME()->noOneLoopNoBorn();
return res;
}
Energy Merger::fillProjector( int pjs ){
// in the shower handler the scale is multiplied
// by DSH()->hardScaleFactor() so here we need
// to devide by the factor.
double xiQSh = history.begin()->node->legsize() == N0()?DSH()->hardScaleFactor():1.;
if( pjs == 0 ){
return ( history.size() == 1?1.:( 1./xiQSh ) )*history.back().scale;
}
for( auto const & hs : history )
if ( isProjectorStage( hs.node , pjs )&&pjs != 0 ){
currentNode()->xcomb()->lastProjector( hs.node->xcomb() );
return ( hs.node == history[0].node?1.:( 1./xiQSh ) )*hs.scale;
}
throw Exception() << "Could not fill projector." << Exception::abortnow;
return ZERO;
}
double Merger::pdfReweight(){ // TODO factorization scale inside
double res = 1.;
double facfactor = ( history[0].node->legsize() == N0()? currentME()->renFac():DSH()->facFac() );
for( int side : {0 , 1} ){
if( history[0].node->xcomb()->mePartonData()[side]->coloured() ){
for ( auto const & hs : history ){
//pdf-ratio only to the last step
if( !hs.node->parent() )continue;
if ( hs.node == history.back().node )continue;
if( !hs.node->dipole() ){
throw Exception() << "\nMerger: pdfReweight: history step has no dipol. "<< Exception::abortnow;
return 0.;
}
res *= pdfratio( hs.node , facfactor*hs.scale , DSH()->facFac()*( hs.node->pT() ) , side );
facfactor = DSH()->facFac();
}
res *= pdfratio( history.back().node , facfactor*history.back().scale , history[0].scale* currentME()->renFac() , side );
}
}
if ( std::isnan( res ) )
generator()->logWarning(Exception()
<< "Merger: pdfReweight is nan."
<< Exception::warning);
return res;
}
double Merger::cmwAlphaS(Energy q)const{
double alphas=1.;
using Constants::pi;
// No cmw-scheme
if (theCMWScheme==0) {
alphas = as( q );
}
// Linear cmw-scheme
else if(theCMWScheme==1){
alphas = as( q );
alphas *=1.+(3.*(67./18.-1./6.*sqr(pi))
-5./9.*Nf(q))* alphas/2./pi;
}
// cmw-scheme as factor in argument.
else if(theCMWScheme==2){
double kg=exp(-(67.-3.*sqr(pi)-10/3*Nf(q))
/( 2. *(33.-2.*Nf(q))));
//Note factor 2 since we here dealing with Energy
alphas = as(max(kg*q,1_GeV));
}else{
throw Exception()
<< "This CMW-Scheme is not implemented."
<< Exception::abortnow;
}
return alphas;
}
double Merger::alphaReweight(){
double res = 1.;
Energy Q_R = ( history[0].node->legsize() == N0()?
currentME()->renFac():
DSH()->renFac() )*
history[0].scale;
using Constants::pi;
const auto Q_qed=history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED();
const auto Oew=history[0].node->nodeME()->orderInAlphaEW();
const auto Oqcd=history[0].node->nodeME()->orderInAlphaS();
if (!history[0].node->children().empty()) {
assert(Oqcd!=0);
}
res *= pow( SM().alphaEMME( Q_qed )/ SM().alphaEMMZ() , Oew );
res *= pow( cmwAlphaS(Q_R) / SM().alphaS() , Oqcd );
for ( auto const & hs : history )
if ( hs.node!= history.back().node ){
Energy q_i = DSH()->renFac()* hs.node->pT();
res *= cmwAlphaS(q_i) / SM().alphaS();
}
if ( std::isnan( res ) )
generator()->logWarning(Exception()
<< "Merger: alphaReweight is nan. "<< Exception::warning);
return res;
}
void Merger::fillHistory( Energy scale , NodePtr begin , NodePtr endNode ){
history.clear();
double sudakov0_n = 1.;
history.push_back( {begin , sudakov0_n , scale} );
double xiQSh = history.begin()->node->legsize() == N0()?
DSH()->hardScaleFactor():1.;
scale *= xiQSh;
if ( begin->parent()||!isUnitarized ) {
while ( begin->parent() && ( begin != endNode ) ) {
if ( !dosudakov( begin , scale , begin->pT() , sudakov0_n ) ){
history.push_back( { begin->parent() , 0. , scale } );
}
if ( std::isnan( sudakov0_n ) )
generator()->logWarning(Exception()
<< "Merger: sudakov"<<scale/GeV<<" "
<<begin->pT()/GeV<<"0_n is nan. "
<< Exception::warning);
scale = begin->pT();
history.push_back( { begin->parent() , sudakov0_n , begin->pT() } );
begin = begin->parent();
}
Energy notunirunning = scale;
if ( !isUnitarized&&N()+N0() > int( currentNode()->legsize() ) ) {
if ( !dosudakov( begin , notunirunning , mergePt() , sudakov0_n ) ){
history.back().weight = 0.;
}else{
history.back().weight = sudakov0_n;
}
}
}
if( history.size() == 1 )scale /= DSH()->hardScaleFactor();
}
double Merger::sumPdfReweightExpansion()const{
double res = 0.;
Energy beam1Scale = history[0].scale*
( history[0].node->legsize() == N0()?
currentME()->renFac():
DSH()->facFac() );
Energy beam2Scale = history[0].scale*
( history[0].node->legsize() == N0()?
currentME()->renFac():
DSH()->facFac() );
for ( auto const & hs : history ){
//pdf expansion only to the last step
if( !hs.node->parent() )continue;
if( hs.node->xcomb()->mePartonData()[0]->coloured()&&
hs.node->nodeME()->lastX1()>0.99 )return 0.;
if( hs.node->xcomb()->mePartonData()[1]->coloured()&&
hs.node->nodeME()->lastX2()>0.99 )return 0.;
if( hs.node->nodeME()->lastX1()<0.00001 )return 0.;
if( hs.node->nodeME()->lastX2()<0.00001 )return 0.;
if ( hs.node->dipole()->bornEmitter() == 0 ){
res += pdfExpansion( hs.node , 0 ,
beam1Scale ,
( hs.node->pT() ) ,
hs.node->nodeME()->lastX1() ,
Nf( history[0].scale ) ,
history[0].scale );
beam1Scale = ( hs.node->pT() )*DSH()->facFac();
}
else if ( hs.node->dipole()->bornEmitter() == 1 ){
res += pdfExpansion( hs.node , 1 ,
beam2Scale ,
( hs.node->pT() ) ,
hs.node->nodeME()->lastX2() ,
Nf( history[0].scale ) ,
history[0].scale );
beam2Scale = ( hs.node->pT() )*DSH()->facFac();
}
// if we're here we know hs.node->dipole()->bornEmitter() > 1
// works only in collinear scheme
else if ( hs.node->dipole()->bornSpectator() == 0 ){
res += pdfExpansion( hs.node , 0 ,
beam1Scale ,
( hs.node->pT() ) ,
hs.node->nodeME()->lastX1() ,
Nf( history[0].scale ) ,
history[0].scale );
beam1Scale = ( hs.node->pT() )*DSH()->facFac();
}
else if ( hs.node->dipole()->bornSpectator() == 1 ){
res += pdfExpansion( hs.node , 1 ,
beam2Scale ,
( hs.node->pT() ) ,
hs.node->nodeME()->lastX2() ,
Nf( history[0].scale ) ,
history[0].scale );
beam2Scale = ( hs.node->pT() )*DSH()->facFac();
}
}
if ( currentNode()->xcomb()->mePartonData()[0]->coloured() ){
res += pdfExpansion( history.back().node , 0 ,
beam1Scale ,
history[0].scale* currentME()->renFac() ,
( history.back() ).node->nodeME()->lastX1() ,
Nf( history[0].scale ) ,
history[0].scale );
}
if ( currentNode()->xcomb()->mePartonData()[1]->coloured() ) {
res += pdfExpansion( history.back().node , 1 ,
beam2Scale ,
history[0].scale* currentME()->renFac() ,
( history.back() ).node->nodeME()->lastX2() ,
Nf( history[0].scale ) ,
history[0].scale );
}
return res;
}
#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
double Merger::pdfExpansion( NodePtr node ,
int side , Energy running ,
Energy next , double x ,
int nlp , Energy fixedScale ) const {
tcPDPtr particle , parton;
tcPDFPtr pdf;
if ( side == 0 ) {
particle = node->nodeME()->lastParticles().first->dataPtr();
parton = node->nodeME()->lastPartons().first->dataPtr();
pdf = node->xcomb()->partonBins().first->pdf();
}else{
assert( side == 1 );
particle = node->nodeME()->lastParticles().second->dataPtr();
parton = node->nodeME()->lastPartons().second->dataPtr();
pdf = node->xcomb()->partonBins().second->pdf();
}
//copied from PKOperator
double res = 0.;
int number = 10;
for ( int nr = 0;nr<number;nr++ ){
double restmp = 0;
using namespace RandomHelpers;
double r = UseRandom::rnd();
double eps = 1e-3;
pair<double , double> zw =
generate( ( piecewise() ,
flat( 0.0 , x ) ,
match( inverse( 0.0 , x , 1.0 ) +
inverse( 1.0+eps , x , 1.0 ) ) ) , r );
double z = zw.first;
double mapz = zw.second;
double PDFxparton = pdf->xfx( particle , parton , sqr( fixedScale ) , x )/x;
double CA = SM().Nc();
double CF = ( SM().Nc()*SM().Nc()-1.0 )/( 2.*SM().Nc() );
if ( abs( parton->id() ) < 7 ) {
double PDFxByzgluon = pdf->xfx( particle ,
getParticleData( ParticleID::g ) ,
sqr( fixedScale ) , x/z )*z/x;
double PDFxByzparton = pdf->xfx( particle , parton ,
sqr( fixedScale ) , x/z )*z/x;
assert( abs( parton->id() ) < 7 );
restmp += CF*( 3./2.+2.*log( 1.-x ) ) * PDFxparton;
if ( z > x ) {
restmp += 0.5 * ( sqr( z ) + sqr( 1.-z ) ) * PDFxByzgluon / z;
restmp += CF*2.*( PDFxByzparton - z*PDFxparton )/( z*( 1.-z ) );
restmp -= CF*PDFxByzparton * ( 1.+z )/z;
}
}else{
assert( parton->id() == ParticleID::g );
double PDFxByzgluon = pdf->xfx( particle ,
getParticleData( ParticleID::g ) ,
sqr( fixedScale ) , x/z )*z/x;
if ( z > x ){
double factor = CF * ( 1. + sqr( 1.-z ) ) / sqr( z );
for ( int f = -nlp; f <= nlp; ++f ) {
if ( f == 0 )
continue;
restmp += pdf->xfx( particle , getParticleData( f ) ,
sqr( fixedScale ) , x/z )*z/x*factor;
}
}
restmp += ( ( 11./6. ) * CA
- ( 1./3. )*Nf( history[0].scale )
+ 2.*CA*log( 1.-x ) ) *PDFxparton;
if ( z > x ) {
restmp += 2. * CA * ( PDFxByzgluon - z*PDFxparton ) / ( z*( 1.-z ) );
restmp += 2.* CA *( ( 1.-z )/z - 1. + z*( 1.-z ) ) * PDFxByzgluon / z;
}
}
if ( PDFxparton<1e-8 )restmp = 0.;
res += 1*restmp*log( sqr( running/next ) )/PDFxparton*mapz;
}
return res/number;
}
double Merger::sumAlphaSReweightExpansion()const{
double res = 0.;
if ( !( history[0].node->children().empty() ) ){
res += alphasExpansion( history[0].scale* DSH()->renFac() ,
history[0].scale* currentME()->renFac() )*
int( history[0].node->legsize()-N0() );
}
// dsig is calculated with fixed alpha_s
for ( auto const & hs : history ){
//expansion only to the last step
if( !hs.node->parent() )continue;
res += alphasExpansion( hs.node->pT()*DSH()->renFac() , history[0].scale );
}
return res;
}
double Merger::sumFillHistoryExpansion(){
double res = 0.;
double xiQSh = history[0].node->legsize() == N0()?DSH()->hardScaleFactor():1.;
for ( auto const & hs : history ){
if( !hs.node->parent() )continue;
doHistExpansion( hs.node , ( hs.node == history[0].node?xiQSh:1. )*hs.scale ,
hs.node->pT() , history[0].scale , res );
}
return res;
}
MergingFactoryPtr Merger::treefactory()const{return theTreeFactory;}
void Merger::doinit(){
if ( !DSH()->hardScaleIsMuF() ) {
throw Exception()
<< "Merger: Merging is currently only sensible "
<< "if we are using the hardScale as MuF."
<< Exception::abortnow;
}
}
CrossSection Merger::MergingDSigDR() {
history.clear();
assert(currentNode()==theFirstNodeMap[ currentME()]);
if(DSH()->doesSplitHardProcess()){
throw Exception()
<< "Merger: The splithardprocess option is currently not supported."
<< Exception::abortnow;
}
DSH()->eventHandler( generator()->eventHandler() );
CrossSection res = ZERO;
if( currentNode()->subtractedReal() ){
res = MergingDSigDRRealStandard();
theCurrentMaxLegs = maxLegsNLO();
}else if( currentNode()->virtualContribution() ){
res = MergingDSigDRVirtualStandard();
theCurrentMaxLegs = maxLegsNLO();
}else if( theGamma!= 1. ){
res = MergingDSigDRBornGamma();
theCurrentMaxLegs = maxLegsLO();
}else{
res = MergingDSigDRBornStandard();
theCurrentMaxLegs = maxLegsLO();
}
auto lxc= currentME()->lastXCombPtr();
lxc->lastShowerScale( sqr( currentNode()->runningPt() ) );
auto lp= currentME()->lastXCombPtr()->lastProjector();
if( lp ){
lp->lastShowerScale( sqr( currentNode()->runningPt() ) );
}
if ( res == ZERO ){
history.clear();
return ZERO;
}
cleanup( currentNode() );
lxc->subProcess( SubProPtr() );
history.clear();
if( !std::isfinite( double( res/nanobarn ) ) ){
generator()->logWarning(Exception()
<< "Merger weight is " << res/nanobarn<< " -> setting to 0"
<< Exception::warning);
return ZERO;
}
return res;
}
#include "Herwig/PDF/HwRemDecayer.h"
void Merger::CKKW_PrepareSudakov( NodePtr node , Energy running ){
//cleanup( node );
tSubProPtr sub = node->xcomb()->construct();
const auto pb=node->xcomb()->partonBins();
DSH()->resetPDFs( {pb.first->pdf(),pb.second->pdf() },pb );
DSH()->setCurrentHandler();
DSH()->currentHandler()->generator()->currentEventHandler(
currentNode()->xcomb()->eventHandlerPtr() );
DSH()->currentHandler()->remnantDecayer()->setHadronContent( currentNode()->xcomb()->lastParticles() );
DSH()->eventRecord().clear();
DSH()->eventRecord().slimprepare( sub , dynamic_ptr_cast<tStdXCombPtr>( node->xcomb() ) , DSH()->pdfs() , currentNode()->xcomb()->lastParticles() );
DSH()->hardScales( sqr( running ) );
}
Energy Merger::CKKW_StartScale( NodePtr node ) const {
Energy res = generator()->maximumCMEnergy();;
if( !node->children().empty() ){
for ( int i = 0;i<node->legsize();i++ ){
if ( !node->nodeME()->mePartonData()[i]->coloured() )continue;
for ( int j = i+1;j<node->legsize();j++ ){
if ( i == j||!node->nodeME()->mePartonData()[j]->coloured() )continue;
res = min( res , sqrt( 2.*node->nodeME()->lastMEMomenta()[i]*node->nodeME()->lastMEMomenta()[j] ) );
}
}
}else{
node->nodeME()->factory()->scaleChoice()->setXComb( node->xcomb() );
res = sqrt( node->nodeME()->factory()->scaleChoice()->renormalizationScale() );
}
node->nodeME()->factory()->scaleChoice()->setXComb( node->xcomb() );
res = max( res , sqrt( node->nodeME()->factory()->scaleChoice()->renormalizationScale() ) );
return res;
}
double Merger::alphasExpansion( Energy next , Energy fixedScale ) const {
double betaZero = ( 11./6. )*SM().Nc() - ( 1./3. )*Nf( history[0].scale );
// TODO factorize the return statement, turn ?: into if
return ( betaZero*log( sqr( fixedScale/next ) ) )+
( theCMWScheme>0?( 3.*( 67./18.-1./6.*Constants::pi*Constants::pi )
-5./9.*Nf( history[0].scale ) ):0. );
}
double Merger::pdfratio( NodePtr node , Energy numerator_scale , Energy denominator_scale , int side ){
StdXCombPtr bXc = node->xcomb();
if( !bXc->mePartonData()[side]->coloured() )
throw Exception()
<< "Merger: pdf-ratio required for non-coloured particle."
<< Exception::abortnow;
double from , to;
if ( side == 0 ){
if ( denominator_scale == numerator_scale ) {
return 1.;
}
from = node->nodeME()->pdf1( sqr( denominator_scale ) );
to = node->nodeME()->pdf1( sqr( numerator_scale ) );
if ( ( to < 1e-8||from < 1e-8 )&&( to/from>10000000. ) ){
generator()->logWarning(Exception()
<< "Merger: pdfratio to = " << to << " from = " << from
<< Exception::warning);
return 0.;
}
}
else{
if ( denominator_scale == numerator_scale ) {
return 1.;
}
from = node->nodeME()->pdf2( sqr( denominator_scale ) );
to = node->nodeME()->pdf2( sqr( numerator_scale ) );
if ( ( to < 1e-8||from < 1e-8 )&&( to/from>10000000. ) ){
generator()->logWarning(Exception()
<< "Merger: pdfratio to = " << to << " from = " << from
<< Exception::warning);
return 0.;}
}
return to/from;
}
bool Merger::dosudakov( NodePtr node , Energy running , Energy next , double& sudakov0_n ) {
CKKW_PrepareSudakov( node , running );
for( DipoleChain const & chain : DSH()->eventRecord().chains() ){
for( Dipole const & dip : chain.dipoles() ){
sudakov0_n *= singlesudakov( dip , next , running , { true , false } );
sudakov0_n *= singlesudakov( dip , next , running , { false , true } );
if ( sudakov0_n == 0.0 ){
cleanup( node );
return false;
}
}
}
cleanup( node );
return true;
}
bool Merger::doHistExpansion( NodePtr node , Energy running , Energy next , Energy fixedScale , double& histExpansion ) {
CKKW_PrepareSudakov( node , running );
for( DipoleChain const & chain : DSH()->eventRecord().chains() ){
for( Dipole const & dip : chain.dipoles() ){
histExpansion += singleHistExpansion( dip , next , running , fixedScale , { true , false } );;
histExpansion += singleHistExpansion( dip , next , running , fixedScale , { false , true } );
}
}
cleanup( node );
return true;
}
bool Merger::isProjectorStage( NodePtr node , int pjs )const{
return ( pjs == int( ( currentNode()->legsize() - node->legsize() ) ) );
}
void Merger::cleanup( NodePtr node ) {
DSH()->eventRecord().clear();
if( !node->xcomb()->subProcess() )return;
ParticleVector vecfirst = node->xcomb()->subProcess()->incoming().first->children();
for( auto const & particle : vecfirst )
node->xcomb()->subProcess()->incoming().first->abandonChild( particle );
ParticleVector vecsecond = node->xcomb()->subProcess()->incoming().second->children();
for( auto const & particle : vecsecond )
node->xcomb()->subProcess()->incoming().second->abandonChild( particle );
node->xcomb()->subProcess( SubProPtr() );
}
double Merger::singlesudakov( Dipole dip , Energy next , Energy running , pair<bool , bool> conf ){
double res = 1.;
tPPtr emitter = dip.emitter( conf );
tPPtr spectator = dip.spectator( conf );
DipoleSplittingInfo candidate( dip.index( conf ) , conf ,
dip.emitterX( conf ) ,
dip.spectatorX( conf ) ,
emitter , spectator );
if ( DSH()->generators().find( candidate.index() ) == DSH()->generators().end() ) DSH()->getGenerators( candidate.index() );
auto const & gens = DSH()->generators().equal_range( candidate.index() );
for ( auto gen = gens.first; gen != gens.second; ++gen ) {
if ( !( gen->first == candidate.index() ) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale( emitter->momentum() , spectator->momentum() );
candidate.scale( dScale );
candidate.continuesEvolving();
Energy ptMax = gen->second->splittingKinematics()->ptMax( candidate.scale() , candidate.emitterX() , candidate.spectatorX() ,
candidate.index() , *gen->second->splittingKernel() );
candidate.hardPt( min( running , ptMax ) );
if ( candidate.hardPt()>next ){
res *= gen->second->sudakov( candidate , next );
}
}
return res;
}
double Merger::singleHistExpansion( Dipole dip , Energy next , Energy running , Energy fixedScale , pair<bool , bool> conf ){
double res = 0.;
tPPtr emitter = dip.emitter( conf );
tPPtr spectator = dip.spectator( conf );
DipoleSplittingInfo candidate( dip.index( conf ) ,
conf , dip.emitterX( conf ) ,
dip.spectatorX( conf ) ,
emitter , spectator );
if ( DSH()->generators().find( candidate.index() ) ==
DSH()->generators().end() )
DSH()->getGenerators( candidate.index() );
auto const & gens = DSH()->generators().equal_range( candidate.index() );
for ( auto gen = gens.first; gen != gens.second; ++gen ) {
if ( !( gen->first == candidate.index() ) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(
emitter->momentum() , spectator->momentum() );
candidate.scale( dScale );
candidate.continuesEvolving();
Energy ptMax = gen->second->
splittingKinematics()->ptMax(
candidate.scale() , candidate.emitterX() ,
candidate.spectatorX() , candidate.index() ,
*gen->second->splittingKernel() );
candidate.hardPt( min( running , ptMax ) );
if ( candidate.hardPt()>next ){
res += gen->second->sudakovExpansion( candidate , next , fixedScale );
}
}
return res;
}
void Merger::firstNodeMap( MatchboxMEBasePtr a , NodePtr b ){
theFirstNodeMap.insert( { a , b } );
}
map<MatchboxMEBasePtr , NodePtr> Merger::firstNodeMap()const{return theFirstNodeMap;}
void Merger::setXComb( tStdXCombPtr xc ){
currentNode()->setXComb( xc );
}
void Merger::setKinematics( ){
currentNode()->setKinematics();
}
void Merger::clearKinematics( ){
currentNode()->clearKinematics();
}
void Merger::flushCaches(){
if (currentNode()&&currentNode()->xcomb()->lastParticles().first) {
currentNode()->flushCaches();
}
}
bool Merger::generateKinematics( const double * r ){
return currentNode()->firstgenerateKinematics( r , ! currentNode()->subtractedReal() );
}
bool Merger::matrixElementRegion( PVector incoming , PVector outgoing , Energy winnerScale , Energy cutscale )const{
Energy ptx = Constants::MaxEnergy;
bool foundwinnerpt = false;
using namespace boost;
//FF
for( auto const & em : outgoing ){ if ( ! em->coloured() ) continue;
for( auto const & emm : outgoing ){ if ( !emm->coloured() ) continue; if ( em == emm ) continue;
for( auto const & spe : outgoing ){ if ( !spe->coloured() ) continue; if ( em == spe||emm == spe ) continue;
if ( !( em->id() == -emm->id()||emm->id()>6 ) )continue;
Energy pt = ZERO;
if ( em->momentum().m()<= 10_MeV &&
emm->momentum().m()<= 10_MeV &&
spe->momentum().m()<= 10_MeV ) {
pt = FFLTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
}else{
pt = FFMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
}
if ( abs( pt-winnerScale ) < 10_MeV ) {
foundwinnerpt = true;
}
ptx = min( ptx , pt );
}
}
}
//FI
for( auto const & spe : incoming ){ if ( ! spe->coloured() ) continue;
for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue;
for( auto const & em : outgoing ){ if ( ! em->coloured() ) continue; if ( em == emm ) continue;
if ( !( em->id() == -emm->id() || emm->id()>6 ) )continue;
Energy pt = ZERO;
if ( em->momentum().m()<= 10_MeV &&
emm->momentum().m()<= 10_MeV &&
spe->momentum().m()<= 10_MeV ) {
pt = FILTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
}else{
pt = FIMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
}
if ( abs( pt-winnerScale )<10_MeV ) {
foundwinnerpt = true;
}
if( pt > ZERO )
ptx = min( ptx , pt );
}
}
}
//IF
for( auto const & em : incoming ){ if ( ! em->coloured() ) continue;
for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue;
for( auto const & spe : outgoing ){ if ( ! spe->coloured() ) continue; if ( emm == spe ) continue;
if ( !( em->id()>6|| em->id() == emm->id() ||emm->id()>6 ) )continue;
Energy pt = ZERO;
if ( em->momentum().m()<= 10_MeV &&
emm->momentum().m()<= 10_MeV &&
spe->momentum().m()<= 10_MeV ) {
//massless
pt = IFLTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
}else{
//massiv
pt = IFMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
}
if ( abs( pt-winnerScale )< 10_MeV ) {
foundwinnerpt = true;
}
ptx = min( ptx , pt );
}
}
}
//II
for( auto const & em : incoming ){ if ( ! em->coloured() ) continue;
for( auto const & spe : incoming ){ if ( ! spe->coloured() ) continue; if ( em == spe ) continue;
for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue;
if ( !( em->id()>6||em->id() == emm->id() ||emm->id()>6 ) )continue;
Energy pt = IILTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
if ( abs( pt-winnerScale )< 10_MeV ) {
foundwinnerpt = true;
}
ptx = min( ptx , pt );
}
}
}
if( !foundwinnerpt ){
generator()->logWarning( Exception()
<< "Merger: Could not find winner with pt."
- << "Run with -d2 to get phase space points. "
+ << "Run with -d3 to get phase space points. "
<< Exception::warning );
if( Debug::level > 2 ) {
generator()->log() << "\nWarning: could not find winner with pt: " << winnerScale/GeV;
for( auto const & emm : incoming ){
generator()->log() << "\n" << emm->id() << " " << emm->momentum()/GeV << " " << emm->momentum().m()/GeV << flush;
}
for( auto const & emm : outgoing ){
generator()->log() <<"\n" << emm->id() << " " << emm->momentum()/GeV << " " << emm->momentum().m()/GeV << flush;
}
}
}
return ( ptx>cutscale ) ;
}
bool Merger::notOnlyMulti() const {
return ( treefactory()->onlymulti() != -1&&
treefactory()->onlymulti() !=
int( currentNode()->legsize()-( projected ? 1 : 0 ) ) );
}
int Merger::M()const{return theTreeFactory->M();}
int Merger::N()const{return theTreeFactory->N();}
void Merger::debugVirt(double weight,double w1,double w2,double w3,
CrossSection matrixElement,double ww1,double ww2,
double ww3, NodePtr node,CrossSection bornWeight)const{
Energy minPT = Constants::MaxEnergy;
for( auto const & no :currentNode()->children() )minPT = min( no->pT() , minPT );
generator()->log() << "\nVIRT " << minPT/GeV << " " << weight << " " << w1;
generator()->log() << " " << w2;
generator()->log() << " " << w3;
generator()->log() << " " << ( matrixElement/nanobarn ) << " " << ww1 << " " << ww2 << " " << ww3 << " " << node->pT()/GeV << " " << node->nodeME()->mePartonData()[3]->mass()/GeV << " " << ( bornWeight*SM().alphaS()/( 2.*ThePEG::Constants::pi )/nanobarn );
}
void Merger::debugReal(string realstr, double weight,
CrossSection me, CrossSection dip)const {
Energy minPT = Constants::MaxEnergy;
for( auto const & no :currentNode()->children() )minPT = min( no->pT() , minPT );
generator()->log() << "\n"<<realstr <<" "<< minPT/GeV << " " << weight << " " << ( me-dip )/nanobarn << " " << me/nanobarn << " " << dip/nanobarn;
}
// If needed , insert default implementations of virtual function defined
// in the InterfacedBase class here ( using ThePEG-interfaced-impl in Emacs ).
#include "ThePEG/Persistency/PersistentOStream.h"
void Merger::persistentOutput( PersistentOStream & os ) const {
os << theShowerExpansionWeights << theCMWScheme << projected <<
isUnitarized << isNLOUnitarized <<
theOpenInitialStateZ << theChooseHistory <<
theN0 << theOnlyN << weight <<
weightCB << theGamma << theSmearing << ounit( theIRSafePT , GeV ) <<
ounit( theMergePt , GeV ) << ounit( theCentralMergePt , GeV )
<< theLargeNBasis << FFLTK << FILTK <<
IFLTK << IILTK << FFMTK << FIMTK << IFMTK <<
theDipoleShowerHandler << theTreeFactory << theFirstNodeMap<<theRealSubtractionRatio;
}
#include "ThePEG/Persistency/PersistentIStream.h"
void Merger::persistentInput( PersistentIStream & is , int ) {
is >> theShowerExpansionWeights >> theCMWScheme >> projected >>
isUnitarized >> isNLOUnitarized >>
theOpenInitialStateZ >>
theChooseHistory >> theN0 >> theOnlyN >>
weight >> weightCB >>
theGamma >> theSmearing >> iunit( theIRSafePT , GeV ) >>
iunit( theMergePt , GeV ) >>
iunit( theCentralMergePt , GeV ) >> theLargeNBasis >>
FFLTK >> FILTK >> IFLTK >>
IILTK >> FFMTK >> FIMTK >>
IFMTK >> theDipoleShowerHandler >>
theTreeFactory >> theFirstNodeMap>>theRealSubtractionRatio ;
}
// *** 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 ).
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<Merger , Herwig::MergerBase>
describeHerwigMerger( "Herwig::Merger" , "HwDipoleShower.so" );
//ClassDescription<Merger> Merger::initMerger;
// Definition of the static class description member.
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
void Merger::Init() {
static ClassDocumentation<Merger> documentation
( "The Merger class takes care of merging multiple LO & NLO cross sections." );
static Reference<Merger , DipoleShowerHandler> interfaceShowerHandler
( "DipoleShowerHandler" ,
"Set the pointer to the Dipole Shower Handler" ,
&Merger::theDipoleShowerHandler ,
false , false ,
true , true , false );
static Switch<Merger , bool>
interfaceShowerExpansionWeights
( "ShowerExpansionWeights" ,
"Calculate the expansions of the shower history to be NLO accurate" ,
&Merger::theShowerExpansionWeights ,
true , false , false );
static SwitchOption
interfaceShowerExpansionWeightsYes
( interfaceShowerExpansionWeights ,
"Yes" ,
"Calculate shower expansions" ,
true );
static SwitchOption
interfaceShowerExpansionWeightsNo
( interfaceShowerExpansionWeights ,
"No" ,
"Do not calculate expansions" ,
false );
static Switch<Merger , unsigned int>
interfacetheCMWScheme
( "CMWScheme" ,
"Use CMW-Scheme to calculate the alpha_s for the shower expressions." ,
&Merger::theCMWScheme ,
0 ,
false , false );
static SwitchOption interfacetheCMWSchemeOff
(interfacetheCMWScheme,
"Off",
"No CMW-Scheme",
0);
static SwitchOption interfacetheCMWSchemeLinear
(interfacetheCMWScheme,
"Linear",
"Linear CMW multiplication: alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi )",
1);
static SwitchOption interfacetheCMWSchemeFactor
(interfacetheCMWScheme,
"Factor",
"Use factor in alpha_s argument: alpha_s(q) -> alpha_s(fac*q) with fac=exp(-(67-3pi^2-10/3*Nf)/(33-2Nf)) ",
2);
static Parameter<Merger , Energy> interfaceMergerScale
( "MergingScale" ,
"The MergingScale." ,
&Merger::theCentralMergePt ,
GeV , 20.0*GeV , 0.0*GeV , 0*GeV ,
false , false , Interface::lowerlim );
static Reference<Merger , MergingFactory> interfaceMergerHelper
( "MergingFactory" ,
"Set a pointer to the MergingFactory" ,
&Merger::theTreeFactory ,
false , false ,
true , true , false );
static Parameter<Merger , double> interfacegamma
( "gamma" ,
"gamma parameter." ,
&Merger::theGamma , 1.0 , 0.0 , 0 ,
false , false , Interface::lowerlim );
static Reference<Merger , ColourBasis> interfaceLargeNBasis
( "LargeNBasis" ,
"Set the large-N colour basis implementation." ,
&Merger::theLargeNBasis ,
false , false ,
true , true , false );
static Switch<Merger , bool>
interfaceOpenInitialSateZ
( "OpenInitialStateZ" , "" , &Merger::theOpenInitialStateZ , false , false , false );
static SwitchOption interfaceOpenInitialSateZYes
( interfaceOpenInitialSateZ , "Yes" , "" , true );
static SwitchOption interfaceOpenInitialSateZNo
( interfaceOpenInitialSateZ , "No" , "" , false );
static Parameter<Merger , Energy>
interfaceIRSafePT
( "IRSafePT" ,
"Set the pt for which a matrixelement should be treated as IR-safe." ,
&Merger::theIRSafePT ,
GeV , 0.0 * GeV , ZERO , Constants::MaxEnergy ,
true , false , Interface::limited );
interfaceIRSafePT.setHasDefault( false );
static Parameter<Merger , double>
interfacemergePtsmearing(
"MergingScaleSmearing" ,
"Set the percentage the merging pt should be smeared." ,
&Merger::theSmearing ,
0. , 0. , 0.0 , 0.5 ,
true , false , Interface::limited );
static Parameter<Merger , int>
interfacechooseHistory(
"chooseHistory" ,
"different ways of choosing the history" ,
&Merger::theChooseHistory ,
3 , -1 , 0 ,
false , false , Interface::lowerlim );
static Parameter<Merger , double>
interfacetheRealSubtractionRatio(
"RealSubtractionRatio" ,
"Set the percentage the merging pt should be smeared." ,
&Merger::theRealSubtractionRatio , 0. , 0. ,
1. , 10.0 ,
true , false , Interface::limited );
}
diff --git a/Shower/Dipole/Merging/Node.cc b/Shower/Dipole/Merging/Node.cc
--- a/Shower/Dipole/Merging/Node.cc
+++ b/Shower/Dipole/Merging/Node.cc
@@ -1,454 +1,461 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Node class.
//
#include "Node.h"
#include "MergingFactory.h"
#include "Merger.h"
using namespace Herwig;
Node::Node(MatchboxMEBasePtr nodeME, int cutstage, MergerPtr mh)
:Interfaced(),
thenodeMEPtr(nodeME),
thedipol(),
theparent(),
theCutStage(cutstage),
isOrdered(true),
theSubtractedReal(false),
theVirtualContribution(false),
theMergingHelper(mh)
{
nodeME->maxMultCKKW(1);
nodeME->minMultCKKW(0);
}
Node::Node(NodePtr deephead,
NodePtr head,
SubtractionDipolePtr dipol,
MatchboxMEBasePtr nodeME,
int cutstage)
:Interfaced(), thenodeMEPtr(nodeME),
thedipol(dipol),
theparent(head),
theDeepHead(deephead),
theCutStage(cutstage),
isOrdered(true),
theSubtractedReal(false),
theVirtualContribution(false),
theMergingHelper() //The subnodes have no merging helper
{
}
Node::~Node() { }
SubtractionDipolePtr Node::dipole() const {
return thedipol;
}
/** returns the matrix element pointer */
const MatchboxMEBasePtr Node::nodeME() const {
return thenodeMEPtr;
}
/** access the matrix element pointer */
MatchboxMEBasePtr Node::nodeME() {
return thenodeMEPtr;
}
pair<PVector , PVector> Node::getInOut( ){
PVector in;
const auto me= nodeME();
const auto pd=me->mePartonData();
for( auto i : {0 , 1} )
in.push_back(pd[i]->produceParticle( me->lastMEMomenta()[i] ) );
PVector out;
for ( size_t i = 2;i< pd.size();i++ ){
PPtr p = pd[i]->produceParticle( me->lastMEMomenta()[i] );
out.push_back( p );
}
return { in , out };
}
int Node::legsize() const {return nodeME()->legsize();}
NodePtr Node::randomChild() {
return thechildren[UseRandom::irnd(thechildren.size())];
}
bool Node::allAbove(Energy pt) {
for (NodePtr child : thechildren)
if ( child->pT() < pt )
return false;
return true;
}
bool Node::isInHistoryOf(NodePtr other) {
while (other->parent()) {
if (other == this)
return true;
other = other->parent();
}
return false;
}
void Node::flushCaches() {
+ if (didflush) return;
+ didflush=true;
for ( auto const & ch: thechildren) {
ch->xcomb()->clean();
ch->nodeME()->flushCaches();
ch->flushCaches();
}
}
void Node::setKinematics() {
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
ch->dipole()->setKinematics();
ch->nodeME()->setKinematics();
ch->setKinematics();
}
}
void Node::clearKinematics() {
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
ch->nodeME()->clearKinematics();
ch->dipole()->clearKinematics();
ch->clearKinematics();
}
}
bool Node::generateKinematics(const double *r, bool directCut) {
+ didflush=false;
// If there are no children to the child process we are done.
if(children().empty()) return true;
assert(parent());
if ( ! directCut && pT() < deepHead()->MH()->mergePt()) {
// Real emission:
// If there are children to the child process,
// we now require that all subsequent children
// with pt < merging scale are in their ME region.
// Since the possible children of the real emission
// contribution are now in their ME region,
// it is clear that a second clustering is possible
// -- modulo phase space restrictions.
// Therefore the real emission contribution
// are unitarised and the cross section is
// hardly modified.
auto inOutPair = getInOut();
NodePtr rc = randomChild();
rc->dipole()->setXComb(rc->xcomb());
if(!rc->dipole()->generateKinematics(r))assert(false);
// If not in ME -> return false
if(!deepHead()->MH()->matrixElementRegion( inOutPair.first ,
inOutPair.second ,
rc->pT() ,
deepHead()->MH()->mergePt() ) )return false;
}
for (auto & ch : children() ) {
ch->dipole()->setXComb(ch->xcomb());
if ( !ch->dipole()->generateKinematics(r) ) { assert(false); }
ch->generateKinematics( r, true);
}
return true;
}
bool Node::firstgenerateKinematics(const double *r, bool directCut) {
+ didflush=false;
// This is called form the merging helper for the first node. So:
assert(!parent());
assert(xcomb());
+
+ ///// This should not be needed!!!
+ ///// ( Warning inMerger::matrixElementRegion gets triggered.)
flushCaches();
//Set here the new merge Pt for the next phase space point.( Smearing!!!)
MH()->smearMergePt();
// If there are no children to this node, we are done here:
if (children().empty())
return true;
// directCut is for born and for virtual contributions.
// if directCut is true, then cut on the first ME region.
// call recursiv generate kinematics for subsequent nodes.
if ( directCut ){
auto inOutPair = getInOut();
NodePtr rc = randomChild();
rc->dipole()->setXComb(rc->xcomb());
if ( !rc->dipole()->generateKinematics(r) ) { assert(false); }
rc->nodeME()->setXComb(rc->xcomb());
if(MH()->gamma() == 1.){
if(!MH()->matrixElementRegion( inOutPair.first ,
inOutPair.second ,
rc->pT() ,
MH()->mergePt() ) ){
return false;
}
}else{
// Different treatment if gamma is not 1.
// Since the dipoles need to be calculated always
// their alpha region is touched.
// AlphaRegion != MERegion !!!
bool inAlphaPS = false;
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
if ( !ch->dipole()->generateKinematics(r) )
assert(false);
MH()->treefactory()->setAlphaParameter( MH()->gamma() );
inAlphaPS |= ch->dipole()->aboveAlpha();
MH()->treefactory()->setAlphaParameter( 1. );
}
NodePtr rc = randomChild();
if(!inAlphaPS&&
!MH()->matrixElementRegion( inOutPair.first ,
inOutPair.second ,
rc->pT() ,
MH()->mergePt() ) )
return false;
}
}
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
if ( !ch->dipole()->generateKinematics(r) ) { assert(false); }
if( ! ch->generateKinematics(r,directCut) )return false;
}
return true;
}
StdXCombPtr Node::xcomb() const {
assert(thexcomb);
return thexcomb;
}
StdXCombPtr Node::xcomb(){
if(thexcomb)return thexcomb;
assert(parent());
thexcomb=dipole()->makeBornXComb(parent()->xcomb());
xcomb()->head(parent()->xcomb());
dipole()->setXComb(thexcomb);
return thexcomb;
}
void Node::setXComb(tStdXCombPtr xc) {
assert ( !parent() );
thexcomb=xc;
assert(thexcomb->lastParticles().first);
}
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
void Node::birth(const vector<MatchboxMEBasePtr> & vec) {
// produce the children
vector<SubtractionDipolePtr> dipoles =
nodeME()->getDipoles(DipoleRepository::dipoles(
nodeME()->factory()->dipoleSet()), vec, true);
for ( auto const & dip : dipoles ) {
dip->doSubtraction();
NodePtr node = new_ptr(Node(theDeepHead,
this,
dip,
dip->underlyingBornME(),
theDeepHead->cutStage()));
thechildren.push_back(node);
}
}
vector<NodePtr> Node::getNextOrderedNodes(bool normal, double hardScaleFactor) const {
vector<NodePtr> temp = children();
vector<NodePtr> res;
for (NodePtr const & child : children()) {
if(deepHead()->MH()->mergePt()>child->pT()) {
res.clear();
return res;
}
}
for (NodePtr const & child: children()) {
if (parent()&& normal) {
if ( child->pT() < pT() ) {
continue;
}
}
if ( child->children().size() != 0 ) {
for (NodePtr itChild: child->children()) {
if( itChild->pT() > child->pT()&&child->inShowerPS(itChild->pT()) ) {
res.push_back(child);
break;
}
}
}
else {
const auto sc=child->nodeME()->factory()->scaleChoice();
sc->setXComb(child->xcomb());
if ( sqr(hardScaleFactor)*
sc->renormalizationScale() >= sqr(child->pT()) &&
child->inShowerPS(hardScaleFactor*sqrt(sc->renormalizationScale()))) {
res.push_back(child);
}
}
}
return res;
}
bool Node::inShowerPS(Energy hardpT)const {
// Here we decide if the current phase space
// point can be reached from the underlying Node.
double z_ = dipole()->lastZ();
// II
if( dipole()->bornEmitter()<2&&
dipole()->bornSpectator()<2&&
deepHead()->MH()->openInitialStateZ()) return true;
// IF
if( dipole()->bornEmitter()<2&&
dipole()->bornSpectator() >= 2&&
deepHead()->MH()->openInitialStateZ())
return true;
pair<double, double> zbounds =
dipole()->tildeKinematics()->zBounds(pT(), hardpT);
return (zbounds.first<z_&&z_<zbounds.second);
}
NodePtr Node::getHistory(bool normal, double hardScaleFactor) {
NodePtr res = this;
vector<NodePtr> temp = getNextOrderedNodes(normal, hardScaleFactor);
Energy minpt = Constants::MaxEnergy;
Selector<NodePtr> subprosel;
while (temp.size() != 0) {
minpt = Constants::MaxEnergy;
subprosel.clear();
for (NodePtr const & child : temp) {
if( child->dipole()->underlyingBornME()->largeNColourCorrelatedME2(
{child->dipole()->bornEmitter(),
child->dipole()->bornSpectator()},
deepHead()->MH()->largeNBasis()) != 0.
) {
double weight = abs(child->dipole()->dSigHatDR()/nanobarn);
if(weight != 0.) {
subprosel.insert(weight , child);
minpt = min(minpt, child->pT());
}
//TODO choosehistories
}
}
if (subprosel.empty())
return res;
res = subprosel.select(UseRandom::rnd());
temp = res->getNextOrderedNodes(true, hardScaleFactor);
}
return res;
}
pair<CrossSection, CrossSection> Node::calcDipandPS(Energy scale)const {
return dipole()->dipandPs(sqr(scale), deepHead()->MH()->largeNBasis());
}
CrossSection Node::calcPs(Energy scale)const {
return dipole()->ps(sqr(scale), deepHead()->MH()->largeNBasis());
}
CrossSection Node::calcDip(Energy scale)const {
return dipole()->dip(sqr(scale));
}
IBPtr Node::clone() const {
return new_ptr(*this);
}
IBPtr Node::fullclone() const {
return new_ptr(*this);
}
#include "ThePEG/Persistency/PersistentOStream.h"
void Node::persistentOutput(PersistentOStream & os) const {
os <<
thexcomb<<
thenodeMEPtr<<
thedipol<<
thechildren<<
theparent<<
theProjector<<
theDeepHead<<
theCutStage<<
ounit(theRunningPt, GeV)<<
theSubtractedReal<<
theVirtualContribution<<
theMergingHelper;
}
#include "ThePEG/Persistency/PersistentIStream.h"
void Node::persistentInput(PersistentIStream & is, int) {
is >>
thexcomb>>
thenodeMEPtr>>
thedipol>>
thechildren>>
theparent>>
theProjector>>
theDeepHead>>
theCutStage>>
iunit(theRunningPt, GeV)>>
theSubtractedReal>>
theVirtualContribution>>
theMergingHelper;
}
// *** 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).
#include "ThePEG/Utilities/DescribeClass.h"
DescribeClass<Node, Interfaced>
describeHerwigNode("Herwig::Node", "HwDipoleShower.so");
void Node::Init() {
static ClassDocumentation<Node>
documentation("There is no documentation for the Node class");
}
diff --git a/Shower/Dipole/Merging/Node.h b/Shower/Dipole/Merging/Node.h
--- a/Shower/Dipole/Merging/Node.h
+++ b/Shower/Dipole/Merging/Node.h
@@ -1,238 +1,235 @@
// -*- C++ -*-
#ifndef Herwig_Node_H
#define Herwig_Node_H
//
// This is the declaration of the Node class.
//
#include "Node.fh"
#include "MergingFactory.fh"
#include "Merger.h"
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Config/std.h"
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include <vector>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the Node class.
*
* @see \ref NodeInterfaces "The interfaces"
* defined for Node.
*/
- /**
- * Define the SafeClusterMap type map< tuple<emitter, emission, spectator >
- * , pair<first-clustering, second-clustering> >
- */
- typedef map< std::tuple<int, int, int> , pair<bool, bool> > SafeClusterMap;
class Node : public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
Node (){};
// constructor for first nodes
Node(MatchboxMEBasePtr nodeME , int cutstage , MergerPtr mh );
// another constructor for underlying nodes
Node(NodePtr deephead,
NodePtr head,
SubtractionDipolePtr dipol,
MatchboxMEBasePtr nodeME,
int cutstage);
/// The destructor.
virtual ~Node();
//@}
public:
// get children from vector<MatchboxMEBasePtr>
void birth(const vector<MatchboxMEBasePtr> & vec);
/// recursive setXComb. proStage is the number of clusterings
/// before the projectors get filled.
void setXComb(tStdXCombPtr xc);
/// calculate the dipole and ps approximation
pair<CrossSection, CrossSection> calcDipandPS(Energy scale)const;
/// calculate the ps approximation
CrossSection calcPs(Energy scale)const;
/// calculate the dipole
CrossSection calcDip(Energy scale)const;
/// recursive flush caches and clean up XCombs.
void flushCaches();
/// recursive clearKinematics
void clearKinematics();
/// recursive setKinematics
void setKinematics();
/// recursive generateKinematics using tilde kinematics of the dipoles
bool generateKinematics(const double *r, bool directCut);
/// generate the kinamatics of the first node
bool firstgenerateKinematics(const double *r, bool directCut);
//return the ME
const MatchboxMEBasePtr nodeME() const;
//return the node ME
MatchboxMEBasePtr nodeME();
//return the parent Node
NodePtr parent() const {return theparent;}
/// vector of children nodes created in birth
vector< NodePtr > children() const {return thechildren;}
//pick a random child (flat)
NodePtr randomChild();
/// true if all children show scales above pt
bool allAbove(Energy pt);
/// true if the node is in the history of other.
bool isInHistoryOf(NodePtr other);
/// legsize of the node ME
int legsize() const;
/// set the first node (first men). only use in factory
void deepHead(NodePtr deephead) {theDeepHead = deephead;}
/// return the first node
NodePtr deepHead() const {return theDeepHead;}
/// returns the dipol of the node.
SubtractionDipolePtr dipole() const;
/// return the xcomb
StdXCombPtr xcomb() const;
/// return the xcomb (if not created, create one from head)
StdXCombPtr xcomb() ;
/// return the current running pt
Energy runningPt() const { return theRunningPt; }
/// set the current running pt
void runningPt(Energy x) { theRunningPt=x; }
/// return the cut stage to cut on merging pt in generate kinematics
int cutStage() const { return theCutStage; }
/// get a vector of the next nodes, ordered in pt (and in parton shower phace space)
vector<NodePtr> getNextOrderedNodes(bool normal=true, double hardscalefactor=1.) const;
//true if the node is in shower history for a given pt
bool inShowerPS(Energy hardpt)const;
//get the history
NodePtr getHistory(bool normal=true, double hardscalefactor=1.);
//true if node correspond to a subtracted real.
bool subtractedReal() const {return theSubtractedReal;}
/// set if node correspont to a subtracted real.
void subtractedReal(bool x) { theSubtractedReal = x;}
//true if node correspond to a virtual contribution.
bool virtualContribution() const { return theVirtualContribution ;}
/// set if node correspont to a virtual contribution.
void virtualContribution(bool x) {theVirtualContribution = x;}
//pointer to the merging helper
MergerPtr MH()const{return theMergingHelper;}
/// set the merging helper
void MH(MergerPtr a){theMergingHelper=a;}
/// pT of the dipole
Energy pT()const{return dipole()->lastPt();}
/// get incoming and outgoing particles (TODO: expensive)
pair<PVector , PVector> getInOut();
private:
/// the Matrixelement representing this node.
MatchboxMEBasePtr thenodeMEPtr;
/// the dipol used to substract
/// and generate kinematics using tilde kinematics
SubtractionDipolePtr thedipol;
/// the parent node
NodePtr theparent;
/// The godfather node of whole tree.(Firstnode)
NodePtr theDeepHead;
/**
* The CutStage is number of clusterings which are possible without
* introducing a merging scale to cut away singularities.
* -> subtracted MEs have the CutStage 1.
* -> virtual and normal tree level ME get 0.
*/
int theCutStage;
/// tell if node belongs to an ordered history
bool isOrdered;
/// flag to tell if node is subtracted real
bool theSubtractedReal;
/// flag to tell if node is virtual contribution
bool theVirtualContribution;
/// the merging helper
MergerPtr theMergingHelper;
//the xcomb of the node
StdXCombPtr thexcomb;
/// vector of the children node
vector< NodePtr > thechildren;
/// the current running pt
Energy theRunningPt;
/// The nodes of the projection stage.
NodePtr theProjector;
+ /// flag not to enter infinite loop. (There should be a better solution...)
+ bool didflush=false;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Node & operator=(const Node &) = delete;
};
}
#endif /* Herwig_Node_H */
diff --git a/src/Merging/3LO-0NLO-LHC-Z.in b/src/Merging/3LO-0NLO-LHC-Z.in
--- a/src/Merging/3LO-0NLO-LHC-Z.in
+++ b/src/Merging/3LO-0NLO-LHC-Z.in
@@ -1,170 +1,170 @@
# -*- ThePEG-repository -*-
##################################################
## Herwig/Merging example input file
##################################################
##################################################
## Collider type
##################################################
read Matchbox/PPCollider.in
read Merging/Merging.in
##################################################
## Beam energy sqrt(s)
##################################################
cd /Herwig/EventHandlers
set EventHandler:LuminosityFunction:Energy 7000*GeV
##################################################
## Process selection
##################################################
## Note that event generation may fail if no matching matrix element has
## been found. Coupling orders are with respect to the Born process,
## i.e. NLO QCD does not require an additional power of alphas.
## Model assumptions
read Matchbox/StandardModelLike.in
read Matchbox/DiagonalCKM.in
## Set the order of the couplings
cd /Herwig/Merging
set MergingFactory:OrderInAlphaS 0
set MergingFactory:OrderInAlphaEW 2
## Select the process
## You may use identifiers such as p, pbar, j, l, mu+, h0 etc.
do MergingFactory:Process p p -> e+ e- [ j j]
set MergingFactory:NLOProcesses 0
set Merger:MergingScale 10.*GeV
set Merger:MergingScaleSmearing 0.1
set Merger:gamma 1.
set MPreWeight:HTPower 0
set MPreWeight:MaxPTPower 0
set MPreWeight:OnlyColoured False
read Matchbox/PQCDLevel.in
set /Herwig/Samplers/Sampler:BinSampler:Kappa 0.02
## Special settings required for on-shell production of unstable particles
## enable for on-shell top production
# read Matchbox/OnShellTopProduction.in
## enable for on-shell W, Z or h production
# read Matchbox/OnShellWProduction.in
# read Matchbox/OnShellZProduction.in
# read Matchbox/OnShellHProduction.in
# Special settings for the VBF approximation
# read Matchbox/VBFDiagramsOnly.in
##################################################
## Matrix element library selection
##################################################
## Select a generic tree/loop combination or a
## specialized NLO package
# read Matchbox/MadGraph-GoSam.in
-# read Matchbox/MadGraph-MadGraph.in
+ read Matchbox/MadGraph-MadGraph.in
# read Matchbox/MadGraph-NJet.in
# read Matchbox/MadGraph-OpenLoops.in
# read Matchbox/HJets.in
# read Matchbox/VBFNLO.in
## Uncomment this to use ggh effective couplings
## currently only supported by MadGraph-GoSam
# read Matchbox/HiggsEffective.in
##################################################
## Cut selection
## See the documentation for more options
##################################################
set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV
set /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV
cd /Herwig/MatrixElements/Matchbox/Utility
insert DiagramGenerator:ExcludeInternal 0 /Herwig/Particles/gamma
## cuts on additional jets
# read Matchbox/DefaultPPJets.in
# insert JetCuts:JetRegions 0 FirstJet
# insert JetCuts:JetRegions 1 SecondJet
# insert JetCuts:JetRegions 2 ThirdJet
# insert JetCuts:JetRegions 3 FourthJet
##################################################
## Scale choice
## See the documentation for more options
##################################################
cd /Herwig/MatrixElements/Matchbox/Scales/
set /Herwig/Merging/MergingFactory:ScaleChoice LeptonPairMassScale
##################################################
## Scale uncertainties
##################################################
# read Matchbox/MuDown.in
# read Matchbox/MuUp.in
##################################################
## Shower scale uncertainties
##################################################
# read Matchbox/MuQDown.in
# read Matchbox/MuQUp.in
##################################################
## PDF choice
##################################################
read Matchbox/FiveFlavourNoBMassScheme.in
read Matchbox/MMHT2014.in
##################################################
## CMW - Scheme
##################################################
### Use factor in alpha_s argument: alpha_s(q) -> alpha_s(fac*q)
### with fac=exp(-(67-3pi^2-10/3*Nf)/(33-2Nf))
read Merging/FactorCMWScheme.in
### Linear CMW multiplication:
### alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi )
# read Merging/LinearCMWScheme.in
##################################################
## Analyses
##################################################
cd /Herwig/Analysis
insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 Rivet
# insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 HepMC
read Merging/LHC7-Z-Analysis.in
##################################################
## Save the generator
##################################################
do /Herwig/Merging/MergingFactory:ProductionMode
set /Herwig/Generators/EventGenerator:IntermediateOutput Yes
cd /Herwig/Generators
saverun 3LO-0NLO-LHC-Z EventGenerator
diff --git a/src/Merging/4LO-0NLO-LHC-Z.in b/src/Merging/4LO-0NLO-LHC-Z.in
--- a/src/Merging/4LO-0NLO-LHC-Z.in
+++ b/src/Merging/4LO-0NLO-LHC-Z.in
@@ -1,170 +1,170 @@
# -*- ThePEG-repository -*-
##################################################
## Herwig/Merging example input file
##################################################
##################################################
## Collider type
##################################################
read Matchbox/PPCollider.in
read Merging/Merging.in
##################################################
## Beam energy sqrt(s)
##################################################
cd /Herwig/EventHandlers
set EventHandler:LuminosityFunction:Energy 7000*GeV
##################################################
## Process selection
##################################################
## Note that event generation may fail if no matching matrix element has
## been found. Coupling orders are with respect to the Born process,
## i.e. NLO QCD does not require an additional power of alphas.
## Model assumptions
read Matchbox/StandardModelLike.in
read Matchbox/DiagonalCKM.in
## Set the order of the couplings
cd /Herwig/Merging
set MergingFactory:OrderInAlphaS 0
set MergingFactory:OrderInAlphaEW 2
## Select the process
## You may use identifiers such as p, pbar, j, l, mu+, h0 etc.
do MergingFactory:Process p p -> e+ e- [ j j j]
set MergingFactory:NLOProcesses 0
set Merger:MergingScale 10.*GeV
set Merger:MergingScaleSmearing 0.1
set Merger:gamma 1.
set MPreWeight:HTPower 0
set MPreWeight:MaxPTPower 0
set MPreWeight:OnlyColoured False
read Matchbox/PQCDLevel.in
## Special settings required for on-shell production of unstable particles
## enable for on-shell top production
# read Matchbox/OnShellTopProduction.in
## enable for on-shell W, Z or h production
# read Matchbox/OnShellWProduction.in
# read Matchbox/OnShellZProduction.in
# read Matchbox/OnShellHProduction.in
# Special settings for the VBF approximation
# read Matchbox/VBFDiagramsOnly.in
##################################################
## Matrix element library selection
##################################################
## Select a generic tree/loop combination or a
## specialized NLO package
# read Matchbox/MadGraph-GoSam.in
-# read Matchbox/MadGraph-MadGraph.in
+ read Matchbox/MadGraph-MadGraph.in
# read Matchbox/MadGraph-NJet.in
-read Matchbox/MadGraph-OpenLoops.in
+# read Matchbox/MadGraph-OpenLoops.in
# read Matchbox/HJets.in
# read Matchbox/VBFNLO.in
## Uncomment this to use ggh effective couplings
## currently only supported by MadGraph-GoSam
# read Matchbox/HiggsEffective.in
##################################################
## Cut selection
## See the documentation for more options
##################################################
set /Herwig/Cuts/ChargedLeptonPairMassCut:MinMass 60*GeV
set /Herwig/Cuts/ChargedLeptonPairMassCut:MaxMass 120*GeV
cd /Herwig/MatrixElements/Matchbox/Utility
insert DiagramGenerator:ExcludeInternal 0 /Herwig/Particles/gamma
## cuts on additional jets
# read Matchbox/DefaultPPJets.in
# insert JetCuts:JetRegions 0 FirstJet
# insert JetCuts:JetRegions 1 SecondJet
# insert JetCuts:JetRegions 2 ThirdJet
# insert JetCuts:JetRegions 3 FourthJet
##################################################
## Scale choice
## See the documentation for more options
##################################################
cd /Herwig/MatrixElements/Matchbox/Scales/
set /Herwig/Merging/MergingFactory:ScaleChoice LeptonPairMassScale
##################################################
## Scale uncertainties
##################################################
# read Matchbox/MuDown.in
# read Matchbox/MuUp.in
##################################################
## Shower scale uncertainties
##################################################
# read Matchbox/MuQDown.in
# read Matchbox/MuQUp.in
##################################################
## PDF choice
##################################################
read Matchbox/FiveFlavourNoBMassScheme.in
read Matchbox/MMHT2014.in
##################################################
## CMW - Scheme
##################################################
### Use factor in alpha_s argument: alpha_s(q) -> alpha_s(fac*q)
### with fac=exp(-(67-3pi^2-10/3*Nf)/(33-2Nf))
read Merging/FactorCMWScheme.in
### Linear CMW multiplication:
### alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi )
# read Merging/LinearCMWScheme.in
##################################################
## Analyses
##################################################
cd /Herwig/Analysis
insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 Rivet
# insert /Herwig/Generators/EventGenerator:AnalysisHandlers 0 HepMC
read Merging/LHC7-Z-Analysis.in
##################################################
## Save the generator
##################################################
do /Herwig/Merging/MergingFactory:ProductionMode
set /Herwig/Generators/EventGenerator:IntermediateOutput Yes
cd /Herwig/Generators
saverun 4LO-0NLO-LHC-Z EventGenerator

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:55 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3799605
Default Alt Text
(144 KB)

Event Timeline