No OneTemporary
Actions

Size

38 KB

Subscribers

None

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.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<PowhegSplittingGenerator,Energy> interfaceFFScreeningScale
	("FFScreeningScale",
	"Set the screening scale",
	&PowhegSplittingGenerator::theFFScreeningScale, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<PowhegSplittingGenerator,Energy> interfaceFIPtCut
	("FIPtCut",
	"Set the pt infrared cutoff",
	&PowhegSplittingGenerator::theFIPtCut, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<PowhegSplittingGenerator,Energy> interfaceFIScreeningScale
	("FIScreeningScale",
	"Set the screening scale",
	&PowhegSplittingGenerator::theFIScreeningScale, GeV, 0.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<PowhegSplittingGenerator,Energy> interfaceIIPtCut
	("IIPtCut",
	"Set the pt infrared cutoff",
	&PowhegSplittingGenerator::theIIPtCut, GeV, 1.0GeV, 0.0GeV, 0*GeV,
	false, false, Interface::lowerlim);

	static Parameter<PowhegSplittingGenerator,Energy> interfaceIIScreeningScale
	("IIScreeningScale",
	"Set the screening scale",
	&PowhegSplittingGenerator::theIIScreeningScale, GeV, 0.0GeV, 0.0GeV, 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, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805843
Default Alt Text: (38 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions