Page MenuHomeHEPForge
Files F7879093

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Powheg/ME2byDipoles.cc b/MatrixElement/Matchbox/Powheg/ME2byDipoles.cc
--- a/MatrixElement/Matchbox/Powheg/ME2byDipoles.cc
+++ b/MatrixElement/Matchbox/Powheg/ME2byDipoles.cc
@@ -1,342 +1,342 @@
// -*- C++ -*-
//
// ME2byDipoles.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 ME2byDipoles class.
//
#include "ME2byDipoles.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
ME2byDipoles::ME2byDipoles()
: MatchboxReweightBase() {}
ME2byDipoles::~ME2byDipoles() {}
IBPtr ME2byDipoles::clone() const {
return new_ptr(*this);
}
IBPtr ME2byDipoles::fullclone() const {
return new_ptr(*this);
}
void ME2byDipoles::setup(Ptr<SubtractedME>::tptr sub) {
theRealME = dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(sub->head());
assert(theRealME);
Ptr<MatchboxMEBase>::ptr nreal = theRealME->cloneMe();
ostringstream pname;
pname << fullName() << "/" << nreal->name();
if ( ! (generator()->preinitRegister(nreal,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
nreal->cloneDependencies();
theRealME = nreal;
vector<Ptr<SubtractionDipole>::ptr> inDipoles = sub->dipoles();
theDipoles.clear();
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
inDipoles.begin(); d != inDipoles.end(); ++d ) {
Ptr<SubtractionDipole>::ptr ndipole = (**d).cloneMe();
ostringstream dname;
dname << fullName() << "/" << (**d).name();
if ( ! (generator()->preinitRegister(ndipole,dname.str())) )
throw InitException() << "Dipole '" << dname.str() << "' already existing.";
ndipole->cloneDependencies();
theDipoles.push_back(ndipole);
}
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
(**d).doSubtraction();
}
void ME2byDipoles::setup(Ptr<SubtractionDipole>::tptr dip,
Ptr<SubtractedME>::tptr sub) {
vector<Ptr<SubtractionDipole>::ptr> inDipoles = sub->dipoles();
theDipoles.clear();
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
inDipoles.begin(); d != inDipoles.end(); ++d ) {
Ptr<SubtractionDipole>::ptr ndipole = (**d).cloneMe();
ostringstream dname;
dname << fullName() << "/" << (**d).name();
if ( ! (generator()->preinitRegister(ndipole,dname.str())) )
throw InitException() << "Dipole '" << dname.str() << "' already existing.";
ndipole->cloneDependencies();
theDipoles.push_back(ndipole);
if ( *d == dip ) {
projectionDipole(theDipoles.back());
}
}
theRealME = Ptr<MatchboxMEBase>::ptr();
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
(**d).doSubtraction();
}
void ME2byDipoles::setXComb(tStdXCombPtr real) {
assert(real);
theLastXComb = real;
map<StdXCombPtr,vector<StdDependentXCombPtr> >::iterator xcs
= theXCombMap.find(theLastXComb);
if ( xcs == theXCombMap.end() ) {
getXCombs(theLastXComb);
xcs = theXCombMap.find(theLastXComb);
}
if ( theRealME )
theRealME->setXComb(theLastXComb);
vector<StdDependentXCombPtr>::iterator xcit = xcs->second.begin();
vector<Ptr<SubtractionDipole>::ptr>::iterator dip = theDipoles.begin();
- for ( ; xcit != xcs->second.end(); ++xcit, ++dip )
+ for ( ; xcit != xcs->second.end(); ++xcit, ++dip ) {
+ (**xcit).clean();
(**dip).setXComb(*xcit);
+ }
}
double ME2byDipoles::scaledBornScreen() const {
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' evaluating Born screening\n";
Energy scale = projectionDipole()->lastDipoleScale();
Energy pt = projectionDipole()->lastPt();
if ( projectionDipole()->verbose() )
generator()->log() << "from pt/GeV = " << (pt/GeV)
<< " scale/GeV = " << (scale/GeV)
<< "\n" << flush;
return pow(pt/scale,4.);
}
double ME2byDipoles::scaledBorn(Energy2 factorizationScale) const {
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' evaluating scaled Born\n" << flush;
projectionDipole()->underlyingBornME()->setScale();
projectionDipole()->underlyingBornME()->getPDFWeight(factorizationScale);
double me2 = projectionDipole()->underlyingBornME()->me2();
double pdf = projectionDipole()->underlyingBornME()->lastXComb().lastMEPDFWeight();
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' done evaluating scaled Born\n" << flush;
return me2 * pdf;
}
void ME2byDipoles::flushCaches() {
if ( theRealME )
theRealME->flushCaches();
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
(**d).flushCaches();
}
double ME2byDipoles::evaluate(double& sratio) const {
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' ME2byDipoles evaluating\n" << flush;
double den = 0.0;
sratio = 0.;
double numDip = 0.0;
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator dip =
theDipoles.begin(); dip != theDipoles.end(); ++dip ) {
if ( !(**dip).apply() )
continue;
tStdDependentXCombPtr depXComb =
dynamic_ptr_cast<tStdDependentXCombPtr>((**dip).lastXCombPtr());
assert(depXComb);
- depXComb->setProcess();
if ( !(**dip).generateTildeKinematics() )
continue;
- depXComb->remakeIncoming();
depXComb->setIncomingPartons();
(**dip).realEmissionME()->setScale();
(**dip).underlyingBornME()->setScale();
double res = (**dip).me2();
den += abs(res);
if ( depXComb->willPassCuts() )
sratio += res;
if ( *dip == projectionDipole() ) {
numDip = abs(res);
}
}
if ( sratio != 0. ) {
assert(abs(den) != 0.0);
sratio /= den;
}
if ( theRealME ) {
if ( !realME()->lastXCombPtr()->willPassCuts() ) {
if ( projectionDipole()->verbose() )
generator()->log() << "real emission did not pass the cuts\n" << flush;
return 0.0;
}
}
assert(abs(den) != 0.);
double num = theRealME ? realME()->me2() : numDip;
double res = num / den;
if ( projectionDipole()->verbose() ) {
generator()->log() << "'" << name() << "' done evaluating\n"
<< "numerator = " << num << " denominator = "
<< den << "\n" << flush;
}
return res;
}
void ME2byDipoles::getXCombs(tStdXCombPtr xc) {
vector<StdDependentXCombPtr> xcs;
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator dip = theDipoles.begin();
dip != theDipoles.end(); ++dip ) {
StdDependentXCombPtr depxc = (**dip).makeBornXComb(xc);
xcs.push_back(depxc);
}
theXCombMap[xc] = xcs;
}
void ME2byDipoles::print(ostream& os) const {
os << "--- ME2byDipoles setup ---------------------------------------------------------\n";
os << " '" << name() << "'\n"
<< " real emission matrix element '" << theRealME->name() << "'\n"
<< " projection dipole: '"
<< (projectionDipole() ? projectionDipole()->name() : "")
<< "'\n";
os << " associated dipoles are:\n";
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator dip =
theDipoles.begin(); dip != theDipoles.end(); ++dip ) {
os << " '" << (**dip).name() << "'\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void ME2byDipoles::printLastEvent(ostream& os) const {
os << "--- ME2byDipoles last event information ----------------------------------------\n";
os << " for ratio '" << name() << "'\n";
os << " real emission event information:\n";
if ( dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(theRealME) )
dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(theRealME)->printLastEvent(os);
else if ( dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(theRealME) )
dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(theRealME)->printLastEvent(os);
else
os << " unknown MEBase object.\n";
if ( projectionDipole() ) {
os << " projection dipole event information:\n";
projectionDipole()->printLastEvent(os);
}
os << " dipoles event information:\n";
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d = theDipoles.begin();
d != theDipoles.end(); ++d )
(**d).printLastEvent(os);
os << "--- end ME2byDipoles last event information ------------------------------------\n";
os << flush;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void ME2byDipoles::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theRealME << theProjectionDipole << theDipoles;
}
void ME2byDipoles::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theRealME >> theProjectionDipole >> theDipoles;
}
void ME2byDipoles::Init() {
static ClassDocumentation<ME2byDipoles> documentation
("ME2byDipoles");
static Reference<ME2byDipoles,MatchboxMEBase> interfaceRealME
("RealME",
"The real emission matrix element.",
&ME2byDipoles::theRealME, false, false, true, false, false);
static Reference<ME2byDipoles,SubtractionDipole> interfaceProjectionDipole
("ProjectionDipole",
"The projection dipole.",
&ME2byDipoles::theProjectionDipole, false, false, true, false, false);
static RefVector<ME2byDipoles,SubtractionDipole> interfaceDipoles
("Dipoles",
"The dipoles associated to the real emission matrix element.",
&ME2byDipoles::theDipoles, -1, false, false, true, false, false);
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<ME2byDipoles,MatchboxReweightBase>
describeME2byDipoles("Herwig::ME2byDipoles", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Powheg/PowhegSplittingGenerator.cc b/MatrixElement/Matchbox/Powheg/PowhegSplittingGenerator.cc
--- a/MatrixElement/Matchbox/Powheg/PowhegSplittingGenerator.cc
+++ b/MatrixElement/Matchbox/Powheg/PowhegSplittingGenerator.cc
@@ -1,435 +1,436 @@
// -*- C++ -*-
//
// PowhegSplittingGenerator.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PowhegSplittingGenerator class.
//
#include "PowhegSplittingGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "PowhegInclusiveME.h"
#include "ThePEG/PDF/PDF.h"
using namespace Herwig;
PowhegSplittingGenerator::PowhegSplittingGenerator()
: StepHandler(),
theFFPtCut(1.0*GeV), theFFScreeningScale(ZERO),
theFIPtCut(1.0*GeV), theFIScreeningScale(ZERO),
theIIPtCut(1.0*GeV), theIIScreeningScale(ZERO),
discardNoEmissions(false), theVerbose(false),
theDiscardNext(false) {}
PowhegSplittingGenerator::~PowhegSplittingGenerator() {
for ( GeneratorMap::iterator g = theGeneratorMap.begin();
g != theGeneratorMap.end(); ++g )
delete g->second.second;
theGeneratorMap.clear();
}
void PowhegSplittingGenerator::
handle(EventHandler & eh, const tPVector &, const Hint &) {
if ( theVerbose ) {
generator()->log() << "PowhegSplittingGenerator generating real emission off the sub-process\n"
<< (*(eh.lastXCombPtr()->subProcess())) << "\n"
<< "with x1 = " << eh.lastXCombPtr()->lastX1()
<< " x2 = " << eh.lastXCombPtr()->lastX2() << "\n" << flush;
}
if ( !generate(eh) ) {
if ( theVerbose ) {
generator()->log() << "PowhegSplittingGenerator did not select radiation above the IR cutoff\n" << flush;
}
if ( discardNext() ) {
setDiscardNext(false);
if ( theVerbose ) {
generator()->log() << "Splitting kernels have been presampled, will discard this event.\n" << flush;
}
throw Veto();
}
if ( discardNoEmissions )
throw Veto();
veto(eh);
return;
}
if ( theVerbose ) {
generator()->log() << "PowhegSplittingGenerator selected the kernel '"
<< lastSplitting()->name() << "' to generate radiation\n" << flush;
}
if ( discardNext() ) {
setDiscardNext(false);
if ( theVerbose ) {
generator()->log() << "Splitting kernels have been presampled, will discard this event.\n" << flush;
}
throw Veto();
}
SubProPtr oldSub =
lastSplitting()->bornSubProcess();
SubProPtr newSub;
try {
Energy pt = lastSplitting()->projectionDipole()->lastPt();
newSub = lastSplitting()->construct(pt);
} catch(Veto&) {
if ( theVerbose ) {
generator()->log() << "The generated real emission process did not pass the cuts.\n" << flush;
}
veto(eh);
return;
}
tParticleSet firstS = oldSub->incoming().first->siblings();
assert(firstS.empty() || firstS.size() == 1);
if ( !firstS.empty() ) {
eh.currentStep()->removeParticle(*firstS.begin());
}
tParticleSet secondS = oldSub->incoming().second->siblings();
assert(secondS.empty() || secondS.size() == 1);
if ( !secondS.empty() ) {
eh.currentStep()->removeParticle(*secondS.begin());
}
// prevent the colliding particles from disappearing
// in the initial state and appearing
// in the final state when we've cut off all their
// (physical) children; only applies to the case
// where we have a parton extractor not build from
// noPDF, so check wether the incoming particle
// doesnt equal the incoming parton -- this needs fixing in ThePEG
PPtr dummy = new_ptr(Particle(getParticleData(ParticleID::gamma)));
bool usedDummy = false;
if ( eh.currentStep()->incoming().first != oldSub->incoming().first ) {
eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().first,dummy);
usedDummy = true;
}
if ( eh.currentStep()->incoming().second != oldSub->incoming().second ) {
eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().second,dummy);
usedDummy = true;
}
eh.currentStep()->removeSubProcess(oldSub);
eh.currentStep()->addSubProcess(newSub);
// get rid of the dummy
if ( usedDummy ) {
eh.currentStep()->removeParticle(dummy);
}
eh.select(lastSplitting()->lastXCombPtr());
- dynamic_ptr_cast<tStdDependentXCombPtr>(lastSplitting()->lastXCombPtr())->setPartonBinInstances();
+ dynamic_ptr_cast<tStdDependentXCombPtr>(lastSplitting()->lastXCombPtr())->
+ recreatePartonBinInstances(lastSplitting()->lastXCombPtr()->lastScale());
eh.lastExtractor()->constructRemnants(lastSplitting()->lastXCombPtr()->partonBinInstances(),
newSub, eh.currentStep());
if ( theVerbose ) {
generator()->log() << "PowhegSplittingGenerator generated the real emission sub-process\n"
<< (*(eh.lastXCombPtr()->subProcess())) << "\n"
<< "with x1 = " << eh.lastXCombPtr()->lastX1()
<< " x2 = " << eh.lastXCombPtr()->lastX2() << "\n" << flush;
}
}
bool PowhegSplittingGenerator::generate(EventHandler & eh) {
pair<GeneratorMap::iterator,GeneratorMap::iterator> generators
= getGenerators(eh);
Energy winnerPt = 0.*GeV;
Energy pt;
GeneratorMap::iterator winner = theGeneratorMap.end();
for ( GeneratorMap::iterator gen = generators.first; gen != generators.second; ++gen ) {
assert(gen->second.first->lastHeadXCombPtr() == eh.lastXCombPtr());
assert(gen->second.first->projectionDipole()->lastHeadXCombPtr() == eh.lastXCombPtr());
assert(gen->second.first->projectionDipole()->lastXCombPtr() ==
gen->second.first->lastXCombPtr());
pt = generate(gen->second);
if ( pt > winnerPt ) {
winnerPt = pt;
winner = gen;
}
}
if ( winner == theGeneratorMap.end() ) {
theLastSplitting = Ptr<PowhegSplittingKernel>::tptr();
return false;
}
theLastSplitting = winner->second.first;
return true;
}
void PowhegSplittingGenerator::veto(EventHandler & eh) const {
tSubProPtr sub = eh.currentStep()->subProcesses().front();
if ( sub->incoming().first->coloured() ) {
sub->incoming().first->vetoScale(ZERO);
}
if ( sub->incoming().second->coloured() ) {
sub->incoming().first->vetoScale(ZERO);
}
for ( ParticleVector::const_iterator p = sub->outgoing().begin();
p != sub->outgoing().end(); ++p )
if ( (**p).coloured() )
(**p).vetoScale(ZERO);
}
IBPtr PowhegSplittingGenerator::clone() const {
return new_ptr(*this);
}
IBPtr PowhegSplittingGenerator::fullclone() const {
return new_ptr(*this);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void PowhegSplittingGenerator::persistentOutput(PersistentOStream & os) const {
os << ounit(theFFPtCut,GeV) << ounit(theFFScreeningScale,GeV)
<< ounit(theFIPtCut,GeV) << ounit(theFIScreeningScale,GeV)
<< ounit(theIIPtCut,GeV) << ounit(theIIScreeningScale,GeV)
<< discardNoEmissions << theVerbose;
}
void PowhegSplittingGenerator::persistentInput(PersistentIStream & is, int) {
is >> iunit(theFFPtCut,GeV) >> iunit(theFFScreeningScale,GeV)
>> iunit(theFIPtCut,GeV) >> iunit(theFIScreeningScale,GeV)
>> iunit(theIIPtCut,GeV) >> iunit(theIIScreeningScale,GeV)
>> discardNoEmissions >> theVerbose;
}
pair<PowhegSplittingGenerator::GeneratorMap::iterator,
PowhegSplittingGenerator::GeneratorMap::iterator>
PowhegSplittingGenerator::getGenerators(EventHandler& eh) {
tXCombPtr xc = eh.lastXCombPtr();
const ThePEG::StandardXComb& xcref = *dynamic_ptr_cast<tStdXCombPtr>(xc);
if ( theVerbose ) {
generator()->log() << "getting splitting generators for xcomb "
<< xc << " and process ";
generator()->log() << xcref.mePartonData()[0]->PDGName() << " "
<< xcref.mePartonData()[1]->PDGName() << " -> ";
for ( ThePEG::cPDVector::const_iterator pid =
xcref.mePartonData().begin() + 2;
pid != xcref.mePartonData().end(); ++pid )
generator()->log() << (**pid).PDGName() << " ";
generator()->log() << "\n" << flush;
}
pair<GeneratorMap::iterator,GeneratorMap::iterator> res =
theGeneratorMap.equal_range(xc);
if ( res.first != res.second ) {
if ( theVerbose )
generator()->log() << "generators already known\n" << flush;
return res;
}
tStdXCombGroupPtr xcGrp =
dynamic_ptr_cast<tStdXCombGroupPtr>(xc);
if ( !xcGrp ) {
if ( theVerbose )
generator()->log() << "xcomb is not an xcomb group\n" << flush;
return make_pair(theGeneratorMap.end(),theGeneratorMap.end());
}
Ptr<PowhegInclusiveME>::tptr me =
dynamic_ptr_cast<Ptr<PowhegInclusiveME>::tptr>((*xcGrp).matrixElement());
if ( !me ) {
if ( theVerbose )
generator()->log() << "matrix element is not a powheg inclusive me\n" << flush;
return make_pair(theGeneratorMap.end(),theGeneratorMap.end());
}
for ( vector<Ptr<PowhegSplittingKernel>::ptr>::iterator k =
me->splittingKernels().begin(); k != me->splittingKernels().end(); ++k ) {
if ( !(**k).apply() )
continue;
if ( theVerbose )
generator()->log() << "initializing generator for kernel '" << (**k).name() << "'\n" << flush;
assert((**k).lastHeadXCombPtr() == eh.lastXCombPtr());
assert((**k).projectionDipole()->lastHeadXCombPtr() == eh.lastXCombPtr());
assert((**k).lastXCombPtr() == (**k).projectionDipole()->lastXCombPtr());
(**k).splittingGenerator(this);
if ( (**k).projectionDipole()->realEmitter() > 1 &&
(**k).projectionDipole()->realSpectator() > 1 ) {
(**k).ptCut(theFFPtCut);
(**k).screeningScale(theFFScreeningScale);
} else if ( (**k).projectionDipole()->realEmitter() < 2 &&
(**k).projectionDipole()->realSpectator() < 2 ) {
(**k).ptCut(theIIPtCut);
(**k).screeningScale(theIIScreeningScale);
} else {
(**k).ptCut(theFIPtCut);
(**k).screeningScale(theFIScreeningScale);
}
ExponentialGeneratorPtr gen = new ExponentialGenerator();
gen->sampling_parameters().maxtry = (**k).maxtry();
gen->sampling_parameters().presampling_points = (**k).presamplingPoints();
gen->function(*k);
gen->initialize();
theGeneratorMap.insert(make_pair(xc,make_pair(*k,gen)));
}
return getGenerators(eh);
}
Energy PowhegSplittingGenerator::generate(pair<Ptr<PowhegSplittingKernel>::ptr,ExponentialGeneratorPtr>& gen) {
double res = 0.;
gen.first->splittingGenerator(this);
while (true) {
try {
res = gen.second->generate();
} catch (exsample::exponential_regenerate&) {
continue;
} catch (exsample::hit_and_miss_maxtry&) {
throw Veto();
} catch (exsample::selection_maxtry&) {
throw Veto();
}
break;
}
if ( theVerbose ) {
generator()->log() << "Generating splitting from '" << gen.first->projectionDipole()->name() << "'.\n" << flush;
if ( res == 0. )
generator()->log() << "Below infrared cutoff.\n" << flush;
else
generator()->log() << "pt/GeV = " << (gen.first->projectionDipole()->lastPt()/GeV) << ".\n" << flush;
}
if ( res == 0. )
return 0.*GeV;
return gen.first->projectionDipole()->lastPt();
}
void PowhegSplittingGenerator::Init() {
static ClassDocumentation<PowhegSplittingGenerator> documentation
("PowhegSplittingGenerator");
static Parameter<PowhegSplittingGenerator,Energy> interfaceFFPtCut
("FFPtCut",
"Set the pt infrared cutoff",
&PowhegSplittingGenerator::theFFPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<PowhegSplittingGenerator,Energy> interfaceFFScreeningScale
("FFScreeningScale",
"Set the screening scale",
&PowhegSplittingGenerator::theFFScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<PowhegSplittingGenerator,Energy> interfaceFIPtCut
("FIPtCut",
"Set the pt infrared cutoff",
&PowhegSplittingGenerator::theFIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<PowhegSplittingGenerator,Energy> interfaceFIScreeningScale
("FIScreeningScale",
"Set the screening scale",
&PowhegSplittingGenerator::theFIScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<PowhegSplittingGenerator,Energy> interfaceIIPtCut
("IIPtCut",
"Set the pt infrared cutoff",
&PowhegSplittingGenerator::theIIPtCut, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<PowhegSplittingGenerator,Energy> interfaceIIScreeningScale
("IIScreeningScale",
"Set the screening scale",
&PowhegSplittingGenerator::theIIScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<PowhegSplittingGenerator,bool> interfaceVerbose
("Verbose",
"",
&PowhegSplittingGenerator::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"",
false);
static Switch<PowhegSplittingGenerator,bool> interfaceDiscardNoEmissions
("DiscardNoEmissions",
"Discard events without radiation.",
&PowhegSplittingGenerator::discardNoEmissions, false, false, false);
static SwitchOption interfaceDiscardNoEmissionsOn
(interfaceDiscardNoEmissions,
"On",
"Discard events without radiation.",
true);
static SwitchOption interfaceDiscardNoEmissionsOff
(interfaceDiscardNoEmissions,
"Off",
"Do not discard events without radiation.",
false);
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<PowhegSplittingGenerator,StepHandler>
describeHerwigPowhegSplittingGenerator("Herwig::PowhegSplittingGenerator", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Powheg/PowhegSplittingKernel.cc b/MatrixElement/Matchbox/Powheg/PowhegSplittingKernel.cc
--- a/MatrixElement/Matchbox/Powheg/PowhegSplittingKernel.cc
+++ b/MatrixElement/Matchbox/Powheg/PowhegSplittingKernel.cc
@@ -1,371 +1,371 @@
// -*- C++ -*-
//
// PowhegSplittingKernel.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 PowhegSplittingKernel class.
//
#include "PowhegSplittingKernel.h"
#include "PowhegSplittingGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Cuts/Cuts.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxNLOME.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
PowhegSplittingKernel::PowhegSplittingKernel()
: ME2byDipoles(), theBornScreening(true),
thePresamplingPoints(10000), theMaxTry(100000),
theScreeningScale(ZERO),
theBornRandom(0,0), theRadiationRandom(0,0),
thePresampling(false) {}
PowhegSplittingKernel::~PowhegSplittingKernel() {}
IBPtr PowhegSplittingKernel::clone() const {
return new_ptr(*this);
}
IBPtr PowhegSplittingKernel::fullclone() const {
return new_ptr(*this);
}
void PowhegSplittingKernel::splittingGenerator(Ptr<PowhegSplittingGenerator>::tptr gen) {
theGenerator = gen;
}
void PowhegSplittingKernel::setXComb(tStdXCombPtr real) {
if ( real )
ME2byDipoles::setXComb(real);
projectionDipole()->setXComb(real);
if ( !real )
return;
if ( theBornRandom.first == theBornRandom.second ) {
theBornRandom.first =
projectionDipole()->lastHeadXComb().pExtractor()->
nDims(projectionDipole()->lastHeadXComb().partonBins()).first;
theBornRandom.second =
theBornRandom.first + projectionDipole()->underlyingBornME()->nDim();
theRadiationRandom.first = theBornRandom.second +
( projectionDipole()->lastHeadXComb().matrixElement()->nDim() - // includes radiation & the coll remainders
projectionDipole()->underlyingBornME()->nDim() -
projectionDipole()->nDimRadiation() );
theRadiationRandom.second =
theRadiationRandom.first + projectionDipole()->nDimRadiation();
}
if ( thePresamplingXCombs.find(real->head()) == thePresamplingXCombs.end() ) {
tStdXCombPtr bornxc = real->head();
thePresamplingXCombs[bornxc] =
new_ptr(StandardXComb(bornxc->maxEnergy(),bornxc->particles(),
bornxc->eventHandlerPtr(),
const_ptr_cast<tSubHdlPtr>(bornxc->subProcessHandler()),
bornxc->pExtractor(),bornxc->CKKWHandler(),
bornxc->partonBins(),bornxc->cuts(),
projectionDipole()->underlyingBornME(),
bornxc->diagrams(),bornxc->mirror()));
}
}
void PowhegSplittingKernel::flushCaches() {
if ( realME() )
realME()->flushCaches();
if ( projectionDipole() )
projectionDipole()->flushCaches();
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
dipoles().begin(); d != dipoles().end(); ++d )
(**d).flushCaches();
}
tSubProPtr PowhegSplittingKernel::construct(Energy pt) {
lastXCombPtr()->matrixElement()->flushCaches();
tSubProPtr subpro = lastXCombPtr()->construct();
if ( projectionDipole()->realEmitter() == 0 ||
projectionDipole()->realSpectator() == 0 ) {
if ( subpro->incoming().first->vetoScale() < 0.0*GeV2 ||
subpro->incoming().first->vetoScale() > sqr(pt) )
subpro->incoming().first->vetoScale(sqr(pt));
}
if ( projectionDipole()->realEmitter() == 1 ||
projectionDipole()->realSpectator() == 1 ) {
if ( subpro->incoming().second->vetoScale() < 0.0*GeV2 ||
subpro->incoming().second->vetoScale() > sqr(pt) )
subpro->incoming().second->vetoScale(sqr(pt));
}
if ( projectionDipole()->realEmitter() > 1 ) {
if ( subpro->outgoing()[projectionDipole()->realEmitter()-2]->vetoScale() < 0.0*GeV2 ||
subpro->outgoing()[projectionDipole()->realEmitter()-2]->vetoScale() > sqr(pt) )
subpro->outgoing()[projectionDipole()->realEmitter()-2]->vetoScale(sqr(pt));
}
if ( projectionDipole()->realSpectator() > 1 ) {
if ( subpro->outgoing()[projectionDipole()->realSpectator()-2]->vetoScale() < 0.0*GeV2 ||
subpro->outgoing()[projectionDipole()->realSpectator()-2]->vetoScale() > sqr(pt) )
subpro->outgoing()[projectionDipole()->realSpectator()-2]->vetoScale(sqr(pt));
}
if ( subpro->outgoing()[projectionDipole()->realEmission()-2]->vetoScale() < 0.0*GeV2 ||
subpro->outgoing()[projectionDipole()->realEmission()-2]->vetoScale() > sqr(pt) )
subpro->outgoing()[projectionDipole()->realEmission()-2]->vetoScale(sqr(pt));
return subpro;
}
double PowhegSplittingKernel::evaluate() const {
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' evaluating\n";
if ( !projectionDipole()->underlyingBornME()->
lastXCombPtr()->willPassCuts() ) {
if ( projectionDipole()->verbose() )
generator()->log() << "Born did not pass the cuts\n";
return 0.;
}
double fscaleFactor =
projectionDipole()->realEmissionME()->factorizationScaleFactor();
Energy2 fscale =
fscaleFactor * (sqr(projectionDipole()->lastPt()) + sqr(theScreeningScale));
double dummy;
double ratio = ME2byDipoles::evaluate(dummy);
assert(ratio >= 0.);
Energy2 bornSHat =
projectionDipole()->underlyingBornME()->lastXComb().lastSHat();
double bornJacobian =
projectionDipole()->underlyingBornME()->lastXComb().jacobian();
if ( projectionDipole()->verbose() )
generator()->log() << "Born sHat/GeV2 = " << (bornSHat/GeV2)
<< " Born Jacobian = " << bornJacobian << "\n" << flush;
CrossSection born =
sqr(hbarc) * scaledBorn(fscale) * bornJacobian / (2.*bornSHat);
if ( born == ZERO ) {
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' done evaluating\n";
return 0.0;
}
if ( bornScreening() ) {
born += sqr(hbarc) * scaledBornScreen() * bornJacobian / (2.*bornSHat);
}
CrossSection dip = projectionDipole()->dSigHatDR(fscale);
ratio *= ( dip / born );
double rscaleFactor =
projectionDipole()->realEmissionME()->renormalizationScaleFactor();
double runAlpha =
SM().alphaS(rscaleFactor * (sqr(projectionDipole()->lastPt())
+ sqr(theScreeningScale)));
if ( projectionDipole()->verbose() )
generator()->log() << "real emission alpha_s = "
<< projectionDipole()->realEmissionME()->lastXComb().lastAlphaS()
<< " pt running alpha_s = " << runAlpha
<< " from pt/GeV = " << (projectionDipole()->lastPt()/GeV)
<< "\n" << flush;
ratio *= runAlpha / projectionDipole()->realEmissionME()->lastXComb().lastAlphaS();
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' done evaluating\n";
return ratio;
}
int PowhegSplittingKernel::nDim() const {
return
projectionDipole()->lastHeadXComb().pExtractor()->
nDims(projectionDipole()->lastHeadXComb().partonBins()).first +
projectionDipole()->lastHeadXComb().pExtractor()->
nDims(projectionDipole()->lastHeadXComb().partonBins()).second +
projectionDipole()->underlyingBornME()->nDim() +
projectionDipole()->nDimRadiation();
}
const vector<double>& PowhegSplittingKernel::parameterPoint() {
assert(!presampling());
theLastParameterPoint.resize(nDim());
copy(projectionDipole()->lastHeadXComb().lastRandomNumbers().begin(),
projectionDipole()->lastHeadXComb().lastRandomNumbers().begin() + theBornRandom.second,
theLastParameterPoint.begin());
copy(projectionDipole()->lastHeadXComb().lastRandomNumbers().begin() + theRadiationRandom.first,
projectionDipole()->lastHeadXComb().lastRandomNumbers().end(),
theLastParameterPoint.begin() + theBornRandom.second);
theLastParameterPoint[evolutionVariable()] = 1.;
return theLastParameterPoint;
}
int PowhegSplittingKernel::evolutionVariable() const {
return
theBornRandom.second +
projectionDipole()->invertedTildeKinematics()->evolutionVariable();
}
const vector<bool>& PowhegSplittingKernel::sampleFlags() {
if ( !theFlags.empty() )
return theFlags;
theFlags.resize(nDim(),false);
for ( int k = theBornRandom.second;
k < theBornRandom.second + projectionDipole()->nDimRadiation(); ++k )
theFlags[k] = true;
return theFlags;
}
const pair<vector<double>,vector<double> >& PowhegSplittingKernel::support() {
if ( !theSupport.first.empty() )
return theSupport;
vector<double> l(nDim(),0.0);
vector<double> u(nDim(),1.0);
theSupport.first = l;
theSupport.second = u;
return theSupport;
}
void PowhegSplittingKernel::startPresampling() {
thePresampling = true;
theXCombBackup = lastHeadXCombPtr();
lastXCombPtr()->head(thePresamplingXCombs[theXCombBackup]);
ME2byDipoles::setXComb(lastXCombPtr());
projectionDipole()->setXComb(lastXCombPtr());
- thePresamplingXCombs[theXCombBackup]->prepare(theXCombBackup->lastParticles());
}
void PowhegSplittingKernel::stopPresampling() {
thePresampling = false;
lastXCombPtr()->head(theXCombBackup);
projectionDipole()->setXComb(lastXCombPtr());
ME2byDipoles::setXComb(lastXCombPtr());
theGenerator->setDiscardNext();
}
double PowhegSplittingKernel::evaluate(const vector<double>& p) {
flushCaches();
try {
if ( projectionDipole()->verbose() )
generator()->log() << "'" << name() << "' preparing\n";
if ( presampling() ) {
if ( projectionDipole()->verbose() )
generator()->log() << "presampling\n";
int bdim = thePresamplingXCombs[theXCombBackup]->nDim();
thePresamplingPoint.resize(bdim);
copy(p.begin(),p.begin()+theBornRandom.second,
thePresamplingPoint.begin());
copy(p.begin()+theBornRandom.second+projectionDipole()->nDimRadiation(),p.end(),
thePresamplingPoint.begin()+theBornRandom.second);
+ thePresamplingXCombs[theXCombBackup]->prepare(theXCombBackup->lastParticles());
+
if ( thePresamplingXCombs[theXCombBackup]->
dSigDR(make_pair(0.0,0.0), bdim,
&thePresamplingPoint[0]) == ZERO ) {
if ( projectionDipole()->verbose() )
generator()->log() << "Born outside phase space\n";
return 0.0;
}
}
const double * rad = &p[theBornRandom.second];
tStdDependentXCombPtr depXComb =
dynamic_ptr_cast<tStdDependentXCombPtr>(lastXCombPtr());
assert(depXComb);
- depXComb->setProcess();
if ( !projectionDipole()->generateKinematics(rad) ) {
if ( projectionDipole()->verbose() )
generator()->log() << "failed to generate radiation kinematics\n";
return 0.0;
}
- depXComb->remakeIncoming();
+ depXComb->clean();
depXComb->setIncomingPartons();
double res = evaluate();
return res >= 0. ? res : 0.;
} catch (...) {
if ( presampling() )
stopPresampling();
throw;
}
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void PowhegSplittingKernel::persistentOutput(PersistentOStream & os) const {
os << theBornScreening << thePresamplingPoints << theMaxTry << ounit(theScreeningScale,GeV)
<< theBornRandom << theRadiationRandom << theLastParameterPoint
<< thePresamplingPoint << thePresampling;
}
void PowhegSplittingKernel::persistentInput(PersistentIStream & is, int) {
is >> theBornScreening >> thePresamplingPoints >> theMaxTry >> iunit(theScreeningScale,GeV)
>> theBornRandom >> theRadiationRandom >> theLastParameterPoint
>> thePresamplingPoint >> thePresampling;
}
void PowhegSplittingKernel::Init() {
static ClassDocumentation<PowhegSplittingKernel> documentation
("PowhegSplittingKernel");
}
// *** 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<PowhegSplittingKernel,ME2byDipoles>
describeHerwigPowhegSplittingKernel("Herwig::PowhegSplittingKernel", "HwMatchbox.so");

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:33 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805843
Default Alt Text
(38 KB)

Event Timeline