Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7877817
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
68 KB
Subscribers
None
View Options
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,790 +1,793 @@
// -*- 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"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
using namespace Herwig;
SubtractedME::SubtractedME()
: MEGroup(),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
theLoopSimSubtraction(false), theSubProcessGroups(false) {}
SubtractedME::~SubtractedME() {}
Ptr<MatchboxFactory>::tcptr SubtractedME::factory() const { return theFactory; }
void SubtractedME::factory(Ptr<MatchboxFactory>::tcptr f) { theFactory = f; }
bool SubtractedME::subProcessGroups() const {
return
(factory()->subProcessGroups() && !showerApproximation()) ||
factory()->subtractionData() != "" || theSubProcessGroups;
}
Ptr<ShowerApproximation>::tptr SubtractedME::showerApproximation() const { return factory()->showerApproximation(); }
const vector<Ptr<MatchboxMEBase>::ptr>& SubtractedME::borns() const {
return theBorns.empty() ? factory()->bornMEs() : theBorns;
}
bool SubtractedME::verbose() const { return factory()->verbose(); }
bool SubtractedME::initVerbose() const { return factory()->initVerbose(); }
IBPtr SubtractedME::clone() const {
return new_ptr(*this);
}
IBPtr SubtractedME::fullclone() const {
return new_ptr(*this);
}
StdXCombPtr SubtractedME::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allPBins,
tStdXCombPtr newHead,
tMEPtr newME) {
tMEGroupPtr newMEGroup = dynamic_ptr_cast<tMEGroupPtr>(newME);
if ( !newMEGroup )
newMEGroup = this;
Ptr<MatchboxXCombGroup>::ptr res =
new_ptr(MatchboxXCombGroup(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts, newMEGroup,
newDiagrams, mir,
newHead));
res->build(allPBins);
theReal->prepareXComb(*res);
if ( factory()->subtractionData() != "" ) {
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.00001,1000.);
ostringstream fname("");
fname << factory()->subtractionData();
const cPDVector& myproc = SoftSubtractionIndex(*p,i).first;
for (cPDVector::const_iterator pp = myproc.begin(); pp != myproc.end(); ++pp) fname << (**pp).PDGName();
fname << "-" << i << "-" << i << "-scatter.dat";
fnamesSoftSubtraction[SoftSubtractionIndex(*p,i)] = fname.str();
if ( theReal->phasespace() )
res->singularLimits().insert(make_pair(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.00001,1000.);
ostringstream fname("");
fname << factory()->subtractionData();
const cPDVector& myproc = CollinearSubtractionIndex(*p,make_pair(i,j)).first;
for (cPDVector::const_iterator pp = myproc.begin(); pp != myproc.end(); ++pp) fname << (**pp).PDGName();
fname << "-" << i << "-" << j << "-scatter.dat";
fnamesCollinearSubtraction[CollinearSubtractionIndex(*p,make_pair(i,j))] = fname.str();
if ( theReal->phasespace() )
res->singularLimits().insert(make_pair(i,j));
}
}
}
}
return res;
}
void SubtractedME::setXComb(tStdXCombPtr xc) {
MEGroup::setXComb(xc);
lastMatchboxXComb(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 Exception() << "SubtractedME: A dependent matrix element of SubtractedME "
<< "has not been derived from SubtractionDipole. "
<< "Please check the corresponding input file." << Exception::runerror;
}
return dipole->underlyingBornDiagrams(proc);
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::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;
Ptr<MatchboxMEBase>::tptr real =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( borns().empty() || !real )
throw Exception() << "SubtractedME: The SubtractedME '"
<< name() << "' could not generate "
<< "subtraction terms for the real emission "
<< "matrix element '" << real->name() << "'. "
<< "Please check the corresponding input file." << Exception::runerror;
Ptr<MatchboxMEBase>::ptr myRealEmissionME = real->cloneMe();
ostringstream pname;
pname << fullName() << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw Exception() << "SubtractedME: Matrix element " << pname.str() << " already existing." << Exception::runerror;
myRealEmissionME->cloneDependencies(pname.str());
head(myRealEmissionME);
real = myRealEmissionME;
dependent().clear();
vector<Ptr<SubtractionDipole>::ptr> genDipoles
= real->getDipoles(DipoleRepository::dipoles(factory()->dipoleSet()),borns());
if ( factory()->subtractionData() != "" ) {
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
genDipoles.begin(); d != genDipoles.end(); ++d )
(**d).doTestSubtraction();
}
if ( genDipoles.empty() && factory()->initVerbose() ) {
// probably finite real contribution, but warn
generator()->log() << "\nWarning: No subtraction dipoles could be found for the process:\n";
generator()->log() << real->subProcess().legs[0]->PDGName() << " "
<< real->subProcess().legs[1]->PDGName() << " -> ";
for ( PDVector::const_iterator p = real->subProcess().legs.begin() + 2;
p != real->subProcess().legs.end(); ++p )
generator()->log() << (**p).PDGName() << " ";
generator()->log() << "\n" << flush;
generator()->log() << "Assuming finite tree-level O(alphaS) correction.\n";
}
dependent().resize(genDipoles.size());
copy(genDipoles.begin(),genDipoles.end(),dependent().begin());
if ( !factory()->reweighters().empty() ) {
for ( MEVector::const_iterator d = dependent().begin(); d != dependent().end(); ++d ) {
for ( vector<ReweightPtr>::const_iterator rw = factory()->reweighters().begin();
rw != factory()->reweighters().end(); ++rw )
(**d).addReweighter(*rw);
}
}
if ( !factory()->preweighters().empty() ) {
for ( MEVector::const_iterator d = dependent().begin(); d != dependent().end(); ++d ) {
for ( vector<ReweightPtr>::const_iterator rw = factory()->preweighters().begin();
rw != factory()->preweighters().end(); ++rw )
(**d).addPreweighter(*rw);
}
}
}
void SubtractedME::cloneRealME(const string& prefix) {
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theReal ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = theReal->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw Exception() << "SubtractedME: Matrix element " << pname.str() << " already existing." << Exception::runerror;
myRealEmissionME->cloneDependencies(pname.str());
theReal = myRealEmissionME;
}
head(theReal);
}
void SubtractedME::cloneDipoles(const string& prefix) {
MEVector dipMEs;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
Ptr<SubtractionDipole>::ptr cloned = dip->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << cloned->name();
if ( ! (generator()->preinitRegister(cloned,pname.str()) ) )
throw Exception() << "SubtractedME: Subtraction dipole " << pname.str() << " already existing." << Exception::runerror;
cloned->cloneDependencies(pname.str());
dipMEs.push_back(cloned);
}
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::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->showerApproximation(showerApproximation());
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->showerApproximation(showerApproximation());
dip->doVirtualShowerSubtraction();
}
}
void SubtractedME::doLoopSimSubtraction() {
theLoopSimSubtraction = 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->showerApproximation(showerApproximation());
dip->doLoopSimSubtraction();
}
}
void SubtractedME::setVetoScales(tSubProPtr) const {}
void SubtractedME::fillProjectors() {
if ( !virtualShowerSubtraction() && !loopSimSubtraction() )
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).matrixElement()->apply() ||
!(**d).kinematicsGenerated() )
continue;
- if ( (**d).willPassCuts() )
- lastXCombPtr()->projectors().insert((**d).cutWeight(),*d);
+ if ( (**d).willPassCuts() &&
+ (**d).lastMECrossSection()/picobarn != 0.0 ) {
+ lastXCombPtr()->projectors().insert(abs((**d).cutWeight()*(**d).lastMECrossSection()/picobarn),*d);
+ }
}
}
double SubtractedME::reweightHead(const vector<tStdXCombPtr>&) {
if ( showerApproximation() ) {
if ( realShowerSubtraction() )
return 1.;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
return 0.;
}
return 1.;
}
double SubtractedME::reweightDependent(tStdXCombPtr xc, const vector<tStdXCombPtr>& dep) {
if ( showerApproximation() ) {
if ( realShowerSubtraction() )
return 1.0;
if ( virtualShowerSubtraction() || loopSimSubtraction() ) {
if ( !lastXComb().lastProjector() )
return 0.0;
if ( xc != lastXComb().lastProjector() )
return 0.0;
double invPAlpha = 0.;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( !(**d).matrixElement()->apply() ||
!(**d).kinematicsGenerated() )
continue;
- if ( (**d).willPassCuts() ) {
- invPAlpha += (**d).cutWeight();
+ if ( (**d).willPassCuts() &&
+ (**d).lastMECrossSection()/picobarn != 0.0 ) {
+ invPAlpha += abs((**d).cutWeight()*(**d).lastMECrossSection()/picobarn);
}
}
- assert(invPAlpha != 0.0);
- double palpha = xc->cutWeight()/invPAlpha;
+ assert(invPAlpha != 0.0 && xc->cutWeight() != 0.0 && xc->lastMECrossSection()/picobarn != 0.0);
+ double palpha = abs((xc->cutWeight())*(xc->lastMECrossSection()/picobarn))/invPAlpha;
return 1./palpha;
}
}
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();
}
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator b = theBorns.begin();
b != theBorns.end(); ++b )
(**b).init();
if ( initVerbose() )
print(Repository::clog());
MEGroup::doinit();
}
void SubtractedME::doinitrun() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinitrun();
return;
}
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
for ( vector<Ptr<MatchboxMEBase>::ptr>::iterator b = theBorns.begin();
b != theBorns.end(); ++b )
(**b).initrun();
MEGroup::doinitrun();
}
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(factory()->subtractionData(),
factory()->subtractionPlotType(),
factory()->subtractionScatterPlot(),
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(factory()->subtractionData(),
factory()->subtractionPlotType(),
factory()->subtractionScatterPlot(),
b->first.first,
b->first.second,
b->first.second);
}
}
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 ( verbose() )
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::persistentOutput(PersistentOStream& os) const {
os << lower << bins;
}
void SubtractedME::SubtractionHistogram::persistentInput(PersistentIStream& is) {
is >> lower >> bins;
}
void SubtractedME::SubtractionHistogram::
dump(const std::string& prefix,
const int& plottype,
const bool& scatterplot,
const cPDVector& proc,
int i, int j) const {
bool bbmin = true;
double bmin = bins.begin()->first;
double bmax = bins.begin()->first;
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.0 ) {
if ( b != bins.begin() ){
out << bp->first;
if (bbmin){
bmin = bp->first;
bbmin = false;
}
}
else {
out << lower;
if (bbmin){
bmin = lower;
bbmin = false;
}
}
bmax = b->first;
out << " " << b->first
<< " " << b->second.first
<< " " << b->second.second
<< "\n" << flush;
}
}
double xmin = pow(10.0, floor(log10(bmin)));
double xmax = pow(10.0, ceil(log10(bmax)));
ofstream gpout((prefix+fname.str()+".gp").c_str());
gpout << "set terminal epslatex color solid\n"
<< "set output '" << fname.str() << "-plot.tex'\n"
<< "set format x '$10^{%T}$'\n"
<< "set logscale x\n"
<< "set xrange [" << xmin << ":" << xmax << "]\n";
if ( i != j ) {
gpout << "set xlabel '$\\sqrt{s_{" << i << j << "}}/{\\rm GeV}$'\n";
} else {
gpout << "set xlabel '$E_{" << i << "}/{\\rm GeV}$'\n";
}
if (plottype == 1){
gpout << "set size 0.5,0.6\n"
<< "set yrange [0:2]\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 '$";
}
else if (plottype == 2){
gpout << "set key left top Left reverse\n"
<< "set logscale y\n"
<< "set format y '$10^{%T}$'\n"
<< "set size 0.7,0.8\n"
<< "set yrange [1e-6:1e1]\n"
<< "set ylabel '$\\max\\left\\{\\left|\\mathcal{D}-\\mathcal{M}\\right|/\\left|\\mathcal{M}\\right|\\right\\}$'\n"
<< "unset bars\n";
gpout << "plot '";
if (scatterplot) gpout << fname.str() << "-scatter.dat' w points pt 7 ps 0.5 lc rgbcolor \"#00AACC\" not, \\\n'";
gpout << fname.str() << ".dat' u (($1+$2)/2.):4 w lines lw 4 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;
for ( vector<StdXCombPtr>::const_iterator d = xc->dependent().begin();
d != xc->dependent().end(); ++d ) {
if ( !(*d) )
continue;
if ( !(**d).matrixElement()->apply() )
continue;
if ( !(**d).willPassCuts() )
continue;
xcdip += (**d).lastCrossSection();
}
// want a real emission safely above the cut
if ( xc->cutWeight() < 1.0 )
return;
double delta;
if (factory()->subtractionPlotType() == 2) delta = abs(xcdip+xcme2)/abs(xcme2);
else delta = abs(xcdip)/abs(xcme2);
if ( theReal->phasespace() ) {
size_t i = lastSingularLimit()->first;
size_t j = lastSingularLimit()->second;
if ( i == j &&
softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,delta);
if ( factory()->subtractionScatterPlot() ){
ofstream outstream((fnamesSoftSubtraction[SoftSubtractionIndex(head()->mePartonData(),i)]).c_str(),ofstream::app);
outstream << meMomenta()[i].t()/GeV << " " << delta << "\n";
}
}
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,delta);
if ( factory()->subtractionScatterPlot() ){
ofstream outstream((fnamesCollinearSubtraction[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))]).c_str(),ofstream::app);
outstream << s << " " << delta << "\n";
}
}
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,delta);
if ( factory()->subtractionScatterPlot() ){
ofstream outstream((fnamesSoftSubtraction[SoftSubtractionIndex(head()->mePartonData(),i)]).c_str(),ofstream::app);
outstream << meMomenta()[i].t()/GeV << " " << delta << "\n";
}
}
}
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,delta);
if ( factory()->subtractionScatterPlot() ){
ofstream outstream((fnamesCollinearSubtraction[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))]).c_str(),ofstream::app);
outstream << s << " " << delta << "\n";
}
}
}
}
void SubtractedME::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theFactory << theBorns << theReal
<< collinearHistograms << softHistograms
<< fnamesSoftSubtraction
<< theRealShowerSubtraction << theVirtualShowerSubtraction
<< theLoopSimSubtraction
<< theSubProcessGroups;
}
void SubtractedME::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theFactory >> theBorns >> theReal
>> collinearHistograms >> softHistograms
>> fnamesSoftSubtraction
>> theRealShowerSubtraction >> theVirtualShowerSubtraction
>> theLoopSimSubtraction
>> theSubProcessGroups;
lastMatchboxXComb(theLastXComb);
}
void SubtractedME::Init() {
static ClassDocumentation<SubtractedME> documentation
("SubtractedME represents a subtracted real emission matrix element.");
}
// *** 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", "Herwig.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,1221 +1,1218 @@
// -*- C++ -*-
//
// SubtractionDipole.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SubtractionDipole class.
//
#include "SubtractionDipole.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
SubtractionDipole::SubtractionDipole()
: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
theIgnoreCuts(false),
theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
lastRealEmissionKey(realEmissionKey(cPDVector(),-1,-1,-1)),
lastUnderlyingBornKey(underlyingBornKey(cPDVector(),-1,-1)),
theBornEmitter(-1), theBornSpectator(-1),
theLastSubtractionScale(ZERO), theLastSplittingScale(ZERO),
theLastSubtractionPt(ZERO), theLastSplittingPt(ZERO),
theLastSubtractionZ(0.0), theLastSplittingZ(0.0),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
theLoopSimSubtraction(false), theRealEmissionScales(false),
theShowerHardScale(ZERO), theShowerScale(ZERO),
theIsInShowerPhasespace(false), theIsAboveCutoff(false) {}
SubtractionDipole::~SubtractionDipole() {}
double SubtractionDipole::alpha() const{
return factory()->alphaParameter();
}
void SubtractionDipole::clearBookkeeping() {
theRealEmitter = -1;
theRealEmission = -1;
theRealSpectator = -1;
theBornEmitter = -1;
theBornSpectator = -1;
theMergingMap.clear();
theSplittingMap.clear();
theIndexMap.clear();
theUnderlyingBornDiagrams.clear();
theRealEmissionDiagrams.clear();
}
void SubtractionDipole::setupBookkeeping(const map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo>& mergeInfo) {
theMergingMap.clear();
theSplittingMap.clear();
theUnderlyingBornDiagrams.clear();
theRealEmissionDiagrams.clear();
int xemitter = -1;
int xspectator = -1;
map<int,int> mergeLegs;
map<int,int> remapLegs;
map<int,int> realBornMap;
map<int,int> bornRealMap;
for ( map<Ptr<DiagramBase>::ptr,MergeInfo>::const_iterator mit = mergeInfo.begin();
mit != mergeInfo.end(); ++mit ) {
DiagramVector::const_iterator bd =
theUnderlyingBornME->diagrams().end();
map<int,int> theRemapLegs;
for ( DiagramVector::const_iterator b =
theUnderlyingBornME->diagrams().begin();
b != theUnderlyingBornME->diagrams().end(); ++b )
if ( mit->second.diagram->isSame(*b,theRemapLegs) ) {
bd = b; break;
}
// no underlying Born
if ( bd == theUnderlyingBornME->diagrams().end() )
continue;
if ( xemitter == -1 ) {
xemitter = mit->second.emitter;
mergeLegs = mit->second.mergeLegs;
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((*mit->first).partons(),realEmitter(),realEmission(),realSpectator());
UnderlyingBornKey bornKey = underlyingBornKey((**bd).partons(),xemitter,xspectator);
if ( theMergingMap.find(realKey) == theMergingMap.end() )
theMergingMap.insert(make_pair(realKey,make_pair(bornKey,realBornMap)));
RealEmissionInfo realInfo = make_pair(realKey,bornRealMap);
bool gotit = false;
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spIterator;
pair<spIterator,spIterator> range = theSplittingMap.equal_range(bornKey);
for ( ; range.first != range.second; ++range.first )
if ( range.first->second == realInfo ) {
gotit = true;
break;
}
if ( !gotit ) {
theSplittingMap.insert(make_pair(bornKey,realInfo));
theUnderlyingBornDiagrams[process(realKey)].push_back(*bd);
theRealEmissionDiagrams[process(bornKey)].push_back(mit->first);
}
}
if ( theSplittingMap.empty() )
return;
theIndexMap.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator s =
theSplittingMap.begin(); s != theSplittingMap.end(); ++s ) {
theIndexMap[process(s->first)] = make_pair(emitter(s->first),spectator(s->first));
}
}
void SubtractionDipole::subtractionBookkeeping() {
/*
if ( theMergingMap.empty() )
setupBookkeeping();
*/
assert(!theMergingMap.empty());
lastRealEmissionKey =
realEmissionKey(lastHeadXComb().mePartonData(),realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() ) {
theApply = false;
return;
}
theApply = true;
lastUnderlyingBornKey = k->second.first;
bornEmitter(emitter(lastUnderlyingBornKey));
bornSpectator(spectator(lastUnderlyingBornKey));
}
void SubtractionDipole::splittingBookkeeping() {
/*
if ( theMergingMap.empty() )
setupBookkeeping();
*/
assert(!theMergingMap.empty());
map<cPDVector,pair<int,int> >::const_iterator esit =
theIndexMap.find(lastHeadXComb().mePartonData());
if ( esit == theIndexMap.end() ) {
theApply = false;
return;
}
theApply = true;
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(lastHeadXComb().mePartonData(),bornEmitter(),bornSpectator());
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
assert(kr.first != kr.second);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
assert(lastRealEmissionInfo != kr.second);
lastRealEmissionKey = lastRealEmissionInfo->second.first;
realEmitter(emitter(lastRealEmissionKey));
realEmission(emission(lastRealEmissionKey));
realSpectator(spectator(lastRealEmissionKey));
}
StdXCombPtr SubtractionDipole::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allBins,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
return
realEmissionME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
StdXCombPtr SubtractionDipole::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
return
realEmissionME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
StdXCombPtr SubtractionDipole::makeBornXComb(tStdXCombPtr realXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();
lastRealEmissionKey =
realEmissionKey(proc,realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() )
return StdXCombPtr();
PartonPairVec pbs = realXC->pExtractor()->getPartons(realXC->maxEnergy(),
realXC->particles(),
*(realXC->cuts()));
DiagramVector bornDiags = underlyingBornDiagrams(proc);
assert(!bornDiags.empty());
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == bornDiags.front()->partons()[0] &&
ppit->second->parton() == bornDiags.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
return
underlyingBornME()->makeXComb(realXC,*ppit,bornDiags,this);
}
vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();
map<cPDVector,pair<int,int> >::const_iterator esit = theIndexMap.find(proc);
if ( esit == theIndexMap.end() )
return vector<StdXCombPtr>();
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(proc,bornEmitter(),bornSpectator());
if ( theSplittingMap.find(lastUnderlyingBornKey) == theSplittingMap.end() )
return vector<StdXCombPtr>();
PartonPairVec pbs = bornXC->pExtractor()->getPartons(bornXC->maxEnergy(),
bornXC->particles(),
*(bornXC->cuts()));
DiagramVector realDiags = realEmissionDiagrams(proc);
assert(!realDiags.empty());
vector<StdXCombPtr> res;
map<cPDVector,DiagramVector> realProcs;
for ( MEBase::DiagramVector::const_iterator d = realDiags.begin();
d != realDiags.end(); ++d ) {
realProcs[(**d).partons()].push_back(*d);
}
for ( map<cPDVector,DiagramVector>::const_iterator pr =
realProcs.begin(); pr != realProcs.end(); ++pr ) {
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == pr->second.front()->partons()[0] &&
ppit->second->parton() == pr->second.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
StdXCombPtr rxc =
realEmissionME()->makeXComb(bornXC,*ppit,pr->second,this);
res.push_back(rxc);
}
return res;
}
const MEBase::DiagramVector& SubtractionDipole::underlyingBornDiagrams(const cPDVector& real) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theUnderlyingBornDiagrams.find(real);
if (k == theUnderlyingBornDiagrams.end() )
return empty;
return k->second;
}
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()->setXComb(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagrams(dv);
}
Selector<const ColourLines *>
SubtractionDipole::colourGeometries(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->colourGeometries(diag) :
underlyingBornME()->colourGeometries(diag);
}
const ColourLines &
SubtractionDipole::selectColourGeometry(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->selectColourGeometry(diag) :
underlyingBornME()->selectColourGeometry(diag);
}
void SubtractionDipole::flushCaches() {
theUnderlyingBornME->flushCaches();
theRealEmissionME->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
}
void SubtractionDipole::setXComb(tStdXCombPtr xc) {
if ( !xc ) {
theApply = false;
return;
} else {
theApply = true;
}
lastMatchboxXComb(xc);
MEBase::setXComb(xc);
if ( splitting() ) {
realEmissionME()->setXComb(xc);
underlyingBornME()->setXComb(xc->head());
splittingBookkeeping();
} else {
realEmissionME()->setXComb(xc->head());
underlyingBornME()->setXComb(xc);
subtractionBookkeeping();
}
if ( !apply() )
return;
}
void SubtractionDipole::setKinematics() {
MEBase::setKinematics();
if ( splitting() )
realEmissionME()->setKinematics();
else
underlyingBornME()->setKinematics();
}
bool SubtractionDipole::generateKinematics(const double * r) {
if ( lastXCombPtr()->kinematicsGenerated() )
return true;
if ( splitting() ) {
if ( !generateRadiationKinematics(r) )
return false;
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()->lastXCombPtr()->setIncomingPartons();
underlyingBornME()->setScale();
assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
underlyingBornME()->lastXCombPtr()->setIncomingPartons();
// need to have the scale and x's available for checking shower phase space
if ( showerApproximation() &&
lastXCombPtr()->willPassCuts() )
showerApproximation()->getShowerVariables();
lastXCombPtr()->didGenerateKinematics();
return true;
}
int SubtractionDipole::nDim() const {
if ( !splitting() )
return underlyingBornME()->nDim();
return underlyingBornME()->nDim() + nDimRadiation();
}
void SubtractionDipole::clearKinematics() {
MEBase::clearKinematics();
if ( splitting() )
realEmissionME()->clearKinematics();
else
underlyingBornME()->clearKinematics();
}
void SubtractionDipole::tildeKinematics(Ptr<TildeKinematics>::tptr tk) {
theTildeKinematics = tk;
}
bool SubtractionDipole::generateTildeKinematics() {
assert(!splitting());
Ptr<TildeKinematics>::tptr kinematics = theTildeKinematics;
if ( showerApproximation() ) {
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
showerApproximation()->checkCutoff();
if ( showerApproximation()->showerTildeKinematics() &&
isAboveCutoff() &&
realShowerSubtraction() )
kinematics = showerApproximation()->showerTildeKinematics();
}
if ( !kinematics ) {
jacobian(0.0);
return false;
}
kinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());
if ( !kinematics->doMap() ) {
jacobian(0.0);
return false;
}
theLastSubtractionScale = kinematics->lastScale();
theLastSubtractionPt = kinematics->lastPt();
theLastSubtractionZ = kinematics->lastZ();
meMomenta().resize(lastHeadXComb().meMomenta().size() - 1);
assert(mergingMap().find(lastRealEmissionKey) != mergingMap().end());
map<int,int>& trans = theMergingMap[lastRealEmissionKey].second;
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == realEmitter() || k == realEmission() || k == realSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( kinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[bornEmitter()] =
const_cast<const TildeKinematics&>(*kinematics).bornEmitterMomentum();
meMomenta()[bornSpectator()] =
const_cast<const TildeKinematics&>(*kinematics).bornSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).hardProcessMass());
p->rescaleRho();
}
jacobian(realEmissionME()->lastXComb().jacobian());
logGenerateTildeKinematics();
return true;
}
void SubtractionDipole::invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr itk) {
theInvertedTildeKinematics = itk;
}
int SubtractionDipole::nDimRadiation() const {
return invertedTildeKinematics() ?
invertedTildeKinematics()->nDimRadiation() :
0;
}
bool SubtractionDipole::generateRadiationKinematics(const double * r) {
assert(splitting());
Ptr<InvertedTildeKinematics>::tptr kinematics = theInvertedTildeKinematics;
if ( showerApproximation() ) {
showerApproximation()->setBornXComb(lastHeadXCombPtr());
showerApproximation()->setRealXComb(lastXCombPtr());
showerApproximation()->setDipole(this);
if ( showerApproximation()->showerInvertedTildeKinematics() ) {
kinematics = showerApproximation()->showerInvertedTildeKinematics();
}
}
if ( !kinematics ) {
jacobian(0.0);
return false;
}
kinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());
if ( !kinematics->doMap(r) ) {
jacobian(0.0);
return false;
}
theLastSplittingScale = kinematics->lastScale();
theLastSplittingPt = kinematics->lastPt();
theLastSplittingZ = kinematics->lastZ();
meMomenta().resize(lastHeadXComb().meMomenta().size() + 1);
assert(splittingMap().find(lastUnderlyingBornKey) != splittingMap().end());
map<int,int>& trans = const_cast<map<int,int>&>(lastRealEmissionInfo->second.second);
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == bornEmitter() || k == bornSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( kinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[realEmitter()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realEmitterMomentum();
meMomenta()[realEmission()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realEmissionMomentum();
meMomenta()[realSpectator()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).hardProcessMass());
p->rescaleRho();
}
jacobian(underlyingBornME()->lastXComb().jacobian() *
kinematics->jacobian());
logGenerateRadiationKinematics(r);
return true;
}
void SubtractionDipole::ptCut(Energy cut) {
theInvertedTildeKinematics->ptCut(cut);
}
CrossSection SubtractionDipole::dSigHatDR(Energy2 factorizationScale) const {
double pdfweight = 1.;
double jac = jacobian();
if ( splitting() && jac == 0.0 ) {
lastMECrossSection(ZERO);
- lastME2(0.0);
return ZERO;
}
if ( factorizationScale == ZERO ) {
factorizationScale = underlyingBornME()->lastScale();
}
if ( havePDFWeight1() ) {
pdfweight *= realEmissionME()->pdf1(factorizationScale);
}
if ( havePDFWeight2() ) {
pdfweight *= realEmissionME()->pdf2(factorizationScale);
}
lastMEPDFWeight(pdfweight);
bool needTheDipole = true;
CrossSection shower = ZERO;
double lastThetaMu = 1.0;
double showerFactor = 1.;
if ( showerApproximation() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(const_cast<SubtractionDipole*>(this));
if ( !isAboveCutoff() ) {
showerApproximation()->wasBelowCutoff();
lastThetaMu = 0.0;
} else {
lastThetaMu = 1.0;
}
if ( lastThetaMu > 0.0 && isInShowerPhasespace() ) {
if ( realShowerSubtraction() )
shower = showerApproximation()->dSigHatDR()*lastThetaMu;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
shower = -showerApproximation()->dSigHatDR()*lastThetaMu;
if ( virtualShowerSubtraction() &&
isAboveCutoff() &&
showerApproximation()->showerTildeKinematics() ) {
// map shower to dipole kinematics; we are always above the
// cutoff in this case
showerFactor *=
showerApproximation()->showerTildeKinematics()->jacobianRatio();
}
shower *= showerFactor;
}
if ( realShowerSubtraction() && lastThetaMu == 1.0 )
needTheDipole = false;
if ( virtualShowerSubtraction() && lastThetaMu == 0.0 )
needTheDipole = false;
if ( factory()->loopSimCorrections() ||
factory()->meCorrectionsOnly() )
needTheDipole = false;
}
double xme2 = 0.0;
if ( needTheDipole )
xme2 = me2();
if ( factory()->loopSimCorrections() ||
factory()->meCorrectionsOnly() ) {
assert(showerApproximation());
xme2 = realEmissionME()->me2() * showerApproximation()->channelWeight();
double rws =
pow(underlyingBornME()->lastXComb().lastAlphaS()/
realEmissionME()->lastXComb().lastAlphaS(),
realEmissionME()->orderInAlphaS());
xme2 *= rws;
double rwe =
pow(underlyingBornME()->lastXComb().lastAlphaEM()/
realEmissionME()->lastXComb().lastAlphaEM(),
underlyingBornME()->orderInAlphaEW());
xme2 *= rwe;
}
if ( realShowerSubtraction() )
xme2 *= 1. - lastThetaMu;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
xme2 *= lastThetaMu;
double coupl = lastMECouplings();
coupl *= underlyingBornME()->lastXComb().lastAlphaS();
lastMECouplings(coupl);
- lastME2(xme2);
-
CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
if ( !showerApproximation() && xme2 != 0.0 ) {
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(theRealEmissionME->lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
}
lastMECrossSection(-res-shower);
logDSigHatDR(jac);
return lastMECrossSection();
}
void SubtractionDipole::print(ostream& os) const {
os << "--- SubtractionDipole setup ----------------------------------------------------\n";
os << " subtraction '" << name() << "'\n for real emission '"
<< theRealEmissionME->name() << "'\n using underlying Born '"
<< theUnderlyingBornME->name() << "'\n";
os << " tilde kinematics are '"
<< (theTildeKinematics ? theTildeKinematics->name() : "")
<< " '\n inverted tilde kinematics are '"
<< (theInvertedTildeKinematics ? theInvertedTildeKinematics->name() : "") << "'\n";
os << " the following subtraction mappings have been found:\n";
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator m =
theMergingMap.begin(); m != theMergingMap.end(); ++m ) {
os << " " << process(m->second.first)[0]->PDGName() << " "
<< process(m->second.first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->second.first).begin() + 2;
p != process(m->second.first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[" << emitter(m->second.first) << "," << spectator(m->second.first) << "] <=> ";
os << process(m->first)[0]->PDGName() << " "
<< process(m->first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->first).begin() + 2;
p != process(m->first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[(" << emitter(m->first) << "," << emission(m->first) << ")," << spectator(m->first) << "]\n"
<< " non-dipole momenta ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->second << " ";
}
os << ") <=> ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->first << " ";
}
os << ")\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractionDipole::printLastEvent(ostream& os) const {
os << "--- SubtractionDipole last event information -----------------------------------\n";
os << " for dipole '" << name() << "' applying ["
<< bornEmitter() << "," << bornSpectator() << "] <=> [("
<< realEmitter() << "," << realEmission() << ")," << realSpectator() << "]\n"
<< " evaluated the cross section/nb " << (lastMECrossSection()/nanobarn) << "\n"
<< " with subtraction parameters x[0] = " << subtractionParameters()[0]
<< " x[1] = " << subtractionParameters()[1] << "\n";
os << " the last real emission event was:\n";
realEmissionME()->printLastEvent(os);
os << " the last underlying Born event was:\n";
underlyingBornME()->printLastEvent(os);
os << "--- end SubtractionDipole last event information -------------------------------\n";
os << flush;
}
void SubtractionDipole::logME2() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated me2 using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "Born phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = bornxc->meMomenta().begin();
cPDVector::const_iterator dit = bornxc->mePartonData().begin();
for ( ; pit != bornxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << bornxc->lastX1() << " x2 = " << bornxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (bornxc->lastSHat()/GeV2) << "\n";
generator()->log() << "Real emission phase space point (in GeV):\n";
pit = realxc->meMomenta().begin();
dit = realxc->mePartonData().begin();
for ( ; pit != realxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << realxc->lastX1() << " x2 = " << realxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (realxc->lastSHat()/GeV2) << "\n";
}
void SubtractionDipole::logDSigHatDR(double effectiveJac) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated cross section using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n"
<< "Jacobian = " << jacobian()
<< " effective Jacobian = " << effectiveJac << "\n"
<< "Born sHat/GeV2 = " << (bornxc->lastSHat()/GeV2)
<< " real sHat/GeV2 = " << (realxc->lastSHat()/GeV2)
<< " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void SubtractionDipole::logGenerateTildeKinematics() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating tilde kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with real xcomb " << lastHeadXCombPtr() << " born xcomb "
<< lastXCombPtr() << "\n"
<< "from real emission phase space point:\n";
Lorentz5Momentum rSum;
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
size_t count = 0;
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
rSum -= *pr;
} else {
rSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (rSum/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n"
<< "with scale/GeV = " << (theLastSubtractionScale/GeV)
<< "and pt/GeV = " << (theLastSubtractionPt/GeV) << "\n";
generator()->log() << "generated tilde kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
count = 0;
Lorentz5Momentum bSum;
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
bSum -= *pr;
} else {
bSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (bSum/GeV) << "\n";
generator()->log() << "Jacobian = " << jacobian() << "\n" << flush;
}
void SubtractionDipole::logGenerateRadiationKinematics(const double * r) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating radiation kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with born xcomb " << lastHeadXCombPtr() << " real xcomb "
<< lastXCombPtr() << "\n"
<< "from random numbers:\n";
copy(r,r+nDimRadiation(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "and born phase space point:\n";
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n" << flush;
generator()->log() << "scales: scale/GeV = " << (theLastSplittingScale/GeV)
<< " pt/GeV = " << (theLastSplittingPt/GeV) << "\n" << flush;
generator()->log() << "generated real emission kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "Jacobian = "
<< jacobian() << " = "
<< underlyingBornME()->lastXComb().jacobian()
<< "|Born * "
<< invertedTildeKinematics()->jacobian()
<< "|Radiation\n" << flush;
}
void SubtractionDipole::doinit() {
MEBase::doinit();
if ( underlyingBornME() ) {
theUnderlyingBornME->init();
}
if ( realEmissionME() ) {
theRealEmissionME->init();
}
if ( tildeKinematics() ) {
theTildeKinematics->init();
}
if ( invertedTildeKinematics() ) {
theInvertedTildeKinematics->init();
}
if ( showerApproximation() ) {
theShowerApproximation->init();
}
for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
p != thePartners.end(); ++p ) {
(**p).init();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).init();
}
}
void SubtractionDipole::doinitrun() {
MEBase::doinitrun();
if ( underlyingBornME() ) {
theUnderlyingBornME->initrun();
}
if ( realEmissionME() ) {
theRealEmissionME->initrun();
}
if ( tildeKinematics() ) {
theTildeKinematics->initrun();
}
if ( invertedTildeKinematics() ) {
theInvertedTildeKinematics->initrun();
}
if ( showerApproximation() ) {
theShowerApproximation->initrun();
}
for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
p != thePartners.end(); ++p ) {
(**p).initrun();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).initrun();
}
}
void SubtractionDipole::cloneDependencies(const std::string& prefix) {
if ( underlyingBornME() ) {
Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
myUnderlyingBornME->cloneDependencies(pname.str());
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 Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
myRealEmissionME->cloneDependencies(pname.str());
realEmissionME(myRealEmissionME);
}
if ( tildeKinematics() ) {
Ptr<TildeKinematics>::ptr myTildeKinematics = tildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myTildeKinematics->name();
if ( ! (generator()->preinitRegister(myTildeKinematics,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
myTildeKinematics->dipole(this);
tildeKinematics(myTildeKinematics);
}
if ( invertedTildeKinematics() ) {
Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Inverted tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
myInvertedTildeKinematics->dipole(this);
invertedTildeKinematics(myInvertedTildeKinematics);
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror;
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
}
void SubtractionDipole::constructVertex(tSubProPtr sub) {
if ( splitting() )
realEmissionME()->constructVertex(sub);
else
underlyingBornME()->constructVertex(sub);
}
void SubtractionDipole::constructVertex(tSubProPtr sub, const ColourLines* cl) {
if ( splitting() )
realEmissionME()->constructVertex(sub,cl);
else
underlyingBornME()->constructVertex(sub,cl);
}
void SubtractionDipole::generateSubCollision(SubProcess & sub) {
if ( splitting() )
realEmissionME()->generateSubCollision(sub);
else
underlyingBornME()->generateSubCollision(sub);
}
void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theSplitting << theApply << theSubtractionTest
<< theIgnoreCuts << theRealEmissionME << theUnderlyingBornME
<< thePartners << theTildeKinematics << theInvertedTildeKinematics
<< theReweights << theRealEmitter << theRealEmission << theRealSpectator
<< theSubtractionParameters << theMergingMap << theSplittingMap
<< theIndexMap << theUnderlyingBornDiagrams << theRealEmissionDiagrams
<< lastRealEmissionKey << lastUnderlyingBornKey
<< theBornEmitter << theBornSpectator << ounit(theLastSubtractionScale,GeV)
<< ounit(theLastSplittingScale,GeV) << ounit(theLastSubtractionPt,GeV)
<< ounit(theLastSplittingPt,GeV) << theLastSubtractionZ
<< theLastSplittingZ << theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction
<< theLoopSimSubtraction << theRealEmissionScales << theFactory
<< ounit(theShowerHardScale,GeV) << ounit(theShowerScale,GeV)
<< theShowerParameters << theIsInShowerPhasespace << theIsAboveCutoff;
}
void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theSplitting >> theApply >> theSubtractionTest
>> theIgnoreCuts >> theRealEmissionME >> theUnderlyingBornME
>> thePartners >> theTildeKinematics >> theInvertedTildeKinematics
>> theReweights >> theRealEmitter >> theRealEmission >> theRealSpectator
>> theSubtractionParameters >> theMergingMap >> theSplittingMap
>> theIndexMap >> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
>> lastRealEmissionKey >> lastUnderlyingBornKey
>> theBornEmitter >> theBornSpectator >> iunit(theLastSubtractionScale,GeV)
>> iunit(theLastSplittingScale,GeV) >> iunit(theLastSubtractionPt,GeV)
>> iunit(theLastSplittingPt,GeV) >> theLastSubtractionZ
>> theLastSplittingZ >> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction
>> theLoopSimSubtraction >> theRealEmissionScales >> theFactory
>> iunit(theShowerHardScale,GeV) >> iunit(theShowerScale,GeV)
>> theShowerParameters >> theIsInShowerPhasespace >> theIsAboveCutoff;
lastMatchboxXComb(theLastXComb);
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
}
Ptr<MatchboxFactory>::tptr SubtractionDipole::factory() const {
return theFactory;
}
void SubtractionDipole::factory(Ptr<MatchboxFactory>::tptr f) {
theFactory = f;
}
void SubtractionDipole::Init() {
static ClassDocumentation<SubtractionDipole> documentation
("SubtractionDipole represents a dipole subtraction "
"term in the formalism of Catani and Seymour.");
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<SubtractionDipole,MEBase>
describeSubtractionDipole("Herwig::SubtractionDipole", "Herwig.so");
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 4:26 PM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805161
Default Alt Text
(68 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment