Page MenuHomeHEPForge
Files F8310538

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
@@ -1,1330 +1,1331 @@
// -*- C++ -*-
//
// SubtractionDipole.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 SubtractionDipole class.
//
#include "SubtractionDipole.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
SubtractionDipole::SubtractionDipole()
: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
theIgnoreCuts(false),
theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
lastRealEmissionKey(realEmissionKey(cPDVector(),-1,-1,-1)),
lastUnderlyingBornKey(underlyingBornKey(cPDVector(),-1,-1)),
theBornEmitter(-1), theBornSpectator(-1),
theLastSubtractionScale(ZERO), theLastSplittingScale(ZERO),
theLastSubtractionPt(ZERO), theLastSplittingPt(ZERO),
theLastSubtractionZ(0.0), theLastSplittingZ(0.0),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
theLoopSimSubtraction(false), theRealEmissionScales(false),
theShowerHardScale(ZERO), theShowerScale(ZERO),
theIsInShowerPhasespace(false), theIsAboveCutoff(false) {}
SubtractionDipole::~SubtractionDipole() {}
double SubtractionDipole::alpha() const{
return factory()->alphaParameter();
}
void SubtractionDipole::clearBookkeeping() {
theRealEmitter = -1;
theRealEmission = -1;
theRealSpectator = -1;
theBornEmitter = -1;
theBornSpectator = -1;
theMergingMap.clear();
theSplittingMap.clear();
theIndexMap.clear();
theUnderlyingBornDiagrams.clear();
theRealEmissionDiagrams.clear();
theBornToRealDiagrams.clear();
theRealToBornDiagrams.clear();
}
void SubtractionDipole::setupBookkeeping(const map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo>& mergeInfo,bool slim) {
theMergingMap.clear();
theSplittingMap.clear();
theUnderlyingBornDiagrams.clear();
theRealEmissionDiagrams.clear();
theBornToRealDiagrams.clear();
theRealToBornDiagrams.clear();
int xemitter = -1;
int xspectator = -1;
map<int,int> mergeLegs;
map<int,int> remapLegs;
map<int,int> realBornMap;
map<int,int> bornRealMap;
set<Ptr<DiagramBase>::cptr> usedDiagrams;
for ( map<Ptr<DiagramBase>::ptr,MergeInfo>::const_iterator mit = mergeInfo.begin();
mit != mergeInfo.end(); ++mit ) {
DiagramVector::const_iterator bd =
theUnderlyingBornME->diagrams().end();
// work out the most similar underlying Born diagram
map<int,int> xRemapLegs;
int nomapScore = 0;
for ( DiagramVector::const_iterator b =
theUnderlyingBornME->diagrams().begin();
b != theUnderlyingBornME->diagrams().end(); ++b ) {
map<int,int> theRemapLegs;
if ( mit->second.diagram->isSame(*b,theRemapLegs) &&
usedDiagrams.find(*b) == usedDiagrams.end() ) {
int theNomapScore = 0;
for ( map<int,int>::const_iterator m = theRemapLegs.begin();
m != theRemapLegs.end(); ++m )
if ( m->first == m->second )
theNomapScore += 1;
if ( theNomapScore >= nomapScore ) {
nomapScore = theNomapScore;
xRemapLegs = theRemapLegs;
bd = b;
}
}
}
// no underlying Born
if ( bd == theUnderlyingBornME->diagrams().end() )
continue;
// as we deal with one splitting only we now mark this diagram as used
// since we fixed the overall remapping of the process from the first
// occurence, see below. TODO: This confuses this code even more, and
// clearly calls for a cleanup. This is just grown historically and got
// messed up with experiencing different processes and setups.
usedDiagrams.insert(*bd);
if ( xemitter == -1 ) {
xemitter = mit->second.emitter;
mergeLegs = mit->second.mergeLegs;
remapLegs = xRemapLegs;
assert(remapLegs.find(xemitter) != remapLegs.end());
xemitter = remapLegs[xemitter];
// work out the leg remapping real -> born
for ( map<int,int>::const_iterator k = mergeLegs.begin();
k != mergeLegs.end(); ++k ) {
assert(remapLegs.find(k->second) != remapLegs.end());
realBornMap[k->first] = remapLegs[k->second];
}
// work out the leg remapping born -> real
for ( map<int,int>::const_iterator k = realBornMap.begin();
k != realBornMap.end(); ++k ) {
bornRealMap[k->second] = k->first;
}
// work out the spectator
assert(mergeLegs.find(realSpectator()) != mergeLegs.end());
assert(remapLegs.find(mergeLegs[realSpectator()]) != remapLegs.end());
xspectator = realBornMap[realSpectator()];
}
RealEmissionKey realKey = realEmissionKey((*mit->first).partons(),realEmitter(),realEmission(),realSpectator());
UnderlyingBornKey bornKey = underlyingBornKey((**bd).partons(),xemitter,xspectator);
if ( theMergingMap.find(realKey) == theMergingMap.end() )
theMergingMap.insert(make_pair(realKey,make_pair(bornKey,realBornMap)));
RealEmissionInfo realInfo = make_pair(realKey,bornRealMap);
bool gotit = false;
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spIterator;
pair<spIterator,spIterator> range = theSplittingMap.equal_range(bornKey);
for ( ; range.first != range.second; ++range.first )
if ( range.first->second == realInfo ) {
gotit = true;
break;
}
if ( !gotit )
theSplittingMap.insert(make_pair(bornKey,realInfo));
theUnderlyingBornDiagrams[process(realKey)].push_back(*bd);
theRealEmissionDiagrams[process(bornKey)].push_back(mit->first);
theBornToRealDiagrams[*bd] = mit->first;
theRealToBornDiagrams[mit->first] = *bd;
}
if (slim) {
theIndexMap.clear();
theSplittingMap.clear();
theBornToRealDiagrams.clear();
theRealEmissionDiagrams.clear();
}
if ( theSplittingMap.empty() )
return;
theIndexMap.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator s =
theSplittingMap.begin(); s != theSplittingMap.end(); ++s ) {
theIndexMap[process(s->first)] = make_pair(emitter(s->first),spectator(s->first));
}
}
void SubtractionDipole::subtractionBookkeeping() {
/*
if ( theMergingMap.empty() )
setupBookkeeping();
*/
assert(!theMergingMap.empty());
lastRealEmissionKey =
realEmissionKey(lastHeadXComb().mePartonData(),realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() ) {
theApply = false;
return;
}
theApply = true;
lastUnderlyingBornKey = k->second.first;
bornEmitter(emitter(lastUnderlyingBornKey));
bornSpectator(spectator(lastUnderlyingBornKey));
}
void SubtractionDipole::splittingBookkeeping() {
/*
if ( theMergingMap.empty() )
setupBookkeeping();
*/
assert(!theMergingMap.empty());
map<cPDVector,pair<int,int> >::const_iterator esit =
theIndexMap.find(lastHeadXComb().mePartonData());
if ( esit == theIndexMap.end() ) {
theApply = false;
return;
}
theApply = true;
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(lastHeadXComb().mePartonData(),bornEmitter(),bornSpectator());
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
assert(kr.first != kr.second);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
assert(lastRealEmissionInfo != kr.second);
lastRealEmissionKey = lastRealEmissionInfo->second.first;
realEmitter(emitter(lastRealEmissionKey));
realEmission(emission(lastRealEmissionKey));
realSpectator(spectator(lastRealEmissionKey));
}
StdXCombPtr SubtractionDipole::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allBins,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
return
realEmissionME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
StdXCombPtr SubtractionDipole::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
return
realEmissionME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
StdXCombPtr SubtractionDipole::makeBornXComb(tStdXCombPtr realXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();
lastRealEmissionKey =
realEmissionKey(proc,realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() )
return StdXCombPtr();
PartonPairVec pbs = realXC->pExtractor()->getPartons(realXC->maxEnergy(),
realXC->particles(),
*(realXC->cuts()));
DiagramVector bornDiags = underlyingBornDiagrams(proc);
assert(!bornDiags.empty());
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == bornDiags.front()->partons()[0] &&
ppit->second->parton() == bornDiags.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
return
underlyingBornME()->makeXComb(realXC,*ppit,bornDiags,this);
}
vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();
map<cPDVector,pair<int,int> >::const_iterator esit = theIndexMap.find(proc);
if ( esit == theIndexMap.end() )
return vector<StdXCombPtr>();
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(proc,bornEmitter(),bornSpectator());
if ( theSplittingMap.find(lastUnderlyingBornKey) == theSplittingMap.end() )
return vector<StdXCombPtr>();
PartonPairVec pbs = bornXC->pExtractor()->getPartons(bornXC->maxEnergy(),
bornXC->particles(),
*(bornXC->cuts()));
DiagramVector realDiags = realEmissionDiagrams(proc);
assert(!realDiags.empty());
vector<StdXCombPtr> res;
map<cPDVector,DiagramVector> realProcs;
for ( MEBase::DiagramVector::const_iterator d = realDiags.begin();
d != realDiags.end(); ++d ) {
realProcs[(**d).partons()].push_back(*d);
}
for ( map<cPDVector,DiagramVector>::const_iterator pr =
realProcs.begin(); pr != realProcs.end(); ++pr ) {
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == pr->second.front()->partons()[0] &&
ppit->second->parton() == pr->second.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
StdXCombPtr rxc =
realEmissionME()->makeXComb(bornXC,*ppit,pr->second,this);
res.push_back(rxc);
}
return res;
}
const MEBase::DiagramVector& SubtractionDipole::underlyingBornDiagrams(const cPDVector& real) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theUnderlyingBornDiagrams.find(real);
if (k == theUnderlyingBornDiagrams.end() )
return empty;
return k->second;
}
tcDiagPtr SubtractionDipole::underlyingBornDiagram(tcDiagPtr realDiag) const {
map<tcDiagPtr,tcDiagPtr>::const_iterator it = theRealToBornDiagrams.find(realDiag);
assert(it != theRealToBornDiagrams.end());
return it->second;
}
const MEBase::DiagramVector& SubtractionDipole::realEmissionDiagrams(const cPDVector& born) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theRealEmissionDiagrams.find(born);
if ( k == theRealEmissionDiagrams.end() )
return empty;
return k->second;
}
tcDiagPtr SubtractionDipole::realEmissionDiagram(tcDiagPtr bornDiag) const {
map<tcDiagPtr,tcDiagPtr>::const_iterator it = theBornToRealDiagrams.find(bornDiag);
assert(it != theBornToRealDiagrams.end());
return it->second;
}
void SubtractionDipole::getDiagrams() const {
if ( splitting() ) {
realEmissionME()->diagrams();
useDiagrams(realEmissionME());
} else {
underlyingBornME()->diagrams();
useDiagrams(underlyingBornME());
}
}
Selector<MEBase::DiagramIndex> SubtractionDipole::diagrams(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->setXComb(lastXCombPtr());
+ me->phasespace()->clearDiagramWeights();
me->phasespace()->fillDiagramWeights();
}
return
me->diagrams(dv);
}
Selector<const ColourLines *>
SubtractionDipole::colourGeometries(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->colourGeometries(diag) :
underlyingBornME()->colourGeometries(diag);
}
const ColourLines &
SubtractionDipole::selectColourGeometry(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->selectColourGeometry(diag) :
underlyingBornME()->selectColourGeometry(diag);
}
void SubtractionDipole::flushCaches() {
theUnderlyingBornME->flushCaches();
theRealEmissionME->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
}
void SubtractionDipole::setXComb(tStdXCombPtr xc) {
if ( !xc ) {
theApply = false;
return;
} else {
theApply = true;
}
lastMatchboxXComb(xc);
MEBase::setXComb(xc);
if ( splitting() ) {
realEmissionME()->setXComb(xc);
underlyingBornME()->setXComb(xc->head());
splittingBookkeeping();
} else {
realEmissionME()->setXComb(xc->head());
underlyingBornME()->setXComb(xc);
subtractionBookkeeping();
}
if ( !apply() )
return;
}
void SubtractionDipole::setKinematics() {
MEBase::setKinematics();
if ( splitting() )
realEmissionME()->setKinematics();
else
underlyingBornME()->setKinematics();
}
bool SubtractionDipole::generateKinematics(const double * r) {
if ( lastXCombPtr()->kinematicsGenerated() )
return true;
if ( splitting() ) {
if ( !generateRadiationKinematics(r) )
return false;
if( ! realEmissionME()->lastXCombPtr()->setIncomingPartons())
return false;
realEmissionME()->setScale();
double jac = jacobian();
jac *= pow(underlyingBornME()->lastXComb().lastSHat() / realEmissionME()->lastXComb().lastSHat(),
realEmissionME()->lastXComb().mePartonData().size()-4.);
jacobian(jac);
assert(lastXCombPtr() == realEmissionME()->lastXCombPtr());
lastXCombPtr()->didGenerateKinematics();
return true;
}
if ( !generateTildeKinematics() ){ return false;}
if( ! underlyingBornME()->lastXCombPtr()->setIncomingPartons() )
return false;
underlyingBornME()->setScale();
assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
if( ! underlyingBornME()->lastXCombPtr()->setIncomingPartons() )
return false;
// need to have the scale and x's available for checking shower phase space
if ( showerApproximation() &&
lastXCombPtr()->willPassCuts() )
showerApproximation()->getShowerVariables();
lastXCombPtr()->didGenerateKinematics();
return true;
}
int SubtractionDipole::nDim() const {
if ( !splitting() )
return underlyingBornME()->nDim();
return underlyingBornME()->nDim() + nDimRadiation();
}
void SubtractionDipole::clearKinematics() {
MEBase::clearKinematics();
if ( splitting() )
realEmissionME()->clearKinematics();
else
underlyingBornME()->clearKinematics();
}
void SubtractionDipole::tildeKinematics(Ptr<TildeKinematics>::tptr tk) {
theTildeKinematics = tk;
}
bool SubtractionDipole::generateTildeKinematics() {
assert(!splitting());
Ptr<TildeKinematics>::tptr kinematics = theTildeKinematics;
if ( showerApproximation() ) {
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
showerApproximation()->checkCutoff();
if ( showerApproximation()->showerTildeKinematics() &&
isAboveCutoff() &&
realShowerSubtraction() )
kinematics = showerApproximation()->showerTildeKinematics();
}
if ( !kinematics ) {
jacobian(0.0);
return false;
}
kinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());
if ( !kinematics->doMap() ) {
jacobian(0.0);
return false;
}
theLastSubtractionScale = kinematics->lastScale();
theLastSubtractionPt = kinematics->lastPt();
theLastSubtractionZ = kinematics->lastZ();
meMomenta().resize(lastHeadXComb().meMomenta().size() - 1);
assert(mergingMap().find(lastRealEmissionKey) != mergingMap().end());
map<int,int>& trans = theMergingMap[lastRealEmissionKey].second;
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == realEmitter() || k == realEmission() || k == realSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( kinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[bornEmitter()] =
const_cast<const TildeKinematics&>(*kinematics).bornEmitterMomentum();
meMomenta()[bornSpectator()] =
const_cast<const TildeKinematics&>(*kinematics).bornSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).hardProcessMass());
p->rescaleRho();
}
jacobian(realEmissionME()->lastXComb().jacobian());
logGenerateTildeKinematics();
return true;
}
void SubtractionDipole::invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr itk) {
theInvertedTildeKinematics = itk;
}
int SubtractionDipole::nDimRadiation() const {
return invertedTildeKinematics() ?
invertedTildeKinematics()->nDimRadiation() :
0;
}
bool SubtractionDipole::generateRadiationKinematics(const double * r) {
assert(splitting());
Ptr<InvertedTildeKinematics>::tptr kinematics = theInvertedTildeKinematics;
if ( showerApproximation() ) {
showerApproximation()->setBornXComb(lastHeadXCombPtr());
showerApproximation()->setRealXComb(lastXCombPtr());
showerApproximation()->setDipole(this);
if ( showerApproximation()->showerInvertedTildeKinematics() ) {
kinematics = showerApproximation()->showerInvertedTildeKinematics();
}
}
if ( !kinematics ) {
jacobian(0.0);
return false;
}
kinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());
if ( !kinematics->doMap(r) ) {
jacobian(0.0);
return false;
}
theLastSplittingScale = kinematics->lastScale();
theLastSplittingPt = kinematics->lastPt();
theLastSplittingZ = kinematics->lastZ();
meMomenta().resize(lastHeadXComb().meMomenta().size() + 1);
assert(splittingMap().find(lastUnderlyingBornKey) != splittingMap().end());
map<int,int>& trans = const_cast<map<int,int>&>(lastRealEmissionInfo->second.second);
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == bornEmitter() || k == bornSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( kinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[realEmitter()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realEmitterMomentum();
meMomenta()[realEmission()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realEmissionMomentum();
meMomenta()[realSpectator()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).hardProcessMass());
p->rescaleRho();
}
jacobian(underlyingBornME()->lastXComb().jacobian() *
kinematics->jacobian());
logGenerateRadiationKinematics(r);
return true;
}
void SubtractionDipole::ptCut(Energy cut) {
theInvertedTildeKinematics->ptCut(cut);
}
CrossSection SubtractionDipole::dSigHatDR(Energy2 factorizationScale) const {
double pdfweight = 1.;
double jac = jacobian();
if ( splitting() && jac == 0.0 ) {
lastMECrossSection(ZERO);
return ZERO;
}
if ( factorizationScale == ZERO ) {
factorizationScale = underlyingBornME()->lastScale();
}
if ( havePDFWeight1() ) {
pdfweight *= realEmissionME()->pdf1(factorizationScale);
}
if ( havePDFWeight2() ) {
pdfweight *= realEmissionME()->pdf2(factorizationScale);
}
lastMEPDFWeight(pdfweight);
bool needTheDipole = true;
CrossSection shower = ZERO;
double lastThetaMu = 1.0;
double showerFactor = 1.;
if ( showerApproximation() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(const_cast<SubtractionDipole*>(this));
if ( !isAboveCutoff() ) {
showerApproximation()->wasBelowCutoff();
lastThetaMu = 0.0;
} else {
lastThetaMu = 1.0;
}
if ( lastThetaMu > 0.0 && isInShowerPhasespace() ) {
if ( realShowerSubtraction() )
shower = showerApproximation()->dSigHatDR()*lastThetaMu;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
shower = -showerApproximation()->dSigHatDR()*lastThetaMu;
if ( virtualShowerSubtraction() &&
isAboveCutoff() &&
showerApproximation()->showerTildeKinematics() ) {
// map shower to dipole kinematics; we are always above the
// cutoff in this case
showerFactor *=
showerApproximation()->showerTildeKinematics()->jacobianRatio();
}
shower *= showerFactor;
}
if ( realShowerSubtraction() && lastThetaMu == 1.0 )
needTheDipole = false;
if ( virtualShowerSubtraction() && lastThetaMu == 0.0 )
needTheDipole = false;
if ( factory()->loopSimCorrections() ||
factory()->meCorrectionsOnly() )
needTheDipole = false;
}
double xme2 = 0.0;
if ( needTheDipole )
xme2 = me2();
if ( factory()->loopSimCorrections() ||
factory()->meCorrectionsOnly() ) {
assert(showerApproximation());
xme2 = realEmissionME()->me2() * showerApproximation()->channelWeight();
double rws =
pow(underlyingBornME()->lastXComb().lastAlphaS()/
realEmissionME()->lastXComb().lastAlphaS(),
realEmissionME()->orderInAlphaS());
xme2 *= rws;
double rwe =
pow(underlyingBornME()->lastXComb().lastAlphaEM()/
realEmissionME()->lastXComb().lastAlphaEM(),
underlyingBornME()->orderInAlphaEW());
xme2 *= rwe;
}
if ( realShowerSubtraction() )
xme2 *= 1. - lastThetaMu;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
xme2 *= lastThetaMu;
double coupl = lastMECouplings();
coupl *= underlyingBornME()->lastXComb().lastAlphaS();
lastMECouplings(coupl);
CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
if ( !showerApproximation() && xme2 != 0.0 ) {
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(theRealEmissionME->lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
}
lastMECrossSection(-res-shower);
logDSigHatDR(jac);
return lastMECrossSection();
}
bool SubtractionDipole::aboveAlpha() const{return theTildeKinematics->aboveAlpha();}
CrossSection SubtractionDipole::prefactor(Energy2 factorizationScale)const{
const double jac = jacobian();
assert( factorizationScale != ZERO );
assert (! splitting());
double pdfweight = 1.;
if ( havePDFWeight1() ) pdfweight *= realEmissionME()->pdf1(factorizationScale);
if ( havePDFWeight2() ) pdfweight *= realEmissionME()->pdf2(factorizationScale);
return sqr(hbarc) * jac * pdfweight / (2. * realEmissionME()->lastXComb().lastSHat());
}
CrossSection SubtractionDipole::ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
double ccme2 =underlyingBornME()->me2()*
underlyingBornME()->
largeNColourCorrelatedME2(
make_pair(bornEmitter(),bornSpectator()),largeNBasis)/
underlyingBornME()->largeNME2(largeNBasis);
return prefactor(factorizationScale) * me2Avg(ccme2);
}
pair<CrossSection,CrossSection> SubtractionDipole::dipandPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
CrossSection factor= prefactor(factorizationScale);
double ccme2 =underlyingBornME()->me2()*
underlyingBornME()->
largeNColourCorrelatedME2(
make_pair(bornEmitter(),bornSpectator()),largeNBasis)/
underlyingBornME()->largeNME2(largeNBasis);
double ps = me2Avg(ccme2);
double dip = me2();
return make_pair(factor*dip,factor*ps);
}
CrossSection SubtractionDipole::dip(Energy2 factorizationScale) const {
CrossSection factor= prefactor(factorizationScale);
double dip = me2();
return factor*dip;
}
void SubtractionDipole::print(ostream& os) const {
os << "--- SubtractionDipole setup ----------------------------------------------------\n";
os << " subtraction '" << name() << "'\n for real emission '"
<< theRealEmissionME->name() << "'\n using underlying Born '"
<< theUnderlyingBornME->name() << "'\n";
os << " tilde kinematics are '"
<< (theTildeKinematics ? theTildeKinematics->name() : "")
<< " '\n inverted tilde kinematics are '"
<< (theInvertedTildeKinematics ? theInvertedTildeKinematics->name() : "") << "'\n";
os << " the following subtraction mappings have been found:\n";
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator m =
theMergingMap.begin(); m != theMergingMap.end(); ++m ) {
os << " " << process(m->second.first)[0]->PDGName() << " "
<< process(m->second.first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->second.first).begin() + 2;
p != process(m->second.first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[" << emitter(m->second.first) << "," << spectator(m->second.first) << "] <=> ";
os << process(m->first)[0]->PDGName() << " "
<< process(m->first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->first).begin() + 2;
p != process(m->first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[(" << emitter(m->first) << "," << emission(m->first) << ")," << spectator(m->first) << "]\n"
<< " non-dipole momenta ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->second << " ";
}
os << ") <=> ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->first << " ";
}
os << ")\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractionDipole::printLastEvent(ostream& os) const {
os << "--- SubtractionDipole last event information -----------------------------------\n";
os << " for dipole '" << name() << "' applying ["
<< bornEmitter() << "," << bornSpectator() << "] <=> [("
<< realEmitter() << "," << realEmission() << ")," << realSpectator() << "]\n"
<< " evaluated the cross section/nb " << (lastMECrossSection()/nanobarn) << "\n"
<< " with subtraction parameters x[0] = " << subtractionParameters()[0]
<< " x[1] = " << subtractionParameters()[1] << "\n";
os << " the last real emission event was:\n";
realEmissionME()->printLastEvent(os);
os << " the last underlying Born event was:\n";
underlyingBornME()->printLastEvent(os);
os << "--- end SubtractionDipole last event information -------------------------------\n";
os << flush;
}
void SubtractionDipole::logME2() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated me2 using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "Born phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = bornxc->meMomenta().begin();
cPDVector::const_iterator dit = bornxc->mePartonData().begin();
for ( ; pit != bornxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << bornxc->lastX1() << " x2 = " << bornxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (bornxc->lastSHat()/GeV2) << "\n";
generator()->log() << "Real emission phase space point (in GeV):\n";
pit = realxc->meMomenta().begin();
dit = realxc->mePartonData().begin();
for ( ; pit != realxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << realxc->lastX1() << " x2 = " << realxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (realxc->lastSHat()/GeV2) << "\n";
}
void SubtractionDipole::logDSigHatDR(double effectiveJac) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated cross section using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n"
<< "Jacobian = " << jacobian()
<< " effective Jacobian = " << effectiveJac << "\n"
<< "Born sHat/GeV2 = " << (bornxc->lastSHat()/GeV2)
<< " real sHat/GeV2 = " << (realxc->lastSHat()/GeV2)
<< " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void SubtractionDipole::logGenerateTildeKinematics() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating tilde kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with real xcomb " << lastHeadXCombPtr() << " born xcomb "
<< lastXCombPtr() << "\n"
<< "from real emission phase space point:\n";
Lorentz5Momentum rSum;
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
size_t count = 0;
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
rSum -= *pr;
} else {
rSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (rSum/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n"
<< "with scale/GeV = " << (theLastSubtractionScale/GeV)
<< "and pt/GeV = " << (theLastSubtractionPt/GeV) << "\n";
generator()->log() << "generated tilde kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
count = 0;
Lorentz5Momentum bSum;
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
bSum -= *pr;
} else {
bSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (bSum/GeV) << "\n";
generator()->log() << "Jacobian = " << jacobian() << "\n" << flush;
}
void SubtractionDipole::logGenerateRadiationKinematics(const double * r) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating radiation kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with born xcomb " << lastHeadXCombPtr() << " real xcomb "
<< lastXCombPtr() << "\n"
<< "from random numbers:\n";
copy(r,r+nDimRadiation(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "and born phase space point:\n";
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n" << flush;
generator()->log() << "scales: scale/GeV = " << (theLastSplittingScale/GeV)
<< " pt/GeV = " << (theLastSplittingPt/GeV) << "\n" << flush;
generator()->log() << "generated real emission kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "Jacobian = "
<< jacobian() << " = "
<< underlyingBornME()->lastXComb().jacobian()
<< "|Born * "
<< invertedTildeKinematics()->jacobian()
<< "|Radiation\n" << flush;
}
void SubtractionDipole::doinit() {
MEBase::doinit();
if ( underlyingBornME() ) {
theUnderlyingBornME->init();
}
if ( realEmissionME() ) {
theRealEmissionME->init();
}
if ( tildeKinematics() ) {
theTildeKinematics->init();
}
if ( invertedTildeKinematics() ) {
theInvertedTildeKinematics->init();
}
if ( showerApproximation() ) {
theShowerApproximation->init();
}
for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
p != thePartners.end(); ++p ) {
(**p).init();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).init();
}
}
void SubtractionDipole::doinitrun() {
MEBase::doinitrun();
if ( underlyingBornME() ) {
theUnderlyingBornME->initrun();
}
if ( realEmissionME() ) {
theRealEmissionME->initrun();
}
if ( tildeKinematics() ) {
theTildeKinematics->initrun();
}
if ( invertedTildeKinematics() ) {
theInvertedTildeKinematics->initrun();
}
if ( showerApproximation() ) {
theShowerApproximation->initrun();
}
for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
p != thePartners.end(); ++p ) {
(**p).initrun();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).initrun();
}
}
void SubtractionDipole::cloneDependencies(const std::string& prefix,bool slim) {
if ( underlyingBornME() ) {
Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
myUnderlyingBornME->cloneDependencies(pname.str(),slim);
underlyingBornME(myUnderlyingBornME);
}
if ( realEmissionME()&& !slim ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = realEmissionME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
myRealEmissionME->cloneDependencies(pname.str());
realEmissionME(myRealEmissionME);
}
if ( tildeKinematics() ) {
Ptr<TildeKinematics>::ptr myTildeKinematics = tildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myTildeKinematics->name();
if ( ! (generator()->preinitRegister(myTildeKinematics,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
myTildeKinematics->dipole(this);
tildeKinematics(myTildeKinematics);
}
if ( invertedTildeKinematics()&& !slim ) {
Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Inverted tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
myInvertedTildeKinematics->dipole(this);
invertedTildeKinematics(myInvertedTildeKinematics);
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror;
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
}
void SubtractionDipole::constructVertex(tSubProPtr sub) {
if ( splitting() )
realEmissionME()->constructVertex(sub);
else
underlyingBornME()->constructVertex(sub);
}
void SubtractionDipole::constructVertex(tSubProPtr sub, const ColourLines* cl) {
if ( splitting() )
realEmissionME()->constructVertex(sub,cl);
else
underlyingBornME()->constructVertex(sub,cl);
}
void SubtractionDipole::generateSubCollision(SubProcess & sub) {
if ( splitting() )
realEmissionME()->generateSubCollision(sub);
else
underlyingBornME()->generateSubCollision(sub);
}
void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theSplitting << theApply << theSubtractionTest
<< theIgnoreCuts << theRealEmissionME << theUnderlyingBornME
<< thePartners << theTildeKinematics << theInvertedTildeKinematics
<< theReweights << theRealEmitter << theRealEmission << theRealSpectator
<< theSubtractionParameters << theMergingMap << theSplittingMap
<< theIndexMap << theUnderlyingBornDiagrams << theRealEmissionDiagrams
<< theBornToRealDiagrams << theRealToBornDiagrams
<< lastRealEmissionKey << lastUnderlyingBornKey
<< theBornEmitter << theBornSpectator << ounit(theLastSubtractionScale,GeV)
<< ounit(theLastSplittingScale,GeV) << ounit(theLastSubtractionPt,GeV)
<< ounit(theLastSplittingPt,GeV) << theLastSubtractionZ
<< theLastSplittingZ << theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction
<< theLoopSimSubtraction << theRealEmissionScales << theFactory
<< ounit(theShowerHardScale,GeV) << ounit(theShowerScale,GeV)
<< theShowerParameters << theIsInShowerPhasespace << theIsAboveCutoff;
}
void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theSplitting >> theApply >> theSubtractionTest
>> theIgnoreCuts >> theRealEmissionME >> theUnderlyingBornME
>> thePartners >> theTildeKinematics >> theInvertedTildeKinematics
>> theReweights >> theRealEmitter >> theRealEmission >> theRealSpectator
>> theSubtractionParameters >> theMergingMap >> theSplittingMap
>> theIndexMap >> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
>> theBornToRealDiagrams >> theRealToBornDiagrams
>> lastRealEmissionKey >> lastUnderlyingBornKey
>> theBornEmitter >> theBornSpectator >> iunit(theLastSubtractionScale,GeV)
>> iunit(theLastSplittingScale,GeV) >> iunit(theLastSubtractionPt,GeV)
>> iunit(theLastSplittingPt,GeV) >> theLastSubtractionZ
>> theLastSplittingZ >> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction
>> theLoopSimSubtraction >> theRealEmissionScales >> theFactory
>> iunit(theShowerHardScale,GeV) >> iunit(theShowerScale,GeV)
>> theShowerParameters >> theIsInShowerPhasespace >> theIsAboveCutoff;
lastMatchboxXComb(theLastXComb);
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
}
Ptr<MatchboxFactory>::tptr SubtractionDipole::factory() const {
return theFactory;
}
void SubtractionDipole::factory(Ptr<MatchboxFactory>::tptr f) {
theFactory = f;
}
void SubtractionDipole::Init() {
static ClassDocumentation<SubtractionDipole> documentation
("SubtractionDipole represents a dipole subtraction "
"term 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<SubtractionDipole,MEBase>
describeSubtractionDipole("Herwig::SubtractionDipole", "Herwig.so");
diff --git a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
--- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
+++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
@@ -1,277 +1,275 @@
// -*- C++ -*-
//
// VBFNLOPhasespace.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 VBFNLOPhasespace class.
//
#include "VBFNLOPhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "VBFNLO/utilities/BLHAinterface.h"
#define DEFSTR(s) CPPSTR(s)
#define CPPSTR(s) #s
using namespace Herwig;
VBFNLOPhasespace::VBFNLOPhasespace() :
lastSqrtS(0*GeV), needToReshuffle(false), VBFNLOlib_(DEFSTR(VBFNLOLIB))
{}
void VBFNLOPhasespace::loadVBFNLO() {
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") ) {
string error1 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") ) {
string error2 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.so") ) {
string error3 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.dylib") ) {
string error4 = DynamicLoader::lastErrorMessage;
throw Exception() << "VBFNLOPhasespace: failed to load libVBFNLO.so/dylib\n"
<< "Error messages are:\n\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.so:\n"
<< error1 << "\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.dylib:\n"
<< error2 << "\n"
<< "* libVBFNLO.so:\n"
<< error3 << "\n"
<< "* libVBFNLO.dylib:\n"
<< error4 << "\n"
<< Exception::runerror;
}
}
}
}
}
VBFNLOPhasespace::~VBFNLOPhasespace() {}
IBPtr VBFNLOPhasespace::clone() const {
return new_ptr(*this);
}
IBPtr VBFNLOPhasespace::fullclone() const {
return new_ptr(*this);
}
void VBFNLOPhasespace::setXComb(tStdXCombPtr xco) {
MatchboxPhasespace::setXComb(xco);
// test for resuffling
needToReshuffle = false;
if ( xco ) {
for ( cPDVector::const_iterator d = mePartonData().begin();
d != mePartonData().end(); ++d ) {
// Higgs is massive -> does not need reshuffling
if ( ( (**d).id() != ParticleID::h0 ) && ( (**d).hardProcessMass() != ZERO ) ) {
needToReshuffle = true;
break;
}
}
}
// set CMS energy
int pStatus = 0;
double zero = 0.0;
double value = sqrt(lastXCombPtr()->lastS())/GeV;
if (value && (value != lastSqrtS/GeV)) {
lastSqrtS = value*GeV;
string name = "sqrtS";
OLP_SetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
if ( !pStatus )
throw Exception() << "VBFNLOPhasespace::setXComb(): VBFNLO failed to set parameter '"
<< name << "' to " << value << "\n"
<< Exception::runerror;
}
}
double VBFNLOPhasespace::generateTwoToNKinematics(const double* random,
vector<Lorentz5Momentum>& momenta) {
double weight;
int id =
olpId()[ProcessType::oneLoopInterference] ?
olpId()[ProcessType::oneLoopInterference] :
olpId()[ProcessType::treeME2];
double* p = new double[4*momenta.size()];
OLP_PhaseSpacePoint(&id, const_cast<double*>(random), const_cast<double*>(random+1), p, &weight);
if (weight < 0) {
throw Exception() << "VBFNLOPhasespace::generateTwoToNKinematics(): Negative weight in VBFNLOPhaseSpace\n"
<< Exception::runerror;
}
if (weight == 0) {
delete[] p;
return 0;
}
for ( size_t i = 0; i < momenta.size(); ++i ) {
momenta[i].setT(p[4*i] *GeV);
momenta[i].setX(p[4*i+1]*GeV);
momenta[i].setY(p[4*i+2]*GeV);
momenta[i].setZ(p[4*i+3]*GeV);
momenta[i].rescaleMass();
}
delete[] p;
Energy beamenergy = sqrt(lastXCombPtr()->lastS())/2.;
double x1 = momenta[0].e()/beamenergy;
double x2 = momenta[1].e()/beamenergy;
Energy2 thisSHat = (momenta[0] + momenta[1]).m2();
// reshuffle so that particles have correct mass
if ( needToReshuffle ) {
// boost final-state into partonic CMS
Boost toCMS = (momenta[0]+momenta[1]).findBoostToCM();
for ( size_t i = 2; i < momenta.size(); ++i ) {
momenta[i].boost(toCMS);
}
// copied from MatchboxRambo phasespace
double xi;
ReshuffleEquation solve(sqrt(thisSHat),mePartonData().begin()+2,mePartonData().end(),
momenta.begin()+2,momenta.end());
GSLBisection solver(1e-10,1e-8,10000);
try {
xi = solver.value(solve,0.0,1.1);
} catch (GSLBisection::GSLerror) {
return 0.;
} catch (GSLBisection::IntervalError) {
return 0.;
}
weight *= pow(xi,3.*(momenta.size()-3.));
Energy num = ZERO;
Energy den = ZERO;
cPDVector::const_iterator d = mePartonData().begin()+2;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin()+2;
k != momenta.end(); ++k, ++d ) {
num += (*k).vect().mag2()/(*k).t();
Energy q = (*k).t();
(*k).setT(sqrt(sqr((**d).hardProcessMass())+xi*xi*sqr((*k).t())));
(*k).setVect(xi*(*k).vect());
weight *= q/(*k).t();
den += (*k).vect().mag2()/(*k).t();
(*k).setMass((**d).hardProcessMass());
}
// unboost
for ( size_t i = 2; i < momenta.size(); ++i ) {
momenta[i].boost(-toCMS);
}
}
if ( !matchConstraints(momenta) )
return 0.;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat(thisSHat);
weight /= pow(thisSHat/GeV2,momenta.size()-4);
weight /= x1*x2;
- fillDiagramWeights();
-
return weight;
}
int VBFNLOPhasespace::nDimPhasespace(int nFinal) const {
return 3*nFinal;
//get this from within VBFNLO
int pStatus = 0;
double value, zero;
string name = "PSdimension";
OLP_GetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
if ( pStatus != 1) {
throw Exception() << "VBFNLOPhasespace::nDimPhasespace(): Cannot get phasespace dimension in VBFNLOPhaseSpace\n"
<< "error code: " << pStatus << "\n"
<< Exception::runerror;
}
// one additional number (first) needed for channel selection
// one additional number (last) needed for global phi integration
return value+2;
}
Energy VBFNLOPhasespace::ReshuffleEquation::operator() (double xi) const {
cPDVector::const_iterator d = dataBegin;
vector<Lorentz5Momentum>::const_iterator p = momentaBegin;
Energy res = -w;
for ( ; d != dataEnd; ++d, ++p ) {
res += sqrt(sqr((**d).hardProcessMass()) +
xi*xi*sqr(p->t()));
}
return res;
}
void VBFNLOPhasespace::doinit() {
loadVBFNLO();
MatchboxPhasespace::doinit();
}
void VBFNLOPhasespace::doinitrun() {
loadVBFNLO();
MatchboxPhasespace::doinitrun();
}
void VBFNLOPhasespace::persistentOutput(PersistentOStream & os) const {
os << needToReshuffle << theLastXComb;
}
void VBFNLOPhasespace::persistentInput(PersistentIStream & is, int) {
is >> needToReshuffle >> theLastXComb;
}
// *** 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<VBFNLOPhasespace,MatchboxPhasespace>
describeHerwigVBFNLOPhasespace("Herwig::VBFNLOPhasespace", "HwMatchboxVBFNLO.so");
void VBFNLOPhasespace::Init() {
static ClassDocumentation<VBFNLOPhasespace> documentation
("VBFNLOPhasespace is an interface to the internal phasespace generator "
"of VBFNLO. It uses the information passed via the BLHA interface to "
"obtain information on the required channels.");
}
diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertibleLabframePhasespace.cc b/MatrixElement/Matchbox/Phasespace/FlatInvertibleLabframePhasespace.cc
--- a/MatrixElement/Matchbox/Phasespace/FlatInvertibleLabframePhasespace.cc
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertibleLabframePhasespace.cc
@@ -1,173 +1,171 @@
// -*- C++ -*-
//
// FlatInvertiblePhasespace.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 FlatInvertibleLabframePhasespace class.
//
#include "FlatInvertibleLabframePhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FlatInvertibleLabframePhasespace::FlatInvertibleLabframePhasespace()
: theLogSHat(false) {}
FlatInvertibleLabframePhasespace::~FlatInvertibleLabframePhasespace() {}
IBPtr FlatInvertibleLabframePhasespace::clone() const {
return new_ptr(*this);
}
IBPtr FlatInvertibleLabframePhasespace::fullclone() const {
return new_ptr(*this);
}
double FlatInvertibleLabframePhasespace::invertTwoToNKinematics(const vector<Lorentz5Momentum>& momenta,
double* r) const {
double weight = 1.;
Energy finalstatemass = 0*GeV;
for ( vector<Lorentz5Momentum>::const_iterator p =
momenta.begin()+2; p != momenta.end(); ++p )
finalstatemass += p->mass();
Lorentz5Momentum pinitial = momenta[0]+momenta[1];
Energy2 sh = pinitial.m2();
double tau = sh/lastS();
Energy2 shmax = lastCuts().sHatMax();
Energy2 shmin = max(lastCuts().sHatMin(),sqr(finalstatemass));
if (theLogSHat) {
r[0] = log(sh/shmin)/log(shmax/shmin);
weight *= tau*log(shmax/shmin);
} else {
r[0] = (sh-shmin)/(shmax-shmin);
weight *= (shmax-shmin)/lastS();
}
double ltau = log(tau);
r[1] = 0.5 - pinitial.rapidity()/ltau;
weight *= -ltau;
vector<Lorentz5Momentum> Pcms = momenta;
Boost toCMS = pinitial.findBoostToCM();
for ( vector<Lorentz5Momentum>::iterator pit =
Pcms.begin(); pit != Pcms.end(); ++pit )
pit->boost(toCMS);
weight *= FlatInvertiblePhasespace::invertTwoToNKinematics(Pcms, r+2);
return weight;
}
double FlatInvertibleLabframePhasespace::generateTwoToNKinematics(const double* r,
vector<Lorentz5Momentum>& momenta) {
double weight = 1.;
Energy finalstatemass = 0*GeV;
for ( vector<Lorentz5Momentum>::const_iterator p =
momenta.begin()+2; p != momenta.end(); ++p )
finalstatemass += p->mass();
Energy beamenergy = sqrt(lastS())/2.;
Energy2 shmax = lastCuts().sHatMax();
Energy2 shmin = max(lastCuts().sHatMin(),sqr(finalstatemass));
Energy2 sh;
double tau;
if (theLogSHat) {
sh = shmin*pow(shmax/shmin, r[0]);
tau = sh/lastS();
weight *= tau*log(shmax/shmin);
} else {
sh = r[0]*(shmax-shmin)+shmin;
tau = sh/lastS();
weight *= (shmax-shmin)/lastS();
}
double ltau = log(tau);
double y = ltau*(0.5 - r[1]);
weight *= -ltau;
double x1 = sqrt(tau)*exp(y);
double x2 = sqrt(tau)*exp(-y);
momenta[0] = Lorentz5Momentum(0*GeV,0*GeV,+x1*beamenergy,x1*beamenergy);
momenta[1] = Lorentz5Momentum(0*GeV,0*GeV,-x2*beamenergy,x2*beamenergy);
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat(sh);
weight *= FlatInvertiblePhasespace::generateTwoToNKinematics(r+2, momenta);
// find boost to the relevant partonic frame note final state kinematics are
// always generated in the CMS for this phase space algorithm
Boost boostinitial = (momenta[0]+momenta[1]).findBoostToCM();
for ( vector<Lorentz5Momentum>::iterator pit =
momenta.begin()+2; pit != momenta.end(); ++pit )
pit->boost(-boostinitial);
- fillDiagramWeights();
-
return weight;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FlatInvertibleLabframePhasespace::persistentOutput(PersistentOStream & os) const {
os << theLogSHat;
}
void FlatInvertibleLabframePhasespace::persistentInput(PersistentIStream & is, int) {
is >> theLogSHat;
}
// *** 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<FlatInvertibleLabframePhasespace,MatchboxPhasespace>
describeHerwigFlatInvertibleLabframePhasespace("Herwig::FlatInvertibleLabframePhasespace", "Herwig.so");
void FlatInvertibleLabframePhasespace::Init() {
static ClassDocumentation<FlatInvertibleLabframePhasespace> documentation
("FlatInvertibleLabframePhasespace implements flat, invertible phase space generation in the lab frame.");
static Switch<FlatInvertibleLabframePhasespace,bool> interfaceLogSHat
("LogSHat",
"Generate a flat distribution in \\f$\\log(\\hat{s})\\f$.",
&FlatInvertibleLabframePhasespace::theLogSHat, false, false, false);
static SwitchOption interfaceLogSHatOn
(interfaceLogSHat,
"True", "Generate flat in \\f$\\log(\\hat{s})\\f$", true);
static SwitchOption interfaceLogSHatOff
(interfaceLogSHat,
"False", "Generate flat in \\f$\\hat{s}\\f$", false);
interfaceLogSHat.rank(-1);
}
diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
--- a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.cc
@@ -1,307 +1,304 @@
// -*- C++ -*-
//
// FlatInvertiblePhasespace.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 FlatInvertiblePhasespace class.
//
#include "FlatInvertiblePhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
FlatInvertiblePhasespace::FlatInvertiblePhasespace() {}
FlatInvertiblePhasespace::~FlatInvertiblePhasespace() {}
IBPtr FlatInvertiblePhasespace::clone() const {
return new_ptr(*this);
}
IBPtr FlatInvertiblePhasespace::fullclone() const {
return new_ptr(*this);
}
double FlatInvertiblePhasespace::bisect(double v, double n,
double target, double maxLevel) const {
if ( v != 0.0 && v != 1.0 ) {
double level = 0;
double left = 0;
double right = 1;
double checkV = -1.;
double u = -1;
while ( level < maxLevel ) {
u = (left+right)*pow(0.5,level+1.);
checkV =
pow(u,n+1.)*(n+2.-(n+1.)*u);
if ( log10(abs(1.-checkV/v)) <= target )
break;
left *= 2.;
right *= 2.;
if ( v <= checkV ) {
right -= 1.;
++level;
}
if ( v > checkV ) {
left += 1.;
++level;
}
}
return u;
}
return v;
}
-static double flatWeights[7] = {
-
- -1.,-1.,
- 0.039788735772973833942,
- 0.00012598255637968550463,
- 1.3296564302788840628E-7,
- 7.0167897579949011130E-11,
- 2.2217170114046130768E-14
-
-};
-
double FlatInvertiblePhasespace::generateIntermediates(vector<Energy>& K,
const double* r) const {
size_t n = K.size() + 1;
+
+
for ( size_t i = 2; i <= n-1; ++i ) {
double u = bisect(r[i-2],n-1-i);
K[i-1] = sqrt(u*sqr(K[i-2]));
}
-
- return flatWeights[n];
+
+ int kap = K.size() + 1;
+
+ return flatWeights(kap);
}
double FlatInvertiblePhasespace::invertIntermediates(const vector<Energy>& K,
double* r) const {
size_t n = K.size() + 1;
for ( size_t i = 2; i <= n-1; ++i ) {
double u = sqr(K[i-1]/K[i-2]);
r[i-2] = (n+1-i)*pow(u,(double)(n-i)) - (n-i)*pow(u,(double)(n+1-i));
}
- return flatWeights[n];
+ int kap = K.size() + 1;
+ return flatWeights(kap);
}
double FlatInvertiblePhasespace::generateIntermediates(vector<Energy>& M,
const vector<Energy>& m,
const double* r) const {
size_t n = M.size() + 1;
vector<Energy> K = M;
for ( size_t i = 1; i <= n; ++i )
K[0] -= m[i-1];
double w0 = generateIntermediates(K,r);
M = K;
for ( size_t i = 1; i <= n-1; ++i ) {
for ( size_t k = i; k <= n; ++k )
M[i-1] += m[k-1];
}
double weight = 8.*w0*rho(M[n-2],m[n-1],m[n-2]);
for ( size_t i = 2; i <= n-1; ++i ) {
weight *=
(rho(M[i-2],M[i-1],m[i-2])/rho(K[i-2],K[i-1],ZERO)) * (M[i-1]/K[i-1]);
}
weight *= pow(K[0]/M[0],2.*n-4.);
-
+
return weight;
}
double FlatInvertiblePhasespace::invertIntermediates(const vector<Energy>& M,
const vector<Energy>& m,
double* r) const {
size_t n = M.size() + 1;
vector<Energy> K = M;
for ( size_t i = 1; i <= n-1; ++i ) {
for ( size_t k = i; k <= n; ++k )
K[i-1] -= m[k-1];
}
double w0 = invertIntermediates(K,r);
double weight = 8.*w0*rho(M[n-2],m[n-1],m[n-2]);
for ( size_t i = 2; i <= n-1; ++i ) {
weight *=
(rho(M[i-2],M[i-1],m[i-2])/rho(K[i-2],K[i-1],ZERO)) * (M[i-1]/K[i-1]);
}
weight *= pow(K[0]/M[0],2.*n-4.);
return weight;
}
double FlatInvertiblePhasespace::generateKinematics(vector<Lorentz5Momentum>& P,
Energy Ecm,
const double* r) const {
vector<Energy> m;
for ( vector<Lorentz5Momentum>::const_iterator p =
P.begin() + 2; p != P.end(); ++p )
m.push_back(p->mass());
size_t n = P.size() - 2;
vector<Energy> M(n-1);
M[0] = Ecm;
double weight = generateIntermediates(M,m,r);
M.push_back(m.back());
Lorentz5Momentum Q(M[0]);
Lorentz5Momentum nextQ;
for ( size_t i = 2; i <= n; ++i ) {
Energy q = 4.*M[i-2]*rho(M[i-2],M[i-1],m[i-2]);
double c = 2.*r[n-6+2*i]-1.;
double s = sqrt(1.-sqr(c));
double phi = 2.*Constants::pi*r[n-5+2*i];
double cphi = cos(phi);
double sphi = sqrt(1.-sqr(cphi));
if ( phi > Constants::pi )
sphi = -sphi;
P[i].setX(q*cphi*s);
P[i].setY(q*sphi*s);
P[i].setZ(q*c);
P[i].rescaleEnergy();
P[i].boost(Q.boostVector());
P[i].rescaleEnergy();
nextQ = Q - P[i];
nextQ.setMass(M[i-1]);
nextQ.rescaleEnergy();
Q = nextQ;
}
P.back() = Q;
return weight;
}
double FlatInvertiblePhasespace::invertKinematics(const vector<Lorentz5Momentum>& P,
Energy Ecm,
double* r) const {
vector<Energy> m;
for ( vector<Lorentz5Momentum>::const_iterator p =
P.begin() + 2; p != P.end(); ++p )
m.push_back(p->mass());
size_t n = P.size() - 2;
vector<Energy> M(n-1);
M[0] = Ecm;
vector<Lorentz5Momentum> Q(n-1);
Q[0] = Lorentz5Momentum(M[0]);
for ( size_t i = 2; i <= n-1; ++i ) {
for ( size_t k = i; k <= n; ++k )
Q[i-1] += P[k+1];
M[i-1] = Q[i-1].m();
}
double weight = invertIntermediates(M,m,r);
for ( size_t i = 2; i <= n; ++i ) {
Lorentz5Momentum p = P[i];
p.boost(-Q[i-2].boostVector());
r[n-6+2*i] = (p.cosTheta()+1.)/2.;
double phi = p.phi();
if ( phi < 0. )
phi = 2.*Constants::pi + phi;
r[n-5+2*i] = phi/(2.*Constants::pi);
}
return weight;
}
double FlatInvertiblePhasespace::generateTwoToNKinematics(const double* r,
vector<Lorentz5Momentum>& momenta) {
double weight = generateKinematics(momenta,sqrt(lastXCombPtr()->lastSHat()),r);
-
- fillDiagramWeights();
-
return weight;
}
+long double FlatInvertiblePhasespace::flatWeights(int k) const{
+ using Constants::pi;
+ if(k<2) { return -1;
+ } else return pow((pi/2),(k-1)) * pow((2*pi),(4-3*k))/factorial(k-1)/factorial(k-2);
+}
+
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FlatInvertiblePhasespace::persistentOutput(PersistentOStream &) const {}
void FlatInvertiblePhasespace::persistentInput(PersistentIStream &, int) {}
// *** 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<FlatInvertiblePhasespace,MatchboxPhasespace>
describeHerwigFlatInvertiblePhasespace("Herwig::FlatInvertiblePhasespace", "Herwig.so");
void FlatInvertiblePhasespace::Init() {
static ClassDocumentation<FlatInvertiblePhasespace> documentation
("FlatInvertiblePhasespace implements flat, invertible phase space generation.");
}
diff --git a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
--- a/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
+++ b/MatrixElement/Matchbox/Phasespace/FlatInvertiblePhasespace.h
@@ -1,194 +1,204 @@
// -*- C++ -*-
//
// FlatInvertiblePhasespace.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_FlatInvertiblePhasespace_H
#define Herwig_FlatInvertiblePhasespace_H
//
// This is the declaration of the FlatInvertiblePhasespace class.
//
#include "Herwig/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief FlatInvertiblePhasespace implements flat, invertible phase space generation.
*
*/
class FlatInvertiblePhasespace: public MatchboxPhasespace {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FlatInvertiblePhasespace();
/**
* The destructor.
*/
virtual ~FlatInvertiblePhasespace();
//@}
public:
/**
* Generate a phase space point and return its weight.
*/
virtual double generateTwoToNKinematics(const double*,
vector<Lorentz5Momentum>& momenta);
/**
* Return the number of random numbers required to produce a given
* multiplicity final state.
*/
virtual int nDimPhasespace(int nFinal) const {
if ( nFinal == 1 )
return 1;
return 3*nFinal - 4;
}
public:
/**
* Return true, if this phase space generator is invertible
*/
virtual bool isInvertible() const { return true; }
/**
* Invert the given phase space point to the random numbers which
* would have generated it.
*/
virtual double invertTwoToNKinematics(const vector<Lorentz5Momentum>& momenta,
double* r) const {
return invertKinematics(momenta,(momenta[0]+momenta[1]).m(),r);
}
private:
/**
* Solve v = (n+2) * u^(n+1) - (n+1) * u^(n+2) for u
*/
double bisect(double v, double n,
double target = -16., double maxLevel = 80.) const;
/**
* Return rho
*/
double rho(Energy M, Energy N, Energy m) const {
return sqrt((sqr(M)-sqr(N+m))*(sqr(M)-sqr(N-m)))/(8.*sqr(M));
}
/**
* Generate intermediate masses for a massless final state
*/
double generateIntermediates(vector<Energy>& K,
const double* r) const;
/**
* Invert intermediate masses for a massless final state
*/
double invertIntermediates(const vector<Energy>& K,
double* r) const;
/**
* Generate intermediate masses for a massive final state
*/
double generateIntermediates(vector<Energy>& M,
const vector<Energy>& m,
const double* r) const;
/**
* Invert intermediate masses for a massive final state
*/
double invertIntermediates(const vector<Energy>& M,
const vector<Energy>& m,
double* r) const;
/**
* Generate momenta in the CMS
*/
double generateKinematics(vector<Lorentz5Momentum>& P,
Energy Ecm,
const double* r) const;
/**
* Invert momenta in the CMS
*/
double invertKinematics(const vector<Lorentz5Momentum>& P,
Energy Ecm,
double* r) const;
+ /**
+ * Return the appropriate phase space weight,
+ * Eq. 11 in 1308.2922
+ * with the factor (2 pi)^4/(2 pi)^(3n) included
+ * and the SHat of the process divided out to have everything expressed in the units of the ThePEG conventions, i.e.
+ * without the Q^2 factor
+ */
+
+ long double flatWeights(int n) 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();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FlatInvertiblePhasespace & operator=(const FlatInvertiblePhasespace &);
};
}
#endif /* Herwig_FlatInvertiblePhasespace_H */
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc
--- a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc
@@ -1,562 +1,565 @@
// -*- C++ -*-
//
// MatchboxPhasespace.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 MatchboxPhasespace class.
//
#include "MatchboxPhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Command.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 "Herwig/MatrixElement/Matchbox/Utility/ProcessData.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
using namespace Herwig;
MatchboxPhasespace::MatchboxPhasespace()
: singularCutoff(10*GeV), theUseMassGenerators(false),
theLoopParticleIdMin(200001), theLoopParticleIdMax(200100) {}
MatchboxPhasespace::~MatchboxPhasespace() {}
void MatchboxPhasespace::cloneDependencies(const std::string&) {}
Ptr<MatchboxFactory>::tcptr MatchboxPhasespace::factory() const {
return lastMatchboxXComb()->factory();
}
Ptr<ProcessData>::tptr MatchboxPhasespace::processData() const {
return factory()->processData();
}
double MatchboxPhasespace::generateKinematics(const double* r,
vector<Lorentz5Momentum>& momenta) {
+ diagramWeights().clear();
+
cPDVector::const_iterator pd = mePartonData().begin() + 2;
vector<Lorentz5Momentum>::iterator p = momenta.begin() + 2;
double massJacobian = 1.;
Energy summ = ZERO;
if ( useMassGenerators() ) {
Energy gmass = ZERO;
tGenericMassGeneratorPtr mgen;
Energy maxMass =
(!haveX1X2() && momenta.size() > 3) ?
sqrt(lastSHat()) : sqrt(lastS());
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
mgen = processData()->massGenerator(*pd);
if ( mgen && !isInvertible() ) {
Energy massMax = min((**pd).massMax(),maxMass);
Energy massMin = (**pd).massMin();
if ( massMin > massMax )
return 0.0;
gmass = mgen->mass(massJacobian,**pd,massMin,massMax,r[0]);
++r;
} else if ( (**pd).hardProcessWidth() != ZERO ) {
Energy massMax = min((**pd).massMax(),maxMass);
Energy massMin = (**pd).massMin();
if ( massMin > massMax )
return 0.0;
// use a standard Breit Wigner here which we can invert
// see invertKinematics as well
double bwILow =
atan((sqr(massMin)-sqr((**pd).hardProcessMass()))/((**pd).hardProcessMass() * (**pd).hardProcessWidth()));
double bwIUp =
atan((sqr(massMax)-sqr((**pd).hardProcessMass()))/((**pd).hardProcessMass() * (**pd).hardProcessWidth()));
gmass = sqrt(sqr((**pd).hardProcessMass()) +
(**pd).hardProcessMass()*(**pd).hardProcessWidth()*tan(bwILow+r[0]*(bwIUp-bwILow)));
++r;
} else {
gmass = (**pd).hardProcessMass();
}
maxMass -= gmass;
p->setMass(gmass);
summ += gmass;
}
} else {
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
summ += (**pd).hardProcessMass();
p->setMass((**pd).hardProcessMass());
}
}
if ( momenta.size() > 3 && !haveX1X2() ) {
if ( summ > (momenta[0]+momenta[1]).m() )
return 0.0;
}
double weight = momenta.size() > 3 ?
generateTwoToNKinematics(r,momenta) :
generateTwoToOneKinematics(r,momenta);
+ fillDiagramWeights();
+
return weight*massJacobian;
}
double MatchboxPhasespace::generateTwoToOneKinematics(const double* r,
vector<Lorentz5Momentum>& momenta) {
double tau = momenta[2].mass2()/lastXCombPtr()->lastS();
double ltau = log(tau)/2.;
//old: y = ltau - 2.*r[0]*ltau; x1 = sqrt(tau)*exp(y); x2 = sqrt(tau)*exp(-y);
double x1=pow(tau,1.-r[0]);
double x2=pow(tau,r[0]);
// Due to the proton mass and P1.e() + P2.e() == lastS() we multiply here
// with the correction factor abs(P1.e()/P1.z()) to produce incoming
// p1/2 = (e1/2,0,0,+/- e1/2)
Lorentz5Momentum P1 = lastXCombPtr()->lastParticles().first->momentum();
ThreeVector<Energy> p1 = x1 * (P1.vect()) * abs(P1.e()/P1.z());
Lorentz5Momentum P2 = lastXCombPtr()->lastParticles().second->momentum();
ThreeVector<Energy> p2 = x2 * (P2.vect()) * abs(P2.e()/P2.z());
ThreeVector<Energy> q = p1 + p2;
momenta[0] = Lorentz5Momentum(momenta[0].mass(),p1);
momenta[1] = Lorentz5Momentum(momenta[1].mass(),p2);
momenta[2] = Lorentz5Momentum(momenta[2].mass(),q);
// check for energy conservation:
if ((momenta[0]+momenta[1]-momenta[2]).e()>pow(10,-9)*GeV)
generator()->log()
<< "Warning: Momentum conservation in generateTwoToOneKinematics not precise.\n"
<< flush;
lastXCombPtr()->lastX1X2({x1,x2});
lastXCombPtr()->lastSHat((momenta[0]+momenta[1]).m2());
- fillDiagramWeights();
-
return -4.*Constants::pi*ltau;
}
double MatchboxPhasespace::invertKinematics(const vector<Lorentz5Momentum>& momenta,
double* r) const {
if ( useMassGenerators() ) {
Energy gmass = ZERO;
Energy maxMass =
(!haveX1X2() && momenta.size() > 3) ?
sqrt((momenta[0]+momenta[1]).m2()) : sqrt(lastS());
cPDVector::const_iterator pd = mePartonData().begin() + 2;
vector<Lorentz5Momentum>::const_iterator p = momenta.begin() + 2;
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
if ( (**pd).hardProcessWidth() != ZERO ) {
Energy massMax = min((**pd).massMax(),maxMass);
Energy massMin = (**pd).massMin();
if ( massMin > massMax )
return 0.0;
double bwILow =
atan((sqr(massMin)-sqr((**pd).hardProcessMass()))/((**pd).hardProcessMass() * (**pd).hardProcessWidth()));
double bwIUp =
atan((sqr(massMax)-sqr((**pd).hardProcessMass()))/((**pd).hardProcessMass() * (**pd).hardProcessWidth()));
gmass = p->mass();
double bw =
atan((sqr(gmass)-sqr((**pd).hardProcessMass()))/((**pd).hardProcessMass() * (**pd).hardProcessWidth()));
r[0] = (bw-bwILow)/(bwIUp-bwILow);
++r;
} else {
gmass = (**pd).hardProcessMass();
}
maxMass -= gmass;
}
}
return momenta.size() > 3 ?
invertTwoToNKinematics(momenta,r) :
invertTwoToOneKinematics(momenta,r);
}
double MatchboxPhasespace::invertTwoToOneKinematics(const vector<Lorentz5Momentum>& momenta,
double* r) const {
double tau = momenta[2].mass2()/lastXCombPtr()->lastS();
double ltau = log(tau)/2.;
r[0] = (ltau - (momenta[0]+momenta[1]).rapidity())/(2.*ltau);
return -4.*Constants::pi*ltau;
}
void MatchboxPhasespace::setCoupling(long a, long b, long c,
double coupling, bool includeCrossings) {
cPDPtr A = getParticleData(a);
cPDPtr B = getParticleData(b);
cPDPtr C = getParticleData(c);
if ( !A || !B || !C ) {
generator()->log() << "Warning: could not determine particle data for ids "
<< a << " " << b << " " << c << " when setting coupling in MatchboxPhasespace.\n"
<< flush;
return;
}
if ( !includeCrossings ) {
theCouplings->couplings()[LTriple(a,b,c)] = coupling;
return;
}
if ( A->CC() ) {
theCouplings->couplings()[LTriple(-a,b,c)] = coupling;
theCouplings->couplings()[LTriple(-a,c,b)] = coupling;
} else {
theCouplings->couplings()[LTriple(a,b,c)] = coupling;
theCouplings->couplings()[LTriple(a,c,b)] = coupling;
}
if ( B->CC() ) {
theCouplings->couplings()[LTriple(-b,a,c)] = coupling;
theCouplings->couplings()[LTriple(-b,c,a)] = coupling;
} else {
theCouplings->couplings()[LTriple(b,a,c)] = coupling;
theCouplings->couplings()[LTriple(b,c,a)] = coupling;
}
if ( C->CC() ) {
theCouplings->couplings()[LTriple(-c,a,b)] = coupling;
theCouplings->couplings()[LTriple(-c,b,a)] = coupling;
} else {
theCouplings->couplings()[LTriple(c,a,b)] = coupling;
theCouplings->couplings()[LTriple(c,b,a)] = coupling;
}
}
string MatchboxPhasespace::doSetCoupling(string in) {
istringstream is(in);
long a,b,c; double coupling;
is >> a >> b >> c >> coupling;
if ( !is )
return "MatchboxPhasespace: error in setting coupling.";
setCoupling(a,b,c,coupling,true);
return "";
}
string MatchboxPhasespace::doSetPhysicalCoupling(string in) {
istringstream is(in);
long a,b,c; double coupling;
is >> a >> b >> c >> coupling;
if ( !is )
return "MatchboxPhasespace: error in setting coupling.";
setCoupling(a,b,c,coupling,false);
return "";
}
pair<double,Lorentz5Momentum>
MatchboxPhasespace::timeLikeWeight(const Tree2toNDiagram& diag,
int branch, double flatCut) const {
pair<int,int> children = diag.children(branch);
if ( children.first == -1 ) {
return make_pair(1.,meMomenta()[diag.externalId(branch)]);
}
pair<double,Lorentz5Momentum> res
= timeLikeWeight(diag,children.first,flatCut);
pair<double,Lorentz5Momentum> other
= timeLikeWeight(diag,children.second,flatCut);
res.first *= other.first;
res.second += other.second;
LTriple vertexKey(diag.allPartons()[branch]->id(),
diag.allPartons()[children.first]->id(),
diag.allPartons()[children.second]->id());
map<LTriple,double>::const_iterator cit = theCouplings->couplings().find(vertexKey);
if ( cit != theCouplings->couplings().end() ){
res.first *= cit->second;
}
Energy2 mass2 = sqr(diag.allPartons()[branch]->hardProcessMass());
Energy2 width2 = sqr(diag.allPartons()[branch]->hardProcessWidth());
if ( abs(diag.allPartons()[branch]->id()) >= theLoopParticleIdMin
&& abs(diag.allPartons()[branch]->id()) <= theLoopParticleIdMax ) { // "loop particle"
if ( abs((res.second.m2()-mass2)/lastSHat()) > flatCut ) {
res.first /=
abs((res.second.m2()-mass2)/GeV2);
res.first *=
log(abs((res.second.m2()-mass2)/GeV2)); // normal. of the argument in the log?
}
} else {
if ( width2 == ZERO ) {
if ( abs((res.second.m2()-mass2)/lastSHat()) > flatCut )
res.first /=
abs((res.second.m2()-mass2)/GeV2);
} else {
res.first /=
(sqr((res.second.m2()-mass2)/GeV2) +
mass2*width2/sqr(GeV2))/(abs(res.second.m2()/GeV2));
}
}
return res;
}
double MatchboxPhasespace::spaceLikeWeight(const Tree2toNDiagram& diag,
const Lorentz5Momentum& incoming,
int branch, double flatCut) const {
if ( branch == -1 )
return 1.;
pair<int,int> children = diag.children(branch);
pair<double,Lorentz5Momentum> res =
timeLikeWeight(diag,children.second,flatCut);
LTriple vertexKey(diag.allPartons()[branch]->id(),
diag.allPartons()[children.first]->id(),
diag.allPartons()[children.second]->id());
if ( children.first == diag.nSpace() - 1 ) {
if ( diag.allPartons()[children.first]->CC() )
vertexKey = LTriple(diag.allPartons()[branch]->id(),
diag.allPartons()[children.second]->id(),
diag.allPartons()[children.first]->CC()->id());
else
vertexKey = LTriple(diag.allPartons()[branch]->id(),
diag.allPartons()[children.second]->id(),
diag.allPartons()[children.first]->id());
}
map<LTriple,double>::const_iterator cit = theCouplings->couplings().find(vertexKey);
if ( cit != theCouplings->couplings().end() ){
res.first *= cit->second;
}
if ( children.first == diag.nSpace() - 1 ) {
return res.first;
}
res.second = incoming - res.second;
Energy2 mass2 = sqr(diag.allPartons()[children.first]->hardProcessMass());
Energy2 width2 = sqr(diag.allPartons()[children.first]->hardProcessWidth());
if ( abs(diag.allPartons()[children.first]->id()) >= theLoopParticleIdMin
&& (diag.allPartons()[children.first]->id()) <= theLoopParticleIdMax ) { // "loop particle"
if ( abs((res.second.m2()-mass2)/lastSHat()) > flatCut ) {
res.first /=
abs((res.second.m2()-mass2)/GeV2);
res.first *=
log(abs((res.second.m2()-mass2)/GeV2)); // normal. of the argument in the log?
}
} else {
if ( width2 == ZERO ) {
if ( abs((res.second.m2()-mass2)/lastSHat()) > flatCut )
res.first /=
abs((res.second.m2()-mass2)/GeV2);
} else {
res.first /=
(sqr((res.second.m2()-mass2)/GeV2) +
mass2*width2/sqr(GeV2))/(abs(res.second.m2()/GeV2));
}
}
return
res.first * spaceLikeWeight(diag,res.second,children.first,flatCut);
}
void MatchboxPhasespace::fillDiagramWeights(double flatCut) {
- diagramWeights().clear();
+ if ( !diagramWeights().empty() )
+ return;
for ( auto & d : lastXComb().diagrams() ) {
diagramWeights()[d->id()] =
spaceLikeWeight(dynamic_cast<const Tree2toNDiagram&>(*d),meMomenta()[0],0,flatCut);
}
}
Selector<MEBase::DiagramIndex>
MatchboxPhasespace::selectDiagrams(const MEBase::DiagramVector& diags) const {
Selector<MEBase::DiagramIndex> ret;
for ( MEBase::DiagramIndex d = 0; d < diags.size(); ++d ) {
ret.insert(diagramWeight(dynamic_cast<const Tree2toNDiagram&>(*diags[d])),d);
}
return ret;
}
bool MatchboxPhasespace::matchConstraints(const vector<Lorentz5Momentum>& momenta) {
if ( singularLimits().empty() )
return true;
lastSingularLimit() = singularLimits().begin();
for ( ; lastSingularLimit() != singularLimits().end(); ++lastSingularLimit() ) {
if ( lastSingularLimit()->first == lastSingularLimit()->second &&
momenta[lastSingularLimit()->first].t() < singularCutoff )
break;
if ( lastSingularLimit()->first != lastSingularLimit()->second &&
sqrt(momenta[lastSingularLimit()->first]*
momenta[lastSingularLimit()->second]) < singularCutoff ) {
bool match = true;
for ( set<pair<size_t,size_t> >::const_iterator other =
singularLimits().begin(); other != singularLimits().end(); ++other ) {
if ( other == lastSingularLimit() )
continue;
if ( other->first == other->second &&
momenta[other->first].t() < singularCutoff ) {
match = false;
break;
}
if ( other->first != other->second &&
sqrt(momenta[other->first]*
momenta[other->second]) < singularCutoff ) {
match = false;
break;
}
}
if ( match )
break;
}
}
return lastSingularLimit() != singularLimits().end();
}
int MatchboxPhasespace::nDim(const cPDVector& data) const {
int ndimps = nDimPhasespace(data.size()-2);
if ( useMassGenerators() ) {
for ( cPDVector::const_iterator pd = data.begin();
pd != data.end(); ++pd ) {
if ( (**pd).massGenerator() ||
(**pd).hardProcessWidth() != ZERO ) {
++ndimps;
}
}
}
return ndimps;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void MatchboxPhasespace::persistentOutput(PersistentOStream & os) const {
os << theLastXComb
<< ounit(singularCutoff,GeV) << theUseMassGenerators
<< theLoopParticleIdMin << theLoopParticleIdMax
<< theCouplings;
}
void MatchboxPhasespace::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb
>> iunit(singularCutoff,GeV) >> theUseMassGenerators
>> theLoopParticleIdMin >> theLoopParticleIdMax
>> theCouplings;
lastMatchboxXComb(theLastXComb);
}
// *** 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<MatchboxPhasespace,HandlerBase>
describeMatchboxPhasespace("Herwig::MatchboxPhasespace", "Herwig.so");
void MatchboxPhasespace::Init() {
static ClassDocumentation<MatchboxPhasespace> documentation
("MatchboxPhasespace defines an abstract interface to a phase "
"space generator.");
static Parameter<MatchboxPhasespace,Energy> interfaceSingularCutoff
("SingularCutoff",
"[debug] Cutoff below which a region is considered singular.",
&MatchboxPhasespace::singularCutoff, GeV, 10.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
interfaceSingularCutoff.rank(-1);
/*
static Switch<MatchboxPhasespace,bool> interfaceUseMassGenerators
("UseMassGenerators",
"Use mass generators instead of fixed masses.",
&MatchboxPhasespace::theUseMassGenerators, false, false, false);
static SwitchOption interfaceUseMassGeneratorsYes
(interfaceUseMassGenerators,
"Yes",
"Use mass generators.",
true);
static SwitchOption interfaceUseMassGeneratorsNo
(interfaceUseMassGenerators,
"No",
"Do not use mass generators.",
false);
*/
static Command<MatchboxPhasespace> interfaceSetCoupling
("SetCoupling",
"",
&MatchboxPhasespace::doSetCoupling, false);
static Command<MatchboxPhasespace> interfaceSetPhysicalCoupling
("SetPhysicalCoupling",
"",
&MatchboxPhasespace::doSetPhysicalCoupling, false);
static Parameter<MatchboxPhasespace,int> interfaceLoopParticleIdMin
("LoopParticleIdMin",
"First id in a range of id's meant to denote fictitious "
"'ghost' particles to be used by the diagram generator "
"in loop induced processes.",
&MatchboxPhasespace::theLoopParticleIdMin, 200001, 0, 0,
false, false, Interface::lowerlim);
interfaceLoopParticleIdMin.rank(-1);
static Parameter<MatchboxPhasespace,int> interfaceLoopParticleIdMax
("LoopParticleIdMax",
"Last id in a range of id's meant to denote fictitious "
"'ghost' particles to be used by the diagram generator "
"in loop induced processes.",
&MatchboxPhasespace::theLoopParticleIdMax, 200100, 0, 0,
false, false, Interface::lowerlim);
interfaceLoopParticleIdMax.rank(-1);
static Reference<MatchboxPhasespace,PhasespaceCouplings> interfaceCouplingData
("CouplingData",
"Set the storage for the couplings.",
&MatchboxPhasespace::theCouplings, false, false, true, false, false);
interfaceCouplingData.rank(-1);
}
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
--- a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
@@ -1,359 +1,366 @@
// -*- C++ -*-
//
// MatchboxPhasespace.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_MatchboxPhasespace_H
#define HERWIG_MatchboxPhasespace_H
//
// This is the declaration of the MatchboxPhasespace class.
//
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig/MatrixElement/Matchbox/Utility/ProcessData.fh"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
#include "Herwig/MatrixElement/Matchbox/Phasespace/PhasespaceCouplings.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Wrap around a vector of random numbers to behave as a stream
* of those.
*/
struct StreamingRnd {
/**
* The random numbers
*/
const double* numbers;
/**
* The number of random numbers available.
*/
size_t nRnd;
/**
* Default constructor.
*/
StreamingRnd()
: numbers(0), nRnd(0) {}
/**
* Construct from random numbers.
*/
explicit StreamingRnd(const double* newNumbers,
size_t n)
: numbers(newNumbers), nRnd(n) {}
/**
* Return next random number
*/
inline double operator()() {
assert(numbers && nRnd > 0);
const double ret = numbers[0];
++numbers; --nRnd;
return ret;
}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxPhasespace defines an abstract interface to a phase
* space generator.
*
*/
class MatchboxPhasespace:
public HandlerBase,
public LastXCombInfo<StandardXComb>,
public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxPhasespace();
/**
* The destructor.
*/
virtual ~MatchboxPhasespace();
//@}
public:
/**
* Set the XComb object steering the Born matrix
* element this class represents virtual corrections to.
*/
virtual void setXComb(tStdXCombPtr xc) {
theLastXComb = xc;
lastMatchboxXComb(xc);
}
/**
* Return the factory object
*/
Ptr<MatchboxFactory>::tcptr factory() const;
/**
* Return the process data object
*/
Ptr<ProcessData>::tptr processData() const;
/**
* Generate a phase space point and return its weight.
*/
virtual double generateKinematics(const double* r,
vector<Lorentz5Momentum>& momenta);
/**
* Generate a phase space point and return its weight.
*/
virtual double generateTwoToNKinematics(const double*,
vector<Lorentz5Momentum>& momenta) = 0;
/**
* Generate a 2 -> 1 phase space point and return its weight.
*/
virtual double generateTwoToOneKinematics(const double*,
vector<Lorentz5Momentum>& momenta);
/**
* Return the number of random numbers required to produce a given
* multiplicity final state.
*/
virtual int nDim(const cPDVector&) const;
/**
* Return the number of random numbers required to produce a given
* multiplicity final state.
*/
virtual int nDimPhasespace(int nFinal) const = 0;
/**
* Return true, if this phasespace generator will generate incoming
* partons itself.
*/
virtual bool haveX1X2() const { return false; }
/**
* Return true, if this phase space generator expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS() const { return true; }
/**
* True, if mass generators should be used instead of fixed masses
*/
bool useMassGenerators() const { return theUseMassGenerators; }
/**
* Fill a diagram selector for the last phase space point.
*/
virtual Selector<MEBase::DiagramIndex> selectDiagrams(const MEBase::DiagramVector&) const;
/**
* Return the momentum and weight appropriate to the given timelike
* branch of the diagram.
*/
pair<double,Lorentz5Momentum> timeLikeWeight(const Tree2toNDiagram& diag,
int branch, double flatCut) const;
/**
* Return the weight appropriate to the given spacelike branch of
* the diagram.
*/
double spaceLikeWeight(const Tree2toNDiagram& diag,
const Lorentz5Momentum& incoming,
int branch, double flatCut) const;
/**
* Return the weight appropriate to the given diagram.
*/
double diagramWeight(const Tree2toNDiagram& diag) const {
assert( !diagramWeights().empty() );
return diagramWeights().find(diag.id())->second;
}
/**
* Fill the diagram weights.
*/
void fillDiagramWeights(double flatCut = 0.0);
/**
+ * Clear the diagram weights.
+ */
+ void clearDiagramWeights() {
+ diagramWeights().clear();
+ }
+
+ /**
* Clone this phase space generator.
*/
Ptr<MatchboxPhasespace>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxPhasespace>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
virtual void cloneDependencies(const std::string& prefix = "");
public:
/**
* Return true, if this phase space generator is invertible
*/
virtual bool isInvertible() const { return false; }
/**
* Invert the given phase space point to the random numbers which
* would have generated it.
*/
virtual double invertKinematics(const vector<Lorentz5Momentum>& momenta,
double* r) const;
/**
* Invert the given phase space point to the random numbers which
* would have generated it.
*/
virtual double invertTwoToNKinematics(const vector<Lorentz5Momentum>&,
double*) const {
return 0.;
}
/**
* Invert the given 2 -> 1 phase space point to the random numbers which
* would have generated it.
*/
virtual double invertTwoToOneKinematics(const vector<Lorentz5Momentum>&, double*) const;
public:
/**
* Limit phasespace generation to a given collinear or soft limit.
*/
void singularLimit(size_t i, size_t j) {
if ( i > j )
swap(i,j);
singularLimits().insert(make_pair(i,j));
}
/**
* Return the last matched singular limit.
*/
const pair<size_t,size_t>& lastSingularIndices() const {
assert(lastSingularLimit() != singularLimits().end());
return *lastSingularLimit();
}
/**
* Return true, if constraints on phasespace generation have been met.
*/
bool matchConstraints(const vector<Lorentz5Momentum>& momenta);
protected:
/**
* Set a coupling for the given vertex; the convention is that all
* legs are outgoing, and all possible crossings will be taken care
* of. If not set, coupling weights default to one.
*/
void setCoupling(long a, long b, long c,
double coupling, bool includeCrossings = true);
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);
//@}
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();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* A cutoff below which a region is considered singular.
*/
Energy singularCutoff;
/**
* True, if mass generators should be used instead of fixed masses
*/
bool theUseMassGenerators;
/**
* Couplings to be used in diagram weighting
*/
Ptr<PhasespaceCouplings>::ptr theCouplings;
/**
* Interface function to setcoupling
*/
string doSetCoupling(string);
/**
* Interface function to setcoupling
*/
string doSetPhysicalCoupling(string);
/**
* The first id in a range of id's meant to denote fictitious
* 'ghost' particles to be used by the diagram generator
* in loop induced processes.
*/
int theLoopParticleIdMin;
/**
* The last id in a range of id's meant to denote fictitious
* 'ghost' particles to be used by the diagram generator
* in loop induced processes.
*/
int theLoopParticleIdMax;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxPhasespace & operator=(const MatchboxPhasespace &);
};
}
#endif /* HERWIG_MatchboxPhasespace_H */
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc b/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc
--- a/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc
@@ -1,264 +1,263 @@
// -*- C++ -*-
//
// MatchboxRambo.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 MatchboxRambo class.
//
#include "MatchboxRambo.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
MatchboxRambo::MatchboxRambo()
: needToReshuffle(false), theMakeReferenceSample(false),
referenceSample(0) {}
MatchboxRambo::~MatchboxRambo() {}
IBPtr MatchboxRambo::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxRambo::fullclone() const {
return new_ptr(*this);
}
static double weights[7] = {
-1.,-1.,
0.039788735772973833942,
0.00012598255637968550463,
1.3296564302788840628E-7,
7.0167897579949011130E-11,
2.2217170114046130768E-14
};
void MatchboxRambo::setXComb(tStdXCombPtr xc) {
MatchboxPhasespace::setXComb(xc);
needToReshuffle = false;
if ( xc ) {
for ( cPDVector::const_iterator d = mePartonData().begin();
d != mePartonData().end(); ++d ) {
if ( (**d).hardProcessMass() != ZERO ) {
needToReshuffle = true;
break;
}
}
}
}
void MatchboxRambo::dumpReference(const vector<Lorentz5Momentum>& momenta, double weight) const {
*referenceSample << lastX1() << " " << lastX2() << " ";
Boost toLab = (lastPartons().first->momentum() +
lastPartons().second->momentum()).boostVector();
for ( vector<Lorentz5Momentum>::const_iterator p = momenta.begin();
p != momenta.end(); ++p ) {
Lorentz5Momentum pl = *p;
if ( toLab.mag2() > Constants::epsilon )
pl.boost(toLab);
*referenceSample
<< (pl.x()/GeV) << " "
<< (pl.y()/GeV) << " "
<< (pl.z()/GeV) << " "
<< (pl.t()/GeV) << " "
<< (pl.mass()/GeV) << " ";
}
double ymax = lastCuts().yHatMax();
double ymin = lastCuts().yHatMin();
double km = log(lastCuts().sHatMax()/lastCuts().sHatMin());
ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/lastSHat())));
ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/lastSHat())));
*referenceSample << weight*km*(ymax-ymin)/(lastX1()*lastX2()) << "\n" << flush;
}
double MatchboxRambo::generateTwoToNKinematics(const double* r,
vector<Lorentz5Momentum>& momenta) {
if ( theMakeReferenceSample ) {
map<cPDVector,ofstream*>::iterator ref =
referenceSamples.find(mePartonData());
if ( ref == referenceSamples.end() ) {
ostringstream refname;
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
refname << (**p).PDGName();
}
refname << ".rambo";
referenceSamples[mePartonData()] = new ofstream(refname.str().c_str(),std::ios_base::app);
ref = referenceSamples.find(mePartonData());
*(ref->second) << setprecision(26);
}
assert(ref != referenceSamples.end());
referenceSample = ref->second;
}
- size_t offset = dynamic_cast<const Tree2toNDiagram&>(*lastXComb().diagrams().front()).nSpace() > 0 ? 2 : 1;
+ size_t offset = 2;
+ if ( lastXCombPtr() )
+ offset = dynamic_cast<const Tree2toNDiagram&>(*lastXComb().diagrams().front()).nSpace() > 0 ? 2 : 1;
Energy w = sqrt(lastSHat());
size_t count = 0;
Lorentz5Momentum Q;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin() + offset;
k != momenta.end(); ++k ) {
Energy q = -w*log(r[count]*r[count+1]);
double ct = 2.*r[count+2]-1.;
double st = sqrt(1.-sqr(ct));
double phi = 2.*Constants::pi*r[count+3];
double cphi = cos(phi);
double sphi = sqrt(1.-sqr(cphi));
if ( phi > Constants::pi )
sphi = -sphi;
(*k).setMass(ZERO);
(*k).setT(q);
(*k).setX(q*cphi*st);
(*k).setY(q*sphi*st);
(*k).setZ(q*ct);
count += 4;
Q += *k;
}
Energy M = sqrt(Q.m2());
double x = w/M;
Boost beta = -(Q.vect() * (1./M));
double gamma = Q.t()/M;
double a = 1./(1.+gamma);
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin() + offset;
k != momenta.end(); ++k ) {
Energy q = (*k).t();
Energy bq = beta*(*k).vect();
(*k).setT(x*(gamma*q+bq));
(*k).setVect(x*((*k).vect()+(q+a*bq)*beta));
}
size_t n = momenta.size()-offset;
double weight = weights[n];
if ( !needToReshuffle ) {
if ( !matchConstraints(momenta) )
return 0.;
- fillDiagramWeights();
if ( theMakeReferenceSample )
dumpReference(momenta, weight);
return weight;
}
double xi;
ReshuffleEquation solve(w,mePartonData().begin()+offset,mePartonData().end(),
momenta.begin()+2,momenta.end());
GSLBisection solver(1e-10,1e-8,10000);
try {
xi = solver.value(solve,0.0,1.1);
} catch (GSLBisection::GSLerror) {
return 0.;
} catch (GSLBisection::IntervalError) {
return 0.;
}
weight *= pow(xi,3.*(n-1.));
Energy num = ZERO;
Energy den = ZERO;
cPDVector::const_iterator d = mePartonData().begin()+offset;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin()+offset;
k != momenta.end(); ++k, ++d ) {
num += (*k).vect().mag2()/(*k).t();
Energy q = (*k).t();
(*k).setT(sqrt(sqr((**d).hardProcessMass())+xi*xi*sqr((*k).t())));
(*k).setVect(xi*(*k).vect());
weight *= q/(*k).t();
den += (*k).vect().mag2()/(*k).t();
(*k).setMass((**d).hardProcessMass());
}
if ( !matchConstraints(momenta) )
return 0.;
weight *= num/den;
- fillDiagramWeights();
-
if ( theMakeReferenceSample )
dumpReference(momenta, weight);
return weight;
}
Energy MatchboxRambo::ReshuffleEquation::operator() (double xi) const {
cPDVector::const_iterator d = dataBegin;
vector<Lorentz5Momentum>::const_iterator p = momentaBegin;
Energy res = -w;
for ( ; d != dataEnd; ++d, ++p ) {
res += sqrt(sqr((**d).hardProcessMass()) +
xi*xi*sqr(p->t()));
}
return res;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void MatchboxRambo::persistentOutput(PersistentOStream & os) const {
os << needToReshuffle << theMakeReferenceSample;
}
void MatchboxRambo::persistentInput(PersistentIStream & is, int) {
is >> needToReshuffle >> theMakeReferenceSample;
}
// *** 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<MatchboxRambo,MatchboxPhasespace>
describeHerwigMatchboxRambo("Herwig::MatchboxRambo", "Herwig.so");
void MatchboxRambo::Init() {
static ClassDocumentation<MatchboxRambo> documentation
("MatchboxRambo implements RAMBO phase space generation.");
static Switch<MatchboxRambo,bool> interfaceMakeReferenceSample
("MakeReferenceSample",
"Switch on generation of a reference sample of phase space points.",
&MatchboxRambo::theMakeReferenceSample, false, false, false);
static SwitchOption interfaceMakeReferenceSampleOn
(interfaceMakeReferenceSample,
"On",
"Generate a reference sample.",
true);
static SwitchOption interfaceMakeReferenceSampleOff
(interfaceMakeReferenceSample,
"Off",
"Do not generate a reference sample.",
false);
interfaceMakeReferenceSample.rank(-1);
}
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxReference.cc b/MatrixElement/Matchbox/Phasespace/MatchboxReference.cc
--- a/MatrixElement/Matchbox/Phasespace/MatchboxReference.cc
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxReference.cc
@@ -1,107 +1,105 @@
// -*- C++ -*-
//
// MatchboxReference.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 MatchboxReference class.
//
#include "MatchboxReference.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
MatchboxReference::MatchboxReference() {}
MatchboxReference::~MatchboxReference() {}
IBPtr MatchboxReference::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxReference::fullclone() const {
return new_ptr(*this);
}
double MatchboxReference::generateTwoToNKinematics(const double*,
vector<Lorentz5Momentum>& momenta) {
map<cPDVector,ifstream*>::iterator ref =
referenceSamples.find(mePartonData());
if ( ref == referenceSamples.end() ) {
ostringstream refname;
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
refname << (**p).PDGName();
}
refname << ".rambo";
referenceSamples[mePartonData()] = new ifstream(refname.str().c_str());
ref = referenceSamples.find(mePartonData());
}
assert(ref != referenceSamples.end());
ifstream& in = *(ref->second);
assert(in);
double x1,x2;
double x,y,z,t,m;
double weight;
in >> x1 >> x2;
for ( vector<Lorentz5Momentum>::iterator p = momenta.begin();
p != momenta.end(); ++p ) {
in >> x >> y >> z >> t >> m;
*p = Lorentz5Momentum(x*GeV,y*GeV,z*GeV,t*GeV,m*GeV);
}
in >> weight;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat((momenta[0]+momenta[1]).m2());
- fillDiagramWeights();
-
return weight;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void MatchboxReference::persistentOutput(PersistentOStream &) const {}
void MatchboxReference::persistentInput(PersistentIStream &, int) {}
// *** 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<MatchboxReference,MatchboxPhasespace>
describeHerwigMatchboxReference("Herwig::MatchboxReference", "Herwig.so");
void MatchboxReference::Init() {
static ClassDocumentation<MatchboxReference> documentation
("MatchboxReference implements reference sample phase space generation.");
}
diff --git a/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h b/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h
--- a/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h
+++ b/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h
@@ -1,528 +1,530 @@
// -*- C++ -*-
//
// MatchboxXComb.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_LastMatchboxXCombInfo_H
#define Herwig_LastMatchboxXCombInfo_H
//
// This is the declaration of the MatchboxXComb class.
//
#include "Herwig/MatrixElement/Matchbox/Utility/MatchboxXCombData.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Provide easy access to MatchboxXComb XComb extensions
*/
class LastMatchboxXCombInfo {
public:
/**
* Default constructor
*/
LastMatchboxXCombInfo()
: theLastMatchboxXComb(0), theLastHeadMatchboxXComb(0) {}
/**
* Return a pointer to the last selected XComb.
*/
MatchboxXCombData* lastMatchboxXComb() const { return theLastMatchboxXComb; }
/**
* If the last selected XComb object belongs to a
* group of XComb's return a pointer to the head
* XComb object for this group.
*/
MatchboxXCombData* lastHeadMatchboxXComb() const { return theLastHeadMatchboxXComb; }
public:
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
const vector<int>& crossingMap() const { return lastMatchboxXComb()->crossingMap(); }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
const map<size_t,size_t>& amplitudeToColourMap() const { return lastMatchboxXComb()->amplitudeToColourMap(); }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
const map<size_t,size_t>& colourToAmplitudeMap() const { return lastMatchboxXComb()->colourToAmplitudeMap(); }
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
double crossingSign() const { return lastMatchboxXComb()->crossingSign(); }
/**
* The last renormalization scale
*/
Energy2 lastRenormalizationScale() const { return lastMatchboxXComb()->lastRenormalizationScale(); }
/**
* The amplitude parton data.
*/
const cPDVector& amplitudePartonData() const { return lastMatchboxXComb()->amplitudePartonData(); }
/**
* The crossed momenta
*/
const vector<Lorentz5Momentum>& amplitudeMomenta() const { return lastMatchboxXComb()->amplitudeMomenta(); }
/**
* True, if the the tree level amplitudes need to be calculated
*/
bool calculateTreeAmplitudes() const { return lastMatchboxXComb()->calculateTreeAmplitudes(); }
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastAmplitudes() const { return lastMatchboxXComb()->lastAmplitudes(); }
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastLargeNAmplitudes() const { return lastMatchboxXComb()->lastLargeNAmplitudes(); }
/**
* True, if the the one-loop amplitudes need to be calculated
*/
bool calculateOneLoopAmplitudes() const { return lastMatchboxXComb()->calculateOneLoopAmplitudes(); }
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastOneLoopAmplitudes() const { return lastMatchboxXComb()->lastOneLoopAmplitudes(); }
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool calculateTreeME2() const { return lastMatchboxXComb()->calculateTreeME2(); }
/**
* The last tree-level matrix element squared
*/
double lastTreeME2() const { return lastMatchboxXComb()->lastTreeME2(); }
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool calculateLargeNME2() const { return lastMatchboxXComb()->calculateLargeNME2(); }
/**
* The last tree-level matrix element squared
*/
double lastLargeNME2() const { return lastMatchboxXComb()->lastLargeNME2(); }
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool calculateOneLoopInterference() const { return lastMatchboxXComb()->calculateOneLoopInterference(); }
/**
* The last one-loop/tree-level interference.
*/
double lastOneLoopInterference() const { return lastMatchboxXComb()->lastOneLoopInterference(); }
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool calculateOneLoopPoles() const { return lastMatchboxXComb()->calculateOneLoopPoles(); }
/**
* The last one-loop/tree-level interference.
*/
pair<double,double> lastOneLoopPoles() const { return lastMatchboxXComb()->lastOneLoopPoles(); }
/**
* True, if the indexed colour correlated matrix element needs to be
* calculated.
*/
bool calculateColourCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->calculateColourCorrelator(ij); }
/**
* The colour correlated matrix element.
*/
double lastColourCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->lastColourCorrelator(ij); }
/**
* True, if the indexed large-N colour correlated matrix element needs to be
* calculated.
*/
bool calculateLargeNColourCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->calculateLargeNColourCorrelator(ij); }
/**
* The large-N colour correlated matrix element.
*/
double lastLargeNColourCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->lastLargeNColourCorrelator(ij); }
/**
* True, if the indexed colour/spin correlated matrix element needs to be
* calculated.
*/
bool calculateColourSpinCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->calculateColourSpinCorrelator(ij); }
/**
* The colour/spin correlated matrix element.
*/
Complex lastColourSpinCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->lastColourSpinCorrelator(ij); }
/**
* True, if the indexed spin correlated matrix element needs to be
* calculated.
*/
bool calculateSpinCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->calculateSpinCorrelator(ij); }
/**
* The spin correlated matrix element.
*/
Complex lastSpinCorrelator(const pair<int,int>& ij) const { return lastMatchboxXComb()->lastSpinCorrelator(ij); }
/**
* Return the number of light flavours to be considered for this process.
*/
unsigned int nLight() const { return lastMatchboxXComb()->nLight(); }
/**
* Return the vector that contains the PDG ids of
* the light flavours, which are contained in the
* jet particle group.
*/
vector<long> nLightJetVec() const { return lastMatchboxXComb()->nLightJetVec(); }
/**
* Return the vector that contains the PDG ids of
* the heavy flavours, which are contained in the
* jet particle group.
*/
vector<long> nHeavyJetVec() const { return lastMatchboxXComb()->nHeavyJetVec(); }
/**
* Return the vector that contains the PDG ids of
* the light flavours, which are contained in the
* proton particle group.
*/
vector<long> nLightProtonVec() const { return lastMatchboxXComb()->nLightProtonVec(); }
/**
* Get the dimensionality of the colour basis for this process.
*/
size_t colourBasisDim() const { return lastMatchboxXComb()->colourBasisDim(); }
/**
* Return the number of degrees of freedom required by the phase space generator
*/
int nDimPhasespace() const { return lastMatchboxXComb()->nDimPhasespace(); }
/**
* Return the number of degrees of freedom required by the amplitude
*/
int nDimAmplitude() const { return lastMatchboxXComb()->nDimAmplitude(); }
/**
* Return the number of degrees of freedom required by the insertion operators
*/
int nDimInsertions() const { return lastMatchboxXComb()->nDimInsertions(); }
/**
* Get the additional random numbers required by the amplitude
*/
const vector<double>& amplitudeRandomNumbers() const { return lastMatchboxXComb()->amplitudeRandomNumbers(); }
/**
* Get the additional random numbers required by the insertion operator
*/
const vector<double>& insertionRandomNumbers() const { return lastMatchboxXComb()->insertionRandomNumbers(); }
/**
* Return the diagram weights indexed by diagram id.
*/
const map<int,double>& diagramWeights() const { return lastMatchboxXComb()->diagramWeights(); }
/**
* Return the singular limits
*/
const set<pair<size_t,size_t> >& singularLimits() const { return lastMatchboxXComb()->singularLimits(); }
/**
* Return the last matched singular limit.
*/
const set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() const { return lastMatchboxXComb()->lastSingularLimit(); }
/**
* Get the Herwig StandardModel object
*/
Ptr<StandardModel>::tcptr hwStandardModel() const { return lastMatchboxXComb()->hwStandardModel(); }
/**
* Return the symmetry factor
*/
double symmetryFactor() const { return lastMatchboxXComb()->symmetryFactor(); }
/**
* Return the OLP process id
*/
const vector<int>& olpId() const { return lastMatchboxXComb()->olpId(); }
/**
* Return the olp momentum vector
*/
double* olpMomenta() const { return lastMatchboxXComb()->olpMomenta(); }
/**
* Fill the olp momentum vector
*/
void fillOLPMomenta(const vector<Lorentz5Momentum>& mm,
const cPDVector& mePartonData,
const map<long,Energy>& reshuffleMap) const {
lastMatchboxXComb()->fillOLPMomenta(mm,mePartonData,reshuffleMap);
}
protected:
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
vector<int>& crossingMap() { return lastMatchboxXComb()->crossingMap(); }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t>& amplitudeToColourMap() { return lastMatchboxXComb()->amplitudeToColourMap(); }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t>& colourToAmplitudeMap() { return lastMatchboxXComb()->colourToAmplitudeMap(); }
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
void crossingSign(double c) { lastMatchboxXComb()->crossingSign(c); }
/**
* The last renormalization scale
*/
void lastRenormalizationScale(Energy2 lrs) { lastMatchboxXComb()->lastRenormalizationScale(lrs); }
/**
* The amplitude parton data.
*/
cPDVector& amplitudePartonData() { return lastMatchboxXComb()->amplitudePartonData(); }
/**
* The crossed momenta
*/
vector<Lorentz5Momentum>& amplitudeMomenta() { return lastMatchboxXComb()->amplitudeMomenta(); }
/**
* True, if the the tree level amplitudes need to be calculated
*/
void haveTreeAmplitudes(bool f = true) { lastMatchboxXComb()->haveTreeAmplitudes(f); }
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastAmplitudes() { return lastMatchboxXComb()->lastAmplitudes(); }
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastLargeNAmplitudes() { return lastMatchboxXComb()->lastLargeNAmplitudes(); }
/**
* True, if the the one-loop amplitudes need to be calculated
*/
void haveOneLoopAmplitudes(bool f = true) { lastMatchboxXComb()->haveOneLoopAmplitudes(f); }
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastOneLoopAmplitudes() { return lastMatchboxXComb()->lastOneLoopAmplitudes(); }
/**
* The last tree-level matrix element squared
*/
void lastTreeME2(double v) const { lastMatchboxXComb()->lastTreeME2(v); }
/**
* The last tree-level matrix element squared
*/
void lastLargeNME2(double v) const { lastMatchboxXComb()->lastLargeNME2(v); }
/**
* The last one-loop/tree-level interference.
*/
void lastOneLoopInterference(double v) const { lastMatchboxXComb()->lastOneLoopInterference(v); }
/**
* The last one-loop/tree-level interference.
*/
void lastOneLoopPoles(pair<double,double> v) const { lastMatchboxXComb()->lastOneLoopPoles(v); }
/**
* The colour correlated matrix element.
*/
void lastColourCorrelator(const pair<int,int>& ij, double v) const { lastMatchboxXComb()->lastColourCorrelator(ij,v); }
/**
* The large-N colour correlated matrix element.
*/
void lastLargeNColourCorrelator(const pair<int,int>& ij, double v) const { lastMatchboxXComb()->lastLargeNColourCorrelator(ij,v); }
/**
* The colour/spin correlated matrix element.
*/
void lastColourSpinCorrelator(const pair<int,int>& ij, Complex v) const { lastMatchboxXComb()->lastColourSpinCorrelator(ij,v); }
/**
* The spin correlated matrix element.
*/
void lastSpinCorrelator(const pair<int,int>& ij, Complex v) const { lastMatchboxXComb()->lastSpinCorrelator(ij,v); }
/**
* Set the number of light flavours to be considered for this process.
*/
void nLight(unsigned int n) { lastMatchboxXComb()->nLight(n); }
/**
* Set the elements of the vector that contains the PDG
* ids of the light flavours, which are contained in the
* jet particle group.
*/
void nLightJetVec(int n) { lastMatchboxXComb()->nLightJetVec(n); }
/**
* Set the elements of the vector that contains the PDG
* ids of the heavy flavours, which are contained in the
* jet particle group.
*/
void nHeavyJetVec(int n) { lastMatchboxXComb()->nHeavyJetVec(n); }
/**
* Set the elements of the vector that contains the PDG
* ids of the light flavours, which are contained in the
* proton particle group.
*/
void nLightProtonVec(int n) { lastMatchboxXComb()->nLightProtonVec(n); }
/**
* Set the dimensionality of the colour basis for this process.
*/
void colourBasisDim(size_t d) { lastMatchboxXComb()->colourBasisDim(d); }
/**
* Set the number of degrees of freedom required by the phase space generator
*/
void nDimPhasespace(int d) { lastMatchboxXComb()->nDimPhasespace(d); }
/**
* Set the number of degrees of freedom required by the amplitude
*/
void nDimAmplitude(int d) { lastMatchboxXComb()->nDimAmplitude(d); }
/**
* Set the number of degrees of freedom required by the insertion operators
*/
void nDimInsertions(int d) { lastMatchboxXComb()->nDimInsertions(d); }
/**
* Access the additional random numbers required by the amplitude
*/
vector<double>& amplitudeRandomNumbers() { return lastMatchboxXComb()->amplitudeRandomNumbers(); }
/**
* Access the additional random numbers required by the insertion operator
*/
vector<double>& insertionRandomNumbers() { return lastMatchboxXComb()->insertionRandomNumbers(); }
/**
* Access the diagram weights indexed by diagram id.
*/
map<int,double>& diagramWeights() { return lastMatchboxXComb()->diagramWeights(); }
/**
* Access the singular limits
*/
set<pair<size_t,size_t> >& singularLimits() { return lastMatchboxXComb()->singularLimits(); }
/**
* Access the last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() { return lastMatchboxXComb()->lastSingularLimit(); }
/**
* Set the Herwig StandardModel object
*/
void hwStandardModel(Ptr<StandardModel>::tcptr sm) { lastMatchboxXComb()->hwStandardModel(sm); }
/**
* Set the symmetry factor
*/
void symmetryFactor(double f) const { lastMatchboxXComb()->symmetryFactor(f); }
/**
* Set the OLP process id
*/
void olpId(int pType, int id) { lastMatchboxXComb()->olpId(pType,id); }
protected:
/**
* Set the XComb pointer cast to MatchboxXComb
*/
void lastMatchboxXComb(tStdXCombPtr xc) {
- theLastMatchboxXComb = xc ? &dynamic_cast<MatchboxXCombData&>(*xc) : 0;
+ theLastMatchboxXComb = xc ?
+ dynamic_cast<MatchboxXCombData*>(PtrTraits<tStdXCombPtr>::barePointer(xc)) : 0;
theLastHeadMatchboxXComb =
- xc && xc->head() ? &dynamic_cast<MatchboxXCombData&>(*xc->head()) : 0;
+ xc && xc->head() ?
+ dynamic_cast<MatchboxXCombData*>(PtrTraits<tStdXCombPtr>::barePointer(xc->head())) : 0;
}
/**
* The XComb pointer cast to MatchboxXComb
*/
MatchboxXCombData* theLastMatchboxXComb;
/**
* The head XComb pointer cast to MatchboxXComb
*/
MatchboxXCombData* theLastHeadMatchboxXComb;
};
}
#endif // Herwig_LastMatchboxXCombInfo_H

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:45 PM (1 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023860
Default Alt Text
(130 KB)

Event Timeline