diff --git a/Shower/Dipole/Merging/Merger.cc b/Shower/Dipole/Merging/Merger.cc
--- a/Shower/Dipole/Merging/Merger.cc
+++ b/Shower/Dipole/Merging/Merger.cc
@@ -1,1590 +1,1586 @@
// -*- C++ -*-
//
// Merger.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Merger class.
//
#include "Merger.h"
#include "Node.h"
#include "MergingFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "Herwig/Shower/Dipole/DipoleShowerHandler.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/Constants.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "fastjet/ClusterSequence.hh"
using namespace Herwig;
Merger::Merger()
: MergerBase(),
Unlopsweights(true),
theCMWScheme(true),
isUnitarized(true),
isNLOUnitarized(true),
defMERegionByJetAlg(false),
theOpenInitialStateZ(false),
theChooseHistory(0),
theN0(0),
theOnlyN(-1),
theCurrentMaxLegs(-1),
xiRenME(1.),
xiFacME(1.),
xiRenSh(1.),
xiFacSh(1.),
xiQSh(1.),
theGamma(1.),
ee_ycut(-1),
pp_dcut(-1),
theSmearing(0.),
therenormscale(-1.*GeV),
theIRSafePT(1.*GeV),
theMergePt(4.*GeV),
theCentralMergePt(4.*GeV),
theMergingJetFinder(),
theLargeNBasis(),
theDipoleShowerHandler(),
theTreeFactory()
{
FFLTK=new_ptr(FFLightTildeKinematics());
FILTK=new_ptr(FILightTildeKinematics());
IFLTK=new_ptr(IFLightTildeKinematics());
IILTK=new_ptr(IILightTildeKinematics());
FFMTK=new_ptr(FFMassiveTildeKinematics());
FIMTK=new_ptr(FIMassiveTildeKinematics());
IFMTK=new_ptr(IFMassiveTildeKinematics());
}
Merger::~Merger() {}
IBPtr Merger::clone() const {
return new_ptr(*this);
}
IBPtr Merger::fullclone() const {
return new_ptr(*this);
}
pair<PVector,PVector> getInOut(NPtr Node){
PVector incoming;
for(auto const & i : {0,1})
incoming.push_back(Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]));
PVector outgoing;
for (size_t i=2;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
outgoing.push_back(p);
}
return make_pair(incoming,outgoing);
}
CrossSection Merger::MergingDSigDRBornGamma(NPtr Node){
double weightCL=0.;
weight=1.;
int pjs=-1;
if (!Node->children().empty()) {
auto const inoutpair=getInOut(Node);
NPtr rndCh= Node->randomChild(); // Here we make sure that clustering and splitting are invertible
if(!matrixElementRegion(inoutpair.first,inoutpair.second, rndCh->pT(), theMergePt))weight*= 0.;
}
NPtr Born= Node-> getHistory(true,xiQSh);
NPtr CLNode;
NPtr BornCL;
if( !Node->children().empty()){
if (UseRandom::rnd()<.5){
CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
weight*=inhist?(-2.*int(Node->children().size())):0.;
projected=true;
pjs=1;
weightCL=2.*int(Node->children().size());
BornCL=CLNode-> getHistory(false,xiQSh);
}else{
weight=2.;projected=false;
pjs=0;
}
}else{
weight=1.;
projected=false;
pjs=0;
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->legsize()-pjs))
return ZERO;
if(!Node->allAbove(mergePt()-0.1*GeV))weight=0.;
if(CLNode&&!CLNode->allAbove(mergePt()-0.1*GeV))weightCL=0.;
if (!Born->xcomb()->willPassCuts()){
if (Born->parent()) {
//cout<<"\n parent not pass";
}
return ZERO;
}
CrossSection res=ZERO;
bool maxMulti=Node->legsize() == int(maxLegsLO());
if(weight!=0.){
Energy startscale=CKKW_StartScale(Born);
fillHistory( startscale, Born, Node);
Node->runningPt(fillProjector(pjs));
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.&&weightCL==0.)return ZERO;
res+=weight*TreedSigDR(startscale,Node,(!maxMulti&&!projected)?theGamma:1.);
}
if(CLNode&&theGamma!=1.){
Energy startscale=CKKW_StartScale(BornCL);
fillHistory( startscale, BornCL, CLNode);
Node->runningPt(fillProjector(pjs));
weightCL*=history.back().weight*alphaReweight()*pdfReweight();
CrossSection diff=ZERO;
Node->nodeME()->factory()->setAlphaParameter(1.);
diff-=weightCL*CLNode->dipole()->dSigHatDR(sqr(startscale*xiFacME));
Node->nodeME()->factory()->setAlphaParameter(theGamma);
string alp=(CLNode->dipole()->aboveAlpha()?"Above":"Below");
diff+=weightCL*CLNode->dipole()->dSigHatDR(sqr(startscale*xiFacME));
Node->nodeME()->factory()->setAlphaParameter(1.);
res+=diff;
}
return res;
}
CrossSection Merger::MergingDSigDRBornStandard(NPtr Node){
int pjs=-1;
// Check if phase space poing is in ME region
if (!Node->children().empty()) {
auto const & inoutpair=getInOut(Node);
NPtr rndCh= Node->randomChild(); // Here we make sure that clustering and splitting are invertible
if(!matrixElementRegion(inoutpair.first,inoutpair.second, rndCh->pT(), theMergePt))return ZERO;
}
// get the history for the process
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
// at least one history step -> unitarisation
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;
pjs=1;
}else{
weight=2.;projected=false;
pjs=0;
}
}else{
// no ordered history -> no projection
weight=1.;projected=false;
pjs=0;
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->legsize()-(projected?1:0)))
return ZERO;
assert(Node->allAbove(mergePt()-0.1*GeV));
if (!Born->xcomb()->willPassCuts())return ZERO;
// calculate the staring scale
Energy startscale=CKKW_StartScale(Born);
// fill history with caluclation of sudakov supression
fillHistory( startscale, Born, Node);
// fill the projector, the argument gets set to current scale
Node->runningPt(fillProjector(pjs));
// the weight has three components to get shower history weight
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
//calculate the cross section
return weight*TreedSigDR(startscale,Node,1.);
}
CrossSection Merger::MergingDSigDRVirtualStandard(NPtr Node ){
int pjs=-1;
// Check if phase space poing is in ME region
if (!Node->children().empty()) {
auto inoutpair=getInOut(Node);
NPtr rndCh= Node->randomChild(); // Here we make sure that clustering and splitting are invertible
if(!matrixElementRegion(inoutpair.first,inoutpair.second, rndCh->pT(), theMergePt))return ZERO;
}
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;
pjs=1;
}else{
weight=2.;projected=false;
pjs=0;
}
}else{
weight=1.;projected=false;
pjs=0;
}
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
fillHistory( startscale, Born, Node);
Node->runningPt(fillProjector(pjs));
double ww1=history.back().weight;
double ww2=alphaReweight();
double ww3=pdfReweight();
weight*=ww1*ww2*ww3;
if(weight==0.)return ZERO;
CrossSection matrixElement=LoopdSigDR(startscale,Node);
CrossSection Bornweight=Node->nodeME()->dSigHatDRB();
double w1=-sumpdfReweightUnlops();
double w2=-sumalphaReweightUnlops();
double w3=-sumfillHistoryUnlops();
CrossSection unlopsweight =(w1+w2+w3)
*Bornweight
*SM().alphaS()/(2.*ThePEG::Constants::pi);
if (Node->legsize()==5&&Debug::level > 2) {
Energy minPT=1000*GeV;
for(auto const & no :Node->children() )minPT=min(no->pT(),minPT);
cout<<"\nVIRT "<<minPT/GeV<<" "<<weight<<" "<<w1;
cout<<" "<<w2;
cout<<" "<<w3;
cout<<" "<<(matrixElement/nanobarn)<<" "<<ww1<<" "<<ww2<<" "<<ww3<<" "<<Born->pT()/GeV<<" "<<Born->nodeME()->mePartonData()[3]->mass()/GeV<<" "<<(Bornweight*SM().alphaS()/(2.*ThePEG::Constants::pi)/nanobarn);
}
return weight*
as(startscale*xiRenSh)/SM().alphaS()*
(matrixElement+unlopsweight);
}
CrossSection Merger::MergingDSigDRRealStandard(NPtr Node){
bool allAbove=Node->allAbove(mergePt());
//TODO: VW Abgas Skandal
if(!Node->allAbove((Debug::level > 2?0.01:1.)*theIRSafePT))return ZERO;
if (allAbove)return MergingDSigDRRealAllAbove(Node);
if (UseRandom::rnd()<.5)
return 2.*MergingDSigDRRealBelowSubReal( Node );
return 2.*MergingDSigDRRealBelowSubInt( Node);
}
CrossSection Merger::MergingDSigDRRealAllAbove(NPtr Node){
int pjs=-1;
if (Node->children().empty()) {
throw Exception()
<< "Real emission contribution without underlying born."
<< "These are finite contibutions already handled in LO merging."
<< Exception::abortnow;
}
//If all dipoles pts are above, we cluster to this dipole.
NPtr CLNode= Node->randomChild();
// Check if phase space poing is in ME region--> else rm PSP
if (!CLNode->children().empty()) {
auto inoutpair=getInOut(CLNode);
NPtr rndCh= CLNode->randomChild(); // Here we make sure that clustering and splitting are invertible
if(!matrixElementRegion(inoutpair.first,inoutpair.second, rndCh->pT(), theMergePt))return ZERO;
}
// first find the history for the acctual Node
NPtr Born= Node-> getHistory(true,xiQSh);
// check if the randomly choosen CLNode is part of the history.
// If CLNode is not part of the history, dont calculate the Real contribution
// else multiply the real contribution with N (number of children).
// this returns the sudakov suppression according to the clustering of the born parts.
bool inhist=CLNode->isInHistoryOf(Born);
// get the history for the clustered Node.
Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true;
pjs=2;
}
else{
weight= 2.; projected=false;
pjs=1;
}
}else{
weight=1.;
projected=false;
pjs=1;
}
if (!CLNode->allAbove(mergePt()))return ZERO;
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
fillHistory( startscale, Born, CLNode);
Node->runningPt(fillProjector(pjs));
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
CrossSection me=(inhist?TreedSigDR(startscale,Node):ZERO);
CrossSection dip=CLNode->calcDip(startscale*xiFacME);
CrossSection res= weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*(me -dip);
if (Node->legsize()==6&&Debug::level > 2) {
Energy minPT=1000*GeV;
for(auto const & no :Node->children() )minPT=min(no->pT(),minPT);
cout<<"\nRAA "<<minPT/GeV<<" "<<weight<<" "<<(me-dip)/nanobarn<< " "<<me/nanobarn<<" "<<dip/nanobarn;
}
// cout<<"\nCLNode->dipole()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn "<<CLNode->dipole()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
return res;
}
CrossSection Merger::MergingDSigDRRealBelowSubReal(NPtr Node){
int pjs=-1;
NPtr HistNode=Node->randomChild();
if (!HistNode->children().empty()) {
auto inoutpair=getInOut(HistNode);
NPtr rndCh= HistNode->randomChild(); // Here we make sure that clustering and splitting are invertible
if(!matrixElementRegion(inoutpair.first,inoutpair.second, rndCh->pT(), theMergePt))return ZERO;
}
NPtr Born=HistNode-> getHistory(false,xiQSh);
if( Born!= HistNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true;
pjs=1;}
else{
weight= 2.; projected=false;
pjs=0;}
}else{
weight=1.; projected=false;
pjs=0;
}
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
fillHistory( startscale, Born, HistNode);
Node->runningPt(fillProjector(pjs));
weight*=history.back().weight*alphaReweight()*pdfReweight();
//cout<<"\n"<<history.back().weight<<" "<<alphaReweight()<<" "<<pdfReweight();
if(weight==0.)return ZERO;
CrossSection sumPS=ZERO;
for(auto const & child : Node->children()){
if (child->allAbove(mergePt())){
if((child)->pT()>mergePt()/3.)//TODO: this is a dynamical cutoff(see below)
sumPS+=child->calcPs(startscale*xiFacME);
else
sumPS+=child->calcDip(startscale*xiFacME);
}else{
assert(child->pT()>mergePt());
}
}
CrossSection me=TreedSigDR(startscale,Node);
if (Node->legsize()==6&&Debug::level > 2) {
Energy minPT=1000*GeV;
for(auto const & no :Node->children() )minPT=min(no->pT(),minPT);
cout<<"\nRBSR "<<minPT/GeV<< " " <<weight<<" "<<(me-sumPS)/nanobarn<<" "<<me/nanobarn<<" "<<sumPS/nanobarn;
}
//Here we subtract the PS (and below the dynamical cutoff the Dip)
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(me-sumPS);
}
CrossSection Merger::MergingDSigDRRealBelowSubInt(NPtr Node){
int pjs=-1;
if(Node->children().empty())return ZERO;
NPtr CLNode= Node->randomChild();
if(CLNode->pT()<mergePt()/3.)return ZERO;//TODO: this is a dynamical cutoff(see above)
if (!CLNode->children().empty()) {
auto inoutpair=getInOut(CLNode);
NPtr rndCh= CLNode->randomChild(); // Here we make sure that clustering and splitting are invertible
if(!matrixElementRegion(inoutpair.first,inoutpair.second, rndCh->pT(), theMergePt))return ZERO;
}
NPtr Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true;
pjs=2;}
else{
weight= 2.; projected=false;
pjs=1;}
}else{
weight=1.; projected=false;
pjs=1;
}
if (!CLNode->allAbove(mergePt()))return ZERO;
if (!Born->xcomb()->willPassCuts())return ZERO;
Energy startscale=CKKW_StartScale(Born);
fillHistory( startscale, Born, CLNode);
Node->runningPt(fillProjector(pjs));
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return ZERO;
pair<CrossSection,CrossSection> DipAndPs=
CLNode->calcDipandPS(startscale*xiFacME);
if (Node->legsize()==6&&Debug::level > 2) {
Energy minPT=1000*GeV;
for(auto const & no :Node->children() )minPT=min(no->pT(),minPT);
cout<<"\nRBSI "<<CLNode->pT()/GeV<<" "<<weight<<" "<<(DipAndPs.second-DipAndPs.first)/nanobarn<<" "<<DipAndPs.second/nanobarn<<" "<<DipAndPs.first/nanobarn;
}
//Here we add the PS and subtrac the Dip (only above the dynamical cutoff)
return weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*(DipAndPs.second-DipAndPs.first);
}
CrossSection Merger::TreedSigDR(Energy startscale,NPtr Node,double diffAlpha){
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
CrossSection res=DeepHead->nodeME()->dSigHatDRB();
if (diffAlpha!=1.) {
res+=DeepHead->nodeME()->dSigHatDRAlphaDiff(diffAlpha);
}
renormscale(0.0*GeV);
if(std::isnan(double(res/nanobarn))){cout<<"\nTreedSigDR is nan";res=ZERO;};
return res;
}
CrossSection Merger::LoopdSigDR(Energy startscale,NPtr Node){
// The deephead should be calculated here.
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
DeepHead->nodeME()->doOneLoopNoBorn();
CrossSection res=DeepHead->nodeME()->dSigHatDRV()+DeepHead->nodeME()->dSigHatDRI();
DeepHead->nodeME()->noOneLoopNoBorn();
renormscale(0.0*GeV);
return res;
}
Energy Merger::fillProjector(int pjs){
// in the shower handler the scale is multiplied by xiQSh
// so here we need to devide this
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
if(pjs == 0){
return (history.size()==1?1.:(1./xiQShfactor))*history.back().scale;
}
for(auto const & hs : history)
if (isProjectorStage(hs.node,pjs)&&pjs != 0){
history.begin()->node->deepHead()->xcomb()->lastProjector(hs.node->xcomb());
return (hs.node==history[0].node?1.:(1./xiQShfactor))*hs.scale;
}
throw Exception() << "Could not fill projector."<< Exception::abortnow;
return ZERO;
}
double Merger::pdfReweight(){
double res=1.;
double facfactor=(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
for(int side : {0,1}){
if(history[0].node->xcomb()->mePartonData()[side]->coloured()){
for (auto const & hs : history){
//pdf-ratio only to the last step
if(!hs.node->parent())continue;
if (hs.node==history.back().node)continue;
if(!hs.node->dipole()){
cout<<"\nthis should not happen";
return 0.;
}
res *= pdfratio(hs.node, facfactor*hs.scale,xiFacSh*(hs.node->pT()), side);
facfactor=xiFacSh;
}
res*=pdfratio(history.back().node,facfactor*history.back().scale ,history[0].scale*xiFacME, side);
}
}
if (std::isnan(res))cout<<"\nWarning:pdfReweight is nan.";
return res;
}
double Merger::alphaReweight(){
double res=1.;
Energy Q_R=(history[0].node->legsize()==N0()?xiRenME:xiRenSh)*history[0].scale;
res *= pow(as(Q_R) / SM().alphaS(), history[0].node->nodeME()->orderInAlphaS());
res *= pow(SM().alphaEMME(history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED())/ SM().alphaEMMZ(), history[0].node->nodeME()->orderInAlphaEW());
if (!(history[0].node->children().empty())){
res *=pow((theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(Q_R))*as(Q_R))/2./Constants::pi):1.),int(history[0].node->legsize()-N0()));
}
for (auto const & hs : history)
if (hs.node!=history.back().node){
Energy q_i=xiRenSh* hs.node->pT();
res *= as(q_i)/ SM().alphaS()
*(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(q_i))*as(q_i))/2./Constants::pi):1.);
if (std::isnan(res))cout<<"\nWarning:alphaReweight is nan:q_i "<<q_i/GeV<<" "<<Nf(q_i);
}
if (std::isnan(res))cout<<"\nWarning:alphaReweight is nan.";
return res;
}
void Merger::fillHistory(Energy scale, NPtr Begin, NPtr EndNode){
history.clear();
double sudakov0_n=1.;
history.push_back(HistoryStep(Begin,sudakov0_n,scale));
double xiQShfactor=history.begin()->node->legsize()==N0()?xiQSh:1.;
scale*=xiQShfactor;
if (Begin->parent()||!isUnitarized) {
while (Begin->parent() && (Begin != EndNode)) {
if (!dosudakov(Begin,scale, Begin->pT(),sudakov0_n)){
history.push_back(HistoryStep(Begin->parent(),0.,scale));
}
scale=Begin->pT();
if (std::isnan(sudakov0_n))cout<<"\nWarning:sudakov0_n is nan.";
history.push_back(HistoryStep(Begin->parent(),sudakov0_n,Begin->pT()));
Begin = Begin->parent();
}
Energy notunirunning=scale;
if (!isUnitarized&&N()+N0() > int(Begin->deepHead()->legsize())) {
if (!dosudakov(Begin,notunirunning,mergePt(),sudakov0_n)){
history.back().weight=0.;
}else{
history.back().weight=sudakov0_n;
}
}
}
if( history.size()==1)scale/=xiQSh;
}
double Merger::sumpdfReweightUnlops(){
double res=0.;
Energy beam1Scale=history[0].scale*(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
Energy beam2Scale=history[0].scale*(history[0].node->legsize()==N0()?xiFacME:xiFacSh);
for (auto const & hs : history){
//pdf expansion only to the last step
if(!hs.node->parent())continue;
if(hs.node->xcomb()->mePartonData()[0]->coloured()&&hs.node->nodeME()->lastX1()>0.99)return 0.;
if(hs.node->xcomb()->mePartonData()[1]->coloured()&&hs.node->nodeME()->lastX2()>0.99)return 0.;
if(hs.node->nodeME()->lastX1()<0.00001)return 0.;
if(hs.node->nodeME()->lastX2()<0.00001)return 0.;
if (hs.node->dipole()->bornEmitter() == 0 ){
res +=pdfUnlops(hs.node,0,
beam1Scale,
(hs.node->pT()),
hs.node->nodeME()->lastX1(),
Nf(history[0].scale),
history[0].scale);
beam1Scale=(hs.node->pT())*xiFacSh;
}
if (hs.node->dipole()->bornEmitter() == 1 ){
res +=pdfUnlops(hs.node,1,
beam2Scale,
(hs.node->pT()),
hs.node->nodeME()->lastX2(),
Nf(history[0].scale),
history[0].scale);
beam2Scale=(hs.node->pT())*xiFacSh;
}
if (hs.node->dipole()->bornSpectator() == 0 &&hs.node->dipole()->bornEmitter() >1){//
res +=pdfUnlops(hs.node,0,
beam1Scale,
(hs.node->pT()),
hs.node->nodeME()->lastX1(),
Nf(history[0].scale),
history[0].scale);
//pdfratio(hs.node, beam1Scale, sqrt(hs.node->pT()), 1);
beam1Scale=(hs.node->pT())*xiFacSh;
}
if (hs.node->dipole()->bornSpectator() == 1 &&hs.node->dipole()->bornEmitter() >1){//
res +=pdfUnlops(hs.node,1,
beam2Scale,
(hs.node->pT()),
hs.node->nodeME()->lastX2(),
Nf(history[0].scale),
history[0].scale);
//pdfratio(hs.node, beam2Scale , sqrt(hs.node->pT()), 2);
beam2Scale=(hs.node->pT())*xiFacSh;
}
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[0]->coloured()){
res +=pdfUnlops(history.back().node,0,
beam1Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX1(),
Nf(history[0].scale),
history[0].scale);
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[1]->coloured()) {
res +=pdfUnlops(history.back().node,1,
beam2Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX2(),
Nf(history[0].scale),
history[0].scale);
}
return res;
}
double Merger::pdfUnlops(NPtr node ,int side,Energy running, Energy next,double x,int nlp,Energy fixedScale) {
tcPDPtr particle,parton;
tcPDFPtr pdf;
if (side==0) {
particle = node->nodeME()->lastParticles().first->dataPtr();
parton=node->nodeME()->lastPartons().first->dataPtr();
pdf =node->xcomb()->partonBins().first->pdf();
}else{
assert(side==1);
particle = node->nodeME()->lastParticles().second->dataPtr();
parton=node->nodeME()->lastPartons().second->dataPtr();
pdf =node->xcomb()->partonBins().second->pdf();
}
//copied from PKOperator
double res=0.;
int number=10;
for (int nr=0;nr<number;nr++){
double restmp=0;
using namespace RandomHelpers;
double r = UseRandom::rnd();
double eps = 1e-3;
pair<double,double> zw =
generate((piecewise(),
flat(0.0,x),
match(inverse(0.0,x,1.0) + inverse(1.0+eps,x,1.0))),r);
double z = zw.first;
double mapz = zw.second;
double PDFxparton=pdf->xfx(particle,parton,sqr(fixedScale),x)/x;
double CA = SM().Nc();
double CF = (SM().Nc()*SM().Nc()-1.0)/(2.*SM().Nc());
if (abs(parton->id()) < 7) {
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
double PDFxByzparton=pdf->xfx(particle,parton,sqr(fixedScale),x/z)*z/x;
assert(abs(parton->id()) < 7);
restmp += CF*(3./2.+2.*log(1.-x)) * PDFxparton;
if ( z > x ) {
restmp+= 0.5 * ( sqr(z) + sqr(1.-z) ) * PDFxByzgluon / z;
restmp += CF*2.*(PDFxByzparton - z*PDFxparton)/(z*(1.-z));
restmp -= CF*PDFxByzparton * (1.+z)/z;
}
}else{
assert(parton->id() == ParticleID::g);
double PDFxByzgluon=pdf->xfx(particle,getParticleData(ParticleID::g),sqr(fixedScale),x/z)*z/x;
if ( z > x ){
double factor = CF * ( 1. + sqr(1.-z) ) / sqr(z);
for ( int f = -nlp; f <= nlp; ++f ) {
if ( f == 0 )
continue;
restmp += pdf->xfx(particle,getParticleData(f),sqr(fixedScale),x/z)*z/x*factor;
}
}
restmp += ( (11./6.) * CA - (1./3.)*Nf(history[0].scale) + 2.*CA*log(1.-x) ) *PDFxparton;
if ( z > x ) {
restmp += 2. * CA * ( PDFxByzgluon - z*PDFxparton ) / (z*(1.-z));
restmp += 2.* CA *( (1.-z)/z - 1. + z*(1.-z) ) * PDFxByzgluon / z;
}
}
if (PDFxparton<1e-8)restmp= 0.;
res+=1*restmp*log(sqr(running/next))/PDFxparton*mapz;
}
return res/number;
}
double Merger::sumalphaReweightUnlops(){
double res=0.;
if (!(history[0].node->children().empty())){
res +=alphasUnlops(history[0].scale*xiRenSh,
history[0].scale*xiRenME)*int(history[0].node->legsize()-N0());
}
// dsig is calculated with fixed alpha_s
for (auto const & hs : history){
//expansion only to the last step
if(!hs.node->parent())continue;
res +=alphasUnlops(hs.node->pT()*xiRenSh ,history[0].scale);
}
return res;
}
double Merger::sumfillHistoryUnlops(){
double res=0.;
double xiQShfactor=history[0].node->legsize()==N0()?xiQSh:1.;
for (auto const & hs : history){
if(!hs.node->parent())continue;
doUNLOPS(hs.node,(hs.node == history[0].node?xiQShfactor:1.)*hs.scale , hs.node->pT() , history[0].scale, res);
}
return res;
}
Ptr<MergingFactory>::ptr Merger::treefactory(){return theTreeFactory;}
void Merger::doinit(){
if (!DSH()->hardScaleIsMuF()) {
throw Exception()
<< "Merger: Merging is currently only sensible if we are using the hardScale as MuF."
<< Exception::abortnow;
}
}
CrossSection Merger::MergingDSigDR() {
history.clear();
if (theFirstNodeMap.find(theCurrentME)==theFirstNodeMap.end()) {
throw Exception()
<< "Merger: The current matrix element could not be found."
<< Exception::abortnow;
}
NPtr Node = theFirstNodeMap[theCurrentME];
-
xiRenME=theCurrentME->renormalizationScaleFactor();
xiFacME=theCurrentME->factorizationScaleFactor();
xiRenSh=DSH()->renormalizationScaleFactor();
xiFacSh=DSH()->factorizationScaleFactor();
xiQSh=DSH()->hardScaleFactor();
-
DSH()->setCurrentHandler();
-
DSH()->eventHandler(generator()->eventHandler());
-
-
+
CrossSection res=ZERO;
if(Node->deepHead()->subtractedReal()){
res=MergingDSigDRRealStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(Node->deepHead()->virtualContribution()){
res=MergingDSigDRVirtualStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(theGamma!=1.){
res=MergingDSigDRBornGamma(Node);
theCurrentMaxLegs=maxLegsLO();
}else{
res=MergingDSigDRBornStandard(Node);
theCurrentMaxLegs=maxLegsLO();
}
theCurrentME->lastXCombPtr()->lastCentralScale(sqr(Node->runningPt()));
theCurrentME->lastXCombPtr()->lastShowerScale(sqr(Node->runningPt()));
if(theCurrentME->lastXCombPtr()->lastProjector()){
theCurrentME->lastXCombPtr()->lastProjector()->lastCentralScale(sqr(Node->runningPt()));
theCurrentME->lastXCombPtr()->lastProjector()->lastShowerScale(sqr(Node->runningPt()));
}
renormscale(0.0*GeV);
if (res == ZERO){
history.clear();
return res;
}
cleanup(Node);
DSH()->eventRecord().clear();
theCurrentME->lastXCombPtr()->subProcess(SubProPtr());
history.clear();
if(std::isnan(double(res/nanobarn))|| !std::isfinite(double(res/nanobarn))){
cout<<"Warning merger weight is "<<res/nanobarn<<" -> setting to 0";
return ZERO;
}
return res;
}
void Merger::CKKW_PrepareSudakov(NPtr Born,Energy running){
//cleanup(Born);
tSubProPtr sub = Born->xcomb()->construct();
DSH()->resetPDFs(make_pair(Born->xcomb()->partonBins().first->pdf(),
Born->xcomb()->partonBins().second->pdf()),
Born->xcomb()->partonBins());
DSH()->setCurrentHandler();
DSH()->currentHandler()->generator()->currentEventHandler(Born->deepHead()->xcomb()->eventHandlerPtr());
DSH()->currentHandler()->remnantDecayer()->setHadronContent(Born->deepHead()->xcomb()->lastParticles());
DSH()->eventRecord().clear();
DSH()->eventRecord().slimprepare(sub, dynamic_ptr_cast<tStdXCombPtr>(Born->xcomb()), DSH()->pdfs(), Born->deepHead()->xcomb()->lastParticles());
DSH()->hardScales(sqr(running));
}
Energy Merger::CKKW_StartScale(NPtr Born){
Energy res=generator()->maximumCMEnergy();;
if(!Born->children().empty()){
for (int i=0;i<Born->legsize();i++){
if (!Born->nodeME()->mePartonData()[i]->coloured())continue;
for (int j=i+1;j<Born->legsize();j++){
if (i==j||!Born->nodeME()->mePartonData()[j]->coloured())continue;
res= min(res,sqrt(2.*Born->nodeME()->lastMEMomenta()[i]*Born->nodeME()->lastMEMomenta()[j]));
}
}
}else{
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res= sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale());
}
Born->nodeME()->factory()->scaleChoice()->setXComb(Born->xcomb());
res=max(res, sqrt(Born->nodeME()->factory()->scaleChoice()->renormalizationScale()));
return res;
}
double Merger::alphasUnlops( Energy next,Energy fixedScale) {
double betaZero = (11./6.)*SM().Nc() - (1./3.)*Nf(history[0].scale);
return (betaZero*log(sqr(fixedScale/next)))+
(theCMWScheme?(((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*Nf(history[0].scale)))):0.);
}
double Merger::pdfratio(NPtr Born,Energy nominator_scale, Energy denominator_scale,int side){
StdXCombPtr bXc = Born->xcomb();
if(!bXc->mePartonData()[side]->coloured())
throw Exception()
<< "Merger: pdf-ratio required for non-coloured particle."
<< Exception::abortnow;
double from,to;
if (side==0){
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf1(sqr(denominator_scale));
to =Born->nodeME()->pdf1(sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
else{
if (denominator_scale==nominator_scale) {
return 1.;
}
from=Born->nodeME()->pdf2(sqr(denominator_scale));
to=Born->nodeME()->pdf2( sqr( nominator_scale ));
if ( (to < 1e-8||from < 1e-8)&&(to/from>10000000.)){cout<<"\npdfratio to="<<to<<" from="<<from;return 0.;}
}
return to/from;
}
bool Merger::dosudakov(NPtr Born,Energy running, Energy next, double& sudakov0_n) {
CKKW_PrepareSudakov(Born, running);
for(DipoleChain const & chain : DSH()->eventRecord().chains()){
for(Dipole const & dip : chain.dipoles()){
sudakov0_n*=singlesudakov( dip, next,running,make_pair(true,false) );
sudakov0_n*=singlesudakov( dip, next,running,make_pair(false,true) );
if (sudakov0_n==0.0){
cleanup(Born);
return false;
}
}
}
cleanup(Born);
return true;
}
bool Merger::doUNLOPS(NPtr Born,Energy running, Energy next,Energy fixedScale, double& UNLOPS) {
CKKW_PrepareSudakov(Born, running);
for(DipoleChain const & chain : DSH()->eventRecord().chains()){
for(Dipole const & dip : chain.dipoles()){
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(true,false) );;
UNLOPS+=singleUNLOPS( dip, next,running,fixedScale,make_pair(false,true) );
}
}
cleanup(Born);
return true;
}
bool Merger::isProjectorStage(NPtr Born,int pjs){
return (pjs == int((Born->deepHead()->legsize() - Born->legsize())));
}
void Merger::cleanup(NPtr Born) {
DSH()->eventRecord().clear();
if(!Born->xcomb()->subProcess())return;
ParticleVector vecfirst = Born->xcomb()->subProcess()->incoming().first->children();
for(auto const & particle : vecfirst)
Born->xcomb()->subProcess()->incoming().first->abandonChild(particle);
ParticleVector vecsecond = Born->xcomb()->subProcess()->incoming().second->children();
for(auto const & particle : vecsecond)
Born->xcomb()->subProcess()->incoming().second->abandonChild(particle);
Born->xcomb()->subProcess(SubProPtr());
}
double Merger::singlesudakov(Dipole dip ,Energy next,Energy running,pair<bool,bool> conf ){
double res=1.;
tPPtr emitter = dip.emitter(conf);
tPPtr spectator = dip.spectator(conf);
DipoleSplittingInfo candidate(dip.index(conf),conf,
dip.emitterX(conf),
dip.spectatorX(conf),
emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
auto const & gens = DSH()->generators().equal_range(candidate.index());
for ( auto gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=gen->second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res*=gen->second->sudakov(candidate,next);
}
}
return res;
}
double Merger::singleUNLOPS(Dipole dip ,Energy next,Energy running,Energy fixedScale,pair<bool,bool> conf ){
double res=0.;
tPPtr emitter = dip.emitter(conf);
tPPtr spectator = dip.spectator(conf);
DipoleSplittingInfo candidate(dip.index(conf),conf,dip.emitterX(conf),dip.spectatorX(conf),emitter,spectator);
if ( DSH()->generators().find(candidate.index()) == DSH()->generators().end() ) DSH()->getGenerators(candidate.index());
auto const & gens = DSH()->generators().equal_range(candidate.index());
for ( auto gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
candidate.scale(dScale);
candidate.continuesEvolving();
Energy ptMax=gen->second->splittingKinematics()->ptMax(candidate.scale(),candidate.emitterX(), candidate.spectatorX(),
candidate.index(),*gen->second->splittingKernel());
candidate.hardPt(min(running,ptMax));
if (candidate.hardPt()>next){
res+=gen->second->unlops(candidate,next,fixedScale);
}
}
return res;
}
void Merger::firstNodeMap(Ptr<MatchboxMEBase>::ptr a,NPtr b){theFirstNodeMap.insert(make_pair(a,b));}
map<Ptr<MatchboxMEBase>::ptr,NPtr> Merger::firstNodeMap(){return theFirstNodeMap;}
void Merger::setXComb(Ptr<MatchboxMEBase>::ptr me,tStdXCombPtr xc){
theFirstNodeMap[me]->setXComb(xc);
}
void Merger::setKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->setKinematics();
}
void Merger::clearKinematics(Ptr<MatchboxMEBase>::ptr me){
theFirstNodeMap[me]->clearKinematics();
}
bool Merger::generateKinematics(Ptr<MatchboxMEBase>::ptr me,const double * r){
theFirstNodeMap[me]->firstgenerateKinematics(r, 0);
if (theFirstNodeMap[me]->cutStage()==0 ){
bool inAlphaPS=false;
NPtrVec children=theFirstNodeMap[me]->children();
for (Ptr<Node>::ptr const & child: children) {
treefactory()->setAlphaParameter(theGamma);
inAlphaPS|=theGamma!=1.?child->dipole()->aboveAlpha():false;
treefactory()->setAlphaParameter(1.);
}
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for (auto const & cs: temp) {
if (!cs.second.first&&!inAlphaPS)return false;
}
}
if (theFirstNodeMap[me]->cutStage()==1 ){
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for (auto const & sc: temp) {
if (!sc.second.first && !sc.second.second)return false;
}
}
return true;
}
void Merger::fillProjectors(Ptr<MatchboxMEBase>::ptr me){
for (auto const & propair: theFirstNodeMap[me]->Projector()) {
me->lastXCombPtr()->projectors().insert(propair.first,
propair.second->xcomb());
}
}
pair<bool,bool> Merger::clusterSafe(Ptr<MatchboxMEBase>::ptr me,int emit,int emis,int spec){
return theFirstNodeMap[me]->clusterSafe().find(make_pair(make_pair(emit,emis),spec))->second;
}
bool Merger::matrixElementRegion(PVector incoming,PVector outgoing,Energy winnerScale,Energy cutscale){
//cout<<"\nparticles s"<<particles.size()<<" "<<particles[0]<<" "<<particles[1]<<flush;
/*
if (defMERegionByJetAlg && !particles[0]->coloured()&& !particles[1]->coloured()) {
assert(false);
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
input_particles.push_back(fastjet::PseudoJet(em->momentum().x()/GeV,
em->momentum().y()/GeV,
em->momentum().z()/GeV,
em->momentum().e()/GeV));
}
fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets_ycut(ee_ycut);
return n==exclusive_jets.size();
}
if (defMERegionByJetAlg ) {
assert(false);
size_t noncol=0;
vector<fastjet::PseudoJet> input_particles;
for(size_t em=2; em < particles.size();em++){
if (em->coloured())
input_particles.push_back(fastjet::PseudoJet(em->momentum().x()/GeV,
em->momentum().y()/GeV,
em->momentum().z()/GeV,
em->momentum().e()/GeV));
else
noncol++;
}
double Rparam = 1.0;
fastjet::Strategy strategy = fastjet::Best;
fastjet::RecombinationScheme recomb_scheme = fastjet::E_scheme;
fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
size_t n = particles.size()-2-noncol;
vector<fastjet::PseudoJet> exclusive_jets = clust_seq.exclusive_jets(pp_dcut);
// cout<<"\nn= "<<n<<" jets= "<<exclusive_jets.size();
//for (size_t i=0; i<exclusive_jets.size(); i++) {
//cout<<"\nn= "<<n<<" jetspt= "<<exclusive_jets[i].perp();
//}
return n==exclusive_jets.size();
}
*/
Energy ptx=1000000.*GeV;
bool foundwinnerpt=false;
using namespace boost;
//FF
for(auto const & em : outgoing){ if (! em->coloured()) continue;
for(auto const & emm : outgoing){ if (!emm->coloured()) continue; if (em==emm) continue;
for(auto const & spe : outgoing){ if (!spe->coloured()) continue; if (em==spe||emm==spe) continue;
if (!(em->id()==-emm->id()||emm->id()>6))continue;
const Lorentz5Momentum emittermom = em->momentum();
const Lorentz5Momentum emissionmom = emm->momentum();
const Lorentz5Momentum spectatormom = spe->momentum();
Energy pt=0*GeV;
if (emittermom.m()<=0.001*GeV&&emissionmom.m()<=0.001*GeV&&spectatormom.m()<=0.001*GeV) {
pt=FFLTK->lastPt(emittermom,emissionmom,spectatormom);
}else{
pt=FFMTK->lastPt(emittermom,emissionmom,spectatormom);
}
if (abs(pt-winnerScale)<0.001*GeV) {
foundwinnerpt=true;
}
ptx =min(ptx,pt);
}
}
}
//FI
for(auto const & spe : incoming){ if (! spe->coloured()) continue;
for(auto const & emm : outgoing){ if (! emm->coloured()) continue;
for(auto const & em : outgoing){ if (! em->coloured()) continue; if (em==emm) continue;
if (!(em->id()==-emm->id()||emm->id()>6))continue;
const Lorentz5Momentum emittermom = em->momentum();
const Lorentz5Momentum emissionmom = emm->momentum();
const Lorentz5Momentum spectatormom = spe->momentum();
Energy pt=0*GeV;
if (emittermom.m()<=0.001*GeV&&emissionmom.m()<=0.001*GeV&&spectatormom.m()<=0.001*GeV) {
pt=FILTK->lastPt(emittermom,emissionmom,spectatormom);
}else{
pt=FIMTK->lastPt(emittermom,emissionmom,spectatormom);
}
if (abs(pt-winnerScale)<0.001*GeV) {
foundwinnerpt=true;
}
if(pt>0.*GeV)
ptx =min(ptx,pt);
}
}
}
//IF
for(auto const & em : incoming){ if (! em->coloured()) continue;
for(auto const & emm : outgoing){ if (! emm->coloured()) continue;
for(auto const & spe : outgoing){ if (! spe->coloured()) continue; if (emm==spe) continue;
if (!(em->id()>6|| em->id()==emm->id() ||emm->id()>6))continue;
const Lorentz5Momentum emittermom = em->momentum();
const Lorentz5Momentum emissionmom = emm->momentum();
const Lorentz5Momentum spectatormom = spe->momentum();
Energy pt=0*GeV;
if (emittermom.m()<=0.001*GeV&&emissionmom.m()<=0.001*GeV&&spectatormom.m()<=0.001*GeV) {
pt=IFLTK->lastPt(emittermom,emissionmom,spectatormom);
}else{
pt=IFMTK->lastPt(emittermom,emissionmom,spectatormom);
}
if (abs(pt-winnerScale)<0.01*GeV) {
foundwinnerpt=true;
}
ptx =min(ptx,pt);
}
}
}
//II
for(auto const & em : incoming){ if (! em->coloured()) continue;
for(auto const & spe : incoming){ if (! spe->coloured()) continue; if (em==spe) continue;
for(auto const & emm : outgoing){ if (! emm->coloured()) continue;
if (!(em->id()>6||em->id()==emm->id() ||emm->id()>6))continue;
const Lorentz5Momentum emittermom = em->momentum();
const Lorentz5Momentum emissionmom = emm->momentum();
const Lorentz5Momentum spectatormom = spe->momentum();
Energy pt=IILTK->lastPt(emittermom,emissionmom,spectatormom);
if (abs(pt-winnerScale)<0.01*GeV) {
foundwinnerpt=true;
}
ptx =min(ptx, pt);
}
}
}
if(!foundwinnerpt){
generator()->logWarning(Exception()
<< "Merger: Could not find winner with pt."
<< "Run with -d2 to get phase space points. "
<< Exception::warning);
if(Debug::level > 2) {
cout<<"\nWarning: could not find winner with pt: "<<winnerScale/GeV;
for(auto const & emm : incoming){
cout<<"\n"<<emm->id()<<" "<<emm->momentum()/GeV<<" "<<emm->momentum().m()/GeV<<flush;
}
for(auto const & emm : outgoing){
cout<<"\n"<<emm->id()<<" "<<emm->momentum()/GeV<<" "<<emm->momentum().m()/GeV<<flush;
}
}
}
return (ptx>cutscale) ;
}
int Merger::M()const{return theTreeFactory->M();}
int Merger::N()const{return theTreeFactory->N();}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Merger::persistentOutput(PersistentOStream & os) const {
os << Unlopsweights<<
theCMWScheme<< projected<<
isUnitarized<< isNLOUnitarized<<
defMERegionByJetAlg<<theOpenInitialStateZ<<
theChooseHistory<<
theN0<< theOnlyN<<
xiRenME<< xiFacME<<
xiRenSh<< xiFacSh<<
xiQSh<< weight<<weightCB<<theGamma<<ee_ycut<<pp_dcut<<theSmearing<<ounit(therenormscale, GeV)<<ounit(theIRSafePT, GeV)<<ounit(theMergePt, GeV)<<ounit(theCentralMergePt, GeV)<<theMergingJetFinder
<<theLargeNBasis
<< FFLTK
<< FILTK
<< IFLTK
<< IILTK
<< FFMTK
<< FIMTK
<< IFMTK
<<theDipoleShowerHandler<< theTreeFactory<<theFirstNodeMap ;
}
void Merger::persistentInput(PersistentIStream & is, int) {
is >> Unlopsweights>>
theCMWScheme>> projected>>
isUnitarized>> isNLOUnitarized>>
defMERegionByJetAlg>>theOpenInitialStateZ>>
theChooseHistory>>
theN0>> theOnlyN>>
xiRenME>> xiFacME>>
xiRenSh>> xiFacSh>>
xiQSh>>
weight>>weightCB>>theGamma>>ee_ycut>>pp_dcut>>theSmearing>>iunit(therenormscale, GeV)>>iunit(theIRSafePT, GeV)>>iunit(theMergePt, GeV)>>iunit(theCentralMergePt, GeV)>>theMergingJetFinder>>theLargeNBasis>>
FFLTK
>> FILTK
>> IFLTK
>> IILTK
>> FFMTK
>> FIMTK
>> IFMTK>>
theDipoleShowerHandler>> theTreeFactory>>theFirstNodeMap ;
}
// *** 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<Merger, Herwig::MergerBase>
describeHerwigMerger("Herwig::Merger", "HwDipoleShower.so");
//ClassDescription<Merger> Merger::initMerger;
// Definition of the static class description member.
void Merger::Init() {
static ClassDocumentation<Merger> documentation
("Merger.");
static Reference<Merger,DipoleShowerHandler> interfaceShowerHandler
("DipoleShowerHandler",
"",
&Merger::theDipoleShowerHandler, false, false, true, true, false);
static Switch<Merger,bool>
interfaceUnlopsweights ("Unlopsweights","",&Merger::Unlopsweights, false, false, false);
static SwitchOption interfaceUnlopsweightsYes
(interfaceUnlopsweights,"Yes","",true);
static SwitchOption interfaceUnlopsweightsNo
(interfaceUnlopsweights,"No","",false);
static Switch<Merger,bool>
interfacetheCMWScheme ("CMWScheme","",&Merger::theCMWScheme, false, false, false);
static SwitchOption interfacetheCMWSchemeYes
(interfacetheCMWScheme,"Yes","",true);
static SwitchOption interfacetheCMWSchemeNo
(interfacetheCMWScheme,"No","",false);
static Parameter<Merger,Energy> interfaceMergerScale
("MergingScale",
"The MergingScale.",
&Merger::theCentralMergePt, GeV, 20.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Reference<Merger,MergingFactory> interfaceMergerHelper
("MergingFactory",
"",
&Merger::theTreeFactory, false, false, true, true, false);
static Parameter<Merger,double> interfacedcut
("pp_dcut",
"The MergingScale.",
&Merger::pp_dcut, 5.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,double> interfaceycut
("ee_ycut",
"The MergingScale.",
&Merger::ee_ycut, 0.2, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<Merger,double> interfacegamma
("gamma",
"gamma parameter.",
&Merger::theGamma, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Reference<Merger,JetFinder> interfaceMergingJetFinder
("MergingJetFinder",
"A reference to the jet finder used in Merging.",
&Merger::theMergingJetFinder, false, false, true, false, false);
static Reference<Merger,ColourBasis> interfaceLargeNBasis
("LargeNBasis",
"Set the large-N colour basis implementation.",
&Merger::theLargeNBasis, false, false, true, true, false);
static Switch<Merger,bool>
interfacedefMERegionByJetAlg
("MERegionByJetAlg","",&Merger::defMERegionByJetAlg, false, false, false);
static SwitchOption interfacedefMERegionByJetAlgYes
(interfacedefMERegionByJetAlg,"Yes","",true);
static SwitchOption interfacedefMERegionByJetAlgNo
(interfacedefMERegionByJetAlg,"No","",false);
static Switch<Merger,bool>
interfaceOpenInitialSateZ
("OpenInitialStateZ","",&Merger::theOpenInitialStateZ, false, false, false);
static SwitchOption interfaceOpenInitialSateZYes
(interfaceOpenInitialSateZ,"Yes","",true);
static SwitchOption interfaceOpenInitialSateZNo
(interfaceOpenInitialSateZ,"No","",false);
static Parameter<Merger, Energy>
interfaceIRSafePT
("IRSafePT", "Set the pt for which a matrixelement should be treated as IR-safe.",
&Merger::theIRSafePT,
GeV, 0.0 * GeV, ZERO, Constants::MaxEnergy, true, false, Interface::limited);
interfaceIRSafePT.setHasDefault(false);
static Parameter<Merger, double> interfacemergePtsmearing("MergingScaleSmearing", "Set the percentage the merging pt should be smeared.",
&Merger::theSmearing, 0., 0.,
0.0, 0.5, true, false, Interface::limited);
static Parameter<Merger, int> interfacechooseHistory("chooseHistory", "different ways of choosing the history", &Merger::theChooseHistory, 3, -1, 0,
false, false, Interface::lowerlim);
}
diff --git a/Shower/Dipole/Merging/MergingFactory.cc b/Shower/Dipole/Merging/MergingFactory.cc
--- a/Shower/Dipole/Merging/MergingFactory.cc
+++ b/Shower/Dipole/Merging/MergingFactory.cc
@@ -1,671 +1,678 @@
// -*- C++ -*-
//
// MergeboxFactory.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 MergeboxFactory class.
//
#include "MergingFactory.h"
#include "Node.h"
#include "Herwig/MatrixElement/Matchbox/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 "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "ThePEG/MatrixElement/MEBase.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 "ThePEG/Cuts/JetFinder.h"
#include "fastjet/ClusterSequence.hh"
#include <boost/progress.hpp>
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
using std::ostream_iterator;
MergingFactory::MergingFactory() :MatchboxFactory(),
theonlyNLOParts(false),
theonlyvirtualNLOParts(false),
theonlyrealNLOParts(false),
theonlyUnlopsweights(false),
theunitarizeNLOParts(true),
calc_born(true),
calc_virtual(true),
calc_real(true),
unitarized(true),
NLOunitarized(true),
ransetup(false),
theonlymulti(-1),
theonlysub(-1),
divideSub(-1),
divideSubNumber(-1)
{}
MergingFactory::~MergingFactory(){}
IBPtr MergingFactory::clone() const {
return new_ptr(*this);
}
IBPtr MergingFactory::fullclone() const {
return new_ptr(*this);
}
void MergingFactory::fill_amplitudes() {
olpProcesses().clear();
processMap[0] = getProcesses()[0];
if(MH()->M()>=0) setHighestVirt(processMap[0].size()+MH()->M());
MH()->N0(processMap[0].size());
for ( int i = 1 ; i <= MH()->N() ; ++i ) {
processMap[i] = processMap[i - 1];
processMap[i].push_back("j");
}
for ( int i = 0 ; i <= MH()->N() ; ++i ) {
vector<Ptr<MatchboxMEBase>::ptr> ames = makeMEs(processMap[i], orderInAlphaS() + i,i<MH()->M()+1);
copy(ames.begin(), ames.end(), back_inserter(pureMEsMap()[i]));
cout<<"\n"<<processMap[i].size()<<" "<<pureMEsMap()[i].size()<<" "<<orderInAlphaS()<<" "<<orderInAlphaEW()<<" "<<MH()->M()<<" "<< MH()->N0()<<" "<< MH()->N();
}
}
void MergingFactory::prepare_BV(int i) {
// check if we have virtual contributions
bool haveVirtuals = true;
for ( auto born : pureMEsMap()[i]) {
prepareME(born);
if ( born->isOLPTree() ) {
int id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::treeME2);
born->olpProcess(ProcessType::treeME2,id);
id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
born->olpProcess(ProcessType::colourCorrelatedME2,id);
bool haveGluon = false;
for ( const auto p : born->subProcess().legs )
if ( p->id() == 21 ) {
haveGluon = true;
break;
}
if ( haveGluon ) {
id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::spinColourCorrelatedME2);
born->olpProcess(ProcessType::spinColourCorrelatedME2,id);
}
}
if ( born->isOLPLoop() && i <= MH()->M()) {
int id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::oneLoopInterference);
born->olpProcess(ProcessType::oneLoopInterference,id);
if ( !born->onlyOneLoop() && born->needsOLPCorrelators() ) {
id = orderOLPProcess(born->subProcess(),
born->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
born->olpProcess(ProcessType::colourCorrelatedME2,id);
}
}
haveVirtuals &= born->haveOneLoop();
}
// check the additional insertion operators
if ( !theVirtualsMap[i].empty() ) haveVirtuals = true;
for ( auto & virt : theVirtualsMap[i] ) virt->factory(this);
// check for consistent conventions on virtuals, if we are to include MH()->M()
assert(i > MH()->M()||haveVirtuals);
for ( auto & virt : DipoleRepository::insertionIOperators(dipoleSet()) )
virt->factory(this);
for ( auto & virt : DipoleRepository::insertionPKOperators(dipoleSet()) )
virt->factory(this);
}
void MergingFactory::prepare_R(int i) {
for ( auto real : pureMEsMap()[i])
prepareME(real);
}
void MergingFactory::pushB(Ptr<MatchboxMEBase>::ptr born, int i) {
Ptr<MatchboxMEBase>::ptr bornme = born->cloneMe();
bornme->maxMultCKKW(1);
bornme->minMultCKKW(0);
string pname = fullName() + "/" + bornme->name() + ".Born";
if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Born ME " << pname << " already existing.";
bornme->virtuals().clear();
for ( auto virt : theVirtualsMap[i] )
if ( virt->apply(born->diagrams().front()->partons()) )
bornme->virtuals().push_back(virt);
for ( auto I : DipoleRepository::insertionIOperators(dipoleSet()) )
if ( I->apply(born->diagrams().front()->partons()) )
bornme->virtuals().push_back(I);
for ( auto PK : DipoleRepository::insertionPKOperators(dipoleSet()) )
if ( PK->apply(born->diagrams().front()->partons()) )
bornme->virtuals().push_back(PK);
Ptr<Node>::ptr clusternode = new_ptr(Node(bornme, 0,MH()));
+
+ clusternode->deepHead(clusternode);
MH()->firstNodeMap(bornme,clusternode);
bornme->factory(this);
bornme->merger(MH());
vector<Ptr<Node>::ptr> temp;
vector<Ptr<Node>::ptr> temp1;
temp.push_back(clusternode);
unsigned int k = 1;
while (thePureMEsMap[i - k].size() != 0) {
for ( auto tmp : temp){//j
tmp->birth(thePureMEsMap[i - k]);
for ( auto tmpchild : tmp->children()){//m
if ( i <= MH()->M() ) {
for ( auto I : DipoleRepository::insertionIOperators(dipoleSet()) )
if ( I->apply(tmpchild->nodeME()->diagrams().front()->partons()) ){
Ptr<MatchboxInsertionOperator>::ptr myI = I->cloneMe();
tmpchild->nodeME()->virtuals().push_back(myI);
}
for ( auto I : DipoleRepository::insertionPKOperators(dipoleSet()) )
if ( I->apply(tmpchild->nodeME()->diagrams().front()->partons()) ){
Ptr<MatchboxInsertionOperator>::ptr myI = I->cloneMe();
tmpchild->nodeME()->virtuals().push_back(myI);
}
tmpchild->nodeME()->noOneLoop();
}
temp1.push_back(tmpchild);
}
}
temp = temp1;
temp1.clear();
k++;
}
if(MH()->N()>i)
bornme->needsCorrelations();
else bornme->needsNoCorrelations();
bornme->cloneDependencies();
MEs().push_back(bornme);
}
void MergingFactory::pushV(Ptr<MatchboxMEBase>::ptr born, int i) {
Ptr<MatchboxMEBase>::ptr nlo = born->cloneMe();
string pname = fullName() + "/" + nlo->name() + ".Virtual";
if ( !(generator()->preinitRegister(nlo, pname)) ) throw InitException() << "Virtual ME " << pname << " already existing.";
////////////////////////////////////NLO///////////////////////////
nlo->virtuals().clear();
if ( !nlo->onlyOneLoop() ) {
for ( auto OV : theVirtualsMap[i] )
if ( OV->apply(nlo->diagrams().front()->partons()) ){
nlo->virtuals().push_back(OV);
}
for ( auto I : DipoleRepository::insertionIOperators(dipoleSet()) )
if ( I->apply(nlo->diagrams().front()->partons()) ){
nlo->virtuals().push_back(I);
}
for ( auto PK : DipoleRepository::insertionPKOperators(dipoleSet()) )
if ( PK->apply(nlo->diagrams().front()->partons()) ){
nlo->virtuals().push_back(PK);
}
if ( nlo->virtuals().empty() ) throw InitException() << "No insertion operators have been found for " << born->name() << "\n";
}
nlo->doOneLoopNoBorn();
////////////////////////////////////NLO///////////////////////////
Ptr<Node>::ptr clusternode = new_ptr(Node(nlo, 0,MH()));
+
+ clusternode->deepHead(clusternode);
clusternode->virtualContribution(true);
MH()->firstNodeMap(nlo,clusternode);
nlo->merger(MH());
vector<Ptr<Node>::ptr> temp;
vector<Ptr<Node>::ptr> temp1;
temp.push_back(clusternode);
unsigned int k = 1;
while (thePureMEsMap[i - k].size() != 0) {
for ( auto tmp : temp ){
tmp->birth(thePureMEsMap[i - k]);
for ( auto tmpchild : tmp->children())
temp1.push_back(tmpchild);
}
temp = temp1;
k++;
}
if ( nlo->isOLPLoop() ) {
int id = orderOLPProcess(nlo->subProcess(),
born->matchboxAmplitude(),
ProcessType::oneLoopInterference);
nlo->olpProcess(ProcessType::oneLoopInterference,id);
if ( !nlo->onlyOneLoop() && nlo->needsOLPCorrelators() ) {
id = orderOLPProcess(nlo->subProcess(),
born->matchboxAmplitude(),
ProcessType::colourCorrelatedME2);
nlo->olpProcess(ProcessType::colourCorrelatedME2,id);
}
}
nlo->needsCorrelations();
nlo->cloneDependencies();
MEs().push_back(nlo);
}
void MergingFactory::pushProR(Ptr<MatchboxMEBase>::ptr born, int i) {
Ptr<MatchboxMEBase>::ptr bornme = born->cloneMe();
string pname = fullName() + "/" + bornme->name() + ".Real";
if ( !(generator()->preinitRegister(bornme, pname)) ) throw InitException() << "Subtracted ME " << pname << " already existing.";
Ptr<Node>::ptr clusternode = new_ptr(Node(bornme, 1,MH()));
+ clusternode->deepHead(clusternode);
clusternode->subtractedReal(true);
MH()->firstNodeMap(bornme,clusternode);
bornme->merger(MH());
+
+
vector<Ptr<Node>::ptr> temp;
vector<Ptr<Node>::ptr> temp1;
temp.push_back(clusternode);
unsigned int k = 1;
while (thePureMEsMap[i - k].size() != 0) {
for ( auto tmp : temp ){
tmp->birth(thePureMEsMap[i - k]);
for ( auto tmpchild : tmp->children()) temp1.push_back(tmpchild);
}
temp = temp1;
k++;
}
if(MH()->N()>i)
bornme->needsCorrelations();
else bornme->needsNoCorrelations();
bornme->cloneDependencies(pname);
MEs().push_back(bornme);
}
void MergingFactory::orderOLPs() {
}
vector<string> MergingFactory::parseProcess(string in) {
vector<string> process = StringUtils::split(in);
if ( process.size() < 3 )
throw Exception() << "MatchboxFactory: Invalid process."
<< Exception::runerror;
for ( string & p : process) {
p = StringUtils::stripws(p);
}
theN=0;
bool prodprocess=true;
vector<string> pprocess;
for (string const p : process) {
if ( p == "->" )
continue;
if (p=="[") {
prodprocess=false;
}else if (p=="]") {
prodprocess=false;
break;
}else if (p=="[j") {
prodprocess=false;
theN++;
}else if (p=="j"&&!prodprocess) {
theN++;
prodprocess=false;
}else if (p=="j]") {
theN++;
prodprocess=false;
break;
}else if(prodprocess){
pprocess.push_back(p);
}else{
cout<<"\nWarning: "<<p<<" in the brackets of the process definition is not recognized.\n Only j as jets are recognized.";
}
}
return pprocess;
}
void MergingFactory::setup() {
useMe();
if(!ransetup){
olpProcesses().clear();
externalAmplitudes().clear();
setFixedCouplings(true);
setFixedQEDCouplings(true);
const PDVector& partons = particleGroups()["j"];
unsigned int nl = 0;
// rebind the particle data objects
for ( auto & g : particleGroups()) {
for ( auto & p : g.second) {
p = getParticleData(p->id());
}
}
for ( const auto p : partons ) {
if ( abs(p->id()) < 7 && p->hardProcessMass() == ZERO )
++nl;
if ( p->id() > 0 && p->id() < 7 && p->hardProcessMass() == ZERO )
nLightJetVec( p->id() );
if ( p->id() > 0 && p->id() < 7 && p->hardProcessMass() != ZERO )
nHeavyJetVec( p->id() );
}
nLight(nl/2);
const PDVector& partonsInP = particleGroups()["p"];
for ( const auto pip : partonsInP )
if ( pip->id() > 0 && pip->id() < 7 && pip->hardProcessMass() == ZERO )
nLightProtonVec( pip->id() );
for ( auto & amp: amplitudes() ) amp->factory(this);
MH()->largeNBasis()->factory(this);
assert(!(divideSub!=-1&÷SubNumber==-1)||!(divideSub==-1&÷SubNumber!=-1));
assert(!subProcessGroups());
//fill the amplitudes
if ( !amplitudes().empty() ) fill_amplitudes();
// prepare the Born and virtual matrix elements
for ( int i = 0 ; i <= max(0, MH()->N()) ; ++i ) prepare_BV(i);
// prepare the real emission matrix elements
for ( int i = 0 ; i <= MH()->N() ; ++i ) prepare_R(i);
orderOLPs();
// start creating matrix elements
MEs().clear();
int onlysubcounter=0;
int i = 0 ;
for (; i <= max(0, MH()->N()) ; ++i )
for ( auto born : thePureMEsMap[i])
if (calc_born && !theonlyUnlopsweights){
if(((theonlysub==-1||theonlysub==onlysubcounter)&÷Sub==-1)
||(divideSub!=-1&&onlysubcounter%divideSub==divideSubNumber))
pushB(born, i);
onlysubcounter++;
}
i = 0 ;
for (; i <=max(0, MH()->N()); ++i )
for ( auto virt : thePureMEsMap[i])
if ( calc_virtual && i <= MH()->M()){
if(((theonlysub==-1||theonlysub==onlysubcounter)&÷Sub==-1)
||(divideSub!=-1&&onlysubcounter%divideSub==divideSubNumber))
pushV(virt, i);
onlysubcounter++;
}
i = 0 ;
for (; i <= max(0, MH()->N()) ; ++i )
for ( auto real : thePureMEsMap[i] )
if (calc_real&& i <= MH()->M() + 1 && i > 0 && !theonlyvirtualNLOParts&&!theonlyUnlopsweights){
if(((theonlysub==-1||theonlysub==onlysubcounter)&÷Sub==-1)
||(divideSub!=-1&&onlysubcounter%divideSub==divideSubNumber))
pushProR(real, i);
onlysubcounter++;
}
if ( !externalAmplitudes().empty() ) {
generator()->log() << "Initializing external amplitudes.\n" << flush;
boost::progress_display * progressBar =
new boost::progress_display(externalAmplitudes().size(),generator()->log());
for ( const auto ext : externalAmplitudes()) {
if ( ! ext->initializeExternal() ) {
throw InitException()
<< "error: failed to initialize amplitude '" << ext->name() << "'\n";
}
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n"
<< flush;
}
if ( !olpProcesses().empty() ) {
generator()->log() << "Initializing one-loop provider(s).\n" << flush;
map<Ptr<MatchboxAmplitude>::tptr, map<pair<Process, int>, int> > olps;
for (const auto oit : olpProcesses()) {
olps[oit.first] = oit.second;
}
boost::progress_display * progressBar = new boost::progress_display(olps.size(), generator()->log());
for ( const auto olpit : olps ) {
if ( !olpit.first->startOLP(olpit.second) ) {
throw InitException() << "error: failed to start OLP for amplitude '" << olpit.first->name() << "'\n";
}
++(*progressBar);
}
delete progressBar;
generator()->log() << "--------------------------------------------------------------------------------\n" << flush;
}
cout<<"\n Generated "<<MEs().size()<<" Subprocesses."<<flush;
if(theonlysub!=-1)cout<<" ( "<<theonlysub<<"/"<<onlysubcounter<<" )"<<flush;
cout<<"\n"<<flush;
generator()->log() << "process setup finished.\n" << flush;
ransetup=true;
}
}
void MergingFactory::persistentOutput(PersistentOStream & os) const {
os
<< calc_born << calc_virtual << calc_real
<< theonlyUnlopsweights << theonlymulti
<< divideSub << divideSubNumber
<< theonlysub << ransetup
<< processMap << theMergingHelper <<theM<<theN;
}
void MergingFactory::persistentInput(PersistentIStream & is, int) {
is
>> calc_born >> calc_virtual >> calc_real
>> theonlyUnlopsweights >> theonlymulti
>> divideSub >> divideSubNumber
>> theonlysub >> ransetup
>> processMap >> theMergingHelper >>theM>>theN;
}
void MergingFactory::Init() {
static Parameter<MergingFactory, int> interfaceonlymulti("onlymulti", "calculate only the ME with k additional partons.", &MergingFactory::theonlymulti, -1, -1, 0,
false, false, Interface::lowerlim);
static Parameter<MergingFactory, int> interfaceonlysub("onlysub", "calculate only one subProcess. this is for building grids.", &MergingFactory::theonlysub, -1, -1, 0,
false, false, Interface::lowerlim);
static Parameter<MergingFactory, int> interfacedivideSub("divideSub", "calculate only one subProcess. this is for building grids.", &MergingFactory::divideSub, -1, -1, 0,
false, false, Interface::lowerlim);
static Parameter<MergingFactory, int> interfacedivideSubNumber("divideSubNumber", "calculate only one subProcess. this is for building grids.", &MergingFactory::divideSubNumber, -1, -1, 0,
false, false, Interface::lowerlim);
static Switch<MergingFactory, bool> interface_calc_born("calc_born", "[debug] Switch on or off the born contribution.", &MergingFactory::calc_born, true,
false, false);
static SwitchOption interface_calc_bornOn(interface_calc_born, "On", "Switch on calculation of born.", true);
static SwitchOption interface_calc_bornOff(interface_calc_born, "Off", "Switch off calculation of born.", false);
static Switch<MergingFactory, bool> interface_calc_virtual("calc_virtual", "[debug] Switch on or off the virtual contribution.",
&MergingFactory::calc_virtual, true, false, false);
static SwitchOption interface_calc_virtualOn(interface_calc_virtual, "On", "Switch on calculation of virtual.", true);
static SwitchOption interface_calc_virtualOff(interface_calc_virtual, "Off", "Switch off calculation of virtual.", false);
static Switch<MergingFactory, bool> interface_calc_real("calc_real", "[debug] Switch on or off the real contribution.", &MergingFactory::calc_real, true,
false, false);
static SwitchOption interface_calc_realOn(interface_calc_real, "On", "Switch on calculation of real.", true);
static SwitchOption interface_calc_realOff(interface_calc_real, "Off", "Switch off calculation of real.", false);
static Switch<MergingFactory, bool> interface_theonlyNLOParts("onlyNLOParts", "Switch on or off the onlyNLOParts.", &MergingFactory::theonlyNLOParts, true, false,
false);
static SwitchOption interface_theonlyNLOPartsOn(interface_theonlyNLOParts, "On", "Switch on the theonlyNLOParts.", true);
static SwitchOption interface_theonlyNLOPartsOff(interface_theonlyNLOParts, "Off", "Switch off the theonlyNLOParts.", false);
// theonlyvirtualNLOParts;
// theonlyrealNLOParts;
static Switch<MergingFactory, bool> interface_theonlyvirtualNLOParts("onlyvirtualNLOParts", "Switch on or off the onlyvirtualNLOParts.", &MergingFactory::theonlyvirtualNLOParts, true, false,
false);
static SwitchOption interface_theonlyvirtualNLOPartsOn(interface_theonlyvirtualNLOParts, "On", "Switch on the theonlyvirtualNLOParts.", true);
static SwitchOption interface_theonlyvirtualNLOPartsOff(interface_theonlyvirtualNLOParts, "Off", "Switch off the theonlyvirtualNLOParts.", false);
static Switch<MergingFactory, bool> interface_theonlyrealNLOParts("onlyrealNLOParts", "Switch on or off the onlyrealNLOParts.", &MergingFactory::theonlyrealNLOParts, true, false,
false);
static SwitchOption interface_theonlyrealNLOPartsOn(interface_theonlyrealNLOParts, "On", "Switch on the theonlyrealNLOParts.", true);
static SwitchOption interface_theonlyrealNLOPartsOff(interface_theonlyrealNLOParts, "Off", "Switch off the theonlyrealNLOParts.", false);
static Switch<MergingFactory, bool> interface_theunitarizeNLOParts("unitarizeNLOParts", "Switch on or off the unitarizeNLOParts.", &MergingFactory::theunitarizeNLOParts, true, false,
false);
static SwitchOption interface_theunitarizeNLOPartsOn(interface_theunitarizeNLOParts, "On", "Switch on the unitarizeNLOParts.", true);
static SwitchOption interface_theunitarizeNLOPartsOff(interface_theunitarizeNLOParts, "Off", "Switch off the unitarizeNLOParts.", false);
static Switch<MergingFactory, bool> interface_theonlyUnlopsweights("onlyUnlopsweights", "Switch on or off the onlyUnlopsweights.", &MergingFactory::theonlyUnlopsweights, true, false,
false);
static SwitchOption interface_theonlyUnlopsweightsOn(interface_theonlyUnlopsweights, "On", "Switch on the onlyUnlopsweights.", true);
static SwitchOption interface_theonlyUnlopsweightsOff(interface_theonlyUnlopsweights, "Off", "Switch off the onlyUnlopsweights.", false);
static Switch<MergingFactory, bool> interface_Unitarized("Unitarized", "Unitarize the cross section.", &MergingFactory::unitarized, true, false, false);
static SwitchOption interface_UnitarizedOn(interface_Unitarized, "On", "Switch on the unitarized cross section.", true);
static SwitchOption interface_UnitarizedOff(interface_Unitarized, "Off", "Switch off the unitarized cross section.", false);
static Switch<MergingFactory, bool> interface_NLOUnitarized("NLOUnitarized", "Unitarize the cross section.", &MergingFactory::NLOunitarized, true, false, false);
static SwitchOption interface_NLOUnitarizedOn(interface_NLOUnitarized, "On", "Switch on the unitarized NLO cross section.", true);
static SwitchOption interface_NLOUnitarizedOff(interface_NLOUnitarized, "Off", "Switch off the unitarized NLO cross section.", false);
static Reference<MergingFactory,Merger> interfaceMergingHelper
("MergingHelper",
"",
&MergingFactory::theMergingHelper, false, false, true, true, false);
static Parameter<MergingFactory, int> interfaceaddNLOLegs("NLOProcesses",
"Set the number of virtual corrections to consider. 0 is default for no virtual correction.", &MergingFactory::theM, 0, 0, 0, false, false,
Interface::lowerlim);
}
// *** 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<MergingFactory, Herwig::MatchboxFactory> describeHerwigMergingFactory("Herwig::MergingFactory", "HwDipoleShower.so");
diff --git a/Shower/Dipole/Merging/Node.cc b/Shower/Dipole/Merging/Node.cc
--- a/Shower/Dipole/Merging/Node.cc
+++ b/Shower/Dipole/Merging/Node.cc
@@ -1,405 +1,403 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Node class.
//
#include "Node.h"
#include "MergingFactory.h"
#include "Merger.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "Herwig/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
Node::Node() {
}
-bool NodeDebug=false;
Node::Node(Ptr<MatchboxMEBase>::ptr nodeME, int cutstage,Ptr<Merger>::ptr mh)
:Interfaced(),
thenodeMEPtr(nodeME),
thedipol(),
- thechildren(),
theparent(),
- theDeepHead(this),
+
theCutStage(cutstage),
isOrdered(true),
theSubtractedReal(false),
theVirtualContribution(false),
theMergingHelper(mh)
{
nodeME->maxMultCKKW(1);
nodeME->minMultCKKW(0);
}
Node::Node(Ptr<Node>::ptr deephead,
Ptr<Node>::ptr head,
Ptr<SubtractionDipole>::ptr dipol,
Ptr<MatchboxMEBase>::ptr nodeME,
int cutstage)
:Interfaced(),thenodeMEPtr(nodeME),
thedipol(dipol),
theparent(head),
theDeepHead(deephead),
theCutStage(cutstage),
isOrdered(true),
theSubtractedReal(false),
-theVirtualContribution(false),
-theMergingHelper() //The subnodes have no merging helper
+theVirtualContribution(false)
+ //The subnodes have no merging helper
{
}
Node::~Node() { }
Ptr<SubtractionDipole>::ptr Node::dipole() const {
return thedipol;
}
/** returns the matrix element pointer */
const Ptr<MatchboxMEBase>::ptr Node::nodeME() const {
return thenodeMEPtr;
}
/** access the matrix element pointer */
Ptr<MatchboxMEBase>::ptr Node::nodeME() {
return thenodeMEPtr;
}
int Node::legsize() const {return nodeME()->legsize();}
Ptr<Node>::ptr Node::randomChild() {
return thechildren[(int)(UseRandom::rnd() * thechildren.size())];
}
bool Node::allAbove(Energy pt){
for (Ptr<Node>::ptr child : thechildren)
if(child->pT()<pt)return false;
return true;
}
bool Node::isInHistoryOf(Ptr<Node>::ptr other){
while (other->parent()) {
if(other==this)return true;
other=other->parent();
}
return false;
}
void Node::flushCaches() {
this->theProjectors.clear();
for ( auto const & ch: thechildren) {
ch->dipole()->underlyingBornME()->flushCaches();
for (Ptr<MatchboxReweightBase>::ptr r : ch->dipole()->reweights())
r->flushCaches();
if ( ch->xcomb() ) ch->xcomb()->clean();
ch->nodeME()->flushCaches();
ch->flushCaches();
}
}
void Node::setKinematics() {
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
ch->dipole()->setKinematics();
ch->nodeME()->setKinematics();
ch->setKinematics();
}
}
void Node::clearKinematics() {
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
ch->nodeME()->clearKinematics();
ch->dipole()->clearKinematics();
ch->clearKinematics();
}
}
bool Node::generateKinematics(const double *r, int stage, Energy2 ) {
bool isthissafe=true;
for (auto const & ch: thechildren) {
ch->dipole()->setXComb(ch->xcomb());
if ( !ch->dipole()->generateKinematics(r) ) cout << "stop";
ch->generateKinematics(r, stage + 1, ch->xcomb()->lastSHat());
isthissafe = (isthissafe && ch->pT() >=deepHead()->MH()->mergePt());
}
return isthissafe;
}
void Node::firstgenerateKinematics(const double *r, int stage) {
flushCaches();
MH()->smeareMergePt();
//Set here the new merge Pt for the next phase space point.( Smearing!!!)
clustersafer.clear();
for (auto const & ch: thechildren) {
bool ifirst = true;
bool isecond = true;
ch->dipole()->setXComb(ch->xcomb());
if ( !ch->dipole()->generateKinematics(r) ) cout << "stop";
isecond = ch->generateKinematics(r, stage + 1, ch->xcomb()->lastSHat());
ifirst = (ch->pT() >= deepHead()->MH()->mergePt());
pair<pair<int, int>, int> EmitEmisSpec =
make_pair(make_pair(ch->dipole()->realEmitter(),
ch->dipole()->realEmission()),
ch->dipole()->realSpectator());
clustersafer.insert(make_pair(EmitEmisSpec, make_pair(ifirst, isecond)));
}
}
void Node::setXComb(tStdXCombPtr xc) {
if ( !parent() ) this->xcomb(xc);
for (auto const & ch: thechildren) {
if ( !ch->xcomb() ) {
ch->xcomb(ch->dipole()->makeBornXComb(xc));
ch->xcomb()->head(xc);
if ( !ch->dipole()->lastXCombPtr() ) {
ch->dipole()->setXComb(ch->xcomb());
}
ch->setXComb(ch->xcomb());
} else {
if ( !(ch->dipole()->lastXCombPtr()->lastScale() == ch->xcomb()->lastScale()) ) {
ch->dipole()->setXComb(ch->xcomb());
}
if ( ch->xcomb()->head() != xc ) ch->xcomb()->head(xc);
ch->setXComb(ch->xcomb());
}
}
}
void Node::birth(vector<Ptr<MatchboxMEBase>::ptr> vec) {
vector<Ptr<SubtractionDipole>::ptr> dipoles =
nodeME()->getDipoles(DipoleRepository::dipoles(
nodeME()->factory()->dipoleSet()), vec,true);
for ( auto const & dip : dipoles ) {
dip->doSubtraction();
Ptr<Node>::ptr node = new_ptr(Node(theDeepHead,
this,
dip,
dip->underlyingBornME(),
theDeepHead->cutStage()));
thechildren.push_back(node);
}
}
vector<Ptr<Node>::ptr> Node::getNextOrderedNodes(bool normal,double hardScaleFactor) {
vector<Ptr<Node>::ptr> temp = children();
vector<Ptr<Node>::ptr> res;
for (Ptr<Node>::ptr const & child : children()) {
if(deepHead()->MH()->mergePt()>child->pT()){
res.clear();
return res;
}
}
for (Ptr<Node>::ptr const & child: children()) {
if (parent()&& normal){
if ( child->pT() < pT() ){
continue;
}
}
if ( child->children().size() != 0 ){
for (Ptr<Node>::ptr itChild: child->children()) {
if( itChild->pT() > child->pT()&&child->inShowerPS(itChild->pT()) ){
res.push_back(child);
break;
}
}
}
else {
child->nodeME()->factory()->scaleChoice()->setXComb(child->xcomb());
if ( sqr(hardScaleFactor)*child->nodeME()->factory()->scaleChoice()->renormalizationScale()
>= sqr(child->pT()) &&
child->inShowerPS(hardScaleFactor*sqrt(child->nodeME()->factory()->scaleChoice()->renormalizationScale()))) {
res.push_back(child);
}
}
}
return res;
}
bool Node::inShowerPS(Energy hardpT){
//Here we decide if the current phase space point can be reached from the underlying Node.
double z_=dipole()->lastZ();
// II
if( dipole()->bornEmitter()<2&&dipole()->bornSpectator()<2&&deepHead()->MH()->openInitialStateZ()) return true;
// IF
if( dipole()->bornEmitter()<2&&dipole()->bornSpectator()>=2&&deepHead()->MH()->openInitialStateZ())
return true;
pair<double,double> zbounds=
dipole()->tildeKinematics()->zBounds(pT(),hardpT);
return (zbounds.first<z_&&z_<zbounds.second);
}
Ptr<Node>::ptr Node::getHistory(bool normal,double hardScaleFactor) {
Ptr<Node>::ptr res=this;
//cout<<"\nstart get next"<<flush;
vector<Ptr<Node>::ptr> temp = getNextOrderedNodes(normal,hardScaleFactor);
Energy minpt=100000.*GeV;
Selector<Ptr<Node>::ptr> subprosel;
while (temp.size()!=0){
minpt=100000.*GeV;
subprosel.clear();
for (Ptr<Node>::ptr const & child : temp) {
if( child->dipole()->underlyingBornME()->largeNColourCorrelatedME2(
make_pair(child->dipole()->bornEmitter(),child->dipole()->bornSpectator()),
deepHead()->MH()->largeNBasis())!=0.
){
double weight=abs(child->dipole()->dSigHatDR()/nanobarn);
if(weight!=0.){
subprosel.insert(weight , child);
minpt=min(minpt,child->pT());
}
/*
if((*it)->nodeME()->dSigHatDR()/nanobarn!=0.){
subprosel.insert((abs((*it)->dipole()->dSigHatDR() /
(*it)->nodeME()->dSigHatDR()*deepHead()->MH()->as((*it)->pT()))), (*it));
minpt=min(minpt,(*it)->pT());
}
*/
//TODO choosehistories
}
}
if (subprosel.empty())
return res;
res = subprosel.select(UseRandom::rnd());
temp = res->getNextOrderedNodes(true,hardScaleFactor);
}
return res;
}
pair<CrossSection,CrossSection> Node::calcDipandPS(Energy scale){
return dipole()->dipandPs(sqr(scale),deepHead()->MH()->largeNBasis());
}
CrossSection Node::calcPs(Energy scale){
return dipole()->ps(sqr(scale),deepHead()->MH()->largeNBasis());
}
CrossSection Node::calcDip(Energy scale){
return dipole()->dip(sqr(scale));
}
IBPtr Node::clone() const {
return new_ptr(*this);
}
IBPtr Node::fullclone() const {
return new_ptr(*this);
}
void Node::persistentOutput(PersistentOStream & os) const {
os <<
thexcomb<<
thenodeMEPtr<<
thedipol<<
thechildren<<
theparent<<
theProjectors<<
theDeepHead<<
theCutStage<<
clustersafer<<
ounit(theRunningPt,GeV)<<
theSubtractedReal<<
theVirtualContribution<<
theMergingHelper;
}
void Node::persistentInput(PersistentIStream & is, int) {
is >>
thexcomb>>
thenodeMEPtr>>
thedipol>>
thechildren>>
theparent>>
theProjectors>>
theDeepHead>>
theCutStage>>
clustersafer>>
iunit(theRunningPt,GeV)>>
theSubtractedReal>>
theVirtualContribution>>
theMergingHelper;
}
// *** 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<Node,Interfaced> describeHerwigNode("Herwig::Node", "HwDipoleShower.so");
void Node::Init() {
static ClassDocumentation<Node> documentation("There is no documentation for the Node class");
}
diff --git a/Shower/Dipole/Merging/Node.h b/Shower/Dipole/Merging/Node.h
--- a/Shower/Dipole/Merging/Node.h
+++ b/Shower/Dipole/Merging/Node.h
@@ -1,242 +1,241 @@
// -*- C++ -*-
#ifndef Herwig_Node_H
#define Herwig_Node_H
//
// This is the declaration of the Node class.
//
//#include "Node.fh"
#include "MergingFactory.fh"
#include "Merger.h"
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/Shower/Dipole/Base/DipoleEventRecord.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Config/Pointers.h"
#include <vector>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the Node class.
*
* @see \ref NodeInterfaces "The interfaces"
* defined for Node.
*/
/**
* Define the SafeClusterMap type map<pair<pair<emitter, emmision>, spectator >
* , pair<first-clustering, second-clustering> >
*/
typedef map<pair<pair<int, int>, int >, pair<bool, bool> > SafeClusterMap;
class Node : public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
///The default constructor. Do not use!
Node();
// another constructor for first nodes
Node(Ptr<MatchboxMEBase>::ptr nodeME, int cutstage, Ptr<Merger>::ptr mh);
// another constructor for underlying nodes
Node(Ptr<Node>::ptr deephead,
Ptr<Node>::ptr head,
Ptr<SubtractionDipole>::ptr dipol,
Ptr<MatchboxMEBase>::ptr nodeME,
int cutstage);
/// The destructor.
virtual ~Node();
//@}
public:
// get children from vector<Ptr<MatchboxMEBase>::ptr>
void birth(vector<Ptr<MatchboxMEBase>::ptr> vec);
/// recursive setXComb. proStage is the number of clusterings
/// before the projectors get filled.
void setXComb(tStdXCombPtr xc);
/// calculate the dipole and ps approximation
pair<CrossSection, CrossSection> calcDipandPS(Energy scale);
/// calculate the ps approximation
CrossSection calcPs(Energy scale);
/// calculate the dipole
CrossSection calcDip(Energy scale);
/// recursive flush caches and clean up XCombs.
void flushCaches();
/// recursive clearKinematics
void clearKinematics();
/// recursive setKinematics
void setKinematics();
/// recursive generateKinematics using tilde kinematics of the dipoles
bool generateKinematics(const double *r, int stage, Energy2 shat);
/// generate the kinamatics of the first node
void firstgenerateKinematics(const double *r, int stage);
//return the ME
const Ptr<MatchboxMEBase>::ptr nodeME()const;
//return the node ME
Ptr<MatchboxMEBase>::ptr nodeME();
//return the parent Node
Ptr<Node>::ptr parent()const {return theparent;}
/// vector of children nodes created in birth
vector< Ptr<Node>::ptr > children()const {return thechildren;}
//pick a random child (flat)
Ptr<Node>::ptr randomChild();
/// true if all children show scales above pt
bool allAbove(Energy pt);
/// true if the node is in the history of other.
bool isInHistoryOf(Ptr<Node>::ptr other);
/// legsize of the node ME
int legsize() const;
/// set the first node (first men). only use in factory
void deepHead(Ptr<Node>::ptr deephead) {theDeepHead = deephead;}
/// return the first node
Ptr<Node>::ptr deepHead() const {return theDeepHead;}
/// insert nodes to projector vector
void Projector(double a, Ptr<Node>::ptr pro) {
pair<double, Ptr<Node>::ptr> p;
p.first = a;
p.second = pro;
theProjectors.push_back(p);
}
/// insert nodes to projector vector
vector< pair <double , Ptr<Node>::ptr > > Projector() {return theProjectors;}
/// returns the dipol of the node.
Ptr<SubtractionDipole>::ptr dipole() const;
/// set the xcomb of the node
void xcomb(StdXCombPtr xc) { thexcomb = xc;}
/// return the xcomb
StdXCombPtr xcomb() const {return thexcomb;}
/// return the current running pt
Energy runningPt(){return theRunningPt;}
/// set the current running pt
void runningPt(Energy x){theRunningPt=x;}
/// return the cut stage to cut on merging pt in generate kinematics
int cutStage() const {return theCutStage;}
/// get the clustersafe map for this node
SafeClusterMap clusterSafe() const {return clustersafer;}
/// get a vector of the next nodes, ordered in pt (and in parton shower phace space)
vector<Ptr<Node>::ptr> getNextOrderedNodes(bool normal=true, double hardscalefactor=1.);
//true if the node is in shower history for a given pt
bool inShowerPS(Energy hardpt);
//get the history
Ptr<Node>::ptr getHistory(bool normal=true, double hardscalefactor=1.);
//true if node correspond to a subtracted real.
- bool subtractedReal() {return theSubtractedReal;}
+ bool subtractedReal() const {return theSubtractedReal;}
/// set if node correspont to a subtracted real.
void subtractedReal(bool x) { theSubtractedReal = x;}
//true if node correspond to a virtual contribution.
- bool virtualContribution() { return theVirtualContribution ;}
+ bool virtualContribution() const { return theVirtualContribution ;}
/// set if node correspont to a virtual contribution.
void virtualContribution(bool x) {theVirtualContribution = x;}
//pointer to the merging helper
Ptr<Merger>::ptr MH()const{return theMergingHelper;}
/// set the merging helper
void MH(Ptr<Merger>::ptr a){theMergingHelper=a;}
Energy pT()const{return dipole()->lastPt();}
private:
- //the xcomb of the node
- StdXCombPtr thexcomb;
/// the Matrixelement representing this node.
Ptr<MatchboxMEBase>::ptr thenodeMEPtr;
/// the dipol used to substract
/// and generate kinematics using tilde kinematics
Ptr<SubtractionDipole>::ptr thedipol;
- /// vector of the children node
- vector< Ptr<Node>::ptr > thechildren;
/// the parent node
Ptr<Node>::ptr theparent;
- /// The nodes of the projection stage.
- vector< pair <double, Ptr<Node>::ptr> > theProjectors;
/// The godfather node of whole tree.(Firstnode)
Ptr<Node>::ptr theDeepHead;
/**
* The CutStage is number of clusterings which are possible without
* introducing a merging scale to cut away singularities.
* -> subtracted MEs have the CutStage 1.
* -> virtual and normal tree level ME get 0.
*/
int theCutStage;
- /// For [[Emitter, Emission], Spectator] the mapped pair gives
- /// information if the first and the second cluster is safe.
- SafeClusterMap clustersafer;
- /// the current running pt
- Energy theRunningPt;
/// tell if node belongs to an ordered history
bool isOrdered;
/// flag to tell if node is subtracted real
bool theSubtractedReal;
/// flag to tell if node is virtual contribution
bool theVirtualContribution;
/// the merging helper
Ptr<Merger>::ptr theMergingHelper;
-
+ //the xcomb of the node
+ StdXCombPtr thexcomb;
+ /// vector of the children node
+ vector< Ptr<Node>::ptr > thechildren;
+ /// the current running pt
+ Energy theRunningPt;
+ /// The nodes of the projection stage.
+ vector< pair <double, Ptr<Node>::ptr> > theProjectors;
+ /// For [[Emitter, Emission], Spectator] the mapped pair gives
+ /// information if the first and the second cluster is safe.
+ SafeClusterMap clustersafer;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Node & operator=(const Node &);
};
}
#endif /* Herwig_Node_H */