Page MenuHomeHEPForge
Files F7878340

No OneTemporary
Actions

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/DipoleShower/Base/DipoleSplittingGenerator.cc b/DipoleShower/Base/DipoleSplittingGenerator.cc
--- a/DipoleShower/Base/DipoleSplittingGenerator.cc
+++ b/DipoleShower/Base/DipoleSplittingGenerator.cc
@@ -1,555 +1,553 @@
// -*- C++ -*-
//
// DipoleSplittingGenerator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleSplittingGenerator class.
//
#include "DipoleSplittingGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/DipoleShower/DipoleShowerHandler.h"
using namespace Herwig;
DipoleSplittingGenerator::DipoleSplittingGenerator()
: HandlerBase(),
theExponentialGenerator(0), prepared(false), presampling(false) {
if ( ShowerHandler::currentHandler() )
setGenerator(ShowerHandler::currentHandler()->generator());
}
DipoleSplittingGenerator::~DipoleSplittingGenerator() {
if ( theExponentialGenerator ) {
delete theExponentialGenerator;
theExponentialGenerator = 0;
}
}
IBPtr DipoleSplittingGenerator::clone() const {
return new_ptr(*this);
}
IBPtr DipoleSplittingGenerator::fullclone() const {
return new_ptr(*this);
}
void DipoleSplittingGenerator::wrap(Ptr<DipoleSplittingGenerator>::ptr other) {
assert(!prepared);
theOtherGenerator = other;
}
void DipoleSplittingGenerator::prepare(const DipoleSplittingInfo& sp) {
generatedSplitting = sp;
generatedSplitting.splittingKinematics(splittingKernel()->splittingKinematics());
generatedSplitting.splittingParameters().resize(splittingKernel()->nDimAdditional());
if ( wrapping() ) {
generatedSplitting.emitterData(theSplittingKernel->emitter(generatedSplitting.index()));
generatedSplitting.spectatorData(theSplittingKernel->spectator(generatedSplitting.index()));
generatedSplitting.emissionData(theSplittingKernel->emission(generatedSplitting.index()));
parameters.resize(theOtherGenerator->nDim());
prepared = true;
return;
}
generatedSplitting.emitterData(splittingKernel()->emitter(generatedSplitting.index()));
generatedSplitting.spectatorData(splittingKernel()->spectator(generatedSplitting.index()));
generatedSplitting.emissionData(splittingKernel()->emission(generatedSplitting.index()));
presampledSplitting = generatedSplitting;
prepared = true;
parameters.resize(nDim());
theExponentialGenerator =
new exsample::exponential_generator<DipoleSplittingGenerator,UseRandom>();
theExponentialGenerator->sampling_parameters().maxtry = maxtry();
theExponentialGenerator->sampling_parameters().presampling_points = presamplingPoints();
theExponentialGenerator->function(this);
theExponentialGenerator->initialize();
}
void DipoleSplittingGenerator::fixParameters(const DipoleSplittingInfo& sp) {
assert(generator());
assert(!presampling);
assert(prepared);
assert(sp.index() == generatedSplitting.index());
generatedSplitting.scale(sp.scale());
parameters[3] = sp.scale()/generator()->maximumCMEnergy();
Energy maxPossible =
generatedSplitting.splittingKinematics()->ptMax(sp.scale(),
sp.emitterX(), sp.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
if ( maxPossible >= sp.hardPt() )
generatedSplitting.hardPt(sp.hardPt());
else
generatedSplitting.hardPt(maxPossible);
parameters[0] = splittingKinematics()->ptToRandom(generatedSplitting.hardPt(),
sp.scale(),generatedSplitting.index());
size_t shift = 4;
if ( generatedSplitting.index().emitterPDF().pdf() &&
generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.emitterX(sp.emitterX());
generatedSplitting.spectatorX(sp.spectatorX());
parameters[4] = sp.emitterX();
parameters[5] = sp.spectatorX();
shift += 2;
}
if ( generatedSplitting.index().emitterPDF().pdf() &&
!generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.emitterX(sp.emitterX());
parameters[4] = sp.emitterX();
++shift;
}
if ( !generatedSplitting.index().emitterPDF().pdf() &&
generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.spectatorX(sp.spectatorX());
parameters[4] = sp.spectatorX();
++shift;
}
if ( splittingReweight() ) {
parameters[shift] = splittingReweight()->evaluate(sp);
++shift;
}
if ( splittingKernel()->nDimAdditional() )
copy(sp.lastSplittingParameters().begin(),sp.lastSplittingParameters().end(),parameters.begin()+shift);
if ( sp.emitter() )
generatedSplitting.emitter(sp.emitter());
if ( sp.spectator() )
generatedSplitting.spectator(sp.spectator());
}
int DipoleSplittingGenerator::nDim() const {
assert(!wrapping());
assert(prepared);
int ret = 4; // 0 pt, 1 z, 2 phi, 3 scale, 4/5 xs + parameters
if ( generatedSplitting.index().emitterPDF().pdf() ) {
++ret;
}
if ( generatedSplitting.index().spectatorPDF().pdf() ) {
++ret;
}
if ( splittingReweight() )
++ret;
ret += splittingKernel()->nDimAdditional();
return ret;
}
const vector<bool>& DipoleSplittingGenerator::sampleFlags() {
assert(!wrapping());
if ( !theFlags.empty() )
return theFlags;
theFlags.resize(nDim(),false);
theFlags[0] = true; theFlags[1] = true; theFlags[2] = true; // 0 pt, 1 z, 2 phi
return theFlags;
}
const pair<vector<double>,vector<double> >& DipoleSplittingGenerator::support() {
assert(!wrapping());
if ( !theSupport.first.empty() )
return theSupport;
vector<double> lower(nDim(),0.);
vector<double> upper(nDim(),1.);
pair<double,double> kSupport =
generatedSplitting.splittingKinematics()->kappaSupport(generatedSplitting);
pair<double,double> xSupport =
generatedSplitting.splittingKinematics()->xiSupport(generatedSplitting);
lower[0] = kSupport.first;
lower[1] = xSupport.first;
upper[0] = kSupport.second;
upper[1] = xSupport.second;
if ( splittingReweight() ) {
pair<double,double> bounds = splittingReweight()->reweightBounds(generatedSplitting.index());
int pos = 4;
if ( generatedSplitting.index().emitterPDF().pdf() ) {
++pos;
}
if ( generatedSplitting.index().spectatorPDF().pdf() ) {
++pos;
}
lower[pos] = bounds.first;
upper[pos] = bounds.second;
}
theSupport.first = lower;
theSupport.second = upper;
return theSupport;
}
void DipoleSplittingGenerator::startPresampling() {
assert(!wrapping());
splittingKernel()->startPresampling(generatedSplitting.index());
presampling = true;
}
void DipoleSplittingGenerator::stopPresampling() {
assert(!wrapping());
splittingKernel()->stopPresampling(generatedSplitting.index());
presampling = false;
}
bool DipoleSplittingGenerator::haveOverestimate() const {
assert(!wrapping());
assert(prepared);
return
generatedSplitting.splittingKinematics()->haveOverestimate() &&
splittingKernel()->haveOverestimate(generatedSplitting);
}
bool DipoleSplittingGenerator::overestimate(const vector<double>& point) {
assert(!wrapping());
assert(prepared);
assert(!presampling);
assert(haveOverestimate());
if ( ! generatedSplitting.splittingKinematics()->generateSplitting(point[0],point[1],point[2],
generatedSplitting) )
return 0.;
generatedSplitting.splittingKinematics()->prepareSplitting(generatedSplitting);
return
( generatedSplitting.splittingKinematics()->jacobianOverestimate() *
splittingKernel()->overestimate(generatedSplitting) *
(splittingReweight() ? splittingReweight()->evaluate(generatedSplitting) : 1.) );
}
double DipoleSplittingGenerator::invertOverestimateIntegral(double value) const {
assert(!wrapping());
assert(prepared);
assert(!presampling);
assert(haveOverestimate());
return
splittingKernel()->invertOverestimateIntegral(generatedSplitting,value);
}
double DipoleSplittingGenerator::evaluate(const vector<double>& point) {
assert(!wrapping());
assert(prepared);
assert(generator());
DipoleSplittingInfo& split =
( !presampling ? generatedSplitting : presampledSplitting );
split.continuesEvolving();
size_t shift = 4;
if ( presampling ) {
split.scale(point[3] * generator()->maximumCMEnergy());
if ( split.index().emitterPDF().pdf() &&
split.index().spectatorPDF().pdf() ) {
split.emitterX(point[4]);
split.spectatorX(point[5]);
shift += 2;
}
if ( split.index().emitterPDF().pdf() &&
!split.index().spectatorPDF().pdf() ) {
split.emitterX(point[4]);
++shift;
}
if ( !split.index().emitterPDF().pdf() &&
split.index().spectatorPDF().pdf() ) {
split.spectatorX(point[4]);
++shift;
}
if ( splittingReweight() )
++shift;
if ( splittingKernel()->nDimAdditional() )
copy(point.begin()+shift,point.end(),split.splittingParameters().begin());
split.hardPt(split.splittingKinematics()->ptMax(split.scale(),
split.emitterX(),
split.spectatorX(),
split.index(),
*splittingKernel()));
}
if ( ! split.splittingKinematics()->generateSplitting(point[0],point[1],point[2],split) ) {
split.lastValue(0.);
return 0.;
}
split.splittingKinematics()->prepareSplitting(split);
if ( split.stoppedEvolving() ) {
split.lastValue(0.);
return 0.;
}
double kernel = splittingKernel()->evaluate(split);
if ( splittingReweight() ) {
if ( !presampling )
kernel *= splittingReweight()->evaluate(split);
else
kernel *= point[shift-1];
}
double jac = split.splittingKinematics()->jacobian();
split.lastValue( abs(jac) * kernel );
if ( kernel < 0. )
return 0.;
return split.lastValue();
}
void DipoleSplittingGenerator::doGenerate() {
assert(!wrapping());
double res = 0.;
Energy startPt = generatedSplitting.hardPt();
while (true) {
try {
res = theExponentialGenerator->generate();
} catch (exsample::exponential_regenerate&) {
generatedSplitting.hardPt(startPt);
continue;
} catch (exsample::hit_and_miss_maxtry&) {
throw DipoleShowerHandler::RedoShower();
} catch (exsample::selection_maxtry&) {
throw DipoleShowerHandler::RedoShower();
}
break;
}
if ( res == 0. ) {
generatedSplitting.lastPt(0.0*GeV);
- generatedSplitting.hardPt(0.0*GeV);
generatedSplitting.didStopEvolving();
} else {
- generatedSplitting.hardPt(generatedSplitting.lastPt());
generatedSplitting.continuesEvolving();
if ( theMCCheck )
theMCCheck->book(generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.scale(),
startPt,
generatedSplitting.lastPt(),
generatedSplitting.lastZ(),
1.);
}
}
Energy DipoleSplittingGenerator::generate(const DipoleSplittingInfo& split) {
fixParameters(split);
if ( wrapping() ) {
return theOtherGenerator->generateWrapped(generatedSplitting);
}
doGenerate();
return generatedSplitting.lastPt();
}
Energy DipoleSplittingGenerator::generateWrapped(DipoleSplittingInfo& split) {
assert(!wrapping());
DipoleSplittingInfo backup = generatedSplitting;
generatedSplitting = split;
fixParameters(split);
try {
doGenerate();
} catch (...) {
split = generatedSplitting;
generatedSplitting = backup;
throw;
}
Energy pt = generatedSplitting.lastPt();
split = generatedSplitting;
generatedSplitting = backup;
return pt;
}
void DipoleSplittingGenerator::completeSplitting(DipoleSplittingInfo& sp) const {
pair<bool,bool> conf = sp.configuration();
sp = generatedSplitting;
sp.configuration(conf);
}
Ptr<DipoleSplittingKernel>::tptr DipoleSplittingGenerator::splittingKernel() const {
if ( wrapping() )
return theOtherGenerator->splittingKernel();
return theSplittingKernel;
}
Ptr<DipoleSplittingReweight>::tptr DipoleSplittingGenerator::splittingReweight() const {
if ( wrapping() )
return theOtherGenerator->splittingReweight();
return theSplittingReweight;
}
Ptr<DipoleSplittingKinematics>::tptr DipoleSplittingGenerator::splittingKinematics() const {
if ( wrapping() )
return theOtherGenerator->splittingKinematics();
return theSplittingKernel->splittingKinematics();
}
void DipoleSplittingGenerator::splittingKernel(Ptr<DipoleSplittingKernel>::tptr sp) {
theSplittingKernel = sp;
if ( theSplittingKernel->mcCheck() )
theMCCheck = theSplittingKernel->mcCheck();
}
void DipoleSplittingGenerator::splittingReweight(Ptr<DipoleSplittingReweight>::tptr sp) {
theSplittingReweight = sp;
}
void DipoleSplittingGenerator::debugGenerator(ostream& os) const {
os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";
os << " generating splittings using\n"
<< " splittingKernel = " << splittingKernel()->name()
<< " splittingKinematics = " << generatedSplitting.splittingKinematics()->name() << "\n"
<< " to sample splittings of type:\n";
os << generatedSplitting;
os << "--------------------------------------------------------------------------------\n";
}
void DipoleSplittingGenerator::debugLastEvent(ostream& os) const {
os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";
os << " last generated event:\n";
os << generatedSplitting;
os << "--------------------------------------------------------------------------------\n";
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSplittingGenerator::persistentOutput(PersistentOStream & os) const {
os << theOtherGenerator << theSplittingKernel << theSplittingReweight << theMCCheck;
}
void DipoleSplittingGenerator::persistentInput(PersistentIStream & is, int) {
is >> theOtherGenerator >> theSplittingKernel >> theSplittingReweight >> theMCCheck;
}
ClassDescription<DipoleSplittingGenerator> DipoleSplittingGenerator::initDipoleSplittingGenerator;
// Definition of the static class description member.
void DipoleSplittingGenerator::Init() {
static ClassDocumentation<DipoleSplittingGenerator> documentation
("DipoleSplittingGenerator is used by the dipole shower "
"to sample splittings from a given dipole splitting kernel.");
static Reference<DipoleSplittingGenerator,DipoleSplittingKernel> interfaceSplittingKernel
("SplittingKernel",
"Set the splitting kernel to sample from.",
&DipoleSplittingGenerator::theSplittingKernel, false, false, true, false, false);
static Reference<DipoleSplittingGenerator,DipoleSplittingReweight> interfaceSplittingReweight
("SplittingReweight",
"Set the splitting reweight.",
&DipoleSplittingGenerator::theSplittingReweight, false, false, true, true, false);
static Reference<DipoleSplittingGenerator,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingGenerator::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
}
diff --git a/DipoleShower/DipoleShowerHandler.cc b/DipoleShower/DipoleShowerHandler.cc
--- a/DipoleShower/DipoleShowerHandler.cc
+++ b/DipoleShower/DipoleShowerHandler.cc
@@ -1,833 +1,851 @@
// -*- C++ -*-
//
// DipoleShowerHandler.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleShowerHandler class.
//
#include "DipoleShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
// include theses to have complete types
#include "Herwig++/Shower/Base/Evolver.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "Herwig++/PDF/MPIPDF.h"
#include "Herwig++/PDF/MinBiasPDF.h"
#include "Herwig++/Shower/Base/ShowerTree.h"
#include "Herwig++/Shower/Base/KinematicsReconstructor.h"
#include "Herwig++/Shower/Base/PartnerFinder.h"
#include "Herwig++/PDF/HwRemDecayer.h"
#include "Herwig++/DipoleShower/Utility/DipolePartonSplitter.h"
using namespace Herwig;
DipoleShowerHandler::DipoleShowerHandler()
: ShowerHandler(), chainOrderVetoScales(true),
nEmissions(0), discardNoEmissions(false),
doFSR(true), doISR(true), realignmentScheme(0),
+ hardFirstEmission(false),
verbosity(0), printEvent(0), nTries(0),
didRadiate(false), didRealign(false),
theFactorizationScaleFactor(1.0),
theRenormalizationScaleFactor(1.0) {}
DipoleShowerHandler::~DipoleShowerHandler() {}
IBPtr DipoleShowerHandler::clone() const {
return new_ptr(*this);
}
IBPtr DipoleShowerHandler::fullclone() const {
return new_ptr(*this);
}
tPPair DipoleShowerHandler::cascade(tSubProPtr sub, XCPtr) {
prepareCascade(sub);
if ( !doFSR && ! doISR )
return sub->incoming();
eventRecord().clear();
eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs());
if ( eventRecord().outgoing().empty() && !doISR )
return sub->incoming();
if ( !eventRecord().incoming().first->coloured() &&
!eventRecord().incoming().second->coloured() &&
!doFSR )
return sub->incoming();
nTries = 0;
while ( true ) {
try {
didRadiate = false;
didRealign = false;
hardScales();
if ( verbosity > 1 ) {
generator()->log() << "DipoleShowerHandler starting off:\n";
eventRecord().debugLastEvent(generator()->log());
generator()->log() << flush;
}
unsigned int nEmitted = 0;
doCascade(nEmitted);
if ( discardNoEmissions ) {
if ( !didRadiate )
throw Veto();
if ( nEmissions )
- if ( nEmissions != nEmitted )
+ if ( nEmissions < nEmitted )
throw Veto();
}
if ( intrinsicPtGenerator ) {
if ( eventRecord().incoming().first->coloured() &&
eventRecord().incoming().second->coloured() ) {
SpinOneLorentzRotation rot =
intrinsicPtGenerator->kick(eventRecord().incoming(),
eventRecord().intermediates());
eventRecord().transform(rot);
}
}
didRealign = realign();
constituentReshuffle();
break;
} catch (RedoShower&) {
if ( ++nTries > maxtry() )
throw ShowerTriesVeto(maxtry());
eventRecord().clear();
eventRecord().prepare(sub,dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr()),pdfs());
continue;
} catch (...) {
throw;
}
}
return eventRecord().fillEventRecord(newStep(),firstInteraction(),didRealign);
}
void DipoleShowerHandler::constituentReshuffle() {
if ( constituentReshuffler ) {
constituentReshuffler->reshuffle(eventRecord().outgoing(),
eventRecord().incoming(),
eventRecord().intermediates());
}
}
void DipoleShowerHandler::hardScales() {
Energy maxPt = generator()->maximumCMEnergy();
- if ( eventRecord().incoming().first->coloured() &&
- eventRecord().incoming().second->coloured() ) {
+ if ( (eventRecord().incoming().first->coloured() ||
+ eventRecord().incoming().second->coloured()) &&
+ !hardFirstEmission ) {
for ( PList::const_iterator p = eventRecord().outgoing().begin();
p != eventRecord().outgoing().end(); ++p )
maxPt = min(maxPt,(**p).momentum().perp());
}
for ( list<DipoleChain>::iterator ch = eventRecord().chains().begin();
ch != eventRecord().chains().end(); ++ch ) {
Energy minVetoScale = -1.*GeV;
for ( list<Dipole>::iterator dip = ch->dipoles().begin();
dip != ch->dipoles().end(); ++dip ) {
// max scale per config
Energy maxFirst = 0.0*GeV;
Energy maxSecond = 0.0*GeV;
for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k =
kernels.begin(); k != kernels.end(); ++k ) {
pair<bool,bool> conf = make_pair(true,false);
if ( (**k).canHandle(dip->index(conf)) ) {
Energy scale =
evolutionOrdering()->hardScale(dip->emitter(conf),dip->spectator(conf),
dip->emitterX(conf),dip->spectatorX(conf),
**k,dip->index(conf));
maxFirst = max(maxFirst,scale);
}
conf = make_pair(false,true);
if ( (**k).canHandle(dip->index(conf)) ) {
Energy scale =
evolutionOrdering()->hardScale(dip->emitter(conf),dip->spectator(conf),
dip->emitterX(conf),dip->spectatorX(conf),
**k,dip->index(conf));
maxSecond = max(maxSecond,scale);
}
}
if ( dip->leftParticle()->vetoScale() >= ZERO ) {
maxFirst = min(maxFirst,sqrt(dip->leftParticle()->vetoScale()));
if ( minVetoScale >= ZERO )
minVetoScale = min(minVetoScale,sqrt(dip->leftParticle()->vetoScale()));
else
minVetoScale = sqrt(dip->leftParticle()->vetoScale());
}
if ( dip->rightParticle()->vetoScale() >= ZERO ) {
maxSecond = min(maxSecond,sqrt(dip->rightParticle()->vetoScale()));
if ( minVetoScale >= ZERO )
minVetoScale = min(minVetoScale,sqrt(dip->rightParticle()->vetoScale()));
else
minVetoScale = sqrt(dip->rightParticle()->vetoScale());
}
maxFirst = min(maxPt,maxFirst);
dip->emitterScale(make_pair(true,false),maxFirst);
maxSecond = min(maxPt,maxSecond);
dip->emitterScale(make_pair(false,true),maxSecond);
}
if ( !evolutionOrdering()->independentDipoles() &&
chainOrderVetoScales &&
minVetoScale >= ZERO ) {
for ( list<Dipole>::iterator dip = ch->dipoles().begin();
dip != ch->dipoles().end(); ++dip ) {
dip->leftScale(min(dip->leftScale(),minVetoScale));
dip->rightScale(min(dip->rightScale(),minVetoScale));
}
}
}
}
Energy DipoleShowerHandler::getWinner(DipoleSplittingInfo& winner,
const Dipole& dip,
pair<bool,bool> conf) {
if ( !dip.index(conf).initialStateEmitter() &&
!doFSR ) {
winner.didStopEvolving();
return 0.0*GeV;
}
if ( dip.index(conf).initialStateEmitter() &&
!doISR ) {
winner.didStopEvolving();
return 0.0*GeV;
}
DipoleSplittingInfo candidate;
candidate.index(dip.index(conf));
candidate.configuration(conf);
candidate.emitterX(dip.emitterX(conf));
candidate.spectatorX(dip.spectatorX(conf));
if ( generators().find(candidate.index()) == generators().end() )
getGenerators(candidate.index());
//
// NOTE -- needs proper fixing at some point
//
// For some very strange reason, equal_range gives back
// key ranges it hasn't been asked for. This particularly
// happens e.g. for FI dipoles of the same kind, but different
// PDF (hard vs MPI PDF). I can't see a reason for this,
// as DipoleIndex properly implements comparison for equality
// and (lexicographic) ordering; for the time being, we
// use equal_range, extented by an explicit check for wether
// the key is indeed what we wanted. See line after (*) comment
// below.
//
pair<GeneratorMap::iterator,GeneratorMap::iterator> gens
= generators().equal_range(candidate.index());
tPPtr emitter = dip.emitter(conf);
tPPtr spectator = dip.spectator(conf);
Energy startScale = dip.emitterScale(conf);
Energy winnerScale = 0.0*GeV;
GeneratorMap::iterator winnerGen = generators().end();
for ( GeneratorMap::iterator gen = gens.first; gen != gens.second; ++gen ) {
// (*) see NOTE above
if ( !(gen->first == candidate.index()) )
continue;
if ( startScale <= gen->second->splittingKinematics()->IRCutoff() )
continue;
Energy dScale =
gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),
spectator->momentum());
// in very exceptional cases happening in DIS
if ( isnan(dScale/GeV ) )
throw RedoShower();
candidate.scale(dScale);
candidate.continuesEvolving();
candidate.hardPt(evolutionOrdering()->maxPt(startScale,candidate,*(gen->second->splittingKernel())));
gen->second->generate(candidate);
Energy nextScale = evolutionOrdering()->evolutionScale(gen->second->lastSplitting(),*(gen->second->splittingKernel()));
if ( nextScale > winnerScale ) {
winner = candidate;
gen->second->completeSplitting(winner);
winnerGen = gen;
winnerScale = nextScale;
}
}
if ( winnerGen == generators().end() ) {
winner.didStopEvolving();
return 0.0*GeV;
}
if ( winner.stoppedEvolving() )
return 0.0*GeV;
return winnerScale;
}
void DipoleShowerHandler::doCascade(unsigned int& emDone) {
if ( nEmissions )
if ( emDone == nEmissions )
return;
DipoleSplittingInfo winner;
DipoleSplittingInfo dipoleWinner;
while ( eventRecord().haveChain() ) {
if ( verbosity > 2 ) {
generator()->log() << "DipoleShowerHandler selecting splittings for the chain:\n"
<< eventRecord().currentChain() << flush;
}
list<Dipole>::iterator winnerDip = eventRecord().currentChain().dipoles().end();
Energy winnerScale = 0.0*GeV;
Energy nextLeftScale = 0.0*GeV;
Energy nextRightScale = 0.0*GeV;
for ( list<Dipole>::iterator dip = eventRecord().currentChain().dipoles().begin();
dip != eventRecord().currentChain().dipoles().end(); ++dip ) {
nextLeftScale = getWinner(dipoleWinner,*dip,make_pair(true,false));
if ( nextLeftScale > winnerScale ) {
winnerScale = nextLeftScale;
winner = dipoleWinner;
winnerDip = dip;
}
nextRightScale = getWinner(dipoleWinner,*dip,make_pair(false,true));
if ( nextRightScale > winnerScale ) {
winnerScale = nextRightScale;
winner = dipoleWinner;
winnerDip = dip;
}
if ( evolutionOrdering()->independentDipoles() ) {
Energy dipScale = max(nextLeftScale,nextRightScale);
if ( dip->leftScale() > dipScale )
dip->leftScale(dipScale);
if ( dip->rightScale() > dipScale )
dip->rightScale(dipScale);
}
}
if ( verbosity > 1 ) {
if ( winnerDip != eventRecord().currentChain().dipoles().end() )
generator()->log() << "DipoleShowerHandler selected the splitting:\n"
<< winner << " for the dipole\n"
<< (*winnerDip) << flush;
else
generator()->log() << "DipoleShowerHandler could not select a splitting above the IR cutoff\n"
<< flush;
}
// pop the chain if no dipole did radiate
if ( winnerDip == eventRecord().currentChain().dipoles().end() ) {
eventRecord().popChain();
continue;
}
// otherwise perform the splitting
didRadiate = true;
pair<list<Dipole>::iterator,list<Dipole>::iterator> children;
DipoleChain* firstChain = 0;
DipoleChain* secondChain = 0;
eventRecord().split(winnerDip,winner,children,firstChain,secondChain);
assert(firstChain && secondChain);
evolutionOrdering()->setEvolutionScale(winnerScale,winner,*firstChain,children);
if ( !secondChain->dipoles().empty() )
evolutionOrdering()->setEvolutionScale(winnerScale,winner,*secondChain,children);
if ( verbosity > 1 ) {
generator()->log() << "DipoleShowerHandler did split the last selected dipole into:\n"
<< (*children.first) << (*children.second) << flush;
}
if ( verbosity > 2 ) {
generator()->log() << "After splitting the last selected dipole, "
<< "DipoleShowerHandler encountered the following chains:\n"
<< (*firstChain) << (*secondChain) << flush;
}
if ( nEmissions )
if ( ++emDone == nEmissions )
return;
}
}
bool DipoleShowerHandler::realign() {
if ( !didRadiate && !intrinsicPtGenerator )
return false;
if ( eventRecord().incoming().first->coloured() ||
eventRecord().incoming().second->coloured() ) {
if ( eventRecord().incoming().first->momentum().perp2()/GeV2 < 1e-10 &&
eventRecord().incoming().second->momentum().perp2()/GeV2 < 1e-10 )
return false;
pair<Lorentz5Momentum,Lorentz5Momentum> inMomenta
(eventRecord().incoming().first->momentum(),
eventRecord().incoming().second->momentum());
SpinOneLorentzRotation transform((inMomenta.first+inMomenta.second).findBoostToCM());
Axis dir = (transform * inMomenta.first).vect().unit();
Axis rot (-dir.y(),dir.x(),0);
double theta = dir.theta();
if ( lastParticles().first->momentum().z() < ZERO )
theta = -theta;
transform.rotate(-theta,rot);
inMomenta.first = transform*inMomenta.first;
inMomenta.second = transform*inMomenta.second;
assert(inMomenta.first.z() > ZERO &&
inMomenta.second.z() < ZERO);
Energy2 sHat =
(eventRecord().incoming().first->momentum() +
eventRecord().incoming().second->momentum()).m2();
+ Energy q = inMomenta.first.z();
pair<Energy,Energy> masses(eventRecord().incoming().first->mass(),
eventRecord().incoming().second->mass());
pair<Energy,Energy> qs;
if ( !eventRecord().incoming().first->coloured() ) {
assert(masses.second == ZERO);
qs.first = eventRecord().incoming().first->momentum().z();
qs.second = (sHat-sqr(masses.first))/(2.*(qs.first+sqrt(sqr(masses.first)+sqr(qs.first))));
} else if ( !eventRecord().incoming().second->coloured() ) {
assert(masses.first == ZERO);
qs.second = eventRecord().incoming().second->momentum().z();
qs.first = (sHat-sqr(masses.second))/(2.*(qs.second+sqrt(sqr(masses.second)+sqr(qs.second))));
} else {
assert(masses.first == ZERO && masses.second == ZERO);
if ( realignmentScheme == 0 ) {
double yX = eventRecord().pX().rapidity();
double yInt = (transform*eventRecord().pX()).rapidity();
double dy = yX-yInt;
qs.first = (sqrt(sHat)/2.)*exp(dy);
qs.second = (sqrt(sHat)/2.)*exp(-dy);
} else if ( realignmentScheme == 1 ) {
Energy sS = sqrt((lastParticles().first->momentum() +
lastParticles().second->momentum()).m2());
qs.first = eventRecord().fractions().first * sS / 2.;
qs.second = eventRecord().fractions().second * sS / 2.;
}
}
double beta =
(qs.first-qs.second) /
( sqrt(sqr(masses.first)+sqr(qs.first)) +
sqrt(sqr(masses.second)+sqr(qs.second)) );
transform.boostZ(beta);
Lorentz5Momentum tmp;
if ( eventRecord().incoming().first->coloured() ) {
tmp = eventRecord().incoming().first->momentum();
tmp = transform * tmp;
eventRecord().incoming().first->set5Momentum(tmp);
}
if ( eventRecord().incoming().second->coloured() ) {
tmp = eventRecord().incoming().second->momentum();
tmp = transform * tmp;
eventRecord().incoming().second->set5Momentum(tmp);
}
eventRecord().transform(transform);
return true;
}
return false;
}
void DipoleShowerHandler::resetAlphaS(Ptr<AlphaSBase>::tptr as) {
for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k = kernels.begin();
k != kernels.end(); ++k )
(**k).alphaS(as);
// clear the generators to be rebuild
// actually, there shouldn't be any generators
// when this happens.
generators().clear();
}
void DipoleShowerHandler::resetReweight(Ptr<DipoleSplittingReweight>::tptr rw) {
for ( GeneratorMap::iterator k = generators().begin();
k != generators().end(); ++k )
k->second->splittingReweight(rw);
}
void DipoleShowerHandler::getGenerators(const DipoleIndex& ind,
Ptr<DipoleSplittingReweight>::tptr rw) {
bool gotone = false;
for ( vector<Ptr<DipoleSplittingKernel>::ptr>::iterator k =
kernels.begin(); k != kernels.end(); ++k ) {
if ( (**k).canHandle(ind) ) {
if ( verbosity > 0 ) {
generator()->log() << "DipoleShowerHandler encountered the dipole configuration\n"
<< ind << " in event number "
<< eventHandler()->currentEvent()->number()
<< "\nwhich can be handled by the splitting kernel '"
<< (**k).name() << "'.\n" << flush;
}
gotone = true;
Ptr<DipoleSplittingGenerator>::ptr nGenerator =
new_ptr(DipoleSplittingGenerator());
nGenerator->splittingKernel(*k);
nGenerator->splittingKernel()->renormalizationScaleFactor(theRenormalizationScaleFactor);
nGenerator->splittingKernel()->factorizationScaleFactor(theFactorizationScaleFactor);
GeneratorMap::const_iterator equivalent = generators().end();
for ( GeneratorMap::const_iterator eq = generators().begin();
eq != generators().end(); ++eq ) {
if ( !eq->second->wrapping() )
if ( (**k).canHandleEquivalent(ind,*(eq->second->splittingKernel()),eq->first) ) {
equivalent = eq;
if ( verbosity > 0 ) {
generator()->log() << "The dipole configuration "
<< ind
<< " can equivalently be handled by the existing\n"
<< "generator for configuration "
<< eq->first << " using the kernel '"
<< eq->second->splittingKernel()->name()
<< "'\n" << flush;
}
break;
}
}
if ( equivalent != generators().end() ) {
nGenerator->wrap(equivalent->second);
}
DipoleSplittingInfo dummy;
dummy.index(ind);
nGenerator->splittingReweight(rw);
nGenerator->prepare(dummy);
generators().insert(make_pair(ind,nGenerator));
}
}
if ( !gotone ) {
generator()->logWarning(Exception()
<< "DipoleShowerHandler could not "
<< "find a splitting kernel which is able "
<< "to handle splittings off the dipole "
<< ind << ".\n"
<< "Please check the input files."
<< Exception::warning);
}
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleShowerHandler::doinit() {
ShowerHandler::doinit();
if ( theGlobalAlphaS )
resetAlphaS(theGlobalAlphaS);
}
void DipoleShowerHandler::dofinish() {
ShowerHandler::dofinish();
}
void DipoleShowerHandler::doinitrun() {
ShowerHandler::doinitrun();
}
void DipoleShowerHandler::persistentOutput(PersistentOStream & os) const {
os << kernels << theEvolutionOrdering
<< constituentReshuffler << intrinsicPtGenerator
<< theGlobalAlphaS << chainOrderVetoScales
<< nEmissions << discardNoEmissions << doFSR << doISR
- << realignmentScheme << verbosity << printEvent
+ << realignmentScheme << hardFirstEmission << verbosity << printEvent
<< theFactorizationScaleFactor << theRenormalizationScaleFactor;
}
void DipoleShowerHandler::persistentInput(PersistentIStream & is, int) {
is >> kernels >> theEvolutionOrdering
>> constituentReshuffler >> intrinsicPtGenerator
>> theGlobalAlphaS >> chainOrderVetoScales
>> nEmissions >> discardNoEmissions >> doFSR >> doISR
- >> realignmentScheme >> verbosity >> printEvent
+ >> realignmentScheme >> hardFirstEmission >> verbosity >> printEvent
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor;
}
ClassDescription<DipoleShowerHandler> DipoleShowerHandler::initDipoleShowerHandler;
// Definition of the static class description member.
void DipoleShowerHandler::Init() {
static ClassDocumentation<DipoleShowerHandler> documentation
("The DipoleShowerHandler class manages the showering using "
"the dipole shower algorithm.",
"The shower evolution was performed using the algorithm described in "
"\\cite{Platzer:2009jq} and \\cite{Platzer:2011bc}.",
"%\\cite{Platzer:2009jq}\n"
"\\bibitem{Platzer:2009jq}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Coherent Parton Showers with Local Recoils,''\n"
" JHEP {\\bf 1101}, 024 (2011)\n"
"arXiv:0909.5593 [hep-ph].\n"
"%%CITATION = ARXIV:0909.5593;%%\n"
"%\\cite{Platzer:2011bc}\n"
"\\bibitem{Platzer:2011bc}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Dipole Showers and Automated NLO Matching in Herwig++,''\n"
"arXiv:1109.6256 [hep-ph].\n"
"%%CITATION = ARXIV:1109.6256;%%");
static RefVector<DipoleShowerHandler,DipoleSplittingKernel> interfaceKernels
("Kernels",
"Set the splitting kernels to be used by the dipole shower.",
&DipoleShowerHandler::kernels, -1, false, false, true, false, false);
static Reference<DipoleShowerHandler,DipoleEvolutionOrdering> interfaceEvolutionOrdering
("EvolutionOrdering",
"Set the evolution ordering to be used.",
&DipoleShowerHandler::theEvolutionOrdering, false, false, true, false, false);
static Reference<DipoleShowerHandler,ConstituentReshuffler> interfaceConstituentReshuffler
("ConstituentReshuffler",
"The object to be used to reshuffle partons to their constitutent mass shells.",
&DipoleShowerHandler::constituentReshuffler, false, false, true, true, false);
static Reference<DipoleShowerHandler,IntrinsicPtGenerator> interfaceIntrinsicPtGenerator
("IntrinsicPtGenerator",
"Set the object in charge to generate intrinsic pt for incoming partons.",
&DipoleShowerHandler::intrinsicPtGenerator, false, false, true, true, false);
static Reference<DipoleShowerHandler,AlphaSBase> interfaceGlobalAlphaS
("GlobalAlphaS",
"Set a global strong coupling for all splitting kernels.",
&DipoleShowerHandler::theGlobalAlphaS, false, false, true, true, false);
static Switch<DipoleShowerHandler,bool> interfaceDoFSR
("DoFSR",
"Switch on or off final state radiation.",
&DipoleShowerHandler::doFSR, true, false, false);
static SwitchOption interfaceDoFSROn
(interfaceDoFSR,
"On",
"Switch on final state radiation.",
true);
static SwitchOption interfaceDoFSROff
(interfaceDoFSR,
"Off",
"Switch off final state radiation.",
false);
static Switch<DipoleShowerHandler,bool> interfaceDoISR
("DoISR",
"Switch on or off initial state radiation.",
&DipoleShowerHandler::doISR, true, false, false);
static SwitchOption interfaceDoISROn
(interfaceDoISR,
"On",
"Switch on initial state radiation.",
true);
static SwitchOption interfaceDoISROff
(interfaceDoISR,
"Off",
"Switch off initial state radiation.",
false);
+ static Switch<DipoleShowerHandler,bool> interfaceHardFirstEmission
+ ("HardFirstEmission",
+ "Switch on or off hard first emission.",
+ &DipoleShowerHandler::hardFirstEmission, false, false, false);
+ static SwitchOption interfaceHardFirstEmissionOn
+ (interfaceHardFirstEmission,
+ "On",
+ "Switch on hard first emission.",
+ true);
+ static SwitchOption interfaceHardFirstEmissionOff
+ (interfaceHardFirstEmission,
+ "Off",
+ "Switch off hard first emission.",
+ false);
+
static Switch<DipoleShowerHandler,int> interfaceRealignmentScheme
("RealignmentScheme",
"The realignment scheme to use.",
&DipoleShowerHandler::realignmentScheme, 0, false, false);
static SwitchOption interfaceRealignmentSchemePreserveRapidity
(interfaceRealignmentScheme,
"PreserveRapidity",
"Preserve the rapidity of non-coloured outgoing system.",
0);
static SwitchOption interfaceRealignmentSchemeEvolutionFractions
(interfaceRealignmentScheme,
"EvolutionFractions",
"Use momentum fractions as generated by the evolution.",
1);
static SwitchOption interfaceRealignmentSchemeCollisionFrame
(interfaceRealignmentScheme,
"CollisionFrame",
"Determine realignment from collision frame.",
2);
static Switch<DipoleShowerHandler,bool> interfaceChainOrderVetoScales
("ChainOrderVetoScales",
"[experimental] Switch on or off the chain ordering for veto scales.",
&DipoleShowerHandler::chainOrderVetoScales, true, false, false);
static SwitchOption interfaceChainOrderVetoScalesOn
(interfaceChainOrderVetoScales,
"On",
"Switch on chain ordering for veto scales.",
true);
static SwitchOption interfaceChainOrderVetoScalesOff
(interfaceChainOrderVetoScales,
"Off",
"Switch off chain ordering for veto scales.",
false);
interfaceChainOrderVetoScales.rank(-1);
static Parameter<DipoleShowerHandler,unsigned int> interfaceNEmissions
("NEmissions",
"[debug option] Limit the number of emissions to be generated. Zero does not limit the number of emissions.",
&DipoleShowerHandler::nEmissions, 0, 0, 0,
false, false, Interface::lowerlim);
interfaceNEmissions.rank(-1);
static Switch<DipoleShowerHandler,bool> interfaceDiscardNoEmissions
("DiscardNoEmissions",
"[debug option] Discard events without radiation.",
&DipoleShowerHandler::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);
interfaceDiscardNoEmissions.rank(-1);
static Parameter<DipoleShowerHandler,int> interfaceVerbosity
("Verbosity",
"[debug option] Set the level of debug information provided.",
&DipoleShowerHandler::verbosity, 0, 0, 0,
false, false, Interface::lowerlim);
interfaceVerbosity.rank(-1);
static Parameter<DipoleShowerHandler,int> interfacePrintEvent
("PrintEvent",
"[debug option] The number of events for which debugging information should be provided.",
&DipoleShowerHandler::printEvent, 0, 0, 0,
false, false, Interface::lowerlim);
interfacePrintEvent.rank(-1);
static Parameter<DipoleShowerHandler,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&DipoleShowerHandler::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleShowerHandler,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&DipoleShowerHandler::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
}
diff --git a/DipoleShower/DipoleShowerHandler.h b/DipoleShower/DipoleShowerHandler.h
--- a/DipoleShower/DipoleShowerHandler.h
+++ b/DipoleShower/DipoleShowerHandler.h
@@ -1,399 +1,404 @@
// -*- C++ -*-
//
// DipoleShowerHandler.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleShowerHandler_H
#define HERWIG_DipoleShowerHandler_H
//
// This is the declaration of the DipoleShowerHandler class.
//
#include "Herwig++/Shower/ShowerHandler.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingReweight.h"
#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleEventRecord.h"
#include "Herwig++/DipoleShower/Base/DipoleEvolutionOrdering.h"
#include "Herwig++/DipoleShower/Utility/ConstituentReshuffler.h"
#include "Herwig++/DipoleShower/Utility/IntrinsicPtGenerator.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief The DipoleShowerHandler class manages the showering using
* the dipole shower algorithm.
*
* @see \ref DipoleShowerHandlerInterfaces "The interfaces"
* defined for DipoleShowerHandler.
*/
class DipoleShowerHandler: public ShowerHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleShowerHandler();
/**
* The destructor.
*/
virtual ~DipoleShowerHandler();
//@}
public:
/**
* Indicate a problem in the shower.
*/
struct RedoShower {};
/**
* Insert an additional splitting kernel.
*/
void addSplitting(Ptr<DipoleSplittingKernel>::ptr sp) {
kernels.push_back(sp);
}
/**
* Reset the alpha_s for all splitting kernels.
*/
void resetAlphaS(Ptr<AlphaSBase>::tptr);
/**
* Reset the splitting reweight for all splitting kernels.
*/
void resetReweight(Ptr<DipoleSplittingReweight>::tptr);
protected:
typedef multimap<DipoleIndex,Ptr<DipoleSplittingGenerator>::ptr> GeneratorMap;
/**
* The main method which manages the showering of a subprocess.
*/
virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb);
/**
* Build splitting generators for the given
* dipole index.
*/
void getGenerators(const DipoleIndex&,
Ptr<DipoleSplittingReweight>::tptr rw =
Ptr<DipoleSplittingReweight>::tptr());
/**
* Setup the hard scales.
*/
void hardScales();
/**
* Return the evolution ordering
*/
Ptr<DipoleEvolutionOrdering>::tptr evolutionOrdering() const { return theEvolutionOrdering; }
/**
* Reshuffle to constituent mass shells
*/
void constituentReshuffle();
/**
* Access the generator map
*/
GeneratorMap& generators() { return theGenerators; }
/**
* Access the event record
*/
DipoleEventRecord& eventRecord() { return theEventRecord; }
/**
* Return the event record
*/
const DipoleEventRecord& eventRecord() const { return theEventRecord; }
/**
* Return the splitting kernels.
*/
const vector<Ptr<DipoleSplittingKernel>::ptr>& splittingKernels() const {
return kernels;
}
/**
* Realign the event such as to have the incoming partons along thre
* beam axes.
*/
bool realign();
private:
/**
* Perform the cascade.
*/
void doCascade(unsigned int& emDone);
/**
* Get the winning splitting for the
* given dipole and configuration.
*/
Energy getWinner(DipoleSplittingInfo& winner,
const Dipole& dip,
pair<bool,bool> conf);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The splitting kernels to be used.
*/
vector<Ptr<DipoleSplittingKernel>::ptr> kernels;
/**
* The evolution ordering considered
*/
Ptr<DipoleEvolutionOrdering>::ptr theEvolutionOrdering;
/**
* The ConstituentReshuffler to be used
*/
Ptr<ConstituentReshuffler>::ptr constituentReshuffler;
/**
* The intrinsic pt generator to be used.
*/
Ptr<IntrinsicPtGenerator>::ptr intrinsicPtGenerator;
/**
* A global alpha_s to be used for all splitting kernels.
*/
Ptr<AlphaSBase>::ptr theGlobalAlphaS;
/**
* Apply chain ordering to events from matrix
* element corrections.
*/
bool chainOrderVetoScales;
/**
* Limit the number of emissions.
* Limit applied if > 0.
*/
unsigned int nEmissions;
/**
* Discard events which did not radiate.
*/
bool discardNoEmissions;
/**
* Switch on or off final state radiation.
*/
bool doFSR;
/**
* Switch on or off initial state radiation.
*/
bool doISR;
/**
* The realignment scheme
*/
int realignmentScheme;
+ /**
+ * True, if first emission should use the available phase space
+ */
+ bool hardFirstEmission;
+
private:
/**
* The verbosity level.
* 0 - print no info
* 1 - print diagnostic information on setting up
* splitting generators etc.
* 2 - print detailed event information for up to
* printEvent events.
* 3 - print dipole chains after each splitting.
*/
int verbosity;
/**
* See verbosity.
*/
int printEvent;
private:
/**
* The splitting generators indexed by the dipole
* indices they can work on.
*/
GeneratorMap theGenerators;
/**
* The evnt record used.
*/
DipoleEventRecord theEventRecord;
/**
* The number of shoer tries so far.
*/
unsigned int nTries;
/**
* Whether or not we did radiate anything
*/
bool didRadiate;
/**
* Whether or not we did realign the event
*/
bool didRealign;
private:
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<DipoleShowerHandler> initDipoleShowerHandler;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleShowerHandler & operator=(const DipoleShowerHandler &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleShowerHandler. */
template <>
struct BaseClassTrait<Herwig::DipoleShowerHandler,1> {
/** Typedef of the first base class of DipoleShowerHandler. */
typedef Herwig::ShowerHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleShowerHandler class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleShowerHandler>
: public ClassTraitsBase<Herwig::DipoleShowerHandler> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleShowerHandler"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleShowerHandler is implemented. It may also include several, space-separated,
* libraries if the class DipoleShowerHandler depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleShowerHandler_H */
diff --git a/DipoleShower/Kernels/FFgx2qqxDipoleKernel.cc b/DipoleShower/Kernels/FFgx2qqxDipoleKernel.cc
--- a/DipoleShower/Kernels/FFgx2qqxDipoleKernel.cc
+++ b/DipoleShower/Kernels/FFgx2qqxDipoleKernel.cc
@@ -1,100 +1,101 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FFgx2qqxDipoleKernel class.
//
#include "FFgx2qqxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FFgx2qqxDipoleKernel::FFgx2qqxDipoleKernel()
: DipoleSplittingKernel() {}
FFgx2qqxDipoleKernel::~FFgx2qqxDipoleKernel() {}
IBPtr FFgx2qqxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FFgx2qqxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FFgx2qqxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
+ flavour()->mass() == ZERO &&
!ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool FFgx2qqxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() + sk.emission(b)->id() == 0 &&
abs(sk.emitter(b)->id()) < 6 &&
sk.emitter(b)->mass() == ZERO;
}
tcPDPtr FFgx2qqxDipoleKernel::emitter(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour();
}
tcPDPtr FFgx2qqxDipoleKernel::emission(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour()->CC();
}
tcPDPtr FFgx2qqxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FFgx2qqxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
- ret *= .5 * ( 1. - 2.*z*(1.-z) );
+ ret *= .25 * ( 1. - 2.*z*(1.-z) );
return ret;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFgx2qqxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void FFgx2qqxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FFgx2qqxDipoleKernel> FFgx2qqxDipoleKernel::initFFgx2qqxDipoleKernel;
// Definition of the static class description member.
void FFgx2qqxDipoleKernel::Init() {
static ClassDocumentation<FFgx2qqxDipoleKernel> documentation
("FFgx2qqxDipoleKernel");
}
diff --git a/DipoleShower/Kernels/FIgx2qqxDipoleKernel.cc b/DipoleShower/Kernels/FIgx2qqxDipoleKernel.cc
--- a/DipoleShower/Kernels/FIgx2qqxDipoleKernel.cc
+++ b/DipoleShower/Kernels/FIgx2qqxDipoleKernel.cc
@@ -1,101 +1,102 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FIgx2qqxDipoleKernel class.
//
#include "FIgx2qqxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FIgx2qqxDipoleKernel::FIgx2qqxDipoleKernel()
: DipoleSplittingKernel() {}
FIgx2qqxDipoleKernel::~FIgx2qqxDipoleKernel() {}
IBPtr FIgx2qqxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr FIgx2qqxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool FIgx2qqxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
+ flavour()->mass() == ZERO &&
!ind.initialStateEmitter() && ind.initialStateSpectator();
}
bool FIgx2qqxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() + sk.emission(b)->id() == 0 &&
abs(sk.emitter(b)->id()) < 6 &&
- sk.emitter(b)->mass() == ZERO &&
+ // sk.emitter(b)->mass() == ZERO &&
a.spectatorPDF() == b.spectatorPDF();
}
tcPDPtr FIgx2qqxDipoleKernel::emitter(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour();
}
tcPDPtr FIgx2qqxDipoleKernel::emission(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour()->CC();
}
tcPDPtr FIgx2qqxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double FIgx2qqxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
- ret *= .5 * (1.-z*(1.-z));
+ ret *= .25 * (1.-2.*z*(1.-z));
return ret;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIgx2qqxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void FIgx2qqxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FIgx2qqxDipoleKernel> FIgx2qqxDipoleKernel::initFIgx2qqxDipoleKernel;
// Definition of the static class description member.
void FIgx2qqxDipoleKernel::Init() {
static ClassDocumentation<FIgx2qqxDipoleKernel> documentation
("FIgx2qqxDipoleKernel");
}
diff --git a/DipoleShower/Kernels/IFgx2ggxDipoleKernel.cc b/DipoleShower/Kernels/IFgx2ggxDipoleKernel.cc
--- a/DipoleShower/Kernels/IFgx2ggxDipoleKernel.cc
+++ b/DipoleShower/Kernels/IFgx2ggxDipoleKernel.cc
@@ -1,101 +1,101 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IFgx2ggxDipoleKernel class.
//
#include "IFgx2ggxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFgx2ggxDipoleKernel::IFgx2ggxDipoleKernel()
: DipoleSplittingKernel() {}
IFgx2ggxDipoleKernel::~IFgx2ggxDipoleKernel() {}
IBPtr IFgx2ggxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr IFgx2ggxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool IFgx2ggxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool IFgx2ggxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() == ParticleID::g &&
sk.emission(b)->id() == ParticleID::g &&
a.emitterPDF() == b.emitterPDF();
}
tcPDPtr IFgx2ggxDipoleKernel::emitter(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr IFgx2ggxDipoleKernel::emission(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr IFgx2ggxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double IFgx2ggxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
double ratio = sqr(split.lastPt()/split.scale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z - ratio );
double u = ratio / (1.-z);
//double rho = (1.+x+u)/x;
//ret *= 6. * ( 1./(1.-x+u) + (1.-x)/x - 1. + x*(1.-x) +u/rho );
- ret *= 6. * ( 1./(1.-x+u) + (1.-x)/x - 1. + x*(1.-x) );
+ ret *= 3. * ( 1./(1.-x+u) + (1.-x)/x - 1. + x*(1.-x) );
return ret;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFgx2ggxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void IFgx2ggxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<IFgx2ggxDipoleKernel> IFgx2ggxDipoleKernel::initIFgx2ggxDipoleKernel;
// Definition of the static class description member.
void IFgx2ggxDipoleKernel::Init() {
static ClassDocumentation<IFgx2ggxDipoleKernel> documentation
("IFgx2ggxDipoleKernel");
}
diff --git a/DipoleShower/Kernels/IFgx2qqxDipoleKernel.cc b/DipoleShower/Kernels/IFgx2qqxDipoleKernel.cc
--- a/DipoleShower/Kernels/IFgx2qqxDipoleKernel.cc
+++ b/DipoleShower/Kernels/IFgx2qqxDipoleKernel.cc
@@ -1,103 +1,103 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IFgx2qqxDipoleKernel class.
//
#include "IFgx2qqxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFgx2qqxDipoleKernel::IFgx2qqxDipoleKernel()
: DipoleSplittingKernel() {}
IFgx2qqxDipoleKernel::~IFgx2qqxDipoleKernel() {}
IBPtr IFgx2qqxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr IFgx2qqxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool IFgx2qqxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
ind.initialStateEmitter() && !ind.initialStateSpectator();
}
bool IFgx2qqxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
flavour() == sk.flavour() &&
abs(flavour()->id()) < 6 &&
flavour()->mass() == ZERO &&
a.emitterPDF() == b.emitterPDF();
}
tcPDPtr IFgx2qqxDipoleKernel::emitter(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour();
}
tcPDPtr IFgx2qqxDipoleKernel::emission(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour();
}
tcPDPtr IFgx2qqxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double IFgx2qqxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
double ratio = sqr(split.lastPt()/split.scale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z - ratio );
- ret *= (!strictLargeN() ? 4./3. : 3./2.) * ( x + 2.*(1.-x)/x );
+ ret *= 0.5 * (!strictLargeN() ? 4./3. : 3./2.) * ( x + 2.*(1.-x)/x );
return ret;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFgx2qqxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void IFgx2qqxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<IFgx2qqxDipoleKernel> IFgx2qqxDipoleKernel::initIFgx2qqxDipoleKernel;
// Definition of the static class description member.
void IFgx2qqxDipoleKernel::Init() {
static ClassDocumentation<IFgx2qqxDipoleKernel> documentation
("IFgx2qqxDipoleKernel");
}
diff --git a/DipoleShower/Kernels/IIgx2ggxDipoleKernel.cc b/DipoleShower/Kernels/IIgx2ggxDipoleKernel.cc
--- a/DipoleShower/Kernels/IIgx2ggxDipoleKernel.cc
+++ b/DipoleShower/Kernels/IIgx2ggxDipoleKernel.cc
@@ -1,100 +1,100 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IIgx2ggxDipoleKernel class.
//
#include "IIgx2ggxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IIgx2ggxDipoleKernel::IIgx2ggxDipoleKernel()
: DipoleSplittingKernel() {}
IIgx2ggxDipoleKernel::~IIgx2ggxDipoleKernel() {}
IBPtr IIgx2ggxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr IIgx2ggxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool IIgx2ggxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
ind.initialStateEmitter() && ind.initialStateSpectator();
}
bool IIgx2ggxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
sk.emitter(b)->id() == ParticleID::g &&
sk.emission(b)->id() == ParticleID::g &&
a.emitterPDF() == b.emitterPDF() &&
a.spectatorData() == b.spectatorData() &&
a.spectatorPDF() == b.spectatorPDF();
}
tcPDPtr IIgx2ggxDipoleKernel::emitter(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr IIgx2ggxDipoleKernel::emission(const DipoleIndex&) const {
return getParticleData(ParticleID::g);
}
tcPDPtr IIgx2ggxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double IIgx2ggxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
double ratio = sqr(split.lastPt()/split.scale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z );
- ret *= 6. * ( x/(1.-x) + (1.-x)/x +x*(1.-x) );
+ ret *= 3. * ( x/(1.-x) + (1.-x)/x +x*(1.-x) );
return ret;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IIgx2ggxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void IIgx2ggxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<IIgx2ggxDipoleKernel> IIgx2ggxDipoleKernel::initIIgx2ggxDipoleKernel;
// Definition of the static class description member.
void IIgx2ggxDipoleKernel::Init() {
static ClassDocumentation<IIgx2ggxDipoleKernel> documentation
("IIgx2ggxDipoleKernel");
}
diff --git a/DipoleShower/Kernels/IIgx2qqxDipoleKernel.cc b/DipoleShower/Kernels/IIgx2qqxDipoleKernel.cc
--- a/DipoleShower/Kernels/IIgx2qqxDipoleKernel.cc
+++ b/DipoleShower/Kernels/IIgx2qqxDipoleKernel.cc
@@ -1,105 +1,105 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IIgx2qqxDipoleKernel class.
//
#include "IIgx2qqxDipoleKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IIgx2qqxDipoleKernel::IIgx2qqxDipoleKernel()
: DipoleSplittingKernel() {}
IIgx2qqxDipoleKernel::~IIgx2qqxDipoleKernel() {}
IBPtr IIgx2qqxDipoleKernel::clone() const {
return new_ptr(*this);
}
IBPtr IIgx2qqxDipoleKernel::fullclone() const {
return new_ptr(*this);
}
bool IIgx2qqxDipoleKernel::canHandle(const DipoleIndex& ind) const {
return
ind.emitterData()->id() == ParticleID::g &&
ind.spectatorData()->mass() == ZERO &&
ind.initialStateEmitter() && ind.initialStateSpectator();
}
bool IIgx2qqxDipoleKernel::canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const {
assert(canHandle(a));
if ( !canHandle(b) )
return false;
return
flavour() == sk.flavour() &&
abs(flavour()->id()) < 6 &&
flavour()->mass() == ZERO &&
a.emitterPDF() == b.emitterPDF() &&
a.spectatorData() == b.spectatorData() &&
a.spectatorPDF() == b.spectatorPDF();
}
tcPDPtr IIgx2qqxDipoleKernel::emitter(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour();
}
tcPDPtr IIgx2qqxDipoleKernel::emission(const DipoleIndex&) const {
assert(flavour());
assert(abs(flavour()->id()) < 6 && flavour()->mass() == ZERO);
return flavour();
}
tcPDPtr IIgx2qqxDipoleKernel::spectator(const DipoleIndex& ind) const {
return ind.spectatorData();
}
double IIgx2qqxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const {
double ret = alphaPDF(split);
double z = split.lastZ();
double ratio = sqr(split.lastPt()/split.scale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z );
- ret *= (!strictLargeN() ? 4./3. : 3./2.) * ( 1./x +sqr(1.-x)/x );
+ ret *= 0.5 * (!strictLargeN() ? 4./3. : 3./2.) * ( 1./x +sqr(1.-x)/x );
return ret;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IIgx2qqxDipoleKernel::persistentOutput(PersistentOStream & ) const {
}
void IIgx2qqxDipoleKernel::persistentInput(PersistentIStream & , int) {
}
ClassDescription<IIgx2qqxDipoleKernel> IIgx2qqxDipoleKernel::initIIgx2qqxDipoleKernel;
// Definition of the static class description member.
void IIgx2qqxDipoleKernel::Init() {
static ClassDocumentation<IIgx2qqxDipoleKernel> documentation
("IIgx2qqxDipoleKernel");
}
diff --git a/DipoleShower/Kinematics/DipoleSplittingKinematics.cc b/DipoleShower/Kinematics/DipoleSplittingKinematics.cc
--- a/DipoleShower/Kinematics/DipoleSplittingKinematics.cc
+++ b/DipoleShower/Kinematics/DipoleSplittingKinematics.cc
@@ -1,141 +1,167 @@
// -*- C++ -*-
//
// DipoleSplittingKinematics.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleSplittingKinematics class.
//
#include "DipoleSplittingKinematics.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <limits>
using namespace Herwig;
DipoleSplittingKinematics::DipoleSplittingKinematics()
: HandlerBase(), theIRCutoff(1.0*GeV),
theXMin(1.e-5), theJacobian(0.0),
theLastPt(0.0*GeV), theLastZ(0.0), theLastPhi(0.0),
theLastEmitterZ(1.0), theLastSpectatorZ(1.0),
theLastSplittingParameters() {}
DipoleSplittingKinematics::~DipoleSplittingKinematics() {}
void DipoleSplittingKinematics::persistentOutput(PersistentOStream & os) const {
os << ounit(theIRCutoff,GeV) << theXMin << theMCCheck;
}
void DipoleSplittingKinematics::persistentInput(PersistentIStream & is, int) {
is >> iunit(theIRCutoff,GeV) >> theXMin >> theMCCheck;
}
void DipoleSplittingKinematics::prepareSplitting(DipoleSplittingInfo& dInfo) {
dInfo.splittingKinematics(this);
if ( lastPt() > IRCutoff() )
dInfo.lastPt(lastPt());
else {
dInfo.lastPt(0.0*GeV);
dInfo.didStopEvolving();
}
dInfo.lastZ(lastZ());
dInfo.lastPhi(lastPhi());
dInfo.lastEmitterZ(lastEmitterZ());
dInfo.lastSpectatorZ(lastSpectatorZ());
dInfo.splittingParameters().resize(lastSplittingParameters().size());
copy(lastSplittingParameters().begin(),lastSplittingParameters().end(),
dInfo.splittingParameters().begin());
}
-Lorentz5Momentum DipoleSplittingKinematics::getKt (const Lorentz5Momentum& p1,
- const Lorentz5Momentum& p2,
- Energy pt,
- double phi) const {
+Lorentz5Momentum DipoleSplittingKinematics::getKt(const Lorentz5Momentum& p1,
+ const Lorentz5Momentum& p2,
+ Energy pt,
+ double phi,
+ bool spacelike) const {
- Boost beta = (p1+p2).findBoostToCM();
-
- Lorentz5Momentum p1c = p1;
+ Lorentz5Momentum P;
+ if ( !spacelike )
+ P = p1 + p2;
+ else
+ P = p1 - p2;
- if (beta.mag2() > Constants::epsilon) {
- p1c.boost(beta);
+ Energy2 Q2 = abs(P.m2());
+
+ Lorentz5Momentum Q =
+ !spacelike ?
+ Lorentz5Momentum(ZERO,ZERO,ZERO,sqrt(Q2),sqrt(Q2)) :
+ Lorentz5Momentum(ZERO,ZERO,sqrt(Q2),ZERO,-sqrt(Q2));
+
+ if ( spacelike && Q.z() < P.z() )
+ Q.setZ(-Q.z());
+
+ bool boost =
+ abs((P-Q).vect().mag2()/GeV2) > 1e-10 ||
+ abs((P-Q).t()/GeV) > 1e-5;
+
+ Lorentz5Momentum inFrame1;
+ if ( boost )
+ inFrame1 = p1 + ((P*p1-Q*p1)/(P*Q-Q.mass2()))*(P-Q);
+ else
+ inFrame1 = p1;
+
+ Energy ptx = inFrame1.x();
+ Energy pty = inFrame1.y();
+ Energy q = 2.*inFrame1.z();
+
+ Energy Qp = sqrt(4.*(sqr(ptx)+sqr(pty))+sqr(q));
+ Energy Qy = sqrt(4.*sqr(pty)+sqr(q));
+
+ double cPhi = cos(phi);
+ double sPhi = sqrt(1.-sqr(cPhi));
+ if ( phi > Constants::pi )
+ sPhi = -sPhi;
+
+ Lorentz5Momentum kt;
+
+ if ( !spacelike ) {
+ kt.setT(ZERO);
+ kt.setX(pt*Qy*cPhi/Qp);
+ kt.setY(-pt*(4*ptx*pty*cPhi/Qp+q*sPhi)/Qy);
+ kt.setZ(2.*pt*(-ptx*q*cPhi/Qp + pty*sPhi)/Qy);
+ } else {
+ kt.setT(2.*pt*(ptx*q*cPhi+pty*Qp*sPhi)/(q*Qy));
+ kt.setX(pt*(Qp*q*cPhi+4.*ptx*pty*sPhi)/(q*Qy));
+ kt.setY(pt*Qy*sPhi/q);
+ kt.setZ(ZERO);
}
-
- Lorentz5Momentum k (0.*GeV,0.*GeV,0.*GeV,0.*GeV);
-
- double ct = p1c.vect().unit().z();
- double st = sqrt(1.-ct*ct);
-
- double cphi = cos(phi);
- double sphi = sqrt(1.-cphi*cphi);
- if (phi > Constants::pi) sphi = -sphi;
-
- if (st > Constants::epsilon) {
- double cchi = p1c.vect().unit().x()/st;
- double schi = p1c.vect().unit().y()/st;
- k.setX((cphi*cchi*ct-sphi*schi)*pt);
- k.setY((cphi*schi*ct+sphi*cchi)*pt);
- k.setZ(-cphi*st*pt);
- } else {
- k.setX(pt*cphi);
- k.setY(pt*sphi);
- k.setZ(0.*GeV);
- }
-
- if (beta.mag2() > Constants::epsilon)
- k.boost(-beta);
-
- return k;
+
+ if ( boost )
+ kt = kt + ((P*kt-Q*kt)/(P*Q-Q.mass2()))*(P-Q);
+ kt.setMass(-pt);
+ kt.rescaleRho();
+
+ return kt;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
AbstractClassDescription<DipoleSplittingKinematics> DipoleSplittingKinematics::initDipoleSplittingKinematics;
// Definition of the static class description member.
void DipoleSplittingKinematics::Init() {
static ClassDocumentation<DipoleSplittingKinematics> documentation
("DipoleSplittingKinematics is the base class for dipole splittings "
"as performed in the dipole shower.");
static Parameter<DipoleSplittingKinematics,Energy> interfaceIRCutoff
("IRCutoff",
"The IR cutoff to be used by this splitting kinematics.",
&DipoleSplittingKinematics::theIRCutoff, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKinematics,double> interfaceXMin
("XMin",
"The minimum momentum fraction for incoming partons",
&DipoleSplittingKinematics::theXMin, 1.0e-5, 0.0, 1.0,
false, false, Interface::limited);
static Reference<DipoleSplittingKinematics,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingKinematics::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
}
diff --git a/DipoleShower/Kinematics/DipoleSplittingKinematics.h b/DipoleShower/Kinematics/DipoleSplittingKinematics.h
--- a/DipoleShower/Kinematics/DipoleSplittingKinematics.h
+++ b/DipoleShower/Kinematics/DipoleSplittingKinematics.h
@@ -1,475 +1,476 @@
// -*- C++ -*-
//
// DipoleSplittingKinematics.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleSplittingKinematics_H
#define HERWIG_DipoleSplittingKinematics_H
//
// This is the declaration of the DipoleSplittingKinematics class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Vectors/Lorentz5Vector.h"
#include "Herwig++/DipoleShower/Utility/DipoleMCCheck.h"
namespace Herwig {
using namespace ThePEG;
class DipoleIndex;
class DipoleSplittingInfo;
class DipoleSplittingKernel;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief DipoleSplittingKinematics is the base class for dipole splittings
* as performed in the dipole shower.
*
* @see \ref DipoleSplittingKinematicsInterfaces "The interfaces"
* defined for DipoleSplittingKinematics.
*/
class DipoleSplittingKinematics: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleSplittingKinematics();
/**
* The destructor.
*/
virtual ~DipoleSplittingKinematics();
//@}
public:
/**
* Return the boundaries in between the evolution
* variable random number is to be sampled; the lower
* cuoff is assumed to correspond to the infrared cutoff.
*/
virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const = 0;
/**
* Return the boundaries in between the momentum
* fraction random number is to be sampled.
*/
virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const = 0;
/**
* Return the dipole scale associated to the
* given pair of emitter and spectator. This
* should be the invariant mass or absolute value
* final/final or initial/initial and the absolute
* value of the momentum transfer for intial/final or
* final/initial dipoles.
*/
virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const = 0;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel& split) const = 0;
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex) const = 0;
/**
* Return the pt given a virtuality.
*/
virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const = 0;
/**
* Return the virtuality given a pt.
*/
virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const = 0;
/**
* Return the infrared cutoff.
*/
virtual Energy IRCutoff() const { return theIRCutoff; }
/**
* Return the minimum momentum fraction for
* incoming partons
*/
double xMin() const { return theXMin; }
/**
* Return the random number associated to
* the given pt.
*/
virtual double ptToRandom(Energy pt, Energy dScale,
const DipoleIndex& dIndex) const = 0;
/**
* Generate splitting variables given three random numbers
* and the momentum fractions of the emitter and spectator.
* Return true on success.
*/
virtual bool generateSplitting(double kappa, double xi, double phi,
DipoleSplittingInfo& dIndex) = 0;
/**
* For the splitting products present in the given dipole splitting
* info object calculate the kinematics parameters and return the
* propagator factor.
*/
virtual InvEnergy2 setKinematics(DipoleSplittingInfo&) const = 0;
/**
* For the splitting parameters given in the dipole splitting info
* object, calculate the phasespace Jacobian times the propagator
* factor.
*/
virtual double jacobianTimesPropagator(const DipoleSplittingInfo&,
Energy) const = 0;
/**
* Get the splitting phasespace weight associated to
* the last call to generateSplitting. This is taken to
* be the single particle phasespace times 16 \pi^2 divided
* by the relevant propagator invariant.
*/
double jacobian() const { return theJacobian; }
/**
* Return true, if this splitting kinematics
* class is capable of delivering an overestimate
* to the jacobian.
*/
virtual bool haveOverestimate() const { return false; }
/**
* Return the overestimated jacobian for the
* last generated parameters.
*/
virtual double jacobianOverestimate() const { return -1.; }
/**
* Return the last generated pt
*/
Energy lastPt() const { return theLastPt; }
/**
* Return the last generated momentum fraction.
*/
double lastZ() const { return theLastZ; }
/**
* Return the last generated azimuthal angle.
*/
double lastPhi() const { return theLastPhi; }
/**
* Return the momentum fraction, by which the emitter's
* momentum fraction should be divided after the splitting.
*/
double lastEmitterZ() const { return theLastEmitterZ; }
/**
* Return the momentum fraction, by which the spectator's
* momentum fraction should be divided after the splitting.
*/
double lastSpectatorZ() const { return theLastSpectatorZ; }
/**
* Return any additional parameters needed to
* evaluate the splitting kernel or to generate the
* full splitting.
*/
const vector<double>& lastSplittingParameters() const { return theLastSplittingParameters; }
/**
* Complete a DipoleSplittingInfo object with
* the parameters generated by the last call to
* generateSplitting()
*/
void prepareSplitting(DipoleSplittingInfo& dInfo);
public:
/**
* Generate the full kinematics given emitter and
* spectator momentum and a previously completeted
* DipoleSplittingInfo object.
*/
virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) = 0;
/**
* Return the emitter's momentum after the splitting.
*/
const Lorentz5Momentum& lastEmitterMomentum() const { return theEmitterMomentum; }
/**
* Return the spectator's momentum after the splitting.
*/
const Lorentz5Momentum& lastSpectatorMomentum() const { return theSpectatorMomentum; }
/**
* Return the emission's momentum.
*/
const Lorentz5Momentum& lastEmissionMomentum() const { return theEmissionMomentum; }
/*
* Return true, if there is a transformation which should
* be applied to all other final state particles except the ones
* involved in the splitting after having performed the splitting
*/
virtual bool doesTransform () const { return false; }
/*
* perform the transformation, if existing
*/
virtual Lorentz5Momentum transform (const Lorentz5Momentum& p) const { return p; }
protected:
/**
* Calculate a transverse momentum for the given momenta,
* invariant pt and azimuth.
*/
- Lorentz5Momentum getKt (const Lorentz5Momentum& p1,
- const Lorentz5Momentum& p2,
- Energy pt,
- double phi) const;
+ Lorentz5Momentum getKt(const Lorentz5Momentum& p1,
+ const Lorentz5Momentum& p2,
+ Energy pt,
+ double phi,
+ bool spacelike = false) const;
/**
* Set the splitting phasespace weight associated to
* the last call to generateSplitting. This is taken to
* be the single particle phasespace times 16 \pi^2 divided
* by the relevant propagator invariant.
*/
void jacobian(double w) { theJacobian = w; }
/**
* Set the last generated pt
*/
void lastPt(Energy p) { theLastPt = p; }
/**
* Set the last generated momentum fraction.
*/
void lastZ(double z) { theLastZ = z; }
/**
* Set the last generated azimuthal angle.
*/
void lastPhi(double p) { theLastPhi = p; }
/**
* Set the momentum fraction, by which the emitter's
* momentum fraction should be divided after the splitting.
*/
void lastEmitterZ(double z) { theLastEmitterZ = z; }
/**
* Set the momentum fraction, by which the spectator's
* momentum fraction should be divided after the splitting.
*/
void lastSpectatorZ(double z) { theLastSpectatorZ = z; }
/**
* Access any additional parameters needed to
* evaluate the splitting kernel or to generate the
* full splitting.
*/
vector<double>& splittingParameters() { return theLastSplittingParameters; }
/**
* Set the emitter's momentum after the splitting.
*/
void emitterMomentum(const Lorentz5Momentum& p) { theEmitterMomentum = p; }
/**
* Set the spectator's momentum after the splitting.
*/
void spectatorMomentum(const Lorentz5Momentum& p) { theSpectatorMomentum = p; }
/**
* Set the emission's momentum.
*/
void emissionMomentum(const Lorentz5Momentum& p) { theEmissionMomentum = p; }
public:
/**
* 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();
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);
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The infrared cutoff associated to this
* splitting kinematics.
*/
Energy theIRCutoff;
/**
* The minimum momentum fraction for
* incoming partons
*/
double theXMin;
/**
* The last calculated splitting phase space weight.
*/
double theJacobian;
/**
* The last generated pt
*/
Energy theLastPt;
/**
* The last generated momentum fraction.
*/
double theLastZ;
/**
* The last generated azimuthal angle.
*/
double theLastPhi;
/**
* The momentum fraction, by which the emitter's
* momentum fraction should be divided after the splitting.
*/
double theLastEmitterZ;
/**
* The momentum fraction, by which the spectator's
* momentum fraction should be divided after the splitting.
*/
double theLastSpectatorZ;
/**
* Any additional parameters needed to
* evaluate the splitting kernel or to generate the
* full splitting.
*/
vector<double> theLastSplittingParameters;
/**
* The emitter's momentum after the splitting.
*/
Lorentz5Momentum theEmitterMomentum;
/**
* The emission's momentum after the splitting.
*/
Lorentz5Momentum theEmissionMomentum;
/**
* The spectator's momentum after the splitting.
*/
Lorentz5Momentum theSpectatorMomentum;
protected:
/**
* Pointer to a check histogram object
*/
Ptr<DipoleMCCheck>::ptr theMCCheck;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class.
*/
static AbstractClassDescription<DipoleSplittingKinematics> initDipoleSplittingKinematics;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleSplittingKinematics & operator=(const DipoleSplittingKinematics &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleSplittingKinematics. */
template <>
struct BaseClassTrait<Herwig::DipoleSplittingKinematics,1> {
/** Typedef of the first base class of DipoleSplittingKinematics. */
typedef HandlerBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleSplittingKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleSplittingKinematics>
: public ClassTraitsBase<Herwig::DipoleSplittingKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleSplittingKinematics"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleSplittingKinematics is implemented. It may also include several, space-separated,
* libraries if the class DipoleSplittingKinematics depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleSplittingKinematics_H */
diff --git a/DipoleShower/Kinematics/FILightKinematics.cc b/DipoleShower/Kinematics/FILightKinematics.cc
--- a/DipoleShower/Kinematics/FILightKinematics.cc
+++ b/DipoleShower/Kinematics/FILightKinematics.cc
@@ -1,242 +1,242 @@
// -*- C++ -*-
//
// FILightKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FILightKinematics class.
//
#include "FILightKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
using namespace Herwig;
FILightKinematics::FILightKinematics()
: DipoleSplittingKinematics() {}
FILightKinematics::~FILightKinematics() {}
IBPtr FILightKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FILightKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> FILightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> FILightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emissionData()->id() != ParticleID::g )
return make_pair(0.5*(1.-c),0.5*(1.+c));
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
}
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
Energy FILightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return sqrt(2.*(pEmitter*pSpectator));
}
Energy FILightKinematics::ptMax(Energy dScale,
double, double specX,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
return dScale * sqrt((1.-specX)/specX) /2.;
}
Energy FILightKinematics::QMax(Energy dScale,
double, double specX,
const DipoleIndex&) const {
return dScale * sqrt((1.-specX)/specX);
}
Energy FILightKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
double z = split.lastZ();
return scale*sqrt(z*(1.-z));
}
Energy FILightKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
double z = split.lastZ();
return scale/sqrt(z*(1.-z));
}
double FILightKinematics::ptToRandom(Energy pt, Energy,
const DipoleIndex&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
bool FILightKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info) {
if ( info.spectatorX() < xMin() ) {
jacobian(0.0);
return false;
}
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
if ( pt > info.hardPt() ) {
jacobian(0.0);
return false;
}
double z;
double mapZJacobian;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
z = xi;
mapZJacobian = 1.;
} else {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
}
} else {
z = 1.-exp(-xi);
mapZJacobian = 1.-z;
}
double s = z*(1.-z);
double xs = info.spectatorX();
double x = 1. / ( 1. + sqr(pt/info.scale()) / s );
double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));
if ( pt < IRCutoff() ||
pt > info.hardPt() ||
z > zp || z < zm ||
x < xs ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
jacobian( 2. * mapZJacobian * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
lastPt(pt);
lastZ(z);
lastPhi(phi);
lastSpectatorZ(x);
if ( theMCCheck )
theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
InvEnergy2 FILightKinematics::setKinematics(DipoleSplittingInfo& split) const {
Lorentz5Momentum emitter = split.splitEmitter()->momentum();
Lorentz5Momentum emission = split.emission()->momentum();
Lorentz5Momentum spectator = split.splitSpectator()->momentum();
split.splittingKinematics(const_cast<FILightKinematics*>(this));
Energy2 scale = 2.*(- emission*emitter + emission*spectator + emitter*spectator);
split.scale(sqrt(scale));
double x =
scale / (2.*(emitter*spectator + emission*spectator));
double z = emitter*spectator / (emitter*spectator + emission*spectator);
split.lastPt(split.scale() * sqrt(z*(1.-z)*(1.-x)/x));
split.lastZ(z);
split.hardPt(split.lastPt());
if ( split.hardPt() > IRCutoff() ) {
split.continuesEvolving();
} else {
split.didStopEvolving();
}
return 1./(2.*x*(emitter*emission));
}
double FILightKinematics::
jacobianTimesPropagator(const DipoleSplittingInfo&,
Energy) const {
assert(false && "implementation missing");
return 0.;
}
void FILightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
Energy pt = dInfo.lastPt();
double z = dInfo.lastZ();
Lorentz5Momentum kt =
- getKt (pEmitter, pSpectator, pt, dInfo.lastPhi());
+ getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);
double ratio = sqr(pt/(-(pEmitter-pSpectator).m()));
double xInv = (1.+ratio/(z*(1.-z)));
Lorentz5Momentum em = z*pEmitter + (ratio/z)*pSpectator + kt;
em.setMass(0.*GeV);
em.rescaleEnergy();
Lorentz5Momentum emm = (1.-z)*pEmitter + (ratio/(1.-z))*pSpectator - kt;
emm.setMass(0.*GeV);
emm.rescaleEnergy();
Lorentz5Momentum spe = xInv*pSpectator;
spe.setMass(0.*GeV);
spe.rescaleEnergy();
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FILightKinematics::persistentOutput(PersistentOStream & ) const {
}
void FILightKinematics::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FILightKinematics> FILightKinematics::initFILightKinematics;
// Definition of the static class description member.
void FILightKinematics::Init() {
static ClassDocumentation<FILightKinematics> documentation
("FILightKinematics implements massless splittings "
"off a final-initial dipole.");
}
diff --git a/DipoleShower/Kinematics/IFLightKinematics.cc b/DipoleShower/Kinematics/IFLightKinematics.cc
--- a/DipoleShower/Kinematics/IFLightKinematics.cc
+++ b/DipoleShower/Kinematics/IFLightKinematics.cc
@@ -1,314 +1,324 @@
// -*- C++ -*-
//
// IFLightKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IFLightKinematics class.
//
#include "IFLightKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
using namespace Herwig;
IFLightKinematics::IFLightKinematics()
: DipoleSplittingKinematics(), theCollinearScheme(false) {}
IFLightKinematics::~IFLightKinematics() {}
IBPtr IFLightKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IFLightKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> IFLightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> IFLightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emitterData()->id() == ParticleID::g ) {
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
} else {
return make_pair(log(0.5*(1.-c)),log(0.5*(1.+c)));
}
}
if ( split.index().emitterData()->id() != ParticleID::g &&
split.emitterData()->id() != ParticleID::g ) {
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
return make_pair(0.5*(1.-c),0.5*(1.+c));
}
Energy IFLightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return sqrt(2.*(pEmitter*pSpectator));
}
Energy IFLightKinematics::ptMax(Energy dScale,
double emX, double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
- return dScale * sqrt(1.-emX) /2.;
+ return dScale * sqrt((1.-emX)/emX) /2.;
}
Energy IFLightKinematics::QMax(Energy,
double, double,
const DipoleIndex&) const {
assert(false && "add this");
return 0.0*GeV;
}
Energy IFLightKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
assert(false && "add this");
return 0.0*GeV;
}
Energy IFLightKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
assert(false && "add this");
return 0.0*GeV;
}
double IFLightKinematics::ptToRandom(Energy pt, Energy,
const DipoleIndex&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
bool IFLightKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info) {
if ( info.emitterX() < xMin() ) {
jacobian(0.0);
return false;
}
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
- if ( pt > info.hardPt() ) {
+ if ( sqr(pt) > sqr(info.hardPt())/(1.+4.*sqr(info.hardPt()/info.scale())) ) {
jacobian(0.0);
return false;
}
double z = 0.;
double mapZJacobian = 0.;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emitterData()->id() == ParticleID::g ) {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
} else {
z = exp(xi);
mapZJacobian = z;
}
}
if ( info.index().emitterData()->id() != ParticleID::g ) {
if ( info.emitterData()->id() != ParticleID::g ) {
z = 1.-exp(-xi);
mapZJacobian = 1.-z;
} else {
z = xi;
mapZJacobian = 1.;
}
}
double ratio = sqr(pt/info.scale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z - ratio );
double u = ratio/(1.-z);
if ( x < 0. || x > 1. || u > 1. ) {
jacobian(0.0);
return false;
}
double xe = info.emitterX();
- double zp = 0.5*( 1.+ xe +
- (1.-xe)*sqrt(1.-sqr(pt/info.hardPt()) ) );
- double zm = 0.5*( 1.+ xe -
- (1.-xe)*sqrt(1.-sqr(pt/info.hardPt()) ) );
+ double zpx = 0.5*( 1.+ xe +
+ (1.-xe)*sqrt(1.-sqr(2.*pt/info.scale())/(1.-xe) ) );
+ double zmx = 0.5*( 1.+ xe -
+ (1.-xe)*sqrt(1.-sqr(2.*pt/info.scale())/(1.-xe) ) );
+
+ double xq = sqr(pt/info.hardPt());
+
+ double zpq = 0.5*( 1.+ xq +
+ (1.-xq)*sqrt(1.-sqr(2.*pt/info.scale())/(1.-xq) ) );
+ double zmq = 0.5*( 1.+ xq -
+ (1.-xq)*sqrt(1.-sqr(2.*pt/info.scale())/(1.-xq) ) );
+
+ double zp = min(zpx,zpq);
+ double zm = max(zmx,zmq);
if ( pt < IRCutoff() ||
pt > info.hardPt() ||
z > zp || z < zm ||
x < xe ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
jacobian(2. * mapZJacobian * (1.-z)/(z*(1.-z)-ratio) * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
lastPt(pt);
lastZ(z);
lastPhi(phi);
lastEmitterZ(x);
if ( theMCCheck )
theMCCheck->book(info.emitterX(),1.,info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
InvEnergy2 IFLightKinematics::setKinematics(DipoleSplittingInfo& split) const {
Lorentz5Momentum emitter = split.splitEmitter()->momentum();
Lorentz5Momentum emission = split.emission()->momentum();
Lorentz5Momentum spectator = split.splitSpectator()->momentum();
split.splittingKinematics(const_cast<IFLightKinematics*>(this));
Energy2 scale = 2.*(emission*emitter - emission*spectator + emitter*spectator);
split.scale(sqrt(scale));
double x =
scale / (2.*(emitter*emission + emitter*spectator));
double u = emitter*emission / (emitter*emission + emitter*spectator);
split.lastPt(split.scale() * sqrt(u*(1.-u)*(1.-x)));
split.lastZ(u+x*(1.-u));
split.hardPt(split.lastPt());
if ( split.hardPt() > IRCutoff() ) {
split.continuesEvolving();
} else {
split.didStopEvolving();
}
return 1./(2.*x*(emitter*emission));
}
double IFLightKinematics::
jacobianTimesPropagator(const DipoleSplittingInfo&,
Energy) const {
assert(false && "implementation missing");
return 0.;
}
void IFLightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
Energy pt = dInfo.lastPt();
double ratio = sqr(pt/-(pEmitter-pSpectator).m());
double z = dInfo.lastZ();
double x = (z*(1.-z)-ratio)/(1.-z-ratio);
double u = ratio / (1.-z);
pt = -(pEmitter-pSpectator).m()*sqrt(u*(1.-u)*(1.-x)/x);
Lorentz5Momentum kt =
- getKt (pEmitter, pSpectator, pt, dInfo.lastPhi());
+ getKt (pEmitter, pSpectator, pt, dInfo.lastPhi(),true);
Lorentz5Momentum em;
Lorentz5Momentum emm;
Lorentz5Momentum spe;
if ( !theCollinearScheme &&
x > u && (1.-x)/(x-u) < 1. ) {
em =
((1.-u)/(x-u))*pEmitter + ((u/x)*(1.-x)/(x-u))*pSpectator - kt/(x-u);
em.setMass(0.*GeV);
em.rescaleEnergy();
emm =
((1.-x)/(x-u))*pEmitter + ((u/x)*(1.-u)/(x-u))*pSpectator - kt/(x-u);
emm.setMass(0.*GeV);
emm.rescaleEnergy();
spe =
(1.-u/x)*pSpectator;
spe.setMass(0.*GeV);
spe.rescaleEnergy();
} else {
em = (1./x)*pEmitter;
emm = ((1.-x)*(1.-u)/x)*pEmitter + u*pSpectator + kt;
emm.setMass(0.*GeV);
emm.rescaleEnergy();
spe = ((1.-x)*u/x)*pEmitter + (1.-u)*pSpectator - kt;
spe.setMass(0.*GeV);
spe.rescaleEnergy();
}
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFLightKinematics::persistentOutput(PersistentOStream & os) const {
os << theCollinearScheme;
}
void IFLightKinematics::persistentInput(PersistentIStream & is, int) {
is >> theCollinearScheme;
}
ClassDescription<IFLightKinematics> IFLightKinematics::initIFLightKinematics;
// Definition of the static class description member.
void IFLightKinematics::Init() {
static ClassDocumentation<IFLightKinematics> documentation
("IFLightKinematics implements massless splittings "
"off a initial-final dipole.");
static Switch<IFLightKinematics,bool> interfaceCollinearScheme
("CollinearScheme",
"[experimental] Switch on or off the collinear scheme",
&IFLightKinematics::theCollinearScheme, false, false, false);
static SwitchOption interfaceCollinearSchemeOn
(interfaceCollinearScheme,
"On",
"Switch on the collinear scheme.",
true);
static SwitchOption interfaceCollinearSchemeOff
(interfaceCollinearScheme,
"Off",
"Switch off the collinear scheme",
false);
interfaceCollinearScheme.rank(-1);
}
diff --git a/DipoleShower/Kinematics/IILightKinematics.cc b/DipoleShower/Kinematics/IILightKinematics.cc
--- a/DipoleShower/Kinematics/IILightKinematics.cc
+++ b/DipoleShower/Kinematics/IILightKinematics.cc
@@ -1,353 +1,364 @@
// -*- C++ -*-
//
// IILightKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IILightKinematics class.
//
#include "IILightKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
using namespace Herwig;
IILightKinematics::IILightKinematics()
: DipoleSplittingKinematics(), theCollinearScheme(false), didCollinear(false) {}
IILightKinematics::~IILightKinematics() {}
IBPtr IILightKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IILightKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> IILightKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> IILightKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emitterData()->id() == ParticleID::g ) {
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
} else {
return make_pair(log(0.5*(1.-c)),log(0.5*(1.+c)));
}
}
if ( split.index().emitterData()->id() != ParticleID::g &&
split.emitterData()->id() != ParticleID::g ) {
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
return make_pair(0.5*(1.-c),0.5*(1.+c));
}
Energy IILightKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return sqrt(2.*(pEmitter*pSpectator));
}
Energy IILightKinematics::ptMax(Energy dScale,
double emX, double specX,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
double tau = emX*specX;
- return (1.-tau) * dScale / 2.;
+ return (1.-tau) * dScale / (2.*sqrt(tau));
}
Energy IILightKinematics::QMax(Energy,
double, double,
const DipoleIndex&) const {
assert(false && "add this");
return 0.0*GeV;
}
Energy IILightKinematics::PtFromQ(Energy, const DipoleSplittingInfo&) const {
assert(false && "add this");
return 0.0*GeV;
}
Energy IILightKinematics::QFromPt(Energy, const DipoleSplittingInfo&) const {
assert(false && "add this");
return 0.0*GeV;
}
double IILightKinematics::ptToRandom(Energy pt, Energy,
const DipoleIndex&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
bool IILightKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info) {
if ( info.emitterX() < xMin() ||
info.spectatorX() < xMin() ) {
jacobian(0.0);
return false;
}
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
- if ( pt > info.hardPt() ) {
+ double r = sqr(info.hardPt()/info.scale());
+ if ( sqr(pt) > sqr(info.hardPt())*r*((2.+1./r)-2.*sqrt(1.+1./r)) ) {
jacobian(0.0);
return false;
}
double z = 0.;
double mapZJacobian = 0.;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emitterData()->id() == ParticleID::g ) {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
} else {
z = exp(xi);
mapZJacobian = z;
}
}
if ( info.index().emitterData()->id() != ParticleID::g ) {
if ( info.emitterData()->id() != ParticleID::g ) {
z = 1.-exp(-xi);
mapZJacobian = (1.-z);
} else {
z = xi;
mapZJacobian = 1.;
}
}
double ratio = sqr(pt/info.scale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z );
double v = ratio / (1.-z);
if ( x < 0. || x > 1. || v > 1. || v > 1.-x ) {
jacobian(0.0);
return false;
}
double tau = info.emitterX()*info.spectatorX();
- double zp = 0.5*( 1.+ tau +
- (1.-tau)*sqrt(1.-sqr(pt/info.hardPt()) ) );
- double zm = 0.5*( 1.+ tau -
- (1.-tau)*sqrt(1.-sqr(pt/info.hardPt()) ) );
+ double zpx = 0.5*( 1.+ tau +
+ (1.-tau)*sqrt(1.-sqr(2.*pt/((1.-tau)*info.scale()))) );
+ double zmx = 0.5*( 1.+ tau -
+ (1.-tau)*sqrt(1.-sqr(2.*pt/((1.-tau)*info.scale()))) );
+
+ double xq = sqr(pt/info.hardPt());
+
+ double zpq = 0.5*( 1.+ xq +
+ (1.-xq)*sqrt(1.-sqr(2.*pt/((1.-xq)*info.scale()))) );
+ double zmq = 0.5*( 1.+ xq -
+ (1.-xq)*sqrt(1.-sqr(2.*pt/((1.-xq)*info.scale()))) );
+
+ double zp = min(zpx,zpq);
+ double zm = max(zmx,zmq);
if ( pt < IRCutoff() ||
pt > info.hardPt() ||
z > zp || z < zm ) {
jacobian(0.0);
return false;
}
if ( !theCollinearScheme &&
(1.-v-x)/(v+x) < 1. ) {
if ( (x+v) < info.emitterX() ||
x/(x+v) < info.spectatorX() ) {
jacobian(0.0);
return false;
}
} else {
if ( x < info.emitterX() ) {
jacobian(0.0);
return false;
}
}
double phi = 2.*Constants::pi*rphi;
jacobian(2. * mapZJacobian * (1.-z)/(z*(1.-z)-ratio) * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
lastPt(pt);
lastZ(z);
lastPhi(phi);
if ( !theCollinearScheme &&
(1.-v-x)/(v+x) < 1. ) {
lastEmitterZ(x+v);
lastSpectatorZ(x/(x+v));
} else {
lastEmitterZ(x);
lastSpectatorZ(1.);
}
if ( theMCCheck )
theMCCheck->book(info.emitterX(),info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
InvEnergy2 IILightKinematics::setKinematics(DipoleSplittingInfo& split) const {
Lorentz5Momentum emitter = split.splitEmitter()->momentum();
Lorentz5Momentum emission = split.emission()->momentum();
Lorentz5Momentum spectator = split.splitSpectator()->momentum();
split.splittingKinematics(const_cast<IILightKinematics*>(this));
Energy2 scale = 2.*(-emission*emitter - emission*spectator + emitter*spectator);
split.scale(sqrt(scale));
double x = scale/(2.*(emitter*spectator));
double v = (emitter*emission)/(emitter*spectator);
split.lastZ(v+x);
split.lastPt(split.scale() * sqrt(v*(1.-x-v)));
split.hardPt(split.lastPt());
if ( split.hardPt() > IRCutoff() ) {
split.continuesEvolving();
} else {
split.didStopEvolving();
}
return 1./(2.*x*(emitter*emission));
}
double IILightKinematics::
jacobianTimesPropagator(const DipoleSplittingInfo&,
Energy) const {
assert(false && "implementation missing");
return 0.;
}
void IILightKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
Energy pt = dInfo.lastPt();
double z = dInfo.lastZ();
double ratio = sqr(pt/(pEmitter+pSpectator).m());
double x = ( z*(1.-z) - ratio ) / ( 1. - z );
double v = ratio / (1.-z);
pt = sqrt(v*(1.-x-v)/x) * (pEmitter+pSpectator).m();
Lorentz5Momentum kt =
getKt (pEmitter, pSpectator, pt, dInfo.lastPhi());
if ( !theCollinearScheme &&
(1.-v-x)/(v+x) < 1. ) {
Lorentz5Momentum em =
(1./(v+x))*pEmitter+(v*(1.-v-x)/(x*(x+v)))*pSpectator+kt/(x+v);
em.setMass(0.*GeV);
em.rescaleEnergy();
Lorentz5Momentum emm =
((1.-v-x)/(v+x))*pEmitter+(v/(x*(x+v)))*pSpectator+kt/(x+v);
emm.setMass(0.*GeV);
emm.rescaleEnergy();
Lorentz5Momentum spe =
(1.+v/x)*pSpectator;
spe.setMass(0.*GeV);
spe.rescaleEnergy();
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
didCollinear = false;
} else {
Lorentz5Momentum em =
(1./x)*pEmitter;
em.setMass(0.*GeV);
em.rescaleEnergy();
Lorentz5Momentum emm =
((1.-x-v)/x)*pEmitter+v*pSpectator+kt;
emm.setMass(0.*GeV);
emm.rescaleEnergy();
Lorentz5Momentum spe =
pSpectator;
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
K = em + spe - emm;
K2 = K.m2();
Ktilde = pEmitter + pSpectator;
KplusKtilde = K + Ktilde;
KplusKtilde2 = KplusKtilde.m2();
didCollinear = true;
}
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IILightKinematics::persistentOutput(PersistentOStream & os) const {
os << theCollinearScheme;
}
void IILightKinematics::persistentInput(PersistentIStream & is, int) {
is >> theCollinearScheme;
}
ClassDescription<IILightKinematics> IILightKinematics::initIILightKinematics;
// Definition of the static class description member.
void IILightKinematics::Init() {
static ClassDocumentation<IILightKinematics> documentation
("IILightKinematics implements massless splittings "
"off an initial-initial dipole.");
static Switch<IILightKinematics,bool> interfaceCollinearScheme
("CollinearScheme",
"[experimental] Switch on or off the collinear scheme",
&IILightKinematics::theCollinearScheme, false, false, false);
static SwitchOption interfaceCollinearSchemeOn
(interfaceCollinearScheme,
"On",
"Switch on the collinear scheme.",
true);
static SwitchOption interfaceCollinearSchemeOff
(interfaceCollinearScheme,
"Off",
"Switch off the collinear scheme",
false);
interfaceCollinearScheme.rank(-1);
}
diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
@@ -1,391 +1,393 @@
// -*- C++ -*-
//
// MatchboxAmplitude.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 MatchboxAmplitude class.
//
#include "MatchboxAmplitude.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinorHelicity.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SU2Helper.h"
#include "MatchboxMEBase.h"
using namespace Herwig;
MatchboxAmplitude::MatchboxAmplitude()
: Amplitude(), theColourBasisDim(0),
calculateTrees(true), calculateLoops(true) {}
MatchboxAmplitude::~MatchboxAmplitude() {}
void MatchboxAmplitude::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theColourBasis << theColourBasisDim
<< theCrossingMap << theColourMap << theCrossingSigns
<< theAmplitudePartonData << theNLight;
}
void MatchboxAmplitude::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theColourBasis >> theColourBasisDim
>> theCrossingMap >> theColourMap >> theCrossingSigns
>> theAmplitudePartonData >> theNLight;
}
void MatchboxAmplitude::cloneDependencies(const std::string&) {
}
Ptr<MatchboxMEBase>::ptr MatchboxAmplitude::makeME(const vector<PDVector>&) const {
return new_ptr(MatchboxMEBase());
}
void MatchboxAmplitude::dumpInfo(const string& prefix) const {
generator()->log() << prefix << fullName()
<< " [" << this << "]\n";
generator()->log() << prefix << " | XComb " << lastXCombPtr()
<< " for ";
if ( lastXCombPtr() ) {
for ( cPDVector::const_iterator p = lastXComb().mePartonData().begin();
p != lastXComb().mePartonData().end(); ++p ) {
generator()->log() << (**p).PDGName() << " ";
}
}
generator()->log() << "\n";
}
struct orderPartonData {
bool operator()(const pair<tcPDPtr,int>& a,
const pair<tcPDPtr,int>& b) const {
if ( a.first == b.first )
return a.second < b.second;
int acolour = a.first->iColour();
int bcolour = b.first->iColour();
if ( abs(acolour) != abs(bcolour) )
return abs(acolour) < abs(bcolour);
if ( a.first->iSpin() != b.first->iSpin() )
return a.first->iSpin() < b.first->iSpin();
int acharge = a.first->iCharge();
int bcharge = b.first->iCharge();
if ( abs(acharge) != abs(bcharge) )
return abs(acharge) < abs(bcharge);
if ( abs(a.first->id()) != abs(b.first->id()) )
return abs(a.first->id()) < abs(b.first->id());
return a.first->id() > b.first->id();
}
};
void MatchboxAmplitude::fillCrossingMap(size_t shift) {
theLastCrossingMap = crossingMap().find(lastXCombPtr());
if ( theLastCrossingMap != crossingMap().end() ) {
assert(amplitudePartonData().find(lastXCombPtr()) != amplitudePartonData().end());
assert(colourMap().find(lastXCombPtr()) != colourMap().end());
theLastAmplitudePartonData = amplitudePartonData().find(lastXCombPtr());
theLastColourMap = colourMap().find(lastXCombPtr());
assert(crossingSigns().find(lastXCombPtr()) != crossingSigns().end());
theLastCrossingSign = crossingSigns().find(lastXCombPtr())->second;
return;
} else {
crossingMap()[lastXCombPtr()] = vector<int>(mePartonData().size());
amplitudePartonData()[lastXCombPtr()] = cPDVector(mePartonData().size());
colourMap()[lastXCombPtr()] = map<size_t,size_t>();
theLastCrossingMap = crossingMap().find(lastXCombPtr());
theLastAmplitudePartonData = amplitudePartonData().find(lastXCombPtr());
theLastColourMap = colourMap().find(lastXCombPtr());
}
double csign = 1.;
set<pair<tcPDPtr,int>,orderPartonData > processLegs;
for ( unsigned int l = 0; l < mePartonData().size(); ++l ) {
if ( l > 1 )
processLegs.insert(make_pair(mePartonData()[l],l));
else {
if ( mePartonData()[l]->CC() ) {
processLegs.insert(make_pair(mePartonData()[l]->CC(),l));
if ( mePartonData()[l]->iSpin() == PDT::Spin1Half )
csign *= -1.;
} else {
processLegs.insert(make_pair(mePartonData()[l],l));
}
}
}
lastCrossingSign(csign);
crossingSigns()[lastXCombPtr()] = csign;
set<pair<tcPDPtr,int> > amplitudeLegs;
int ampCount = 0;
// process legs are already sorted, we only need to arrange for
// adjacent particles and anti-particles
while ( !processLegs.empty() ) {
set<pair<tcPDPtr,int>,orderPartonData >::iterator next
= processLegs.begin();
while ( next->first->id() < 0 ) {
if ( ++next == processLegs.end() )
break;
}
assert(next != processLegs.end());
lastCrossingMap()[ampCount] = next->second - shift;
amplitudeLegs.insert(make_pair(next->first,ampCount));
tcPDPtr check = next->first;
processLegs.erase(next);
++ampCount;
if ( check->CC() ) {
set<pair<tcPDPtr,int>,orderPartonData>::iterator checkcc
= processLegs.end();
for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
c != processLegs.end(); ++c ) {
if ( c->first == check->CC() ) {
checkcc = c; break;
}
}
if ( checkcc == processLegs.end() )
for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
c != processLegs.end(); ++c ) {
assert(SU2Helper::SU2CC(check)->CC());
if ( c->first == SU2Helper::SU2CC(check)->CC() ) {
checkcc = c; break;
}
}
assert(checkcc != processLegs.end());
lastCrossingMap()[ampCount] = checkcc->second - shift;
amplitudeLegs.insert(make_pair(checkcc->first,ampCount));
processLegs.erase(checkcc);
++ampCount;
}
}
for ( set<pair<tcPDPtr,int> >::const_iterator l = amplitudeLegs.begin();
l != amplitudeLegs.end(); ++l )
lastAmplitudePartonData()[l->second] = l->first;
if ( colourBasis() ) {
assert(colourBasis()->indexMap().find(mePartonData()) !=
colourBasis()->indexMap().end());
const map<size_t,size_t> colourCross =
colourBasis()->indexMap().find(mePartonData())->second;
for ( size_t k = 0; k < lastCrossingMap().size(); ++k ) {
if ( colourCross.find(lastCrossingMap()[k]) !=
colourCross.end() ) {
lastColourMap()[k] = colourCross.find(lastCrossingMap()[k])->second;
}
}
}
}
Lorentz5Momentum MatchboxAmplitude::amplitudeMomentum(int i) const {
int iCrossed = lastCrossingMap()[i];
Lorentz5Momentum res = meMomenta()[iCrossed];
return iCrossed > 1 ? res : -res;
}
set<vector<int> > MatchboxAmplitude::generateHelicities() const {
set<vector<int> > res;
vector<int> current(lastAmplitudePartonData().size());
doGenerateHelicities(res,current,0);
return res;
}
void MatchboxAmplitude::doGenerateHelicities(set<vector<int> >& res,
vector<int>& current,
size_t pos) const {
if ( pos == lastAmplitudePartonData().size() ) {
res.insert(current);
return;
}
if ( lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin0 ||
( lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin1 &&
lastAmplitudePartonData()[pos]->mass() != ZERO ) ) {
current[pos] = 0;
doGenerateHelicities(res,current,pos+1);
} else if ( lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin1Half ||
lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin1 ) {
current[pos] = 1;
doGenerateHelicities(res,current,pos+1);
current[pos] = -1;
doGenerateHelicities(res,current,pos+1);
}
}
void MatchboxAmplitude::prepareAmplitudes() {
if ( !calculateTrees )
return;
if ( lastAmplitudes().empty() ) {
set<vector<int> > helicities = generateHelicities();
for ( set<vector<int> >::const_iterator h = helicities.begin();
h != helicities.end(); ++h )
lastAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
lastLargeNAmplitudes() = lastAmplitudes();
}
AmplitudeIterator amp = lastAmplitudes().begin();
AmplitudeIterator lamp = lastLargeNAmplitudes().begin();
for ( ;amp != lastAmplitudes().end(); ++amp, ++lamp ) {
for ( size_t k = 0; k < colourBasisDim(); ++k )
amp->second(k) = evaluate(k,amp->first,lamp->second(k));
}
calculateTrees = false;
}
void MatchboxAmplitude::prepareOneLoopAmplitudes() {
if ( !calculateLoops )
return;
if ( lastOneLoopAmplitudes().empty() ) {
set<vector<int> > helicities = generateHelicities();
for ( set<vector<int> >::const_iterator h = helicities.begin();
h != helicities.end(); ++h )
lastOneLoopAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
}
for ( AmplitudeIterator amp = lastOneLoopAmplitudes().begin();
amp != lastOneLoopAmplitudes().end(); ++amp ) {
for ( size_t k = 0; k < colourBasisDim(); ++k )
amp->second(k) = evaluateOneLoop(k,amp->first);
}
calculateLoops = false;
}
Complex MatchboxAmplitude::value(const tcPDVector&,
const vector<Lorentz5Momentum>&,
const vector<int>&) {
assert(false && "ThePEG::Amplitude interface is not sufficient at the moment.");
throw Exception() << "ThePEG::Amplitude interface is not sufficient at the moment."
<< Exception::abortnow;
return 0.;
}
double MatchboxAmplitude::colourCorrelatedME2(pair<int,int> ij) const {
double Nc = generator()->standardModel()->Nc();
double cfac = mePartonData()[ij.first]->id() == ParticleID::g ? Nc : (sqr(Nc)-1.)/(2.*Nc);
- return colourBasis()->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes())/cfac;
+ return
+ lastCrossingSign()*colourBasis()->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes())/cfac;
}
// compare int vectors modulo certain element
// which needs to differe between the two
bool equalsModulo(unsigned int i, const vector<int>& a, const vector<int>& b) {
assert(a.size()==b.size());
if ( a[i] == b[i] )
return false;
for ( unsigned int k = 0; k < a.size(); ++k ) {
if ( k == i )
continue;
if ( a[k] != b[k] )
return false;
}
return true;
}
double MatchboxAmplitude::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
using namespace SpinorHelicity;
Lorentz5Momentum p = meMomenta()[ij.first];
Lorentz5Momentum n = meMomenta()[ij.second];
LorentzVector<complex<Energy> > num =
PlusSpinorCurrent(PlusConjugateSpinor(n),MinusSpinor(p)).eval();
complex<Energy> den =
sqrt(2.)*PlusSpinorProduct(PlusConjugateSpinor(n),PlusSpinor(p)).eval();
LorentzVector<Complex> polarization(num.x()/den,num.y()/den,num.z()/den,num.t()/den);
Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale()));
double avg =
colourCorrelatedME2(ij)*(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)*norm(pFactor));
Complex corr = 0.;
int iCrossed = -1;
for ( unsigned int k = 0; k < lastCrossingMap().size(); ++k )
if ( lastCrossingMap()[k] == ij.first ) {
iCrossed = k;
break;
}
assert(iCrossed >= 0);
set<const CVector*> done;
for ( AmplitudeConstIterator a = theLastAmplitudes.begin();
a != theLastAmplitudes.end(); ++a ) {
if ( done.find(&(a->second)) != done.end() )
continue;
AmplitudeConstIterator b = theLastAmplitudes.begin();
while ( !equalsModulo(iCrossed,a->first,b->first) )
if ( ++b == theLastAmplitudes.end() )
break;
if ( b == theLastAmplitudes.end() || done.find(&(b->second)) != done.end() )
continue;
done.insert(&(a->second)); done.insert(&(b->second));
if ( a->first[iCrossed] == 1 )
swap(a,b);
corr += colourBasis()->interference(mePartonData(),a->second,b->second);
}
double Nc = generator()->standardModel()->Nc();
double cfac = mePartonData()[ij.first]->id() == ParticleID::g ? Nc : (sqr(Nc)-1.)/(2.*Nc);
- return avg + 2.*(c.scale() > ZERO ? 1. : -1.)*real(corr*sqr(pFactor))/cfac;
+ return
+ avg + 2.*lastCrossingSign()*(c.scale() > ZERO ? 1. : -1.)*real(corr*sqr(pFactor))/cfac;
}
void MatchboxAmplitude::Init() {
static ClassDocumentation<MatchboxAmplitude> documentation
("MatchboxAmplitude is the base class for amplitude "
"implementations inside Matchbox.");
static Reference<MatchboxAmplitude,ColourBasis> interfaceColourBasis
("ColourBasis",
"Set the colour basis implementation.",
&MatchboxAmplitude::theColourBasis, false, false, true, true, 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<MatchboxAmplitude,Amplitude>
describeMatchboxAmplitude("Herwig::MatchboxAmplitude", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
@@ -1,569 +1,571 @@
// -*- C++ -*-
//
// MatchboxAmplitude.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_MatchboxAmplitude_H
#define HERWIG_MatchboxAmplitude_H
//
// This is the declaration of the MatchboxAmplitude class.
//
#include "ThePEG/MatrixElement/Amplitude.h"
#include "ThePEG/Handlers/LastXCombInfo.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/ColourBasis.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
namespace Herwig {
using namespace ThePEG;
class MatchboxMEBase;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxAmplitude is the base class for amplitude
* implementations inside Matchbox.
*
* @see \ref MatchboxAmplitudeInterfaces "The interfaces"
* defined for MatchboxAmplitude.
*/
class MatchboxAmplitude: public Amplitude, public LastXCombInfo<StandardXComb> {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxAmplitude();
/**
* The destructor.
*/
virtual ~MatchboxAmplitude();
//@}
public:
typedef map<vector<int>,CVector> AmplitudeMap;
typedef map<vector<int>,CVector>::iterator AmplitudeIterator;
typedef map<vector<int>,CVector>::const_iterator AmplitudeConstIterator;
/**
* Return the amplitude. Needs to be implemented from
* ThePEG::Amplitude but is actually ill-defined, as colours of the
* external particles are not specified. To this extent, this
* implementation just asserts.
*/
virtual Complex value(const tcPDVector & particles,
const vector<Lorentz5Momentum> & momenta,
const vector<int> & helicities);
/** @name Subprocess information */
//@{
/**
* Return true, if this amplitude can handle the given process.
*/
virtual bool canHandle(const PDVector&) const { return false; }
/**
* Return a ME instance appropriate for this amplitude and the given
* subprocesses
*/
Ptr<MatchboxMEBase>::ptr makeME(const vector<PDVector>&) const;
/**
* Return the amplitude parton data.
*/
const cPDVector& lastAmplitudePartonData() const { return theLastAmplitudePartonData->second; }
/**
* Access the amplitude parton data.
*/
cPDVector& lastAmplitudePartonData() { return theLastAmplitudePartonData->second; }
/**
* Access the amplitude parton data.
*/
map<tStdXCombPtr,cPDVector>& amplitudePartonData() { return theAmplitudePartonData; }
/**
* Return the number of light flavours
*/
unsigned int nLight() const { return theNLight; }
/**
* Set the number of light flavours
*/
void nLight(unsigned int n) { theNLight = n; }
/**
* Set the (tree-level) order in \f$g_S\f$ in which this matrix
* element should be evaluated.
*/
virtual void orderInGs(unsigned int) {}
/**
* Return the (tree-level) order in \f$g_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGs() const = 0;
/**
* Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element should be evaluated.
*/
virtual void orderInGem(unsigned int) {}
/**
* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInGem() const = 0;
/**
* Return the Herwig++ StandardModel object
*/
Ptr<StandardModel>::tcptr standardModel() {
if ( !theStandardModel )
theStandardModel =
dynamic_ptr_cast<Ptr<StandardModel>::tcptr>(HandlerBase::standardModel());
return theStandardModel;
}
//@}
/** @name Colour basis. */
//@{
/**
* Return the colour basis.
*/
Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }
/**
* Set the colour basis dimensionality.
*/
void colourBasisDim(size_t dim) { theColourBasisDim = dim; }
/**
* Get the colour basis dimensionality.
*/
size_t colourBasisDim() const { return theColourBasisDim; }
/**
* Return true, if this amplitude will not require colour correlations.
*/
virtual bool noCorrelations() const { return !haveOneLoop(); }
/**
* Return true, if the colour basis is capable of assigning colour
* flows.
*/
virtual bool haveColourFlows() const {
return colourBasis() ? colourBasis()->haveColourFlows() : false;
}
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
*/
virtual Selector<const ColourLines *> colourGeometries(tcDiagPtr diag) const {
return
haveColourFlows() ?
theColourBasis->colourGeometries(diag,lastLargeNAmplitudes()) :
Selector<const ColourLines *>();
}
/**
* Return the colour crossing information as filled by the last call to
* fillCrossingMap(...), mapping amplitude ids to colour basis ids.
*/
const map<size_t,size_t>& lastColourMap() const { return theLastColourMap->second; }
/**
* Access the colour crossing information.
*/
map<size_t,size_t>& lastColourMap() { return theLastColourMap->second; }
/**
* Access the colour crossing information.
*/
map<tStdXCombPtr,map<size_t,size_t> >& colourMap() { return theColourMap; }
//@}
/** @name Phasespace point, crossing and helicities */
//@{
/**
* Set the xcomb object.
*/
virtual void setXComb(tStdXCombPtr xc) {
theLastXComb = xc;
fillCrossingMap();
}
/**
* Return the momentum as crossed appropriate for this amplitude.
*/
Lorentz5Momentum amplitudeMomentum(int) const;
/**
* Perform a normal ordering of external legs and fill the
* crossing information as. This default implementation sorts
* lexicographically in (abs(colour)/spin/abs(charge)), putting pairs
* of particles/anti-particles where possible.
*/
virtual void fillCrossingMap(size_t shift = 0);
/**
* Return the crossing sign.
*/
double lastCrossingSign() const { return theLastCrossingSign; }
/**
* Set the crossing sign.
*/
void lastCrossingSign(double s) { theLastCrossingSign = s; }
/**
* Return the crossing information as filled by the last call to
* fillCrossingMap(...), mapping amplitude ids to process ids.
*/
const vector<int>& lastCrossingMap() const { return theLastCrossingMap->second; }
/**
* Access the crossing information.
*/
vector<int>& lastCrossingMap() { return theLastCrossingMap->second; }
/**
* Access the crossing information.
*/
map<tStdXCombPtr,vector<int> >& crossingMap() { return theCrossingMap; }
/**
* Access the crossing signs.
*/
map<tStdXCombPtr,double>& crossingSigns() { return theCrossingSigns; }
/**
* Generate the helicity combinations.
*/
virtual set<vector<int> > generateHelicities() const;
//@}
/** @name Tree-level amplitudes */
//@{
/**
* Calculate the tree level amplitudes for the phasespace point
* stored in lastXComb.
*/
virtual void prepareAmplitudes();
/**
* Return last evaluated helicity amplitudes.
*/
const AmplitudeMap& lastAmplitudes() const { return theLastAmplitudes; }
/**
* Access the last evaluated helicity amplitudes.
*/
AmplitudeMap& lastAmplitudes() { return theLastAmplitudes; }
/**
* Return last evaluated, leading colour helicity amplitudes.
*/
const AmplitudeMap& lastLargeNAmplitudes() const { return theLastLargeNAmplitudes; }
/**
* Access the last evaluated, leading colour helicity amplitudes.
*/
AmplitudeMap& lastLargeNAmplitudes() { return theLastLargeNAmplitudes; }
/**
* Return the matrix element squared.
*/
virtual double me2() const {
- return colourBasis()->me2(mePartonData(),lastAmplitudes());
+ return
+ lastCrossingSign()*colourBasis()->me2(mePartonData(),lastAmplitudes());
}
/**
* Return the colour correlated matrix element.
*/
virtual double colourCorrelatedME2(pair<int,int> ij) const;
/**
* Return the colour and spin correlated matrix element.
*/
virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
const SpinCorrelationTensor& c) const;
/**
* Evaluate the amplitude for the given colour tensor id and
* helicity assignment
*/
virtual Complex evaluate(size_t, const vector<int>&, Complex&) { return 0.; }
//@}
/** @name One-loop amplitudes */
//@{
/**
* Return true, if this amplitude is capable of calculating one-loop
* (QCD) corrections.
*/
virtual bool haveOneLoop() const { return false; }
/**
* Return true, if this amplitude only provides
* one-loop (QCD) corrections.
*/
virtual bool onlyOneLoop() const { return false; }
/**
* Return true, if one loop corrections have been calculated in
* dimensional reduction. Otherwise conventional dimensional
* regularization is assumed. Note that renormalization is always
* assumed to be MSbar.
*/
bool isDR() const { return false; }
/**
* Return true, if one loop corrections are given in the conventions
* of the integrated dipoles.
*/
bool isCS() const { return false; }
/**
* Return the value of the dimensional regularization
* parameter. Note that renormalization scale dependence is fully
* restored in DipoleIOperator.
*/
virtual Energy2 mu2() const { return 0.*GeV2; }
/**
* Calculate the one-loop amplitudes for the phasespace point
* stored in lastXComb, if provided.
*/
virtual void prepareOneLoopAmplitudes();
/**
* Return last evaluated one-loop helicity amplitudes.
*/
const AmplitudeMap& lastOneLoopAmplitudes() const { return theLastOneLoopAmplitudes; }
/**
* Access the last evaluated one-loop helicity amplitudes.
*/
AmplitudeMap& lastOneLoopAmplitudes() { return theLastOneLoopAmplitudes; }
/**
* Return the one-loop/tree interference.
*/
virtual double oneLoopInterference() const {
- return colourBasis()->interference(mePartonData(),
- lastOneLoopAmplitudes(),lastAmplitudes());
+ return
+ lastCrossingSign()*colourBasis()->interference(mePartonData(),
+ lastOneLoopAmplitudes(),lastAmplitudes());
}
/**
* Evaluate the amplitude for the given colour tensor id and
* helicity assignment
*/
virtual Complex evaluateOneLoop(size_t, const vector<int>&) { return 0.; }
//@}
/** @name Caching and helpers to setup amplitude objects. */
//@{
/**
* Flush all cashes.
*/
virtual void flushCaches() {
calculateTrees = true;
calculateLoops = true;
}
/**
* Clone this amplitude.
*/
Ptr<MatchboxAmplitude>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
virtual void cloneDependencies(const std::string& prefix = "");
//@}
/** @name Diagnostic information */
//@{
/**
* Dump xcomb hierarchies.
*/
void dumpInfo(const string& prefix = "") const;
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* Recursively generate helicities
*/
void doGenerateHelicities(set<vector<int> >& res,
vector<int>& current,
size_t pos) const;
/**
* The Herwig++ StandardModel object
*/
Ptr<StandardModel>::tcptr theStandardModel;
/**
* The number of light flavours to be used.
*/
unsigned int theNLight;
/**
* The colour basis implementation to be used.
*/
Ptr<ColourBasis>::ptr theColourBasis;
/**
* The dimensionality of the colour basis for the processes covered
* by the colour basis.
*/
size_t theColourBasisDim;
/**
* References to the amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastAmplitudes;
/**
* References to the leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastLargeNAmplitudes;
/**
* References to the one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastOneLoopAmplitudes;
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
map<tStdXCombPtr,vector<int> > theCrossingMap;
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<tStdXCombPtr,map<size_t,size_t> > theColourMap;
/**
* The crossing signs as filled by the last call to
* fillCrossingMap()
*/
map<tStdXCombPtr,double> theCrossingSigns;
/**
* The amplitude parton data.
*/
map<tStdXCombPtr,cPDVector> theAmplitudePartonData;
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
map<tStdXCombPtr,vector<int> >::iterator theLastCrossingMap;
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<tStdXCombPtr,map<size_t,size_t> >::iterator theLastColourMap;
/**
* The amplitude parton data.
*/
map<tStdXCombPtr,cPDVector>::iterator theLastAmplitudePartonData;
/**
* The crossing sign.
*/
double theLastCrossingSign;
/**
* True, if tree amplitudes need to be recalculated.
*/
bool calculateTrees;
/**
* True, if loop amplitudes need to be recalculated.
*/
bool calculateLoops;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxAmplitude & operator=(const MatchboxAmplitude &);
};
}
#endif /* HERWIG_MatchboxAmplitude_H */
diff --git a/MatrixElement/Matchbox/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc
--- a/MatrixElement/Matchbox/MatchboxFactory.cc
+++ b/MatrixElement/Matchbox/MatchboxFactory.cc
@@ -1,754 +1,780 @@
// -*- C++ -*-
//
// MatchboxFactory.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 MatchboxFactory class.
//
#include "MatchboxFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
using namespace Herwig;
using std::ostream_iterator;
MatchboxFactory::MatchboxFactory()
: SubProcessHandler(), theNLight(0),
theOrderInAlphaS(0), theOrderInAlphaEW(0),
theBornContributions(true), theVirtualContributions(true),
theRealContributions(true), theSubProcessGroups(false),
theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0),
theFixedCouplings(false), theVetoScales(false),
theVerbose(false), theSubtractionData("") {}
MatchboxFactory::~MatchboxFactory() {}
IBPtr MatchboxFactory::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxFactory::fullclone() const {
return new_ptr(*this);
}
void MatchboxFactory::prepareME(Ptr<MatchboxMEBase>::ptr me) const {
Ptr<MatchboxAmplitude>::ptr amp =
dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>((*me).amplitude());
me->matchboxAmplitude(amp);
if ( diagramGenerator() && !me->diagramGenerator() )
me->diagramGenerator(diagramGenerator());
if ( me->nLight() == 0 )
me->nLight(nLight());
if ( phasespace() && !me->phasespace() )
me->phasespace(phasespace());
if ( scaleChoice() && !me->scaleChoice() )
me->scaleChoice(scaleChoice());
if ( me->factorizationScaleFactor() == 1.0 )
me->factorizationScaleFactor(factorizationScaleFactor());
if ( me->renormalizationScaleFactor() == 1.0 )
me->renormalizationScaleFactor(renormalizationScaleFactor());
if ( fixedCouplings() )
me->setFixedCouplings();
if ( cache() && !me->cache() )
me->cache(cache());
if ( verbose() )
me->setVerbose();
}
struct ProjectQN {
inline pair<int,pair<int,int> > operator()(PDPtr p) const {
return
pair<int,pair<int,int> >((*p).iSpin(),pair<int,int>((*p).iCharge(),(*p).iColour()));
}
};
string pid(const vector<pair<int,pair<int,int> > >& key) {
ostringstream res;
for ( vector<pair<int,pair<int,int> > >::const_iterator k =
key.begin(); k != key.end(); ++k )
res << k->first << k->second.first
<< k->second.second;
return res.str();
}
vector<Ptr<MatchboxMEBase>::ptr> MatchboxFactory::
makeMEs(const vector<string>& proc, unsigned int orderas) const {
typedef vector<pair<int,pair<int,int> > > QNKey;
map<Ptr<MatchboxAmplitude>::ptr,map<QNKey,vector<PDVector> > > ampProcs;
set<PDVector> processes = makeSubProcesses(proc);
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
+ (**amp).orderInGs(orderas);
+ (**amp).orderInGem(orderInAlphaEW());
if ( (**amp).orderInGs() != orderas ||
(**amp).orderInGem() != orderInAlphaEW() )
continue;
for ( set<PDVector>::const_iterator p = processes.begin();
p != processes.end(); ++p ) {
if ( !(**amp).canHandle(*p) )
continue;
QNKey key;
transform(p->begin(),p->end(),back_inserter(key),ProjectQN());
ampProcs[*amp][key].push_back(*p);
}
}
vector<Ptr<MatchboxMEBase>::ptr> res;
for ( map<Ptr<MatchboxAmplitude>::ptr,map<QNKey,vector<PDVector> > >::const_iterator
ap = ampProcs.begin(); ap != ampProcs.end(); ++ap ) {
for ( map<QNKey,vector<PDVector> >::const_iterator m = ap->second.begin();
m != ap->second.end(); ++m ) {
- Ptr<MatchboxMEBase>::ptr me = new_ptr(MatchboxMEBase());
+ Ptr<MatchboxMEBase>::ptr me = ap->first->makeME(m->second);
me->subProcesses() = m->second;
me->amplitude(ap->first);
string pname = "ME" + ap->first->name() + pid(m->first);
if ( ! (generator()->preinitRegister(me,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
res.push_back(me);
}
}
return res;
}
void MatchboxFactory::setup() {
if ( !amplitudes().empty() ) {
if ( particleGroups().find("j") == particleGroups().end() )
throw InitException() << "Could not find a jet particle group named 'j'";
// rebind the particle data objects
for ( map<string,PDVector>::iterator g = particleGroups().begin();
g != particleGroups().end(); ++g )
for ( PDVector::iterator p = g->second.begin();
p != g->second.end(); ++p ) {
long checkid = (**p).id();
*p = getParticleData((**p).id());
assert((**p).id() == checkid);
}
nLight(particleGroups()["j"].size());
- bornMEs() = makeMEs(process,orderInAlphaS());
+ vector<Ptr<MatchboxMEBase>::ptr> ames = makeMEs(process,orderInAlphaS());
+ copy(ames.begin(),ames.end(),back_inserter(bornMEs()));
if ( realContributions() ) {
vector<string> rproc = process;
rproc.push_back("j");
- realEmissionMEs() = makeMEs(rproc,orderInAlphaS()+1);
+ ames = makeMEs(rproc,orderInAlphaS()+1);
+ copy(ames.begin(),ames.end(),back_inserter(realEmissionMEs()));
}
}
// check if we have virtual contributions
bool haveVirtuals = true;
// check DR conventions of virtual contributions
bool virtualsAreDR = false;
bool virtualsAreCDR = false;
// check finite term conventions of virtual contributions
bool virtualsAreCS = false;
bool virtualsAreStandard = false;
// check and prepare the Born and virtual matrix elements
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
prepareME(*born);
haveVirtuals &= (**born).haveOneLoop();
if ( (**born).haveOneLoop() ) {
virtualsAreDR |= (**born).isDR();
virtualsAreCDR |= !(**born).isDR();
virtualsAreCS |= (**born).isCS();
virtualsAreStandard |= !(**born).isCS();
}
}
// check the additional insertion operators
if ( !virtuals().empty() )
haveVirtuals = true;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
virtualsAreDR |= (**virt).isDR();
virtualsAreCDR |= !(**virt).isDR();
virtualsAreCS |= (**virt).isCS();
virtualsAreStandard |= !(**virt).isCS();
}
// check for consistent conventions on virtuals, if we are to include them
if ( virtualContributions() ) {
if ( virtualsAreDR && virtualsAreCDR ) {
throw InitException() << "Virtual corrections use inconsistent regularization schemes.\n";
}
if ( virtualsAreCS && virtualsAreStandard ) {
throw InitException() << "Virtual corrections use inconsistent conventions on finite terms.\n";
}
if ( !haveVirtuals ) {
throw InitException() << "Could not find amplitudes for all virtual contributions needed.\n";
}
}
// prepare dipole insertion operators
if ( virtualContributions() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( virtualsAreDR )
(**virt).useDR();
else
(**virt).useCDR();
if ( virtualsAreCS )
(**virt).useCS();
else
(**virt).useNonCS();
}
}
// prepare the real emission matrix elements
if ( realContributions() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
prepareME(*real);
}
}
// start creating matrix elements
MEs().clear();
// setup born and virtual contributions
if ( !bornContributions() && virtualContributions() ) {
throw InitException() << "Virtual corrections without Born contributions not yet supported.\n";
}
if ( bornContributions() && !virtualContributions() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
if ( (**born).onlyOneLoop() )
continue;
Ptr<MatchboxMEBase>::ptr bornme = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( ! (generator()->preinitRegister(bornme,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
bornme->cloneDependencies();
MEs().push_back(bornme);
}
}
if ( bornContributions() && virtualContributions() ) {
bornVirtualMEs().clear();
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
Ptr<MatchboxNLOME>::ptr nlo = new_ptr(MatchboxNLOME());
string pname = fullName() + "/" + (**born).name();
if ( ! (generator()->preinitRegister(nlo,pname) ) )
throw InitException() << "NLO ME " << pname << " already existing.";
nlo->matrixElement(*born);
nlo->virtuals().clear();
if ( !nlo->matrixElement()->onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
if ( nlo->virtuals().empty() )
throw InitException() << "No insertion operators have been found for "
<< (**born).name() << "\n";
}
nlo->cloneDependencies();
bornVirtualMEs().push_back(nlo);
MEs().push_back(nlo);
}
}
if ( realContributions() ) {
if ( theSubtractionData != "" )
if ( theSubtractionData[theSubtractionData.size()-1] != '/' )
theSubtractionData += "/";
subtractedMEs().clear();
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
Ptr<SubtractedME>::ptr sub = new_ptr(SubtractedME());
string pname = fullName() + "/" + (**real).name();
if ( ! (generator()->preinitRegister(sub,pname) ) )
throw InitException() << "Subtracted ME " << pname << " already existing.";
sub->borns() = bornMEs();
sub->head(*real);
sub->allDipoles().clear();
sub->dependent().clear();
sub->getDipoles();
if ( verbose() )
sub->setVerbose();
if ( subProcessGroups() )
sub->setSubProcessGroups();
if ( vetoScales() )
sub->doVetoScales();
if ( subtractionData() != "" )
sub->subtractionData(subtractionData());
subtractedMEs().push_back(sub);
MEs().push_back(sub);
}
}
}
void MatchboxFactory::print(ostream& os) const {
os << "--- MatchboxFactory setup -----------------------------------------------------------\n";
if ( !amplitudes().empty() ) {
os << " generated Born matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = bornMEs().begin();
m != bornMEs().end(); ++m ) {
- os << " '" << (**m).name() << "'\n";
+ os << " '" << (**m).name() << "' for subprocesses:\n";
+ for ( vector<PDVector>::const_iterator p = (**m).subProcesses().begin();
+ p != (**m).subProcesses().end(); ++p ) {
+ os << " ";
+ for ( PDVector::const_iterator pp = p->begin();
+ pp != p->end(); ++pp ) {
+ os << (**pp).PDGName() << " ";
+ if ( pp == p->begin() + 1 )
+ os << "-> ";
+ }
+ os << "\n";
+ }
}
os << flush;
os << " generated real emission matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = realEmissionMEs().begin();
m != realEmissionMEs().end(); ++m ) {
- os << " '" << (**m).name() << "'\n";
+ os << " '" << (**m).name() << "' for subprocesses:\n";
+ for ( vector<PDVector>::const_iterator p = (**m).subProcesses().begin();
+ p != (**m).subProcesses().end(); ++p ) {
+ os << " ";
+ for ( PDVector::const_iterator pp = p->begin();
+ pp != p->end(); ++pp ) {
+ os << (**pp).PDGName() << " ";
+ if ( pp == p->begin() + 1 )
+ os << "-> ";
+ }
+ os << "\n";
+ }
}
os << flush;
}
os << " generated Born+virtual matrix elements:\n";
for ( vector<Ptr<MatchboxNLOME>::ptr>::const_iterator bv
= bornVirtualMEs().begin(); bv != bornVirtualMEs().end(); ++bv ) {
(**bv).print(os);
}
os << " generated subtracted matrix elements:\n";
for ( vector<Ptr<SubtractedME>::ptr>::const_iterator sub
= subtractedMEs().begin(); sub != subtractedMEs().end(); ++sub ) {
os << " '" << (**sub).name() << "'\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxFactory::doinit() {
setup();
if ( theVerbose )
print(Repository::clog());
SubProcessHandler::doinit();
}
void MatchboxFactory::persistentOutput(PersistentOStream & os) const {
os << theDiagramGenerator
<< theNLight << theOrderInAlphaS << theOrderInAlphaEW
<< theBornContributions << theVirtualContributions
<< theRealContributions << theSubProcessGroups
<< thePhasespace << theScaleChoice
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theFixedCouplings << theVetoScales
<< theAmplitudes << theCache
<< theBornMEs << theVirtuals << theRealEmissionMEs
<< theBornVirtualMEs << theSubtractedMEs
<< theVerbose << theSubtractionData
<< theParticleGroups << process;
}
void MatchboxFactory::persistentInput(PersistentIStream & is, int) {
is >> theDiagramGenerator
>> theNLight >> theOrderInAlphaS >> theOrderInAlphaEW
>> theBornContributions >> theVirtualContributions
>> theRealContributions >> theSubProcessGroups
>> thePhasespace >> theScaleChoice
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theFixedCouplings >> theVetoScales
>> theAmplitudes >> theCache
>> theBornMEs >> theVirtuals >> theRealEmissionMEs
>> theBornVirtualMEs >> theSubtractedMEs
>> theVerbose >> theSubtractionData
>> theParticleGroups >> process;
}
string MatchboxFactory::startParticleGroup(string name) {
particleGroupName = StringUtils::stripws(name);
particleGroup.clear();
return "";
}
string MatchboxFactory::endParticleGroup(string) {
if ( particleGroup.empty() )
throw InitException() << "Empty particle group.";
particleGroups()[particleGroupName] = particleGroup;
particleGroup.clear();
return "";
}
string MatchboxFactory::doProcess(string in) {
process = StringUtils::split(in);
if ( process.size() < 3 )
throw InitException() << "Invalid process.";
for ( vector<string>::iterator p = process.begin();
p != process.end(); ++p ) {
*p = StringUtils::stripws(*p);
}
return "";
}
struct SortPID {
inline bool operator()(PDPtr a, PDPtr b) const {
return a->id() < b->id();
}
};
set<PDVector> MatchboxFactory::
makeSubProcesses(const vector<string>& proc) const {
if ( proc.empty() )
throw InitException() << "No process specified.";
vector<PDVector> allProcs(1);
size_t pos = 0;
typedef map<string,PDVector>::const_iterator GroupIterator;
while ( pos < proc.size() ) {
GroupIterator git =
particleGroups().find(proc[pos]);
if ( git == particleGroups().end() ) {
throw InitException() << "particle group '"
<< proc[pos] << "' not defined.";
}
vector<PDVector> mine;
for ( vector<PDVector>::const_iterator i = allProcs.begin();
i != allProcs.end(); ++i ) {
for ( PDVector::const_iterator p = git->second.begin();
p != git->second.end(); ++p ) {
PDVector v = *i;
v.push_back(*p);
mine.push_back(v);
}
}
allProcs = mine;
++pos;
}
set<PDVector> allCheckedProcs;
for ( vector<PDVector>::const_iterator p = allProcs.begin();
p != allProcs.end(); ++p ) {
int charge = -(*p)[0]->iCharge() -(*p)[1]->iCharge();
for ( size_t k = 2; k < (*p).size(); ++k )
charge += (*p)[k]->iCharge();
if ( charge != 0 )
continue;
PDVector pr = *p;
sort(pr.begin()+2,pr.end(),SortPID());
allCheckedProcs.insert(pr);
}
return allCheckedProcs;
}
void MatchboxFactory::Init() {
static ClassDocumentation<MatchboxFactory> documentation
("MatchboxFactory",
"NLO QCD corrections have been calculated "
"using Matchbox \\cite{Platzer:2011bc}",
"%\\cite{Platzer:2011bc}\n"
"\\bibitem{Platzer:2011bc}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Dipole Showers and Automated NLO Matching in Herwig++,''\n"
"arXiv:1109.6256 [hep-ph].\n"
"%%CITATION = ARXIV:1109.6256;%%");
static Reference<MatchboxFactory,Tree2toNGenerator> interfaceDiagramGenerator
("DiagramGenerator",
"Set the diagram generator.",
&MatchboxFactory::theDiagramGenerator, false, false, true, true, false);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaS
("OrderInAlphaS",
"The order in alpha_s to consider.",
&MatchboxFactory::theOrderInAlphaS, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaEW
("OrderInAlphaEW",
"The order in alpha_EW",
&MatchboxFactory::theOrderInAlphaEW, 2, 0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceBornContributions
("BornContributions",
"Switch on or off the Born contributions.",
&MatchboxFactory::theBornContributions, true, false, false);
static SwitchOption interfaceBornContributionsOn
(interfaceBornContributions,
"On",
"Switch on Born contributions.",
true);
static SwitchOption interfaceBornContributionsOff
(interfaceBornContributions,
"Off",
"Switch off Born contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceVirtualContributions
("VirtualContributions",
"Switch on or off the virtual contributions.",
&MatchboxFactory::theVirtualContributions, true, false, false);
static SwitchOption interfaceVirtualContributionsOn
(interfaceVirtualContributions,
"On",
"Switch on virtual contributions.",
true);
static SwitchOption interfaceVirtualContributionsOff
(interfaceVirtualContributions,
"Off",
"Switch off virtual contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceRealContributions
("RealContributions",
"Switch on or off the real contributions.",
&MatchboxFactory::theRealContributions, true, false, false);
static SwitchOption interfaceRealContributionsOn
(interfaceRealContributions,
"On",
"Switch on real contributions.",
true);
static SwitchOption interfaceRealContributionsOff
(interfaceRealContributions,
"Off",
"Switch off real contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceSubProcessGroups
("SubProcessGroups",
"Switch on or off production of sub-process groups.",
&MatchboxFactory::theSubProcessGroups, false, false, false);
static SwitchOption interfaceSubProcessGroupsOn
(interfaceSubProcessGroups,
"On",
"On",
true);
static SwitchOption interfaceSubProcessGroupsOff
(interfaceSubProcessGroups,
"Off",
"Off",
false);
static Reference<MatchboxFactory,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"Set the phasespace generator.",
&MatchboxFactory::thePhasespace, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxScaleChoice> interfaceScaleChoice
("ScaleChoice",
"Set the scale choice object.",
&MatchboxFactory::theScaleChoice, false, false, true, true, false);
static Parameter<MatchboxFactory,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&MatchboxFactory::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&MatchboxFactory::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceFixedCouplings
("FixedCouplings",
"Switch on or off fixed couplings.",
&MatchboxFactory::theFixedCouplings, true, false, false);
static SwitchOption interfaceFixedCouplingsOn
(interfaceFixedCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedCouplingsOff
(interfaceFixedCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceVetoScales
("VetoScales",
"Switch on or setting veto scales.",
&MatchboxFactory::theVetoScales, true, false, false);
static SwitchOption interfaceVetoScalesOn
(interfaceVetoScales,
"On",
"On",
true);
static SwitchOption interfaceVetoScalesOff
(interfaceVetoScales,
"Off",
"Off",
false);
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceAmplitudes
("Amplitudes",
"The amplitude objects.",
&MatchboxFactory::theAmplitudes, -1, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxMECache> interfaceCache
("Cache",
"Set the matrix element cache object.",
&MatchboxFactory::theCache, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceBornMEs
("BornMEs",
"The Born matrix elements to be used",
&MatchboxFactory::theBornMEs, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxInsertionOperator> interfaceVirtuals
("Virtuals",
"The virtual corrections to include",
&MatchboxFactory::theVirtuals, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceRealEmissionMEs
("RealEmissionMEs",
"The RealEmission matrix elements to be used",
&MatchboxFactory::theRealEmissionMEs, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxNLOME> interfaceBornVirtuals
("BornVirtualMEs",
"The generated Born/virtual contributions",
&MatchboxFactory::theBornVirtualMEs, -1, false, true, true, true, false);
static RefVector<MatchboxFactory,SubtractedME> interfaceSubtractedMEs
("SubtractedMEs",
"The generated Born/virtual contributions",
&MatchboxFactory::theSubtractedMEs, -1, false, true, true, true, false);
static Switch<MatchboxFactory,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&MatchboxFactory::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"On",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"Off",
false);
static Parameter<MatchboxFactory,string> interfaceSubtractionData
("SubtractionData",
"Prefix for subtraction check data.",
&MatchboxFactory::theSubtractionData, "",
false, false);
static RefVector<MatchboxFactory,ParticleData> interfaceParticleGroup
("ParticleGroup",
"The particle group just started.",
&MatchboxFactory::particleGroup, -1, false, false, true, false, false);
static Command<MatchboxFactory> interfaceStartParticleGroup
("StartParticleGroup",
"Start a particle group.",
&MatchboxFactory::startParticleGroup, false);
static Command<MatchboxFactory> interfaceEndParticleGroup
("EndParticleGroup",
"End a particle group.",
&MatchboxFactory::endParticleGroup, false);
static Command<MatchboxFactory> interfaceProcess
("Process",
"Set the process to consider.",
&MatchboxFactory::doProcess, 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<MatchboxFactory,SubProcessHandler>
describeHerwigMatchboxFactory("Herwig::MatchboxFactory", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.cc
@@ -1,132 +1,132 @@
// -*- C++ -*-
//
// FILightInvertedTildeKinematics.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 FILightInvertedTildeKinematics class.
//
#include "FILightInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FILightInvertedTildeKinematics::FILightInvertedTildeKinematics() {}
FILightInvertedTildeKinematics::~FILightInvertedTildeKinematics() {}
IBPtr FILightInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FILightInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool FILightInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
double y = sqr(pt/lastScale())/(z*(1.-z));
double x = 1./(1.+y);
if ( x < spectatorX() ) {
jacobian(0.0);
return false;
}
mapping /= z*(1.-z);
jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
double phi = 2.*Constants::pi*r[2];
Lorentz5Momentum kt
- = getKt(emitter,spectator,pt,phi);
+ = getKt(spectator,emitter,pt,phi,true);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = z;
realEmitterMomentum() = z*emitter + y*(1.-z)*spectator + kt;
realEmissionMomentum() = (1.-z)*emitter + y*z*spectator - kt;
realSpectatorMomentum() = (1.+y)*spectator;
realEmitterMomentum().setMass(ZERO);
realEmitterMomentum().rescaleEnergy();
realEmissionMomentum().setMass(ZERO);
realEmissionMomentum().rescaleEnergy();
realSpectatorMomentum().setMass(ZERO);
realSpectatorMomentum().rescaleEnergy();
return true;
}
Energy FILightInvertedTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
return scale * sqrt(z*(1.-z)*(1.-x)/x);
}
Energy FILightInvertedTildeKinematics::ptMax() const {
double x = spectatorX();
return sqrt((1.-x)/x)*lastScale()/2.;
}
pair<double,double> FILightInvertedTildeKinematics::zBounds(Energy pt) const {
double s = sqrt(1.-sqr(pt/ptMax()));
return make_pair(0.5*(1.-s),0.5*(1.+s));
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FILightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void FILightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void FILightInvertedTildeKinematics::Init() {
static ClassDocumentation<FILightInvertedTildeKinematics> documentation
("FILightInvertedTildeKinematics inverts the final-initial tilde "
"kinematics.");
}
// *** 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<FILightInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigFILightInvertedTildeKinematics("Herwig::FILightInvertedTildeKinematics", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.cc
@@ -1,137 +1,137 @@
// -*- C++ -*-
//
// IFLightInvertedTildeKinematics.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 IFLightInvertedTildeKinematics class.
//
#include "IFLightInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IFLightInvertedTildeKinematics::IFLightInvertedTildeKinematics() {}
IFLightInvertedTildeKinematics::~IFLightInvertedTildeKinematics() {}
IBPtr IFLightInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IFLightInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool IFLightInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
double ratio = sqr(pt/lastScale());
double x = ( z*(1.-z) - ratio ) / ( 1. - z - ratio );
double u = ratio / (1.-z);
pt = lastScale()*sqrt(u*(1.-u)*(1.-x)/x);
if ( x < emitterX() || x > 1. || u > 1. ) {
jacobian(0.0);
return false;
}
mapping /= sqr(z*(1.-z)-ratio)/(1.-z-ratio);
jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
double phi = 2.*Constants::pi*r[2];
- Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi);
+ Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi,true);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = u;
realEmitterMomentum() = (1./x)*emitter;
realEmissionMomentum() = ((1.-x)*(1.-u)/x)*emitter + u*spectator + kt;
realSpectatorMomentum() = ((1.-x)*u/x)*emitter + (1.-u)*spectator - kt;
realEmitterMomentum().setMass(ZERO);
realEmitterMomentum().rescaleEnergy();
realEmissionMomentum().setMass(ZERO);
realEmissionMomentum().rescaleEnergy();
realSpectatorMomentum().setMass(ZERO);
realSpectatorMomentum().rescaleEnergy();
return true;
}
Energy IFLightInvertedTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
return scale * sqrt(u*(1.-u)*(1.-x));
}
Energy IFLightInvertedTildeKinematics::ptMax() const {
double x = emitterX();
return sqrt(1.-x)*lastScale()/2.;
}
pair<double,double> IFLightInvertedTildeKinematics::zBounds(Energy pt) const {
double s = sqrt(1.-sqr(pt/ptMax()));
double x = emitterX();
return make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void IFLightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
}
void IFLightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
}
void IFLightInvertedTildeKinematics::Init() {
static ClassDocumentation<IFLightInvertedTildeKinematics> documentation
("IFLightInvertedTildeKinematics inverts the initial-final tilde "
"kinematics.");
}
// *** 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<IFLightInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigIFLightInvertedTildeKinematics("Herwig::IFLightInvertedTildeKinematics", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.cc
@@ -1,173 +1,199 @@
// -*- C++ -*-
//
// InvertedTildeKinematics.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 InvertedTildeKinematics class.
//
#include <limits>
#include "InvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
using namespace Herwig;
InvertedTildeKinematics::InvertedTildeKinematics()
: HandlerBase(), theJacobian(0.0), thePtCut(0.0*GeV) {}
InvertedTildeKinematics::~InvertedTildeKinematics() {}
void InvertedTildeKinematics::dumpInfo(const string& prefix) const {
generator()->log() << prefix << fullName()
<< " [" << this << "]\n";
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
-Lorentz5Momentum InvertedTildeKinematics::getKt (const Lorentz5Momentum& p1,
- const Lorentz5Momentum& p2,
- Energy pt,
- double phi) const {
+Lorentz5Momentum InvertedTildeKinematics::getKt(const Lorentz5Momentum& p1,
+ const Lorentz5Momentum& p2,
+ Energy pt,
+ double phi,
+ bool spacelike) const {
- Boost beta = (p1+p2).findBoostToCM();
-
- Lorentz5Momentum p1c = p1;
+ Lorentz5Momentum P;
+ if ( !spacelike )
+ P = p1 + p2;
+ else
+ P = p1 - p2;
- if (beta.mag2() > Constants::epsilon) {
- p1c.boost(beta);
+ Energy2 Q2 = abs(P.m2());
+
+ Lorentz5Momentum Q =
+ !spacelike ?
+ Lorentz5Momentum(ZERO,ZERO,ZERO,sqrt(Q2),sqrt(Q2)) :
+ Lorentz5Momentum(ZERO,ZERO,sqrt(Q2),ZERO,-sqrt(Q2));
+
+ if ( spacelike && Q.z() < P.z() )
+ Q.setZ(-Q.z());
+
+ bool boost =
+ abs((P-Q).vect().mag2()/GeV2) > 1e-10 ||
+ abs((P-Q).t()/GeV) > 1e-5;
+
+ Lorentz5Momentum inFrame1;
+ if ( boost )
+ inFrame1 = p1 + ((P*p1-Q*p1)/(P*Q-Q.mass2()))*(P-Q);
+ else
+ inFrame1 = p1;
+
+ Energy ptx = inFrame1.x();
+ Energy pty = inFrame1.y();
+ Energy q = 2.*inFrame1.z();
+
+ Energy Qp = sqrt(4.*(sqr(ptx)+sqr(pty))+sqr(q));
+ Energy Qy = sqrt(4.*sqr(pty)+sqr(q));
+
+ double cPhi = cos(phi);
+ double sPhi = sqrt(1.-sqr(cPhi));
+ if ( phi > Constants::pi )
+ sPhi = -sPhi;
+
+ Lorentz5Momentum kt;
+
+ if ( !spacelike ) {
+ kt.setT(ZERO);
+ kt.setX(pt*Qy*cPhi/Qp);
+ kt.setY(-pt*(4*ptx*pty*cPhi/Qp+q*sPhi)/Qy);
+ kt.setZ(2.*pt*(-ptx*q*cPhi/Qp + pty*sPhi)/Qy);
+ } else {
+ kt.setT(2.*pt*(ptx*q*cPhi+pty*Qp*sPhi)/(q*Qy));
+ kt.setX(pt*(Qp*q*cPhi+4.*ptx*pty*sPhi)/(q*Qy));
+ kt.setY(pt*Qy*sPhi/q);
+ kt.setZ(ZERO);
}
-
- Lorentz5Momentum k (0.*GeV,0.*GeV,0.*GeV,0.*GeV);
-
- double ct = p1c.vect().unit().z();
- double st = sqrt(1.-ct*ct);
-
- double cphi = cos(phi);
- double sphi = sqrt(1.-cphi*cphi);
- if (phi > Constants::pi) sphi = -sphi;
-
- if (st > Constants::epsilon) {
- double cchi = p1c.vect().unit().x()/st;
- double schi = p1c.vect().unit().y()/st;
- k.setX((cphi*cchi*ct-sphi*schi)*pt);
- k.setY((cphi*schi*ct+sphi*cchi)*pt);
- k.setZ(-cphi*st*pt);
- } else {
- k.setX(pt*cphi);
- k.setY(pt*sphi);
- k.setZ(0.*GeV);
- }
-
- if (beta.mag2() > Constants::epsilon)
- k.boost(-beta);
-
- return k;
+
+ if ( boost )
+ kt = kt + ((P*kt-Q*kt)/(P*Q-Q.mass2()))*(P-Q);
+ kt.setMass(-pt);
+ kt.rescaleRho();
+
+ return kt;
}
Energy InvertedTildeKinematics::lastScale() const {
if ( ( theDipole->bornEmitter() < 2 && theDipole->bornSpectator() > 1 ) ||
( theDipole->bornEmitter() > 1 && theDipole->bornSpectator() < 2 ) ) {
return -(bornEmitterMomentum()-bornSpectatorMomentum()).m();
}
return (bornEmitterMomentum()+bornSpectatorMomentum()).m();
}
pair<Energy,double> InvertedTildeKinematics::generatePtZ(double& jac, const double * r) const {
double kappaMin =
ptCut() != ZERO ?
sqr(ptCut()/ptMax()) :
sqr(0.1*GeV/GeV);
double kappa;
using namespace RandomHelpers;
if ( ptCut() > ZERO ) {
pair<double,double> kw =
generate(inverse(0.,kappaMin,1.),r[0]);
kappa = kw.first;
jac *= kw.second;
} else {
pair<double,double> kw =
generate((piecewise(),
flat(1e-4,kappaMin),
match(inverse(0.,kappaMin,1.))),r[0]);
kappa = kw.first;
jac *= kw.second;
}
Energy pt = sqrt(kappa)*ptMax();
pair<double,double> zLims = zBounds(pt);
pair<double,double> zw =
generate(inverse(0.,zLims.first,zLims.second)+
inverse(1.,zLims.first,zLims.second),r[1]);
double z = zw.first;
jac *= zw.second;
jac *= sqr(ptMax()/lastScale());
return make_pair(pt,z);
}
void InvertedTildeKinematics::rebind(const TranslationMap & trans) {
theDipole = trans.translate(theDipole);
HandlerBase::rebind(trans);
}
IVector InvertedTildeKinematics::getReferences() {
IVector ret = HandlerBase::getReferences();
ret.push_back(theDipole);
return ret;
}
void InvertedTildeKinematics::persistentOutput(PersistentOStream & os) const {
os << theDipole << theRealXComb << theBornXComb
<< ounit(theRealEmitterMomentum,GeV) << ounit(theRealEmissionMomentum,GeV)
<< ounit(theRealSpectatorMomentum,GeV) << theJacobian
<< ounit(thePtCut,GeV);
}
void InvertedTildeKinematics::persistentInput(PersistentIStream & is, int) {
is >> theDipole >> theRealXComb >> theBornXComb
>> iunit(theRealEmitterMomentum,GeV) >> iunit(theRealEmissionMomentum,GeV)
>> iunit(theRealSpectatorMomentum,GeV) >> theJacobian
>> iunit(thePtCut,GeV);
}
void InvertedTildeKinematics::Init() {
static ClassDocumentation<InvertedTildeKinematics> documentation
("InvertedTildeKinematics is the base class for the inverted 'tilde' "
"kinematics being used for subtraction terms in the "
"formalism of Catani and Seymour.");
}
// *** 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<InvertedTildeKinematics,HandlerBase>
describeInvertedTildeKinematics("Herwig::InvertedTildeKinematics", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h b/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h
--- a/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h
+++ b/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h
@@ -1,437 +1,438 @@
// -*- C++ -*-
//
// InvertedTildeKinematics.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_InvertedTildeKinematics_H
#define HERWIG_InvertedTildeKinematics_H
//
// This is the declaration of the InvertedTildeKinematics class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief InvertedTildeKinematics is the base class for the inverted 'tilde'
* kinematics being used for subtraction terms in the
* formalism of Catani and Seymour.
*
*/
class InvertedTildeKinematics: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
InvertedTildeKinematics();
/**
* The destructor.
*/
virtual ~InvertedTildeKinematics();
//@}
public:
/** @name Access to kinematic quantities. */
//@{
/**
* Return the momentum of the emitter in the real emission process
*/
const Lorentz5Momentum& realEmitterMomentum() const { return theRealEmitterMomentum; }
/**
* Return the momentum of the emission in the real emission process
*/
const Lorentz5Momentum& realEmissionMomentum() const { return theRealEmissionMomentum; }
/**
* Return the momentum of the spectator in the real emission process
*/
const Lorentz5Momentum& realSpectatorMomentum() const { return theRealSpectatorMomentum; }
/**
* Return the momentum of the emitter in the underlying Born process
*/
const Lorentz5Momentum& bornEmitterMomentum() const {
return theBornXComb->meMomenta()[theDipole->bornEmitter()];
}
/**
* Return the momentum of the spectator in the underlying Born process
*/
const Lorentz5Momentum& bornSpectatorMomentum() const {
return theBornXComb->meMomenta()[theDipole->bornSpectator()];
}
/**
* Return the momentum fraction of the emitter
*/
double emitterX() const {
return
theDipole->bornEmitter() == 0 ?
theBornXComb->lastX1() :
theBornXComb->lastX2();
}
/**
* Return the momentum fraction of the spectator
*/
double spectatorX() const {
return
theDipole->bornSpectator() == 0 ?
theBornXComb->lastX1() :
theBornXComb->lastX2();
}
/**
* Return the vector of dimensionless variables calculated
*/
const vector<double>& subtractionParameters() const { return theDipole->subtractionParameters(); }
/**
* Return true, if this InvertedTildeKinematics object needs to transform
* all other particles in the process except the emitter, emission and spectator
*/
virtual bool doesTransform() const { return false; }
/**
* If this InvertedTildeKinematics object needs to transform all other particles
* in the process except the emitter, emission and spectator, return the transformed
* momentum.
*/
virtual Lorentz5Momentum transform(const Lorentz5Momentum& p) const { return p; }
/**
* Return the centre of mass energy for the underlying Born configuration
*/
Energy2 sHat() const { return theBornXComb->lastSHat(); }
//@}
public:
/**
* Clone this object
*/
Ptr<InvertedTildeKinematics>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<InvertedTildeKinematics>::ptr>(clone());
}
/**
* Dump xcomb hierarchies.
*/
void dumpInfo(const string& prefix = "") const;
/** @name Access to process data. */
//@{
/**
* Prepare given a dipole, and XCombs describing the real emission
* and underlying Born processes, respectively.
*/
void prepare(tcStdXCombPtr newRealXComb,
tcStdXCombPtr newBornXComb) {
theRealXComb = newRealXComb; theBornXComb = newBornXComb;
}
/**
* Set the current dipole
*/
void dipole(Ptr<SubtractionDipole>::tptr dip) { theDipole = dip; }
/**
* Return the current dipole
*/
Ptr<SubtractionDipole>::tptr dipole() { return theDipole; }
/**
* Return the current dipole
*/
Ptr<SubtractionDipole>::tcptr dipole() const { return theDipole; }
/**
* Return the number of random numbers needed to generate
* a real emission configuration off the underlying Born
* configuration.
*/
virtual int nDimRadiation() const { return 3; }
/**
* Perform the mapping of the tilde kinematics for the
* last selected process and store all dimensionless
* variables in the subtractionParameters() vector.
* Return false, if the calculation of the real
* kinematics was impossible for the selected configuration
* and true on success.
*/
virtual bool doMap(const double *) = 0;
/**
* Set an optional cutoff on the emission's
* transverse momentum.
*/
void ptCut(Energy pt) { thePtCut = pt; }
/**
* Return the optional cutoff on the emission's
* transverse momentum.
*/
Energy ptCut() const { return thePtCut; }
/**
* Return the random number index
* corresponding to the evolution variable.
*/
virtual int evolutionVariable() const { return 0; }
/**
* Return the cutoff on the evolution
* random number corresponding to the pt cut.
*/
virtual double evolutionCutoff() const { return 0.0; }
/**
* Return the pt associated to the last generated splitting.
*/
virtual Energy lastPt() const = 0;
/**
* Return the relevant dipole scale
*/
virtual Energy lastScale() const;
/**
* Return the upper bound on pt
*/
virtual Energy ptMax() const = 0;
/**
* Given a pt, return the boundaries on z
*/
virtual pair<double,double> zBounds(Energy pt) const = 0;
/**
* Generate pt and z
*/
virtual pair<Energy,double> generatePtZ(double& jac, const double * r) const;
/**
* Return the single particle phasespace weight in units
* of sHat() for the last selected configuration.
*/
double jacobian() const { return theJacobian; }
/**
* Return the particle type of the emitter in the real emission process
*/
cPDPtr realEmitterData() const {
return
(theDipole && theRealXComb) ?
theRealXComb->mePartonData()[theDipole->realEmitter()] :
cPDPtr();
}
/**
* Return the particle type of the emission in the real emission process
*/
cPDPtr realEmissionData() const {
return
(theDipole && theRealXComb) ?
theRealXComb->mePartonData()[theDipole->realEmission()] :
cPDPtr();
}
/**
* Return the particle type of the spectator in the real emission process
*/
cPDPtr realSpectatorData() const {
return
(theDipole && theRealXComb) ?
theRealXComb->mePartonData()[theDipole->realSpectator()] :
cPDPtr();
}
/**
* Return the particle type of the emitter in the underlying Born process
*/
cPDPtr bornEmitterData() const {
return
(theDipole && theBornXComb) ?
theBornXComb->mePartonData()[theDipole->bornEmitter()] :
cPDPtr();
}
/**
* Return the particle type of the spectator in the underlying Born process
*/
cPDPtr bornSpectatorData() const {
return
(theDipole && theBornXComb) ?
theBornXComb->mePartonData()[theDipole->bornSpectator()] :
cPDPtr();
}
//@}
protected:
/**
* Access the momentum of the emitter in the real emission process
*/
Lorentz5Momentum& realEmitterMomentum() { return theRealEmitterMomentum; }
/**
* Access the momentum of the emission in the real emission process
*/
Lorentz5Momentum& realEmissionMomentum() { return theRealEmissionMomentum; }
/**
* Access the momentum of the spectator in the real emission process
*/
Lorentz5Momentum& realSpectatorMomentum() { return theRealSpectatorMomentum; }
/**
* Access the vector of dimensionless variables calculated
*/
vector<double>& subtractionParameters() { return theDipole->subtractionParameters(); }
/**
* Set the single particle phasespace weight in units
* of sHat() for the last selected configuration.
*/
void jacobian(double w) { theJacobian = w; }
/**
* Calculate a transverse momentum for the given momenta,
* invariant pt and azimuth.
*/
- Lorentz5Momentum getKt (const Lorentz5Momentum& p1,
- const Lorentz5Momentum& p2,
- Energy pt,
- double phi) const;
+ Lorentz5Momentum getKt(const Lorentz5Momentum& p1,
+ const Lorentz5Momentum& p2,
+ Energy pt,
+ double phi,
+ bool spacelike = false) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans);
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
//@}
private:
/**
* The last dipole this InvertedTildeKinematics has been selected for
*/
Ptr<SubtractionDipole>::tptr theDipole;
/**
* The XComb object describing the real emission process
*/
tcStdXCombPtr theRealXComb;
/**
* The XComb object describing the underlying Born process
*/
tcStdXCombPtr theBornXComb;
/**
* The momentum of the emitter in the real emission process
*/
Lorentz5Momentum theRealEmitterMomentum;
/**
* The momentum of the emission in the real emission process
*/
Lorentz5Momentum theRealEmissionMomentum;
/**
* The momentum of the spectator in the real emission process
*/
Lorentz5Momentum theRealSpectatorMomentum;
/**
* Return the single particle phasespace weight in units
* of sHat() for the last selected configuration.
*/
double theJacobian;
/**
* The optional cutoff on the emission's
* transverse momentum.
*/
Energy thePtCut;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
InvertedTildeKinematics & operator=(const InvertedTildeKinematics &);
};
}
#endif /* HERWIG_InvertedTildeKinematics_H */
diff --git a/src/DIS-Matchbox.in b/src/DIS-Matchbox.in
new file mode 100644
--- /dev/null
+++ b/src/DIS-Matchbox.in
@@ -0,0 +1,91 @@
+read Matchbox.in
+
+################################################################################
+# setup the sampler
+################################################################################
+
+set /Herwig/EventHandlers/DISHandler:Sampler /Herwig/MatrixElements/Matchbox/Samplers/Sampler
+
+################################################################################
+# setup the shower
+#
+# use -LO.in or -NLO.in for LO/NLO simulations, respectively.
+#
+################################################################################
+
+read DipoleShower.in
+# read DipoleShowerParameters-LO.in
+# read DipoleShowerParameters-NLO.in
+
+################################################################################
+# setup matrix element and matching
+#
+# uncomment the desired calculation
+#
+################################################################################
+
+cd /Herwig/EventHandlers
+
+set DISHandler:BeamA /Herwig/Particles/e+
+set DISLuminosity:BeamEMaxA 27.5.*GeV
+set DISHandler:BeamB /Herwig/Particles/p+
+set DISLuminosity:BeamEMaxB 820.*GeV
+
+# the only infrared safe cut in this case is a Q^2 cut
+# leave all other cuts as set in here
+
+set /Herwig/Cuts/NeutralCurrentCut:MinQ2 20.0*GeV2
+set /Herwig/Cuts/NeutralCurrentCut:MinW2 0.0*GeV2
+set /Herwig/Cuts/DISCuts:MHatMin 0.0*GeV
+
+cd /Herwig/MatrixElements/Matchbox/Builtin
+
+################################################################################
+# leading order
+################################################################################
+
+# insert /Herwig/Generators/DISGenerator:EventHandler:SubProcessHandlers[0] MElP2lJetLO
+
+################################################################################
+# MC@NLO-type next-to-leading order
+################################################################################
+
+# set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+# set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+# insert /Herwig/Generators/DISGenerator:EventHandler:SubProcessHandlers[0] MElP2lJetNLO
+# uncomment to run plain NLO calculation (no consistent analysis yet)
+# set MElP2lJetNLO:SubProcessGroups On
+
+################################################################################
+# POWHEG-type next-to-leading order
+################################################################################
+
+# set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+# set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+# insert /Herwig/Generators/DISGenerator:EventHandler:SubProcessHandlers[0] MElP2lJetNLOInclusive
+# set MElP2lJetNLOInclusive:BornScreening Off
+
+################################################################################
+# setup generator and analysis
+################################################################################
+
+cd /Herwig/Generators
+
+set DISGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
+set DISGenerator:EventHandler:CascadeHandler:MPIHandler NULL
+
+set DISGenerator:NumberOfEvents 100000000
+set DISGenerator:RandomNumberGenerator:Seed 31122001
+set DISGenerator:DebugLevel 1
+set DISGenerator:PrintEvent 10
+set DISGenerator:MaxErrors 1000000
+
+cd /Herwig/Generators
+
+#insert DISGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
+set /Herwig/Analysis/HepMCFile:PrintEvent 2000000
+set /Herwig/Analysis/HepMCFile:Format GenEvent
+set /Herwig/Analysis/HepMCFile:Units GeV_mm
+set /Herwig/Analysis/HepMCFile:Filename events.fifo
+
+saverun DIS DISGenerator
diff --git a/src/LEP-Matchbox.in b/src/LEP-Matchbox.in
new file mode 100644
--- /dev/null
+++ b/src/LEP-Matchbox.in
@@ -0,0 +1,76 @@
+read Matchbox.in
+
+################################################################################
+# setup the sampler
+################################################################################
+
+set /Herwig/EventHandlers/LEPHandler:Sampler /Herwig/MatrixElements/Matchbox/Samplers/Sampler
+
+################################################################################
+# setup the shower
+#
+# use -LO.in or -NLO.in for LO/NLO simulations, respectively.
+#
+################################################################################
+
+read DipoleShower.in
+# read DipoleShowerParameters-LO.in
+# read DipoleShowerParameters-NLO.in
+
+################################################################################
+# setup matrix element and matching
+#
+# uncomment the desired calculation
+#
+################################################################################
+
+cd /Herwig/MatrixElements/Matchbox/Builtin
+cp /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler MECorrectionHandler
+
+################################################################################
+# leading order
+################################################################################
+
+# insert /Herwig/Generators/LEPGenerator:EventHandler:SubProcessHandlers[0] MEllbar2JetJetLO
+
+################################################################################
+# MC@NLO-type next-to-leading order
+################################################################################
+
+# insert /Herwig/Generators/LEPGenerator:EventHandler:SubProcessHandlers[0] MEllbar2JetJetNLO
+# uncomment to run plain NLO calculation (no consistent analysis yet)
+# set MEllbar2JetJetNLO:SubProcessGroups On
+
+################################################################################
+# POWHEG-type next-to-leading order
+################################################################################
+
+# insert /Herwig/Generators/LEPGenerator:EventHandler:SubProcessHandlers[0] MEllbar2JetJetNLOInclusive
+# insert /Herwig/Generators/LEPGenerator:EventHandler:PostSubProcessHandlers[0] MECorrectionHandler
+
+################################################################################
+# setup generator and analysis
+################################################################################
+
+cd /Herwig/Generators
+
+set LEPGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
+set LEPGenerator:EventHandler:CascadeHandler:MPIHandler NULL
+
+set LEPGenerator:NumberOfEvents 100000000
+set LEPGenerator:RandomNumberGenerator:Seed 31122001
+set LEPGenerator:DebugLevel 1
+set LEPGenerator:PrintEvent 10
+set LEPGenerator:MaxErrors 10000
+
+set LEPGenerator:EventHandler:LuminosityFunction:Energy 91.2
+
+cd /Herwig/Generators
+
+#insert LEPGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
+set /Herwig/Analysis/HepMCFile:PrintEvent 1000000
+set /Herwig/Analysis/HepMCFile:Format GenEvent
+set /Herwig/Analysis/HepMCFile:Units GeV_mm
+set /Herwig/Analysis/HepMCFile:Filename events.fifo
+
+saverun LEP LEPGenerator
diff --git a/src/LHC-Matchbox.in b/src/LHC-Matchbox.in
new file mode 100644
--- /dev/null
+++ b/src/LHC-Matchbox.in
@@ -0,0 +1,88 @@
+read Matchbox.in
+
+################################################################################
+# setup the sampler
+################################################################################
+
+set /Herwig/EventHandlers/LHCHandler:Sampler /Herwig/MatrixElements/Matchbox/Samplers/Sampler
+
+################################################################################
+# setup the shower
+#
+# use -LO.in or -NLO.in for LO/NLO simulations, respectively.
+#
+################################################################################
+
+read DipoleShower.in
+# read DipoleShowerParameters-LO.in
+# read DipoleShowerParameters-NLO.in
+
+################################################################################
+# setup matrix element and matching
+#
+# uncomment the desired calculation
+#
+################################################################################
+
+cd /Herwig/EventHandlers
+
+set LHCHandler:LuminosityFunction:Energy 7000.0*GeV
+
+# the only infrared safe cut in this case is a cut on the lepton pair mass
+# leave all other cuts as set in here
+
+set /Herwig/Cuts/QCDCuts:MHatMin 0.0*GeV
+set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
+
+cd /Herwig/MatrixElements/Matchbox/Builtin
+cp /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler MECorrectionHandler
+
+################################################################################
+# leading order
+################################################################################
+
+# insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] MEPP2llbarLO
+
+################################################################################
+# MC@NLO-type next-to-leading order
+################################################################################
+
+# set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+# set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+# insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] MEPP2llbarNLO
+# uncomment to run plain NLO calculation (no consistent analysis yet)
+# set MElP2lJetNLO:SubProcessGroups On
+
+################################################################################
+# POWHEG-type next-to-leading order
+################################################################################
+
+# set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+# set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+# insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] MEPP2llbarNLOInclusive
+# insert /Herwig/Generators/LHCGenerator:EventHandler:PostSubProcessHandlers[0] MECorrectionHandler
+
+################################################################################
+# setup generator and analysis
+################################################################################
+
+cd /Herwig/Generators
+
+set LHCGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
+set LHCGenerator:EventHandler:CascadeHandler:MPIHandler NULL
+
+set LHCGenerator:NumberOfEvents 100000000
+set LHCGenerator:RandomNumberGenerator:Seed 31122001
+set LHCGenerator:DebugLevel 1
+set LHCGenerator:PrintEvent 10
+set LHCGenerator:MaxErrors 4000000
+
+cd /Herwig/Generators
+
+#insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
+set /Herwig/Analysis/HepMCFile:PrintEvent 10000000
+set /Herwig/Analysis/HepMCFile:Format GenEvent
+set /Herwig/Analysis/HepMCFile:Units GeV_mm
+set /Herwig/Analysis/HepMCFile:Filename events.fifo
+
+saverun LHC LHCGenerator
diff --git a/src/Makefile.am b/src/Makefile.am
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,167 +1,171 @@
SUBDIRS = defaults
AUTOMAKE_OPTIONS = -Wno-portability
defaultsdir = ${pkgdatadir}/defaults
bin_PROGRAMS = Herwig++
Herwig___SOURCES = Herwig++.cc herwigopts.c herwigopts.h
BUILT_SOURCES = herwigopts.c herwigopts.h
Herwig___LDFLAGS = $(AM_LDFLAGS) -export-dynamic $(THEPEGLDFLAGS)
Herwig___LDADD = $(THEPEGLIB) -ldl
Herwig___CPPFLAGS = $(AM_CPPFLAGS) \
-DHERWIG_PKGDATADIR="\"$(pkgdatadir)\"" \
-DHERWIG_PKGLIBDIR="\"$(pkglibdir)\"" \
-DTHEPEG_PKGLIBDIR="\"$(THEPEGPATH)/lib/ThePEG\""
bin_SCRIPTS = herwig-config
HELPERFILES = CMSSM40.1.1.slha RPV3.1.slha NMSSM.spc \
ADD.model \
Leptoquark.model \
MSSM.model \
MUED.model \
NMSSM.model \
RS.model \
Sextet.model \
TTBA.model \
Zprime.model
INPUTFILES = \
DIS.in \
+DIS-Matchbox.in \
GammaGamma.in \
ILC.in \
ILC-MSSM.in \
ILC-MUED.in \
ILC-RS.in \
LEP.in \
+LEP-Matchbox.in \
LHC-ADD.in \
LHC-diffractive.in \
LHC-GammaGamma.in \
LHC.in \
+LHC-Matchbox.in \
LHC-LQ.in \
LHC-MSSM.in \
LHC-MU900-2.in \
LHC-MUED.in \
LHC-NMSSM.in \
LHC-Powheg.in \
LHC-RPV.in \
LHC-RS.in \
LHC-Sextet.in \
LHC-TRP.in \
LHC-TTBA.in \
LHC-MB.in \
LHC-ZP.in \
TVT.in \
+TVT-Matchbox.in \
TVT-Powheg.in \
TVT-TTBA.in
dist_pkgdata_DATA = $(INPUTFILES) $(HELPERFILES)
pkgdata_DATA = Makefile-UserModules
CLEANFILES = HerwigDefaults.rpo \
*.run *.log *.out *.tex \
multi.test *.output probs.test chisq.value \
LHC-RS-BR.spc LHC-MSSM-BR.spc LHC-RPV-BR.spc
## checking targets ##
HerwigDefaults.rpo: Herwig++ $(srcdir)/defaults/*.in defaults/PDF.in defaults/Analysis.in $(top_builddir)/lib/*.so
./Herwig++ init -L$(top_builddir)/lib -i defaults/HerwigDefaults.in -D
check_BSM_Full=
check_BSM=
if WANT_MSSM
check_BSM += check-LHC-RPV
check_BSM_Full += check-LHC-RPV check-LHC-MSSM
endif
if WANT_NMSSM
check_BSM_Full += check-LHC-NMSSM
endif
if WANT_UED
check_BSM += check-LHC-MUED
check_BSM_Full += check-LHC-MUED
endif
if WANT_RS
check_BSM += check-LHC-RS
check_BSM_Full += check-LHC-RS
endif
if WANT_TRP
check_BSM_Full += check-LHC-TRP
endif
if WANT_ADD
check_BSM_Full += check-LHC-ADD
endif
if WANT_TTBA
check_BSM_Full += check-LHC-TTBA
endif
if WANT_ZPRIME
check_BSM_Full += check-LHC-ZP
endif
check-local: check-LHC check-LEP check-DIS check-ILC check-GammaGamma $(check_BSM) check-LHC-Powheg
check-Powheg: check-LHC-Powheg check-TVT-Powheg
check-BSM: $(check_BSM_Full)
link-helper-files:
@for i in $(HELPERFILES); do \
if test -f $(srcdir)/$$i -a ! -e $$i; then \
$(LN_S) -f $(srcdir)/$$i; fi; done
check-%: $(srcdir)/%.in HerwigDefaults.rpo link-helper-files
./Herwig++ read $< -D
./Herwig++ run $(notdir $(subst .in,.run,$<)) -N500 -d1 -D
## valgrind targets ##
VALGRIND=valgrind --leak-check=full --num-callers=25 --freelist-vol=100000000 --leak-resolution=med --trace-children=yes
valgrind: valgrind-init valgrind-read valgrind-run
valgrind-init:
$(VALGRIND) ./Herwig++ init -d1 -D -L$(top_builddir)/lib -i defaults/HerwigDefaults.in \
&> /tmp/valgrind-init.log
valgrind-read:
$(VALGRIND) ./Herwig++ read -d1 -D LHC.in &> /tmp/valgrind-read.log
valgrind-run:
$(VALGRIND) ./Herwig++ run -d1 -D -N5 LHC.run &> /tmp/valgrind-run.log
SETUPTHEPEG=$(THEPEGPATH)/bin/setupThePEG
THEPEGREPO=$(THEPEGPATH)/lib/ThePEG/ThePEGDefaults.rpo
install-data-hook:
@echo Creating repository
@./Herwig++ init -L$(DESTDIR)$(pkglibdir) -i $(DESTDIR)$(defaultsdir)/HerwigDefaults.in -r $(DESTDIR)$(pkgdatadir)/HerwigDefaults.rpo
uninstall-hook:
rm -f $(DESTDIR)$(pkgdatadir)/HerwigDefaults.rpo
register: register-with-thepeg-repo
register-with-thepeg-repo:
@if test -x "$(SETUPTHEPEG)" -a -w "$(THEPEGREPO)"; \
then echo Registering with ThePEG; \
"$(SETUPTHEPEG)" --init \
$(DESTDIR)$(defaultsdir)/HerwigDefaults.in \
-r "$(THEPEGREPO)" -o "$(THEPEGREPO)" \
-l$(DESTDIR)$(pkglibdir) ; \
fi
unregister : unregister-from-thepeg-repo
unregister-from-thepeg-repo:
@if test -x "$(SETUPTHEPEG)" -a -w "$(THEPEGREPO)"; \
then echo Unregistering with ThePEG; \
"$(SETUPTHEPEG)" --init defaults/HerwigCleanup.in \
-r "$(THEPEGREPO)" -o "$(THEPEGREPO)" \
-l$(DESTDIR)$(pkglibdir) ; \
fi
EXTRA_DIST = herwigopts.ggo
GENGETOPT = gengetopt
%opts.h %opts.c : %opts.ggo
$(GENGETOPT) < $<
diff --git a/src/TVT-Matchbox.in b/src/TVT-Matchbox.in
new file mode 100644
--- /dev/null
+++ b/src/TVT-Matchbox.in
@@ -0,0 +1,89 @@
+read Matchbox.in
+
+################################################################################
+# setup the sampler
+################################################################################
+
+set /Herwig/EventHandlers/LHCHandler:Sampler /Herwig/MatrixElements/Matchbox/Samplers/Sampler
+
+################################################################################
+# setup the shower
+#
+# use -LO.in or -NLO.in for LO/NLO simulations, respectively.
+#
+################################################################################
+
+read DipoleShower.in
+# read DipoleShowerParameters-LO.in
+# read DipoleShowerParameters-NLO.in
+
+################################################################################
+# setup matrix element and matching
+#
+# uncomment the desired calculation
+#
+################################################################################
+
+cd /Herwig/EventHandlers
+
+set LHCHandler:LuminosityFunction:Energy 1800.0*GeV
+set LHCHandler:BeamB /Herwig/Particles/pbar-
+
+# the only infrared safe cut in this case is a cut on the lepton mass
+# leave all other cuts as set in here
+
+set /Herwig/Cuts/QCDCuts:MHatMin 0.0*GeV
+set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
+
+cd /Herwig/MatrixElements/Matchbox/Builtin
+cp /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler MECorrectionHandler
+
+################################################################################
+# leading order
+################################################################################
+
+# insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] MEPP2llbarLO
+
+################################################################################
+# MC@NLO-type next-to-leading order
+################################################################################
+
+# set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+# set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+# insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] MEPP2llbarNLO
+# uncomment to run plain NLO calculation (no consistent analysis yet)
+# set MElP2lJetNLO:SubProcessGroups On
+
+################################################################################
+# POWHEG-type next-to-leading order
+################################################################################
+
+# set /Herwig/Particles/p+:PDF /Herwig/Partons/MRST-NLO
+# set /Herwig/Particles/pbar-:PDF /Herwig/Partons/MRST-NLO
+# insert /Herwig/Generators/LHCGenerator:EventHandler:SubProcessHandlers[0] MEPP2llbarNLOInclusive
+# insert /Herwig/Generators/LHCGenerator:EventHandler:PostSubProcessHandlers[0] MECorrectionHandler
+
+################################################################################
+# setup generator and analysis
+################################################################################
+
+cd /Herwig/Generators
+
+set LHCGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
+set LHCGenerator:EventHandler:CascadeHandler:MPIHandler NULL
+
+set LHCGenerator:NumberOfEvents 100000000
+set LHCGenerator:RandomNumberGenerator:Seed 31122001
+set LHCGenerator:DebugLevel 1
+set LHCGenerator:PrintEvent 10
+set LHCGenerator:MaxErrors 4000000
+
+cd /Herwig/Generators
+
+#insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
+set /Herwig/Analysis/HepMCFile:PrintEvent 10000000
+set /Herwig/Analysis/HepMCFile:Format GenEvent
+set /Herwig/Analysis/HepMCFile:Units GeV_mm
+set /Herwig/Analysis/HepMCFile:Filename events.fifo
+
+saverun TVT LHCGenerator
diff --git a/src/defaults/DipoleShower.in b/src/defaults/DipoleShower.in
--- a/src/defaults/DipoleShower.in
+++ b/src/defaults/DipoleShower.in
@@ -1,300 +1,316 @@
################################################################################
# Setup the dipole shower
################################################################################
library HwDipoleShower.so
mkdir /Herwig/DipoleShower
cd /Herwig/DipoleShower
create Herwig::DipoleShowerHandler DipoleShowerHandler
################################################################################
#
# /!\ Nothing to be done below here, except you exactly know
# what you're doing.
#
# Really.
#
################################################################################
# .....
# ....
# ...
# ..
# .
#
################################################################################
# zero quark masses for the time being
################################################################################
set /Herwig/Particles/d:NominalMass 0.*GeV
set /Herwig/Particles/dbar:NominalMass 0.*GeV
set /Herwig/Particles/u:NominalMass 0.*GeV
set /Herwig/Particles/ubar:NominalMass 0.*GeV
set /Herwig/Particles/c:NominalMass 0.*GeV
set /Herwig/Particles/cbar:NominalMass 0.*GeV
set /Herwig/Particles/s:NominalMass 0.*GeV
set /Herwig/Particles/sbar:NominalMass 0.*GeV
set /Herwig/Particles/b:NominalMass 0.*GeV
set /Herwig/Particles/bbar:NominalMass 0.*GeV
################################################################################
# Setup the underlying event and fix missing reference.
################################################################################
set DipoleShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
set DipoleShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
set DipoleShowerHandler:Evolver /Herwig/Shower/Evolver
################################################################################
# Setup the ordering.
################################################################################
create Herwig::DipoleChainOrdering ChainPtOrdering
################################################################################
# Setup the reshuffler.
################################################################################
create Herwig::ConstituentReshuffler ConstituentReshuffler
set DipoleShowerHandler:ConstituentReshuffler ConstituentReshuffler
################################################################################
# Setup intrinsic pt.
################################################################################
create Herwig::IntrinsicPtGenerator IntrinsicPtGenerator
set DipoleShowerHandler:IntrinsicPtGenerator IntrinsicPtGenerator
################################################################################
# Setup the alphas
################################################################################
library HwDipoleShowerAlphaS.so
create matchbox::lo_alpha_s LOAlphaS
set LOAlphaS:min_active_flavours 5
create matchbox::nlo_alpha_s NLOAlphaS
set NLOAlphaS:exact_evaluation large_scale
set NLOAlphaS:min_active_flavours 5
################################################################################
# Setup the splitting kinematics.
################################################################################
mkdir /Herwig/DipoleShower/Kinematics
cd /Herwig/DipoleShower/Kinematics
create Herwig::FFLightKinematics FFLightKinematics
create Herwig::FILightKinematics FILightKinematics
create Herwig::IFLightKinematics IFLightKinematics
create Herwig::IILightKinematics IILightKinematics
################################################################################
# Setup the kernels.
################################################################################
mkdir /Herwig/DipoleShower/Kernels
cd /Herwig/DipoleShower/Kernels
################################################################################
# FF
################################################################################
create Herwig::FFgx2ggxDipoleKernel FFgx2ggxDipoleKernel
set FFgx2ggxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/FFLightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ggxDipoleKernel
create Herwig::FFqx2qgxDipoleKernel FFqx2qgxDipoleKernel
set FFqx2qgxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/FFLightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFqx2qgxDipoleKernel
create Herwig::FFgx2qqxDipoleKernel FFgx2qqxDipoleKernel
set FFgx2qqxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/FFLightKinematics
cp FFgx2qqxDipoleKernel FFgx2ddxDipoleKernel
set FFgx2ddxDipoleKernel:Flavour /Herwig/Particles/d
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ddxDipoleKernel
cp FFgx2qqxDipoleKernel FFgx2uuxDipoleKernel
set FFgx2uuxDipoleKernel:Flavour /Herwig/Particles/u
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2uuxDipoleKernel
cp FFgx2qqxDipoleKernel FFgx2ccxDipoleKernel
set FFgx2ccxDipoleKernel:Flavour /Herwig/Particles/c
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ccxDipoleKernel
cp FFgx2qqxDipoleKernel FFgx2ssxDipoleKernel
set FFgx2ssxDipoleKernel:Flavour /Herwig/Particles/s
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ssxDipoleKernel
cp FFgx2qqxDipoleKernel FFgx2bbxDipoleKernel
set FFgx2bbxDipoleKernel:Flavour /Herwig/Particles/b
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2bbxDipoleKernel
################################################################################
# create the pdf ratio object
################################################################################
create Herwig::PDFRatio PDFRatio
################################################################################
# FI
################################################################################
create Herwig::FIgx2ggxDipoleKernel FIgx2ggxDipoleKernel
set FIgx2ggxDipoleKernel:PDFRatio PDFRatio
set FIgx2ggxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/FILightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIgx2ggxDipoleKernel
create Herwig::FIqx2qgxDipoleKernel FIqx2qgxDipoleKernel
set FIqx2qgxDipoleKernel:PDFRatio PDFRatio
set FIqx2qgxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/FILightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIqx2qgxDipoleKernel
create Herwig::FIgx2qqxDipoleKernel FIgx2qqxDipoleKernel
set FIgx2qqxDipoleKernel:PDFRatio PDFRatio
set FIgx2qqxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/FILightKinematics
cp FIgx2qqxDipoleKernel FIgx2ddxDipoleKernel
set FIgx2ddxDipoleKernel:Flavour /Herwig/Particles/d
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIgx2ddxDipoleKernel
cp FIgx2qqxDipoleKernel FIgx2uuxDipoleKernel
set FIgx2uuxDipoleKernel:Flavour /Herwig/Particles/u
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIgx2uuxDipoleKernel
cp FIgx2qqxDipoleKernel FIgx2ccxDipoleKernel
set FIgx2ccxDipoleKernel:Flavour /Herwig/Particles/c
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIgx2ccxDipoleKernel
cp FIgx2qqxDipoleKernel FIgx2ssxDipoleKernel
set FIgx2ssxDipoleKernel:Flavour /Herwig/Particles/s
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIgx2ssxDipoleKernel
cp FIgx2qqxDipoleKernel FIgx2bbxDipoleKernel
set FIgx2bbxDipoleKernel:Flavour /Herwig/Particles/b
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FIgx2bbxDipoleKernel
################################################################################
# IF
################################################################################
create Herwig::IFgx2ggxDipoleKernel IFgx2ggxDipoleKernel
set IFgx2ggxDipoleKernel:PDFRatio PDFRatio
set IFgx2ggxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IFLightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2ggxDipoleKernel
create Herwig::IFqx2qgxDipoleKernel IFqx2qgxDipoleKernel
set IFqx2qgxDipoleKernel:PDFRatio PDFRatio
set IFqx2qgxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IFLightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFqx2qgxDipoleKernel
create Herwig::IFqx2gqxDipoleKernel IFqx2gqxDipoleKernel
set IFqx2gqxDipoleKernel:PDFRatio PDFRatio
set IFqx2gqxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IFLightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFqx2gqxDipoleKernel
create Herwig::IFgx2qqxDipoleKernel IFgx2qqxDipoleKernel
set IFgx2qqxDipoleKernel:PDFRatio PDFRatio
set IFgx2qqxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IFLightKinematics
cp IFgx2qqxDipoleKernel IFgx2ddbarxDipoleKernel
set IFgx2ddbarxDipoleKernel:Flavour /Herwig/Particles/d
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2ddbarxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2dbardxDipoleKernel
set IFgx2dbardxDipoleKernel:Flavour /Herwig/Particles/dbar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2dbardxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2uubarxDipoleKernel
set IFgx2uubarxDipoleKernel:Flavour /Herwig/Particles/u
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2uubarxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2ubaruxDipoleKernel
set IFgx2ubaruxDipoleKernel:Flavour /Herwig/Particles/ubar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2ubaruxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2ccbarxDipoleKernel
set IFgx2ccbarxDipoleKernel:Flavour /Herwig/Particles/c
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2ccbarxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2cbarcxDipoleKernel
set IFgx2cbarcxDipoleKernel:Flavour /Herwig/Particles/cbar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2cbarcxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2ssbarxDipoleKernel
set IFgx2ssbarxDipoleKernel:Flavour /Herwig/Particles/s
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2ssbarxDipoleKernel
cp IFgx2qqxDipoleKernel IFgx2sbarsxDipoleKernel
set IFgx2sbarsxDipoleKernel:Flavour /Herwig/Particles/sbar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2sbarsxDipoleKernel
+cp IFgx2qqxDipoleKernel IFgx2bbbarxDipoleKernel
+set IFgx2bbbarxDipoleKernel:Flavour /Herwig/Particles/b
+insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2bbbarxDipoleKernel
+
+cp IFgx2qqxDipoleKernel IFgx2bbarbxDipoleKernel
+set IFgx2bbarbxDipoleKernel:Flavour /Herwig/Particles/bbar
+insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IFgx2bbarbxDipoleKernel
+
################################################################################
# II
################################################################################
create Herwig::IIgx2ggxDipoleKernel IIgx2ggxDipoleKernel
set IIgx2ggxDipoleKernel:PDFRatio PDFRatio
set IIgx2ggxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IILightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2ggxDipoleKernel
create Herwig::IIqx2qgxDipoleKernel IIqx2qgxDipoleKernel
set IIqx2qgxDipoleKernel:PDFRatio PDFRatio
set IIqx2qgxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IILightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIqx2qgxDipoleKernel
create Herwig::IIqx2gqxDipoleKernel IIqx2gqxDipoleKernel
set IIqx2gqxDipoleKernel:PDFRatio PDFRatio
set IIqx2gqxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IILightKinematics
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIqx2gqxDipoleKernel
create Herwig::IIgx2qqxDipoleKernel IIgx2qqxDipoleKernel
set IIgx2qqxDipoleKernel:PDFRatio PDFRatio
set IIgx2qqxDipoleKernel:SplittingKinematics /Herwig/DipoleShower/Kinematics/IILightKinematics
cp IIgx2qqxDipoleKernel IIgx2ddbarxDipoleKernel
set IIgx2ddbarxDipoleKernel:Flavour /Herwig/Particles/d
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2ddbarxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2dbardxDipoleKernel
set IIgx2dbardxDipoleKernel:Flavour /Herwig/Particles/dbar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2dbardxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2uubarxDipoleKernel
set IIgx2uubarxDipoleKernel:Flavour /Herwig/Particles/u
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2uubarxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2ubaruxDipoleKernel
set IIgx2ubaruxDipoleKernel:Flavour /Herwig/Particles/ubar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2ubaruxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2ccbarxDipoleKernel
set IIgx2ccbarxDipoleKernel:Flavour /Herwig/Particles/c
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2ccbarxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2cbarcxDipoleKernel
set IIgx2cbarcxDipoleKernel:Flavour /Herwig/Particles/cbar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2cbarcxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2ssbarxDipoleKernel
set IIgx2ssbarxDipoleKernel:Flavour /Herwig/Particles/s
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2ssbarxDipoleKernel
cp IIgx2qqxDipoleKernel IIgx2sbarsxDipoleKernel
set IIgx2sbarsxDipoleKernel:Flavour /Herwig/Particles/sbar
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2sbarsxDipoleKernel
+cp IIgx2qqxDipoleKernel IIgx2bbbarxDipoleKernel
+set IIgx2bbbarxDipoleKernel:Flavour /Herwig/Particles/b
+insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2bbbarxDipoleKernel
+
+cp IIgx2qqxDipoleKernel IIgx2bbarbxDipoleKernel
+set IIgx2bbarbxDipoleKernel:Flavour /Herwig/Particles/bbar
+insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 IIgx2bbarbxDipoleKernel
+
cd /
diff --git a/src/defaults/DipoleShowerParameters-LO.in b/src/defaults/DipoleShowerParameters-LO.in
--- a/src/defaults/DipoleShowerParameters-LO.in
+++ b/src/defaults/DipoleShowerParameters-LO.in
@@ -1,125 +1,131 @@
################################################################################
# Setup the dipole shower parameters
################################################################################
cd /Herwig/DipoleShower
################################################################################
# setup alphas
################################################################################
cp LOAlphaS AlphaS
set AlphaS:input_alpha_s 0.1131785
set AlphaS:freezing_scale 0.0*GeV
+set /Herwig/Model:QCD/RunningAlphaS AlphaS
+
################################################################################
# shower parameters
################################################################################
set DipoleShowerHandler:GlobalAlphaS AlphaS
set DipoleShowerHandler:EvolutionOrdering ChainPtOrdering
set IntrinsicPtGenerator:ValenceIntrinsicPtScale 1.68463*GeV
set IntrinsicPtGenerator:SeaIntrinsicPtScale 1.29001*GeV
cd /Herwig/DipoleShower/Kinematics
set FFLightKinematics:IRCutoff 1.416023*GeV
set FILightKinematics:IRCutoff 0.796205*GeV
set IFLightKinematics:IRCutoff 0.796205*GeV
set IILightKinematics:IRCutoff 0.367359*GeV
################################################################################
# shower parameters at boundary to non-perturbative domain
################################################################################
set /Herwig/Particles/g:ConstituentMass 1.080386*GeV
cd /Herwig/DipoleShower/Kernels
set FFgx2ggxDipoleKernel:ScreeningScale 0.2427254*GeV
set FFqx2qgxDipoleKernel:ScreeningScale 0.2427254*GeV
set FFgx2uuxDipoleKernel:ScreeningScale 0.2427254*GeV
set FFgx2ddxDipoleKernel:ScreeningScale 0.2427254*GeV
set FFgx2ssxDipoleKernel:ScreeningScale 0.2427254*GeV
set FFgx2ccxDipoleKernel:ScreeningScale 0.2427254*GeV
set FFgx2bbxDipoleKernel:ScreeningScale 0.2427254*GeV
set FIgx2ggxDipoleKernel:ScreeningScale 1.355894*GeV
set FIqx2qgxDipoleKernel:ScreeningScale 1.355894*GeV
set FIgx2uuxDipoleKernel:ScreeningScale 1.355894*GeV
set FIgx2ddxDipoleKernel:ScreeningScale 1.355894*GeV
set FIgx2ssxDipoleKernel:ScreeningScale 1.355894*GeV
set FIgx2ccxDipoleKernel:ScreeningScale 1.355894*GeV
set FIgx2bbxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2ggxDipoleKernel:ScreeningScale 1.355894*GeV
set IFqx2qgxDipoleKernel:ScreeningScale 1.355894*GeV
set IFqx2gqxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2uubarxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2ubaruxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2ddbarxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2dbardxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2ssbarxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2sbarsxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2ccbarxDipoleKernel:ScreeningScale 1.355894*GeV
set IFgx2cbarcxDipoleKernel:ScreeningScale 1.355894*GeV
+set IFgx2bbbarxDipoleKernel:ScreeningScale 1.355894*GeV
+set IFgx2bbarbxDipoleKernel:ScreeningScale 1.355894*GeV
set IIgx2ggxDipoleKernel:ScreeningScale 0.205854*GeV
set IIqx2qgxDipoleKernel:ScreeningScale 0.205854*GeV
set IIqx2gqxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2uubarxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2ubaruxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2ddbarxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2dbardxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2ssbarxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2sbarsxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2ccbarxDipoleKernel:ScreeningScale 0.205854*GeV
set IIgx2cbarcxDipoleKernel:ScreeningScale 0.205854*GeV
+set IIgx2bbbarxDipoleKernel:ScreeningScale 0.205854*GeV
+set IIgx2bbarbxDipoleKernel:ScreeningScale 0.205854*GeV
################################################################################
# hadronization parameters
################################################################################
cd /Herwig/Hadronization
set ClusterFissioner:ClMaxLight 4.17032*GeV
set ClusterFissioner:ClPowLight 5.734681
set ClusterFissioner:PSplitLight 0.7651726
set ClusterDecayer:ClDirLight 1
set ClusterDecayer:ClSmrLight 4.548755
set ClusterFissioner:ClMaxCharm 4.17032*GeV
set ClusterFissioner:ClPowCharm 5.734681
set ClusterFissioner:PSplitCharm 0.7651726
set ClusterDecayer:ClDirCharm 1
set ClusterDecayer:ClSmrCharm 4.548755
set LightClusterDecayer:SingleHadronLimitCharm 0.0
set ClusterFissioner:ClMaxBottom 4.17032*GeV
set ClusterFissioner:ClPowBottom 5.734681
set ClusterFissioner:PSplitBottom 0.7651726
set ClusterDecayer:ClDirBottom 1
set ClusterDecayer:ClSmrBottom 4.548755
set LightClusterDecayer:SingleHadronLimitBottom 0.0
set HadronSelector:PwtUquark 1.0
set HadronSelector:PwtDquark 1.0
set HadronSelector:PwtSquark 1.0
set HadronSelector:PwtCquark 1.0
set HadronSelector:PwtBquark 1.0
set HadronSelector:PwtDIquark 1.0
set HadronSelector:SngWt 1.0
set HadronSelector:DecWt 1.0
cd /
diff --git a/src/defaults/DipoleShowerParameters-NLO.in b/src/defaults/DipoleShowerParameters-NLO.in
--- a/src/defaults/DipoleShowerParameters-NLO.in
+++ b/src/defaults/DipoleShowerParameters-NLO.in
@@ -1,140 +1,144 @@
################################################################################
# Setup the dipole shower parameters
################################################################################
cd /Herwig/DipoleShower
################################################################################
# setup alphas
################################################################################
cp NLOAlphaS AlphaS
set AlphaS:input_alpha_s 0.1175576
set AlphaS:freezing_scale 0.0*GeV
+set /Herwig/Model:QCD/RunningAlphaS AlphaS
+
################################################################################
# matching parameters
################################################################################
-set /Herwig/Model:QCD/RunningAlphaS AlphaS
+set /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler:FFPtCut 1.245924*GeV
+set /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler:FFScreeningScale 0.0*GeV
-set MECorrectionHandler:FFPtCut 1.245924*GeV
-set MECorrectionHandler:FFScreeningScale 0.0*GeV
+set /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler:FIPtCut 0.706986*GeV
+set /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler:FIScreeningScale 1.017780*GeV
-set MECorrectionHandler:FIPtCut 0.706986*GeV
-set MECorrectionHandler:FIScreeningScale 1.017780*GeV
-
-set MECorrectionHandler:IIPtCut 0.275894*GeV
-set MECorrectionHandler:IIScreeningScale 0.254028*GeV
+set /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler:IIPtCut 0.275894*GeV
+set /Herwig/MatrixElements/Matchbox/Utility/MECorrectionHandler:IIScreeningScale 0.254028*GeV
################################################################################
# shower parameters
################################################################################
set DipoleShowerHandler:GlobalAlphaS AlphaS
set DipoleShowerHandler:EvolutionOrdering ChainPtOrdering
set IntrinsicPtGenerator:ValenceIntrinsicPtScale 1.26905*GeV
set IntrinsicPtGenerator:SeaIntrinsicPtScale 1.1613*GeV
cd /Herwig/DipoleShower/Kinematics
set FFLightKinematics:IRCutoff 1.245924*GeV
set FILightKinematics:IRCutoff 0.706986*GeV
set IFLightKinematics:IRCutoff 0.706986*GeV
set IILightKinematics:IRCutoff 0.275894*GeV
################################################################################
# shower parameters at boundary to non-perturbative domain
################################################################################
set /Herwig/Particles/g:ConstituentMass 1.006752*GeV
cd /Herwig/DipoleShower/Kernels
set FFgx2ggxDipoleKernel:ScreeningScale 0.0*GeV
set FFqx2qgxDipoleKernel:ScreeningScale 0.0*GeV
set FFgx2uuxDipoleKernel:ScreeningScale 0.0*GeV
set FFgx2ddxDipoleKernel:ScreeningScale 0.0*GeV
set FFgx2ssxDipoleKernel:ScreeningScale 0.0*GeV
set FFgx2ccxDipoleKernel:ScreeningScale 0.0*GeV
set FFgx2bbxDipoleKernel:ScreeningScale 0.0*GeV
set FIgx2ggxDipoleKernel:ScreeningScale 1.017780*GeV
set FIqx2qgxDipoleKernel:ScreeningScale 1.017780*GeV
set FIgx2uuxDipoleKernel:ScreeningScale 1.017780*GeV
set FIgx2ddxDipoleKernel:ScreeningScale 1.017780*GeV
set FIgx2ssxDipoleKernel:ScreeningScale 1.017780*GeV
set FIgx2ccxDipoleKernel:ScreeningScale 1.017780*GeV
set FIgx2bbxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2ggxDipoleKernel:ScreeningScale 1.017780*GeV
set IFqx2qgxDipoleKernel:ScreeningScale 1.017780*GeV
set IFqx2gqxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2uubarxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2ubaruxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2ddbarxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2dbardxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2ssbarxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2sbarsxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2ccbarxDipoleKernel:ScreeningScale 1.017780*GeV
set IFgx2cbarcxDipoleKernel:ScreeningScale 1.017780*GeV
+set IFgx2bbbarxDipoleKernel:ScreeningScale 1.017780*GeV
+set IFgx2bbarbxDipoleKernel:ScreeningScale 1.017780*GeV
set IIgx2ggxDipoleKernel:ScreeningScale 0.254028*GeV
set IIqx2qgxDipoleKernel:ScreeningScale 0.254028*GeV
set IIqx2gqxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2uubarxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2ubaruxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2ddbarxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2dbardxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2ssbarxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2sbarsxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2ccbarxDipoleKernel:ScreeningScale 0.254028*GeV
set IIgx2cbarcxDipoleKernel:ScreeningScale 0.254028*GeV
+set IIgx2bbbarxDipoleKernel:ScreeningScale 0.254028*GeV
+set IIgx2bbarbxDipoleKernel:ScreeningScale 0.254028*GeV
################################################################################
# hadronization parameters
################################################################################
cd /Herwig/Hadronization
set ClusterFissioner:ClMaxLight 3.66643*GeV
set ClusterFissioner:ClPowLight 5.682173
set ClusterFissioner:PSplitLight 0.7712919
set ClusterDecayer:ClDirLight 1
set ClusterDecayer:ClSmrLight 3.118342
set ClusterFissioner:ClMaxCharm 3.66643*GeV
set ClusterFissioner:ClPowCharm 5.682173
set ClusterFissioner:PSplitCharm 0.7712919
set ClusterDecayer:ClDirCharm 1
set ClusterDecayer:ClSmrCharm 3.118342
set LightClusterDecayer:SingleHadronLimitCharm 0.0
set ClusterFissioner:ClMaxBottom 3.66643*GeV
set ClusterFissioner:ClPowBottom 5.682173
set ClusterFissioner:PSplitBottom 0.7712919
set ClusterDecayer:ClDirBottom 1
set ClusterDecayer:ClSmrBottom 3.118342
set LightClusterDecayer:SingleHadronLimitBottom 0.0
set HadronSelector:PwtUquark 1.0
set HadronSelector:PwtDquark 1.0
set HadronSelector:PwtSquark 1.0
set HadronSelector:PwtCquark 1.0
set HadronSelector:PwtBquark 1.0
set HadronSelector:PwtDIquark 1.0
set HadronSelector:SngWt 1.0
set HadronSelector:DecWt 1.0
cd /
diff --git a/src/defaults/Matchbox.in b/src/defaults/Matchbox.in
--- a/src/defaults/Matchbox.in
+++ b/src/defaults/Matchbox.in
@@ -1,324 +1,446 @@
################################################################################
# Input file to setup Matchbox NLO matrix elements
-#
-# The Matchbox NLO matrix elements are provided as SubProcessHandlers
-# containing a complete and consistent NLO setup.
-#
-# Important notes on usage and matching issues:
-# -------------------------------------------------------------------------
-#
-# * Matrix elements with names ending in ...NLO can be directly run
-# with the dipole shower in case of which a MC@NLO-type matching
-# is carried out
-#
-# * Matrix elements with names ending in ...NLOInclusive serve
-# as fixed order calculations for POWHEG Matching, where
-# inserting /Herwig/DipoleShower/MECorrectionHandler as
-# PostSubprocessHandler is mandatory to generate the real
-# emission piece.
-#
-# * Be aware that any other usage will produce inconsistencies.
-#
-# * In addition, Matrix elements with names ending in ...NLO can be run
-# as plain NLO using dedicated analysis handlers without any showering
-# or hadronization. This featrue is currently not completely supported.
-#
-# Matrix elements currently provided are:
-# -------------------------------------------------------------------------
-#
-# * e+e- to two jets
-#
-# /Herwig/MatrixElements/Matchbox/LeptonLepton/MEllbar2JetJetLO
-# /Herwig/MatrixElements/Matchbox/LeptonLepton/MEllbar2JetJetNLO
-# /Herwig/MatrixElements/Matchbox/LeptonLepton/MEllbar2JetJetNLOInclusive
-#
-# * DIS (charged lepton, neutral current)
-#
-# /Herwig/MatrixElements/Matchbox/LeptonHadron/MElP2lJetLO
-# /Herwig/MatrixElements/Matchbox/LeptonHadron/MElP2lJetNLO
-# /Herwig/MatrixElements/Matchbox/LeptonHadron/MElP2lJetNLOInclusive
-#
-# * Drell-Yan (Z/gamma)
-#
-# /Herwig/MatrixElements/Matchbox/HadronHadron/MEPP2llbarLO
-# /Herwig/MatrixElements/Matchbox/HadronHadron/MEPP2llbarNLO
-# /Herwig/MatrixElements/Matchbox/HadronHadron/MEPP2llbarNLOInclusive
-#
-#
################################################################################
################################################################################
#
# /!\ Nothing to be done below here, except you exactly know
# what you're doing.
#
# Really.
#
################################################################################
# .....
# ....
# ...
# ..
# .
#
-
+library HwMatchbox.so
+mkdir /Herwig/MatrixElements/Matchbox
+cd /Herwig/MatrixElements/Matchbox
################################################################################
# For the moment, we can only deal with massless partons
################################################################################
set /Herwig/Particles/d:NominalMass 0*GeV
set /Herwig/Particles/dbar:NominalMass 0*GeV
set /Herwig/Particles/u:NominalMass 0*GeV
set /Herwig/Particles/ubar:NominalMass 0*GeV
set /Herwig/Particles/s:NominalMass 0*GeV
set /Herwig/Particles/sbar:NominalMass 0*GeV
set /Herwig/Particles/c:NominalMass 0*GeV
set /Herwig/Particles/cbar:NominalMass 0*GeV
set /Herwig/Particles/b:NominalMass 0*GeV
set /Herwig/Particles/bbar:NominalMass 0*GeV
-library HwMatchbox.so
-
-mkdir /Herwig/MatrixElements/Matchbox
-cd /Herwig/MatrixElements/Matchbox
-
################################################################################
# Phasespace generators
################################################################################
+
+cd /Herwig/MatrixElements/Matchbox
mkdir Phasespace
cd Phasespace
create Herwig::MatchboxRambo Rambo
create Herwig::TreePhasespace TreePhasespace
+################################################################################
+# Scale choices
+################################################################################
+
cd /Herwig/MatrixElements/Matchbox
-mkdir LeptonLepton
-cd LeptonLepton
+mkdir Scales
+cd Scales
+
+create Herwig::MatchboxScaleChoice SHatScale
+cp SHatScale FixedScale
+set FixedScale:FixedScale 100.*GeV
+create Herwig::MatchboxPtScale MaxPtScale
+create Herwig::MatchboxLeptonMassScale LeptonMassScale
+
+################################################################################
+# AlphaS
+################################################################################
+
+cd /Herwig/MatrixElements/Matchbox
+mkdir AlphaS
+cd AlphaS
+
+library HwDipoleShowerAlphaS.so
+create matchbox::lo_alpha_s LOAlphaS
+set LOAlphaS:min_active_flavours 5
+
+create matchbox::nlo_alpha_s NLOAlphaS
+set NLOAlphaS:exact_evaluation large_scale
+set NLOAlphaS:min_active_flavours 5
+
+################################################################################
+# Utilities
+################################################################################
+
+cd /Herwig/MatrixElements/Matchbox
+mkdir Utility
+cd Utility
+
+create Herwig::PowhegSplittingGenerator MECorrectionHandler
+
+set MECorrectionHandler:FFPtCut 1.0*GeV
+set MECorrectionHandler:FFScreeningScale 0.0*GeV
+
+set MECorrectionHandler:FIPtCut 1.0*GeV
+set MECorrectionHandler:FIScreeningScale 0.0*GeV
+
+set MECorrectionHandler:IIPtCut 1.0*GeV
+set MECorrectionHandler:IIScreeningScale 0.0*GeV
+
+create Herwig::SimpleColourBasis SimpleColourBasis
+
+create Herwig::Tree2toNGenerator DiagramGenerator
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFGVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFHVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFPVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFWVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/FFZVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/GGGVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HGGVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HHHVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/HPPVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/WWHVertex
+insert DiagramGenerator:Vertices 0 /Herwig/Vertices/WWWVertex
+
+create Herwig::MatchboxMECache Cache
+
+create Herwig::MatchboxFactory Factory
+
+set Factory:DiagramGenerator /Herwig/MatrixElements/Matchbox/Utility/DiagramGenerator
+set Factory:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
+set Factory:ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/SHatScale
+set Factory:Cache /Herwig/MatrixElements/Matchbox/Utility/Cache
+
+do Factory:StartParticleGroup p
+insert Factory:ParticleGroup 0 /Herwig/Particles/u
+insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
+insert Factory:ParticleGroup 0 /Herwig/Particles/d
+insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/s
+insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/c
+insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/b
+insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/g
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup j
+insert Factory:ParticleGroup 0 /Herwig/Particles/u
+insert Factory:ParticleGroup 0 /Herwig/Particles/ubar
+insert Factory:ParticleGroup 0 /Herwig/Particles/d
+insert Factory:ParticleGroup 0 /Herwig/Particles/dbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/s
+insert Factory:ParticleGroup 0 /Herwig/Particles/sbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/c
+insert Factory:ParticleGroup 0 /Herwig/Particles/cbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/b
+insert Factory:ParticleGroup 0 /Herwig/Particles/bbar
+insert Factory:ParticleGroup 0 /Herwig/Particles/g
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup e+
+insert Factory:ParticleGroup 0 /Herwig/Particles/e+
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup e-
+insert Factory:ParticleGroup 0 /Herwig/Particles/e-
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup mu+
+insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup mu-
+insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup l+
+insert Factory:ParticleGroup 0 /Herwig/Particles/e+
+insert Factory:ParticleGroup 0 /Herwig/Particles/mu+
+do Factory:EndParticleGroup
+
+do Factory:StartParticleGroup l-
+insert Factory:ParticleGroup 0 /Herwig/Particles/e-
+insert Factory:ParticleGroup 0 /Herwig/Particles/mu-
+do Factory:EndParticleGroup
+
+################################################################################
+# Integrators and samplers
+################################################################################
+
+cd /Herwig/MatrixElements/Matchbox
+mkdir Samplers
+cd Samplers
+
+library HwExsample2.so
+
+create Herwig::ExSampler ExSampler
+create Herwig::BinSampler FlatSampler
+create Herwig::ProjectingSampler ProjectingSampler
+
+create Herwig::GeneralSampler Sampler
+set Sampler:BinSampler ExSampler
+
+################################################################################
+# Builtin simple matrix elements
+################################################################################
+
+cd /Herwig/MatrixElements/Matchbox
+mkdir Builtin
+cd Builtin
################################################################################
# e+e- to jj
################################################################################
create Herwig::MatchboxMEllbar2qqbar MEllbar2qqbar
set MEllbar2qqbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MEllbar2qqbar:LeptonFlavours[0] /Herwig/Particles/e-
insert MEllbar2qqbar:QuarkFlavours[0] /Herwig/Particles/d
insert MEllbar2qqbar:QuarkFlavours[1] /Herwig/Particles/u
insert MEllbar2qqbar:QuarkFlavours[2] /Herwig/Particles/s
insert MEllbar2qqbar:QuarkFlavours[3] /Herwig/Particles/c
insert MEllbar2qqbar:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMEllbar2qqbarg MEllbar2qqbarg
+set MEllbar2qqbarg:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MEllbar2qqbarg:LeptonFlavours[0] /Herwig/Particles/e-
insert MEllbar2qqbarg:QuarkFlavours[0] /Herwig/Particles/d
insert MEllbar2qqbarg:QuarkFlavours[1] /Herwig/Particles/u
insert MEllbar2qqbarg:QuarkFlavours[2] /Herwig/Particles/s
insert MEllbar2qqbarg:QuarkFlavours[3] /Herwig/Particles/c
insert MEllbar2qqbarg:QuarkFlavours[4] /Herwig/Particles/b
#-------------------------------------------------------------------------------
create ThePEG::SubProcessHandler MEllbar2JetJetLO
set MEllbar2JetJetLO:PartonExtractor /Herwig/Partons/EEExtractor
set MEllbar2JetJetLO:Cuts /Herwig/Cuts/EECuts
insert MEllbar2JetJetLO:MatrixElements[0] MEllbar2qqbar
+create ThePEG::SubProcessHandler MEllbar2JetJetJetLO
+set MEllbar2JetJetJetLO:PartonExtractor /Herwig/Partons/EEExtractor
+set MEllbar2JetJetJetLO:Cuts /Herwig/Cuts/EECuts
+insert MEllbar2JetJetJetLO:MatrixElements[0] MEllbar2qqbarg
+
create Herwig::MatchboxFactory MEllbar2JetJetNLO
set MEllbar2JetJetNLO:PartonExtractor /Herwig/Partons/EEExtractor
set MEllbar2JetJetNLO:Cuts /Herwig/Cuts/EECuts
insert MEllbar2JetJetNLO:BornMEs[0] MEllbar2qqbar
insert MEllbar2JetJetNLO:RealEmissionMEs[0] MEllbar2qqbarg
create Herwig::PowhegFactory MEllbar2JetJetNLOInclusive
set MEllbar2JetJetNLOInclusive:PartonExtractor /Herwig/Partons/EEExtractor
set MEllbar2JetJetNLOInclusive:Cuts /Herwig/Cuts/EECuts
set MEllbar2JetJetNLOInclusive:MatchboxFactory MEllbar2JetJetNLO
set MEllbar2JetJetNLOInclusive:BornScreening Off
-cd /Herwig/MatrixElements/Matchbox
-mkdir LeptonHadron
-cd LeptonHadron
-
################################################################################
# DIS
################################################################################
create Herwig::MatchboxMElP2lJet MElq2lq
set MElq2lq:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MElq2lq:LeptonFlavours[0] /Herwig/Particles/e+
insert MElq2lq:LeptonFlavours[1] /Herwig/Particles/e-
insert MElq2lq:QuarkFlavours[0] /Herwig/Particles/d
insert MElq2lq:QuarkFlavours[1] /Herwig/Particles/u
insert MElq2lq:QuarkFlavours[2] /Herwig/Particles/s
insert MElq2lq:QuarkFlavours[3] /Herwig/Particles/c
+insert MElq2lq:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMElP2lJet MElqbar2lqbar
set MElqbar2lqbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MElqbar2lqbar:LeptonFlavours[0] /Herwig/Particles/e+
insert MElqbar2lqbar:LeptonFlavours[1] /Herwig/Particles/e-
insert MElqbar2lqbar:QuarkFlavours[0] /Herwig/Particles/dbar
insert MElqbar2lqbar:QuarkFlavours[1] /Herwig/Particles/ubar
insert MElqbar2lqbar:QuarkFlavours[2] /Herwig/Particles/sbar
insert MElqbar2lqbar:QuarkFlavours[3] /Herwig/Particles/cbar
+insert MElqbar2lqbar:QuarkFlavours[4] /Herwig/Particles/bbar
create Herwig::MatchboxMElq2lqg MElq2lqg
+set MElq2lqg:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MElq2lqg:LeptonFlavours[0] /Herwig/Particles/e+
insert MElq2lqg:LeptonFlavours[1] /Herwig/Particles/e-
insert MElq2lqg:QuarkFlavours[0] /Herwig/Particles/d
insert MElq2lqg:QuarkFlavours[1] /Herwig/Particles/u
insert MElq2lqg:QuarkFlavours[2] /Herwig/Particles/s
insert MElq2lqg:QuarkFlavours[3] /Herwig/Particles/c
+insert MElq2lqg:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMElq2lqg MElqbar2lqbarg
+set MElqbar2lqbarg:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MElqbar2lqbarg:LeptonFlavours[0] /Herwig/Particles/e+
insert MElqbar2lqbarg:LeptonFlavours[1] /Herwig/Particles/e-
insert MElqbar2lqbarg:QuarkFlavours[0] /Herwig/Particles/dbar
insert MElqbar2lqbarg:QuarkFlavours[1] /Herwig/Particles/ubar
insert MElqbar2lqbarg:QuarkFlavours[2] /Herwig/Particles/sbar
insert MElqbar2lqbarg:QuarkFlavours[3] /Herwig/Particles/cbar
+insert MElqbar2lqbarg:QuarkFlavours[4] /Herwig/Particles/bbar
create Herwig::MatchboxMElg2lqqbar MElg2lqqbar
+set MElg2lqqbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MElg2lqqbar:LeptonFlavours[0] /Herwig/Particles/e+
insert MElg2lqqbar:LeptonFlavours[1] /Herwig/Particles/e-
insert MElg2lqqbar:QuarkFlavours[0] /Herwig/Particles/d
insert MElg2lqqbar:QuarkFlavours[1] /Herwig/Particles/u
insert MElg2lqqbar:QuarkFlavours[2] /Herwig/Particles/s
insert MElg2lqqbar:QuarkFlavours[3] /Herwig/Particles/c
+insert MElg2lqqbar:QuarkFlavours[4] /Herwig/Particles/b
#-------------------------------------------------------------------------------
create ThePEG::SubProcessHandler MElP2lJetLO
set MElP2lJetLO:PartonExtractor /Herwig/Partons/DISExtractor
set MElP2lJetLO:Cuts /Herwig/Cuts/DISCuts
insert MElP2lJetLO:MatrixElements[0] MElq2lq
insert MElP2lJetLO:MatrixElements[1] MElqbar2lqbar
create Herwig::MatchboxFactory MElP2lJetNLO
set MElP2lJetNLO:PartonExtractor /Herwig/Partons/DISExtractor
set MElP2lJetNLO:Cuts /Herwig/Cuts/DISCuts
insert MElP2lJetNLO:BornMEs[0] MElq2lq
insert MElP2lJetNLO:BornMEs[1] MElqbar2lqbar
insert MElP2lJetNLO:RealEmissionMEs[0] MElq2lqg
insert MElP2lJetNLO:RealEmissionMEs[1] MElqbar2lqbarg
insert MElP2lJetNLO:RealEmissionMEs[2] MElg2lqqbar
create Herwig::PowhegFactory MElP2lJetNLOInclusive
set MElP2lJetNLOInclusive:PartonExtractor /Herwig/Partons/DISExtractor
set MElP2lJetNLOInclusive:Cuts /Herwig/Cuts/DISCuts
set MElP2lJetNLOInclusive:MatchboxFactory MElP2lJetNLO
-cd /Herwig/MatrixElements/Matchbox
-mkdir HadronHadron
-cd HadronHadron
-
################################################################################
# DY
################################################################################
create Herwig::MatchboxMEPP2llbar MEqqbar2llbar
set MEqqbar2llbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MEqqbar2llbar:LeptonFlavours[0] /Herwig/Particles/e-
insert MEqqbar2llbar:QuarkFlavours[0] /Herwig/Particles/d
insert MEqqbar2llbar:QuarkFlavours[1] /Herwig/Particles/u
insert MEqqbar2llbar:QuarkFlavours[2] /Herwig/Particles/s
insert MEqqbar2llbar:QuarkFlavours[3] /Herwig/Particles/c
+insert MEqqbar2llbar:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMEPP2llbar MEqbarq2llbar
set MEqbarq2llbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MEqbarq2llbar:LeptonFlavours[0] /Herwig/Particles/e-
insert MEqbarq2llbar:QuarkFlavours[0] /Herwig/Particles/dbar
insert MEqbarq2llbar:QuarkFlavours[1] /Herwig/Particles/ubar
insert MEqbarq2llbar:QuarkFlavours[2] /Herwig/Particles/sbar
insert MEqbarq2llbar:QuarkFlavours[3] /Herwig/Particles/cbar
+insert MEqbarq2llbar:QuarkFlavours[4] /Herwig/Particles/bbar
create Herwig::MatchboxMEqqbar2llbarg MEqqbar2llbarg
+set MEqqbar2llbarg:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MEqqbar2llbarg:LeptonFlavours[0] /Herwig/Particles/e-
insert MEqqbar2llbarg:QuarkFlavours[0] /Herwig/Particles/d
insert MEqqbar2llbarg:QuarkFlavours[1] /Herwig/Particles/u
insert MEqqbar2llbarg:QuarkFlavours[2] /Herwig/Particles/s
insert MEqqbar2llbarg:QuarkFlavours[3] /Herwig/Particles/c
+insert MEqqbar2llbarg:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMEqqbar2llbarg MEqbarq2llbarg
+set MEqbarq2llbarg:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
insert MEqbarq2llbarg:LeptonFlavours[0] /Herwig/Particles/e-
insert MEqbarq2llbarg:QuarkFlavours[0] /Herwig/Particles/dbar
insert MEqbarq2llbarg:QuarkFlavours[1] /Herwig/Particles/ubar
insert MEqbarq2llbarg:QuarkFlavours[2] /Herwig/Particles/sbar
insert MEqbarq2llbarg:QuarkFlavours[3] /Herwig/Particles/cbar
+insert MEqbarq2llbarg:QuarkFlavours[4] /Herwig/Particles/bbar
create Herwig::MatchboxMEqg2llbarq MEqg2llbarq
+set MEqg2llbarq:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
set MEqg2llbarq:WhichGluon 1
insert MEqg2llbarq:LeptonFlavours[0] /Herwig/Particles/e-
insert MEqg2llbarq:QuarkFlavours[0] /Herwig/Particles/d
insert MEqg2llbarq:QuarkFlavours[1] /Herwig/Particles/u
insert MEqg2llbarq:QuarkFlavours[2] /Herwig/Particles/s
insert MEqg2llbarq:QuarkFlavours[3] /Herwig/Particles/c
+insert MEqg2llbarq:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMEqg2llbarq MEgq2llbarq
+set MEgq2llbarq:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
set MEgq2llbarq:WhichGluon 0
insert MEgq2llbarq:LeptonFlavours[0] /Herwig/Particles/e-
insert MEgq2llbarq:QuarkFlavours[0] /Herwig/Particles/d
insert MEgq2llbarq:QuarkFlavours[1] /Herwig/Particles/u
insert MEgq2llbarq:QuarkFlavours[2] /Herwig/Particles/s
insert MEgq2llbarq:QuarkFlavours[3] /Herwig/Particles/c
+insert MEgq2llbarq:QuarkFlavours[4] /Herwig/Particles/b
create Herwig::MatchboxMEqg2llbarq MEqbarg2llbarqbar
+set MEqbarg2llbarqbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
set MEqbarg2llbarqbar:WhichGluon 1
insert MEqbarg2llbarqbar:LeptonFlavours[0] /Herwig/Particles/e-
insert MEqbarg2llbarqbar:QuarkFlavours[0] /Herwig/Particles/dbar
insert MEqbarg2llbarqbar:QuarkFlavours[1] /Herwig/Particles/ubar
insert MEqbarg2llbarqbar:QuarkFlavours[2] /Herwig/Particles/sbar
insert MEqbarg2llbarqbar:QuarkFlavours[3] /Herwig/Particles/cbar
+insert MEqbarg2llbarqbar:QuarkFlavours[4] /Herwig/Particles/bbar
create Herwig::MatchboxMEqg2llbarq MEgqbar2llbarqbar
+set MEgqbar2llbarqbar:Phasespace /Herwig/MatrixElements/Matchbox/Phasespace/TreePhasespace
set MEgqbar2llbarqbar:WhichGluon 0
insert MEgqbar2llbarqbar:LeptonFlavours[0] /Herwig/Particles/e-
insert MEgqbar2llbarqbar:QuarkFlavours[0] /Herwig/Particles/dbar
insert MEgqbar2llbarqbar:QuarkFlavours[1] /Herwig/Particles/ubar
insert MEgqbar2llbarqbar:QuarkFlavours[2] /Herwig/Particles/sbar
insert MEgqbar2llbarqbar:QuarkFlavours[3] /Herwig/Particles/cbar
+insert MEgqbar2llbarqbar:QuarkFlavours[4] /Herwig/Particles/bbar
#-------------------------------------------------------------------------------
create ThePEG::SubProcessHandler MEPP2llbarLO
set MEPP2llbarLO:PartonExtractor /Herwig/Partons/QCDExtractor
set MEPP2llbarLO:Cuts /Herwig/Cuts/QCDCuts
insert MEPP2llbarLO:MatrixElements[0] MEqqbar2llbar
insert MEPP2llbarLO:MatrixElements[1] MEqbarq2llbar
create Herwig::MatchboxFactory MEPP2llbarNLO
set MEPP2llbarNLO:PartonExtractor /Herwig/Partons/QCDExtractor
set MEPP2llbarNLO:Cuts /Herwig/Cuts/QCDCuts
insert MEPP2llbarNLO:BornMEs[0] MEqqbar2llbar
insert MEPP2llbarNLO:BornMEs[1] MEqbarq2llbar
insert MEPP2llbarNLO:RealEmissionMEs[0] MEqqbar2llbarg
insert MEPP2llbarNLO:RealEmissionMEs[1] MEqbarq2llbarg
insert MEPP2llbarNLO:RealEmissionMEs[2] MEqg2llbarq
insert MEPP2llbarNLO:RealEmissionMEs[3] MEgq2llbarq
insert MEPP2llbarNLO:RealEmissionMEs[4] MEqbarg2llbarqbar
insert MEPP2llbarNLO:RealEmissionMEs[5] MEgqbar2llbarqbar
create Herwig::PowhegFactory MEPP2llbarNLOInclusive
set MEPP2llbarNLOInclusive:PartonExtractor /Herwig/Partons/QCDExtractor
set MEPP2llbarNLOInclusive:Cuts /Herwig/Cuts/QCDCuts
set MEPP2llbarNLOInclusive:MatchboxFactory MEPP2llbarNLO
cd /

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:53 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805459
Default Alt Text
(262 KB)

Event Timeline