No OneTemporary
Actions

Size

81 KB

Subscribers

None

View Options

	diff --git a/DipoleShower/Merging2/Merger.cc b/DipoleShower/Merging2/Merger.cc
	--- a/DipoleShower/Merging2/Merger.cc
	+++ b/DipoleShower/Merging2/Merger.cc
	@@ -1,1498 +1,1482 @@
	// -- 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;
	theKImproved=true;
	theSmearing=0.;
	theDipoleShowerHandler=
	Ptr<DipoleShowerHandler>::ptr();
	isUnitarized=true;
	isNLOUnitarized=true;
	defMERegionByJetAlg=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-> getLongestHistory_simple(true,xiQSh);//Longest irreführend
	+ 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-> getLongestHistory_simple(false,xiQSh);
	+ 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().weightalphaReweight()*pdfReweight();
	if(weight==0.&&weightCL==0.)return 0.;
	- //Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);

	res+=weight*matrixElementWeight(startscale,Node,(!maxMulti&&!projected)?theGamma:1.);
	}

	-
	-
	- //cout<<"\n"<<((!maxMulti&&!projected)?gamma:1.)
	- //<<history.back().weight<<" "<<alphaReweight()<<" "<<pdfReweight();;
	-
	-
	if(CLNode&&theGamma!=1.){
	- //assert(false);
	Energy startscale=CKKW_StartScale(BornCL);
	Energy projectedscale=startscale;
	fillHistory( startscale, BornCL, CLNode,fast);
	if (!fillProjector(projectedscale))return 0.;
	Node->runningPt(projectedscale);
	weightCL=history.back().weightalphaReweight()*pdfReweight();
	double diff=0;
	Node->nodeME()->factory()->setAlphaParameter(1.);
	diff-=weightCLCLNode->dipol()->dSigHatDR(sqr(startscalexiFacME))/nanobarn;
	Node->nodeME()->factory()->setAlphaParameter(theGamma);

	string alp=(CLNode->dipol()->aboveAlpha()?"Above":"Below");

	diff+=weightCLCLNode->dipol()->dSigHatDR(sqr(startscalexiFacME))/nanobarn;
	Node->nodeME()->factory()->setAlphaParameter(1.);
	- //cout<<"\n diff "<<diff;
	- //cout<<"\ndiff "<<diff<<" "<<history.back().weight<<" "<<alphaReweight()<<" "<<alp;
	+
	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-> getLongestHistory_simple(true,xiQSh);//Longest irreführend
	+ 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()){
	-
	- if (Born->parent()) {
	- //cout<<"\n parent not pass";
	- }
	- return 0.;
	-
	+ 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().weightalphaReweight()*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-> getLongestHistory_simple(true,xiQSh);
	+ 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().weightalphaReweight()*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(startscalexiRenSh)/SM().alphaS()
	(matrixElement+unlopsweight);
	}


	double Merger::reweightCKKWRealStandard(NPtr Node,bool fast){
	bool allAbove=Node->allAbove(mergePt());
	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-> getLongestHistory_simple(true,xiQSh);
	+ NPtr Born= Node-> getHistory(true,xiQSh);
	NPtr CLNode= Node->randomChild();
	bool inhist=CLNode->isInHistoryOf(Born);
	- Born=CLNode-> getLongestHistory_simple(false,xiQSh);
	+ 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().weightalphaReweight()*pdfReweight();
	if(weight==0.)return 0.;
	bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
	Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);

	double res= weightas(startscalexiRenSh)/SM().alphaS()*
	(double)Node->children().size()*
	((inhist?matrixElementWeight(startscale,Node):0.)
	+CLNode->dipol()->dSigHatDR(sqr(startscale*xiFacME))/nanobarn);
	// cout<<"\nCLNode->dipol()->dSigHatDR(sqr(startscalexiFacME))/nanobarn "<<CLNode->dipol()->dSigHatDR(sqr(startscalexiFacME))/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-> getLongestHistory_simple(false,xiQSh);
	+ 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().weightalphaReweight()*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(startscalexiFacME);
	}
	//double me=matrixElementWeight(startscale,Node);
	//cout<<"\nreweightCKKWRealBelowSubReal "<<me<<" "<<sumPS<<" "<<me/sumPS;
	return weightas(startscalexiRenSh)/SM().alphaS()*
	(matrixElementWeight(startscale,Node)-sumPS);
	}



	double Merger::reweightCKKWRealBelowSubInt(NPtr Node,bool fast){
	NPtr CLNode= Node->randomChild();
	- NPtr Born=CLNode-> getLongestHistory_simple(false,xiQSh);
	+ 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().weightalphaReweight()*pdfReweight();
	if(weight==0.)return 0.;
	bool maxMulti=CLNode->xcomb()->meMomenta().size() == maxLegsNLO();
	Node->vetoPt((projected&&maxMulti)?mergePt():history.back().scale);

	double res=0.;
	if (Born==CLNode&&!CLNode->children().empty())
	res=-1.CLNode->dipol()->dSigHatDR(sqr(startscalexiFacME))/nanobarn;
	else
	res=CLNode->calcPsMinusDip(startscale*xiFacME);

	return weightas(startscalexiRenSh)/SM().alphaS()*
	(double)Node->children().size()*res;
	}










	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].scalexiFacME, 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()
	(theKImproved?(1.+((3.(67./18.-1./6.Constants::piConstants::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 += CF2.(PDFxByzparton - zPDFxparton)/(z(1.-z));
	restmp -= CFPDFxByzparton (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/xfactor;
	}
	}

	restmp += ( (11./6.) * CA - (1./3.)theNf + 2.CAlog(1.-x) ) PDFxparton;
	if ( z > x ) {
	restmp += 2. * CA * ( PDFxByzgluon - zPDFxparton ) / (z(1.-z));
	restmp += 2.* CA ( (1.-z)/z - 1. + z(1.-z) ) * PDFxByzgluon / z;
	}

	}
	if (PDFxparton<1e-8)restmp= 0.;
	res+=1restmplog(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());


	- //cout<<"\nfast ";//<<fast<<" HS "<<history.size()<<" stnode "<<StartingBorn<<" "<<eH->didEstimate();
	-
	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);
	}else if(ME->oneLoopNoBorn()){
	res*=reweightCKKWVirtualStandard(Node);
	}else if(theGamma!=1.){
	res*=reweightCKKWBornGamma(Node,fast);
	}else{

	res*=reweightCKKWBornStandard(Node,fast);
	}



	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)))+
	(theKImproved?(((3.(67./18.-1./6.Constants::piConstants::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){
	Born->deepHead()->nodeME()->mePartonData().size();
	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));};



	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.(emissionmomemittermom + emissionmomspectatormom + emittermomspectatormom));
	double y = emissionmomemittermom / (emissionmomemittermom + emissionmomspectatormom + emittermomspectatormom);
	double z = emittermomspectatormom / (emittermomspectatormom + emissionmom*spectatormom);
	if (abs(scale * sqrt(yz(1.-z))-winnerScale)<0.0001*GeV) {
	foundwinnerpt=true;
	}
	// if(scale * sqrt(yz(1.-z))<optVeto&&winnerScale>optVeto)cout<<"\nFF "<<(scale * sqrt(yz(1.-z))/GeV);
	ptx =min(ptx,scale * sqrt(yz(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.EmspEmsp);
	double x = (- emissionmomemittermom + emissionmomspectatormom + emittermomspectatormom) /(emittermomspectatormom + emissionmom*spectatormom);
	double z = emittermomspectatormom / (emittermomspectatormom + 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.0001GeV) {
	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.EmspEmsp);
	double x = (-emissionmomspectatormom + emittermomspectatormom + emittermomemissionmom) /(emittermomemissionmom + emittermom*spectatormom);
	double u = emittermomemissionmom / (emittermomemissionmom + emittermom*spectatormom);
	if (abs(scale sqrt(u(1.-u)(1.-x)/x)-winnerScale)<0.0001GeV) {
	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 = (emittermomspectatormom - emittermomemissionmom - spectatormomemissionmom)/(emittermomspectatormom);
	double v = (emittermomemissionmom)/(emittermomspectatormom);
	if (abs(scale * sqrt(v(1.-x-v)/x)-winnerScale)<0.0001GeV) {
	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<<
	theKImproved<< projected<<
	isUnitarized<< isNLOUnitarized<<
	defMERegionByJetAlg<<
	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>>
	theKImproved>> projected>>
	isUnitarized>> isNLOUnitarized>>
	defMERegionByJetAlg>>
	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>
	interfacetheKImproved ("KImproved","",&Merger::theKImproved, false, false, false);
	static SwitchOption interfacetheKImprovedYes
	(interfacetheKImproved,"Yes","",true);
	static SwitchOption interfacetheKImprovedNo
	(interfacetheKImproved,"No","",false);








	static Parameter<Merger,Energy> interfaceMergerScale
	("MergingScale",
	"The MergingScale.",
	&Merger::theCentralMergePt, GeV, 20.0GeV, 0.0GeV, 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 Parameter<Merger, Energy> interfaceIRSafePT("IRSafePT", "Set the pt for which a matrix element 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> interfaceadditionalN("additionalN", "Set the number of additional jets to consider.", &Merger::theN, 0, 0,
	0, false, false, Interface::lowerlim);

	static Parameter<Merger, int> interfacevirtualM("virtualM",
	"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/Merging2/Node.cc b/DipoleShower/Merging2/Node.cc
	--- a/DipoleShower/Merging2/Node.cc
	+++ b/DipoleShower/Merging2/Node.cc
	@@ -1,607 +1,607 @@
	// -- 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())];
	}

	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(hardScaleFactorsqrt((*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
	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;
	}
	// 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::getLongestHistory_simple(bool normal,double hardScaleFactor) {
	+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 = getLongestHistory_simple(true,hardScaleFactor);
	+ 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;
	}

	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/Merging2/Node.h b/DipoleShower/Merging2/Node.h
	--- a/DipoleShower/Merging2/Node.h
	+++ b/DipoleShower/Merging2/Node.h
	@@ -1,434 +1,434 @@
	// -- 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);
	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();


	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 getLongestHistory_simple(bool normal=true,double hardscalefactor=1.);
	+ 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 */

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 1:25 PM (19 h, 36 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022948
Default Alt Text: (81 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions