Page MenuHomeHEPForge
Files F7878338

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/Merging/Merger.cc b/DipoleShower/Merging/Merger.cc
--- a/DipoleShower/Merging/Merger.cc
+++ b/DipoleShower/Merging/Merger.cc
@@ -1,1504 +1,1506 @@
// -*- 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 "MFactory.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "Herwig/DipoleShower/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() {
theNf=5;
minusL=false;
Unlopsweights=true;
theCMWScheme=true;
theSmearing=0.;
theDipoleShowerHandler=
Ptr<DipoleShowerHandler>::ptr();
isUnitarized=true;
isNLOUnitarized=true;
defMERegionByJetAlg=false;
theOpenInitialStateZ=false;
theChooseHistory=8;
xiRenME=1.;
xiFacME=1.;
xiRenSh=1.;
xiFacSh=1.;
xiQSh=1.;
theGamma=1.;
ee_ycut=1000000;
pp_dcut=1000000;
therenormscale=0.*GeV;
theIRSafePT=1.*GeV;
theMergePt=3.*GeV;
theCentralMergePt=3.*GeV;
}
Merger::~Merger() {}
IBPtr Merger::clone() const {
return new_ptr(*this);
}
IBPtr Merger::fullclone() const {
return new_ptr(*this);
}
double Merger::reweightCKKWBornGamma(NPtr Node,bool fast){
double weightCL=0.;
weight=1.;
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), 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;Node->nodeME()->projectorStage(1);
weightCL=2.*int(Node->children().size());
BornCL=CLNode-> getHistory(false,xiQSh);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
return 0.;
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 0.;
}
double res=0.;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsNLO();
if(weight!=0.){
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node,fast);
if (!fillProjector(projectedscale))weight=0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.&&weightCL==0.)return 0.;
res+=weight*matrixElementWeight(startscale,Node,(!maxMulti&&!projected)?theGamma:1.);
}
if(CLNode&&theGamma!=1.){
Energy startscale=CKKW_StartScale(BornCL);
Energy projectedscale=startscale;
fillHistory( startscale, BornCL, CLNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weightCL*=history.back().weight*alphaReweight()*pdfReweight();
double diff=0;
Node->nodeME()->factory()->setAlphaParameter(1.);
diff-=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
Node->nodeME()->factory()->setAlphaParameter(theGamma);
string alp=(CLNode->dipol()->aboveAlpha()?"Above":"Below");
diff+=weightCL*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
Node->nodeME()->factory()->setAlphaParameter(1.);
res+=diff;
}
return res;
}
double Merger::reweightCKKWBornStandard(NPtr Node,bool fast){
if (!Node->children().empty()) {
PVector particles;
for (size_t i=0;i<Node->nodeME()->mePartonData().size();i++){
Ptr<ThePEG::Particle>::ptr p =Node->nodeME()->mePartonData()[i]->produceParticle(Node->nodeME()->lastMEMomenta()[i]);
particles.push_back(p);
}
if(!matrixElementRegion(particles, Node->children()[0]->dipol()->lastPt(), theMergePt))return 0.;
}
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (treefactory()->onlymulti()!=-1&&
treefactory()->onlymulti()!=int(Node->xcomb()->mePartonData().size()-Node->nodeME()->projectorStage()))
return 0.;
assert(Node->allAbove(mergePt()-0.1*GeV));
if (!Born->xcomb()->willPassCuts()){
return 0.;
}
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
return weight*matrixElementWeight(startscale,Node,1.);
}
double Merger::reweightCKKWVirtualStandard(NPtr Node,bool fast){
NPtr Born= Node-> getHistory(true,xiQSh);
if( Born!= Node){
if (UseRandom::rnd()<.5){
weight=-2.;projected=true;Node->nodeME()->projectorStage(1);
}else{
weight=2.;projected=false;Node->nodeME()->projectorStage(0);
}
}else{
weight=1.;projected=false;Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, Node,fast);
if (history.size()==1&&Node->children().size()!=0) {
cout<<"\n1-->"<<startscale/GeV<<" "<<weight;
}
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=Node->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double matrixElement=matrixElementWeight(startscale,Node);
double Bornweight=Node->nodeME()->lastBorndSigHatDR();
double unlopsweight =(-sumpdfReweightUnlops()
-sumalphaReweightUnlops()
-sumfillHistoryUnlops())
*Bornweight
*SM().alphaS()/(2.*ThePEG::Constants::pi);
//assert(unlopsweight==0.||history.size()!=1);
if (history.size()==1&&Node->children().size()!=0) {
cout<<"\n unlopsweight "<<unlopsweight<<" "<<matrixElement;
}
return weight*
as(startscale*xiRenSh)/SM().alphaS()*
(matrixElement+unlopsweight);
}
double Merger::reweightCKKWRealStandard(NPtr Node,bool fast){
bool allAbove=Node->allAbove(mergePt());
+ if(!Node->allAbove(theIRSafePT))return 0.;
if (allAbove)return reweightCKKWRealAllAbove(Node, fast);
if (UseRandom::rnd()<.5)
return 2.*reweightCKKWRealBelowSubReal( Node, fast);
return 2.*reweightCKKWRealBelowSubInt( Node, fast);
}
double Merger::reweightCKKWRealAllAbove(NPtr Node,bool fast){
NPtr Born= Node-> getHistory(true,xiQSh);
NPtr CLNode= Node->randomChild();
bool inhist=CLNode->isInHistoryOf(Born);
Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
if (!CLNode->allAbove(mergePt()))return 0.;
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double res= weight*as(startscale*xiRenSh)/SM().alphaS()*
(double)Node->children().size()*
((inhist?matrixElementWeight(startscale,Node):0.)
+CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn);
// cout<<"\nCLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn "<<CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
return res;
}
double Merger::reweightCKKWRealBelowSubReal(NPtr Node,bool fast){
NPtrVec children=Node->children();
Selector<NPtr> HistNodeSel;
Energy minScale=generator()->maximumCMEnergy();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->dipol()->lastPt()<minScale) {
minScale=(*child)->dipol()->lastPt();
HistNodeSel.insert(1.,*child);
}
}
NPtr HistNode=HistNodeSel.select(UseRandom::rnd());
NPtr Born=HistNode-> getHistory(false,xiQSh);
if( Born!= HistNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(1);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(0);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(0);
}
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, HistNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=HistNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
double sumPS=0;
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
if ((*child)->allAbove(mergePt()))
sumPS-=(*child)->calcPs(startscale*xiFacME);
}
- //double me=matrixElementWeight(startscale,Node);
- //cout<<"\nreweightCKKWRealBelowSubReal "<<me<<" "<<sumPS<<" "<<me/sumPS;
+ double me=matrixElementWeight(startscale,Node);
+ //cout<<"\nSubReal "<<Node->miniPt()/GeV<<" "<<me<<" "<<sumPS<<" "<<me/sumPS;
+
return weight*as(startscale*xiRenSh)/SM().alphaS()*
- (matrixElementWeight(startscale,Node)-sumPS);
+ (me-sumPS);
}
double Merger::reweightCKKWRealBelowSubInt(NPtr Node,bool fast){
NPtr CLNode= Node->randomChild();
NPtr Born=CLNode-> getHistory(false,xiQSh);
if( Born!= CLNode){
if (UseRandom::rnd()<.5){
weight=-2.; projected=true; Node->nodeME()->projectorStage(2);}
else{
weight= 2.; projected=false; Node->nodeME()->projectorStage(1);}
}else{
weight=1.; projected=false; Node->nodeME()->projectorStage(1);
}
if (!CLNode->allAbove(mergePt()))return 0.;
if (!Born->xcomb()->willPassCuts())return 0.;
Energy startscale=CKKW_StartScale(Born);
Energy projectedscale=startscale;
fillHistory( startscale, Born, CLNode,fast);
if (!fillProjector(projectedscale))return 0.;
Node->runningPt(projectedscale);
weight*=history.back().weight*alphaReweight()*pdfReweight();
if(weight==0.)return 0.;
bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);
- double res=0.;
+ pair<double,double> DipAndPs=make_pair(0.,0.);
if (Born==CLNode&&!CLNode->children().empty())
- res=-1.*CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn;
+ DipAndPs=make_pair(CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn,0.);
else
- res=CLNode->calcPsMinusDip(startscale*xiFacME);
+ DipAndPs=CLNode->calcDipandPS(startscale*xiFacME);
return weight*as(startscale*xiRenSh)/SM().alphaS()*
- (double)Node->children().size()*res;
+ (double)Node->children().size()*(DipAndPs.second-DipAndPs.first);
}
double Merger::matrixElementWeight(Energy startscale,NPtr Node,double diffAlpha){
double res;
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
res=DeepHead->nodeME()->dSigHatDR(false,diffAlpha)/nanobarn;
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
double Merger::matrixElementWeightWithLoops(Energy startscale,NPtr Node,bool fast){
double res=0.;
return res;
// The deephead should be calculated here.
NPtr DeepHead=Node;//->deepHead();
renormscale(startscale);
DeepHead->nodeME()->factory()->scaleChoice()->setXComb(DeepHead->xcomb());
DeepHead->nodeME()->setScale(sqr(startscale),sqr(startscale));
theCalculateInNode=false;
DeepHead->nodeME()->doOneLoopNoBorn();
res=DeepHead->nodeME()->dSigHatDR(fast)/nanobarn;
DeepHead->nodeME()->noOneLoopNoBorn();
renormscale(0.0*GeV);
theCalculateInNode=true;
return res;
}
bool Merger::fillProjector(Energy& prerunning){
if(history.begin()->node->deepHead()->nodeME()->projectorStage() == 0){
prerunning=(history.size()==1?xiQSh:1.)*history.back().scale;
return true;
}
for (History::iterator it=history.begin();it!=history.end();it++){
if (projectorStage((*it).node)&&history.begin()->node->deepHead()->nodeME()->projectorStage() != 0){
history.begin()->node->deepHead()->xcomb()->lastProjector((*it).node->xcomb());
prerunning=(it==history.begin()?xiQSh:1.)*(*it).scale;
return true;
}
}
return false;
}
double Merger::pdfReweight(){
double res=1.;
for(int side=0;side!=2;side++){
if(history[0].node->xcomb()->mePartonData()[side]->coloured()){
History::iterator it=history.begin();
for (;it+1!=history.end();it++){
res *= pdfratio((*it).node, (*it).scale,xiFacSh*((*it).node->dipol()->lastPt()), side);
}
res*=pdfratio(history.back().node,history.back().scale ,history[0].scale*xiFacME, side);
}
}
return res;
}
double Merger::alphaReweight(){
double res=1.;
Energy Q_R=xiRenME*history[0].scale;
res *= pow(as(Q_R) / SM().alphaS(), history[0].node->nodeME()->orderInAlphaS());
res *= pow(history[0].node->deepHead()->xcomb()->eventHandler().SM().alphaEMPtr()->value(history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED())/ SM().alphaEMMZ(), history[0].node->nodeME()->orderInAlphaEW());
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
if ((*it).node->parent()){
Energy q_i=xiRenSh* (*it).node->dipol()->lastPt();
res *= as(q_i)/ SM().alphaS()
*(theCMWScheme?(1.+((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*5.)*as(q_i))/2./Constants::pi):1.);
}
}
return res;
}
void Merger::fillHistory(Energy scale, NPtr Begin, NPtr EndNode,bool fast){
history.clear();
double sudakov0_n=1.;
history.push_back(HistoryStep(Begin,sudakov0_n,scale));
scale*=xiQSh;
if (Begin->parent()||!isUnitarized) {
while (Begin->parent() && (Begin != EndNode)) {
if (!dosudakov(Begin,scale, Begin->dipol()->lastPt(),sudakov0_n,fast)){
history.push_back(HistoryStep(Begin->parent(),0.,scale));
}
scale=Begin->dipol()->lastPt();
history.push_back(HistoryStep(Begin->parent(),sudakov0_n,Begin->dipol()->lastPt()));
Begin = Begin->parent();
}
Energy notunirunning=scale;
if (!isUnitarized&&N()+N0() > Begin->deepHead()->nodeME()->lastMEMomenta().size()) {
if (!dosudakov(Begin,notunirunning,mergePt(),sudakov0_n,fast)){
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*xiFacME;
Energy beam2Scale=history[0].scale*xiFacME;
History::iterator it=history.begin();
for (;it!=history.end();it++){
History::iterator ittmp=it;
ittmp++;
if((*it).node->xcomb()->mePartonData()[0]->coloured()&&(*it).node->nodeME()->lastX1()>0.99)return 0.;
if((*it).node->xcomb()->mePartonData()[1]->coloured()&&(*it).node->nodeME()->lastX2()>0.99)return 0.;
if (ittmp!=history.end()){
if((*it).node->nodeME()->lastX1()<0.00001)return 0.;
if((*it).node->nodeME()->lastX2()<0.00001)return 0.;
if ((*it).node->dipol()->bornEmitter() == 0 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
theNf,
history[0].scale);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornEmitter() == 1 ){
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
theNf,
history[0].scale);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 0 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().first->dataPtr(),
(*it).node->nodeME()->lastPartons().first->dataPtr(),
(*it).node->xcomb()->partonBins().first->pdf(),
beam1Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX1(),
theNf,
history[0].scale);
//pdfratio((*it).node, beam1Scale, sqrt((*it).node->dipol()->lastPt()), 1);
beam1Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
if ((*it).node->dipol()->bornSpectator() == 1 &&(*it).node->dipol()->bornEmitter() >1){//
res +=pdfUnlops((*it).node->nodeME()->lastParticles().second->dataPtr(),
(*it).node->nodeME()->lastPartons().second->dataPtr(),
(*it).node->xcomb()->partonBins().second->pdf(),
beam2Scale,
((*it).node->dipol()->lastPt()),
(*it).node->nodeME()->lastX2(),
theNf,
history[0].scale);
//pdfratio((*it).node, beam2Scale , sqrt((*it).node->dipol()->lastPt()), 2);
beam2Scale=((*it).node->dipol()->lastPt())*xiFacSh;
}
}
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[0]->coloured()){
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().first->dataPtr(),
history.back().node->nodeME()->lastPartons().first->dataPtr(),
history.back().node->xcomb()->partonBins().first->pdf(),
beam1Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX1(),
theNf,
history[0].scale);
}
if (history[0].node->deepHead()->xcomb()->mePartonData()[1]->coloured()) {
res +=pdfUnlops(history.back().node->nodeME()->lastParticles().second->dataPtr(),
history.back().node->nodeME()->lastPartons().second->dataPtr(),
history.back().node->xcomb()->partonBins().second->pdf(),
beam2Scale,
history[0].scale*xiFacME,
(history.back()).node->nodeME()->lastX2(),
theNf,
history[0].scale);
//history[0].node->deepHead()->nodeME()->pdf2(sqr(beam2Scale))/history[0].node->deepHead()->nodeME()->pdf2(sqr(history[0].scale));
}
return res;
}
double Merger::pdfUnlops(tcPDPtr particle,tcPDPtr parton,tcPDFPtr pdf,Energy running, Energy next,double x,int nlp,Energy fixedScale) {
//copied from PKOperator
double res=0.;
int number=10;//*int((max(running/next,next/running)));
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.)*theNf + 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,
history[0].scale);
}
// dsig is calculated with fixed alpha_s
for (History::iterator it=history.begin();(it+1)!=history.end();it++){
assert((*it).node->parent());
res +=alphasUnlops((*it).node->dipol()->lastPt()*xiRenSh ,history[0].scale*xiRenME);
}
return res;
}
double Merger::sumfillHistoryUnlops(){
double res=0.;
for (History::iterator it = history.begin(); (it+1) != history.end();it++){
doUNLOPS((*it).node,(it == history.begin()?xiQSh:1.)*(*it).scale , (*it).node->dipol()->lastPt() , history[0].scale, res);
}
return res;
}
Ptr<MFactory>::ptr Merger::treefactory(){return theTreeFactory;}
bool Merger::reweightCKKWSingle(Ptr<MatchboxXComb>::ptr SX, double & res,bool fast) {
assert(!fast);
Ptr<StandardEventHandler>::ptr eH =
dynamic_ptr_cast<Ptr<StandardEventHandler>::ptr>(generator()->eventHandler());
if (!eH->didEstimate()||fast) {
history.clear();
}
theNf=5;//TODO
if (!SX) return true;
assert(SX->eventHandlerPtr());
Ptr<MatchboxMEBase>::ptr ME = dynamic_ptr_cast<Ptr<MatchboxMEBase>::ptr>(SX->matchboxME());
if (!ME) return true;
if (theFirstNodeMap.find(ME)==theFirstNodeMap.end()) {
cout<<"\nnot in map:"<<theFirstNodeMap.size();
}
NPtr Node = theFirstNodeMap[SX->matchboxME()];
Ptr<MFactory>::ptr factory = dynamic_ptr_cast<Ptr<MFactory>::ptr>(Node->nodeME()->factory());
assert(factory);
if (!Node) return true;
NPtr MENode = Node;
Ptr<AlphaEMBase>::transient_pointer alphaEM = SX->eventHandler().SM().alphaEMPtr();
assert(DSH()->hardScaleIsMuF());
xiRenME=ME->renormalizationScaleFactor();
xiFacME=ME->factorizationScaleFactor();
xiRenSh=DSH()->renormalizationScaleFactor();
xiFacSh=DSH()->factorizationScaleFactor();
xiQSh=DSH()->hardScaleFactor();
if(Node->deepHead()->subtractedReal()){
res*=reweightCKKWRealStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(ME->oneLoopNoBorn()){
res*=reweightCKKWVirtualStandard(Node);
theCurrentMaxLegs=maxLegsNLO();
}else if(theGamma!=1.){
res*=reweightCKKWBornGamma(Node,fast);
theCurrentMaxLegs=maxLegsLO();
}else{
res*=reweightCKKWBornStandard(Node,fast);
theCurrentMaxLegs=maxLegsLO();
}
SX->lastCentralScale(sqr(Node->runningPt()));
SX->lastShowerScale(sqr(Node->runningPt()));
if(SX->lastProjector()){
SX->lastProjector()->lastCentralScale(sqr(Node->runningPt()));
SX->lastProjector()->lastShowerScale(sqr(Node->runningPt()));
}
renormscale(0.0*GeV);
if (res == 0.){
history.clear();
return false;
}
cleanup(MENode);
cleanup(Node);
DSH()->eventRecord().clear();
SX->subProcess(SubProPtr());
if (!fast) {
history.clear();
}
return true;
}
//----------------------------------Reviewed--------------------------------------//
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().prepare(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 (size_t i=0;i<Born->nodeME()->mePartonData().size();i++){
if (!Born->nodeME()->mePartonData()[i]->coloured())continue;
for (size_t j=i+1;j<Born->nodeME()->mePartonData().size();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.)*theNf;
return (betaZero*log(sqr(fixedScale/next)))+
(theCMWScheme?(((3.*(67./18.-1./6.*Constants::pi*Constants::pi)-5./9.*theNf))):0.);
}
double Merger::pdfratio(NPtr Born,Energy & nominator_scale, Energy denominator_scale,int side){
StdXCombPtr bXc = Born->xcomb();
assert (bXc->mePartonData()[side]->coloured());
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,bool fast) {
CKKW_PrepareSudakov(Born, running);
for ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
sudakov0_n*=singlesudakov( dip, next,running,make_pair(true,false), fast );
sudakov0_n*=singlesudakov( dip, next,running,make_pair(false,true), fast );
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 ( list<DipoleChain>::iterator chain = DSH()->eventRecord().chains().begin() ;
chain != DSH()->eventRecord().chains().end() ; chain++ ) {
for ( list<Dipole>::iterator dip = (*chain).dipoles().begin() ; dip != (*chain).dipoles().end() ; ++dip ) {
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::projectorStage(NPtr Born){
return (Born->deepHead()->nodeME()->projectorStage() ==
int((Born->deepHead()->nodeME()->mePartonData().size()
- Born->nodeME()->mePartonData().size())));
}
void Merger::cleanup(NPtr Born) {
DSH()->eventRecord().clear();
if(!Born->xcomb()->subProcess())return;
ParticleVector vecfirst = Born->xcomb()->subProcess()->incoming().first->children();
for ( ParticleVector::const_iterator it = vecfirst.begin() ; it != vecfirst.end() ; it++ )
Born->xcomb()->subProcess()->incoming().first->abandonChild(*it);
ParticleVector vecsecond = Born->xcomb()->subProcess()->incoming().second->children();
for ( ParticleVector::const_iterator it = vecsecond.begin() ; it != vecsecond.end() ; it++ )
Born->xcomb()->subProcess()->incoming().second->abandonChild(*it);
Born->xcomb()->subProcess(SubProPtr());
}
double Merger::singlesudakov(list<Dipole>::iterator dip ,Energy next,Energy running,pair<bool,bool> conf ,bool fast){
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());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
assert (! isnan(dScale/GeV ) );
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,fast);
}
}
return res;
}
double Merger::singleUNLOPS(list<Dipole>::iterator 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());
pair<GeneratorMap2::iterator,GeneratorMap2::iterator> gens = DSH()->generators().equal_range(candidate.index());
for ( GeneratorMap2::iterator gen = gens.first; gen != gens.second; ++gen ) {
if ( !(gen->first == candidate.index()) )
continue;
Energy dScale = gen->second->splittingKinematics()->dipoleScale(emitter->momentum(),spectator->momentum());
assert (! isnan(dScale/GeV ) );
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));};
+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,int st){
theFirstNodeMap[me]->setXComb(xc, st);
}
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,me->lastXCombPtr()->lastSHat());
if (theFirstNodeMap[me]->cutStage()==0 ){
bool inAlphaPS=false;
NPtrVec children=theFirstNodeMap[me]->children();
for (NPtrVec::iterator child = children.begin();
child != children.end(); child++){
treefactory()->setAlphaParameter(theGamma);
inAlphaPS|=(*child)->dipol()->aboveAlpha();
treefactory()->setAlphaParameter(1.);
}
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first&&!inAlphaPS)return false;
}
}
if (theFirstNodeMap[me]->cutStage()==1 ){
SafeClusterMap temp=theFirstNodeMap[me]->clusterSafe();
for(SafeClusterMap::iterator
it=temp.begin();
it!=temp.end();++it){
if (!it->second.first && !it->second.second)return false;
}
}
return true;
}
bool Merger::calculateInNode() const{
return theCalculateInNode;
}
void Merger::fillProjectors(Ptr<MatchboxMEBase>::ptr me){
for (unsigned int i = 0; i < (theFirstNodeMap[me]->Projector()).size(); ++i) {
me->lastXCombPtr()->projectors().insert(
(theFirstNodeMap[me]->Projector())[i].first,
(theFirstNodeMap[me]->Projector())[i].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 particles,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(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[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 (particles[em]->coloured())
input_particles.push_back(fastjet::PseudoJet(particles[em]->momentum().x()/GeV,
particles[em]->momentum().y()/GeV,
particles[em]->momentum().z()/GeV,
particles[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();
}
*/
//assert(false);
Energy ptx=1000000.*GeV;
bool foundwinnerpt=false;
//FF
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (em==spe||emm==spe) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Energy scale = sqrt(2.*(emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom));
double y = emissionmom*emittermom / (emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom);
double z = emittermom*spectatormom / (emittermom*spectatormom + emissionmom*spectatormom);
if (abs(scale * sqrt(y*z*(1.-z))-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
// if(scale * sqrt(y*z*(1.-z))<optVeto&&winnerScale>optVeto)cout<<"\nFF "<<(scale * sqrt(y*z*(1.-z))/GeV);
ptx =min(ptx,scale * sqrt(y*z*(1.-z)));
}
}
}
//FI
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t em=2; em < particles.size();em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==emm) continue;
if (!(particles[em]->id()==-particles[emm]->id()||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Lorentz5Momentum Emsp=emittermom+emissionmom-spectatormom;
Energy scale = sqrt(-1.*Emsp*Emsp);
double x = (- emissionmom*emittermom + emissionmom*spectatormom + emittermom*spectatormom) /(emittermom*spectatormom + emissionmom*spectatormom);
double z = emittermom*spectatormom / (emittermom*spectatormom + emissionmom*spectatormom);
// if(scale * sqrt(z*(1.-z)*(1.-x)/x)<optVeto&&winnerScale>optVeto)cout<<"\nFI "<<(scale * sqrt(z*(1.-z)*(1.-x)/x)/GeV);
if (abs(scale *sqrt(z*(1.-z)*(1.-x)/x)-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
if((z*(1.-z)*(1.-x)/x)>0.)
ptx =min(ptx,scale * sqrt(z*(1.-z)*(1.-x)/x));
}
}
}
//IF
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
for(size_t spe=2; spe < particles.size();spe++){
if (!particles[spe]->coloured()) continue;
if (emm==spe) continue;
if (!(particles[em]->id()>6|| particles[em]->id()==particles[emm]->id() ||particles[emm]->id()>6))continue;
// assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Lorentz5Momentum Emsp=-emittermom+emissionmom+spectatormom;
Energy scale = sqrt(-1.*Emsp*Emsp);
double x = (-emissionmom*spectatormom + emittermom*spectatormom + emittermom*emissionmom) /(emittermom*emissionmom + emittermom*spectatormom);
double u = emittermom*emissionmom / (emittermom*emissionmom + emittermom*spectatormom);
if (abs(scale *sqrt(u*(1.-u)*(1.-x)/x)-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
if((u*(1.-u)*(1.-x)/x)>0.)
ptx =min(ptx,scale * sqrt(u*(1.-u)*(1.-x)/x));
}
}
}
//II
for(size_t em=0; em < 2;em++){
if (!particles[em]->coloured()) continue;
for(size_t spe=0; spe < 2;spe++){
if (!particles[spe]->coloured()) continue;
for(size_t emm=2; emm < particles.size();emm++){
if (!particles[emm]->coloured()) continue;
if (em==spe) continue;
if (!(particles[em]->id()>6||
particles[em]->id()==particles[emm]->id() ||
particles[emm]->id()>6))continue;
//assert(false);
Lorentz5Momentum emittermom = particles[em]->momentum();
Lorentz5Momentum emissionmom = particles[emm]->momentum();
Lorentz5Momentum spectatormom = particles[spe]->momentum();
Lorentz5Momentum Emsp=emittermom-emissionmom+spectatormom;
Energy scale = sqrt(Emsp*Emsp);
double x = (emittermom*spectatormom - emittermom*emissionmom - spectatormom*emissionmom)/(emittermom*spectatormom);
double v = (emittermom*emissionmom)/(emittermom*spectatormom);
if (abs(scale * sqrt(v*(1.-x-v)/x)-winnerScale)<0.0001*GeV) {
foundwinnerpt=true;
}
if ((v*(1.-x-v)/x)>0.)
ptx =min(ptx, scale * sqrt(v*(1.-x-v)/x));
}
}
}
// cout<<"\n"<<cutscale/GeV<< " "<<foundwinnerpt<<" "<<ptx/GeV;
assert(foundwinnerpt);
return (ptx>cutscale) ;
}
// 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 << minusL<< Unlopsweights<<
theCMWScheme<< projected<<
isUnitarized<< isNLOUnitarized<<
defMERegionByJetAlg<<theOpenInitialStateZ<<
theChooseHistory<<
theN0<< theOnlyN<<
theN<< theM<<
xiRenME<< xiFacME<<
xiRenSh<< xiFacSh<<
xiQSh<< theNf<<
weight<<weightCB<<theGamma<<ee_ycut<<pp_dcut<<theSmearing<<ounit(therenormscale, GeV)<<ounit(theIRSafePT, GeV)<<ounit(theMergePt, GeV)<<ounit(theCentralMergePt, GeV)<<theMergingJetFinder<<theLargeNBasis<<theDipoleShowerHandler<< theTreeFactory<<theFirstNodeMap ;
}
void Merger::persistentInput(PersistentIStream & is, int) {
is >> minusL>> Unlopsweights>>
theCMWScheme>> projected>>
isUnitarized>> isNLOUnitarized>>
defMERegionByJetAlg>>theOpenInitialStateZ>>
theChooseHistory>>
theN0>> theOnlyN>>
theN>> theM>>
xiRenME>> xiFacME>>
xiRenSh>> xiFacSh>>
xiQSh>> theNf>>
weight>>weightCB>>theGamma>>ee_ycut>>pp_dcut>>theSmearing>>iunit(therenormscale, GeV)>>iunit(theIRSafePT, GeV)>>iunit(theMergePt, GeV)>>iunit(theCentralMergePt, GeV)>>theMergingJetFinder>>theLargeNBasis>>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>
interfaceminusL ("minusL","",&Merger::minusL, false, false, false);
static SwitchOption interfaceminusLYes
(interfaceminusL,"Yes","",true);
static SwitchOption interfaceminusLNo
(interfaceminusL,"No","",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,MFactory> interfaceMergerHelper
("MFactory",
"",
&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);
static Parameter<Merger, int> interfaceaddLOLegs("addLOLegs", "Set the number of additional jets to consider.", &Merger::theN, 0, 0,
0, false, false, Interface::lowerlim);
static Parameter<Merger, int> interfaceaddNLOLegs("addNLOLegs",
"Set the number of virtual corrections to consider. -1 is default for no virtual correction.", &Merger::theM, -1, -1, 0, false, false,
Interface::lowerlim);
}
diff --git a/DipoleShower/Merging/Node.cc b/DipoleShower/Merging/Node.cc
--- a/DipoleShower/Merging/Node.cc
+++ b/DipoleShower/Merging/Node.cc
@@ -1,613 +1,627 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Node class.
//
#include "Node.h"
#include "MFactory.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 deepprostage, int cutstage) {
nodeME->maxMultCKKW(1);
nodeME->minMultCKKW(0);
theDeepHead = this;
thenodeMEPtr = nodeME;
theDeepProStage = deepprostage;
theCutStage = cutstage;
theSubtractedReal=false;
theVirtualContribution=false;
}
Node::Node(Ptr<Node>::ptr deephead, Ptr<Node>::ptr head, Ptr<SubtractionDipole>::ptr dipol, Ptr<MatchboxMEBase>::ptr nodeME,
int deepprostage, int cutstage) {
theDeepHead = deephead;
theparent = head;
thedipol = dipol;
thenodeMEPtr = nodeME;
theDeepProStage = deepprostage;
theCutStage = cutstage;
theSubtractedReal=false;
theVirtualContribution=false;
}
Node::~Node() { }
Ptr<SubtractionDipole>::ptr Node::dipol() 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;
}
Ptr<Node>::ptr Node::randomChild() {
return thechildren[(int)(UseRandom::rnd() * thechildren.size())];
}
+Energy Node::miniPt() const{
+ Energy res=1000000000*GeV;
+ for (vector<Ptr<Node>::ptr>::const_iterator it = thechildren.begin(); it != thechildren.end(); it++) {
+ res=min(res,(*it)->dipol()->lastPt());
+ }
+ return res;
+
+}
+
+
+
bool Node::allAbove(Energy pt){
for (vector<Ptr<Node>::ptr>::iterator it = thechildren.begin(); it != thechildren.end(); it++) {
if((*it)->dipol()->lastPt()<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 ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->underlyingBornME()->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r = thechildren[i]->dipol()->reweights().begin() ; r != thechildren[i]->dipol()->reweights().end() ;
++r ) {
(**r).flushCaches();
}
if ( thechildren[i]->xcomb() ) thechildren[i]->xcomb()->clean();
thechildren[i]->nodeME()->flushCaches();
thechildren[i]->flushCaches();
}
}
void Node::setKinematics() {
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
thechildren[i]->dipol()->setKinematics();
thechildren[i]->nodeME()->setKinematics();
thechildren[i]->setKinematics();
}
}
void Node::clearKinematics() {
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
thechildren[i]->nodeME()->clearKinematics();
thechildren[i]->dipol()->clearKinematics();
thechildren[i]->clearKinematics();
}
}
bool Node::generateKinematics(const double *r, int stage, Energy2 shat) {
isOrdered = false;
isthissafe=true;
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
if ( !thechildren[i]->dipol()->generateKinematics(r) ) cout << "stop";
if ( dipol()->lastPt() < thechildren[i]->dipol()->lastPt() && parent()->ordered() ) {
isOrdered = true;
deepHead()->setOrderedSteps(stage + 1);
}
thechildren[i]->generateKinematics(r, stage + 1, thechildren[i]->xcomb()->lastSHat());
isthissafe = (isthissafe && thechildren[i]->dipol()->lastPt() >=deepHead()->MH()->mergePt());
}
return isthissafe;
}
void Node::firstgenerateKinematics(const double *r, int stage, Energy2 shat) {
flushCaches();
//Set here the new merge Pt for the next phase space point.( Smearing!!!)
MH()->smeareMergePt();
isOrdered = true;
// if we consider the hard scale to play a role in the steps we must change here to 0.
setOrderedSteps(1);
this->orderedSteps();
clustersafer.clear();
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
bool ifirst = true;
bool isecond = true;
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
if ( !thechildren[i]->dipol()->generateKinematics(r) ) cout << "stop";
isecond = thechildren[i]->generateKinematics(r, stage + 1, thechildren[i]->xcomb()->lastSHat());
ifirst = (thechildren[i]->dipol()->lastPt() >= deepHead()->MH()->mergePt());
pair<pair<int, int>, int> EmitEmisSpec = make_pair(make_pair(thechildren[i]->dipol()->realEmitter(), thechildren[i]->dipol()->realEmission()),
thechildren[i]->dipol()->realSpectator());
clustersafer.insert(make_pair(EmitEmisSpec, make_pair(ifirst, isecond)));
}
}
void Node::setXComb(tStdXCombPtr xc, int proStage) {
if ( !parent() ) this->xcomb(xc);
for ( unsigned int i = 0 ; i < thechildren.size() ; ++i ) {
if ( !thechildren[i]->xcomb() ) {
thechildren[i]->xcomb(thechildren[i]->dipol()->makeBornXComb(xc));
assert(thechildren[i]->xcomb());
thechildren[i]->xcomb()->head(xc);
if ( !thechildren[i]->dipol()->lastXCombPtr() ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
}
thechildren[i]->setXComb(thechildren[i]->xcomb(), (proStage - 1));
} else {
if ( !(thechildren[i]->dipol()->lastXCombPtr()->lastScale() == thechildren[i]->xcomb()->lastScale()) ) {
thechildren[i]->dipol()->setXComb(thechildren[i]->xcomb());
}
if ( thechildren[i]->xcomb()->head() != xc ) thechildren[i]->xcomb()->head(xc);
thechildren[i]->setXComb(thechildren[i]->xcomb(), (proStage - 1));
}
}
}
void Node::birth(vector<Ptr<MatchboxMEBase>::ptr> vec) {
vector<Ptr<SubtractionDipole>::ptr> dipoles =
nodeME()->getDipoles(DipoleRepository::dipoles(
nodeME()->factory()->dipoleSet()), vec,true);
for ( unsigned int j = 0 ; j < dipoles.size() ; ++j ) {
dipoles[j]->doSubtraction();
Ptr<Node>::ptr node = new_ptr(Node(theDeepHead,this, dipoles[j],
dipoles[j]->underlyingBornME(),
theDeepHead->DeepProStage(),
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 ( vector<Ptr<Node>::ptr>::const_iterator it = temp.begin() ; it != temp.end() ; ++it ) {
if(deepHead()->MH()->mergePt()>(*it)->dipol()->lastPt()){
res.clear();
// if any of the nodes is below the merging scale return empty vector
return res;
continue;
}
if (parent()&& normal){
if ( (*it)->dipol()->lastPt() < dipol()->lastPt() ){
continue;
}
}
if ( (*it)->children().size() != 0 ){
vector<Ptr<Node>::ptr> tempdown = (*it)->children();
for ( vector<Ptr<Node>::ptr>::const_iterator itd = tempdown.begin() ; itd != tempdown.end() ; ++itd ) {
if( (*itd)->dipol()->lastPt() > (*it)->dipol()->lastPt()&&(*it)->inShowerPS((*itd)->dipol()->lastPt()) ){
res.push_back(*it);
break;
}
}
}
else {
(*it)->nodeME()->factory()->scaleChoice()->setXComb((*it)->xcomb());
if ( sqr(hardScaleFactor)*(*it)->nodeME()->factory()->scaleChoice()->renormalizationScale()
>= (*it)->dipol()->lastPt() * (*it)->dipol()->lastPt()&&
(*it)->inShowerPS(hardScaleFactor*sqrt((*it)->nodeME()->factory()->scaleChoice()->renormalizationScale()))) {
res.push_back(*it);
}
}
}
return res;
}
bool Node::inShowerPS(Energy hardpT){
assert(deepHead()->MH()->largeNBasis());
double x=0.;
double z_=dipol()->lastZ();
// if (dipol()->lastPt()>50*GeV) {
// return false;
// }
string type;
// II
if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()<2){
type="II";
x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
double ratio = sqr(dipol()->lastPt()/dipol()->lastDipoleScale());
double x__ = z_*(1.-z_)/(1.-z_+ratio);
double v = ratio*z_ /(1.-z_+ratio);
if (dipol()->lastPt()>(1.-x) * dipol()->lastDipoleScale()/ (2.*sqrt(x)))return false;
assert(v< 1.-x__&&x > 0. && x < 1. && v > 0.);
if (deepHead()->MH()->openInitialStateZ()) {
return true;
}
}
// IF
if( dipol()->bornEmitter()<2&&dipol()->bornSpectator()>=2){
type="IF";
x =dipol()->bornEmitter()==0?xcomb()->lastX1():xcomb()->lastX2();
if (dipol()->lastPt()>dipol()->lastDipoleScale()* sqrt((1.- x)/x) /2.)return false;
if (deepHead()->MH()->openInitialStateZ()) {
return true;
}
}
// FI
if( dipol()->bornEmitter()>=2&&dipol()->bornSpectator()<2){
type="FI";
double lastx=dipol()->bornSpectator()==0?xcomb()->lastX1():xcomb()->lastX2();
if (dipol()->lastPt()>dipol()->lastDipoleScale()*sqrt((1.-lastx)/lastx) /2.)return false;
}
// FF
if( dipol()->bornEmitter()>=2&&dipol()->bornSpectator()>=2){
type="FF";
if (dipol()->lastPt()>dipol()->lastDipoleScale()/2.)return false;
}
double kappa=sqr(dipol()->lastPt()/hardpT);
double s = sqrt(1.-kappa);
pair<double,double> zbounds= make_pair(0.5*(1.+x-(1.-x)*s),0.5*(1.+x+(1.-x)*s));
return (zbounds.first<z_&&z_<zbounds.second);
}
Ptr<Node>::ptr Node::getHistory(bool normal,double hardScaleFactor) {
Ptr<Node>::ptr res=this;
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 (vector<Ptr<Node>::ptr>::iterator it=temp.begin();it!=temp.end();it++){
if( (*it)->dipol()->underlyingBornME()->largeNColourCorrelatedME2(
make_pair((*it)->dipol()->bornEmitter(),(*it)->dipol()->bornSpectator()),
deepHead()->MH()->largeNBasis())!=0.
){
if((*it)->nodeME()->dSigHatDR()/nanobarn!=0.){
subprosel.insert((abs((*it)->dipol()->dSigHatDR() /
(*it)->nodeME()->dSigHatDR()*deepHead()->MH()->as((*it)->dipol()->lastPt()))), (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
//TODO choosehistories
}
}
if (subprosel.empty())
return res;
res = subprosel.select(UseRandom::rnd());
temp = res->getNextOrderedNodes(true,hardScaleFactor);
}
return res;
}
bool Node::headContribution(double hardScaleFactor){
bool allabove=true;
vector<Ptr<Node>::ptr> temp2 = children();
for (vector<Ptr<Node>::ptr>::iterator it = temp2.begin(); it != temp2.end(); it++) {
allabove&=(*it)->dipol()->lastPt()>deepHead()->MH()->mergePt();
}
if(allabove){
Ptr<Node>::ptr tmpBorn = getHistory(true,hardScaleFactor);
if(!tmpBorn->parent())return false;
}
return true;
}
bool Node::DipolesAboveMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
double Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()) {
assert((*it)->dipol()->clustersafe());
minpt=min(minpt,(*it)->dipol()->lastPt());
assert((*it)->dipol()->clustersafe());
sum+=Di;
first_subpro.insert(1., (*it));
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
-double Node::calcPsMinusDip(Energy scale){
- return -1.* dipol()->dipMinusPs(sqr(scale),deepHead()->MH()->largeNBasis())/nanobarn;
+pair<double,double> Node::calcDipandPS(Energy scale){
+ return dipol()->dipandPs(sqr(scale),deepHead()->MH()->largeNBasis());
}
double Node::calcPs(Energy scale){
return dipol()->ps(sqr(scale),deepHead()->MH()->largeNBasis())/nanobarn;
}
-
+/*
bool Node::diffPsDipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
bool isInSomePS=false;
Energy maxx=0.*GeV;
Energy minn=100000*GeV;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
isInSomePS|=(*it)->inShowerPS((*it2)->dipol()->lastPt());
maxx=max(maxx,(*it)->dipol()->lastPt());
minn=min(minn,(*it)->dipol()->lastPt());
}
double Di=0.;
if(isInSomePS||(tmp2.empty())){
Di=-1.* (*it)->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->MH()->largeNBasis())/nanobarn;
}else{
Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
}
if ((Di!=0.)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<MH()->mergePt())) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++) assert(((*it2)->dipol()->lastPt()>deepHead()->MH()->mergePt()));
if (Di!=0.) {
first_subpro.insert(1., (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
vector<Ptr<Node>::ptr> tmp2=selectedNode->children();
bool isInSomePS=false;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
isInSomePS|=selectedNode->inShowerPS((*it2)->dipol()->lastPt());
}
if((isInSomePS||(tmp2.empty()))){//selectedNode->dipol()->lastPt()<MH()->mergePt()&&
sum=-1.* selectedNode->dipol()->dipMinusPs(sqr(10.*GeV),deepHead()->MH()->largeNBasis())/nanobarn;
}else{
sum=-1.* selectedNode->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
}
return true;
}
return false;
}
+*/
-
-
+/*
bool Node::psBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
double Di=-1.* (*it)->dipol()->ps(sqr(10.*GeV),deepHead()->MH()->largeNBasis())/nanobarn;
if ((Di!=0)&&(*it)->xcomb()->willPassCuts()){//&&((*it)->dipol()->lastPt()<MH()->mergePt())) {
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
bool isInSomePS=false;
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
assert(((*it2)->dipol()->lastPt()>deepHead()->MH()->mergePt())||((*it)->dipol()->lastPt()>deepHead()->MH()->mergePt()));
isInSomePS|=(*it)->inShowerPS((*it2)->dipol()->lastPt());
}
if(!isInSomePS&&!(tmp2.empty()))continue;
sum+=Di;
if (Di!=0) {
first_subpro.insert(1., (*it));
minpt=min(minpt,(*it)->dipol()->lastPt());
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
+*/
+/*
bool Node::dipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number){
sum=0.;
Selector<Ptr<Node>::ptr> first_subpro;
vector<Ptr<Node>::ptr> tmp=children();
for (vector<Ptr<Node>::ptr>::iterator it = tmp.begin(); it != tmp.end(); it++) {
if (true
//(*it)->dipol()->clustersafe() &&
// (*it)->xcomb()->willPassCuts()
) {
bool calcdip=true;
vector<Ptr<Node>::ptr> tmp2=(*it)->children();
for (vector<Ptr<Node>::ptr>::iterator it2 = tmp2.begin(); it2 != tmp2.end(); it2++){
assert(((*it2)->dipol()->lastPt()>deepHead()->MH()->mergePt())||((*it)->dipol()->lastPt()>deepHead()->MH()->mergePt()));
if(((*it2)->dipol()->lastPt()<deepHead()->MH()->mergePt())){
if((*it)->dipol()->lastPt()<deepHead()->MH()->mergePt()){
sum=0;
return false;
}
calcdip=false;
}
}
if(calcdip){
double Di=-1.* (*it)->dipol()->dSigHatDR(sqr(10.*GeV))/nanobarn;
sum+=Di;
minpt=min(minpt,(*it)->dipol()->lastPt());
if (Di!=0) {
first_subpro.insert(1., (*it));
}
}
}
}
number=int(first_subpro.size());
if (number!=0) {
selectedNode=first_subpro.select(UseRandom::rnd());
return true;
}
return false;
}
+*/
IBPtr Node::clone() const {
return new_ptr(*this);
}
IBPtr Node::fullclone() const {
return new_ptr(*this);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void Node::persistentOutput(PersistentOStream & os) const {
os <<
theheadxcomb << thexcomb << thenodeMEPtr <<
thedipol << thechildren <<theparent <<
theDeepHead << theCutStage <<
theDeepProStage << clustersafer << ounit(theVetoPt, GeV) <<
ounit(theRunningPt, GeV) << theSubtractedReal << theVirtualContribution<<theMergingHelper;
// *** ATTENTION *** os << ; // Add all member variable which should be written persistently here.
}
void Node::persistentInput(PersistentIStream & is, int) {
is >>
theheadxcomb >> thexcomb >> thenodeMEPtr >>
thedipol >> thechildren >> theparent >> theDeepHead >> theCutStage >>
theDeepProStage >> clustersafer >> iunit(theVetoPt, GeV) >>
iunit(theRunningPt, GeV) >> theSubtractedReal >> theVirtualContribution
>>theMergingHelper;
// *** ATTENTION *** is >> ; // Add all member variable which should be read persistently here.
}
// *** 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/DipoleShower/Merging/Node.h b/DipoleShower/Merging/Node.h
--- a/DipoleShower/Merging/Node.h
+++ b/DipoleShower/Merging/Node.h
@@ -1,434 +1,436 @@
// -*- C++ -*-
#ifndef Herwig_Node_H
#define Herwig_Node_H
//
// This is the declaration of the Node class.
//
//#include "Node.fh"
#include "MFactory.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/DipoleShower/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();
Node(Ptr<MatchboxMEBase>::ptr nodeME, int deepprostage,int cutstage);
Node(Ptr<Node>::ptr deephead,
Ptr<Node>::ptr head,
Ptr<SubtractionDipole>::ptr dipol,
Ptr<MatchboxMEBase>::ptr nodeME,
int deepprostage, 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, int proStage);
bool headContribution(double hardscalefactor);
bool DipolesAboveMergeingScale(Ptr<Node>::ptr& selectedNode,double& sum,Energy& minpt,int& number);
-
+
+ /*
bool diffPsDipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number);
bool psBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number);
bool dipBelowMergeingScale(Ptr<Node>::ptr& selectedNode,double & sum,Energy& minpt,int& number);
+ */
-
- double calcPsMinusDip(Energy scale);
+ pair<double,double> calcDipandPS(Energy scale);
double calcPs(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);
void firstgenerateKinematics(const double *r, int stage,Energy2 shat);
/** returns the matrix element pointer */
const Ptr<MatchboxMEBase>::ptr nodeME()const;
/** access the matrix element pointer */
Ptr<MatchboxMEBase>::ptr nodeME();
/** the parent node */
Ptr<Node>::ptr parent()const {
return theparent;
}
/** map of children nodes*/
vector< Ptr<Node>::ptr > children() {
return thechildren;
}
Ptr<Node>::ptr randomChild();
-
+
+ Energy miniPt()const ;
bool allAbove(Energy pt);
bool isInHistoryOf(Ptr<Node>::ptr other);
/** set the first node (godfather). 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);
//theProjectors.push_back(make_pair(a,pro));
}
/** insert nodes to projector vector */
vector< pair <double , Ptr<Node>::ptr > > Projector() {
return theProjectors;
}
/** returns the dipol of the node. */
Ptr<SubtractionDipole>::ptr dipol() const;
/** set the xcomb of the node */
void xcomb(StdXCombPtr xc) {
thexcomb = xc;
}
/** return the xcomb */
StdXCombPtr xcomb() const {
return thexcomb;
}
/** set the headxcomb of the node */
void headxcomb(StdXCombPtr xc) {
thexcomb = xc;
}
/** return the xcomb */
StdXCombPtr headxcomb() const {
return theheadxcomb;
}
Energy vetoPt() const {
return theVetoPt;
}
void vetoPt(Energy x) {
theVetoPt=x;
}
Energy runningPt(){return theRunningPt;}
void runningPt(Energy x){theRunningPt=x;}
int cutStage() const {
return theCutStage;
}
unsigned int DeepProStage() const {
return theDeepProStage;
}
SafeClusterMap clusterSafe() const {
return clustersafer;
}
bool ordered() {
return isOrdered;
}
int orderedSteps() const {
return theOrderdSteps;
}
void setOrderedSteps(int x) {
theOrderdSteps = x;
}
bool isThisSafe() const {
return isthissafe;
}
vector<Ptr<Node>::ptr> getNextOrderedNodes(bool normal=true,double hardscalefactor=1.);
bool inShowerPS(Energy hardpt);
vector<Ptr<Node>::ptr> getNextFullOrderedNodes();
bool hasFullHistory();
Ptr<Node>::ptr getHistory(bool normal=true,double hardscalefactor=1.);
bool subtractedReal() {
return theSubtractedReal;
}
void subtractedReal(bool x) {
theSubtractedReal = x;
}
bool virtualContribution() {
return theVirtualContribution ;
}
void virtualContribution(bool x) {
theVirtualContribution = x;
}
Ptr<Merger>::ptr MH()const{return theMergingHelper;}
void MH(Ptr<Merger>::ptr a){theMergingHelper=a;}
private:
StdXCombPtr theheadxcomb;
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;
/**
* all nodes downstream have pt over merged pt.
*/
bool isthissafe;
/**
* The DeepProStage is the number of additional partons of the DeepHead
* compared to the lowest multiplicity.
*/
int theDeepProStage;
/**
* For [[Emitter,Emission],Spectator] the mapped pair gives
* information if the first and the second cluster is safe.
*/
SafeClusterMap clustersafer;
bool isOrdered;
int theOrderdSteps;
bool theNeedFullOrderedHistory;
unsigned int theSudakovSteps;
Energy theVetoPt;
Energy theRunningPt;
bool theSubtractedReal;
bool theVirtualContribution;
Ptr<Merger>::ptr theMergingHelper;
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 */
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,1358 +1,1361 @@
// -*- C++ -*-
//
// SubtractionDipole.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SubtractionDipole class.
//
#include "SubtractionDipole.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "Herwig/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
SubtractionDipole::SubtractionDipole()
: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
theIgnoreCuts(false),
theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
lastRealEmissionKey(realEmissionKey(cPDVector(),-1,-1,-1)),
lastUnderlyingBornKey(underlyingBornKey(cPDVector(),-1,-1)),
theBornEmitter(-1), theBornSpectator(-1),
theLastSubtractionScale(ZERO), theLastSplittingScale(ZERO),
theLastSubtractionPt(ZERO), theLastSplittingPt(ZERO),
theLastSubtractionZ(0.0), theLastSplittingZ(0.0),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
theLoopSimSubtraction(false), theRealEmissionScales(false),
theShowerHardScale(ZERO), theShowerScale(ZERO),
theIsInShowerPhasespace(false), theIsAboveCutoff(false) {}
SubtractionDipole::~SubtractionDipole() {}
double SubtractionDipole::alpha() const{
return factory()->alphaParameter();
}
void SubtractionDipole::clearBookkeeping() {
theRealEmitter = -1;
theRealEmission = -1;
theRealSpectator = -1;
theBornEmitter = -1;
theBornSpectator = -1;
theMergingMap.clear();
theSplittingMap.clear();
theIndexMap.clear();
theUnderlyingBornDiagrams.clear();
theRealEmissionDiagrams.clear();
theBornToRealDiagrams.clear();
theRealToBornDiagrams.clear();
}
void SubtractionDipole::setupBookkeeping(const map<Ptr<DiagramBase>::ptr,SubtractionDipole::MergeInfo>& mergeInfo,bool slim) {
theMergingMap.clear();
theSplittingMap.clear();
theUnderlyingBornDiagrams.clear();
theRealEmissionDiagrams.clear();
theBornToRealDiagrams.clear();
theRealToBornDiagrams.clear();
int xemitter = -1;
int xspectator = -1;
map<int,int> mergeLegs;
map<int,int> remapLegs;
map<int,int> realBornMap;
map<int,int> bornRealMap;
set<Ptr<DiagramBase>::cptr> usedDiagrams;
for ( map<Ptr<DiagramBase>::ptr,MergeInfo>::const_iterator mit = mergeInfo.begin();
mit != mergeInfo.end(); ++mit ) {
DiagramVector::const_iterator bd =
theUnderlyingBornME->diagrams().end();
// work out the most similar underlying Born diagram
map<int,int> xRemapLegs;
int nomapScore = 0;
for ( DiagramVector::const_iterator b =
theUnderlyingBornME->diagrams().begin();
b != theUnderlyingBornME->diagrams().end(); ++b ) {
map<int,int> theRemapLegs;
if ( mit->second.diagram->isSame(*b,theRemapLegs) &&
usedDiagrams.find(*b) == usedDiagrams.end() ) {
int theNomapScore = 0;
for ( map<int,int>::const_iterator m = theRemapLegs.begin();
m != theRemapLegs.end(); ++m )
if ( m->first == m->second )
theNomapScore += 1;
if ( theNomapScore >= nomapScore ) {
nomapScore = theNomapScore;
xRemapLegs = theRemapLegs;
bd = b;
}
}
}
// no underlying Born
if ( bd == theUnderlyingBornME->diagrams().end() )
continue;
// as we deal with one splitting only we now mark this diagram as used
// since we fixed the overall remapping of the process from the first
// occurence, see below. TODO: This confuses this code even more, and
// clearly calls for a cleanup. This is just grown historically and got
// messed up with experiencing different processes and setups.
usedDiagrams.insert(*bd);
if ( xemitter == -1 ) {
xemitter = mit->second.emitter;
mergeLegs = mit->second.mergeLegs;
remapLegs = xRemapLegs;
assert(remapLegs.find(xemitter) != remapLegs.end());
xemitter = remapLegs[xemitter];
// work out the leg remapping real -> born
for ( map<int,int>::const_iterator k = mergeLegs.begin();
k != mergeLegs.end(); ++k ) {
assert(remapLegs.find(k->second) != remapLegs.end());
realBornMap[k->first] = remapLegs[k->second];
}
// work out the leg remapping born -> real
for ( map<int,int>::const_iterator k = realBornMap.begin();
k != realBornMap.end(); ++k ) {
bornRealMap[k->second] = k->first;
}
// work out the spectator
assert(mergeLegs.find(realSpectator()) != mergeLegs.end());
assert(remapLegs.find(mergeLegs[realSpectator()]) != remapLegs.end());
xspectator = realBornMap[realSpectator()];
}
RealEmissionKey realKey = realEmissionKey((*mit->first).partons(),realEmitter(),realEmission(),realSpectator());
UnderlyingBornKey bornKey = underlyingBornKey((**bd).partons(),xemitter,xspectator);
if ( theMergingMap.find(realKey) == theMergingMap.end() )
theMergingMap.insert(make_pair(realKey,make_pair(bornKey,realBornMap)));
RealEmissionInfo realInfo = make_pair(realKey,bornRealMap);
bool gotit = false;
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spIterator;
pair<spIterator,spIterator> range = theSplittingMap.equal_range(bornKey);
for ( ; range.first != range.second; ++range.first )
if ( range.first->second == realInfo ) {
gotit = true;
break;
}
if ( !gotit )
theSplittingMap.insert(make_pair(bornKey,realInfo));
theUnderlyingBornDiagrams[process(realKey)].push_back(*bd);
theRealEmissionDiagrams[process(bornKey)].push_back(mit->first);
theBornToRealDiagrams[*bd] = mit->first;
theRealToBornDiagrams[mit->first] = *bd;
}
if (slim) {
theIndexMap.clear();
theSplittingMap.clear();
theBornToRealDiagrams.clear();
theRealEmissionDiagrams.clear();
}
if ( theSplittingMap.empty() )
return;
theIndexMap.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator s =
theSplittingMap.begin(); s != theSplittingMap.end(); ++s ) {
theIndexMap[process(s->first)] = make_pair(emitter(s->first),spectator(s->first));
}
}
double SubtractionDipole::jacobianMerging(Energy2 shatBorn,Energy2 shatReal,int n) const {return theTildeKinematics->jacobian(shatBorn,shatReal,n);}
void SubtractionDipole::subtractionBookkeeping() {
/*
if ( theMergingMap.empty() )
setupBookkeeping();
*/
assert(!theMergingMap.empty());
lastRealEmissionKey =
realEmissionKey(lastHeadXComb().mePartonData(),realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() ) {
theApply = false;
return;
}
theApply = true;
lastUnderlyingBornKey = k->second.first;
bornEmitter(emitter(lastUnderlyingBornKey));
bornSpectator(spectator(lastUnderlyingBornKey));
}
void SubtractionDipole::splittingBookkeeping() {
/*
if ( theMergingMap.empty() )
setupBookkeeping();
*/
assert(!theMergingMap.empty());
map<cPDVector,pair<int,int> >::const_iterator esit =
theIndexMap.find(lastHeadXComb().mePartonData());
if ( esit == theIndexMap.end() ) {
theApply = false;
return;
}
theApply = true;
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(lastHeadXComb().mePartonData(),bornEmitter(),bornSpectator());
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
assert(kr.first != kr.second);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
assert(lastRealEmissionInfo != kr.second);
lastRealEmissionKey = lastRealEmissionInfo->second.first;
realEmitter(emitter(lastRealEmissionKey));
realEmission(emission(lastRealEmissionKey));
realSpectator(spectator(lastRealEmissionKey));
}
StdXCombPtr SubtractionDipole::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allBins,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
return
realEmissionME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
StdXCombPtr SubtractionDipole::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
return
realEmissionME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
StdXCombPtr SubtractionDipole::makeBornXComb(tStdXCombPtr realXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();
lastRealEmissionKey =
realEmissionKey(proc,realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() )
return StdXCombPtr();
PartonPairVec pbs = realXC->pExtractor()->getPartons(realXC->maxEnergy(),
realXC->particles(),
*(realXC->cuts()));
DiagramVector bornDiags = underlyingBornDiagrams(proc);
assert(!bornDiags.empty());
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == bornDiags.front()->partons()[0] &&
ppit->second->parton() == bornDiags.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
return
underlyingBornME()->makeXComb(realXC,*ppit,bornDiags,this);
}
vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();
map<cPDVector,pair<int,int> >::const_iterator esit = theIndexMap.find(proc);
if ( esit == theIndexMap.end() )
return vector<StdXCombPtr>();
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(proc,bornEmitter(),bornSpectator());
if ( theSplittingMap.find(lastUnderlyingBornKey) == theSplittingMap.end() )
return vector<StdXCombPtr>();
PartonPairVec pbs = bornXC->pExtractor()->getPartons(bornXC->maxEnergy(),
bornXC->particles(),
*(bornXC->cuts()));
DiagramVector realDiags = realEmissionDiagrams(proc);
assert(!realDiags.empty());
vector<StdXCombPtr> res;
map<cPDVector,DiagramVector> realProcs;
for ( MEBase::DiagramVector::const_iterator d = realDiags.begin();
d != realDiags.end(); ++d ) {
realProcs[(**d).partons()].push_back(*d);
}
for ( map<cPDVector,DiagramVector>::const_iterator pr =
realProcs.begin(); pr != realProcs.end(); ++pr ) {
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == pr->second.front()->partons()[0] &&
ppit->second->parton() == pr->second.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
StdXCombPtr rxc =
realEmissionME()->makeXComb(bornXC,*ppit,pr->second,this);
res.push_back(rxc);
}
return res;
}
const MEBase::DiagramVector& SubtractionDipole::underlyingBornDiagrams(const cPDVector& real) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theUnderlyingBornDiagrams.find(real);
if (k == theUnderlyingBornDiagrams.end() )
return empty;
return k->second;
}
tcDiagPtr SubtractionDipole::underlyingBornDiagram(tcDiagPtr realDiag) const {
map<tcDiagPtr,tcDiagPtr>::const_iterator it = theRealToBornDiagrams.find(realDiag);
assert(it != theRealToBornDiagrams.end());
return it->second;
}
const MEBase::DiagramVector& SubtractionDipole::realEmissionDiagrams(const cPDVector& born) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theRealEmissionDiagrams.find(born);
if ( k == theRealEmissionDiagrams.end() )
return empty;
return k->second;
}
tcDiagPtr SubtractionDipole::realEmissionDiagram(tcDiagPtr bornDiag) const {
map<tcDiagPtr,tcDiagPtr>::const_iterator it = theBornToRealDiagrams.find(bornDiag);
assert(it != theBornToRealDiagrams.end());
return it->second;
}
void SubtractionDipole::getDiagrams() const {
if ( splitting() ) {
realEmissionME()->diagrams();
useDiagrams(realEmissionME());
} else {
underlyingBornME()->diagrams();
useDiagrams(underlyingBornME());
}
}
Selector<MEBase::DiagramIndex> SubtractionDipole::diagrams(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->setXComb(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagrams(dv);
}
Selector<const ColourLines *>
SubtractionDipole::colourGeometries(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->colourGeometries(diag) :
underlyingBornME()->colourGeometries(diag);
}
const ColourLines &
SubtractionDipole::selectColourGeometry(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->selectColourGeometry(diag) :
underlyingBornME()->selectColourGeometry(diag);
}
void SubtractionDipole::flushCaches() {
theUnderlyingBornME->flushCaches();
theRealEmissionME->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
}
void SubtractionDipole::setXComb(tStdXCombPtr xc) {
if ( !xc ) {
theApply = false;
return;
} else {
theApply = true;
}
lastMatchboxXComb(xc);
MEBase::setXComb(xc);
if ( splitting() ) {
realEmissionME()->setXComb(xc);
underlyingBornME()->setXComb(xc->head());
splittingBookkeeping();
} else {
realEmissionME()->setXComb(xc->head());
underlyingBornME()->setXComb(xc);
subtractionBookkeeping();
}
if ( !apply() )
return;
}
void SubtractionDipole::setKinematics() {
MEBase::setKinematics();
if ( splitting() )
realEmissionME()->setKinematics();
else
underlyingBornME()->setKinematics();
}
bool SubtractionDipole::generateKinematics(const double * r) {
if ( lastXCombPtr()->kinematicsGenerated() )
return true;
if ( splitting() ) {
if ( !generateRadiationKinematics(r) )
return false;
if( ! realEmissionME()->lastXCombPtr()->setIncomingPartons())
return false;
realEmissionME()->setScale();
double jac = jacobian();
jac *= pow(underlyingBornME()->lastXComb().lastSHat() / realEmissionME()->lastXComb().lastSHat(),
realEmissionME()->lastXComb().mePartonData().size()-4.);
jacobian(jac);
assert(lastXCombPtr() == realEmissionME()->lastXCombPtr());
lastXCombPtr()->didGenerateKinematics();
return true;
}
if ( !generateTildeKinematics() )
{cout<<"\nxx 1 xx"; return false;}
if( ! underlyingBornME()->lastXCombPtr()->setIncomingPartons() )
{cout<<"\nxx 2 xx"; return false;}
underlyingBornME()->setScale();
assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
if( ! underlyingBornME()->lastXCombPtr()->setIncomingPartons() )
{cout<<"\nxx 3 xx"; return false;}
// need to have the scale and x's available for checking shower phase space
if ( showerApproximation() &&
lastXCombPtr()->willPassCuts() )
showerApproximation()->getShowerVariables();
lastXCombPtr()->didGenerateKinematics();
return true;
}
int SubtractionDipole::nDim() const {
if ( !splitting() )
return underlyingBornME()->nDim();
return underlyingBornME()->nDim() + nDimRadiation();
}
void SubtractionDipole::clearKinematics() {
MEBase::clearKinematics();
if ( splitting() )
realEmissionME()->clearKinematics();
else
underlyingBornME()->clearKinematics();
}
void SubtractionDipole::tildeKinematics(Ptr<TildeKinematics>::tptr tk) {
theTildeKinematics = tk;
}
bool SubtractionDipole::generateTildeKinematics() {
assert(!splitting());
Ptr<TildeKinematics>::tptr kinematics = theTildeKinematics;
if ( showerApproximation() ) {
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
showerApproximation()->checkCutoff();
if ( showerApproximation()->showerTildeKinematics() &&
isAboveCutoff() &&
realShowerSubtraction() )
kinematics = showerApproximation()->showerTildeKinematics();
}
if ( !kinematics ) {
jacobian(0.0);
return false;
}
kinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());
if ( !kinematics->doMap() ) {
jacobian(0.0);
return false;
}
theLastSubtractionScale = kinematics->lastScale();
theLastSubtractionPt = kinematics->lastPt();
theLastSubtractionZ = kinematics->lastZ();
theLastRealR=kinematics->lastRealR();
theLastBornR=kinematics->lastBornR();
meMomenta().resize(lastHeadXComb().meMomenta().size() - 1);
assert(mergingMap().find(lastRealEmissionKey) != mergingMap().end());
map<int,int>& trans = theMergingMap[lastRealEmissionKey].second;
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == realEmitter() || k == realEmission() || k == realSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( kinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[bornEmitter()] =
const_cast<const TildeKinematics&>(*kinematics).bornEmitterMomentum();
meMomenta()[bornSpectator()] =
const_cast<const TildeKinematics&>(*kinematics).bornSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).hardProcessMass());
p->rescaleRho();
}
jacobian(realEmissionME()->lastXComb().jacobian());
logGenerateTildeKinematics();
return true;
}
void SubtractionDipole::invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr itk) {
theInvertedTildeKinematics = itk;
}
int SubtractionDipole::nDimRadiation() const {
return invertedTildeKinematics() ?
invertedTildeKinematics()->nDimRadiation() :
0;
}
bool SubtractionDipole::generateRadiationKinematics(const double * r) {
assert(splitting());
Ptr<InvertedTildeKinematics>::tptr kinematics = theInvertedTildeKinematics;
if ( showerApproximation() ) {
showerApproximation()->setBornXComb(lastHeadXCombPtr());
showerApproximation()->setRealXComb(lastXCombPtr());
showerApproximation()->setDipole(this);
if ( showerApproximation()->showerInvertedTildeKinematics() ) {
kinematics = showerApproximation()->showerInvertedTildeKinematics();
}
}
if ( !kinematics ) {
jacobian(0.0);
return false;
}
kinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());
if ( !kinematics->doMap(r) ) {
jacobian(0.0);
return false;
}
theLastSplittingScale = kinematics->lastScale();
theLastSplittingPt = kinematics->lastPt();
theLastSplittingZ = kinematics->lastZ();
meMomenta().resize(lastHeadXComb().meMomenta().size() + 1);
assert(splittingMap().find(lastUnderlyingBornKey) != splittingMap().end());
map<int,int>& trans = const_cast<map<int,int>&>(lastRealEmissionInfo->second.second);
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == bornEmitter() || k == bornSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( kinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = kinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[realEmitter()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realEmitterMomentum();
meMomenta()[realEmission()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realEmissionMomentum();
meMomenta()[realSpectator()] =
const_cast<const InvertedTildeKinematics&>(*kinematics).realSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).hardProcessMass());
p->rescaleRho();
}
jacobian(underlyingBornME()->lastXComb().jacobian() *
kinematics->jacobian());
logGenerateRadiationKinematics(r);
return true;
}
void SubtractionDipole::ptCut(Energy cut) {
theInvertedTildeKinematics->ptCut(cut);
}
CrossSection SubtractionDipole::dSigHatDR(Energy2 factorizationScale) const {
double pdfweight = 1.;
double jac = jacobian();
if ( splitting() && jac == 0.0 ) {
lastMECrossSection(ZERO);
return ZERO;
}
if ( factorizationScale == ZERO ) {
factorizationScale = underlyingBornME()->lastScale();
}
if ( havePDFWeight1() ) {
pdfweight *= realEmissionME()->pdf1(factorizationScale);
}
if ( havePDFWeight2() ) {
pdfweight *= realEmissionME()->pdf2(factorizationScale);
}
lastMEPDFWeight(pdfweight);
bool needTheDipole = true;
CrossSection shower = ZERO;
double lastThetaMu = 1.0;
double showerFactor = 1.;
if ( showerApproximation() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(const_cast<SubtractionDipole*>(this));
if ( !isAboveCutoff() ) {
showerApproximation()->wasBelowCutoff();
lastThetaMu = 0.0;
} else {
lastThetaMu = 1.0;
}
if ( lastThetaMu > 0.0 && isInShowerPhasespace() ) {
if ( realShowerSubtraction() )
shower = showerApproximation()->dSigHatDR()*lastThetaMu;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
shower = -showerApproximation()->dSigHatDR()*lastThetaMu;
if ( virtualShowerSubtraction() &&
isAboveCutoff() &&
showerApproximation()->showerTildeKinematics() ) {
// map shower to dipole kinematics; we are always above the
// cutoff in this case
showerFactor *=
showerApproximation()->showerTildeKinematics()->jacobianRatio();
}
shower *= showerFactor;
}
if ( realShowerSubtraction() && lastThetaMu == 1.0 )
needTheDipole = false;
if ( virtualShowerSubtraction() && lastThetaMu == 0.0 )
needTheDipole = false;
if ( factory()->loopSimCorrections() ||
factory()->meCorrectionsOnly() )
needTheDipole = false;
}
double xme2 = 0.0;
if ( needTheDipole )
xme2 = me2();
if ( factory()->loopSimCorrections() ||
factory()->meCorrectionsOnly() ) {
assert(showerApproximation());
xme2 = realEmissionME()->me2() * showerApproximation()->channelWeight();
double rws =
pow(underlyingBornME()->lastXComb().lastAlphaS()/
realEmissionME()->lastXComb().lastAlphaS(),
realEmissionME()->orderInAlphaS());
xme2 *= rws;
double rwe =
pow(underlyingBornME()->lastXComb().lastAlphaEM()/
realEmissionME()->lastXComb().lastAlphaEM(),
underlyingBornME()->orderInAlphaEW());
xme2 *= rwe;
}
if ( realShowerSubtraction() )
xme2 *= 1. - lastThetaMu;
if ( virtualShowerSubtraction() || loopSimSubtraction() )
xme2 *= lastThetaMu;
double coupl = lastMECouplings();
coupl *= underlyingBornME()->lastXComb().lastAlphaS();
lastMECouplings(coupl);
CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
if ( !showerApproximation() && xme2 != 0.0 ) {
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(theRealEmissionME->lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
}
lastMECrossSection(-res-shower);
logDSigHatDR(jac);
return lastMECrossSection();
}
CrossSection SubtractionDipole::ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
double jac = jacobian();
assert( factorizationScale != ZERO );
assert (! splitting());
if(!theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second || jac == 0.0 ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
double pdfweight = 1.;
if ( havePDFWeight1() ) pdfweight *= realEmissionME()->pdf1(factorizationScale);
if ( havePDFWeight2() ) pdfweight *= realEmissionME()->pdf2(factorizationScale);
double ccme2 =underlyingBornME()->me2()*
underlyingBornME()->largeNColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),largeNBasis)
/underlyingBornME()->largeNME2(largeNBasis);
//double ccme2 =underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
return -sqr(hbarc) * jac * pdfweight * (me2Avg(ccme2)) / (2. * realEmissionME()->lastXComb().lastSHat());
}
-CrossSection SubtractionDipole::dipMinusPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
+pair<double,double> SubtractionDipole::dipandPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const {
double jac = jacobian();
assert( factorizationScale != ZERO );
assert (! splitting());
if(!theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second || jac == 0.0 ) {
lastMECrossSection(ZERO);
lastME2(0.0);
- return ZERO;
+ return make_pair(0.,0.);
}
double pdfweight = 1.;
if ( havePDFWeight1() ) pdfweight *= realEmissionME()->pdf1(factorizationScale);
if ( havePDFWeight2() ) pdfweight *= realEmissionME()->pdf2(factorizationScale);
double ccme2 =underlyingBornME()->me2()
*
underlyingBornME()->largeNColourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()),largeNBasis)
/underlyingBornME()->largeNME2(largeNBasis);
// cout<<"\nln "<<ccme2;
//ccme2 =underlyingBornME()->colourCorrelatedME2(make_pair(bornEmitter(),bornSpectator()));
// cout<<"\nan "<<ccme2;
double ps = me2Avg(ccme2);
double dip = me2();
+ double factor= sqr(hbarc) * jac * pdfweight / (2. * realEmissionME()->lastXComb().lastSHat())/nanobarn;
- return -sqr(hbarc) * jac * pdfweight * (dip-ps) / (2. * realEmissionME()->lastXComb().lastSHat());
+ pair<double,double> res=make_pair(factor*dip,factor*ps);
+
+ return res ;
}
bool SubtractionDipole::clustersafe()const {
return (theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second);
}
bool SubtractionDipole::clustersafe(){
return (theRealEmissionME->clustersafe(realEmitter(),realEmission(),realSpectator()).second);
}
void SubtractionDipole::print(ostream& os) const {
os << "--- SubtractionDipole setup ----------------------------------------------------\n";
os << " subtraction '" << name() << "'\n for real emission '"
<< theRealEmissionME->name() << "'\n using underlying Born '"
<< theUnderlyingBornME->name() << "'\n";
os << " tilde kinematics are '"
<< (theTildeKinematics ? theTildeKinematics->name() : "")
<< " '\n inverted tilde kinematics are '"
<< (theInvertedTildeKinematics ? theInvertedTildeKinematics->name() : "") << "'\n";
os << " the following subtraction mappings have been found:\n";
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator m =
theMergingMap.begin(); m != theMergingMap.end(); ++m ) {
os << " " << process(m->second.first)[0]->PDGName() << " "
<< process(m->second.first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->second.first).begin() + 2;
p != process(m->second.first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[" << emitter(m->second.first) << "," << spectator(m->second.first) << "] <=> ";
os << process(m->first)[0]->PDGName() << " "
<< process(m->first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->first).begin() + 2;
p != process(m->first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[(" << emitter(m->first) << "," << emission(m->first) << ")," << spectator(m->first) << "]\n"
<< " non-dipole momenta ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->second << " ";
}
os << ") <=> ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->first << " ";
}
os << ")\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractionDipole::printLastEvent(ostream& os) const {
os << "--- SubtractionDipole last event information -----------------------------------\n";
os << " for dipole '" << name() << "' applying ["
<< bornEmitter() << "," << bornSpectator() << "] <=> [("
<< realEmitter() << "," << realEmission() << ")," << realSpectator() << "]\n"
<< " evaluated the cross section/nb " << (lastMECrossSection()/nanobarn) << "\n"
<< " with subtraction parameters x[0] = " << subtractionParameters()[0]
<< " x[1] = " << subtractionParameters()[1] << "\n";
os << " the last real emission event was:\n";
realEmissionME()->printLastEvent(os);
os << " the last underlying Born event was:\n";
underlyingBornME()->printLastEvent(os);
os << "--- end SubtractionDipole last event information -------------------------------\n";
os << flush;
}
void SubtractionDipole::logME2() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated me2 using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "Born phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = bornxc->meMomenta().begin();
cPDVector::const_iterator dit = bornxc->mePartonData().begin();
for ( ; pit != bornxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << bornxc->lastX1() << " x2 = " << bornxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (bornxc->lastSHat()/GeV2) << "\n";
generator()->log() << "Real emission phase space point (in GeV):\n";
pit = realxc->meMomenta().begin();
dit = realxc->mePartonData().begin();
for ( ; pit != realxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << realxc->lastX1() << " x2 = " << realxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (realxc->lastSHat()/GeV2) << "\n";
}
void SubtractionDipole::logDSigHatDR(double effectiveJac) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated cross section using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n"
<< "Jacobian = " << jacobian()
<< " effective Jacobian = " << effectiveJac << "\n"
<< "Born sHat/GeV2 = " << (bornxc->lastSHat()/GeV2)
<< " real sHat/GeV2 = " << (realxc->lastSHat()/GeV2)
<< " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void SubtractionDipole::logGenerateTildeKinematics() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating tilde kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with real xcomb " << lastHeadXCombPtr() << " born xcomb "
<< lastXCombPtr() << "\n"
<< "from real emission phase space point:\n";
Lorentz5Momentum rSum;
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
size_t count = 0;
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
rSum -= *pr;
} else {
rSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (rSum/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n"
<< "with scale/GeV = " << (theLastSubtractionScale/GeV)
<< "and pt/GeV = " << (theLastSubtractionPt/GeV) << "\n";
generator()->log() << "generated tilde kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
count = 0;
Lorentz5Momentum bSum;
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
bSum -= *pr;
} else {
bSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (bSum/GeV) << "\n";
generator()->log() << "Jacobian = " << jacobian() << "\n" << flush;
}
void SubtractionDipole::logGenerateRadiationKinematics(const double * r) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating radiation kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with born xcomb " << lastHeadXCombPtr() << " real xcomb "
<< lastXCombPtr() << "\n"
<< "from random numbers:\n";
copy(r,r+nDimRadiation(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "and born phase space point:\n";
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n" << flush;
generator()->log() << "scales: scale/GeV = " << (theLastSplittingScale/GeV)
<< " pt/GeV = " << (theLastSplittingPt/GeV) << "\n" << flush;
generator()->log() << "generated real emission kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "Jacobian = "
<< jacobian() << " = "
<< underlyingBornME()->lastXComb().jacobian()
<< "|Born * "
<< invertedTildeKinematics()->jacobian()
<< "|Radiation\n" << flush;
}
void SubtractionDipole::doinit() {
MEBase::doinit();
if ( underlyingBornME() ) {
theUnderlyingBornME->init();
}
if ( realEmissionME() ) {
theRealEmissionME->init();
}
if ( tildeKinematics() ) {
theTildeKinematics->init();
}
if ( invertedTildeKinematics() ) {
theInvertedTildeKinematics->init();
}
if ( showerApproximation() ) {
theShowerApproximation->init();
}
for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
p != thePartners.end(); ++p ) {
(**p).init();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).init();
}
}
void SubtractionDipole::doinitrun() {
MEBase::doinitrun();
if ( underlyingBornME() ) {
theUnderlyingBornME->initrun();
}
if ( realEmissionME() ) {
theRealEmissionME->initrun();
}
if ( tildeKinematics() ) {
theTildeKinematics->initrun();
}
if ( invertedTildeKinematics() ) {
theInvertedTildeKinematics->initrun();
}
if ( showerApproximation() ) {
theShowerApproximation->initrun();
}
for ( vector<Ptr<SubtractionDipole>::tptr>::iterator p = thePartners.begin();
p != thePartners.end(); ++p ) {
(**p).initrun();
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).initrun();
}
}
void SubtractionDipole::cloneDependencies(const std::string& prefix,bool slim) {
if ( underlyingBornME() ) {
Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
myUnderlyingBornME->cloneDependencies(pname.str(),slim);
underlyingBornME(myUnderlyingBornME);
}
if ( realEmissionME()&& !slim ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = realEmissionME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Matrix element " << pname.str() << " already existing." << Exception::runerror;
myRealEmissionME->cloneDependencies(pname.str());
realEmissionME(myRealEmissionME);
}
if ( tildeKinematics() ) {
Ptr<TildeKinematics>::ptr myTildeKinematics = tildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myTildeKinematics->name();
if ( ! (generator()->preinitRegister(myTildeKinematics,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
myTildeKinematics->dipole(this);
tildeKinematics(myTildeKinematics);
}
if ( invertedTildeKinematics()&& !slim ) {
Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Inverted tilde kinematics " << pname.str() << " already existing." << Exception::runerror;
myInvertedTildeKinematics->dipole(this);
invertedTildeKinematics(myInvertedTildeKinematics);
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw Exception() << "SubtractionDipole::cloneDependencies(): Reweight " << pname.str() << " already existing." << Exception::runerror;
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
}
void SubtractionDipole::constructVertex(tSubProPtr sub) {
if ( splitting() )
realEmissionME()->constructVertex(sub);
else
underlyingBornME()->constructVertex(sub);
}
void SubtractionDipole::constructVertex(tSubProPtr sub, const ColourLines* cl) {
if ( splitting() )
realEmissionME()->constructVertex(sub,cl);
else
underlyingBornME()->constructVertex(sub,cl);
}
void SubtractionDipole::generateSubCollision(SubProcess & sub) {
if ( splitting() )
realEmissionME()->generateSubCollision(sub);
else
underlyingBornME()->generateSubCollision(sub);
}
void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theSplitting << theApply << theSubtractionTest
<< theIgnoreCuts << theRealEmissionME << theUnderlyingBornME
<< thePartners << theTildeKinematics << theInvertedTildeKinematics
<< theReweights << theRealEmitter << theRealEmission << theRealSpectator
<< theSubtractionParameters << theMergingMap << theSplittingMap
<< theIndexMap << theUnderlyingBornDiagrams << theRealEmissionDiagrams
<< theBornToRealDiagrams << theRealToBornDiagrams
<< lastRealEmissionKey << lastUnderlyingBornKey
<< theBornEmitter << theBornSpectator << ounit(theLastSubtractionScale,GeV)
<< ounit(theLastSplittingScale,GeV) << ounit(theLastSubtractionPt,GeV)
<< ounit(theLastSplittingPt,GeV) << theLastSubtractionZ
<< theLastSplittingZ << theShowerApproximation
<< theRealShowerSubtraction << theVirtualShowerSubtraction
<< theLoopSimSubtraction << theRealEmissionScales << theFactory
<< ounit(theShowerHardScale,GeV) << ounit(theShowerScale,GeV)
<< theShowerParameters << theIsInShowerPhasespace << theIsAboveCutoff;
}
void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theSplitting >> theApply >> theSubtractionTest
>> theIgnoreCuts >> theRealEmissionME >> theUnderlyingBornME
>> thePartners >> theTildeKinematics >> theInvertedTildeKinematics
>> theReweights >> theRealEmitter >> theRealEmission >> theRealSpectator
>> theSubtractionParameters >> theMergingMap >> theSplittingMap
>> theIndexMap >> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
>> theBornToRealDiagrams >> theRealToBornDiagrams
>> lastRealEmissionKey >> lastUnderlyingBornKey
>> theBornEmitter >> theBornSpectator >> iunit(theLastSubtractionScale,GeV)
>> iunit(theLastSplittingScale,GeV) >> iunit(theLastSubtractionPt,GeV)
>> iunit(theLastSplittingPt,GeV) >> theLastSubtractionZ
>> theLastSplittingZ >> theShowerApproximation
>> theRealShowerSubtraction >> theVirtualShowerSubtraction
>> theLoopSimSubtraction >> theRealEmissionScales >> theFactory
>> iunit(theShowerHardScale,GeV) >> iunit(theShowerScale,GeV)
>> theShowerParameters >> theIsInShowerPhasespace >> theIsAboveCutoff;
lastMatchboxXComb(theLastXComb);
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
}
Ptr<MatchboxFactory>::tptr SubtractionDipole::factory() const {
return theFactory;
}
void SubtractionDipole::factory(Ptr<MatchboxFactory>::tptr f) {
theFactory = f;
}
void SubtractionDipole::Init() {
static ClassDocumentation<SubtractionDipole> documentation
("SubtractionDipole represents a dipole subtraction "
"term in the formalism of Catani and Seymour.");
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<SubtractionDipole,MEBase>
describeSubtractionDipole("Herwig::SubtractionDipole", "Herwig.so");
diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h
@@ -1,1265 +1,1265 @@
// -*- C++ -*-
//
// SubtractionDipole.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SubtractionDipole_H
#define HERWIG_SubtractionDipole_H
//
// This is the declaration of the SubtractionDipole class.
//
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig/MatrixElement/Matchbox/Phasespace/TildeKinematics.fh"
#include "Herwig/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.fh"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief SubtractionDipole represents a dipole subtraction
* term in the formalism of Catani and Seymour.
*
*/
class SubtractionDipole:
public MEBase, public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SubtractionDipole();
/**
* The destructor.
*/
virtual ~SubtractionDipole();
//@}
public:
/**
* Return the factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tptr factory() const;
/**
* Set the factory which produced this matrix element
*/
void factory(Ptr<MatchboxFactory>::tptr f);
/** @name Subprocess and diagram information. */
//@{
/**
* A helpre struct to communicate diagram merging and remapping
* information
*/
struct MergeInfo {
/**
* The merged emitter
*/
int emitter;
/**
* The Born diagram
*/
Ptr<Tree2toNDiagram>::ptr diagram;
/**
* The merging map
*/
map<int,int> mergeLegs;
};
/**
* Return true, if this dipole can possibly handle the indicated
* emitter.
*/
virtual bool canHandleEmitter(const cPDVector& partons, int emitter) const = 0;
/**
* Return true, if this dipole can possibly handle the indicated
* splitting.
*/
virtual bool canHandleSplitting(const cPDVector& partons, int emitter, int emission) const = 0;
/**
* Return true, if this dipole can possibly handle the indicated
* spectator.
*/
virtual bool canHandleSpectator(const cPDVector& partons, int spectator) const = 0;
/**
* Return true, if this dipole applies to the selected
* configuration.
*/
virtual bool canHandle(const cPDVector& partons,
int emitter, int emission, int spectator) const = 0;
/**
* Return true, if this dipole is symmetric with respect to emitter
* and emission.
*/
virtual bool isSymmetric() const { return false; }
/**
* If this is a dependent matrix element in a ME group, return true,
* if it applies to the process set in lastXComb()
*/
virtual bool apply() const { return theApply; }
/**
* Clear the bookkeeping
*/
void clearBookkeeping();
/**
* Setup bookkeeping maps.
*/
void setupBookkeeping(const map<Ptr<DiagramBase>::ptr,MergeInfo>& mergeInfo,bool slim);
/**
* Get bookkeeping information for the given
* real emission diagram
*/
void subtractionBookkeeping();
/**
* Determine bookkeeping information for
* the underlying Born process supplied through
* the lastHeadXComb() object.
*/
void splittingBookkeeping();
/**
* For the given event generation setup return a xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr 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 = tStdXCombPtr(),
tMEPtr newME = tMEPtr());
/**
* For the given event generation setup return a dependent xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME = tMEPtr());
/**
* Create a dependent xcomb object for the underlying
* Born process, given a XComb driving the real emission
*/
StdXCombPtr makeBornXComb(tStdXCombPtr realXC);
/**
* Create dependent xcomb objects for the real emission process,
* given a XComb driving the underlying Born
*/
vector<StdXCombPtr> makeRealXCombs(tStdXCombPtr bornXC);
/**
* Return true, if bookkeeping did not find a non-trivial setup.
*/
bool empty() const { return theSplittingMap.empty()&&theMergingMap.empty(); }
/**
* Return the emitter as referred to by the real emission
* matrix element.
*/
int realEmitter() const { return theRealEmitter; }
/**
* Set the emitter as referred to by the real emission
* matrix element.
*/
void realEmitter(int id) { theRealEmitter = id; }
/**
* Return the emission as referred to by the real emission
* matrix element.
*/
int realEmission() const { return theRealEmission; }
/**
* Set the emission as referred to by the real emission
* matrix element.
*/
void realEmission(int id) { theRealEmission = id; }
/**
* Return the spectator as referred to by the real emission
* matrix element.
*/
int realSpectator() const { return theRealSpectator; }
/**
* Set the spectator as referred to by the real emission
* matrix element.
*/
void realSpectator(int id) { theRealSpectator = id; }
/**
* Return the emitter as referred to by the underlying
* Born process.
*/
int bornEmitter() const { return theBornEmitter; }
/**
* Set the emitter as referred to by the underlying
* Born process.
*/
void bornEmitter(int id) { theBornEmitter = id; }
/**
* Return the spectator as referred to by the underlying
* Born process.
*/
int bornSpectator() const { return theBornSpectator; }
/**
* Set the spectator as referred to by the underlying
* Born process.
*/
void bornSpectator(int id) { theBornSpectator = id; }
/**
* Define the real emission key type
*/
typedef pair<pair<cPDVector,int>,pair<int,int> > RealEmissionKey;
/**
* Create a real emission key
*/
static RealEmissionKey realEmissionKey(const cPDVector& proc,
int em, int emm, int sp) {
return make_pair(make_pair(proc,emm),make_pair(em,sp));
}
/**
* Return the diagram of a real emission key
*/
static const cPDVector& process(const RealEmissionKey& key) {
return key.first.first;
}
/**
* Return the emission id of a real emission key
*/
static int emission(const RealEmissionKey& key) {
return key.first.second;
}
/**
* Return the emitter id of a real emission key
*/
static int emitter(const RealEmissionKey& key) {
return key.second.first;
}
/**
* Return the spectator id of a real emission key
*/
static int spectator(const RealEmissionKey& key) {
return key.second.second;
}
/**
* Define the underlying Born key type
*/
typedef pair<cPDVector,pair<int,int> > UnderlyingBornKey;
/**
* Create a underlying Born key
*/
static UnderlyingBornKey underlyingBornKey(const cPDVector& proc,
int em, int sp) {
return make_pair(proc,make_pair(em,sp));
}
/**
* Return the diagram of a underlying Born key
*/
static const cPDVector& process(const UnderlyingBornKey& key) {
return key.first;
}
/**
* Return the emitter id of a underlying Born key
*/
static int emitter(const UnderlyingBornKey& key) {
return key.second.first;
}
/**
* Return the spectator id of a underlying Born key
*/
static int spectator(const UnderlyingBornKey& key) {
return key.second.second;
}
/**
* Define real emission key and index dictionary
* for partons not involved in the given dipole.
*/
typedef pair<RealEmissionKey,map<int,int> > RealEmissionInfo;
/**
* Define underlying Born key and index dictionary
* for partons not involved in the given dipole.
*/
typedef pair<UnderlyingBornKey,map<int,int> > UnderlyingBornInfo;
/**
* Return the merging map
*/
const map<RealEmissionKey,UnderlyingBornInfo>& mergingMap() const { return theMergingMap; }
/**
* Return the splitting map
*/
const multimap<UnderlyingBornKey,RealEmissionInfo>& splittingMap() const { return theSplittingMap; }
/**
* Return the underlying Born diagrams to be considered
* for the given real emission process.
*/
const DiagramVector& underlyingBornDiagrams(const cPDVector& real) const;
/**
* Find the underlying Born diagram for the given real emission diagram
*/
tcDiagPtr underlyingBornDiagram(tcDiagPtr realDiag) const;
/**
* Return the real emission diagrams to be considered
* for the given Born process.
*/
const DiagramVector& realEmissionDiagrams(const cPDVector& born) const;
/**
* Find the real emission diagram for the given underlying Born diagram
*/
tcDiagPtr realEmissionDiagram(tcDiagPtr bornDiag) const;
/**
* Add all possible diagrams with the add() function.
*/
virtual void getDiagrams() const;
/**
* Return true, if this matrix element does not want to
* make use of mirroring processes; in this case all
* possible partonic subprocesses with a fixed assignment
* of incoming particles need to be provided through the diagrams
* added with the add(...) method.
*/
virtual bool noMirror () const { return true; }
/**
* With the information previously supplied with the
* setKinematics(...) method, a derived class may optionally
* override this method to weight the given diagrams with their
* (although certainly not physical) relative probabilities.
*/
virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
/**
* Return a Selector with possible colour geometries for the selected
* diagram weighted by their relative probabilities.
*/
virtual Selector<const ColourLines *>
colourGeometries(tcDiagPtr diag) const;
/**
* Select a ColpurLines geometry. The default version returns a
* colour geometry selected among the ones returned from
* colourGeometries(tcDiagPtr).
*/
virtual const ColourLines &
selectColourGeometry(tcDiagPtr diag) const;
/**
* Return the order in \f$\alpha_S\f$ in which this matrix element
* is given.
*/
virtual unsigned int orderInAlphaS() const { return realEmissionME()->orderInAlphaS(); }
/**
* Return the order in \f$\alpha_{EM}\f$ in which this matrix
* element is given. Returns 0.
*/
virtual unsigned int orderInAlphaEW() const { return underlyingBornME()->orderInAlphaEW(); }
//@}
/** @name Phasespace generation */
//@{
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr xc);
/**
* Set the typed and momenta of the incoming and outgoing partons to
* be used in subsequent calls to me() and colourGeometries()
* according to the associated XComb object. If the function is
* overridden in a sub class the new function must call the base
* class one first.
*/
virtual void setKinematics();
/**
* Generate internal degrees of freedom given nDim() uniform random
* numbers in the interval ]0,1[. To help the phase space generator,
* the 'dSigHatDR' should be a smooth function of these numbers,
* although this is not strictly necessary. The return value should
* be true of the generation succeeded. If so the generated momenta
* should be stored in the meMomenta() vector.
*/
virtual bool generateKinematics(const double * r);
/**
* The number of internal degreed of freedom used in the matrix
* element. This default version returns 0;
*/
virtual int nDim() const;
/**
* Return true, if this matrix element expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS () const { return realEmissionME()->wantCMS(); }
/**
* Clear the information previously provided by a call to
* setKinematics(...).
*/
virtual void clearKinematics();
/**
* If this is a dependent matrix element in a ME group, return true,
* if cuts should be ignored.
*/
virtual bool ignoreCuts() const { return theIgnoreCuts; }
/**
* Indicate that cuts should be ignored
*/
void doIgnoreCuts(bool is = true) { theIgnoreCuts = is; }
//@}
/** @name Tilde kinematics */
//@{
/**
* Return the TildeKinematics object used
*/
Ptr<TildeKinematics>::tcptr tildeKinematics() const { return theTildeKinematics; }
/**
* Set the TildeKinematics object used
*/
void tildeKinematics(Ptr<TildeKinematics>::tptr);
/**
* Generate the tilde kinematics from real emission
* kinematics accessible through the XComb's
* head object and store it in meMomenta(). This default
* implemenation uses the tildeKinematics() object.
*/
virtual bool generateTildeKinematics();
/**
* Return the InvertedTildeKinematics object used
*/
Ptr<InvertedTildeKinematics>::tcptr invertedTildeKinematics() const { return theInvertedTildeKinematics; }
/**
* Set the InvertedTildeKinematics object used
*/
void invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr);
/**
* Return the number of additional random numbers
* needed to generate real emission kinematics off
* the tilde kinematics previously supplied through
* the XComb object. This default implementation
* returns invertedTildeKinematics()->nDimRadiation()
*/
virtual int nDimRadiation() const;
/**
* Generate the real emission kinematics
* off the Born kinematics accessible through the XComb's
* head object and store it in meMomenta(); store
* the single particle phasespace in units of lastHeadXComb()->lastSHat()
* in jacobian(). This default
* implemenation uses the invertedTildeKinematics() object
*/
virtual bool generateRadiationKinematics(const double *);
/**
* Set a pt cut when splitting
*/
void ptCut(Energy cut);
/**
* Return the relevant dipole scale
*/
Energy lastDipoleScale() const {
return splitting() ? theLastSplittingScale : theLastSubtractionScale;
}
/**
* Return the relevant pt
*/
Energy lastPt() const {
return splitting() ? theLastSplittingPt : theLastSubtractionPt;
}
/**
* Return the relevant momentum fractions
*/
double lastZ() const {
return splitting() ? theLastSplittingZ : theLastSubtractionZ;
}
double lastRealR() const{
return theLastRealR;
}
double lastBornR() const{
return theLastBornR;
}
/**
* Return true, if this dipole acts in splitting mode.
*/
bool splitting() const { return theSplitting; }
/**
* Switch on splitting mode for this dipole.
*/
void doSplitting() { theSplitting = true; }
/**
* Switch off splitting mode for this dipole.
*/
void doSubtraction() { theSplitting = false; }
/**
* Return the subtraction parameters.
*/
const vector<double>& subtractionParameters() const { return theSubtractionParameters; }
/**
* Access the subtraction parameters.
*/
vector<double>& subtractionParameters() { return theSubtractionParameters; }
/**
* Return the shower hard scale encountered
*/
Energy showerHardScale() const { return theShowerHardScale; }
/**
* Set the shower hard scale encountered
*/
void showerHardScale(Energy s) { theShowerHardScale = s; }
/**
* Return the shower evolution scale encountered
*/
Energy showerScale() const { return theShowerScale; }
/**
* Set the shower evolution scale encountered
*/
void showerScale(Energy s) { theShowerScale = s; }
/**
* Return the shower splitting variables encountered
*/
const vector<double>& showerParameters() const { return theShowerParameters; }
/**
* Access the shower splitting variables encountered
*/
vector<double>& showerParameters() { return theShowerParameters; }
/**
* Return true, if this configuration is in the shower phase space
*/
bool isInShowerPhasespace() const { return theIsInShowerPhasespace; }
/**
* Indicate whether this configuration is in the shower phase space
*/
void isInShowerPhasespace(bool yes) { theIsInShowerPhasespace = yes; }
double jacobianMerging(Energy2,Energy2,int n) const;
/**
* Return true, if this configuration is above the shower infrared cutoff
*/
bool isAboveCutoff() const { return theIsAboveCutoff; }
/**
* Indicate whether this configuration is above the shower infrared cutoff
*/
void isAboveCutoff(bool yes) { theIsAboveCutoff = yes; }
//@}
/** @name Scale choices, couplings and PDFs */
//@{
/**
* Return true, if scales should be calculated from real emission kinematics
*/
bool realEmissionScales() const { return theRealEmissionScales; }
/**
* Switch on or off that scales should be calculated from real emission kinematics
*/
void doRealEmissionScales(bool on = true) { theRealEmissionScales = on; }
/**
* Return the scale associated with the phase space point provided
* by the last call to setKinematics().
*/
virtual Energy2 scale() const {
return realEmissionScales() ?
realEmissionME()->scale() :
underlyingBornME()->scale();
}
/**
* Return the value of \f$\alpha_S\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaS(scale()).
*/
virtual double alphaS() const {
return realEmissionScales() ?
realEmissionME()->alphaS() :
underlyingBornME()->alphaS();
}
/**
* Return the value of \f$\alpha_EM\f$ associated with the phase
* space point provided by the last call to setKinematics(). This
* versions returns SM().alphaEM(scale()).
*/
virtual double alphaEM() const {
return realEmissionScales() ?
realEmissionME()->alphaEM() :
underlyingBornME()->alphaEM();
}
/**
* Return true, if this matrix element provides the PDF
* weight for the first incoming parton itself.
*/
virtual bool havePDFWeight1() const { return realEmissionME()->havePDFWeight1(); }
/**
* Return true, if this matrix element provides the PDF
* weight for the second incoming parton itself.
*/
virtual bool havePDFWeight2() const { return realEmissionME()->havePDFWeight2(); }
/**
* How to sample the z-distribution.
* FlatZ = 1
* OneOverZ = 2
* OneOverOneMinusZ = 3
* OneOverZOneMinusZ = 4
*/
virtual int samplingZ() const {return 4;}
//@}
/** @name Matrix elements and evaluation */
//@{
/**
* Return the real emission matrix element
*/
Ptr<MatchboxMEBase>::tcptr realEmissionME() const {
return theRealEmissionME;
}
/**
* Return the real emission matrix element
*/
Ptr<MatchboxMEBase>::tptr realEmissionME() {
return theRealEmissionME;
}
/**
* Set the real emission matrix element
*/
void realEmissionME(Ptr<MatchboxMEBase>::tptr me) { theRealEmissionME = me; }
/**
* Return the underlying Born matrix element
*/
Ptr<MatchboxMEBase>::tcptr underlyingBornME() const {
return theUnderlyingBornME;
}
/**
* Return the underlying Born matrix element
*/
Ptr<MatchboxMEBase>::tptr underlyingBornME() {
return theUnderlyingBornME;
}
/**
* Set the underlying Born matrix element
*/
void underlyingBornME(Ptr<MatchboxMEBase>::tptr me) { theUnderlyingBornME = me; }
/**
* Set the dipoles which have been found along with this dipole
*/
void partnerDipoles(const vector<Ptr<SubtractionDipole>::tptr>& p) {
thePartners = p;
}
/**
* Return the dipoles which have been found along with this dipole
*/
const vector<Ptr<SubtractionDipole>::tptr>& partnerDipoles() const {
return thePartners;
}
/**
* Return the matrix element averaged over spin correlations.
*/
virtual double me2Avg(double ccme2) const = 0;
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR(Energy2 factorizationScale) const;
/**
* Return the matrix element squared differential in the variables
* given by the last call to generateKinematics().
*/
virtual CrossSection dSigHatDR() const { return dSigHatDR(ZERO); }
//TODO
bool clustersafe() const;
bool clustersafe();
CrossSection ps(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;
- CrossSection dipMinusPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;
+ pair<double,double> dipandPs(Energy2 factorizationScale,Ptr<ColourBasis>::tptr largeNBasis) const;
//@}
/** @name Methods relevant to matching */
//@{
/**
* Set the shower approximation.
*/
void showerApproximation(Ptr<ShowerApproximation>::tptr app) {
theShowerApproximation = app;
}
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const { return theShowerApproximation; }
/**
* Indicate that the shower real emission contribution should be subtracted.
*/
void doRealShowerSubtraction() { theRealShowerSubtraction = true; }
/**
* Return true, if the shower real emission contribution should be subtracted.
*/
bool realShowerSubtraction() const { return theRealShowerSubtraction; }
/**
* Indicate that the shower virtual contribution should be subtracted.
*/
void doVirtualShowerSubtraction() { theVirtualShowerSubtraction = true; }
/**
* Return true, if the shower virtual contribution should be subtracted.
*/
bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
/**
* Indicate that the loopsim matched virtual contribution should be subtracted.
*/
void doLoopSimSubtraction() { theLoopSimSubtraction = true; }
/**
* Return true, if the loopsim matched virtual contribution should be subtracted.
*/
bool loopSimSubtraction() const { return theLoopSimSubtraction; }
//@}
/** @name Caching and diagnostic information */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches();
/**
* Indicate that the subtraction is being tested.
*/
void doTestSubtraction() { theSubtractionTest = true; }
/**
* Return true, if the subtraction is being tested.
*/
bool testSubtraction() const { return theSubtractionTest; }
/**
* Return true, if verbose
*/
bool verbose() const { return realEmissionME()->verbose() || underlyingBornME()->verbose(); }
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Print debug information on the last event
*/
virtual void printLastEvent(ostream&) const;
/**
* Write out diagnostic information for
* generateTildeKinematics
*/
void logGenerateTildeKinematics() const;
/**
* Write out diagnostic information for
* generateRadiationKinematics
*/
void logGenerateRadiationKinematics(const double * r) const;
/**
* Write out diagnostic information for
* me2 evaluation
*/
void logME2() const;
/**
* Write out diagnostic information
* for dsigdr evaluation
*/
void logDSigHatDR(double effectiveJac) const;
//@}
/** @name Reweight objects */
//@{
/**
* Insert a reweight object
*/
void addReweight(Ptr<MatchboxReweightBase>::ptr rw) { theReweights.push_back(rw); }
/**
* Return the reweight objects
*/
const vector<Ptr<MatchboxReweightBase>::ptr>& reweights() const { return theReweights; }
/**
* Access the reweight objects
*/
vector<Ptr<MatchboxReweightBase>::ptr>& reweights() { return theReweights; }
//@}
/** @name Methods used to setup SubtractionDipole objects */
//@{
/**
* Clone this dipole.
*/
Ptr<SubtractionDipole>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
void cloneDependencies(const std::string& prefix = ""){cloneDependencies(prefix,false);};
void cloneDependencies(const std::string& prefix , bool slim);
//@}
/** @name Methods required to setup the event record */
//@{
/**
* construct the spin information for the interaction
*/
virtual void constructVertex(tSubProPtr sub);
/**
* construct the spin information for the interaction
*/
virtual void constructVertex(tSubProPtr sub, const ColourLines* cl);
/**
* Comlete a SubProcess object using the internal degrees of freedom
* generated in the last generateKinematics() (and possible other
* degrees of freedom which was intergated over in dSigHatDR(). This
* default version does nothing. Will be made purely virtual in the
* future.
*/
virtual void generateSubCollision(SubProcess & sub);
/**
* Alpha parameter as in Nagy
* (http://arxiv.org/pdf/hep-ph/0307268v2.pdf) to restrict dipole
* phase space
*/
double alpha() const;
/*
* True if phase space point is above the alpha cut for this dipole.
*/
virtual bool aboveAlpha() const{return true;}
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
/**
* The factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tptr theFactory;
/**
* Wether or not this dipole acts in splitting mode.
*/
bool theSplitting;
/**
* True, if should apply to process in the xcomb.
*/
bool theApply;
/**
* True, if the subtraction is being tested.
*/
bool theSubtractionTest;
/**
* True if cuts should be ignored
*/
bool theIgnoreCuts;
/**
* The real emission matrix element to be considered
*/
Ptr<MatchboxMEBase>::ptr theRealEmissionME;
/**
* The underlying Born matrix element
*/
Ptr<MatchboxMEBase>::ptr theUnderlyingBornME;
/**
* The dipoles which have been found along with this dipole
*/
vector<Ptr<SubtractionDipole>::tptr> thePartners;
/**
* The TildeKinematics to be used.
*/
Ptr<TildeKinematics>::ptr theTildeKinematics;
/**
* The InvertedTildeKinematics to be used.
*/
Ptr<InvertedTildeKinematics>::ptr theInvertedTildeKinematics;
/**
* A vector of reweight objects the sum of which
* should be applied to reweight this dipole
*/
vector<Ptr<MatchboxReweightBase>::ptr> theReweights;
/**
* The emitter as referred to by the real emission
* matrix element.
*/
int theRealEmitter;
/**
* The emission as referred to by the real emission
* matrix element.
*/
int theRealEmission;
/**
* The spectator as referred to by the real emission
* matrix element.
*/
int theRealSpectator;
/**
* The subtraction parameters
*/
vector<double> theSubtractionParameters;
/**
* Map real emission diagrams to underlying Born diagrams
* and tilde emitter/spectator.
*/
map<RealEmissionKey,UnderlyingBornInfo> theMergingMap;
/**
* Map underlying Born diagrams and tilde emitter/spectator
* to real emission diagram containing the splitting.
*/
multimap<UnderlyingBornKey,RealEmissionInfo> theSplittingMap;
/**
* Map underlying Born diagrams to emitter/spectator pairs
*/
map<cPDVector,pair<int,int> > theIndexMap;
/**
* Map real emission processes to Born diagrams
*/
map<cPDVector,DiagramVector> theUnderlyingBornDiagrams;
/**
* Map Born processes to real emission diagrams
*/
map<cPDVector,DiagramVector> theRealEmissionDiagrams;
/**
* Map underlying Born diagrams to real emission diagrams.
*/
map<tcDiagPtr,tcDiagPtr> theBornToRealDiagrams;
/**
* Map real emission diagrams to underlying Born diagrams.
*/
map<tcDiagPtr,tcDiagPtr> theRealToBornDiagrams;
/**
* The last real emission key encountered
*/
RealEmissionKey lastRealEmissionKey;
/**
* The last underlying Born key encountered
*/
UnderlyingBornKey lastUnderlyingBornKey;
/**
* The last real emission info encountered
*/
multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator lastRealEmissionInfo;
/**
* The emitter as referred to by the underlying Born
* matrix element.
*/
int theBornEmitter;
/**
* The spectator as referred to by the underlying Born
* matrix element.
*/
int theBornSpectator;
/**
* The last scale as generated from the tilde mapping
*/
Energy theLastSubtractionScale;
/**
* The last scale as generated from the splitting mapping
*/
Energy theLastSplittingScale;
/**
* The last pt as generated from the tilde mapping
*/
Energy theLastSubtractionPt;
/**
* The last pt as generated from the splitting mapping
*/
Energy theLastSplittingPt;
/**
* The last z as generated from the tilde mapping
*/
double theLastSubtractionZ;
/**
* The last z as generated from the splitting mapping
*/
double theLastSplittingZ;
double theLastRealR;
double theLastBornR;
/**
* The shower approximation.
*/
Ptr<ShowerApproximation>::ptr theShowerApproximation;
/**
* True, if the shower real emission contribution should be subtracted.
*/
bool theRealShowerSubtraction;
/**
* True, if the shower virtual contribution should be subtracted.
*/
bool theVirtualShowerSubtraction;
/**
* True, if the loopsim matched virtual contribution should be subtracted.
*/
bool theLoopSimSubtraction;
/**
* True, if scales should be calculated from real emission kinematics
*/
bool theRealEmissionScales;
/**
* Return the shower hard scale encountered
*/
Energy theShowerHardScale;
/**
* The shower evolution scale encountered
*/
Energy theShowerScale;
/**
* The shower splitting variables encountered
*/
vector<double> theShowerParameters;
/**
* True, if this configuration is in the shower phase space
*/
bool theIsInShowerPhasespace;
/**
* True, if this configuration is above the shower infrared cutoff
*/
bool theIsAboveCutoff;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SubtractionDipole & operator=(const SubtractionDipole &);
};
}
#endif /* HERWIG_SubtractionDipole_H */
diff --git a/src/Merging/MergingLEP.in b/src/Merging/MergingLEP.in
--- a/src/Merging/MergingLEP.in
+++ b/src/Merging/MergingLEP.in
@@ -1,191 +1,191 @@
set /Herwig/Particles/d:NominalMass 0*GeV
set /Herwig/Particles/dbar:NominalMass 0*GeV
set /Herwig/Particles/u:NominalMass 0*GeV
set /Herwig/Particles/ubar:NominalMass 0*GeV
set /Herwig/Particles/s:NominalMass 0*GeV
set /Herwig/Particles/sbar:NominalMass 0*GeV
set /Herwig/Particles/c:NominalMass 0*GeV
set /Herwig/Particles/cbar:NominalMass 0*GeV
set /Herwig/Particles/b:NominalMass 0*GeV
set /Herwig/Particles/bbar:NominalMass 0*GeV
set /Herwig/Particles/e+:PDF /Herwig/Partons/NoPDF
set /Herwig/Particles/e-:PDF /Herwig/Partons/NoPDF
set /Herwig/DipoleShower/Kinematics/FFLightKinematics:IRCutoff 0.78
set /Herwig/Couplings/NLOAlphaS:min_active_flavours 4
set /Herwig/Couplings/NLOAlphaS:input_alpha_s 0.118
set /Herwig/Model:QCD/RunningAlphaS /Herwig/Couplings/NLOAlphaS
set /Herwig/DipoleShower/DipoleShowerHandler:GlobalAlphaS /Herwig/Couplings/NLOAlphaS
cd /Herwig/Merging
create Herwig::Merger Merger
set Merger:DipoleShowerHandler /Herwig/DipoleShower/DipoleShowerHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper Merger
cd /Herwig/DipoleShower/Kernels
clear /Herwig/DipoleShower/DipoleShowerHandler:Kernels
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ggxDipoleKernel
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFqx2qgxDipoleKernel
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ddxDipoleKernel
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2uuxDipoleKernel
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ccxDipoleKernel
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2ssxDipoleKernel
insert /Herwig/DipoleShower/DipoleShowerHandler:Kernels 0 FFgx2bbxDipoleKernel
cd /Herwig/EventHandlers
set /Herwig/Generators/LEPGenerator:EventHandler:LuminosityFunction:Energy 91.2*GeV
cd /Herwig/Merging
insert /Herwig/Generators/LEPGenerator:EventHandler:SubProcessHandlers[0] EEMFactory
set EEMFactory:OrderInAlphaS 0
set EEMFactory:OrderInAlphaEW 2
do EEMFactory:Process e- e+ j j
create Herwig::SimpleColourBasis LargeNColourBasis
set LargeNColourBasis:LargeN On
set Merger:LargeNBasis LargeNColourBasis
cd /Herwig/Generators
set LEPGenerator:EventHandler:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler
set LEPGenerator:EventHandler:CascadeHandler:MPIHandler NULL
set LEPGenerator:EventHandler:MultipleInteractionHandler NULL
set /Herwig/Samplers/Sampler:Verbose On
set /Herwig/Cuts/JetKtCut:MinKT 0.0*GeV
set /Herwig/Generators/LEPGenerator:EventHandler:CollisionCuts Off
create Herwig::MonacoSampler /Herwig/Samplers/Monaco
set /Herwig/Samplers/Sampler:BinSampler /Herwig/Samplers/Monaco
cd /Herwig/MatrixElements/Matchbox/Phasespace
set TreePhasespace:M0 0.001*GeV
set TreePhasespace:MC 0.00001*GeV
set /Herwig/Merging/EEMFactory:Phasespace TreePhasespace
cd /Herwig/EventHandlers
set LEPHandler:Sampler /Herwig/Samplers/Sampler
cd /Herwig/MatrixElements/Matchbox
cd Utility
insert DiagramGenerator:ExcludeInternal 0 /Herwig/Particles/gamma
cd /Herwig/Samplers
set Monaco:EnhancementFactor 1.2
set Monaco:InitialPoints 1000
set Monaco:LuminosityMapperBins 0
set Monaco:NIterations 4
set Monaco:RemapChannelDimension Yes
set Monaco:RemapperMinSelection 0.0001
set Monaco:RemapperPoints 1000
set Monaco:UseAllIterations No
cd /Herwig/EventHandlers
cd /Herwig/Samplers
set Sampler:UpdateAfter 3000
set Sampler:AddUpSamplers Off
set Sampler:GlobalMaximumWeight Off
set Sampler:FlatSubprocesses Off
set Sampler:MinSelection 0.00001
set Sampler:AlmostUnweighted On
-set Sampler:BinSampler:Kappa 0.001
+set Sampler:BinSampler:Kappa 0.1
set Sampler:RunCombinationData Off
set Sampler:Verbose On
set /Herwig/EventHandlers/LEPHandler:Sampler /Herwig/Samplers/Sampler
cd /Herwig/MatrixElements/Matchbox/Scales/
set /Herwig/Merging/MScale:ScaleChoice FixedScale
set FixedScale:FixedScale 91.2*GeV
cd /Herwig/Generators
#set LEPGenerator:EventHandler:DecayHandler NULL
#set LEPGenerator:EventHandler:HadronizationHandler NULL
cd /Herwig/Analysis
clear /Herwig/Generators/LEPGenerator:AnalysisHandlers
create ThePEG::RivetAnalysis Rivet RivetAnalysis.so
insert /Herwig/Generators/LEPGenerator:AnalysisHandlers 0 Rivet
insert Rivet:Analyses 0 ALEPH_2004_S5765862
insert Rivet:Analyses 0 ALEPH_1991_S2435284
insert Rivet:Analyses 0 ALEPH_1996_S3486095
insert Rivet:Analyses 0 DELPHI_1995_S3137023
insert Rivet:Analyses 0 DELPHI_1996_S3430090
insert Rivet:Analyses 0 ALEPH_2001_S4656318
insert Rivet:Analyses 0 DELPHI_2002_069_CONF_603
insert Rivet:Analyses 0 JADE_OPAL_2000_S4300807
insert Rivet:Analyses 0 OPAL_1998_S3780481
insert Rivet:Analyses 0 OPAL_2001_S4553896
insert Rivet:Analyses 0 OPAL_2002_S5361494
insert Rivet:Analyses 0 OPAL_2004_S6132243
insert Rivet:Analyses 0 SLD_1996_S3398250
insert Rivet:Analyses 0 MC_XS
cd /Herwig/Generators
set /Herwig/Generators/LEPGenerator:IntermediateOutput Yes
set /Herwig/Generators/LEPGenerator:EventHandler /Herwig/EventHandlers/LEPHandler
set /Herwig/DipoleShower/DipoleShowerHandler:MaxPtIsMuF Yes
cd /Herwig/Cuts
set EECuts:Fuzzy FuzzyTheta
set MatchboxJetMatcher:Factory /Herwig/Merging/EEMFactory
cd /Herwig/Merging
set EEMFactory:MergingHelper Merger
set Merger:MFactory EEMFactory
set MScale:MergingHelper Merger
set Merger:MergingJetFinder /Herwig/Cuts/JetFinder
cd /Herwig/Generators
set LEPGenerator:DebugLevel 1
set LEPGenerator:PrintEvent 10
set /Herwig/DipoleShower/DipoleShowerHandler:MergingHelper:minusL No

File Metadata

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

Event Timeline