No OneTemporary
Actions

Size

57 KB

Subscribers

None

View Options

	diff --git a/Shower/Dipole/Merging/Merger.cc b/Shower/Dipole/Merging/Merger.cc
	--- a/Shower/Dipole/Merging/Merger.cc
	+++ b/Shower/Dipole/Merging/Merger.cc
	@@ -1,1636 +1,1636 @@
	// -- C++ --
	//
	// Merger.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright ( C ) 2002-2007 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL , see COPYING for details.
	// Please respect the MCnet academic guidelines , see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined , non-templated member
	// functions of the Merger class.
	//

	#include "Merger.h"
	#include "Node.h"
	#include "MergingFactory.h"
	// other includes when needed below.
	using namespace Herwig;

	IBPtr Merger::clone() const {
	return new_ptr( *this );
	}

	IBPtr Merger::fullclone() const {
	return new_ptr( *this );
	}





	namespace {
	double decideClustering(const NodePtr sub,const NodePtr head,bool& pro){
	if( sub != head ){// at least one history step -> unitarisation
	if ( UseRandom::rndbool() ){ pro = true; return -2.; }
	else{ pro = false; return 2.; }
	} // no ordered history -> no projection
	else{ pro = false; return 1.; }
	}
	}

	CrossSection Merger::MergingDSigDRBornStandard( ){
	// get the history for the process
	const NodePtr productionNode =
	currentNode()-> getHistory( true, DSH()->hardScaleFactor() );
	// decide if to cluster
	weight = decideClustering(productionNode, currentNode(), projected);
	// check if we only want to calculate the current multiplicity.
	if(notOnlyMulti()) return ZERO;
	// Check for cuts on the production proces.
	if ( !productionNode->xcomb()->willPassCuts() ) return ZERO;
	// calculate the staring scale for the production node
	Energy startscale = CKKW_StartScale( productionNode );
	// fill history with caluclation of sudakov supression
	fillHistory( startscale , productionNode , currentNode() );
	// fill the projector -> return the scale of the last splitting
	currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
	// the weight has three components to get shower history weight
	weight = history.back().weight // Sudakov suppression
	alphaReweight()* // alpha_s reweight
	pdfReweight(); // pdf reweight
	// If weight is zero return.
	if( weight == ZERO ) return ZERO;
	//calculate the cross section
	return weight*TreedSigDR( startscale , 1. );
	}




	CrossSection Merger::MergingDSigDRVirtualStandard( ){
	// get the history for the process
	const NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() );
	// decide if to cluster
	weight = decideClustering(productionNode,currentNode(),projected);
	// Check for cuts on the production proces.
	if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
	// calculate the staring scale
	Energy startscale = CKKW_StartScale( productionNode );
	// fill history with caluclation of sudakov supression
	fillHistory( startscale , productionNode , currentNode() );
	// fill the projector -> return the scale of the last splitting
	currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
	// the weight has three components to get shower history weight
	double ww1 = history.back().weight;
	double ww2 = alphaReweight();
	double ww3 = pdfReweight();
	weight = ww1ww2*ww3;
	// If weight is zero return.
	if( weight == 0. )return ZERO;
	// calculate the cross section for virtual contribution
	// and various insertion operators.
	CrossSection matrixElement = LoopdSigDR( startscale );
	// Now calculate the expansion of the shower history.
	// first: the born contibution:
	CrossSection bornWeight = currentME()->dSigHatDRB();
	// second: expansion of pdf ,alpha_s-ratio and sudakov suppression.
	double w1 = -sumPdfReweightExpansion();
	double w2 = -sumAlphaSReweightExpansion();
	double w3 = -sumFillHistoryExpansion();
	// put together the expansion weights.
	CrossSection expansionweight = ( w1+w2+w3 )*(theShowerExpansionWeights?1.:0.)
	*bornWeight
	SM().alphaS()/( 2.ThePEG::Constants::pi );
	// [ DEBUG ]
	if ( currentNode()->legsize() == 5 && Debug::level > 2 )
	debugVirt(weight,w1,w2,w3,matrixElement,ww1,ww2,ww3,productionNode,bornWeight);
	// return with correction that ME was calculated with fixed alpha_s
	return weight* as( startscale*DSH()->renFac() )/
	SM().alphaS()* ( matrixElement+expansionweight );
	}




	CrossSection Merger::MergingDSigDRRealStandard(){
	if ( currentNode()->children().empty() ) {
	throw Exception()
	<< "Real emission contribution without underlying born."
	<< "These are finite contibutions already handled in LO merging."
	<< Exception::abortnow;
	}
	// check for IR Safe Cutoff
	if( !currentNode()->allAbove( theIRSafePT ) )return ZERO;

	auto inOutPair = currentNode()->getInOut();
	NodePtr randomChild = currentNode()->randomChild();
	bool meRegion =matrixElementRegion(
	inOutPair.first ,
	inOutPair.second ,
	randomChild->pT() ,
	theMergePt );


	if ( meRegion )return MergingDSigDRRealAllAbove( );
	if ( UseRandom::rndbool() )
	return 2.*MergingDSigDRRealBelowSubReal( );
	return 2.*MergingDSigDRRealBelowSubInt( );
	}

	CrossSection Merger::MergingDSigDRRealAllAbove( ){
	//If all dipoles pts are above , we cluster to this dipole.
	NodePtr CLNode = currentNode()->randomChild();
	// Check if phase space poing is in ME region--> else rm PSP
	if ( !CLNode->children().empty() ) {
	auto inOutPair = CLNode->getInOut();
	NodePtr randomChild = CLNode->randomChild();
	if( !matrixElementRegion( inOutPair.first ,
	inOutPair.second ,
	randomChild->pT() ,
	theMergePt ) )return ZERO;
	}

	// first find the history for the acctual Node
	NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() );
	// check if the randomly choosen CLNode is part of the history.
	// If CLNode is not part of the history , dont calculate the Real contribution
	// else multiply the real contribution with N ( number of children ).
	// this returns the sudakov suppression according to
	// the clustering of the born parts.
	bool inhist = CLNode->isInHistoryOf( productionNode );
	// get the history for the clustered Node.
	productionNode = CLNode-> getHistory( false , DSH()->hardScaleFactor() );
	// decide if to cluster
	weight = decideClustering(productionNode,CLNode,projected);
	// Check for cuts on the production process.
	if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
	// calculate the staring scale
	Energy startscale = CKKW_StartScale( productionNode );
	// fill history with caluclation of sudakov supression
	fillHistory( startscale , productionNode , CLNode );
	// fill the projector -> return the scale of the last splitting
	currentNode()->runningPt( fillProjector( projected ? 2 : 1 ) );
	// the weight has three components to get shower history weight
	weight = history.back().weightalphaReweight()*pdfReweight();
	if( weight == 0. )return ZERO;
	// The inhist flag produces the correct cluster sensity.
	CrossSection me = ( inhist?TreedSigDR( startscale ):ZERO );
	// calculate the dipole
	CrossSection dip = CLNode->calcDip( startscale* currentME()->renFac() );


	CrossSection res = weightas( startscaleDSH()->renFac() )/SM().alphaS()*
	currentNode()->children().size()*( me -dip );
	// [ DEBUG ]
	if ( currentNode()->legsize() == 6&&Debug::level > 2 )
	debugReal("RAA",weight,me,dip);

	return res;
	}

	CrossSection Merger::MergingDSigDRRealBelowSubReal( ){
	NodePtr HistNode = currentNode()->randomChild();
	if ( !HistNode->children().empty() ) {
	auto inOutPair = HistNode->getInOut();
	NodePtr randomChild = HistNode->randomChild(); // Here we make sure that clustering and splitting are invertible
	if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO;
	}

	const NodePtr productionNode = HistNode-> getHistory( false , DSH()->hardScaleFactor() );

	weight = decideClustering(productionNode,HistNode,projected);

	if ( !productionNode->xcomb()->willPassCuts() )return ZERO;

	Energy startscale = CKKW_StartScale( productionNode );
	fillHistory( startscale , productionNode , HistNode );
	currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
	weight = history.back().weightalphaReweight()*pdfReweight();

	if( weight == 0. )return ZERO;

	CrossSection sumPS = ZERO;

	for( auto const & child : currentNode()->children() ){
	if ( child->allAbove( mergePt() ) ){
	if( ( child )->pT()>mergePt()/theRealSubtractionRatio )
	sumPS += child->calcPs( startscale* currentME()->renFac() );
	else
	sumPS += child->calcDip( startscale* currentME()->renFac() );
	}else{
	assert( child->pT()>mergePt() );
	}
	}

	CrossSection me = TreedSigDR( startscale );
	// [ DEBUG ]
	if ( currentNode()->legsize() == 6&&Debug::level > 2 )
	debugReal("RBSR",weight,me,sumPS);
	//Here we subtract the PS ( and below the dynamical cutoff the Dip )
	return weightas( startscaleDSH()->renFac() )/SM().alphaS()*
	( me-sumPS );
	}



	CrossSection Merger::MergingDSigDRRealBelowSubInt( ){

	if( currentNode()->children().empty() )return ZERO;
	NodePtr CLNode = currentNode()->randomChild();
	if( CLNode->pT()<mergePt()/theRealSubtractionRatio )return ZERO;

	if ( !CLNode->children().empty() ) {
	auto inOutPair = CLNode->getInOut( );
	NodePtr randomChild = CLNode->randomChild(); // Here we make sure that clustering and splitting are invertible
	if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO;
	}


	const NodePtr productionNode = CLNode-> getHistory( false , DSH()->hardScaleFactor() );

	weight = decideClustering(productionNode,CLNode,projected);

	if ( !CLNode->allAbove( mergePt() ) )return ZERO;

	if ( !productionNode->xcomb()->willPassCuts() )return ZERO;

	Energy startscale = CKKW_StartScale( productionNode );

	fillHistory( startscale , productionNode , CLNode );

	currentNode()->runningPt( fillProjector( projected ? 2 : 1 ) );

	weight = history.back().weightalphaReweight()*pdfReweight();

	if( weight == 0. )return ZERO;


	pair<CrossSection , CrossSection> DipAndPs =
	CLNode->calcDipandPS( startscale* currentME()->renFac() );
	// [ DEBUG ]
	if ( currentNode()->legsize() == 6&&Debug::level > 2 )
	debugReal("RBSI",weight,DipAndPs.second,DipAndPs.first);
	//Here we add the PS and subtrac the Dip ( only above the dynamical cutoff )
	return weightas( startscaleDSH()->renFac() )/SM().alphaS()*
	currentNode()->children().size()*( DipAndPs.second-DipAndPs.first );
	}


	CrossSection Merger::MergingDSigDRBornGamma( ){

	double weightCL = 0.;
	weight = 1.;

	if ( !currentNode()->children().empty() ) {
	auto const inOutPair = currentNode()->getInOut();
	// Here we make sure that clustering and splitting are invertible.
	NodePtr randomChild = currentNode()->randomChild();
	// Check if point is part of the ME region.
	if( !matrixElementRegion( inOutPair.first ,
	inOutPair.second ,
	randomChild->pT() ,
	theMergePt ) )weight *= 0.;
	}

	const NodePtr productionNode = currentNode()->getHistory( true , DSH()->hardScaleFactor() );
	NodePtr CLNode;
	NodePtr BornCL;


	if( !currentNode()->children().empty() ){
	if ( UseRandom::rndbool() ){
	CLNode = currentNode()->randomChild();
	bool inhist = CLNode->isInHistoryOf( productionNode );
	weight = inhist?( -2.int( currentNode()->children().size() ) ):0.;
	projected = true;
	weightCL = 2.*int( currentNode()->children().size() );
	BornCL = CLNode-> getHistory( false , DSH()->hardScaleFactor() );
	}else{
	weight = 2.;
	projected = false;
	}
	}else{
	weight = 1.;
	projected = false;
	}


	if ( treefactory()->onlymulti() != -1&&
	treefactory()->onlymulti() !=
	int( currentNode()->legsize()-(projected ? 1 : 0) ) )
	return ZERO;


	if( !currentNode()->allAbove( mergePt()*(1.-1e-6) ) )weight = 0.;
	if( CLNode&&!CLNode->allAbove( mergePt()*(1.-1e-6) ) )weightCL = 0.;
	if ( !productionNode->xcomb()->willPassCuts() ){
	return ZERO;
	}

	CrossSection res = ZERO;
	bool maxMulti = currentNode()->legsize() == int( maxLegsLO() );


	if( weight != 0. ){
	Energy startscale = CKKW_StartScale( productionNode );
	fillHistory( startscale , productionNode , currentNode() );
	currentNode()->runningPt( fillProjector( (projected ? 1 : 0) ) );
	weight = history.back().weightalphaReweight()*pdfReweight();
	if( weight == 0.&&weightCL == 0. )return ZERO;

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

	if( CLNode&&theGamma != 1. ){
	Energy startscale = CKKW_StartScale( BornCL );
	fillHistory( startscale , BornCL , CLNode );
	currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
	weightCL = history.back().weightalphaReweight()*pdfReweight();
	CrossSection diff = ZERO;
	currentME()->factory()->setAlphaParameter( 1. );
	diff -= weightCLCLNode->dipole()->dSigHatDR( sqr( startscale currentME()->renFac() ) );
	currentME()->factory()->setAlphaParameter( theGamma );

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

	diff += weightCLCLNode->dipole()->dSigHatDR( sqr( startscale currentME()->renFac() ) );
	currentME()->factory()->setAlphaParameter( 1. );

	res += diff;
	}


	return res;
	}





	/// calculate the history weighted born cross section
	#include "ThePEG/Handlers/SamplerBase.h"
	CrossSection Merger::MergingDSigDRBornCheapME( ){
	// Check if phase space poing is in ME region
	if ( !currentNode()->children().empty() ) {
	auto const & inOutPair = currentNode()->getInOut();
	NodePtr randomChild = currentNode()->randomChild(); // Here we make sure that clustering and splitting are invertible
	if( !matrixElementRegion( inOutPair.first , inOutPair.second , randomChild->pT() , theMergePt ) )return ZERO;
	}

	CrossSection meweight=(SamplerBase::runLevel() == SamplerBase::RunMode)?TreedSigDR( 60.GeV, 1. ):1.nanobarn;

	if (SamplerBase::runLevel() == SamplerBase::RunMode && theHighMeWeightMap.count(currentNode()) ) {
	if (abs(meweight)>abs(theHighMeWeightMap[currentNode()])) {
	theHighMeWeightMap[currentNode()]=abs(meweight);
	}else{
	if (meweight/theHighMeWeightMap[currentNode()]<UseRandom::rnd()) {
	return ZERO;
	}
	}
	}else{
	theHighMeWeightMap.insert({currentNode(),abs(meweight)});
	}



	// get the history for the process
	const NodePtr productionNode = currentNode()-> getHistory( true , DSH()->hardScaleFactor() );

	// decide if to cluster
	weight = decideClustering(productionNode,currentNode(),projected);

	if ( treefactory()->onlymulti() != -1&&
	treefactory()->onlymulti() != int( currentNode()->legsize()-( projected ? 1 : 0 ) ) )
	return ZERO;
	// Check for cuts on the production proces.
	if ( !productionNode->xcomb()->willPassCuts() )return ZERO;
	// calculate the staring scale
	Energy startscale = CKKW_StartScale( productionNode );

	// fill history with caluclation of sudakov supression
	fillHistory( startscale , productionNode , currentNode() );
	// fill the projector -> return the scale of the last splitting
	currentNode()->runningPt( fillProjector( projected ? 1 : 0 ) );
	// the weight has three components to get shower history weight
	weight = history.back().weightalphaReweight()*pdfReweight();
	if( weight == 0. )return ZERO;
	//calculate the cross section
	return weightTreedSigDR( startscale , 1. )theHighMeWeightMap[currentNode()]/meweight;
	}








	CrossSection Merger::TreedSigDR( Energy startscale , double diffAlpha ){

	currentME()->setScale( sqr( startscale ) , sqr( startscale ) );
	CrossSection res = currentME()->dSigHatDRB();
	- bool useDipolesForME=true;
	+ bool useDipolesForME=false;
	if (useDipolesForME && !currentNode()->children().empty()){
	res=ZERO;
	for (auto const & child : currentNode()->children() )
	res-=child->dipole()->dSigHatDR(sqr( startscale ));
	}
	if ( projected && emitDipoleMEDiff ) {
	CrossSection resDip=ZERO;
	for (auto const & child : currentNode()->children() )
	resDip-=child->dipole()->dSigHatDR(sqr( startscale ));
	setEmissionProbability(1.-min(resDip/res,res/resDip));
	}else{
	setEmissionProbability(0.);
	}

	if ( diffAlpha != 1. ) {
	res += currentME()->dSigHatDRAlphaDiff( diffAlpha );
	}
	if( std::isnan( double( res/nanobarn ) ) ){
	generator()->logWarning(Exception()
	<< "Merger: TreedSigDR is nan"
	<< Exception::warning);
	res = ZERO;};
	return res;
	}

	CrossSection Merger::LoopdSigDR( Energy startscale ){
	currentME()->setScale( sqr( startscale ) , sqr( startscale ) );
	currentME()->doOneLoopNoBorn();
	CrossSection res = currentME()->dSigHatDRV()+
	currentME()->dSigHatDRI();
	currentME()->noOneLoopNoBorn();
	return res;
	}

	Energy Merger::fillProjector( int pjs ){
	// in the shower handler the scale is multiplied
	// by DSH()->hardScaleFactor() so here we need
	// to devide by the factor.
	double xiQSh = history.begin()->node->legsize() == N0()?DSH()->hardScaleFactor():1.;
	if( pjs == 0 ){
	return ( history.size() == 1?1.:( 1./xiQSh ) )*history.back().scale;
	}
	for( auto const & hs : history )
	if ( isProjectorStage( hs.node , pjs )&&pjs != 0 ){
	currentNode()->xcomb()->lastProjector( hs.node->xcomb() );
	return ( hs.node == history[0].node?1.:( 1./xiQSh ) )*hs.scale;
	}

	throw Exception() << "Could not fill projector." << Exception::abortnow;
	return ZERO;
	}

	double Merger::pdfReweight(){ // TODO factorization scale inside

	double res = 1.;
	// consider both sides.
	for( int side : {0 , 1} ){
	// The side scale defines the scale that the leg is changig thou the history.
	// We start reweighting at the seed process.
	// We only need to calculate the pdf if the emission whould change the leg,
	// otherwise the leg remains the same.
	// If no emission is prduced from this leg the pdf ratios from this leg are always 1.
	Energy sidescale=history[0].scale*(
	history[0].node->legsize() == N0() ?
	currentME()->facFac():
	DSH()->facFac());
	bool sidechanged=false;
	// only if the incoming parton is coloured.
	if( history[0].node->xcomb()->mePartonData()[side]->coloured() ){
	// go though the history.
	for ( auto const & hs : history ){
	//pdf-ratio only to the last step
	if ( !hs.node->parent() ) continue;
	if ( hs.node == history.back().node ) continue;
	if ( !hs.node->dipole() ){
	throw Exception()
	<< "\nMerger: pdfReweight: history step has no dipol. "
	<< Exception::abortnow;
	return 0.;
	}
	// if the emitter is the side the emission changes the momentum fraction.
	if(!(hs.node->dipole()->bornEmitter()==side\|\|
	// II dipoles dont change the momentum fraction of the spectator only FI
	( hs.node->dipole()->bornSpectator()==side && hs.node->dipole()->bornEmitter()>1)))
	continue;
	const bool fromIsME = false;
	const bool toIsME = history[0].node == hs.node;
	res *= pdfratio( hs.node,
	//numerator
	sidescale,
	// denominator
	DSH()->facFac()*hs.node->pT(),
	side,
	fromIsME,
	toIsME);
	sidescale= DSH()->facFac()*hs.node->pT();
	sidechanged=true;
	}

	const bool fromIsME = true;
	const bool toIsME = !sidechanged && history[0].node->legsize() == N0();

	res *= pdfratio( history.back().node,
	sidescale,
	history[0].scale * currentME()->facFac(),
	side,
	fromIsME,
	toIsME );
	}
	}
	if ( std::isnan( res ) )
	generator()->logWarning(Exception()
	<< "Merger: pdfReweight is nan."
	<< Exception::warning);

	return res;
	}


	double Merger::cmwAlphaS(Energy q)const{
	double alphas=1.;
	using Constants::pi;
	// No cmw-scheme
	if (theCMWScheme==0) {
	alphas = as( q );
	}
	// Linear cmw-scheme
	else if(theCMWScheme==1){
	alphas = as( q );
	alphas =1.+(3.(67./18.-1./6.*sqr(pi))
	-5./9.Nf(q)) alphas/2./pi;
	}
	// cmw-scheme as factor in argument.
	else if(theCMWScheme==2){
	double kg=exp(-(67.-3.sqr(pi)-10/3Nf(q))
	/( 2. (33.-2.Nf(q))));
	//Note factor 2 since we here dealing with Energy
	alphas = as(max(kg*q,1_GeV));
	}else{
	throw Exception()
	<< "This CMW-Scheme is not implemented."
	<< Exception::abortnow;

	}
	return alphas;

	}




	double Merger::alphaReweight(){
	double res = 1.;

	Energy Q_R = ( history[0].node->legsize() == N0()?
	currentME()->renFac():
	DSH()->renFac() )*
	history[0].scale;

	using Constants::pi;
	const auto Q_qed=history[0].node->nodeME()->factory()->scaleChoice()->renormalizationScaleQED();
	const auto Oew=history[0].node->nodeME()->orderInAlphaEW();
	const auto Oqcd=history[0].node->nodeME()->orderInAlphaS();

	if (!history[0].node->children().empty()) {
	assert(Oqcd!=0);
	}

	res *= pow( SM().alphaEMME( Q_qed )/ SM().alphaEMMZ() , Oew );

	res *= pow( as(Q_R) / SM().alphaS() , Oqcd );





	for ( auto const & hs : history )
	if ( hs.node!= history.back().node ){
	Energy q_i = DSH()->renFac()* hs.node->pT();
	res *= cmwAlphaS(q_i) / SM().alphaS();
	}

	if ( std::isnan( res ) )
	generator()->logWarning(Exception()
	<< "Merger: alphaReweight is nan. "<< Exception::warning);
	return res;
	}

	void Merger::fillHistory( Energy scale , NodePtr begin , NodePtr endNode ){

	history.clear();
	double sudakov0_n = 1.;

	history.push_back( {begin , sudakov0_n , scale} );


	double xiQSh = history.begin()->node->legsize() == N0()?
	DSH()->hardScaleFactor():1.;

	scale *= xiQSh;
	if ( begin->parent()\|\|!isUnitarized ) {
	while ( begin->parent() && ( begin != endNode ) ) {
	if ( !dosudakov( begin , scale , begin->pT() , sudakov0_n ) ){
	history.push_back( { begin->parent() , 0. , scale } );
	}

	if ( std::isnan( sudakov0_n ) )
	generator()->logWarning(Exception()
	<< "Merger: sudakov"<<scale/GeV<<" "
	<<begin->pT()/GeV<<"0_n is nan. "
	<< Exception::warning);

	scale = begin->pT();
	history.push_back( { begin->parent() , sudakov0_n , begin->pT() } );
	begin = begin->parent();
	}

	Energy notunirunning = scale;

	if ( !isUnitarized&&N()+N0() > int( currentNode()->legsize() ) ) {
	if ( !dosudakov( begin , notunirunning , mergePt() , sudakov0_n ) ){
	history.back().weight = 0.;
	}else{
	history.back().weight = sudakov0_n;
	}
	}
	}
	if( history.size() == 1 )scale /= DSH()->hardScaleFactor();
	}




	double Merger::sumPdfReweightExpansion()const{
	double res = 0.;
	Energy beam1Scale = history[0].scale*
	( history[0].node->legsize() == N0()?
	currentME()->renFac():
	DSH()->facFac() );
	Energy beam2Scale = history[0].scale*
	( history[0].node->legsize() == N0()?
	currentME()->renFac():
	DSH()->facFac() );

	for ( auto const & hs : history ){
	//pdf expansion only to the last step
	if( !hs.node->parent() )continue;
	if( hs.node->xcomb()->mePartonData()[0]->coloured()&&
	hs.node->nodeME()->lastX1()>0.99 )return 0.;
	if( hs.node->xcomb()->mePartonData()[1]->coloured()&&
	hs.node->nodeME()->lastX2()>0.99 )return 0.;

	if( hs.node->nodeME()->lastX1()<0.00001 )return 0.;
	if( hs.node->nodeME()->lastX2()<0.00001 )return 0.;

	if ( hs.node->dipole()->bornEmitter() == 0 ){
	res += pdfExpansion( hs.node , 0 ,
	beam1Scale ,
	( hs.node->pT() ) ,
	hs.node->nodeME()->lastX1() ,
	Nf( history[0].scale ) ,
	history[0].scale );
	beam1Scale = ( hs.node->pT() )*DSH()->facFac();
	}
	else if ( hs.node->dipole()->bornEmitter() == 1 ){
	res += pdfExpansion( hs.node , 1 ,
	beam2Scale ,
	( hs.node->pT() ) ,
	hs.node->nodeME()->lastX2() ,
	Nf( history[0].scale ) ,
	history[0].scale );
	beam2Scale = ( hs.node->pT() )*DSH()->facFac();
	}
	// if we're here we know hs.node->dipole()->bornEmitter() > 1
	// works only in collinear scheme
	else if ( hs.node->dipole()->bornSpectator() == 0 ){
	res += pdfExpansion( hs.node , 0 ,
	beam1Scale ,
	( hs.node->pT() ) ,
	hs.node->nodeME()->lastX1() ,
	Nf( history[0].scale ) ,
	history[0].scale );
	beam1Scale = ( hs.node->pT() )*DSH()->facFac();
	}
	else if ( hs.node->dipole()->bornSpectator() == 1 ){
	res += pdfExpansion( hs.node , 1 ,
	beam2Scale ,
	( hs.node->pT() ) ,
	hs.node->nodeME()->lastX2() ,
	Nf( history[0].scale ) ,
	history[0].scale );
	beam2Scale = ( hs.node->pT() )*DSH()->facFac();
	}
	}

	if ( currentNode()->xcomb()->mePartonData()[0]->coloured() ){
	res += pdfExpansion( history.back().node , 0 ,
	beam1Scale ,
	history[0].scale* currentME()->renFac() ,
	( history.back() ).node->nodeME()->lastX1() ,
	Nf( history[0].scale ) ,
	history[0].scale );

	}
	if ( currentNode()->xcomb()->mePartonData()[1]->coloured() ) {
	res += pdfExpansion( history.back().node , 1 ,
	beam2Scale ,
	history[0].scale* currentME()->renFac() ,
	( history.back() ).node->nodeME()->lastX2() ,
	Nf( history[0].scale ) ,
	history[0].scale );
	}
	return res;
	}





	#include "Herwig/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
	double Merger::pdfExpansion( NodePtr node ,
	int side , Energy running ,
	Energy next , double x ,
	int nlp , Energy fixedScale ) const {

	tcPDPtr particle , parton;
	tcPDFPtr pdf;
	if ( side == 0 ) {
	particle = node->nodeME()->lastParticles().first->dataPtr();
	parton = node->nodeME()->lastPartons().first->dataPtr();
	pdf = node->xcomb()->partonBins().first->pdf();
	}else{
	assert( side == 1 );
	particle = node->nodeME()->lastParticles().second->dataPtr();
	parton = node->nodeME()->lastPartons().second->dataPtr();
	pdf = node->xcomb()->partonBins().second->pdf();
	}

	//copied from PKOperator
	double res = 0.;
	int number = 10;
	for ( int nr = 0;nr<number;nr++ ){
	double restmp = 0;

	using namespace RandomHelpers;
	double r = UseRandom::rnd();
	double eps = 1e-3;

	pair<double , double> zw =
	generate( ( piecewise() ,
	flat( 0.0 , x ) ,
	match( inverse( 0.0 , x , 1.0 ) +
	inverse( 1.0+eps , x , 1.0 ) ) ) , r );

	double z = zw.first;
	double mapz = zw.second;
	double PDFxparton = pdf->xfx( particle , parton , sqr( fixedScale ) , x )/x;
	double CA = SM().Nc();
	double CF = ( SM().Nc()SM().Nc()-1.0 )/( 2.SM().Nc() );


	if ( abs( parton->id() ) < 7 ) {

	double PDFxByzgluon = pdf->xfx( particle ,
	getParticleData( ParticleID::g ) ,
	sqr( fixedScale ) , x/z )*z/x;
	double PDFxByzparton = pdf->xfx( particle , parton ,
	sqr( fixedScale ) , x/z )*z/x;
	assert( abs( parton->id() ) < 7 );

	restmp += CF( 3./2.+2.log( 1.-x ) ) * PDFxparton;
	if ( z > x ) {
	restmp += 0.5 * ( sqr( z ) + sqr( 1.-z ) ) * PDFxByzgluon / z;
	restmp += 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. )*Nf( history[0].scale )
	+ 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::sumAlphaSReweightExpansion()const{
	double res = 0.;
	if ( !( history[0].node->children().empty() ) ){
	res += alphasExpansion( history[0].scale* DSH()->renFac() ,
	history[0].scale* currentME()->renFac() )*
	int( history[0].node->legsize()-N0() );
	}
	// dsig is calculated with fixed alpha_s
	for ( auto const & hs : history ){
	//expansion only to the last step
	if( !hs.node->parent() )continue;
	res += alphasExpansion( hs.node->pT()*DSH()->renFac() , history[0].scale );
	}
	return res;
	}

	double Merger::sumFillHistoryExpansion(){
	double res = 0.;
	double xiQSh = history[0].node->legsize() == N0()?DSH()->hardScaleFactor():1.;
	for ( auto const & hs : history ){
	if( !hs.node->parent() )continue;
	doHistExpansion( hs.node , ( hs.node == history[0].node?xiQSh:1. )*hs.scale ,
	hs.node->pT() , history[0].scale , res );
	}
	return res;
	}


	MergingFactoryPtr Merger::treefactory()const{return theTreeFactory;}


	void Merger::doinit(){
	if ( !DSH()->hardScaleIsMuF() ) {
	throw Exception()
	<< "Merger: Merging is currently only sensible "
	<< "if we are using the hardScale as MuF."
	<< Exception::abortnow;
	}
	}


	namespace{
	void setXCombScales(StdXCombPtr xc,Energy2 scale){
	xc->lastShowerScale ( scale );
	xc->lastCentralScale ( scale );
	xc->lastScale ( scale );
	xc->partonBinInstances().first->scale ( scale );
	xc->partonBinInstances().second->scale ( scale );
	xc->lastPartons().first->scale ( scale );
	xc->lastPartons().second->scale ( scale );
	}
	}

	CrossSection Merger::MergingDSigDR() {

	history.clear();
	assert(currentNode()==theFirstNodeMap[ currentME()]);

	if(DSH()->doesSplitHardProcess()){
	throw Exception()
	<< "Merger: The splithardprocess option is currently not supported."
	<< Exception::abortnow;
	}


	//get the PDF's (from ShowerHandler.cc)
	if (!DSH()->getPDFA()\|\|!DSH()->firstPDF().particle()){
	tSubProPtr sub = currentNode()->xcomb()->construct();
	const auto pb=currentNode()->xcomb()->partonBins();
	tcPDFPtr first = DSH()->getPDFA() ?
	tcPDFPtr(DSH()->getPDFA()) : DSH()->firstPDF().pdf();
	tcPDFPtr second = DSH()->getPDFB() ?
	tcPDFPtr(DSH()->getPDFB()) : DSH()->secondPDF().pdf();
	DSH()->resetPDFs( {first,second },pb );
	}



	DSH()->eventHandler( generator()->eventHandler() );

	CrossSection res = ZERO;
	if( currentNode()->subtractedReal() ){
	res = MergingDSigDRRealStandard();
	theCurrentMaxLegs = maxLegsNLO();
	}else if( currentNode()->virtualContribution() ){
	res = MergingDSigDRVirtualStandard();
	theCurrentMaxLegs = maxLegsNLO();
	}else if( theGamma!= 1. ){
	res = MergingDSigDRBornGamma();
	theCurrentMaxLegs = maxLegsLO();
	}else{
	res = MergingDSigDRBornStandard();
	theCurrentMaxLegs = maxLegsLO();
	}

	auto lxc= currentME()->lastXCombPtr();
	setXCombScales(lxc,sqr( currentNode()->runningPt()));


	auto lp= currentME()->lastXCombPtr()->lastProjector();

	if( lp )
	setXCombScales( lp, sqr( currentNode()->runningPt()));

	if ( res == ZERO ){
	history.clear();
	return ZERO;
	}


	cleanup( currentNode() );
	lxc->subProcess( SubProPtr() );

	history.clear();

	if( !std::isfinite( double( res/nanobarn ) ) ){
	generator()->logWarning(Exception()
	<< "Merger weight is " << res/nanobarn<< " -> setting to 0"
	<< Exception::warning);
	return ZERO;
	}

	return res;

	}




	#include "Herwig/PDF/HwRemDecayer.h"
	void Merger::CKKW_PrepareSudakov( NodePtr node , Energy running ){
	//cleanup( node );
	node->xcomb()->lastPartons().first->scale(sqr(running));
	node->xcomb()->lastPartons().second->scale(sqr(running));
	node->xcomb()->partonBinInstances().first->scale(sqr(running));
	node->xcomb()->partonBinInstances().second->scale(sqr(running));
	tSubProPtr sub = node->xcomb()->construct();
	const auto pb=node->xcomb()->partonBins();

	DSH()->setCurrentHandler();

	DSH()->currentHandler()->generator()->currentEventHandler(
	currentNode()->xcomb()->eventHandlerPtr() );

	DSH()->currentHandler()->remnantDecayer()->setHadronContent( currentNode()->xcomb()->lastParticles() );
	DSH()->eventRecord().clear();
	DSH()->eventRecord().slimprepare( sub , dynamic_ptr_cast<tStdXCombPtr>( node->xcomb() ) , DSH()->pdfs() , currentNode()->xcomb()->lastParticles() );
	DSH()->hardScales( sqr( running ) );
	}


	Energy Merger::CKKW_StartScale( NodePtr node ) const {
	Energy res = generator()->maximumCMEnergy();
	const int N=node->legsize();
	const auto & data=node->nodeME()->mePartonData();
	const auto & momenta=node->nodeME()->lastMEMomenta();
	if( !node->children().empty() ){
	for ( int i = 0; i<N ; i++ ){ if ( !data[i]->coloured() )continue;
	for ( int j = 2; j<N ; j++ ){ if ( i == j\|\|!data[j]->coloured() )continue;
	for ( int k = 0; k<N ; k++ ){ if ( !data[k]->coloured() \|\| i == k \|\| j == k )continue;
	if(i<2){
	if(k<2) {res = min( res , IILTK->lastPt( momenta[i] , momenta[j] , momenta[k] ));}
	else {res = min( res , IFLTK->lastPt( momenta[i] , momenta[j] , momenta[k] ));}
	}else{
	if(k<2) {res = min( res , FILTK->lastPt( momenta[i] , momenta[j] , momenta[k] ));}
	else {res = min( res , FFLTK->lastPt( momenta[i] , momenta[j] , momenta[k] ));}
	}
	}
	}
	}
	}else{
	node->nodeME()->factory()->scaleChoice()->setXComb( node->xcomb() );
	res = sqrt( node->nodeME()->factory()->scaleChoice()->renormalizationScale() );
	}
	node->nodeME()->factory()->scaleChoice()->setXComb( node->xcomb() );
	res = max( res , sqrt( node->nodeME()->factory()->scaleChoice()->renormalizationScale() ) );
	return res;
	}




	double Merger::alphasExpansion( Energy next , Energy fixedScale ) const {
	double betaZero = ( 11./6. )SM().Nc() - ( 1./3. )Nf( history[0].scale );
	// TODO factorize the return statement, turn ?: into if
	return ( betaZero*log( sqr( fixedScale/next ) ) )+
	( theCMWScheme>0?( 3.( 67./18.-1./6.Constants::pi*Constants::pi )
	-5./9.*Nf( history[0].scale ) ):0. );
	}


	double Merger::pdfratio( NodePtr node ,
	Energy numerator_scale ,
	Energy denominator_scale ,
	int side ,
	bool fromIsME,
	bool toIsME ){
	StdXCombPtr bXc = node->xcomb();
	if( !bXc->mePartonData()[side]->coloured() )
	throw Exception()
	<< "Merger: pdf-ratio required for non-coloured particle."
	<< Exception::abortnow;


	double from = 1.;
	double to = 1.;
	if ( side == 0 ){
	if ( denominator_scale == numerator_scale && fromIsME==toIsME ) {
	return 1.;
	}


	if (fromIsME) {
	from = node->nodeME()->pdf1( sqr( denominator_scale ) );
	}else{
	from = DSH()->firstPDF().xfx(node->xcomb()->lastPartons().first->dataPtr(),
	sqr( denominator_scale ),
	node->xcomb()->lastX1())/node->xcomb()->lastX1();
	}


	if (toIsME) {
	to = node->nodeME()->pdf1( sqr( numerator_scale ) );
	}else{
	to = DSH()->firstPDF().xfx(node->xcomb()->lastPartons().first->dataPtr(),
	sqr( numerator_scale ),
	node->xcomb()->lastX1())/node->xcomb()->lastX1();
	}

	if ( ( to < 1e-8\|\|from < 1e-8 )&&( to/from>10000000. ) ){
	generator()->logWarning(Exception()
	<< "Merger: pdfratio to = " << to << " from = " << from
	<< Exception::warning);
	return 0.;
	}
	}
	else{
	if ( denominator_scale == numerator_scale && fromIsME==toIsME ) {
	return 1.;
	}

	if (fromIsME) {
	from = node->nodeME()->pdf2( sqr( denominator_scale ) );
	}else{
	from =DSH()->secondPDF().xfx(node->xcomb()->lastPartons().second->dataPtr(),
	sqr( denominator_scale ),
	node->xcomb()->lastX2())/node->xcomb()->lastX2();
	}


	if (toIsME) {
	to = node->nodeME()->pdf2( sqr( numerator_scale ) );
	}else{
	to = DSH()->secondPDF().xfx(node->xcomb()->lastPartons().second->dataPtr(),
	sqr( numerator_scale ),
	node->xcomb()->lastX2())/node->xcomb()->lastX2();
	}

	if ( ( to < 1e-8\|\|from < 1e-8 )&&( to/from>10000000. ) ){
	generator()->logWarning(Exception()
	<< "Merger: pdfratio to = "
	<< to << " from = " << from
	<< Exception::warning);
	return 0.;}
	}
	return to/from;
	}



	bool Merger::dosudakov( NodePtr node , Energy running , Energy next , double& sudakov0_n ) {
	CKKW_PrepareSudakov( node , running );
	for( DipoleChain const & chain : DSH()->eventRecord().chains() ){
	for( Dipole const & dip : chain.dipoles() ){
	sudakov0_n *= singlesudakov( dip , next , running , { true , false } );
	sudakov0_n *= singlesudakov( dip , next , running , { false , true } );
	if ( sudakov0_n == 0.0 ){
	cleanup( node );
	return false;
	}
	}
	}
	cleanup( node );
	return true;
	}

	bool Merger::doHistExpansion( NodePtr node , Energy running , Energy next , Energy fixedScale , double& histExpansion ) {
	CKKW_PrepareSudakov( node , running );
	for( DipoleChain const & chain : DSH()->eventRecord().chains() ){
	for( Dipole const & dip : chain.dipoles() ){
	histExpansion += singleHistExpansion( dip , next , running , fixedScale , { true , false } );;
	histExpansion += singleHistExpansion( dip , next , running , fixedScale , { false , true } );
	}
	}
	cleanup( node );
	return true;
	}



	bool Merger::isProjectorStage( NodePtr node , int pjs )const{
	return ( pjs == int( ( currentNode()->legsize() - node->legsize() ) ) );
	}

	void Merger::cleanup( NodePtr node ) {
	DSH()->eventRecord().clear();
	if( !node->xcomb()->subProcess() )return;
	ParticleVector vecfirst = node->xcomb()->subProcess()->incoming().first->children();
	for( auto const & particle : vecfirst )
	node->xcomb()->subProcess()->incoming().first->abandonChild( particle );

	ParticleVector vecsecond = node->xcomb()->subProcess()->incoming().second->children();
	for( auto const & particle : vecsecond )
	node->xcomb()->subProcess()->incoming().second->abandonChild( particle );
	node->xcomb()->subProcess( SubProPtr() );
	}

	double Merger::singlesudakov( Dipole dip , Energy next , Energy running , pair<bool , bool> conf ){

	double res = 1.;
	tPPtr emitter = dip.emitter( conf );
	tPPtr spectator = dip.spectator( conf );
	DipoleSplittingInfo candidate( dip.index( conf ) , conf ,
	dip.emitterX( conf ) ,
	dip.spectatorX( conf ) ,
	emitter , spectator );


	if ( DSH()->generators().find( candidate.index() ) == DSH()->generators().end() ) DSH()->getGenerators( candidate.index() );

	auto const & gens = DSH()->generators().equal_range( candidate.index() );

	for ( auto gen = gens.first; gen != gens.second; ++gen ) {
	if ( !( gen->first == candidate.index() ) )
	continue;

	Energy dScale = gen->second->splittingKinematics()->dipoleScale( emitter->momentum() , spectator->momentum() );
	candidate.scale( dScale );
	candidate.continuesEvolving();
	Energy ptMax = gen->second->splittingKinematics()->ptMax( candidate.scale() , candidate.emitterX() , candidate.spectatorX() ,
	candidate.index() , *gen->second->splittingKernel() );

	candidate.hardPt( min( running , ptMax ) );

	if ( candidate.hardPt()>next ){
	res *= gen->second->sudakov( candidate , next );
	}
	}

	return res;
	}


	double Merger::singleHistExpansion( Dipole dip , Energy next , Energy running , Energy fixedScale , pair<bool , bool> conf ){

	double res = 0.;
	tPPtr emitter = dip.emitter( conf );
	tPPtr spectator = dip.spectator( conf );
	DipoleSplittingInfo candidate( dip.index( conf ) ,
	conf , dip.emitterX( conf ) ,
	dip.spectatorX( conf ) ,
	emitter , spectator );

	if ( DSH()->generators().find( candidate.index() ) ==
	DSH()->generators().end() )
	DSH()->getGenerators( candidate.index() );

	auto const & gens = DSH()->generators().equal_range( candidate.index() );

	for ( auto gen = gens.first; gen != gens.second; ++gen ) {
	if ( !( gen->first == candidate.index() ) )
	continue;
	Energy dScale = gen->second->splittingKinematics()->dipoleScale(
	emitter->momentum() , spectator->momentum() );
	candidate.scale( dScale );
	candidate.continuesEvolving();
	Energy ptMax = gen->second->
	splittingKinematics()->ptMax(
	candidate.scale() , candidate.emitterX() ,
	candidate.spectatorX() , candidate.index() ,
	*gen->second->splittingKernel() );

	candidate.hardPt( min( running , ptMax ) );
	if ( candidate.hardPt()>next ){
	res += gen->second->sudakovExpansion( candidate , next , fixedScale );
	}
	}

	return res;
	}




	void Merger::firstNodeMap( MatchboxMEBasePtr a , NodePtr b ){
	theFirstNodeMap.insert( { a , b } );
	}



	map<MatchboxMEBasePtr , NodePtr> Merger::firstNodeMap()const{return theFirstNodeMap;}




	void Merger::setXComb( tStdXCombPtr xc ){
	currentNode()->setXComb( xc );
	}
	void Merger::setKinematics( ){
	currentNode()->setKinematics();
	}
	void Merger::clearKinematics( ){
	currentNode()->clearKinematics();
	}

	void Merger::flushCaches(){
	if (currentNode()&&currentNode()->xcomb()->lastParticles().first) {
	currentNode()->flushCaches();
	}
	}

	bool Merger::generateKinematics( const double * r ){
	return currentNode()->firstgenerateKinematics( r , ! currentNode()->subtractedReal() );
	}

	bool Merger::matrixElementRegion( PVector incoming , PVector outgoing , Energy winnerScale , Energy cutscale )const{


	Energy ptx = Constants::MaxEnergy;
	bool foundwinnerpt = false;
	using namespace boost;
	//FF

	for( auto const & em : outgoing ){ if ( ! em->coloured() ) continue;
	for( auto const & emm : outgoing ){ if ( !emm->coloured() ) continue; if ( em == emm ) continue;
	for( auto const & spe : outgoing ){ if ( !spe->coloured() ) continue; if ( em == spe\|\|emm == spe ) continue;

	if ( !( em->id() == -emm->id()\|\|emm->id()>6 ) )continue;

	Energy pt = ZERO;
	if ( em->momentum().m()<= 10_MeV &&
	emm->momentum().m()<= 10_MeV &&
	spe->momentum().m()<= 10_MeV ) {
	pt = FFLTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
	}else{
	pt = FFMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
	}

	if ( abs( pt-winnerScale ) < 10_MeV ) {
	foundwinnerpt = true;
	}
	ptx = min( ptx , pt );
	}
	}
	}

	//FI
	for( auto const & spe : incoming ){ if ( ! spe->coloured() ) continue;
	for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue;
	for( auto const & em : outgoing ){ if ( ! em->coloured() ) continue; if ( em == emm ) continue;

	if ( !( em->id() == -emm->id() \|\| emm->id()>6 ) )continue;


	Energy pt = ZERO;
	if ( em->momentum().m()<= 10_MeV &&
	emm->momentum().m()<= 10_MeV &&
	spe->momentum().m()<= 10_MeV ) {
	pt = FILTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
	}else{
	pt = FIMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
	}


	if ( abs( pt-winnerScale )<10_MeV ) {
	foundwinnerpt = true;
	}

	if( pt > ZERO )
	ptx = min( ptx , pt );
	}
	}
	}

	//IF
	for( auto const & em : incoming ){ if ( ! em->coloured() ) continue;
	for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue;
	for( auto const & spe : outgoing ){ if ( ! spe->coloured() ) continue; if ( emm == spe ) continue;

	if ( !( em->id()>6\|\| em->id() == emm->id() \|\|emm->id()>6 ) )continue;

	Energy pt = ZERO;

	if ( em->momentum().m()<= 10_MeV &&
	emm->momentum().m()<= 10_MeV &&
	spe->momentum().m()<= 10_MeV ) {
	//massless
	pt = IFLTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
	}else{
	//massiv
	pt = IFMTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );
	}


	if ( abs( pt-winnerScale )< 10_MeV ) {
	foundwinnerpt = true;
	}
	ptx = min( ptx , pt );
	}
	}
	}

	//II
	for( auto const & em : incoming ){ if ( ! em->coloured() ) continue;
	for( auto const & spe : incoming ){ if ( ! spe->coloured() ) continue; if ( em == spe ) continue;
	for( auto const & emm : outgoing ){ if ( ! emm->coloured() ) continue;
	if ( !( em->id()>6\|\|em->id() == emm->id() \|\|emm->id()>6 ) )continue;

	Energy pt = IILTK->lastPt( em->momentum() , emm->momentum() , spe->momentum() );

	if ( abs( pt-winnerScale )< 10_MeV ) {
	foundwinnerpt = true;
	}
	ptx = min( ptx , pt );
	}
	}
	}

	if( !foundwinnerpt ){
	generator()->logWarning( Exception()
	<< "Merger: Could not find winner with pt."
	<< "Run with -d3 to get phase space points. "
	<< Exception::warning );

	if( Debug::level > 2 ) {
	generator()->log() << "\nWarning: could not find winner with pt: " << winnerScale/GeV;
	for( auto const & emm : incoming ){
	generator()->log() << "\n" << emm->id() << " " << emm->momentum()/GeV << " " << emm->momentum().m()/GeV << flush;
	}
	for( auto const & emm : outgoing ){
	generator()->log() <<"\n" << emm->id() << " " << emm->momentum()/GeV << " " << emm->momentum().m()/GeV << flush;
	}
	}

	}

	return ( ptx>cutscale ) ;

	}

	bool Merger::notOnlyMulti() const {
	return ( treefactory()->onlymulti() != -1&&
	treefactory()->onlymulti() !=
	int( currentNode()->legsize()-( projected ? 1 : 0 ) ) );
	}

	int Merger::M()const{return theTreeFactory->M();}

	int Merger::N()const{return theTreeFactory->N();}

	void Merger::debugVirt(double weight,double w1,double w2,double w3,
	CrossSection matrixElement,double ww1,double ww2,
	double ww3, NodePtr node,CrossSection bornWeight)const{
	Energy minPT = Constants::MaxEnergy;
	for( auto const & no :currentNode()->children() )minPT = min( no->pT() , minPT );

	generator()->log() << "\nVIRT " << minPT/GeV << " " << weight << " " << w1;
	generator()->log() << " " << w2;
	generator()->log() << " " << w3;
	generator()->log() << " " << ( matrixElement/nanobarn ) << " " << ww1 << " " << ww2 << " " << ww3 << " " << node->pT()/GeV << " " << node->nodeME()->mePartonData()[3]->mass()/GeV << " " << ( bornWeightSM().alphaS()/( 2.ThePEG::Constants::pi )/nanobarn );
	}



	void Merger::debugReal(string realstr, double weight,
	CrossSection me, CrossSection dip)const {
	Energy minPT = Constants::MaxEnergy;
	for( auto const & no :currentNode()->children() )minPT = min( no->pT() , minPT );

	generator()->log() << "\n"<<realstr <<" "<< minPT/GeV << " " << weight << " " << ( me-dip )/nanobarn << " " << me/nanobarn << " " << dip/nanobarn;

	}

	// If needed , insert default implementations of virtual function defined
	// in the InterfacedBase class here ( using ThePEG-interfaced-impl in Emacs ).

	#include "ThePEG/Persistency/PersistentOStream.h"
	void Merger::persistentOutput( PersistentOStream & os ) const {
	os << theShowerExpansionWeights << theCMWScheme << projected <<
	isUnitarized << isNLOUnitarized <<
	theOpenInitialStateZ << theChooseHistory <<
	theN0 << theOnlyN << weight <<
	weightCB << theGamma << theSmearing << ounit( theIRSafePT , GeV ) <<
	ounit( theMergePt , GeV ) << ounit( theCentralMergePt , GeV )
	<< theLargeNBasis << FFLTK << FILTK <<
	IFLTK << IILTK << FFMTK << FIMTK << IFMTK <<
	theDipoleShowerHandler << theTreeFactory <<
	theFirstNodeMap<<theRealSubtractionRatio << emitDipoleMEDiff;

	}
	#include "ThePEG/Persistency/PersistentIStream.h"
	void Merger::persistentInput( PersistentIStream & is , int ) {
	is >> theShowerExpansionWeights >> theCMWScheme >> projected >>
	isUnitarized >> isNLOUnitarized >>
	theOpenInitialStateZ >>
	theChooseHistory >> theN0 >> theOnlyN >>
	weight >> weightCB >>
	theGamma >> theSmearing >> iunit( theIRSafePT , GeV ) >>
	iunit( theMergePt , GeV ) >>
	iunit( theCentralMergePt , GeV ) >> theLargeNBasis >>
	FFLTK >> FILTK >> IFLTK >>
	IILTK >> FFMTK >> FIMTK >>
	IFMTK >> theDipoleShowerHandler >> theTreeFactory >>
	theFirstNodeMap >> theRealSubtractionRatio >> emitDipoleMEDiff;
	}

	// * 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 ).
	#include "ThePEG/Utilities/DescribeClass.h"
	DescribeClass<Merger , Herwig::MergerBase>
	describeHerwigMerger( "Herwig::Merger" , "HwDipoleShower.so" );

	//ClassDescription<Merger> Merger::initMerger;
	// Definition of the static class description member.

	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Switch.h"

	void Merger::Init() {

	static ClassDocumentation<Merger> documentation
	( "The Merger class takes care of merging multiple LO & NLO cross sections." );


	static Reference<Merger , DipoleShowerHandler> interfaceShowerHandler
	( "DipoleShowerHandler" ,
	"Set the pointer to the Dipole Shower Handler" ,
	&Merger::theDipoleShowerHandler ,
	false , false ,
	true , true , false );

	static Switch<Merger , bool>
	interfaceShowerExpansionWeights
	( "ShowerExpansionWeights" ,
	"Calculate the expansions of the shower history to be NLO accurate" ,
	&Merger::theShowerExpansionWeights ,
	true , false , false );

	static SwitchOption
	interfaceShowerExpansionWeightsYes
	( interfaceShowerExpansionWeights ,
	"Yes" ,
	"Calculate shower expansions" ,
	true );
	static SwitchOption
	interfaceShowerExpansionWeightsNo
	( interfaceShowerExpansionWeights ,
	"No" ,
	"Do not calculate expansions" ,
	false );

	static Switch<Merger , unsigned int>
	interfacetheCMWScheme
	( "CMWScheme" ,
	"Use CMW-Scheme to calculate the alpha_s for the shower expressions." ,
	&Merger::theCMWScheme ,
	0 ,
	false , false );
	static SwitchOption interfacetheCMWSchemeOff
	(interfacetheCMWScheme,
	"Off",
	"No CMW-Scheme",
	0);
	static SwitchOption interfacetheCMWSchemeLinear
	(interfacetheCMWScheme,
	"Linear",
	"Linear CMW multiplication: alpha_s(q) -> alpha_s(q)(1+K_g*alpha_s(q)/2pi )",
	1);
	static SwitchOption interfacetheCMWSchemeFactor
	(interfacetheCMWScheme,
	"Factor",
	"Use factor in alpha_s argument: alpha_s(q) -> alpha_s(facq) with fac=exp(-(67-3pi^2-10/3Nf)/(33-2Nf)) ",
	2);

	static Parameter<Merger , Energy> interfaceMergerScale
	( "MergingScale" ,
	"The MergingScale." ,
	&Merger::theCentralMergePt ,
	GeV , 20.0GeV , 0.0GeV , 0*GeV ,
	false , false , Interface::lowerlim );

	static Reference<Merger , MergingFactory> interfaceMergerHelper
	( "MergingFactory" ,
	"Set a pointer to the MergingFactory" ,
	&Merger::theTreeFactory ,
	false , false ,
	true , true , false );


	static Parameter<Merger , double> interfacegamma
	( "gamma" ,
	"gamma parameter." ,
	&Merger::theGamma , 1.0 , 0.0 , 0 ,
	false , false , Interface::lowerlim );


	static Reference<Merger , ColourBasis> interfaceLargeNBasis
	( "LargeNBasis" ,
	"Set the large-N colour basis implementation." ,
	&Merger::theLargeNBasis ,
	false , false ,
	true , true , 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 Switch<Merger , bool>
	interfaceemitDipoleMEDiff
	( "emitDipoleMEDiff" , "Allow emissions of the unitarisation contribution with prob. 1-min(B_n/sum Dip_n,sum Dip_n/B_n)" ,
	&Merger::emitDipoleMEDiff , false , false , false );
	static SwitchOption interfaceemitDipoleMEDiffYes
	( interfaceemitDipoleMEDiff , "Yes" , "" , true );
	static SwitchOption interfaceemitDipoleMEDiffNo
	( interfaceemitDipoleMEDiff , "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" ,
	"Various ways of choosing the history weights: 0(default):dipole xs, 1:dipole/born, 2:flat, 3: 1/pt_dip" ,
	&Merger::theChooseHistory ,
	3 , -1 , 0 ,
	false , false , Interface::lowerlim );


	static Parameter<Merger , double>
	interfacetheRealSubtractionRatio(
	"RealSubtractionRatio" ,
	"Set the percentage the merging pt should be smeared." ,
	&Merger::theRealSubtractionRatio , 0. , 0. ,
	1. , 10.0 ,
	true , false , Interface::limited );

	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 5:36 PM (12 h, 39 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023662
Default Alt Text: (57 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions