No OneTemporary
Actions

Size

36 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
	--- a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
	+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
	@@ -1,481 +1,486 @@
	// -- C++ --
	//
	// ShowerApproximationGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ShowerApproximationGenerator class.
	//

	#include <config.h>
	#include "ShowerApproximationGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/Cuts/Cuts.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;

	ShowerApproximationGenerator::ShowerApproximationGenerator()
	: thePresamplingPoints(2000), theMaxTry(100000), theFreezeGrid(500000),
	theDoCompensate(false) {}

	ShowerApproximationGenerator::~ShowerApproximationGenerator() {}

	double ShowerApproximationGenerator::generateFraction(tcPDPtr pd, double r, double xmin) const {
	if ( pd->coloured() \|\| pd->id() == ParticleID::gamma ) {
	return pow(xmin,r);
	}
	double x0 = 1.e-5;
	return 1. + x0 - x0*pow((1.+x0)/x0,r);
	}

	double ShowerApproximationGenerator::invertFraction(tcPDPtr pd, double x, double xmin) const {
	if ( pd->coloured() \|\| pd->id() == ParticleID::gamma ) {
	return log(x)/log(xmin);
	}
	double x0 = 1.e-5;
	return log((1.-x+x0)/x0)/log((1.+x0)/x0);
	}

	bool ShowerApproximationGenerator::prepare() {

	tSubProPtr oldSub = lastIncomingXComb->subProcess();

	tcStdXCombPtr cIncomingXC = lastIncomingXComb;

	bool hasFractions =
	thePhasespace->haveX1X2() \|\|
	cIncomingXC->mePartonData().size() == 3;

	theLastMomenta.resize(cIncomingXC->mePartonData().size());

	if ( !hasFractions )
	theLastRandomNumbers.resize(thePhasespace->nDim(cIncomingXC->mePartonData()) + 2);
	else
	theLastRandomNumbers.resize(thePhasespace->nDim(cIncomingXC->mePartonData()));

	if ( !hasFractions ) {

	double x1 =
	oldSub->incoming().first->momentum().plus() /
	lastIncomingXComb->lastParticles().first->momentum().plus();

	theLastRandomNumbers[0] = invertFraction(oldSub->incoming().first->dataPtr(),x1,
	lastIncomingXComb->cuts()->x1Min());

	double x2 =
	oldSub->incoming().second->momentum().minus() /
	lastIncomingXComb->lastParticles().second->momentum().minus();

	theLastRandomNumbers[1] = invertFraction(oldSub->incoming().second->dataPtr(),x2,
	lastIncomingXComb->cuts()->x2Min());

	}

	theLastMomenta = cIncomingXC->meMomenta();

	theLastPartons.first =
	oldSub->incoming().first->data().produceParticle(oldSub->incoming().first->momentum());
	theLastPartons.second =
	oldSub->incoming().second->data().produceParticle(oldSub->incoming().second->momentum());

	thePhasespace->setXComb(lastIncomingXComb);
	thePhasespace->invertKinematics(theLastMomenta,
	!hasFractions ? &theLastRandomNumbers[2] : &theLastRandomNumbers[0]);

	theLastBornXComb->clean();
	theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
	theLastMomenta,theLastRandomNumbers);

	if ( !theLastBornXComb->cuts()->initSubProcess(theLastBornXComb->lastSHat(),
	theLastBornXComb->lastY(),
	theLastBornXComb->mirror()) )
	return false;

	theLastBornME->setXComb(theLastBornXComb);

	if ( !theLastBornME->generateKinematics(!hasFractions ? &theLastRandomNumbers[2] : &theLastRandomNumbers[0]) )
	return false;

	CrossSection bornXS = theLastBornME->dSigHatDR();

	if ( bornXS == ZERO )
	return false;

	return true;

	}

	bool ShowerApproximationGenerator::generate(const vector<double>& r) {

	theLastBornXComb->clean();

	bool hasFractions =
	thePhasespace->haveX1X2() \|\|
	theLastBornXComb->mePartonData().size() == 3;

	if ( !hasFractions ) {

	double x = generateFraction(theLastPartons.first->dataPtr(),r[0],
	lastIncomingXComb->cuts()->x1Min());
	Energy Q = lastIncomingXComb->lastParticles().first->momentum().plus();
	Energy mass = theLastPartons.first->dataPtr()->mass();
	double xi = (sqr(xQ) - sqr(mass))/(sqr(Q)x);
	Lorentz5Momentum p1(ZERO,ZERO,xi*Q/2.);
	p1.setMass(mass); p1.rescaleEnergy();
	theLastPartons.first->set5Momentum(p1);

	x = generateFraction(theLastPartons.second->dataPtr(),r[1],
	lastIncomingXComb->cuts()->x2Min());
	Q = lastIncomingXComb->lastParticles().second->momentum().minus();
	mass = theLastPartons.second->dataPtr()->mass();
	xi = (sqr(xQ) - sqr(mass))/(sqr(Q)x);
	Lorentz5Momentum p2(ZERO,ZERO,-xi*Q/2.);
	p2.setMass(mass); p2.rescaleEnergy();
	theLastPartons.second->set5Momentum(p2);

	} else {

	theLastBornME->setXComb(theLastBornXComb);
	theLastBornXComb->lastParticles(lastIncomingXComb->lastParticles());
	theLastBornXComb->lastP1P2(make_pair(0.0, 0.0));
	theLastBornXComb->lastS(lastIncomingXComb->lastS());

	if ( !theLastBornME->generateKinematics(&r[0]) )
	return false;

	theLastPartons.first->set5Momentum(theLastBornME->lastMEMomenta()[0]);
	theLastPartons.second->set5Momentum(theLastBornME->lastMEMomenta()[1]);

	}

	theLastPresamplingMomenta.resize(theLastMomenta.size());

	Boost toCMS =
	(theLastPartons.first->momentum() +
	theLastPartons.second->momentum()).findBoostToCM();

	theLastPresamplingMomenta[0] = theLastPartons.first->momentum();
	if ( !hasFractions )
	theLastPresamplingMomenta[0].boost(toCMS);
	theLastPresamplingMomenta[1] = theLastPartons.second->momentum();
	if ( !hasFractions )
	theLastPresamplingMomenta[1].boost(toCMS);

	if ( hasFractions ) {
	for ( size_t k = 2; k < theLastBornME->lastMEMomenta().size(); ++k )
	theLastPresamplingMomenta[k] = theLastBornME->lastMEMomenta()[k];
	}

	theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
	theLastPresamplingMomenta,r);

	if ( !theLastBornXComb->cuts()->initSubProcess(theLastBornXComb->lastSHat(),
	theLastBornXComb->lastY(),
	theLastBornXComb->mirror()) )
	return false;

	if ( !hasFractions ) {
	theLastBornME->setXComb(theLastBornXComb);
	if ( !theLastBornME->generateKinematics(&r[2]) )
	return false;
	}

	CrossSection bornXS = theLastBornME->dSigHatDR();

	if ( bornXS == ZERO )
	return false;

	return true;

	}

	void ShowerApproximationGenerator::restore() {

	tSubProPtr oldSub = lastIncomingXComb->subProcess();

	theLastPartons.first->set5Momentum(oldSub->incoming().first->momentum());
	theLastPartons.second->set5Momentum(oldSub->incoming().second->momentum());

	theLastBornXComb->clean();
	theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
	theLastMomenta,theLastRandomNumbers);

	theLastBornME->setXComb(theLastBornXComb);
	theLastBornME->generateKinematics(&theLastRandomNumbers[2]);
	theLastBornME->dSigHatDR();

	}

	void ShowerApproximationGenerator::
	handle(EventHandler & eh, const tPVector &,
	const Hint &) {
	theFactory->setHardTreeEmitter(-1);
	theFactory->setHardTreeSpectator(-1);
	theFactory->setHardTreeSubprocess(SubProPtr());


	lastIncomingXComb = dynamic_ptr_cast<tStdXCombPtr>(eh.lastXCombPtr());
	if ( !lastIncomingXComb )
	throw Exception() << "ShowerApproximationGenerator::handle(): Expecting a standard event handler."
	<< Exception::runerror;

	if ( lastIncomingXComb->lastProjector() )
	lastIncomingXComb = lastIncomingXComb->lastProjector();

	const StandardXComb& xc = *lastIncomingXComb;

	map<cPDVector,set<Ptr<ShowerApproximationKernel>::ptr> >::const_iterator
	kernelit = theKernelMap.find(xc.mePartonData());

	if ( kernelit == theKernelMap.end() ) {
	list<MatchboxFactory::SplittingChannel> channels =
	theFactory->getSplittingChannels(lastIncomingXComb);
	set<Ptr<ShowerApproximationKernel>::ptr> newKernels;
	for ( list<MatchboxFactory::SplittingChannel>::const_iterator c =
	channels.begin(); c != channels.end(); ++c ) {
	Ptr<ShowerApproximationKernel>::ptr kernel =
	new_ptr(ShowerApproximationKernel());
	kernel->setBornXComb(c->bornXComb);
	kernel->setRealXComb(c->realXComb);
	kernel->setTildeXCombs(c->tildeXCombs);
	kernel->setDipole(c->dipole);
	kernel->showerApproximation(theShowerApproximation);
	kernel->presamplingPoints(thePresamplingPoints);
	kernel->maxtry(theMaxTry);
	kernel->freezeGrid(theFreezeGrid);
	kernel->showerApproximationGenerator(this);
	kernel->doCompensate(theDoCompensate);
	if ( kernel->dipole()->bornEmitter() > 1 &&
	kernel->dipole()->bornSpectator() > 1 ) {
	kernel->ptCut(ffPtCut());
	} else if ( ( kernel->dipole()->bornEmitter() > 1 &&
	kernel->dipole()->bornSpectator() < 2 ) \|\|
	( kernel->dipole()->bornEmitter() < 2 &&
	kernel->dipole()->bornSpectator() > 1 ) ) {
	kernel->ptCut(fiPtCut());
	} else {
	assert(kernel->dipole()->bornEmitter() < 2 &&
	kernel->dipole()->bornSpectator() < 2);
	kernel->ptCut(iiPtCut());
	}
	newKernels.insert(kernel);
	}
	theKernelMap[xc.mePartonData()] = newKernels;
	kernelit = theKernelMap.find(xc.mePartonData());
	}

	if ( kernelit->second.empty() )
	return;

	const set<Ptr<ShowerApproximationKernel>::ptr>& kernels = kernelit->second;

	theLastBornME = (**kernels.begin()).dipole()->underlyingBornME();
	theLastBornME->phasespace(thePhasespace);
	theLastBornXComb = (**kernels.begin()).bornXComb();

	if ( !prepare() )
	return;

	Energy winnerPt = ZERO;
	Ptr<ShowerApproximationKernel>::ptr winnerKernel;

	- for ( set<Ptr<ShowerApproximationKernel>::ptr>::const_iterator k =
	- kernels.begin(); k != kernels.end(); ++k ) {
	- if ( (*k).generate() != 0. && (k)->dipole()->lastPt() > winnerPt){
	- winnerKernel = *k;
	- winnerPt = winnerKernel->dipole()->lastPt();
	+ try {
	+ for ( set<Ptr<ShowerApproximationKernel>::ptr>::const_iterator k =
	+ kernels.begin(); k != kernels.end(); ++k ) {
	+ if ( (*k).generate() != 0. && (k)->dipole()->lastPt() > winnerPt){
	+ winnerKernel = *k;
	+ winnerPt = winnerKernel->dipole()->lastPt();
	+ }
	}
	+ } catch(ShowerApproximationKernel::MaxTryException&) {
	+ throw Exception() << "Too many tries needed to generate the matrix element correction in '"
	+ << name() << "'" << Exception::eventerror;
	}

	if ( !winnerKernel \|\| winnerPt == ZERO )
	return;

	//Hardest emission should be this one.

	winnerKernel->realXComb()->lastCentralScale(sqr(winnerPt));
	winnerKernel->bornXComb()->lastCentralScale(sqr(winnerPt));
	lastIncomingXComb->lastCentralScale(sqr(winnerPt));


	SubProPtr oldSub = lastIncomingXComb->subProcess();
	SubProPtr newSub;

	try {
	tcDiagPtr bornDiag = lastIncomingXComb->lastDiagram();
	tcDiagPtr realDiag =
	winnerKernel->dipole()->realEmissionDiagram(bornDiag);
	winnerKernel->realXComb()->externalDiagram(realDiag);
	newSub = winnerKernel->realXComb()->construct();
	} catch(Veto&) {
	return;
	}

	if ( !theShowerApproximation->needsTruncatedShower() ){
	tParticleSet firstS = oldSub->incoming().first->siblings();
	assert(firstS.empty() \|\| firstS.size() == 1);
	if ( !firstS.empty() ) {
	eh.currentStep()->removeParticle(*firstS.begin());
	}

	tParticleSet secondS = oldSub->incoming().second->siblings();
	assert(secondS.empty() \|\| secondS.size() == 1);
	if ( !secondS.empty() ) {
	eh.currentStep()->removeParticle(*secondS.begin());
	}

	// prevent the colliding particles from disappearing
	// in the initial state and appearing
	// in the final state when we've cut off all their
	// (physical) children; only applies to the case
	// where we have a parton extractor not build from
	// noPDF, so check wether the incoming particle
	// doesnt equal the incoming parton -- this needs fixing in ThePEG
	PPtr dummy = new_ptr(Particle(getParticleData(ParticleID::gamma)));
	bool usedDummy = false;
	if ( eh.currentStep()->incoming().first != oldSub->incoming().first ) {
	eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().first,dummy);
	usedDummy = true;
	}
	if ( eh.currentStep()->incoming().second != oldSub->incoming().second ) {
	eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().second,dummy);
	usedDummy = true;
	}

	eh.currentStep()->removeSubProcess(oldSub);
	eh.currentStep()->addSubProcess(newSub);

	// get rid of the dummy
	if ( usedDummy ) {
	eh.currentStep()->removeParticle(dummy);
	}

	eh.select(winnerKernel->realXComb());

	winnerKernel->realXComb()->recreatePartonBinInstances(winnerKernel->realXComb()->lastScale());
	winnerKernel->realXComb()->refillPartonBinInstances(&(xc.lastRandomNumbers()[0]));

	winnerKernel->realXComb()->pExtractor()->constructRemnants(winnerKernel->realXComb()->partonBinInstances(),
	newSub, eh.currentStep());
	}
	else{
	theFactory->setHardTreeSubprocess(newSub);
	theFactory->setHardTreeEmitter(winnerKernel->dipole()->bornEmitter());
	theFactory->setHardTreeSpectator(winnerKernel->dipole()->bornSpectator());
	}

	}

	IBPtr ShowerApproximationGenerator::clone() const {
	return new_ptr(*this);
	}

	IBPtr ShowerApproximationGenerator::fullclone() const {
	return new_ptr(*this);
	}


	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void ShowerApproximationGenerator::persistentOutput(PersistentOStream & os) const {
	os << theShowerApproximation << thePhasespace << theFactory
	<< theKernelMap << thePresamplingPoints << theMaxTry << theFreezeGrid
	<< lastIncomingXComb << theLastBornME << ounit(theLastMomenta,GeV)
	<< ounit(theLastPresamplingMomenta,GeV) << theLastRandomNumbers
	<< theLastBornXComb << theLastPartons << theDoCompensate;
	}

	void ShowerApproximationGenerator::persistentInput(PersistentIStream & is, int) {
	is >> theShowerApproximation >> thePhasespace >> theFactory
	>> theKernelMap >> thePresamplingPoints >> theMaxTry >> theFreezeGrid
	>> lastIncomingXComb >> theLastBornME >> iunit(theLastMomenta,GeV)
	>> iunit(theLastPresamplingMomenta,GeV) >> theLastRandomNumbers
	>> theLastBornXComb >> theLastPartons >> theDoCompensate;
	}


	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<ShowerApproximationGenerator,StepHandler>
	describeHerwigShowerApproximationGenerator("Herwig::ShowerApproximationGenerator", "Herwig.so");

	void ShowerApproximationGenerator::Init() {

	static ClassDocumentation<ShowerApproximationGenerator> documentation
	("ShowerApproximationGenerator generates emissions according to a "
	"shower approximation entering a NLO matching.");


	static Reference<ShowerApproximationGenerator,ShowerApproximation> interfaceShowerApproximation
	("ShowerApproximation",
	"Set the shower approximation to sample.",
	&ShowerApproximationGenerator::theShowerApproximation, false, false, true, false, false);


	static Reference<ShowerApproximationGenerator,MatchboxPhasespace> interfacePhasespace
	("Phasespace",
	"The phase space generator to use.",
	&ShowerApproximationGenerator::thePhasespace, false, false, true, false, false);


	static Reference<ShowerApproximationGenerator,MatchboxFactory> interfaceFactory
	("Factory",
	"The factory object to use.",
	&ShowerApproximationGenerator::theFactory, false, false, true, false, false);

	static Parameter<ShowerApproximationGenerator,unsigned long> interfacePresamplingPoints
	("PresamplingPoints",
	"Set the number of presampling points.",
	&ShowerApproximationGenerator::thePresamplingPoints, 2000, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximationGenerator,unsigned long> interfaceMaxTry
	("MaxTry",
	"Set the number of maximum attempts.",
	&ShowerApproximationGenerator::theMaxTry, 100000, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<ShowerApproximationGenerator,unsigned long> interfaceFreezeGrid
	("FreezeGrid",
	"",
	&ShowerApproximationGenerator::theFreezeGrid, 500000, 1, 0,
	false, false, Interface::lowerlim);

	static Switch<ShowerApproximationGenerator,bool> interfaceDoCompensate
	("DoCompensate",
	"",
	&ShowerApproximationGenerator::theDoCompensate, false, false, false);
	static SwitchOption interfaceDoCompensateYes
	(interfaceDoCompensate,
	"Yes",
	"",
	true);
	static SwitchOption interfaceDoCompensateNo
	(interfaceDoCompensate,
	"No",
	"",
	false);

	}

	diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
	--- a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
	+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
	@@ -1,210 +1,210 @@
	// -- C++ --
	//
	// ShowerApproximationKernel.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the ShowerApproximationKernel class.
	//

	#include <config.h>
	#include "ShowerApproximationKernel.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	#include "ShowerApproximationGenerator.h"

	using namespace Herwig;

	ShowerApproximationKernel::ShowerApproximationKernel()
	: thePresampling(false), thePresamplingPoints(10000),
	theMaxTry(100000), theFreezeGrid(500000),
	sampler(0), theDoCompensate(false) {}

	ShowerApproximationKernel::~ShowerApproximationKernel() {}

	IBPtr ShowerApproximationKernel::clone() const {
	return new_ptr(*this);
	}

	IBPtr ShowerApproximationKernel::fullclone() const {
	return new_ptr(*this);
	}

	void ShowerApproximationKernel::showerApproximationGenerator(Ptr<ShowerApproximationGenerator>::tptr gen) {
	theShowerApproximationGenerator = gen;
	}

	Ptr<ShowerApproximationGenerator>::tptr ShowerApproximationKernel::showerApproximationGenerator() const {
	return theShowerApproximationGenerator;
	}

	const vector<bool>& ShowerApproximationKernel::sampleFlags() {
	if ( !theFlags.empty() )
	return theFlags;
	theFlags.resize(nDim(),false);
	for ( int k = nDimBorn();
	k < nDimBorn() + dipole()->nDimRadiation(); ++k )
	theFlags[k] = true;
	return theFlags;
	}

	const pair<vector<double>,vector<double> >& ShowerApproximationKernel::support() {
	if ( !theSupport.first.empty() )
	return theSupport;
	vector<double> l(nDim(),0.0);
	vector<double> u(nDim(),1.0);
	theSupport.first = l;
	theSupport.second = u;
	return theSupport;
	}

	const vector<double>& ShowerApproximationKernel::parameterPoint() {
	theLastParameterPoint.resize(nDim());
	copy(bornCXComb()->lastRandomNumbers().begin(),
	bornCXComb()->lastRandomNumbers().end(),
	theLastParameterPoint.begin());
	theLastParameterPoint[evolutionVariable()] = 1.;
	return theLastParameterPoint;
	}

	void ShowerApproximationKernel::startPresampling() {
	thePresampling = true;
	}

	void ShowerApproximationKernel::stopPresampling() {
	showerApproximationGenerator()->restore();
	thePresampling = false;
	}

	double ShowerApproximationKernel::evaluate(const vector<double>& r) {

	if ( presampling() ) {

	theLastBornPoint.resize(nDimBorn());
	copy(r.begin(),r.begin()+nDimBorn(),theLastBornPoint.begin());
	if ( !showerApproximationGenerator()->generate(theLastBornPoint) )
	return 0.;

	}

	assert(dipole()->splitting());
	realXComb()->clean();
	dipole()->setXComb(realXComb());
	for ( vector<StdXCombPtr>::const_iterator t = tildeXCombs().begin();
	t != tildeXCombs().end(); ++t ) {
	(**t).clean();
	(*t).matrixElement()->setXComb(t);
	}
	if ( !dipole()->generateKinematics(&r[nDimBorn()]) )
	return 0.;

	double jac =
	showerApproximation()->showerInvertedTildeKinematics() ?
	showerApproximation()->showerInvertedTildeKinematics()->jacobian() :
	dipole()->invertedTildeKinematics()->jacobian();

	showerApproximation()->setBornXComb(bornXComb());
	showerApproximation()->setRealXComb(realXComb());
	showerApproximation()->setTildeXCombs(tildeXCombs());
	showerApproximation()->setDipole(dipole());
	showerApproximation()->checkCutoff();
	showerApproximation()->getShowerVariables();

	if ( !dipole()->isInShowerPhasespace() )
	return 0.;

	return showerApproximation()->me2() * jac;

	}

	double ShowerApproximationKernel::generate() {

	if ( !sampler ) {

	sampler = new ExponentialGenerator();
	sampler->sampling_parameters().maxtry = maxtry();
	sampler->sampling_parameters().presampling_points = presamplingPoints();
	sampler->sampling_parameters().freeze_grid = freezeGrid();
	sampler->docompensate(theDoCompensate);
	sampler->function(this);
	sampler->initialize();

	}

	double res = 0.;

	while (true) {
	try {
	res = sampler->generate();
	} catch (exsample::exponential_regenerate&) {
	continue;
	} catch (exsample::hit_and_miss_maxtry&) {
	- throw Veto();
	+ throw MaxTryException();
	} catch (exsample::selection_maxtry&) {
	- throw Veto();
	+ throw MaxTryException();
	}
	break;
	}

	return res;

	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void ShowerApproximationKernel::persistentOutput(PersistentOStream & os) const {
	os << theDipole << theShowerApproximation << theBornXComb
	<< theRealXComb << theTildeXCombs << thePresampling
	<< thePresamplingPoints << theMaxTry << theFreezeGrid << theFlags
	<< theSupport << theShowerApproximationGenerator
	<< theLastParameterPoint << theLastBornPoint
	<< (sampler ? true : false) << theDoCompensate;
	if ( sampler )
	sampler->put(os);
	}

	void ShowerApproximationKernel::persistentInput(PersistentIStream & is, int) {
	bool haveSampler;
	is >> theDipole >> theShowerApproximation >> theBornXComb
	>> theRealXComb >> theTildeXCombs >> thePresampling
	>> thePresamplingPoints >> theMaxTry >> theFreezeGrid >> theFlags
	>> theSupport >> theShowerApproximationGenerator
	>> theLastParameterPoint >> theLastBornPoint
	>> haveSampler >> theDoCompensate;
	if ( haveSampler ) {
	sampler = new ExponentialGenerator();
	sampler->get(is);
	sampler->function(this);
	sampler->initialize();
	}
	}

	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeClass<ShowerApproximationKernel,HandlerBase>
	describeHerwigShowerApproximationKernel("Herwig::ShowerApproximationKernel", "Herwig.so");

	void ShowerApproximationKernel::Init() {

	static ClassDocumentation<ShowerApproximationKernel> documentation
	("ShowerApproximationKernel generates emissions according to a "
	"shower approximation entering a NLO matching.");

	}

	diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
	--- a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
	+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.h
	@@ -1,455 +1,462 @@
	// -- C++ --
	//
	// ShowerApproximationKernel.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef Herwig_ShowerApproximationKernel_H
	#define Herwig_ShowerApproximationKernel_H
	//
	// This is the declaration of the ShowerApproximationKernel class.
	//

	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
	#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
	#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
	#include "Herwig/Sampling/exsample/exponential_generator.h"

	namespace Herwig {

	using namespace ThePEG;

	class ShowerApproximationGenerator;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief ShowerApproximationKernel generates emissions according to a
	* shower approximation entering a NLO matching.
	*
	*/
	class ShowerApproximationKernel: public HandlerBase {

	public:

	+ /**
	+ * Exception to communicate sampler maxtry events.
	+ */
	+ struct MaxTryException {};
	+
	+public:
	+
	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	ShowerApproximationKernel();

	/**
	* The destructor.
	*/
	virtual ~ShowerApproximationKernel();
	//@}

	public:

	/**
	* Set the XComb object describing the Born process
	*/
	void setBornXComb(tStdXCombPtr xc) { theBornXComb = xc; }

	/**
	* Return the XComb object describing the Born process
	*/
	tcStdXCombPtr bornCXComb() const { return theBornXComb; }

	/**
	* Return the XComb object describing the Born process
	*/
	tStdXCombPtr bornXComb() const { return theBornXComb; }

	/**
	* Set the XComb object describing the real emission process
	*/
	void setRealXComb(tStdXCombPtr xc) { theRealXComb = xc; }

	/**
	* Return the XComb object describing the real emission process
	*/
	tcStdXCombPtr realCXComb() const { return theRealXComb; }

	/**
	* Return the XComb object describing the real emission process
	*/
	tStdXCombPtr realXComb() const { return theRealXComb; }

	/**
	* Set the tilde xcomb objects associated to the real xcomb
	*/
	void setTildeXCombs(const vector<StdXCombPtr>& xc) { theTildeXCombs = xc; }

	/**
	* Return the tilde xcomb objects associated to the real xcomb
	*/
	const vector<StdXCombPtr>& tildeXCombs() const { return theTildeXCombs; }

	/**
	* Set the dipole in charge for the emission
	*/
	void setDipole(Ptr<SubtractionDipole>::tptr dip) { theDipole = dip; }

	/**
	* Return the dipole in charge for the emission
	*/
	Ptr<SubtractionDipole>::tptr dipole() const { return theDipole; }

	/**
	* Set the shower approximation.
	*/
	void showerApproximation(Ptr<ShowerApproximation>::tptr app) { theShowerApproximation = app; }

	/**
	* Return the shower approximation.
	*/
	Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }

	/**
	* Set the shower approximation generator.
	*/
	void showerApproximationGenerator(Ptr<ShowerApproximationGenerator>::tptr);

	/**
	* Return the shower approximation generator.
	*/
	Ptr<ShowerApproximationGenerator>::tptr showerApproximationGenerator() const;

	/**
	* Generate the next emission
	*/
	double generate();

	public:

	/**
	* Set a pt cut on the dipole to generate the radiation
	*/
	void ptCut(Energy pt) { dipole()->ptCut(pt); }

	/**
	* Return the number of random numbers
	* needed to sample this kernel.
	*/
	int nDim() const {
	return
	nDimBorn() +
	dipole()->nDimRadiation();
	}

	/**
	* Return the number of random numbers
	* needed to sample the Born process.
	*/
	int nDimBorn() const {
	return bornCXComb()->lastRandomNumbers().size();
	}

	/**
	* Flag, which variables are free variables.
	*/
	const vector<bool>& sampleFlags();

	/**
	* Return the support of the splitting kernel.
	* The lower bound on the first variable is
	* assumed to correspond to the cutoff on the
	* evolution variable.
	*/
	const pair<vector<double>,vector<double> >& support();

	/**
	* Return the parameter point associated to the splitting
	* previously supplied through fixParameters.
	*/
	const vector<double>& parameterPoint();

	/**
	* Indicate that presampling of this kernel
	* will be performed in the next calls to
	* evaluate until stopPresampling() is called.
	*/
	void startPresampling();

	/**
	* Indicate that presampling of this kernel
	* is done until startPresampling() is called.
	*/
	void stopPresampling();

	/**
	* Return true, if currently being presampled
	*/
	bool presampling() const { return thePresampling; }

	/**
	* Return the number of points to presample this
	* splitting generator.
	*/
	unsigned long presamplingPoints() const { return thePresamplingPoints; }

	/**
	* Return the maximum number of trials
	* to generate a splitting.
	*/
	unsigned long maxtry() const { return theMaxTry; }

	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long freezeGrid() const { return theFreezeGrid; }

	/**
	* Set the number of points to presample this
	* splitting generator.
	*/
	void presamplingPoints(unsigned long p) { thePresamplingPoints = p; }

	/**
	* Set the maximum number of trials
	* to generate a splitting.
	*/
	void maxtry(unsigned long p) { theMaxTry = p; }

	/**
	* Set the number of accepted points after which the grid should
	* be frozen
	*/
	void freezeGrid(unsigned long n) { theFreezeGrid = n; }

	/**
	* Evalute the splitting kernel.
	*/
	double evaluate(const vector<double>&);

	/**
	* Return the index of the random number corresponding
	* to the evolution variable.
	*/
	int evolutionVariable() const {
	return
	nDimBorn() +
	(showerApproximation()->showerInvertedTildeKinematics() ?
	showerApproximation()->showerInvertedTildeKinematics()->evolutionVariable() :
	dipole()->invertedTildeKinematics()->evolutionVariable());
	}

	/**
	* Return the cutoff on the evolution
	* random number corresponding to the pt cut.
	*/
	double evolutionCutoff() const {
	return
	showerApproximation()->showerInvertedTildeKinematics() ?
	showerApproximation()->showerInvertedTildeKinematics()->evolutionCutoff() :
	dipole()->invertedTildeKinematics()->evolutionCutoff();
	}

	/**
	* True, if sampler should apply compensation
	*/
	void doCompensate(bool yes = true) { theDoCompensate = yes; }

	public:

	/*@name Wrap to the exsample2 interface until this is finally cleaned up. /
	//@{

	inline const vector<bool>& variable_flags () {
	return sampleFlags();
	}

	inline size_t evolution_variable () const { return evolutionVariable(); }

	inline double evolution_cutoff () const {
	return evolutionCutoff();
	}

	inline const vector<double>& parameter_point () {
	return parameterPoint();
	}

	inline void start_presampling () {
	startPresampling();
	}

	inline void stop_presampling () {
	stopPresampling();
	}

	inline size_t dimension () const {
	return nDim();
	}

	inline unsigned long presampling_points () const {
	return presamplingPoints();
	}

	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}


	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).


	private:

	/**
	* The dipole in charge of the emission
	*/
	Ptr<SubtractionDipole>::ptr theDipole;

	/**
	* The shower approximation to consider
	*/
	Ptr<ShowerApproximation>::ptr theShowerApproximation;

	/**
	* The XComb off which radiation will be generated
	*/
	StdXCombPtr theBornXComb;

	/**
	* The XComb describing the process after radiation
	*/
	StdXCombPtr theRealXComb;

	/**
	* The tilde xcomb objects associated to the real xcomb
	*/
	vector<StdXCombPtr> theTildeXCombs;

	/**
	* True, if currently being presampled
	*/
	bool thePresampling;

	/**
	* The number of points to presample this
	* splitting generator.
	*/
	unsigned long thePresamplingPoints;

	/**
	* The maximum number of trials
	* to generate a splitting.
	*/
	unsigned long theMaxTry;

	/**
	* Return the number of accepted points after which the grid should
	* be frozen
	*/
	unsigned long theFreezeGrid;

	/**
	* The sampling flags
	*/
	vector<bool> theFlags;

	/**
	* The support.
	*/
	pair<vector<double>,vector<double> > theSupport;

	/**
	* The shower approximation generator.
	*/
	Ptr<ShowerApproximationGenerator>::tptr theShowerApproximationGenerator;

	/**
	* The last parameter point
	*/
	vector<double> theLastParameterPoint;

	/**
	* The last random numbers used for Born sampling
	*/
	vector<double> theLastBornPoint;

	/**
	* Define the Sudakov sampler
	*/
	typedef
	exsample::exponential_generator<ShowerApproximationKernel,UseRandom>
	ExponentialGenerator;

	/**
	* Define a pointer to the Sudakov sampler
	*/
	typedef
	exsample::exponential_generator<ShowerApproximationKernel,UseRandom>*
	ExponentialGeneratorPtr;

	/**
	* The Sudakov sampler
	*/
	ExponentialGeneratorPtr sampler;

	/**
	* True, if sampler should apply compensation
	*/
	bool theDoCompensate;

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	ShowerApproximationKernel & operator=(const ShowerApproximationKernel &);

	};

	}

	#endif /* Herwig_ShowerApproximationKernel_H */

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 8:54 PM (23 h, 2 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242346
Default Alt Text: (36 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions