Page MenuHomeHEPForge
Files F7878469

No OneTemporary
Actions

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.cc b/MatrixElement/Matchbox/Base/SubtractedME.cc
--- a/MatrixElement/Matchbox/Base/SubtractedME.cc
+++ b/MatrixElement/Matchbox/Base/SubtractedME.cc
@@ -1,692 +1,702 @@
// -*- C++ -*-
//
// SubtractedME.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SubtractedME class.
//
#include "SubtractedME.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
using namespace Herwig;
SubtractedME::SubtractedME()
: MEGroup(),
theSubtractionData(""),
theVerbose(false),
theSubProcessGroups(false), theInclusive(false),
theVetoScales(false),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false) {}
SubtractedME::~SubtractedME() {}
IBPtr SubtractedME::clone() const {
return new_ptr(*this);
}
IBPtr SubtractedME::fullclone() const {
return new_ptr(*this);
}
void SubtractedME::setXComb(tStdXCombPtr xc) {
MEGroup::setXComb(xc);
}
MEBase::DiagramVector SubtractedME::dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const {
Ptr<SubtractionDipole>::tptr dipole =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(depME);
if ( !dipole ) {
Throw<InitException>() << "A dependent matrix element of SubtractedME "
<< "has not been derived from SubtractionDipole. "
<< "Please check the corresponding input file.";
}
return dipole->underlyingBornDiagrams(proc);
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::dipoles() {
if ( allDipoles().empty() ) {
allDipoles() = DipoleRepository::dipoles();
}
if ( dependent().empty() )
getDipoles();
vector<Ptr<SubtractionDipole>::ptr> res;
for ( MEVector::const_iterator k = dependent().begin();
k != dependent().end(); ++k )
res.push_back(dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(*k));
return res;
}
void SubtractedME::getDipoles() {
if ( !dependent().empty() )
return;
if ( allDipoles().empty() ) {
allDipoles() = DipoleRepository::dipoles();
}
Ptr<MatchboxMEBase>::tptr real =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( allDipoles().empty() || theBorns.empty() || !real )
Throw<InitException>() << "The SubtractedME '"
<< name() << "' could not generate "
<< "subtraction terms for the real emission "
<< "matrix element '" << real->name() << "'. "
<< "Please check the corresponding input file.";
Ptr<MatchboxMEBase>::ptr myRealEmissionME = real->cloneMe();
ostringstream pname;
pname << fullName() << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myRealEmissionME->cloneDependencies(pname.str());
head(myRealEmissionME);
real = myRealEmissionME;
MEVector dipMEs;
vector<Ptr<SubtractionDipole>::ptr> genDipoles
= real->getDipoles(allDipoles(),theBorns);
if ( theSubtractionData != "" ) {
theSubProcessGroups = true;
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
genDipoles.begin(); d != genDipoles.end(); ++d )
(**d).doTestSubtraction();
}
if ( genDipoles.empty() ) {
// probably finite real contribution, but warn
generator()->log() << "\nWarning: No subtraction dipoles could be found for the processes:\n";
for ( vector<PDVector>::const_iterator s = real->subProcesses().begin();
s != real->subProcesses().end(); ++s ) {
generator()->log() << (*s)[0]->PDGName() << " " << (*s)[1]->PDGName()
<< " -> ";
for ( PDVector::const_iterator p = s->begin() + 2; p != s->end(); ++p )
generator()->log() << (**p).PDGName() << " ";
generator()->log() << "\n" << flush;
}
generator()->log() << "Assuming finite tree-level O(alphaS) correction.\n";
}
dipMEs.resize(genDipoles.size());
copy(genDipoles.begin(),genDipoles.end(),dipMEs.begin());
dependent() = dipMEs;
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::splitDipoles(const cPDVector& born) {
vector<Ptr<SubtractionDipole>::ptr> dips = dipoles();
vector<Ptr<SubtractionDipole>::ptr> res;
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = dips.begin();
d != dips.end(); ++d ) {
for ( DiagramVector::const_iterator p = (**d).underlyingBornME()->diagrams().begin();
p != (**d).underlyingBornME()->diagrams().end(); ++p )
if ( born == (**p).partons() ) {
res.push_back(*d);
break;
}
}
return res;
}
void SubtractedME::showerApproximation(Ptr<ShowerApproximation>::tptr app) {
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->showerApproximation(app);
}
}
+void SubtractedME::doRealEmissionScales() {
+ for ( MEVector::const_iterator m = dependent().begin();
+ m != dependent().end(); ++m ) {
+ Ptr<SubtractionDipole>::tptr dip =
+ dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
+ assert(dip);
+ dip->doRealEmissionScales();
+ }
+}
+
void SubtractedME::doRealShowerSubtraction() {
theRealShowerSubtraction = true;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->doRealShowerSubtraction();
}
}
void SubtractedME::doVirtualShowerSubtraction() {
theVirtualShowerSubtraction = true;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->doVirtualShowerSubtraction();
}
}
void SubtractedME::setVetoScales(tSubProPtr) const {}
void SubtractedME::fillProjectors() {
if ( !inclusive() && !virtualShowerSubtraction() )
return;
Ptr<StdXCombGroup>::tptr group =
dynamic_ptr_cast<Ptr<StdXCombGroup>::tptr>(lastXCombPtr());
for ( vector<StdXCombPtr>::const_iterator d = group->dependent().begin();
d != group->dependent().end(); ++d ) {
if ( (**d).lastCrossSection() != ZERO )
lastXCombPtr()->projectors().insert(1.,*d);
}
}
double SubtractedME::reweightHead(const vector<tStdXCombPtr>& dep) {
if ( inclusive() && !lastXComb().lastProjector() )
return 1.;
if ( virtualShowerSubtraction() && !lastXComb().lastProjector() ) {
assert(!showerApproximation()->belowCutoff());
return 0.;
}
if ( realShowerSubtraction() ) {
assert(showerApproximation());
bool below = showerApproximation()->belowCutoff();
if ( below )
return 0.;
return 1.;
}
if ( virtualShowerSubtraction() || inclusive() ) {
if ( virtualShowerSubtraction() ) {
assert(showerApproximation());
bool above = !showerApproximation()->belowCutoff();
if ( above )
return 0.;
}
double sum = 0.;
size_t n = 0;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( (**d).lastCrossSection() != ZERO ) {
sum += (**d).lastME2();
++n;
}
}
return
n * lastXComb().lastProjector()->lastME2() / sum;
}
return 1.;
}
double SubtractedME::reweightDependent(tStdXCombPtr xc, const vector<tStdXCombPtr>& dep) {
if ( inclusive() && !lastXComb().lastProjector() )
return 0.;
if ( virtualShowerSubtraction() && !lastXComb().lastProjector() ) {
assert(!showerApproximation()->belowCutoff());
return 0.;
}
if ( virtualShowerSubtraction() || inclusive() ) {
if ( xc != lastXComb().lastProjector() )
return 0.;
size_t n = 0;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( (**d).lastCrossSection() != ZERO ) {
++n;
}
}
return n;
}
return 1.;
}
void SubtractedME::doinit() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinit();
return;
}
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theReal ) {
getDipoles();
}
if ( theVerbose )
print(Repository::clog());
MEGroup::doinit();
}
void SubtractedME::doinitrun() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinitrun();
return;
}
MEGroup::doinitrun();
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theSubtractionData != "" ) {
set<cPDVector> procs;
for ( DiagramVector::const_iterator d = head()->diagrams().begin();
d != head()->diagrams().end(); ++d ) {
if ( procs.find((**d).partons()) == procs.end() )
procs.insert((**d).partons());
}
for ( set<cPDVector>::const_iterator p = procs.begin();
p != procs.end(); ++p ) {
for ( size_t i = 0; i < (*p).size(); ++i ) {
if ( !(*p)[i]->coloured() )
continue;
if ( i > 1 &&
(*p)[i]->id() == ParticleID::g ) {
softHistograms[SoftSubtractionIndex(*p,i)] = SubtractionHistogram(0.001,10.);
if ( theReal->phasespace() )
theReal->phasespace()->singularLimit(i,i);
}
for ( size_t j = i+1; j < (*p).size(); ++j ) {
if ( !(*p)[j]->coloured() )
continue;
long iid = (*p)[i]->id();
long jid = (*p)[j]->id();
if ( i < 2 && j < 2 )
continue;
if ( i < 2 && j > 1 ) {
if ( abs(iid) < 7 && abs(jid) < 7 && iid != jid )
continue;
}
if ( i > 1 && j > 1 ) {
if ( abs(iid) < 7 && abs(jid) < 7 && iid + jid != 0 )
continue;
}
bool haveDipole = false;
for ( MEVector::const_iterator k = dependent().begin();
k != dependent().end(); ++k ) {
const SubtractionDipole& dip = dynamic_cast<const SubtractionDipole&>(**k);
if ( ( (size_t)(dip.realEmitter()) == i && (size_t)(dip.realEmission()) == j ) ||
( (size_t)(dip.realEmitter()) == j && (size_t)(dip.realEmission()) == i ) ) {
haveDipole = true;
break;
}
}
if ( !haveDipole )
continue;
collinearHistograms[CollinearSubtractionIndex(*p,make_pair(i,j))] = SubtractionHistogram(0.001,20.);
if ( theReal->phasespace() )
theReal->phasespace()->singularLimit(i,j);
}
}
}
}
}
void SubtractedME::dofinish() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::dofinish();
return;
}
MEGroup::dofinish();
for ( map<CollinearSubtractionIndex,SubtractionHistogram>::
const_iterator b = collinearHistograms.begin();
b != collinearHistograms.end(); ++b ) {
b->second.dump(theSubtractionData,
b->first.first,
b->first.second.first,
b->first.second.second);
}
for ( map<SoftSubtractionIndex,SubtractionHistogram>::
const_iterator b = softHistograms.begin();
b != softHistograms.end(); ++b ) {
b->second.dump(theSubtractionData,
b->first.first,
b->first.second,
b->first.second);
}
}
void SubtractedME::rebind(const TranslationMap & trans) {
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
*d = trans.translate(*d);
MEGroup::rebind(trans);
}
IVector SubtractedME::getReferences() {
IVector ret = MEGroup::getReferences();
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
theDipoles.begin(); d != theDipoles.end(); ++d )
ret.push_back(*d);
return ret;
}
void SubtractedME::print(ostream& os) const {
os << "--- SubtractedME setup ---------------------------------------------------------\n";
os << " '" << name() << "' subtracting real emission\n '"
<< head()->name() << "' using the dipoles:\n";
for ( MEVector::const_iterator d = dependent().begin();
d != dependent().end(); ++d )
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*d)->print(os);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractedME::printLastEvent(ostream& os) const {
os << "--- SubtractedME last event information ----------------------------------------\n";
os << " for subtracted matrix element '" << name() << "'\n";
os << " real emission event information:\n";
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head())->printLastEvent(os);
os << " dipoles event information:\n";
for ( MEVector::const_iterator d = dependent().begin();
d != dependent().end(); ++d )
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*d)->printLastEvent(os);
os << "--- end SubtractedME last event information ------------------------------------\n\n\n";
os << flush;
}
void SubtractedME::lastEventStatistics() {
MEGroup::lastEventStatistics();
if ( !generator() )
return;
/*
if ( theVerbose )
printLastEvent(generator()->log());
*/
if ( !collinearHistograms.empty() )
lastEventSubtraction();
}
SubtractedME::SubtractionHistogram::
SubtractionHistogram(double low,
double up,
unsigned int nbins)
: lower(low) {
nbins = nbins + 1;
double c = log10(up/low) / (nbins-1.);
for ( unsigned int k = 1; k < nbins; ++k ) {
bins[low*pow(10.0,k*c)] = make_pair(Constants::MaxDouble,0.);
}
}
void SubtractedME::SubtractionHistogram::
dump(const std::string& prefix,
const cPDVector& proc,
int i, int j) const {
ostringstream fname("");
for ( cPDVector::const_iterator p = proc.begin();
p != proc.end(); ++p )
fname << (**p).PDGName();
fname << "-" << i << "-" << j;
ofstream out((prefix+fname.str()+".dat").c_str());
for ( map<double,pair<double,double> >::const_iterator b = bins.begin();
b != bins.end(); ++b ) {
map<double,pair<double,double> >::const_iterator bp = b; --bp;
if ( b->second.first != Constants::MaxDouble &&
b->second.second != 0. ) {
if ( b != bins.begin() )
out << bp->first;
else
out << lower;
out << " " << b->first
<< " " << b->second.first
<< " " << b->second.second
<< "\n" << flush;
}
}
ofstream gpout((prefix+fname.str()+".gp").c_str());
gpout << "set terminal epslatex color solid;\n"
<< "set output '" << fname.str() << "-plot.tex';\n"
<< "set log x;\n"
<< "set size 0.5,0.6;\n"
<< "set yrange [0:2];\n"
<< "set xrange [0.001:10];\n";
if ( i != j ) {
gpout << "set xlabel '$\\sqrt{s_{" << i << j << "}}/{\\rm GeV}$'\n";
} else {
gpout << "set xlabel '$E_{" << i << "}/{\\rm GeV}$'\n";
}
gpout << "plot 1 w lines lc rgbcolor \"#DDDDDD\" notitle, '" << fname.str()
<< ".dat' u (($1+$2)/2.):3:($4 < 4. ? $4 : 4.) w filledcurves lc rgbcolor \"#00AACC\" t "
<< "'$";
for ( size_t k = 0; k < proc.size(); k++ ) {
if ( k == 2 )
gpout << "\\to ";
gpout << (proc[k]->id() < 0 ? "\\bar{" : "")
<< (proc[k]->id() < 0 ? proc[k]->CC()->PDGName() : proc[k]->PDGName())
<< (proc[k]->id() < 0 ? "}" : "") << " ";
}
gpout << "$';\n";
gpout << "reset;\n";
}
void SubtractedME::lastEventSubtraction() {
tStdXCombGroupPtr xc = dynamic_ptr_cast<tStdXCombGroupPtr>(lastXCombPtr());
CrossSection xcme2 = xc->lastHeadCrossSection();
CrossSection xcdip = ZERO;
if ( xcme2 == ZERO )
return;
for ( vector<StdXCombPtr>::const_iterator d = xc->dependent().begin();
d != xc->dependent().end(); ++d ) {
if ( !(*d) )
continue;
if ( !(**d).matrixElement()->apply() )
continue;
xcdip += (**d).lastCrossSection();
}
if ( theReal->phasespace() ) {
size_t i = theReal->phasespace()->lastSingularLimit().first;
size_t j = theReal->phasespace()->lastSingularLimit().second;
if ( i == j &&
softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,abs(xcdip)/abs(xcme2));
}
if ( i != j &&
collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j)))
!= collinearHistograms.end() ) {
double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV;
collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))].
book(s,abs(xcdip)/abs(xcme2));
}
return;
}
for ( size_t i = 0; i < meMomenta().size(); ++i ) {
if ( i > 1 ) {
if ( softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,abs(xcdip)/abs(xcme2));
}
}
for ( size_t j = i+1; j < meMomenta().size(); ++j ) {
if ( collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j)))
== collinearHistograms.end() )
continue;
double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV;
collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))].
book(s,abs(xcdip)/abs(xcme2));
}
}
}
void SubtractedME::persistentOutput(PersistentOStream & os) const {
os << theDipoles << theBorns << theSubtractionData
<< theVerbose << theInclusive << theSubProcessGroups << theVetoScales
<< theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction;
}
void SubtractedME::persistentInput(PersistentIStream & is, int) {
is >> theDipoles >> theBorns >> theSubtractionData
>> theVerbose >> theInclusive >> theSubProcessGroups >> theVetoScales
>> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction;
}
void SubtractedME::Init() {
static ClassDocumentation<SubtractedME> documentation
("SubtractedME represents a subtracted real emission matrix element.");
static RefVector<SubtractedME,MatchboxMEBase> interfaceBorns
("Borns",
"The underlying Born matrix elements to be considered",
&SubtractedME::theBorns, -1, false, false, true, true, false);
static Parameter<SubtractedME,string> interfaceSubtractionData
("SubtractionData",
"File to dump subtraction check to.",
&SubtractedME::theSubtractionData, "",
false, false);
static Switch<SubtractedME,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&SubtractedME::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"On",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"Off",
false);
static Switch<SubtractedME,bool> interfaceSubProcessGroups
("SubProcessGroups",
"Switch on or off production of sub-process groups.",
&SubtractedME::theSubProcessGroups, false, false, false);
static SwitchOption interfaceSubProcessGroupsOn
(interfaceSubProcessGroups,
"On",
"On",
true);
static SwitchOption interfaceSubProcessGroupsOff
(interfaceSubProcessGroups,
"Off",
"Off",
false);
static Switch<SubtractedME,bool> interfaceInclusive
("Inclusive",
"Switch on or off production of inclusive cross section.",
&SubtractedME::theInclusive, false, false, false);
static SwitchOption interfaceInclusiveOn
(interfaceInclusive,
"On",
"On",
true);
static SwitchOption interfaceInclusiveOff
(interfaceInclusive,
"Off",
"Off",
false);
static Switch<SubtractedME,bool> interfaceVetoScales
("VetoScales",
"Switch on or off production of sub-process groups.",
&SubtractedME::theVetoScales, false, false, false);
static SwitchOption interfaceVetoScalesOn
(interfaceVetoScales,
"On",
"On",
true);
static SwitchOption interfaceVetoScalesOff
(interfaceVetoScales,
"Off",
"Off",
false);
static Reference<SubtractedME,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&SubtractedME::theShowerApproximation, false, false, true, true, false);
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<SubtractedME,MEGroup>
describeHerwigSubtractedME("Herwig::SubtractedME", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.h b/MatrixElement/Matchbox/Base/SubtractedME.h
--- a/MatrixElement/Matchbox/Base/SubtractedME.h
+++ b/MatrixElement/Matchbox/Base/SubtractedME.h
@@ -1,534 +1,539 @@
// -*- C++ -*-
//
// SubtractedME.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_SubtractedME_H
#define HERWIG_SubtractedME_H
//
// This is the declaration of the SubtractedME class.
//
#include "ThePEG/MatrixElement/MEGroup.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief SubtractedME represents a subtracted real emission matrix element.
*
* @see \ref SubtractedMEInterfaces "The interfaces"
* defined for SubtractedME.
*/
class SubtractedME: public MEGroup {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SubtractedME();
/**
* The destructor.
*/
virtual ~SubtractedME();
//@}
public:
/** @name Phasespace and subprocess information */
//@{
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr);
/**
* Return true, if the same additional random numbers
* should be presented to any of the dependent
* matrix elements.
*/
virtual bool uniformAdditional() const { return true; }
/**
* Return true, if the XComb steering this matrix element
* should keep track of the random numbers used to generate
* the last phase space point
*/
virtual bool keepRandomNumbers() const { return true; }
/**
* Given a process from the head matrix element,
* return a list of diagrams which should be considered for
* the given dependent matrix element.
*/
virtual MEBase::DiagramVector dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const;
/**
* Return true, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
virtual bool subProcessGroups() const { return theSubProcessGroups; }
/**
* Switch on or off producing subprocess groups.
*/
void setSubProcessGroups(bool on = true) { theSubProcessGroups = on; }
/**
* Return true, if one of the dependent subprocesses should be
* constructed in place of the one driven by the head matrix element
* or a full subprocess group.
*/
virtual bool selectDependentSubProcess() const { return theInclusive; }
/**
* Return true, if the integral over the unresolved emission should be
* calculated.
*/
bool inclusive() const { return theInclusive; }
/**
* Switch on or off inclusive mode.
*/
void setInclusive(bool on = true) { theInclusive = on; }
/**
* Fill the projectors object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
virtual void fillProjectors();
/**
* Return true, if this MEGroup will reweight the contributing cross
* sections.
*/
virtual bool groupReweighted() const { return inclusive() || showerApproximation(); }
/**
* Reweight the head cross section
*/
virtual double reweightHead(const vector<tStdXCombPtr>&);
/**
* Reweight the dependent cross section
*/
virtual double reweightDependent(tStdXCombPtr, const vector<tStdXCombPtr>&);
+ /**
+ * Switch on or off that scales should be calculated from real emission kinematics
+ */
+ void doRealEmissionScales();
+
//@}
/** @name Methods relevant to matching */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches() {
MEGroup::flushCaches();
if ( showerApproximation() )
showerApproximation()->resetBelowCutoff();
}
/**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr);
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
/**
* Indicate that the shower real emission contribution should be subtracted.
*/
void doRealShowerSubtraction();
/**
* Return true, if the shower real emission contribution should be subtracted.
*/
bool realShowerSubtraction() const { return theRealShowerSubtraction; }
/**
* Indicate that the shower virtual contribution should be subtracted.
*/
void doVirtualShowerSubtraction();
/**
* Return true, if the shower virtual contribution should be subtracted.
*/
bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
//@}
/** @name Matrix element and dipole information */
//@{
/**
* Return the subtraction dipoles.
*/
vector<Ptr<SubtractionDipole>::ptr> dipoles();
/**
* Return the underlying born matrix elements.
*/
const vector<Ptr<MatchboxMEBase>::ptr>& borns() const { return theBorns; }
/**
* Access the underlying born matrix elements.
*/
vector<Ptr<MatchboxMEBase>::ptr>& borns() { return theBorns; }
/**
* Build up dipoles needed.
*/
void getDipoles();
/**
* Return all dipoles matching the given Born process
*/
vector<Ptr<SubtractionDipole>::ptr> splitDipoles(const cPDVector&);
/**
* Return the subtraction dipoles.
*/
const vector<Ptr<SubtractionDipole>::ptr>& allDipoles() const { return theDipoles; }
/**
* Access the subtraction dipoles.
*/
vector<Ptr<SubtractionDipole>::ptr>& allDipoles() { return theDipoles; }
//@}
/** @name Veto scale settings */
//@{
/**
* Set veto scales on the particles at the given
* SubProcess which has been generated using this
* matrix element.
*/
virtual void setVetoScales(tSubProPtr) const;
/**
* Return true, if veto scales should be set
* for the real emission
*/
bool vetoScales() const { return theVetoScales; }
/**
* Switch on setting veto scales
*/
void doVetoScales() { theVetoScales = true; }
/**
* Switch off setting veto scales
*/
void noVetoScales() { theVetoScales = true; }
//@}
/** @name Diagnostic information */
//@{
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Collect information on the last evaluated phasespace
* point for verification or debugging purposes. This
* only called, if the StdXCombGroup did accumulate
* a non-zero cross section from this ME group.
*/
virtual void lastEventStatistics();
/**
* Print debug information on the last event
*/
void printLastEvent(ostream&) const;
/**
* Check the subtraction for the last event
*/
void lastEventSubtraction();
/**
* Set a file to print subtraction check to
*/
void subtractionData(string n) { theSubtractionData = n; }
/**
* Return true, if verbose
*/
bool verbose() const { return theVerbose; }
/**
* Switch on or off verbosity for this subtracted ME
*/
void setVerbose(bool on = true) { theVerbose = on; }
//@}
/** @name Setup of Subtracted ME objects */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return 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);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans);
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
//@}
private:
/**
* Pointer to the head real emission ME casted to a MatchboxMEBase
* object.
*/
Ptr<MatchboxMEBase>::ptr theReal;
/**
* The dipoles to be considered; the dipoles generated
* can be accessed throught the dependent() matrxi element
* vector, provided the head() is a MatchboxMEBase object.
*/
vector<Ptr<SubtractionDipole>::ptr> theDipoles;
/**
* The underlying Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theBorns;
/**
* File name to dump subtraction check to
*/
string theSubtractionData;
/**
* Simple envelope histogram to keep track of subtraction
*/
struct SubtractionHistogram {
/**
* The lower bound
*/
double lower;
/**
* The bins, indexed by upper bound.
*/
map<double,pair<double,double> > bins;
/**
* Constructor
*/
SubtractionHistogram(double low = 0.001,
double up = 10.,
unsigned int nbins = 100);
/**
* Book an event.
*/
void book(double inv, double ratio) {
map<double,pair<double,double> >::iterator b =
bins.upper_bound(inv);
if ( b == bins.end() ) return;
b->second.first = min(b->second.first,abs(ratio));
b->second.second = max(b->second.second,abs(ratio));
}
/**
* Write to file given name and invariant.
*/
void dump(const std::string& prefix,
const cPDVector& proc,
int i, int j) const;
};
/**
* Define the key for the collinear subtraction data.
*/
typedef pair<cPDVector,pair<size_t, size_t> > CollinearSubtractionIndex;
/**
* subtraction data for collinear limits.
*/
map<CollinearSubtractionIndex,SubtractionHistogram> collinearHistograms;
/**
* Define the key for the soft subtraction data.
*/
typedef pair<cPDVector,size_t> SoftSubtractionIndex;
/**
* subtraction data for soft limits.
*/
map<SoftSubtractionIndex,SubtractionHistogram> softHistograms;
/**
* Switch to print full information on the
* last phase space point.
*/
bool theVerbose;
/**
* True, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool theSubProcessGroups;
/**
* True, if the integral over the unresolved emission should be
* calculated.
*/
bool theInclusive;
/**
* True, if veto scales should be set
* for the real emission
*/
bool theVetoScales;
/**
* The shower approximation.
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* True, if the shower real emission contribution should be subtracted.
*/
bool theRealShowerSubtraction;
/**
* True, if the shower virtual contribution should be subtracted.
*/
bool theVirtualShowerSubtraction;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SubtractedME & operator=(const SubtractedME &);
};
}
#endif /* HERWIG_SubtractedME_H */
diff --git a/MatrixElement/Matchbox/Dipoles/FFMggxDipole.cc b/MatrixElement/Matchbox/Dipoles/FFMggxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FFMggxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FFMggxDipole.cc
@@ -1,134 +1,124 @@
// -*- C++ -*-
//
// FFMggxDipole.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 FFMggxDipole class.
//
#include "FFMggxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h"
using namespace Herwig;
FFMggxDipole::FFMggxDipole()
: SubtractionDipole() {}
FFMggxDipole::~FFMggxDipole() {}
IBPtr FFMggxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FFMggxDipole::fullclone() const {
return new_ptr(*this);
}
bool FFMggxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
partons[emitter]->id() == ParticleID::g &&
partons[spectator]->mass() != ZERO;
}
// TODO: (5.20)
double FFMggxDipole::me2Avg(double) const {
assert(false && "implementation missing");
return 0.;
}
double FFMggxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
// masses, g->gg all masses zero except spectator
double muj2 = sqr( realEmissionME()->lastXComb().mePartonData()[realSpectator()]->mass() / lastDipoleScale() );
// massive extra terms
double vijk = sqrt( sqr(2.*muj2+(1.-muj2)*(1.-y))-4.*muj2 ) / ((1.-muj2)*(1.-y));
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double diag = 1./(1-z*(1-y))+1./(1-(1.-z)*(1.-y))-2./vijk; // kappa=0
double zim = z-0.5*(1.-vijk), zjm = (1.-z)-0.5*(1.-vijk);
Lorentz5Momentum pc =
zim*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
zjm*realEmissionME()->lastXComb().meMomenta()[realEmission()];
SpinCorrelationTensor corr(-diag,pc,prop/2.*vijk);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
// extra mass terms all = 0.
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FFMggxDipole::persistentOutput(PersistentOStream &) const {
}
void FFMggxDipole::persistentInput(PersistentIStream &, int) {
}
void FFMggxDipole::Init() {
static ClassDocumentation<FFMggxDipole> documentation
("FFMggxDipole");
DipoleRepository::registerDipole<FFMggxDipole,FFMassiveTildeKinematics,FFMassiveInvertedTildeKinematics>
("FFMggxDipole","FFMassiveTildeKinematics","FFMassiveInvertedTildeKinematics");
}
// *** 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<FFMggxDipole,SubtractionDipole>
describeHerwigFFMggxDipole("Herwig::FFMggxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FFMqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/FFMqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FFMqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FFMqgxDipole.cc
@@ -1,132 +1,122 @@
// -*- C++ -*-
//
// FFMqgxDipole.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 FFMqgxDipole class.
//
#include "FFMqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h"
using namespace Herwig;
FFMqgxDipole::FFMqgxDipole()
: SubtractionDipole() {}
FFMqgxDipole::~FFMqgxDipole() {}
IBPtr FFMqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FFMqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool FFMqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 7 &&
!(partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
// TODO: (5.16)
double FFMqgxDipole::me2Avg(double) const {
assert(false && "implementation missing");
return 0.;
}
double FFMqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
// masses
double muQ2 = sqr( realEmissionME()->lastXComb().mePartonData()[realEmitter()]->mass() / lastDipoleScale() );
double muj2 = sqr( realEmissionME()->lastXComb().mePartonData()[realSpectator()]->mass() / lastDipoleScale() );
// massive extra terms
double vijk = sqrt( sqr(2.*muj2+(1.-muQ2-muj2)*(1.-y))-4.*muj2 ) / ((1.-muQ2-muj2)*(1.-y));
double vbar = sqrt( 1.+sqr(muQ2)+sqr(muj2)-2.*(muQ2+muj2+muQ2*muj2) ) / (1.-muQ2-muj2);
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
// extra mass terms cancel: mi2+m2-Mi2 = mQ2+0-mQ2
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-z*(1.-y)) - vbar/vijk * ( (1.+z) + muQ2*sqr(lastDipoleScale())*2./prop ) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FFMqgxDipole::persistentOutput(PersistentOStream &) const {
}
void FFMqgxDipole::persistentInput(PersistentIStream &, int) {
}
void FFMqgxDipole::Init() {
static ClassDocumentation<FFMqgxDipole> documentation
("FFMqgxDipole");
DipoleRepository::registerDipole<FFMqgxDipole,FFMassiveTildeKinematics,FFMassiveInvertedTildeKinematics>
("FFMqgxDipole","FFMassiveTildeKinematics","FFMassiveInvertedTildeKinematics");
}
// *** 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<FFMqgxDipole,SubtractionDipole>
describeHerwigFFMqgxDipole("Herwig::FFMqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FFMqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/FFMqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FFMqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FFMqqxDipole.cc
@@ -1,138 +1,128 @@
// -*- C++ -*-
//
// FFMqqxDipole.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 FFMqqxDipole class.
//
#include "FFMqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFMassiveTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFMassiveInvertedTildeKinematics.h"
using namespace Herwig;
FFMqqxDipole::FFMqqxDipole()
: SubtractionDipole() {}
FFMqqxDipole::~FFMqqxDipole() {}
IBPtr FFMqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FFMqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool FFMqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator > 1 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() + partons[emitter]->id() == 0 &&
!(partons[emission]->mass() == ZERO &&
partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
// TODO: (5.18)
double FFMqqxDipole::me2Avg(double) const {
assert(false && "implementation missing");
return 0.;
}
double FFMqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
// masses
double muQ2 = sqr( realEmissionME()->lastXComb().mePartonData()[realEmission()]->mass() / lastDipoleScale() );
double muj2 = sqr( realEmissionME()->lastXComb().mePartonData()[realSpectator()]->mass() / lastDipoleScale() );
Energy2 mQ2 = sqr( realEmissionME()->lastXComb().mePartonData()[realEmission()]->mass() );
// massive extra terms
double vijk = sqrt( sqr(2.*muj2+(1.-2.*muQ2-muj2)*(1.-y))-4.*muj2 ) / ((1.-2.*muQ2-muj2)*(1.-y));
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double zim = z-0.5*(1.-vijk), zjm = (1.-z)-0.5*(1.-vijk);
Lorentz5Momentum pc =
zim*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
zjm*realEmissionME()->lastXComb().meMomenta()[realEmission()];
// kappa=0 -- otherwise: extra diagonal term (instead of just -1.)
SpinCorrelationTensor corr(-1.,pc,-(prop+2.*mQ2)/4.);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/ (prop+2.*mQ2);
+ (underlyingBornME()->lastXComb().lastAlphaS())/ (prop+2.*mQ2);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FFMqqxDipole::persistentOutput(PersistentOStream &) const {
}
void FFMqqxDipole::persistentInput(PersistentIStream &, int) {
}
void FFMqqxDipole::Init() {
static ClassDocumentation<FFMqqxDipole> documentation
("FFMqqxDipole");
DipoleRepository::registerDipole<FFMqqxDipole,FFMassiveTildeKinematics,FFMassiveInvertedTildeKinematics>
("FFMqqxDipole","FFMassiveTildeKinematics","FFMassiveInvertedTildeKinematics");
}
// *** 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<FFMqqxDipole,SubtractionDipole>
describeHerwigFFMqqxDipole("Herwig::FFMqqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FFggxDipole.cc b/MatrixElement/Matchbox/Dipoles/FFggxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FFggxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FFggxDipole.cc
@@ -1,161 +1,141 @@
// -*- C++ -*-
//
// FFggxDipole.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 FFggxDipole class.
//
#include "FFggxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.h"
using namespace Herwig;
FFggxDipole::FFggxDipole()
: SubtractionDipole() {}
FFggxDipole::~FFggxDipole() {}
IBPtr FFggxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FFggxDipole::fullclone() const {
return new_ptr(*this);
}
bool FFggxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
partons[emitter]->id() == ParticleID::g &&
partons[spectator]->mass() == ZERO;
}
double FFggxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double res =
1./(1.-z*(1.-y)) + 1./(1.-(1.-z)*(1.-y)) - 2 + z*(1.-z);
res *= -ccme2;
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FFggxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double diag = 1./(1.-z*(1.-y))+1./(1.-(1.-z)*(1.-y))-2.;
Lorentz5Momentum pc =
z*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
(1.-z)*realEmissionME()->lastXComb().meMomenta()[realEmission()];
SpinCorrelationTensor corr(-diag,pc,prop/2.);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FFggxDipole::persistentOutput(PersistentOStream &) const {
}
void FFggxDipole::persistentInput(PersistentIStream &, int) {
}
void FFggxDipole::Init() {
static ClassDocumentation<FFggxDipole> documentation
("FFggxDipole");
DipoleRepository::registerDipole<FFggxDipole,FFLightTildeKinematics,FFLightInvertedTildeKinematics>
("FFggxDipole","FFLightTildeKinematics","FFLightInvertedTildeKinematics");
}
// *** 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<FFggxDipole,SubtractionDipole>
describeHerwigFFggxDipole("Herwig::FFggxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FFqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/FFqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FFqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FFqgxDipole.cc
@@ -1,160 +1,140 @@
// -*- C++ -*-
//
// FFqgxDipole.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 FFqgxDipole class.
//
#include "FFqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.h"
using namespace Herwig;
FFqgxDipole::FFqgxDipole()
: SubtractionDipole() {}
FFqgxDipole::~FFqgxDipole() {}
IBPtr FFqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FFqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool FFqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 6 &&
partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double FFqgxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-z*(1.-y)) - (1.+z) );
res *= -ccme2;
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FFqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double y = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-z*(1.-y)) - (1.+z) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FFqgxDipole::persistentOutput(PersistentOStream &) const {
}
void FFqgxDipole::persistentInput(PersistentIStream &, int) {
}
void FFqgxDipole::Init() {
static ClassDocumentation<FFqgxDipole> documentation
("FFqgxDipole");
DipoleRepository::registerDipole<FFqgxDipole,FFLightTildeKinematics,FFLightInvertedTildeKinematics>
("FFqgxDipole","FFLightTildeKinematics","FFLightInvertedTildeKinematics");
}
// *** 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<FFqgxDipole,SubtractionDipole>
describeHerwigFFqgxDipole("Herwig::FFqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FFqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/FFqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FFqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FFqqxDipole.cc
@@ -1,160 +1,140 @@
// -*- C++ -*-
//
// FFqqxDipole.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 FFqqxDipole class.
//
#include "FFqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FFLightInvertedTildeKinematics.h"
using namespace Herwig;
FFqqxDipole::FFqqxDipole()
: SubtractionDipole() {}
FFqqxDipole::~FFqqxDipole() {}
IBPtr FFqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FFqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool FFqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator > 1 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() + partons[emitter]->id() == 0 &&
partons[emitter]->mass() == ZERO &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double FFqqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
double res = 1.-2.*z*(1.-z);
res *= -ccme2;
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FFqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]));
Lorentz5Momentum pc =
z*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
(1.-z)*realEmissionME()->lastXComb().meMomenta()[realEmission()];
SpinCorrelationTensor corr(-1.,pc,-prop/4.);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FFqqxDipole::persistentOutput(PersistentOStream &) const {
}
void FFqqxDipole::persistentInput(PersistentIStream &, int) {
}
void FFqqxDipole::Init() {
static ClassDocumentation<FFqqxDipole> documentation
("FFqqxDipole");
DipoleRepository::registerDipole<FFqqxDipole,FFLightTildeKinematics,FFLightInvertedTildeKinematics>
("FFqqxDipole","FFLightTildeKinematics","FFLightInvertedTildeKinematics");
}
// *** 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<FFqqxDipole,SubtractionDipole>
describeHerwigFFqqxDipole("Herwig::FFqqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FIMqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/FIMqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FIMqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FIMqgxDipole.cc
@@ -1,174 +1,154 @@
// -*- C++ -*-
//
// FIqgxDipole.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 FIMqgxDipole class.
//
#include "FIMqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h"
using namespace Herwig;
FIMqgxDipole::FIMqgxDipole()
: SubtractionDipole() {}
FIMqgxDipole::~FIMqgxDipole() {}
IBPtr FIMqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FIMqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool FIMqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator < 2 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 6 &&
!(partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
double FIMqgxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
// extra mass terms cancel
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
// NOTE: extra term taken from FIqgxDipole implementation
res *= ( 2./(1.-z+(1.-x)) -(1.+z) +(1.-x)*(1.+3.*x*z) -
sqr(realEmissionME()->lastXComb().mePartonData()[realEmission()]->mass()) / prop * 2.*x);
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FIMqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
// extra mass terms cancel
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
// NOTE: extra term taken from FIqgxDipole implementation
res *= ( 2./(1.-z+(1.-x)) -(1.+z) +(1.-x)*(1.+3.*x*z) -
sqr(realEmissionME()->lastXComb().mePartonData()[realEmission()]->mass()) / prop * 2.*x);
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FIMqgxDipole::persistentOutput(PersistentOStream &) const {
}
void FIMqgxDipole::persistentInput(PersistentIStream &, int) {
}
void FIMqgxDipole::Init() {
static ClassDocumentation<FIMqgxDipole> documentation
("FIMqgxDipole");
DipoleRepository::registerDipole<FIMqgxDipole,FILightTildeKinematics,FILightInvertedTildeKinematics>
("FIMqgxDipole","FILightTildeKinematics","FILightInvertedTildeKinematics");
}
// *** 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<FIMqgxDipole,SubtractionDipole>
describeHerwigFIMqgxDipole("Herwig::FIMqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FIMqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/FIMqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FIMqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FIMqqxDipole.cc
@@ -1,183 +1,163 @@
// -*- C++ -*-
//
// FIMqqxDipole.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 FIqqxDipole class.
//
#include "FIMqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h"
using namespace Herwig;
FIMqqxDipole::FIMqqxDipole()
: SubtractionDipole() {}
FIMqqxDipole::~FIMqqxDipole() {}
IBPtr FIMqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FIMqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool FIMqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator < 2 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() + partons[emitter]->id() == 0 &&
!(partons[emitter]->mass() == ZERO &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
double FIMqqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
Energy2 mQ2 = sqr(realEmissionME()->lastXComb().mePartonData()[realEmitter()]->mass());
double muQ2 = x * mQ2 /
((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realSpectator()]));
// mu_ij=0, mu_i=mu_j=mu_Q.
double zm = ( 1.-x - sqrt( sqr(1.-x-2.*muQ2) - 4.*muQ2 ) ) /
( 2.*(1.-x) );
double zp = ( 1.-x + sqrt( sqr(1.-x-2.*muQ2) - 4.*muQ2 ) ) /
( 2.*(1.-x) );
double res = 1.-2.*(z-zm)*(zp-z);
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/( prop+2.*mQ2*x );
+ (underlyingBornME()->lastXComb().lastAlphaS())/( prop+2.*mQ2*x );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FIMqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
Energy2 mQ2 = sqr((realEmissionME()->lastXComb().mePartonData()[realEmitter()])->mass());
Lorentz5Momentum pc =
z*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
(1.-z)*realEmissionME()->lastXComb().meMomenta()[realEmission()];
SpinCorrelationTensor corr(-1.,pc,-(prop+2.*mQ2*x)/(4.*x));
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/( prop+2.*mQ2*x );
+ (underlyingBornME()->lastXComb().lastAlphaS())/( prop+2.*mQ2*x );
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FIMqqxDipole::persistentOutput(PersistentOStream &) const {
}
void FIMqqxDipole::persistentInput(PersistentIStream &, int) {
}
void FIMqqxDipole::Init() {
static ClassDocumentation<FIMqqxDipole> documentation
("FIMqqxDipole");
DipoleRepository::registerDipole<FIMqqxDipole,FILightTildeKinematics,FILightInvertedTildeKinematics>
("FIMqqxDipole","FILightTildeKinematics","FILightInvertedTildeKinematics");
}
// *** 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<FIMqqxDipole,SubtractionDipole>
describeHerwigFIMqqxDipole("Herwig::FIMqqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FIggxDipole.cc b/MatrixElement/Matchbox/Dipoles/FIggxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FIggxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FIggxDipole.cc
@@ -1,169 +1,149 @@
// -*- C++ -*-
//
// FIggxDipole.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 FIggxDipole class.
//
#include "FIggxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h"
using namespace Herwig;
FIggxDipole::FIggxDipole()
: SubtractionDipole() {}
FIggxDipole::~FIggxDipole() {}
IBPtr FIggxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FIggxDipole::fullclone() const {
return new_ptr(*this);
}
bool FIggxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator < 2 &&
partons[emission]->id() == ParticleID::g &&
partons[emitter]->id() == ParticleID::g &&
partons[spectator]->mass() == ZERO;
}
double FIggxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res = 1./((1.-z)+(1.-x))+1./(z+(1.-x))-2.+z*(1.-z)+(1.-x)*(1.+x*z*(1.-z));
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FIggxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double diag = 1./(1.-z+1.-x)+1./(z+1.-x)-2.+(1.-x)*(1.+x*z*(1.-z));
Lorentz5Momentum pc =
z*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
(1.-z)*realEmissionME()->lastXComb().meMomenta()[realEmission()];
SpinCorrelationTensor corr(-diag,pc,prop/(2.*x));
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FIggxDipole::persistentOutput(PersistentOStream &) const {
}
void FIggxDipole::persistentInput(PersistentIStream &, int) {
}
void FIggxDipole::Init() {
static ClassDocumentation<FIggxDipole> documentation
("FIggxDipole");
DipoleRepository::registerDipole<FIggxDipole,FILightTildeKinematics,FILightInvertedTildeKinematics>
("FIggxDipole","FILightTildeKinematics","FILightInvertedTildeKinematics");
}
// *** 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<FIggxDipole,SubtractionDipole>
describeHerwigFIggxDipole("Herwig::FIggxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FIqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/FIqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FIqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FIqgxDipole.cc
@@ -1,168 +1,148 @@
// -*- C++ -*-
//
// FIqgxDipole.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 FIqgxDipole class.
//
#include "FIqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h"
using namespace Herwig;
FIqgxDipole::FIqgxDipole()
: SubtractionDipole() {}
FIqgxDipole::~FIqgxDipole() {}
IBPtr FIqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FIqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool FIqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator < 2 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 6 &&
partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double FIqgxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-z+(1.-x)) -(1.+z) +(1.-x)*(1.+3.*x*z) );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FIqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-z+(1.-x)) -(1.+z) +(1.-x)*(1.+3.*x*z) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FIqgxDipole::persistentOutput(PersistentOStream &) const {
}
void FIqgxDipole::persistentInput(PersistentIStream &, int) {
}
void FIqgxDipole::Init() {
static ClassDocumentation<FIqgxDipole> documentation
("FIqgxDipole");
DipoleRepository::registerDipole<FIqgxDipole,FILightTildeKinematics,FILightInvertedTildeKinematics>
("FIqgxDipole","FILightTildeKinematics","FILightInvertedTildeKinematics");
}
// *** 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<FIqgxDipole,SubtractionDipole>
describeHerwigFIqgxDipole("Herwig::FIqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/FIqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/FIqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/FIqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/FIqqxDipole.cc
@@ -1,170 +1,150 @@
// -*- C++ -*-
//
// FIqqxDipole.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 FIqqxDipole class.
//
#include "FIqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/FILightInvertedTildeKinematics.h"
using namespace Herwig;
FIqqxDipole::FIqqxDipole()
: SubtractionDipole() {}
FIqqxDipole::~FIqqxDipole() {}
IBPtr FIqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr FIqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool FIqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter > 1 && spectator < 2 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() + partons[emitter]->id() == 0 &&
partons[emitter]->mass() == ZERO &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double FIqqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res = 1.-2.*z*(1.-z);
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double FIqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double z = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
Lorentz5Momentum pc =
z*realEmissionME()->lastXComb().meMomenta()[realEmitter()] -
(1.-z)*realEmissionME()->lastXComb().meMomenta()[realEmission()];
SpinCorrelationTensor corr(-1.,pc,-prop/(4.*x));
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 4.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void FIqqxDipole::persistentOutput(PersistentOStream &) const {
}
void FIqqxDipole::persistentInput(PersistentIStream &, int) {
}
void FIqqxDipole::Init() {
static ClassDocumentation<FIqqxDipole> documentation
("FIqqxDipole");
DipoleRepository::registerDipole<FIqqxDipole,FILightTildeKinematics,FILightInvertedTildeKinematics>
("FIqqxDipole","FILightTildeKinematics","FILightInvertedTildeKinematics");
}
// *** 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<FIqqxDipole,SubtractionDipole>
describeHerwigFIqqxDipole("Herwig::FIqqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFMggxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFMggxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFMggxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFMggxDipole.cc
@@ -1,178 +1,158 @@
// -*- C++ -*-
//
// IFMggxDipole.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 IFggxDipole class.
//
#include "IFMggxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFMggxDipole::IFMggxDipole()
: SubtractionDipole() {}
IFMggxDipole::~IFMggxDipole() {}
IBPtr IFMggxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFMggxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFMggxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
partons[emitter]->id() == ParticleID::g &&
partons[spectator]->mass() != ZERO;
}
double IFMggxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double muj2 = sqr( (realEmissionME()->lastXComb().mePartonData()[realSpectator()]->mass()) ) /
(2.* (realEmissionME()->lastXComb().meMomenta()[bornSpectator()])*
(realEmissionME()->lastXComb().meMomenta()[realEmitter()]) );
double res = 1./(1.-x+u) + (1.-x)/x - 1. + x*(1.-x) -
muj2/x*u/(1.-u);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFMggxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double diag = 1./(1.-x+u)-1.+x*(1.-x);
Lorentz5Momentum pc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]/u -
realEmissionME()->lastXComb().meMomenta()[realSpectator()]/(1.-u);
Energy2 sc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]*
realEmissionME()->lastXComb().meMomenta()[realSpectator()];
sc /= u*(1.-u)*(1.-x)/x;
SpinCorrelationTensor corr(-diag,pc,sc);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFMggxDipole::persistentOutput(PersistentOStream &) const {
}
void IFMggxDipole::persistentInput(PersistentIStream &, int) {
}
void IFMggxDipole::Init() {
static ClassDocumentation<IFMggxDipole> documentation
("IFMggxDipole");
DipoleRepository::registerDipole<IFMggxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFMggxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFMggxDipole,SubtractionDipole>
describeHerwigIFMggxDipole("Herwig::IFMggxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFMgqxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFMgqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFMgqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFMgqxDipole.cc
@@ -1,162 +1,142 @@
// -*- C++ -*-
//
// IFMgqxDipole.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 IFgqxDipole class.
//
#include "IFMgqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFMgqxDipole::IFMgqxDipole()
: SubtractionDipole() {}
IFMgqxDipole::~IFMgqxDipole() {}
IBPtr IFMgqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFMgqxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFMgqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
partons[emitter]->id() == ParticleID::g &&
abs(partons[emission]->id()) < 6 &&
!(partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
double IFMgqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res =
8.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= .5 * ( 1.-2.*x*(1.-x) );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFMgqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res =
8.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= .5 * ( 1.-2.*x*(1.-x) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFMgqxDipole::persistentOutput(PersistentOStream &) const {
}
void IFMgqxDipole::persistentInput(PersistentIStream &, int) {
}
void IFMgqxDipole::Init() {
static ClassDocumentation<IFMgqxDipole> documentation
("IFMgqxDipole");
DipoleRepository::registerDipole<IFMgqxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFMgqxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFMgqxDipole,SubtractionDipole>
describeHerwigIFMgqxDipole("Herwig::IFMgqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFMqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFMqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFMqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFMqgxDipole.cc
@@ -1,170 +1,150 @@
// -*- C++ -*-
//
// IFMqgxDipole.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 IFqgxDipole class.
//
#include "IFMqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFMqgxDipole::IFMqgxDipole()
: SubtractionDipole() {}
IFMqgxDipole::~IFMqgxDipole() {}
IBPtr IFMqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFMqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFMqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 6 &&
!(partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
double IFMqgxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
// NOTE: extra term same as in IFqgxDipole
res *= ( 2./(1.-x+u) - (1.+x) + u*(1.+3.*x*(1.-u)) );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFMqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
// NOTE: extra term same as in IFqgxDipole
res *= ( 2./(1.-x+u) - (1.+x) + u*(1.+3.*x*(1.-u)) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFMqgxDipole::persistentOutput(PersistentOStream &) const {
}
void IFMqgxDipole::persistentInput(PersistentIStream &, int) {
}
void IFMqgxDipole::Init() {
static ClassDocumentation<IFMqgxDipole> documentation
("IFMqgxDipole");
DipoleRepository::registerDipole<IFMqgxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFMqgxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFMqgxDipole,SubtractionDipole>
describeHerwigIFMqgxDipole("Herwig::IFMqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFMqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFMqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFMqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFMqqxDipole.cc
@@ -1,184 +1,164 @@
// -*- C++ -*-
//
// IFMqqxDipole.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 IFqqxDipole class.
//
#include "IFMqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFMqqxDipole::IFMqqxDipole()
: SubtractionDipole() {}
IFMqqxDipole::~IFMqqxDipole() {}
IBPtr IFMqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFMqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFMqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() - partons[emitter]->id() == 0 &&
!(partons[emitter]->mass() == ZERO &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO);
}
double IFMqqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double muj2 = sqr( (realEmissionME()->lastXComb().mePartonData()[realSpectator()]->mass()) ) /
(2.* (realEmissionME()->lastXComb().meMomenta()[bornSpectator()])*
(realEmissionME()->lastXComb().meMomenta()[realEmitter()]) );
double res = x + 2.*(1.-x)/x -
2.*muj2/x*u/(1.-u);
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
res *= 8.*CF*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFMqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
Lorentz5Momentum pc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]/u -
realEmissionME()->lastXComb().meMomenta()[realSpectator()]/(1.-u);
Energy2 sc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]*
realEmissionME()->lastXComb().meMomenta()[realSpectator()];
sc /= u*(1.-u)*(1.-x)/x;
SpinCorrelationTensor corr(-x,pc,sc/2.);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
res *= 8.*CF*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFMqqxDipole::persistentOutput(PersistentOStream &) const {
}
void IFMqqxDipole::persistentInput(PersistentIStream &, int) {
}
void IFMqqxDipole::Init() {
static ClassDocumentation<IFMqqxDipole> documentation
("IFMqqxDipole");
DipoleRepository::registerDipole<IFMqqxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFMqqxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFMqqxDipole,SubtractionDipole>
describeHerwigIFMqqxDipole("Herwig::IFMqqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFggxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFggxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFggxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFggxDipole.cc
@@ -1,173 +1,153 @@
// -*- C++ -*-
//
// IFggxDipole.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 IFggxDipole class.
//
#include "IFggxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFggxDipole::IFggxDipole()
: SubtractionDipole() {}
IFggxDipole::~IFggxDipole() {}
IBPtr IFggxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFggxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFggxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
partons[emitter]->id() == ParticleID::g &&
partons[spectator]->mass() == ZERO;
}
double IFggxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res = 1./(1.-x+u) + (1.-x)/x - 1. + x*(1.-x);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFggxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double diag = 1./(1.-x+u)-1.+x*(1.-x);
Lorentz5Momentum pc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]/u -
realEmissionME()->lastXComb().meMomenta()[realSpectator()]/(1.-u);
Energy2 sc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]*
realEmissionME()->lastXComb().meMomenta()[realSpectator()];
sc /= u*(1.-u)*(1.-x)/x;
SpinCorrelationTensor corr(-diag,pc,sc);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFggxDipole::persistentOutput(PersistentOStream &) const {
}
void IFggxDipole::persistentInput(PersistentIStream &, int) {
}
void IFggxDipole::Init() {
static ClassDocumentation<IFggxDipole> documentation
("IFggxDipole");
DipoleRepository::registerDipole<IFggxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFggxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFggxDipole,SubtractionDipole>
describeHerwigIFggxDipole("Herwig::IFggxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFgqxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFgqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFgqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFgqxDipole.cc
@@ -1,162 +1,142 @@
// -*- C++ -*-
//
// IFgqxDipole.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 IFgqxDipole class.
//
#include "IFgqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFgqxDipole::IFgqxDipole()
: SubtractionDipole() {}
IFgqxDipole::~IFgqxDipole() {}
IBPtr IFgqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFgqxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFgqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
partons[emitter]->id() == ParticleID::g &&
abs(partons[emission]->id()) < 6 &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double IFgqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res =
8.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= .5 * ( 1.-2.*x*(1.-x) );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFgqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res =
8.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= .5 * ( 1.-2.*x*(1.-x) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFgqxDipole::persistentOutput(PersistentOStream &) const {
}
void IFgqxDipole::persistentInput(PersistentIStream &, int) {
}
void IFgqxDipole::Init() {
static ClassDocumentation<IFgqxDipole> documentation
("IFgqxDipole");
DipoleRepository::registerDipole<IFgqxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFgqxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFgqxDipole,SubtractionDipole>
describeHerwigIFgqxDipole("Herwig::IFgqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFqgxDipole.cc
@@ -1,168 +1,148 @@
// -*- C++ -*-
//
// IFqgxDipole.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 IFqgxDipole class.
//
#include "IFqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFqgxDipole::IFqgxDipole()
: SubtractionDipole() {}
IFqgxDipole::~IFqgxDipole() {}
IBPtr IFqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 6 &&
partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double IFqgxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-x+u) - (1.+x) + u*(1.+3.*x*(1.-u)) );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= ( 2./(1.-x+u) - (1.+x) + u*(1.+3.*x*(1.-u)) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFqgxDipole::persistentOutput(PersistentOStream &) const {
}
void IFqgxDipole::persistentInput(PersistentIStream &, int) {
}
void IFqgxDipole::Init() {
static ClassDocumentation<IFqgxDipole> documentation
("IFqgxDipole");
DipoleRepository::registerDipole<IFqgxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFqgxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFqgxDipole,SubtractionDipole>
describeHerwigIFqgxDipole("Herwig::IFqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IFqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/IFqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IFqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IFqqxDipole.cc
@@ -1,178 +1,158 @@
// -*- C++ -*-
//
// IFqqxDipole.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 IFqqxDipole class.
//
#include "IFqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IFLightInvertedTildeKinematics.h"
using namespace Herwig;
IFqqxDipole::IFqqxDipole()
: SubtractionDipole() {}
IFqqxDipole::~IFqqxDipole() {}
IBPtr IFqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IFqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool IFqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator > 1 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() - partons[emitter]->id() == 0 &&
partons[emitter]->mass() == ZERO &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double IFqqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res = x + 2.*(1.-x)/x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
res *= 8.*CF*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IFqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double u = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
Lorentz5Momentum pc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]/u -
realEmissionME()->lastXComb().meMomenta()[realSpectator()]/(1.-u);
Energy2 sc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]*
realEmissionME()->lastXComb().meMomenta()[realSpectator()];
sc /= u*(1.-u)*(1.-x)/x;
SpinCorrelationTensor corr(-x,pc,sc/2.);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
res *= 8.*CF*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IFqqxDipole::persistentOutput(PersistentOStream &) const {
}
void IFqqxDipole::persistentInput(PersistentIStream &, int) {
}
void IFqqxDipole::Init() {
static ClassDocumentation<IFqqxDipole> documentation
("IFqqxDipole");
DipoleRepository::registerDipole<IFqqxDipole,IFLightTildeKinematics,IFLightInvertedTildeKinematics>
("IFqqxDipole","IFLightTildeKinematics","IFLightInvertedTildeKinematics");
}
// *** 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<IFqqxDipole,SubtractionDipole>
describeHerwigIFqqxDipole("Herwig::IFqqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IIggxDipole.cc b/MatrixElement/Matchbox/Dipoles/IIggxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IIggxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IIggxDipole.cc
@@ -1,172 +1,152 @@
// -*- C++ -*-
//
// IIggxDipole.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 IIggxDipole class.
//
#include "IIggxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.h"
using namespace Herwig;
IIggxDipole::IIggxDipole()
: SubtractionDipole() {}
IIggxDipole::~IIggxDipole() {}
IBPtr IIggxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IIggxDipole::fullclone() const {
return new_ptr(*this);
}
bool IIggxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator < 2 &&
partons[emission]->id() == ParticleID::g &&
partons[emitter]->id() == ParticleID::g &&
partons[spectator]->mass() == ZERO;
}
double IIggxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res = x/(1.-x) + (1.-x)/x + x*(1.-x);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IIggxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double v = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double diag = x/(1.-x)+x*(1.-x);
Lorentz5Momentum pc =
realEmissionME()->lastXComb().meMomenta()[realEmission()] -
v*realEmissionME()->lastXComb().meMomenta()[realSpectator()];
Energy2 sc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]*
realEmissionME()->lastXComb().meMomenta()[realSpectator()];
sc /= (1.-x)/(x*v);
SpinCorrelationTensor corr(-diag,pc,sc);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
res *= 16.*Constants::pi*SM().Nc()*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IIggxDipole::persistentOutput(PersistentOStream &) const {
}
void IIggxDipole::persistentInput(PersistentIStream &, int) {
}
void IIggxDipole::Init() {
static ClassDocumentation<IIggxDipole> documentation
("IIggxDipole");
DipoleRepository::registerDipole<IIggxDipole,IILightTildeKinematics,IILightInvertedTildeKinematics>
("IIggxDipole","IILightTildeKinematics","IILightInvertedTildeKinematics");
}
// *** 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<IIggxDipole,SubtractionDipole>
describeHerwigIIggxDipole("Herwig::IIggxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IIgqxDipole.cc b/MatrixElement/Matchbox/Dipoles/IIgqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IIgqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IIgqxDipole.cc
@@ -1,162 +1,142 @@
// -*- C++ -*-
//
// IIgqxDipole.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 IIgqxDipole class.
//
#include "IIgqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.h"
using namespace Herwig;
IIgqxDipole::IIgqxDipole()
: SubtractionDipole() {}
IIgqxDipole::~IIgqxDipole() {}
IBPtr IIgqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IIgqxDipole::fullclone() const {
return new_ptr(*this);
}
bool IIgqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator < 2 &&
partons[emitter]->id() == ParticleID::g &&
abs(partons[emission]->id()) < 6 &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double IIgqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res =
8.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= .5 * ( 1.-2.*x*(1.-x) );
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IIgqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res =
8.*Constants::pi*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= .5 * ( 1.-2.*x*(1.-x) );
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IIgqxDipole::persistentOutput(PersistentOStream &) const {
}
void IIgqxDipole::persistentInput(PersistentIStream &, int) {
}
void IIgqxDipole::Init() {
static ClassDocumentation<IIgqxDipole> documentation
("IIgqxDipole");
DipoleRepository::registerDipole<IIgqxDipole,IILightTildeKinematics,IILightInvertedTildeKinematics>
("IIgqxDipole","IILightTildeKinematics","IILightInvertedTildeKinematics");
}
// *** 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<IIgqxDipole,SubtractionDipole>
describeHerwigIIgqxDipole("Herwig::IIgqxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IIqgxDipole.cc b/MatrixElement/Matchbox/Dipoles/IIqgxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IIqgxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IIqgxDipole.cc
@@ -1,166 +1,146 @@
// -*- C++ -*-
//
// IIqgxDipole.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 IIqgxDipole class.
//
#include "IIqgxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.h"
using namespace Herwig;
IIqgxDipole::IIqgxDipole()
: SubtractionDipole() {}
IIqgxDipole::~IIqgxDipole() {}
IBPtr IIqgxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IIqgxDipole::fullclone() const {
return new_ptr(*this);
}
bool IIqgxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator < 2 &&
partons[emission]->id() == ParticleID::g &&
abs(partons[emitter]->id()) < 6 &&
partons[emitter]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double IIqgxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= 2./(1.-x) - (1.+x);
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IIqgxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
double res =
8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= 2./(1.-x) - (1.+x);
res *= -underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IIqgxDipole::persistentOutput(PersistentOStream &) const {
}
void IIqgxDipole::persistentInput(PersistentIStream &, int) {
}
void IIqgxDipole::Init() {
static ClassDocumentation<IIqgxDipole> documentation
("IIqgxDipole");
DipoleRepository::registerDipole<IIqgxDipole,IILightTildeKinematics,IILightInvertedTildeKinematics>
("IIqgxDipole","IILightTildeKinematics","IILightInvertedTildeKinematics");
}
// *** 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<IIqgxDipole,SubtractionDipole>
describeHerwigIIqgxDipole("Herwig::IIqgxDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/IIqqxDipole.cc b/MatrixElement/Matchbox/Dipoles/IIqqxDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/IIqqxDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/IIqqxDipole.cc
@@ -1,178 +1,158 @@
// -*- C++ -*-
//
// IIqqxDipole.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 IIqqxDipole class.
//
#include "IIqqxDipole.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.h"
using namespace Herwig;
IIqqxDipole::IIqqxDipole()
: SubtractionDipole() {}
IIqqxDipole::~IIqqxDipole() {}
IBPtr IIqqxDipole::clone() const {
return new_ptr(*this);
}
IBPtr IIqqxDipole::fullclone() const {
return new_ptr(*this);
}
bool IIqqxDipole::canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const {
return
emitter < 2 && spectator < 2 &&
abs(partons[emission]->id()) < 6 &&
abs(partons[emitter]->id()) < 6 &&
partons[emission]->id() - partons[emitter]->id() == 0 &&
partons[emitter]->mass() == ZERO &&
partons[emission]->mass() == ZERO &&
partons[spectator]->mass() == ZERO;
}
double IIqqxDipole::me2Avg(double ccme2) const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
double res = (1.+sqr(1.-x))/x;
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
res *= 8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *= -ccme2;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
return res;
}
double IIqqxDipole::me2() const {
if ( jacobian() == 0.0 )
return 0.0;
double x = subtractionParameters()[0];
double v = subtractionParameters()[1];
Energy2 prop =
2.*((realEmissionME()->lastXComb().meMomenta()[realEmitter()])*
(realEmissionME()->lastXComb().meMomenta()[realEmission()]))*x;
Lorentz5Momentum pc =
realEmissionME()->lastXComb().meMomenta()[realEmission()] -
v*realEmissionME()->lastXComb().meMomenta()[realSpectator()];
Energy2 sc =
realEmissionME()->lastXComb().meMomenta()[realEmission()]*
realEmissionME()->lastXComb().meMomenta()[realSpectator()];
sc /= (1.-x)/(x*v);
SpinCorrelationTensor corr(-x,pc,sc/2.);
double res = -underlyingBornME()->spinColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),
corr);
double CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
res *= 8.*Constants::pi*CF*(realEmissionME()->lastXComb().lastSHat())*
- (realEmissionME()->lastXComb().lastAlphaS())/prop;
+ (underlyingBornME()->lastXComb().lastAlphaS())/prop;
res *=
pow(realEmissionME()->lastXComb().lastSHat() / underlyingBornME()->lastXComb().lastSHat(),
underlyingBornME()->lastXComb().mePartonData().size()-4.);
res *=
realEmissionME()->finalStateSymmetry() /
underlyingBornME()->finalStateSymmetry();
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaS()/
- underlyingBornME()->lastXComb().lastAlphaS(),
- underlyingBornME()->orderInAlphaS());
-
- res *=
- pow(realEmissionME()->lastXComb().lastAlphaEM()/
- underlyingBornME()->lastXComb().lastAlphaEM(),
- underlyingBornME()->orderInAlphaEW());
-
lastME2(res);
logME2();
return res;
}
void IIqqxDipole::persistentOutput(PersistentOStream &) const {
}
void IIqqxDipole::persistentInput(PersistentIStream &, int) {
}
void IIqqxDipole::Init() {
static ClassDocumentation<IIqqxDipole> documentation
("IIqqxDipole");
DipoleRepository::registerDipole<IIqqxDipole,IILightTildeKinematics,IILightInvertedTildeKinematics>
("IIqqxDipole","IILightTildeKinematics","IILightInvertedTildeKinematics");
}
// *** 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<IIqqxDipole,SubtractionDipole>
describeHerwigIIqqxDipole("Herwig::IIqqxDipole", "HwMatchbox.so");
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,1043 +1,1068 @@
// -*- 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 <iterator>
using std::ostream_iterator;
using namespace Herwig;
SubtractionDipole::SubtractionDipole()
: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
theIgnoreCuts(false), theShowerKernel(false),
theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
theBornEmitter(-1), theBornSpectator(-1),
- theRealShowerSubtraction(false), theVirtualShowerSubtraction(false) {}
+ theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
+ theRealEmissionScales(false) {}
SubtractionDipole::~SubtractionDipole() {}
void SubtractionDipole::setupBookkeeping() {
theMergingMap.clear();
theSplittingMap.clear();
int xemitter = -1;
int xspectator = -1;
map<int,int> mergeLegs;
map<int,int> remapLegs;
map<int,int> realBornMap;
map<int,int> bornRealMap;
for ( DiagramVector::const_iterator d =
theRealEmissionME->diagrams().begin();
d != theRealEmissionME->diagrams().end(); ++d ) {
Ptr<Tree2toNDiagram>::ptr check =
new_ptr(*dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d));
map<int,int> theMergeLegs;
map<int,int> theRemapLegs;
for ( unsigned int i = 0; i < check->external().size(); ++i )
theMergeLegs[i] = -1;
int theEmitter = check->mergeEmission(realEmitter(),realEmission(),theMergeLegs);
// no underlying Born
if ( theEmitter == -1 )
continue;
DiagramVector::const_iterator bd =
theUnderlyingBornME->diagrams().end();
for ( DiagramVector::const_iterator b =
theUnderlyingBornME->diagrams().begin();
b != theUnderlyingBornME->diagrams().end(); ++b )
if ( check->isSame(*b,theRemapLegs) ) {
bd = b; break;
}
// no underlying Born
if ( bd == theUnderlyingBornME->diagrams().end() )
continue;
if ( xemitter == -1 ) {
xemitter = theEmitter;
mergeLegs = theMergeLegs;
remapLegs = theRemapLegs;
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((**d).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)));
theSplittingMap.insert(make_pair(bornKey,make_pair(realKey,bornRealMap)));
}
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));
}
theUnderlyingBornDiagrams.clear();
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator r =
theMergingMap.begin(); r != theMergingMap.end(); ++r ) {
DiagramVector borns;
for ( DiagramVector::const_iterator d = theUnderlyingBornME->diagrams().begin();
d != theUnderlyingBornME->diagrams().end(); ++d ) {
if ( (**d).partons() == process(r->second.first) )
borns.push_back(*d);
}
theUnderlyingBornDiagrams[process(r->first)] = borns;
}
theRealEmissionDiagrams.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator r =
theSplittingMap.begin(); r != theSplittingMap.end(); ++r ) {
DiagramVector reals;
for ( DiagramVector::const_iterator d = theRealEmissionME->diagrams().begin();
d != theRealEmissionME->diagrams().end(); ++d ) {
if ( (**d).partons() == process(r->second.first) )
reals.push_back(*d);
}
theRealEmissionDiagrams[process(r->first)] = reals;
}
}
void SubtractionDipole::subtractionBookkeeping() {
if ( theMergingMap.empty() )
setupBookkeeping();
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();
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::makeBornXComb(tStdXCombPtr realXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();
if ( theMergingMap.find(realEmissionKey(proc,realEmitter(),
realEmission(),realSpectator())) ==
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());
StdXCombPtr res =
new_ptr(StandardXComb(realXC,*ppit,this,bornDiags));
return res;
}
vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();
if ( theSplittingMap.find(underlyingBornKey(proc,bornEmitter(),bornSpectator())) ==
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 =
new_ptr(StandardXComb(bornXC,*ppit,this,pr->second));
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;
}
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;
}
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()->prepare(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagrams(dv);
}
MEBase::DiagramIndex SubtractionDipole::diagram(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->prepare(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagram(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;
}
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;
realEmissionME()->lastXCombPtr()->setIncomingPartons();
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;
underlyingBornME()->setScale();
assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
underlyingBornME()->lastXCombPtr()->setIncomingPartons();
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());
if ( !tildeKinematics() ) {
jacobian(0.0);
return false;
}
theTildeKinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());
if ( !theTildeKinematics->doMap() ) {
jacobian(0.0);
return false;
}
theLastSubtractionScale = theTildeKinematics->lastScale();
theLastSubtractionPt = theTildeKinematics->lastPt();
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 ( theTildeKinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = theTildeKinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[bornEmitter()] = tildeKinematics()->bornEmitterMomentum();
meMomenta()[bornSpectator()] = tildeKinematics()->bornSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).mass());
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());
if ( !invertedTildeKinematics() ) {
jacobian(0.0);
return false;
}
theInvertedTildeKinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());
if ( !theInvertedTildeKinematics->doMap(r) ) {
jacobian(0.0);
return false;
}
theLastSplittingScale = theInvertedTildeKinematics->lastScale();
theLastSplittingPt = theInvertedTildeKinematics->lastPt();
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 ( invertedTildeKinematics()->doesTransform() && k > 1 )
meMomenta()[trans[k]] = invertedTildeKinematics()->transform(meMomenta()[trans[k]]);
}
meMomenta()[realEmitter()] = invertedTildeKinematics()->realEmitterMomentum();
meMomenta()[realEmission()] = invertedTildeKinematics()->realEmissionMomentum();
meMomenta()[realSpectator()] = invertedTildeKinematics()->realSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).mass());
p->rescaleRho();
}
jacobian(underlyingBornME()->lastXComb().jacobian() *
invertedTildeKinematics()->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);
lastME2(0.0);
return ZERO;
}
- if ( havePDFWeight1() || havePDFWeight2() ) {
- if ( splitting() )
- realEmissionME()->getPDFWeight(factorizationScale);
- pdfweight = realEmissionME()->lastXComb().lastMEPDFWeight();
- lastMEPDFWeight(pdfweight);
+ if ( factorizationScale == ZERO ) {
+ if ( realEmissionScales() )
+ factorizationScale = realEmissionME()->lastScale();
+ else
+ factorizationScale = underlyingBornME()->lastScale();
}
+ if ( havePDFWeight1() ) {
+ pdfweight *= realEmissionME()->pdf1(factorizationScale);
+ }
+
+ if ( havePDFWeight2() ) {
+ pdfweight *= realEmissionME()->pdf2(factorizationScale);
+ }
+
+ lastMEPDFWeight(pdfweight);
+
if ( showerApproximation() && realShowerSubtraction() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
if ( !showerApproximation()->isAboveCutoff() )
showerApproximation()->wasBelowCutoff();
if ( !showerApproximation()->isInShowerPhasespace() ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
lastMECrossSection(-showerApproximation()->dSigHatDR());
return lastMECrossSection();
}
double xme2 = 0.0;
if ( lastME2() == 0.0 ) {
if ( !showerKernel() )
xme2 = me2();
else
xme2 = me2Avg(-underlyingBornME()->me2());
} else {
xme2 = lastME2();
}
if ( xme2 == 0.0 ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
+ if ( realEmissionScales() ) {
+
+ xme2 *=
+ pow(realEmissionME()->lastXComb().lastAlphaS()/
+ underlyingBornME()->lastXComb().lastAlphaS(),
+ realEmissionME()->orderInAlphaS());
+
+ xme2 *=
+ pow(realEmissionME()->lastXComb().lastAlphaEM()/
+ underlyingBornME()->lastXComb().lastAlphaEM(),
+ underlyingBornME()->orderInAlphaEW());
+
+ }
+
lastME2(xme2);
CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
if ( !showerApproximation() ) {
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;
}
if ( showerApproximation() && virtualShowerSubtraction() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
if ( !showerApproximation()->isAboveCutoff() )
showerApproximation()->wasBelowCutoff();
CrossSection shower = ZERO;
if ( showerApproximation()->isInShowerPhasespace() )
shower = showerApproximation()->dSigHatDR();
res -= shower;
}
lastMECrossSection(-res);
logDSigHatDR(jac);
return lastMECrossSection();
}
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";
generator()->log() << "me2 = " << lastME2() << "\n" << flush;
}
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::cloneDependencies(const std::string& prefix) {
if ( underlyingBornME() ) {
Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myUnderlyingBornME->cloneDependencies(pname.str());
underlyingBornME(myUnderlyingBornME);
}
if ( realEmissionME() ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = realEmissionME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
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 InitException() << "Tilde kinematics " << pname.str() << " already existing.";
myTildeKinematics->dipole(this);
tildeKinematics(myTildeKinematics);
}
if ( invertedTildeKinematics() ) {
Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
throw InitException() << "Inverted tilde kinematics " << pname.str() << " already existing.";
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 InitException() << "Reweight " << pname.str() << " already existing.";
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
}
void SubtractionDipole::constructVertex(tSubProPtr sub) {
if ( splitting() )
realEmissionME()->constructVertex(sub);
else
underlyingBornME()->constructVertex(sub);
}
void SubtractionDipole::generateSubCollision(SubProcess & sub) {
if ( splitting() )
realEmissionME()->generateSubCollision(sub);
else
underlyingBornME()->generateSubCollision(sub);
}
void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
os << theSplitting << theApply << theSubtractionTest << theIgnoreCuts
<< theShowerKernel << theRealEmissionME << theUnderlyingBornME
<< theTildeKinematics << theInvertedTildeKinematics << theReweights
<< theRealEmitter << theRealEmission << theRealSpectator
<< theSubtractionParameters
<< theMergingMap << theSplittingMap << theIndexMap
<< theUnderlyingBornDiagrams << theRealEmissionDiagrams
<< lastRealEmissionKey << lastUnderlyingBornKey
<< theBornEmitter << theBornSpectator
<< theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction
- << thePartners;
+ << thePartners << theRealEmissionScales;
}
void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
is >> theSplitting >> theApply >> theSubtractionTest >> theIgnoreCuts
>> theShowerKernel >> theRealEmissionME >> theUnderlyingBornME
>> theTildeKinematics >> theInvertedTildeKinematics >> theReweights
>> theRealEmitter >> theRealEmission >> theRealSpectator
>> theSubtractionParameters
>> theMergingMap >> theSplittingMap >> theIndexMap
>> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
>> lastRealEmissionKey >> lastUnderlyingBornKey
>> theBornEmitter >> theBornSpectator
>> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction
- >> thePartners;
+ >> thePartners >> theRealEmissionScales;
}
void SubtractionDipole::Init() {
static ClassDocumentation<SubtractionDipole> documentation
("SubtractionDipole represents a dipole subtraction "
"term in the formalism of Catani and Seymour.");
static Reference<SubtractionDipole,MatchboxMEBase> interfaceUnderlyingBornME
("UnderlyingBornME",
"The underlying Born matrix element.",
&SubtractionDipole::theUnderlyingBornME, false, false, true, true, false);
static Reference<SubtractionDipole,MatchboxMEBase> interfaceRealEmissionME
("RealEmissionME",
"The real emission matrix element.",
&SubtractionDipole::theRealEmissionME, false, false, true, true, false);
static RefVector<SubtractionDipole,SubtractionDipole> interfacePartners
("Partners",
"The partner dipoles found along with this dipole.",
&SubtractionDipole::thePartners, -1, false, false, true, false, false);
static Reference<SubtractionDipole,TildeKinematics> interfaceTildeKinematics
("TildeKinematics",
"Set the TildeKinematics object to be used.",
&SubtractionDipole::theTildeKinematics, false, false, true, false, false);
static Reference<SubtractionDipole,InvertedTildeKinematics> interfaceInvertedTildeKinematics
("InvertedTildeKinematics",
"Set the InvertedTildeKinematics object to be used.",
&SubtractionDipole::theInvertedTildeKinematics, false, false, true, true, false);
static RefVector<SubtractionDipole,MatchboxReweightBase> interfaceReweights
("Reweights",
"Reweight objects to be applied to this matrix element.",
&SubtractionDipole::theReweights, -1, false, false, true, true, false);
static Reference<SubtractionDipole,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&SubtractionDipole::theShowerApproximation, false, false, true, true, false);
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<SubtractionDipole,MEBase>
describeSubtractionDipole("Herwig::SubtractionDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
@@ -1,991 +1,1018 @@
// -*- C++ -*-
//
// SubtractionDipole.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_SubtractionDipole_H
#define HERWIG_SubtractionDipole_H
//
// This is the declaration of the SubtractionDipole class.
//
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
namespace Herwig {
using namespace ThePEG;
class TildeKinematics;
class InvertedTildeKinematics;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief SubtractionDipole represents a dipole subtraction
* term in the formalism of Catani and Seymour.
*
*/
class SubtractionDipole: public MEBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SubtractionDipole();
/**
* The destructor.
*/
virtual ~SubtractionDipole();
//@}
public:
/** @name Subprocess and diagram information. */
//@{
/**
* Return true, if this dipole applies to the selected
* configuration.
*/
virtual bool canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const = 0;
/**
* Return true, if this dipole is symmetric with respect to emitter
* and emission.
*/
virtual bool isSymmetric() const { return false; }
/**
* If this is a dependent matrix element in a ME group, return true,
* if it applies to the process set in lastXComb()
*/
virtual bool apply() const { return theApply; }
/**
* Setup bookkeeping maps.
*/
void setupBookkeeping();
/**
* Get bookkeeping information for the given
* real emission diagram
*/
void subtractionBookkeeping();
/**
* Determine bookkeeping information for
* the underlying Born process supplied through
* the lastHeadXComb() object.
*/
void splittingBookkeeping();
/**
* Create a dependent xcomb object for the underlying
* Born process, given a XComb driving the real emission
*/
StdXCombPtr makeBornXComb(tStdXCombPtr realXC);
/**
* Create dependent xcomb objects for the real emission process,
* given a XComb driving the underlying Born
*/
vector<StdXCombPtr> makeRealXCombs(tStdXCombPtr bornXC);
/**
* Return true, if bookkeeping did not find a non-trivial setup.
*/
bool empty() const { return theSplittingMap.empty(); }
/**
* Return the emitter as referred to by the real emission
* matrix element.
*/
int realEmitter() const { return theRealEmitter; }
/**
* Set the emitter as referred to by the real emission
* matrix element.
*/
void realEmitter(int id) { theRealEmitter = id; }
/**
* Return the emission as referred to by the real emission
* matrix element.
*/
int realEmission() const { return theRealEmission; }
/**
* Set the emission as referred to by the real emission
* matrix element.
*/
void realEmission(int id) { theRealEmission = id; }
/**
* Return the spectator as referred to by the real emission
* matrix element.
*/
int realSpectator() const { return theRealSpectator; }
/**
* Set the spectator as referred to by the real emission
* matrix element.
*/
void realSpectator(int id) { theRealSpectator = id; }
/**
* Return the emitter as referred to by the underlying
* Born process.
*/
int bornEmitter() const { return theBornEmitter; }
/**
* Set the emitter as referred to by the underlying
* Born process.
*/
void bornEmitter(int id) { theBornEmitter = id; }
/**
* Return the spectator as referred to by the underlying
* Born process.
*/
int bornSpectator() const { return theBornSpectator; }
/**
* Set the spectator as referred to by the underlying
* Born process.
*/
void bornSpectator(int id) { theBornSpectator = id; }
/**
* Define the real emission key type
*/
typedef pair<pair<cPDVector,int>,pair<int,int> > RealEmissionKey;
/**
* Create a real emission key
*/
static RealEmissionKey realEmissionKey(const cPDVector& proc,
int em, int emm, int sp) {
return make_pair(make_pair(proc,emm),make_pair(em,sp));
}
/**
* Return the diagram of a real emission key
*/
static const cPDVector& process(const RealEmissionKey& key) {
return key.first.first;
}
/**
* Return the emission id of a real emission key
*/
static int emission(const RealEmissionKey& key) {
return key.first.second;
}
/**
* Return the emitter id of a real emission key
*/
static int emitter(const RealEmissionKey& key) {
return key.second.first;
}
/**
* Return the spectator id of a real emission key
*/
static int spectator(const RealEmissionKey& key) {
return key.second.second;
}
/**
* Define the underlying Born key type
*/
typedef pair<cPDVector,pair<int,int> > UnderlyingBornKey;
/**
* Create a underlying Born key
*/
static UnderlyingBornKey underlyingBornKey(const cPDVector& proc,
int em, int sp) {
return make_pair(proc,make_pair(em,sp));
}
/**
* Return the diagram of a underlying Born key
*/
static const cPDVector& process(const UnderlyingBornKey& key) {
return key.first;
}
/**
* Return the emitter id of a underlying Born key
*/
static int emitter(const UnderlyingBornKey& key) {
return key.second.first;
}
/**
* Return the spectator id of a underlying Born key
*/
static int spectator(const UnderlyingBornKey& key) {
return key.second.second;
}
/**
* Define real emission key and index dictionary
* for partons not involved in the given dipole.
*/
typedef pair<RealEmissionKey,map<int,int> > RealEmissionInfo;
/**
* Define underlying Born key and index dictionary
* for partons not involved in the given dipole.
*/
typedef pair<UnderlyingBornKey,map<int,int> > UnderlyingBornInfo;
/**
* Return the merging map
*/
const map<RealEmissionKey,UnderlyingBornInfo>& mergingMap() const { return theMergingMap; }
/**
* Return the splitting map
*/
const multimap<UnderlyingBornKey,RealEmissionInfo>& splittingMap() const { return theSplittingMap; }
/**
* Return the underlying Born diagrams to be considered
* for the given real emission process.
*/
const DiagramVector& underlyingBornDiagrams(const cPDVector& real) const;
/**
* Return the real emission diagrams to be considered
* for the given Born process.
*/
const DiagramVector& realEmissionDiagrams(const cPDVector& born) const;
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Return true, if this matrix element does not want to
* make use of mirroring processes; in this case all
* possible partonic subprocesses with a fixed assignment
* of incoming particles need to be provided through the diagrams
* added with the add(...) method.
*/
virtual bool noMirror () const { return true; }
/**
* With the information previously supplied with the
* setKinematics(...) method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Select a diagram. Default version uses diagrams(const
* DiagramVector &) to select a diagram according to the
* weights. This is the only method used that should be outside of
* MEBase.
*/
virtual DiagramIndex diagram(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Select a ColpurLines geometry. The default version returns a
* colour geometry selected among the ones returned from
* colourGeometries(tcDiagPtr).
*/
virtual const ColourLines &
selectColourGeometry(tcDiagPtr diag) const;
/**
* Return the order in \f$\alpha_S\f$ in which this matrix element
* is given.
*/
virtual unsigned int orderInAlphaS() const { return realEmissionME()->orderInAlphaS(); }
/**
* Return the order in \f$\alpha_{EM}\f$ in which this matrix
* element is given. Returns 0.
*/
virtual unsigned int orderInAlphaEW() const { return underlyingBornME()->orderInAlphaEW(); }
//@}
/** @name Phasespace generation */
//@{
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr xc);
/**
* Set the typed and momenta of the incoming and outgoing partons to
* be used in subsequent calls to me() and colourGeometries()
* according to the associated XComb object. If the function is
* overridden in a sub class the new function must call the base
* class one first.
*/
virtual void setKinematics();
/**
* Generate internal degrees of freedom given nDim() uniform random
* numbers in the interval ]0,1[. To help the phase space generator,
* the 'dSigHatDR' should be a smooth function of these numbers,
* although this is not strictly necessary. The return value should
* be true of the generation succeeded. If so the generated momenta
* should be stored in the meMomenta() vector.
*/
virtual bool generateKinematics(const double * r);
/**
* The number of internal degreed of freedom used in the matrix
* element. This default version returns 0;
*/
virtual int nDim() const;
/**
* Return true, if this matrix element expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS () const { return realEmissionME()->wantCMS(); }
/**
* Clear the information previously provided by a call to
* setKinematics(...).
*/
virtual void clearKinematics();
/**
* If this is a dependent matrix element in a ME group, return true,
* if cuts should be ignored.
*/
virtual bool ignoreCuts() const { return theIgnoreCuts; }
/**
* Indicate that cuts should be ignored
*/
void doIgnoreCuts(bool is = true) { theIgnoreCuts = is; }
//@}
/** @name Tilde kinematics */
//@{
/**
* Return the TildeKinematics object used
*/
Ptr<TildeKinematics>::tcptr tildeKinematics() const { return theTildeKinematics; }
/**
* Set the TildeKinematics object used
*/
void tildeKinematics(Ptr<TildeKinematics>::tptr);
/**
* Generate the tilde kinematics from real emission
* kinematics accessible through the XComb's
* head object and store it in meMomenta(). This default
* implemenation uses the tildeKinematics() object.
*/
virtual bool generateTildeKinematics();
/**
* Return the InvertedTildeKinematics object used
*/
Ptr<InvertedTildeKinematics>::tcptr invertedTildeKinematics() const { return theInvertedTildeKinematics; }
/**
* Set the InvertedTildeKinematics object used
*/
void invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr);
/**
* Return the number of additional random numbers
* needed to generate real emission kinematics off
* the tilde kinematics previously supplied through
* the XComb object. This default implementation
* returns invertedTildeKinematics()->nDimRadiation()
*/
virtual int nDimRadiation() const;
/**
* Generate the real emission kinematics
* off the Born kinematics accessible through the XComb's
* head object and store it in meMomenta(); store
* the single particle phasespace in units of lastHeadXComb()->lastSHat()
* in jacobian(). This default
* implemenation uses the invertedTildeKinematics() object
*/
virtual bool generateRadiationKinematics(const double *);
/**
* Set a pt cut when splitting
*/
void ptCut(Energy cut);
/**
* Return the relevant dipole scale
*/
Energy lastDipoleScale() const {
return splitting() ? theLastSplittingScale : theLastSubtractionScale;
}
/**
* Return the relevant pt
*/
Energy lastPt() const {
return splitting() ? theLastSplittingPt : theLastSubtractionPt;
}
/**
* Return true, if this dipole acts in splitting mode.
*/
bool splitting() const { return theSplitting; }
/**
* Switch on splitting mode for this dipole.
*/
void doSplitting() { theSplitting = true; }
/**
* Switch off splitting mode for this dipole.
*/
void doSubtraction() { theSplitting = false; }
/**
* Return the subtraction parameters.
*/
const vector<double>& subtractionParameters() const { return theSubtractionParameters; }
/**
* Access the subtraction parameters.
*/
vector<double>& subtractionParameters() { return theSubtractionParameters; }
//@}
/** @name Scale choices, couplings and PDFs */
//@{
/**
+ * Return true, if scales should be calculated from real emission kinematics
+ */
+ bool realEmissionScales() const { return theRealEmissionScales; }
+
+ /**
+ * Switch on or off that scales should be calculated from real emission kinematics
+ */
+ void doRealEmissionScales(bool on = true) { theRealEmissionScales = on; }
+
+ /**
* Return the scale associated with the phase space point provided
* by the last call to setKinematics().
*/
- virtual Energy2 scale() const { return realEmissionME()->scale(); }
+ virtual Energy2 scale() const {
+ return realEmissionScales() ?
+ realEmissionME()->scale() :
+ underlyingBornME()->scale();
+ }
/**
* Return the value of \f$\alpha_S\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaS(scale()).
*/
- virtual double alphaS() const { return realEmissionME()->lastAlphaS(); }
+ virtual double alphaS() const {
+ return realEmissionScales() ?
+ realEmissionME()->alphaS() :
+ underlyingBornME()->alphaS();
+ }
/**
* Return the value of \f$\alpha_EM\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaEM(scale()).
*/
- virtual double alphaEM() const { return realEmissionME()->lastAlphaEM(); }
+ virtual double alphaEM() const {
+ return realEmissionScales() ?
+ realEmissionME()->alphaEM() :
+ underlyingBornME()->alphaEM();
+ }
/**
* Return true, if this matrix element provides the PDF
* weight for the first incoming parton itself.
*/
virtual bool havePDFWeight1() const { return realEmissionME()->havePDFWeight1(); }
/**
* Return true, if this matrix element provides the PDF
* weight for the second incoming parton itself.
*/
virtual bool havePDFWeight2() const { return realEmissionME()->havePDFWeight2(); }
//@}
/** @name Matrix elements and evaluation */
//@{
/**
* Return the real emission matrix element
*/
Ptr<MatchboxMEBase>::tcptr realEmissionME() const {
return theRealEmissionME;
}
/**
* Return the real emission matrix element
*/
Ptr<MatchboxMEBase>::tptr realEmissionME() {
return theRealEmissionME;
}
/**
* Set the real emission matrix element
*/
void realEmissionME(Ptr<MatchboxMEBase>::tptr me) { theRealEmissionME = me; }
/**
* Return the underlying Born matrix element
*/
Ptr<MatchboxMEBase>::tcptr underlyingBornME() const {
return theUnderlyingBornME;
}
/**
* Return the underlying Born matrix element
*/
Ptr<MatchboxMEBase>::tptr underlyingBornME() {
return theUnderlyingBornME;
}
/**
* Set the underlying Born matrix element
*/
void underlyingBornME(Ptr<MatchboxMEBase>::tptr me) { theUnderlyingBornME = me; }
/**
* Set the dipoles which have been found along with this dipole
*/
void partnerDipoles(const vector<Ptr<SubtractionDipole>::ptr>& p) {
thePartners = p;
}
/**
* Return the dipoles which have been found along with this dipole
*/
const vector<Ptr<SubtractionDipole>::ptr>& partnerDipoles() const {
return thePartners;
}
/**
* Return the matrix element averaged over spin correlations.
*/
virtual double me2Avg(double ccme2) const = 0;
/**
* Return true, if the cross section should actually return the spin
* averaged splitting function times the Born matrix element squared.
*/
bool showerKernel() const { return theShowerKernel; }
/**
* Indicate that the cross section should actually return the spin
* averaged splitting function times the Born matrix element squared.
*/
void doShowerKernel(bool is = true) { theShowerKernel = is; }
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR(Energy2 factorizationScale) const;
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR() const { return dSigHatDR(ZERO); }
//@}
/** @name Methods relevant to matching */
//@{
/**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr app) {
theShowerApproximation = app;
}
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
/**
* Indicate that the shower real emission contribution should be subtracted.
*/
void doRealShowerSubtraction() { theRealShowerSubtraction = true; }
/**
* Return true, if the shower real emission contribution should be subtracted.
*/
bool realShowerSubtraction() const { return theRealShowerSubtraction; }
/**
* Indicate that the shower virtual contribution should be subtracted.
*/
void doVirtualShowerSubtraction() { theVirtualShowerSubtraction = true; }
/**
* Return true, if the shower virtual contribution should be subtracted.
*/
bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
//@}
/** @name Caching and diagnostic information */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches();
/**
* Indicate that the subtraction is being tested.
*/
void doTestSubtraction() { theSubtractionTest = true; }
/**
* Return true, if the subtraction is being tested.
*/
bool testSubtraction() const { return theSubtractionTest; }
/**
* Return true, if verbose
*/
bool verbose() const { return realEmissionME()->verbose() || underlyingBornME()->verbose(); }
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Print debug information on the last event
*/
virtual void printLastEvent(ostream&) const;
/**
* Write out diagnostic information for
* generateTildeKinematics
*/
void logGenerateTildeKinematics() const;
/**
* Write out diagnostic information for
* generateRadiationKinematics
*/
void logGenerateRadiationKinematics(const double * r) const;
/**
* Write out diagnostic information for
* me2 evaluation
*/
void logME2() const;
/**
* Write out diagnostic information
* for dsigdr evaluation
*/
void logDSigHatDR(double effectiveJac) const;
//@}
/** @name Reweight objects */
//@{
/**
* Insert a reweight object
*/
void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }
/**
* Return the reweight objects
*/
const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }
/**
* Access the reweight objects
*/
vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }
//@}
/** @name Methods used to setup SubtractionDipole objects */
//@{
/**
* Clone this dipole.
*/
Ptr<SubtractionDipole>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
void cloneDependencies(const std::string& prefix = "");
//@}
/** @name Methods required to setup the event record */
//@{
/**
* construct the spin information for the interaction
*/
virtual void constructVertex(tSubProPtr sub);
/**
* Comlete a SubProcess object using the internal degrees of freedom
* generated in the last generateKinematics() (and possible other
* degrees of freedom which was intergated over in dSigHatDR(). This
* default version does nothing. Will be made purely virtual in the
* future.
*/
virtual void generateSubCollision(SubProcess & sub);
//@}
protected:
/**
* Handle integer powers appearing downstream.
*/
double pow(double x, unsigned int p) const {
return std::pow(x,(double)p);
}
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();
private:
/**
* Wether or not this dipole acts in splitting mode.
*/
bool theSplitting;
/**
* True, if should apply to process in the xcomb.
*/
bool theApply;
/**
* True, if the subtraction is being tested.
*/
bool theSubtractionTest;
/**
* True if cuts should be ignored
*/
bool theIgnoreCuts;
/**
* True, if the cross section should actually return the spin
* averaged splitting function times the Born matrix element squared.
*/
bool theShowerKernel;
/**
* The real emission matrix element to be considered
*/
Ptr<MatchboxMEBase>::ptr theRealEmissionME;
/**
* The underlying Born matrix element
*/
Ptr<MatchboxMEBase>::ptr theUnderlyingBornME;
/**
* The dipoles which have been found along with this dipole
*/
vector<Ptr<SubtractionDipole>::ptr> thePartners;
/**
* The TildeKinematics to be used.
*/
Ptr<TildeKinematics>::ptr theTildeKinematics;
/**
* The InvertedTildeKinematics to be used.
*/
Ptr<InvertedTildeKinematics>::ptr theInvertedTildeKinematics;
/**
* A vector of reweight objects the sum of which
* should be applied to reweight this dipole
*/
vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
/**
* The emitter as referred to by the real emission
* matrix element.
*/
int theRealEmitter;
/**
* The emission as referred to by the real emission
* matrix element.
*/
int theRealEmission;
/**
* The spectator as referred to by the real emission
* matrix element.
*/
int theRealSpectator;
/**
* The subtraction parameters
*/
vector<double> theSubtractionParameters;
/**
* Map real emission diagrams to underlying Born diagrams
* and tilde emitter/spectator.
*/
map<RealEmissionKey,UnderlyingBornInfo> theMergingMap;
/**
* Map underlying Born diagrams and tilde emitter/spectator
* to real emission diagram containing the splitting.
*/
multimap<UnderlyingBornKey,RealEmissionInfo> theSplittingMap;
/**
* Map underlying Born diagrams to emitter/spectator pairs
*/
map<cPDVector,pair<int,int> > theIndexMap;
/**
* Map real emission processes to Born diagrams
*/
map<cPDVector,DiagramVector> theUnderlyingBornDiagrams;
/**
* Map Born processes to real emission diagrams
*/
map<cPDVector,DiagramVector> theRealEmissionDiagrams;
/**
* The last real emission key encountered
*/
RealEmissionKey lastRealEmissionKey;
/**
* The last underlying Born key encountered
*/
UnderlyingBornKey lastUnderlyingBornKey;
/**
* The last real emission info encountered
*/
multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator lastRealEmissionInfo;
/**
* The emitter as referred to by the underlying Born
* matrix element.
*/
int theBornEmitter;
/**
* The spectator as referred to by the underlying Born
* matrix element.
*/
int theBornSpectator;
/**
* The last scale as generated from the tilde mapping
*/
Energy theLastSubtractionScale;
/**
* The last scale as generated from the splitting mapping
*/
Energy theLastSplittingScale;
/**
* The last pt as generated from the tilde mapping
*/
Energy theLastSubtractionPt;
/**
* The last pt as generated from the splitting mapping
*/
Energy theLastSplittingPt;
/**
* The shower approximation.
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* True, if the shower real emission contribution should be subtracted.
*/
bool theRealShowerSubtraction;
/**
* True, if the shower virtual contribution should be subtracted.
*/
bool theVirtualShowerSubtraction;
+ /**
+ * True, if scales should be calculated from real emission kinematics
+ */
+ bool theRealEmissionScales;
+
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SubtractionDipole & operator=(const SubtractionDipole &);
};
}
#endif /* HERWIG_SubtractionDipole_H */
diff --git a/MatrixElement/Matchbox/MatchboxFactory.cc b/MatrixElement/Matchbox/MatchboxFactory.cc
--- a/MatrixElement/Matchbox/MatchboxFactory.cc
+++ b/MatrixElement/Matchbox/MatchboxFactory.cc
@@ -1,1182 +1,1204 @@
// -*- C++ -*-
//
// MatchboxFactory.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxFactory class.
//
#include "MatchboxFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/SU2Helper.h"
#include <boost/progress.hpp>
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
using std::ostream_iterator;
MatchboxFactory::MatchboxFactory()
: SubProcessHandler(), theNLight(0),
theOrderInAlphaS(0), theOrderInAlphaEW(0),
theBornContributions(true), theVirtualContributions(true),
theRealContributions(true), theIndependentVirtuals(false),
theSubProcessGroups(false), theInclusive(false),
theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0),
theFixedCouplings(false), theFixedQEDCouplings(false), theVetoScales(false),
- theVerbose(false), theInitVerbose(false), theSubtractionData(""), theCheckPoles(false) {}
+ theVerbose(false), theInitVerbose(false), theSubtractionData(""),
+ theCheckPoles(false), theRealEmissionScales(false) {}
MatchboxFactory::~MatchboxFactory() {}
IBPtr MatchboxFactory::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxFactory::fullclone() const {
return new_ptr(*this);
}
void MatchboxFactory::prepareME(Ptr<MatchboxMEBase>::ptr me) const {
Ptr<MatchboxAmplitude>::ptr amp =
dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>((*me).amplitude());
me->matchboxAmplitude(amp);
if ( diagramGenerator() && !me->diagramGenerator() )
me->diagramGenerator(diagramGenerator());
if ( processData() && !me->processData() )
me->processData(processData());
if ( me->nLight() == 0 )
me->nLight(nLight());
if ( phasespace() && !me->phasespace() )
me->phasespace(phasespace());
if ( scaleChoice() && !me->scaleChoice() )
me->scaleChoice(scaleChoice());
if ( me->factorizationScaleFactor() == 1.0 )
me->factorizationScaleFactor(factorizationScaleFactor());
if ( me->renormalizationScaleFactor() == 1.0 )
me->renormalizationScaleFactor(renormalizationScaleFactor());
if ( fixedCouplings() )
me->setFixedCouplings();
if ( fixedQEDCouplings() )
me->setFixedQEDCouplings();
if ( cache() && !me->cache() )
me->cache(cache());
if ( verbose() )
me->setVerbose();
}
struct LegIndex {
int spin;
int charge;
int colour;
int isSameAs;
int isSameFamilyAs;
inline bool operator==(const LegIndex& other) const {
return
spin == other.spin &&
charge == other.charge &&
colour == other.colour &&
isSameAs == other.isSameAs &&
isSameFamilyAs == other.isSameFamilyAs;
}
inline bool operator<(const LegIndex& other) const {
if ( spin != other.spin )
return spin < other.spin;
if ( charge != other.charge )
return charge < other.charge;
if ( colour != other.colour )
return colour < other.colour;
if ( isSameAs != other.isSameAs )
return isSameAs < other.isSameAs;
if ( isSameFamilyAs != other.isSameFamilyAs )
return isSameFamilyAs < other.isSameFamilyAs;
return false;
}
};
vector<LegIndex> makeIndex(const PDVector& proc) {
map<long,int> idMap;
map<int,int> familyIdMap;
int lastId = 0;
int lastFamilyId = 0;
vector<LegIndex> res;
for ( PDVector::const_iterator p = proc.begin();
p != proc.end(); ++p ) {
int id;
if ( idMap.find((**p).id()) != idMap.end() ) {
id = idMap[(**p).id()];
} else {
id = lastId;
idMap[(**p).id()] = lastId;
++lastId;
}
int familyId;
if ( familyIdMap.find(SU2Helper::family(*p)) != familyIdMap.end() ) {
familyId = familyIdMap[SU2Helper::family(*p)];
} else {
familyId = lastFamilyId;
familyIdMap[SU2Helper::family(*p)] = lastFamilyId;
++lastFamilyId;
}
LegIndex idx;
idx.spin = (**p).iSpin();
idx.charge = (**p).iCharge();
idx.colour = (**p).iColour();
idx.isSameAs = id;
idx.isSameFamilyAs = familyId;
res.push_back(idx);
}
return res;
}
string pid(const vector<LegIndex>& key) {
ostringstream res;
for ( vector<LegIndex>::const_iterator k =
key.begin(); k != key.end(); ++k )
res << k->spin << k->charge
<< k->colour << k->isSameAs
<< k->isSameFamilyAs;
return res.str();
}
vector<Ptr<MatchboxMEBase>::ptr> MatchboxFactory::
makeMEs(const vector<string>& proc, unsigned int orderas) const {
typedef vector<LegIndex> QNKey;
generator()->log() << "determining subprocesses for ";
copy(proc.begin(),proc.end(),ostream_iterator<string>(generator()->log()," "));
generator()->log() << flush;
map<Ptr<MatchboxAmplitude>::ptr,map<QNKey,vector<PDVector> > > ampProcs;
set<PDVector> processes = makeSubProcesses(proc);
vector<Ptr<MatchboxAmplitude>::ptr> matchAmplitudes;
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= amplitudes().begin(); amp != amplitudes().end(); ++amp ) {
(**amp).orderInGs(orderas);
(**amp).orderInGem(orderInAlphaEW());
if ( (**amp).orderInGs() != orderas ||
(**amp).orderInGem() != orderInAlphaEW() ) {
continue;
}
matchAmplitudes.push_back(*amp);
}
size_t combinations = processes.size()*matchAmplitudes.size();
size_t procCount = 0;
boost::progress_display * progressBar =
new boost::progress_display(combinations,generator()->log());
for ( vector<Ptr<MatchboxAmplitude>::ptr>::const_iterator amp
= matchAmplitudes.begin(); amp != matchAmplitudes.end(); ++amp ) {
(**amp).orderInGs(orderas);
(**amp).orderInGem(orderInAlphaEW());
for ( set<PDVector>::const_iterator p = processes.begin();
p != processes.end(); ++p ) {
++(*progressBar);
if ( !(**amp).canHandle(*p) )
continue;
QNKey key = makeIndex(*p);
++procCount;
ampProcs[*amp][key].push_back(*p);
}
}
delete progressBar;
generator()->log() << flush;
vector<Ptr<MatchboxMEBase>::ptr> res;
for ( map<Ptr<MatchboxAmplitude>::ptr,map<QNKey,vector<PDVector> > >::const_iterator
ap = ampProcs.begin(); ap != ampProcs.end(); ++ap ) {
for ( map<QNKey,vector<PDVector> >::const_iterator m = ap->second.begin();
m != ap->second.end(); ++m ) {
Ptr<MatchboxMEBase>::ptr me = ap->first->makeME(m->second);
me->subProcesses() = m->second;
me->amplitude(ap->first);
string pname = "ME" + ap->first->name() + pid(m->first);
if ( ! (generator()->preinitRegister(me,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
res.push_back(me);
}
}
generator()->log() << "created " << res.size()
<< " matrix element objects for "
<< procCount << " subprocesses.\n";
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
return res;
}
void MatchboxFactory::setup() {
if ( !amplitudes().empty() ) {
if ( particleGroups().find("j") == particleGroups().end() )
throw InitException() << "Could not find a jet particle group named 'j'";
// rebind the particle data objects
for ( map<string,PDVector>::iterator g = particleGroups().begin();
g != particleGroups().end(); ++g )
for ( PDVector::iterator p = g->second.begin();
p != g->second.end(); ++p ) {
#ifndef NDEBUG
long checkid = (**p).id();
#endif
*p = getParticleData((**p).id());
assert((**p).id() == checkid);
}
const PDVector& partons = particleGroups()["j"];
unsigned int nl = 0;
for ( PDVector::const_iterator p = partons.begin();
p != partons.end(); ++p )
if ( abs((**p).id()) < 6 )
++nl;
nLight(nl/2);
vector<Ptr<MatchboxMEBase>::ptr> ames = makeMEs(process,orderInAlphaS());
copy(ames.begin(),ames.end(),back_inserter(bornMEs()));
if ( realContributions() ) {
vector<string> rproc = process;
if ( realEmissionProcess.empty() ) {
rproc.push_back("j");
} else {
rproc = realEmissionProcess;
}
ames = makeMEs(rproc,orderInAlphaS()+1);
copy(ames.begin(),ames.end(),back_inserter(realEmissionMEs()));
}
}
// check if we have virtual contributions
bool haveVirtuals = true;
// check DR conventions of virtual contributions
bool virtualsAreDR = false;
bool virtualsAreCDR = false;
// check finite term conventions of virtual contributions
bool virtualsAreCS = false;
bool virtualsAreBDK = false;
bool virtualsAreExpanded = false;
// check and prepare the Born and virtual matrix elements
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
prepareME(*born);
haveVirtuals &= (**born).haveOneLoop();
if ( (**born).haveOneLoop() ) {
virtualsAreDR |= (**born).isDR();
virtualsAreCDR |= !(**born).isDR();
virtualsAreCS |= (**born).isCS();
virtualsAreBDK |= (**born).isBDK();
virtualsAreExpanded |= (**born).isExpanded();
}
}
// check the additional insertion operators
if ( !virtuals().empty() )
haveVirtuals = true;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
virtualsAreDR |= (**virt).isDR();
virtualsAreCDR |= !(**virt).isDR();
virtualsAreCS |= (**virt).isCS();
virtualsAreBDK |= (**virt).isBDK();
virtualsAreExpanded |= (**virt).isExpanded();
}
// check for consistent conventions on virtuals, if we are to include them
if ( virtualContributions() ) {
if ( virtualsAreDR && virtualsAreCDR ) {
throw InitException() << "Virtual corrections use inconsistent regularization schemes.\n";
}
if ( (virtualsAreCS && virtualsAreBDK) ||
(virtualsAreCS && virtualsAreExpanded) ||
(virtualsAreBDK && virtualsAreExpanded) ||
(!virtualsAreCS && !virtualsAreBDK && !virtualsAreExpanded) ) {
throw InitException() << "Virtual corrections use inconsistent conventions on finite terms.\n";
}
if ( !haveVirtuals ) {
throw InitException() << "Could not find amplitudes for all virtual contributions needed.\n";
}
}
// prepare dipole insertion operators
if ( virtualContributions() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( virtualsAreDR )
(**virt).useDR();
else
(**virt).useCDR();
if ( virtualsAreCS )
(**virt).useCS();
if ( virtualsAreBDK )
(**virt).useBDK();
if ( virtualsAreExpanded )
(**virt).useExpanded();
}
}
// prepare the real emission matrix elements
if ( realContributions() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
prepareME(*real);
}
}
// start creating matrix elements
MEs().clear();
// setup born and virtual contributions
if ( bornContributions() || virtualContributions() ) {
generator()->log() << "preparing Born "
<< (virtualContributions() ? "and virtual" : "")
<< " matrix elements." << flush;
}
if ( (bornContributions() && !virtualContributions()) ||
(bornContributions() && virtualContributions() && independentVirtuals()) ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
if ( (**born).onlyOneLoop() )
continue;
Ptr<MatchboxMEBase>::ptr bornme = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( independentVirtuals() )
pname += ".Born";
if ( ! (generator()->preinitRegister(bornme,pname) ) )
throw InitException() << "Matrix element " << pname << " already existing.";
bornme->cloneDependencies();
MEs().push_back(bornme);
}
}
if ( virtualContributions() ) {
bornVirtualMEs().clear();
boost::progress_display * progressBar =
new boost::progress_display(bornMEs().size(),generator()->log());
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born ) {
Ptr<MatchboxMEBase>::ptr nlo = (**born).cloneMe();
string pname = fullName() + "/" + (**born).name();
if ( !independentVirtuals() )
pname += ".BornVirtual";
else
pname += ".Virtual";
if ( ! (generator()->preinitRegister(nlo,pname) ) )
throw InitException() << "NLO ME " << pname << " already existing.";
nlo->virtuals().clear();
if ( !nlo->onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= virtuals().begin(); virt != virtuals().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator virt
= DipoleRepository::insertionOperators().begin();
virt != DipoleRepository::insertionOperators().end(); ++virt ) {
if ( (**virt).apply((**born).diagrams().front()->partons()) )
nlo->virtuals().push_back(*virt);
}
if ( nlo->virtuals().empty() )
throw InitException() << "No insertion operators have been found for "
<< (**born).name() << "\n";
if ( checkPoles() ) {
if ( !virtualsAreExpanded ) {
throw InitException() << "Cannot check epsilon poles if virtuals are not in `expanded' convention.\n";
}
nlo->doCheckPoles();
}
}
if ( !bornContributions() || independentVirtuals() ) {
nlo->doOneLoopNoBorn();
} else {
nlo->doOneLoop();
}
nlo->cloneDependencies();
bornVirtualMEs().push_back(nlo);
MEs().push_back(nlo);
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
theSplittingDipoles.clear();
set<cPDVector> bornProcs;
if ( showerApproximation() )
if ( showerApproximation()->needsSplittingGenerator() ) {
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator born
= bornMEs().begin(); born != bornMEs().end(); ++born )
for ( MEBase::DiagramVector::const_iterator d = (**born).diagrams().begin();
d != (**born).diagrams().end(); ++d )
bornProcs.insert((**d).partons());
}
if ( realContributions() ) {
generator()->log() << "preparing real emission matrix elements." << flush;
if ( theSubtractionData != "" )
if ( theSubtractionData[theSubtractionData.size()-1] != '/' )
theSubtractionData += "/";
subtractedMEs().clear();
boost::progress_display * progressBar =
new boost::progress_display(realEmissionMEs().size(),generator()->log());
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator real
= realEmissionMEs().begin(); real != realEmissionMEs().end(); ++real ) {
Ptr<SubtractedME>::ptr sub = new_ptr(SubtractedME());
string pname = fullName() + "/" + (**real).name() + ".Real";
if ( ! (generator()->preinitRegister(sub,pname) ) )
throw InitException() << "Subtracted ME " << pname << " already existing.";
sub->borns() = bornMEs();
sub->head(*real);
sub->allDipoles().clear();
sub->dependent().clear();
if ( subtractionData() != "" )
sub->subtractionData(subtractionData());
sub->getDipoles();
if ( sub->dependent().empty() ) {
// finite real contribution
Ptr<MatchboxMEBase>::ptr fme =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(sub->head())->cloneMe();
string pname = fullName() + "/" + (**real).name() + ".FiniteReal";
if ( ! (generator()->preinitRegister(fme,pname) ) )
throw InitException() << "ME " << pname << " already existing.";
MEs().push_back(fme);
finiteRealMEs().push_back(fme);
sub->head(tMEPtr());
continue;
}
if ( verbose() )
sub->setVerbose();
if ( subProcessGroups() )
sub->setSubProcessGroups();
if ( inclusive() )
sub->setInclusive();
if ( vetoScales() )
sub->doVetoScales();
+ if ( realEmissionScales() )
+ sub->doRealEmissionScales();
+
subtractedMEs().push_back(sub);
MEs().push_back(sub);
if ( showerApproximation() ) {
sub->showerApproximation(showerApproximation());
Ptr<SubtractedME>::ptr subv = new_ptr(*sub);
string vname = sub->fullName() + ".vsub";
if ( ! (generator()->preinitRegister(subv,pname) ) )
throw InitException() << "Subtracted ME " << vname << " already existing.";
sub->doRealShowerSubtraction();
subv->doVirtualShowerSubtraction();
subtractedMEs().push_back(subv);
MEs().push_back(subv);
if ( showerApproximation()->needsSplittingGenerator() )
for ( set<cPDVector>::const_iterator p = bornProcs.begin();
p != bornProcs.end(); ++p ) {
vector<Ptr<SubtractionDipole>::ptr> sdip = sub->splitDipoles(*p);
set<Ptr<SubtractionDipole>::ptr>& dips = theSplittingDipoles[*p];
copy(sdip.begin(),sdip.end(),inserter(dips,dips.begin()));
}
}
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
if ( !theSplittingDipoles.empty() ) {
map<Ptr<SubtractionDipole>::ptr,Ptr<SubtractionDipole>::ptr> cloneMap;
for ( map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::const_iterator sd = theSplittingDipoles.begin();
sd != theSplittingDipoles.end(); ++sd ) {
for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator d = sd->second.begin();
d != sd->second.end(); ++d ) {
cloneMap[*d] = Ptr<SubtractionDipole>::ptr();
}
}
for ( map<Ptr<SubtractionDipole>::ptr,Ptr<SubtractionDipole>::ptr>::iterator cd =
cloneMap.begin(); cd != cloneMap.end(); ++cd ) {
Ptr<SubtractionDipole>::ptr cloned = cd->first->cloneMe();
string dname = cd->first->fullName() + ".splitting";
if ( ! (generator()->preinitRegister(cloned,dname)) )
throw InitException() << "Dipole '" << dname << "' already existing.";
cloned->cloneDependencies();
cloned->showerApproximation(Ptr<ShowerApproximation>::tptr());
cloned->doSplitting();
cd->second = cloned;
}
for ( map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::iterator sd = theSplittingDipoles.begin();
sd != theSplittingDipoles.end(); ++sd ) {
set<Ptr<SubtractionDipole>::ptr> cloned;
for ( set<Ptr<SubtractionDipole>::ptr>::iterator d = sd->second.begin();
d != sd->second.end(); ++d ) {
cloned.insert(cloneMap[*d]);
}
sd->second = cloned;
}
}
generator()->log() << "process setup finished.\n" << flush;
}
list<MatchboxFactory::SplittingChannel>
MatchboxFactory::getSplittingChannels(tStdXCombPtr xcptr) const {
if ( xcptr->lastProjector() )
xcptr = xcptr->lastProjector();
const StandardXComb& xc = *xcptr;
cPDVector proc = xc.mePartonData();
map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >::const_iterator splitEntries
= splittingDipoles().find(proc);
list<SplittingChannel> res;
if ( splitEntries == splittingDipoles().end() )
return res;
const set<Ptr<SubtractionDipole>::ptr>& splitDipoles = splitEntries->second;
for ( set<Ptr<SubtractionDipole>::ptr>::const_iterator sd =
splitDipoles.begin(); sd != splitDipoles.end(); ++sd ) {
Ptr<MatchboxMEBase>::tptr bornME =
const_ptr_cast<Ptr<MatchboxMEBase>::tptr>((**sd).underlyingBornME());
SplittingChannel channel;
channel.dipole = *sd;
channel.bornXComb =
bornME->makeXComb(xc.maxEnergy(),xc.particles(),xc.eventHandlerPtr(),
const_ptr_cast<tSubHdlPtr>(xc.subProcessHandler()),
xc.pExtractor(),xc.CKKWHandler(),
xc.partonBins(),xc.cuts(),xc.diagrams(),xc.mirror(),
PartonPairVec());
vector<StdXCombPtr> realXCombs = (**sd).makeRealXCombs(channel.bornXComb);
for ( vector<StdXCombPtr>::const_iterator rxc = realXCombs.begin();
rxc != realXCombs.end(); ++rxc ) {
channel.realXComb = *rxc;
if ( showerApproximation()->needsTildeXCombs() ) {
channel.tildeXCombs.clear();
assert(!channel.dipole->partnerDipoles().empty());
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator p =
channel.dipole->partnerDipoles().begin();
p != channel.dipole->partnerDipoles().end(); ++p ) {
StdXCombPtr txc = channel.dipole->makeBornXComb(channel.realXComb);
if ( txc )
channel.tildeXCombs.push_back(txc);
}
}
res.push_back(channel);
}
}
return res;
}
void MatchboxFactory::print(ostream& os) const {
os << "--- MatchboxFactory setup -----------------------------------------------------------\n";
if ( !amplitudes().empty() ) {
os << " generated Born matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = bornMEs().begin();
m != bornMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
os << " generated real emission matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator m = realEmissionMEs().begin();
m != realEmissionMEs().end(); ++m ) {
(**m).print(os);
}
os << flush;
}
os << " generated Born+virtual matrix elements:\n";
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator bv
= bornVirtualMEs().begin(); bv != bornVirtualMEs().end(); ++bv ) {
(**bv).print(os);
}
os << " generated subtracted matrix elements:\n";
for ( vector<Ptr<SubtractedME>::ptr>::const_iterator sub
= subtractedMEs().begin(); sub != subtractedMEs().end(); ++sub ) {
os << " '" << (**sub).name() << "'\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxFactory::doinit() {
setup();
if ( initVerbose() )
print(Repository::clog());
SubProcessHandler::doinit();
}
void MatchboxFactory::persistentOutput(PersistentOStream & os) const {
os << theDiagramGenerator << theProcessData
<< theNLight << theOrderInAlphaS << theOrderInAlphaEW
<< theBornContributions << theVirtualContributions
<< theRealContributions << theIndependentVirtuals << theSubProcessGroups << theInclusive
<< thePhasespace << theScaleChoice
<< theFactorizationScaleFactor << theRenormalizationScaleFactor
<< theFixedCouplings << theFixedQEDCouplings << theVetoScales
<< theAmplitudes << theCache
<< theBornMEs << theVirtuals << theRealEmissionMEs
<< theBornVirtualMEs << theSubtractedMEs << theFiniteRealMEs
<< theVerbose << theInitVerbose << theSubtractionData << theCheckPoles
<< theParticleGroups << process << realEmissionProcess
- << theShowerApproximation << theSplittingDipoles;
+ << theShowerApproximation << theSplittingDipoles
+ << theRealEmissionScales;
}
void MatchboxFactory::persistentInput(PersistentIStream & is, int) {
is >> theDiagramGenerator >> theProcessData
>> theNLight >> theOrderInAlphaS >> theOrderInAlphaEW
>> theBornContributions >> theVirtualContributions
>> theRealContributions >> theIndependentVirtuals >> theSubProcessGroups >> theInclusive
>> thePhasespace >> theScaleChoice
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor
>> theFixedCouplings >> theFixedQEDCouplings >> theVetoScales
>> theAmplitudes >> theCache
>> theBornMEs >> theVirtuals >> theRealEmissionMEs
>> theBornVirtualMEs >> theSubtractedMEs >> theFiniteRealMEs
>> theVerbose >> theInitVerbose >> theSubtractionData >> theCheckPoles
>> theParticleGroups >> process >> realEmissionProcess
- >> theShowerApproximation >> theSplittingDipoles;
+ >> theShowerApproximation >> theSplittingDipoles
+ >> theRealEmissionScales;
}
string MatchboxFactory::startParticleGroup(string name) {
particleGroupName = StringUtils::stripws(name);
particleGroup.clear();
return "";
}
string MatchboxFactory::endParticleGroup(string) {
if ( particleGroup.empty() )
throw InitException() << "Empty particle group.";
particleGroups()[particleGroupName] = particleGroup;
particleGroup.clear();
return "";
}
string MatchboxFactory::doProcess(string in) {
process = StringUtils::split(in);
if ( process.size() < 3 )
throw InitException() << "Invalid process.";
for ( vector<string>::iterator p = process.begin();
p != process.end(); ++p ) {
*p = StringUtils::stripws(*p);
}
return "";
}
string MatchboxFactory::doSingleRealProcess(string in) {
realEmissionProcess = StringUtils::split(in);
if ( realEmissionProcess.size() < 3 )
throw InitException() << "Invalid process.";
for ( vector<string>::iterator p = realEmissionProcess.begin();
p != realEmissionProcess.end(); ++p ) {
*p = StringUtils::stripws(*p);
}
return "";
}
struct SortPID {
inline bool operator()(PDPtr a, PDPtr b) const {
return a->id() < b->id();
}
};
set<PDVector> MatchboxFactory::
makeSubProcesses(const vector<string>& proc) const {
if ( proc.empty() )
throw InitException() << "No process specified.";
vector<PDVector> allProcs(1);
size_t pos = 0;
typedef map<string,PDVector>::const_iterator GroupIterator;
while ( pos < proc.size() ) {
GroupIterator git =
particleGroups().find(proc[pos]);
if ( git == particleGroups().end() ) {
throw InitException() << "particle group '"
<< proc[pos] << "' not defined.";
}
vector<PDVector> mine;
for ( vector<PDVector>::const_iterator i = allProcs.begin();
i != allProcs.end(); ++i ) {
for ( PDVector::const_iterator p = git->second.begin();
p != git->second.end(); ++p ) {
PDVector v = *i;
v.push_back(*p);
mine.push_back(v);
}
}
allProcs = mine;
++pos;
}
set<PDVector> allCheckedProcs;
for ( vector<PDVector>::const_iterator p = allProcs.begin();
p != allProcs.end(); ++p ) {
int charge = -(*p)[0]->iCharge() -(*p)[1]->iCharge();
for ( size_t k = 2; k < (*p).size(); ++k )
charge += (*p)[k]->iCharge();
if ( charge != 0 )
continue;
PDVector pr = *p;
sort(pr.begin()+2,pr.end(),SortPID());
allCheckedProcs.insert(pr);
}
return allCheckedProcs;
}
void MatchboxFactory::Init() {
static ClassDocumentation<MatchboxFactory> documentation
("MatchboxFactory",
"NLO QCD corrections have been calculated "
"using Matchbox \\cite{Platzer:2011bc}",
"%\\cite{Platzer:2011bc}\n"
"\\bibitem{Platzer:2011bc}\n"
"S.~Platzer and S.~Gieseke,\n"
"``Dipole Showers and Automated NLO Matching in Herwig++,''\n"
"arXiv:1109.6256 [hep-ph].\n"
"%%CITATION = ARXIV:1109.6256;%%");
static Reference<MatchboxFactory,Tree2toNGenerator> interfaceDiagramGenerator
("DiagramGenerator",
"Set the diagram generator.",
&MatchboxFactory::theDiagramGenerator, false, false, true, true, false);
static Reference<MatchboxFactory,ProcessData> interfaceProcessData
("ProcessData",
"Set the process data object to be used.",
&MatchboxFactory::theProcessData, false, false, true, true, false);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaS
("OrderInAlphaS",
"The order in alpha_s to consider.",
&MatchboxFactory::theOrderInAlphaS, 0, 0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,unsigned int> interfaceOrderInAlphaEW
("OrderInAlphaEW",
"The order in alpha_EW",
&MatchboxFactory::theOrderInAlphaEW, 2, 0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceBornContributions
("BornContributions",
"Switch on or off the Born contributions.",
&MatchboxFactory::theBornContributions, true, false, false);
static SwitchOption interfaceBornContributionsOn
(interfaceBornContributions,
"On",
"Switch on Born contributions.",
true);
static SwitchOption interfaceBornContributionsOff
(interfaceBornContributions,
"Off",
"Switch off Born contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceVirtualContributions
("VirtualContributions",
"Switch on or off the virtual contributions.",
&MatchboxFactory::theVirtualContributions, true, false, false);
static SwitchOption interfaceVirtualContributionsOn
(interfaceVirtualContributions,
"On",
"Switch on virtual contributions.",
true);
static SwitchOption interfaceVirtualContributionsOff
(interfaceVirtualContributions,
"Off",
"Switch off virtual contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceRealContributions
("RealContributions",
"Switch on or off the real contributions.",
&MatchboxFactory::theRealContributions, true, false, false);
static SwitchOption interfaceRealContributionsOn
(interfaceRealContributions,
"On",
"Switch on real contributions.",
true);
static SwitchOption interfaceRealContributionsOff
(interfaceRealContributions,
"Off",
"Switch off real contributions.",
false);
static Switch<MatchboxFactory,bool> interfaceIndependentVirtuals
("IndependentVirtuals",
"Switch on or off virtual contributions as separate subprocesses.",
&MatchboxFactory::theIndependentVirtuals, true, false, false);
static SwitchOption interfaceIndependentVirtualsOn
(interfaceIndependentVirtuals,
"On",
"Switch on virtual contributions as separate subprocesses.",
true);
static SwitchOption interfaceIndependentVirtualsOff
(interfaceIndependentVirtuals,
"Off",
"Switch off virtual contributions as separate subprocesses.",
false);
static Switch<MatchboxFactory,bool> interfaceSubProcessGroups
("SubProcessGroups",
"Switch on or off production of sub-process groups.",
&MatchboxFactory::theSubProcessGroups, false, false, false);
static SwitchOption interfaceSubProcessGroupsOn
(interfaceSubProcessGroups,
"On",
"On",
true);
static SwitchOption interfaceSubProcessGroupsOff
(interfaceSubProcessGroups,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceInclusive
("Inclusive",
"Switch on or off production of inclusive cross section.",
&MatchboxFactory::theInclusive, false, false, false);
static SwitchOption interfaceInclusiveOn
(interfaceInclusive,
"On",
"On",
true);
static SwitchOption interfaceInclusiveOff
(interfaceInclusive,
"Off",
"Off",
false);
static Reference<MatchboxFactory,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"Set the phasespace generator.",
&MatchboxFactory::thePhasespace, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxScaleChoice> interfaceScaleChoice
("ScaleChoice",
"Set the scale choice object.",
&MatchboxFactory::theScaleChoice, false, false, true, true, false);
static Parameter<MatchboxFactory,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&MatchboxFactory::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<MatchboxFactory,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&MatchboxFactory::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Switch<MatchboxFactory,bool> interfaceFixedCouplings
("FixedCouplings",
"Switch on or off fixed couplings.",
&MatchboxFactory::theFixedCouplings, true, false, false);
static SwitchOption interfaceFixedCouplingsOn
(interfaceFixedCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedCouplingsOff
(interfaceFixedCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceFixedQEDCouplings
("FixedQEDCouplings",
"Switch on or off fixed QED couplings.",
&MatchboxFactory::theFixedQEDCouplings, true, false, false);
static SwitchOption interfaceFixedQEDCouplingsOn
(interfaceFixedQEDCouplings,
"On",
"On",
true);
static SwitchOption interfaceFixedQEDCouplingsOff
(interfaceFixedQEDCouplings,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceVetoScales
("VetoScales",
"Switch on or setting veto scales.",
&MatchboxFactory::theVetoScales, false, false, false);
static SwitchOption interfaceVetoScalesOn
(interfaceVetoScales,
"On",
"On",
true);
static SwitchOption interfaceVetoScalesOff
(interfaceVetoScales,
"Off",
"Off",
false);
static RefVector<MatchboxFactory,MatchboxAmplitude> interfaceAmplitudes
("Amplitudes",
"The amplitude objects.",
&MatchboxFactory::theAmplitudes, -1, false, false, true, true, false);
static Reference<MatchboxFactory,MatchboxMECache> interfaceCache
("Cache",
"Set the matrix element cache object.",
&MatchboxFactory::theCache, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceBornMEs
("BornMEs",
"The Born matrix elements to be used",
&MatchboxFactory::theBornMEs, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxInsertionOperator> interfaceVirtuals
("Virtuals",
"The virtual corrections to include",
&MatchboxFactory::theVirtuals, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceRealEmissionMEs
("RealEmissionMEs",
"The RealEmission matrix elements to be used",
&MatchboxFactory::theRealEmissionMEs, -1, false, false, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceBornVirtuals
("BornVirtualMEs",
"The generated Born/virtual contributions",
&MatchboxFactory::theBornVirtualMEs, -1, false, true, true, true, false);
static RefVector<MatchboxFactory,SubtractedME> interfaceSubtractedMEs
("SubtractedMEs",
"The generated subtracted real emission contributions",
&MatchboxFactory::theSubtractedMEs, -1, false, true, true, true, false);
static RefVector<MatchboxFactory,MatchboxMEBase> interfaceFiniteRealMEs
("FiniteRealMEs",
"The generated finite real contributions",
&MatchboxFactory::theFiniteRealMEs, -1, false, true, true, true, false);
static Switch<MatchboxFactory,bool> interfaceVerbose
("Verbose",
"Print full infomation on each evaluated phase space point.",
&MatchboxFactory::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"On",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"Off",
false);
static Switch<MatchboxFactory,bool> interfaceInitVerbose
("InitVerbose",
"Print setup information.",
&MatchboxFactory::theInitVerbose, false, false, false);
static SwitchOption interfaceInitVerboseOn
(interfaceInitVerbose,
"On",
"On",
true);
static SwitchOption interfaceInitVerboseOff
(interfaceInitVerbose,
"Off",
"Off",
false);
static Parameter<MatchboxFactory,string> interfaceSubtractionData
("SubtractionData",
"Prefix for subtraction check data.",
&MatchboxFactory::theSubtractionData, "",
false, false);
static Switch<MatchboxFactory,bool> interfaceCheckPoles
("CheckPoles",
"Switch on or off checks of epsilon poles.",
&MatchboxFactory::theCheckPoles, true, false, false);
static SwitchOption interfaceCheckPolesOn
(interfaceCheckPoles,
"On",
"On",
true);
static SwitchOption interfaceCheckPolesOff
(interfaceCheckPoles,
"Off",
"Off",
false);
static RefVector<MatchboxFactory,ParticleData> interfaceParticleGroup
("ParticleGroup",
"The particle group just started.",
&MatchboxFactory::particleGroup, -1, false, false, true, false, false);
static Command<MatchboxFactory> interfaceStartParticleGroup
("StartParticleGroup",
"Start a particle group.",
&MatchboxFactory::startParticleGroup, false);
static Command<MatchboxFactory> interfaceEndParticleGroup
("EndParticleGroup",
"End a particle group.",
&MatchboxFactory::endParticleGroup, false);
static Command<MatchboxFactory> interfaceProcess
("Process",
"Set the process to consider.",
&MatchboxFactory::doProcess, false);
static Command<MatchboxFactory> interfaceSingleRealProcess
("SingleRealProcess",
"Set the process to consider.",
&MatchboxFactory::doSingleRealProcess, false);
static Reference<MatchboxFactory,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to be considered.",
&MatchboxFactory::theShowerApproximation, false, false, true, true, false);
+ static Switch<MatchboxFactory,bool> interfaceRealEmissionScales
+ ("RealEmissionScales",
+ "Switch on or off calculation of subtraction scales from real emission kinematics.",
+ &MatchboxFactory::theRealEmissionScales, true, false, false);
+ static SwitchOption interfaceRealEmissionScalesOn
+ (interfaceRealEmissionScales,
+ "On",
+ "On",
+ true);
+ static SwitchOption interfaceRealEmissionScalesOff
+ (interfaceRealEmissionScales,
+ "Off",
+ "Off",
+ false);
+
+
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MatchboxFactory,SubProcessHandler>
describeHerwigMatchboxFactory("Herwig::MatchboxFactory", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/MatchboxFactory.h b/MatrixElement/Matchbox/MatchboxFactory.h
--- a/MatrixElement/Matchbox/MatchboxFactory.h
+++ b/MatrixElement/Matchbox/MatchboxFactory.h
@@ -1,781 +1,796 @@
// -*- C++ -*-
//
// MatchboxFactory.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_MatchboxFactory_H
#define HERWIG_MatchboxFactory_H
//
// This is the declaration of the MatchboxFactory class.
//
#include "ThePEG/Handlers/SubProcessHandler.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/Tree2toNGenerator.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/ProcessData.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxScaleChoice.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxMECache.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Base/SubtractedME.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxFactory automatically sets up a NLO
* QCD calculation carried out in dipole subtraction.
*
* @see \ref MatchboxFactoryInterfaces "The interfaces"
* defined for MatchboxFactory.
*/
class MatchboxFactory: public SubProcessHandler {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxFactory();
/**
* The destructor.
*/
virtual ~MatchboxFactory();
//@}
public:
/** @name Process and diagram information */
//@{
/**
* Return the diagram generator.
*/
Ptr<Tree2toNGenerator>::tptr diagramGenerator() const { return theDiagramGenerator; }
/**
* Set the diagram generator.
*/
void diagramGenerator(Ptr<Tree2toNGenerator>::ptr dg) { theDiagramGenerator = dg; }
/**
* Return the process data.
*/
Ptr<ProcessData>::tptr processData() const { return theProcessData; }
/**
* Set the process data.
*/
void processData(Ptr<ProcessData>::ptr pd) { theProcessData = pd; }
/**
* Return the number of light flavours, this matrix
* element is calculated for.
*/
unsigned int nLight() const { return theNLight; }
/**
* Set the number of light flavours, this matrix
* element is calculated for.
*/
void nLight(unsigned int n) { theNLight = n; }
/**
* Return the order in \f$\alpha_S\f$.
*/
unsigned int orderInAlphaS() const { return theOrderInAlphaS; }
/**
* Set the order in \f$\alpha_S\f$.
*/
void orderInAlphaS(unsigned int o) { theOrderInAlphaS = o; }
/**
* Return the order in \f$\alpha_{EM}\f$.
*/
unsigned int orderInAlphaEW() const { return theOrderInAlphaEW; }
/**
* Set the order in \f$\alpha_{EM}\f$.
*/
void orderInAlphaEW(unsigned int o) { theOrderInAlphaEW = o; }
/**
* Return true, if Born contributions should be included.
*/
bool bornContributions() const { return theBornContributions; }
/**
* Switch on or off Born contributions
*/
void setBornContributions(bool on = true) { theBornContributions = on; }
/**
* Return true, if virtual contributions should be included.
*/
bool virtualContributions() const { return theVirtualContributions; }
/**
* Switch on or off virtual contributions
*/
void setVirtualContributions(bool on = true) { theVirtualContributions = on; }
/**
* Return true, if subtracted real emission contributions should be included.
*/
bool realContributions() const { return theRealContributions; }
/**
* Switch on or off subtracted real emission contributions
*/
void setRealContributions(bool on = true) { theRealContributions = on; }
/**
* Return true, if virtual contributions should be treated as independent subprocesses
*/
bool independentVirtuals() const { return theIndependentVirtuals; }
/**
* Switch on/off virtual contributions should be treated as independent subprocesses
*/
void setIndependentVirtuals(bool on = true) { theIndependentVirtuals = on; }
/**
* Return true, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool subProcessGroups() const { return theSubProcessGroups; }
/**
* Switch on or off producing subprocess groups.
*/
void setSubProcessGroups(bool on = true) { theSubProcessGroups = on; }
/**
+ * Return true, if subtraction scales should be caluclated from real emission kinematics
+ */
+ bool realEmissionScales() const { return theRealEmissionScales; }
+
+ /**
+ * Switch on/off that subtraction scales should be caluclated from real emission kinematics
+ */
+ void setRealEmissionScales(bool on = true) { theRealEmissionScales = on; }
+
+ /**
* Return true, if the integral over the unresolved emission should be
* calculated.
*/
bool inclusive() const { return theInclusive; }
/**
* Switch on or off inclusive mode.
*/
void setInclusive(bool on = true) { theInclusive = on; }
/**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr app) { theShowerApproximation = app; }
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
//@}
/** @name Phasespace generation and scale choice */
//@{
/**
* Return the phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::tptr phasespace() const { return thePhasespace; }
/**
* Set the phase space generator to be used.
*/
void phasespace(Ptr<MatchboxPhasespace>::ptr ps) { thePhasespace = ps; }
/**
* Set the scale choice object
*/
void scaleChoice(Ptr<MatchboxScaleChoice>::ptr sc) { theScaleChoice = sc; }
/**
* Return the scale choice object
*/
Ptr<MatchboxScaleChoice>::tptr scaleChoice() const { return theScaleChoice; }
/**
* Get the factorization scale factor
*/
double factorizationScaleFactor() const { return theFactorizationScaleFactor; }
/**
* Set the factorization scale factor
*/
void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }
/**
* Get the renormalization scale factor
*/
double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }
/**
* Set the renormalization scale factor
*/
void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }
/**
* Return true, if fixed couplings are used.
*/
bool fixedCouplings() const { return theFixedCouplings; }
/**
* Switch on fixed couplings.
*/
void setFixedCouplings(bool on = true) { theFixedCouplings = on; }
/**
* Return true, if fixed couplings are used.
*/
bool fixedQEDCouplings() const { return theFixedQEDCouplings; }
/**
* Switch on fixed couplings.
*/
void setFixedQEDCouplings(bool on = true) { theFixedQEDCouplings = on; }
/**
* Return true, if veto scales should be set
* for the real emission
*/
bool vetoScales() const { return theVetoScales; }
/**
* Switch on setting veto scales
*/
void doVetoScales() { theVetoScales = true; }
/**
* Switch off setting veto scales
*/
void noVetoScales() { theVetoScales = true; }
//@}
/** @name Amplitudes and caching */
//@{
/**
* Return the amplitudes to be considered
*/
const vector<Ptr<MatchboxAmplitude>::ptr>& amplitudes() const { return theAmplitudes; }
/**
* Access the amplitudes to be considered
*/
vector<Ptr<MatchboxAmplitude>::ptr>& amplitudes() { return theAmplitudes; }
/**
* Set the ME cache object
*/
void cache(Ptr<MatchboxMECache>::ptr c) { theCache = c; }
/**
* Get the ME cache object
*/
Ptr<MatchboxMECache>::tptr cache() const { return theCache; }
//@}
/** @name Matrix element objects. */
//@{
/**
* Return the Born matrix elements to be considered
*/
const vector<Ptr<MatchboxMEBase>::ptr>& bornMEs() const { return theBornMEs; }
/**
* Access the Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr>& bornMEs() { return theBornMEs; }
/**
* Return the virtual corrections to be considered
*/
const vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() const { return theVirtuals; }
/**
* Access the virtual corrections to be considered
*/
vector<Ptr<MatchboxInsertionOperator>::ptr>& virtuals() { return theVirtuals; }
/**
* Return the produced NLO matrix elements
*/
const vector<Ptr<MatchboxMEBase>::ptr>& bornVirtualMEs() const { return theBornVirtualMEs; }
/**
* Access the produced NLO matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr>& bornVirtualMEs() { return theBornVirtualMEs; }
/**
* Return the real emission matrix elements to be considered
*/
const vector<Ptr<MatchboxMEBase>::ptr>& realEmissionMEs() const { return theRealEmissionMEs; }
/**
* Access the real emission matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr>& realEmissionMEs() { return theRealEmissionMEs; }
/**
* Return the produced subtracted matrix elements
*/
const vector<Ptr<SubtractedME>::ptr>& subtractedMEs() const { return theSubtractedMEs; }
/**
* Access the produced subtracted matrix elements
*/
vector<Ptr<SubtractedME>::ptr>& subtractedMEs() { return theSubtractedMEs; }
/**
* Return the produced finite real emission matrix elements
*/
const vector<Ptr<MatchboxMEBase>::ptr>& finiteRealMEs() const { return theFiniteRealMEs; }
/**
* Access the produced finite real emission elements
*/
vector<Ptr<MatchboxMEBase>::ptr>& finiteRealMEs() { return theFiniteRealMEs; }
/**
* Return the map of Born processes to splitting dipoles
*/
const map<cPDVector,set<Ptr<SubtractionDipole>::ptr> >& splittingDipoles() const {
return theSplittingDipoles;
}
/**
* Identify a splitting channel
*/
struct SplittingChannel {
/**
* The Born XComb
*/
StdXCombPtr bornXComb;
/**
* The real XComb
*/
StdXCombPtr realXComb;
/**
* The set of tilde XCombs to consider for the real xcomb
*/
vector<StdXCombPtr> tildeXCombs;
/**
* The dipole in charge of the splitting
*/
Ptr<SubtractionDipole>::ptr dipole;
};
/**
* Generate all splitting channels for the Born process handled by
* the given XComb
*/
list<SplittingChannel> getSplittingChannels(tStdXCombPtr xc) const;
//@}
/** @name Setup the matrix elements */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return true; }
/**
* Prepare a matrix element.
*/
void prepareME(Ptr<MatchboxMEBase>::ptr) const;
/**
* Setup everything
*/
void setup();
//@}
/** @name Diagnostic information */
//@{
/**
* Return true, if verbose
*/
bool verbose() const { return theVerbose; }
/**
* Switch on diagnostic information.
*/
void setVerbose(bool on = true) { theVerbose = on; }
/**
* Return true, if verbose while initializing
*/
bool initVerbose() const { return theInitVerbose || verbose(); }
/**
* Switch on diagnostic information while initializing
*/
void setInitVerbose(bool on = true) { theInitVerbose = on; }
/**
* Dump the setup
*/
void print(ostream&) const;
/**
* Return the subtraction data prefix.
*/
const string& subtractionData() const { return theSubtractionData; }
/**
* Set the subtraction data prefix.
*/
void subtractionData(const string& s) { theSubtractionData = s; }
/**
* Return true, if cancellationn of epsilon poles should be checked.
*/
bool checkPoles() const { return theCheckPoles; }
/**
* Switch on checking of epsilon pole cancellation.
*/
void doCheckPoles() { theCheckPoles = true; }
//@}
/** @name Process generation */
//@{
/**
* Return the particle groups.
*/
const map<string,PDVector>& particleGroups() const { return theParticleGroups; }
/**
* Access the particle groups.
*/
map<string,PDVector>& particleGroups() { return theParticleGroups; }
//@}
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;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The diagram generator.
*/
Ptr<Tree2toNGenerator>::ptr theDiagramGenerator;
/**
* The process data object to be used
*/
Ptr<ProcessData>::ptr theProcessData;
/**
* The number of light flavours, this matrix
* element is calculated for.
*/
unsigned int theNLight;
/**
* The order in \f$\alpha_S\f$.
*/
unsigned int theOrderInAlphaS;
/**
* The order in \f$\alpha_{EM}\f$.
*/
unsigned int theOrderInAlphaEW;
/**
* Switch on or off Born contributions
*/
bool theBornContributions;
/**
* Switch on or off virtual contributions
*/
bool theVirtualContributions;
/**
* Switch on or off subtracted real emission contributions should be included.
*/
bool theRealContributions;
/**
* True if virtual contributions should be treated as independent subprocesses
*/
bool theIndependentVirtuals;
/**
* True, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
bool theSubProcessGroups;
/**
* True, if the integral over the unresolved emission should be
* calculated.
*/
bool theInclusive;
/**
* The phase space generator to be used.
*/
Ptr<MatchboxPhasespace>::ptr thePhasespace;
/**
* The scale choice object
*/
Ptr<MatchboxScaleChoice>::ptr theScaleChoice;
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
/**
* Use non-running couplings.
*/
bool theFixedCouplings;
/**
* Use non-running couplings.
*/
bool theFixedQEDCouplings;
/**
* True, if veto scales should be set
* for the real emission
*/
bool theVetoScales;
/**
* The amplitudes to be considered
*/
vector<Ptr<MatchboxAmplitude>::ptr> theAmplitudes;
/**
* The ME cache object
*/
Ptr<MatchboxMECache>::ptr theCache;
/**
* The Born matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theBornMEs;
/**
* The virtual corrections to be considered
*/
vector<Ptr<MatchboxInsertionOperator>::ptr> theVirtuals;
/**
* The real emission matrix elements to be considered
*/
vector<Ptr<MatchboxMEBase>::ptr> theRealEmissionMEs;
/**
* The produced NLO matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr> theBornVirtualMEs;
/**
* The produced subtracted matrix elements
*/
vector<Ptr<SubtractedME>::ptr> theSubtractedMEs;
/**
* The produced finite real emission matrix elements
*/
vector<Ptr<MatchboxMEBase>::ptr> theFiniteRealMEs;
/**
* Switch on or off verbosity
*/
bool theVerbose;
/**
* True, if verbose while initializing
*/
bool theInitVerbose;
/**
* Prefix for subtraction data
*/
string theSubtractionData;
/**
* Command to limit the real emission process to be considered.
*/
string doSingleRealProcess(string);
/**
* The real emission process to be included; if empty, all possible
* ones will be considered.
*/
vector<string> realEmissionProcess;
/**
* True, if cancellationn of epsilon poles should be checked.
*/
bool theCheckPoles;
/**
* Particle groups.
*/
map<string,PDVector> theParticleGroups;
/**
* Command to start a particle group.
*/
string startParticleGroup(string);
/**
* The name of the particle group currently edited.
*/
string particleGroupName;
/**
* The particle group currently edited.
*/
PDVector particleGroup;
/**
* Command to end a particle group.
*/
string endParticleGroup(string);
/**
* Command to set the process.
*/
string doProcess(string);
/**
* The process to consider in terms of particle groups.
*/
vector<string> process;
/**
* Generate subprocesses.
*/
set<PDVector> makeSubProcesses(const vector<string>&) const;
/**
* Generate matrix element objects for the given process.
*/
vector<Ptr<MatchboxMEBase>::ptr> makeMEs(const vector<string>&,
unsigned int orderas) const;
/**
* The shower approximation.
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* The map of Born processes to splitting dipoles
*/
map<cPDVector,set<Ptr<SubtractionDipole>::ptr> > theSplittingDipoles;
+ /**
+ * True, if subtraction scales should be caluclated from real emission kinematics
+ */
+ bool theRealEmissionScales;
+
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxFactory & operator=(const MatchboxFactory &);
};
}
#endif /* HERWIG_MatchboxFactory_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:08 PM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805519
Default Alt Text
(280 KB)

Event Timeline