No OneTemporary
Actions

Size

60 KB

Subscribers

None

View Options

	diff --git a/UnderlyingEvent/MPIHandler.cc b/UnderlyingEvent/MPIHandler.cc
	--- a/UnderlyingEvent/MPIHandler.cc
	+++ b/UnderlyingEvent/MPIHandler.cc
	@@ -1,828 +1,849 @@
	// -- C++ --
	//
	// MPIHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 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 MPIHandler class.
	//

	#include "MPIHandler.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/Handlers/SubProcessHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Cuts/JetCuts.h"
	#include "ThePEG/Cuts/SimpleKTCut.h"
	#include "ThePEG/PDF/PartonExtractor.h"

	#include "gsl/gsl_sf_bessel.h"

	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace Herwig;

	MPIHandler * MPIHandler::currentHandler_ = 0;

	bool MPIHandler::beamOK() const {
	return (HadronMatcher::Check(*eventHandler()->incoming().first) &&
	HadronMatcher::Check(*eventHandler()->incoming().second) );
	}


	tStdXCombPtr MPIHandler::generate(unsigned int sel) {
	//generate a certain process
	if(sel+1 > processHandlers().size())
	throw Exception() << "MPIHandler::generate called with argument out of range"
	<< Exception::runerror;

	return processHandlers()[sel]->generate();
	}

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

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

	void MPIHandler::finalize() {
	if( beamOK() ){
	statistics();
	}
	}

	void MPIHandler::initialize() {
	currentHandler_ = this;
	useMe();
	theHandler = generator()->currentEventHandler();
	//stop if the EventHandler is not present:
	assert(theHandler);

	//check if MPI is wanted
	if( !beamOK() ){
	throw Exception() << "You have requested multiple parton-parton scattering,\n"
	<< "but the model is not forseen for the beam setup you chose.\n"
	<< "You should therefore disable that by setting XXXGenerator:EventHandler:"
	<< "CascadeHandler:MPIHandler to NULL"
	<< Exception::runerror;
	}

	numSubProcs_ = subProcesses().size();

	if( numSubProcs_ != cuts().size() )
	throw Exception() << "MPIHandler::each SubProcess needs a Cuts Object"
	<< "ReferenceVectors are not equal in size"
	<< Exception::runerror;

	if( additionalMultiplicities_.size()+1 != numSubProcs_ )
	throw Exception() << "MPIHandler: each additional SubProcess needs "
	<< "a multiplicity assigned. This can be done in with "
	<< "insert MPIHandler:additionalMultiplicities 0 1"
	<< Exception::runerror;

	//identicalToUE_ = 0 hard process is identical to ue, -1 no one
	if( identicalToUE_ > (int)numSubProcs_ \|\| identicalToUE_ < -1 )
	throw Exception() << "MPIHandler:identicalToUE has disallowed value"
	<< Exception::runerror;

	// override the cuts for the additional scatters if energyExtrapolation_ is
	// set
	if (energyExtrapolation_ != 0 ) {
	overrideUECuts();
	}

	tcPDPtr gluon=getParticleData(ParticleID::g);
	//determine ptmin
	Ptmin_ = cuts()[0]->minKT(gluon);

	if(identicalToUE_ == -1){
	algorithm_ = 2;
	}else{
	if(identicalToUE_ == 0){
	//Need to work a bit, in case of LesHouches events for QCD2to2
	Ptr<StandardEventHandler>::pointer eH =
	dynamic_ptr_cast<Ptr<StandardEventHandler>::pointer>(eventHandler());

	if( eH ) {
	PtOfQCDProc_ = eH->cuts()->minKT(gluon);
	// find the jet cut in the new style cuts
	for(unsigned int ix=0;ix<eH->cuts()->multiCuts().size();++ix) {
	Ptr<JetCuts>::pointer jetCuts =
	dynamic_ptr_cast<Ptr<JetCuts>::pointer>(eH->cuts()->multiCuts()[ix]);
	if(jetCuts) {
	Energy ptMin=1e30*GeV;
	for(unsigned int iy=0;iy<jetCuts->jetRegions().size();++iy) {
	ptMin = min(ptMin,jetCuts->jetRegions()[iy]->ptMin());
	}
	if(ptMin<1e29*GeV&&ptMin>PtOfQCDProc_) PtOfQCDProc_ = ptMin;
	}
	}
	}
	else {
	if(PtOfQCDProc_ == -1.0*GeV)
	throw Exception() << "MPIHandler: You need to specify the pt cutoff "
	<< "used to in the LesHouches file for QCD2to2 events"
	<< Exception::runerror;
	}
	}
	else {
	PtOfQCDProc_ = cuts()[identicalToUE_]->minKT(gluon);
	}

	if(PtOfQCDProc_ > 2*Ptmin_)
	algorithm_ = 1;
	else
	algorithm_ = 0;

	if(PtOfQCDProc_ == ZERO)//pure MinBias mode
	algorithm_ = -1;
	}

	//Init all subprocesses
	for(unsigned int i=0; i<numSubProcs_; i++){
	theProcessHandlers.push_back(new_ptr(ProcessHandler()));
	processHandlers().back()->initialize(subProcesses()[i],
	cuts()[i], eventHandler());
	processHandlers().back()->initrun();
	}

	//now calculate the individual Probabilities
	XSVector UEXSecs;
	UEXSecs.push_back(processHandlers()[0]->integratedXSec());
	//save the hard cross section
	hardXSec_ = UEXSecs.front();

	//determine sigma_soft and beta
	if(softInt_){//check that soft ints are requested
	GSLBisection rootFinder;

	if(twoComp_){
	//two component model
	/*
	GSLMultiRoot eqSolver;
	slopeAndTotalXSec eq(this);
	pair<CrossSection, Energy2> res = eqSolver.value(eq, 10millibarn, 0.6GeV2);
	softXSec_ = res.first;
	softMu2_ = res.second;
	*/
	slopeBisection fs(this);
	try{
	softMu2_ = rootFinder.value(fs, 0.3GeV2, 1.GeV2);
	softXSec_ = fs.softXSec();
	}catch(GSLBisection::IntervalError){
	try{
	// Very low energies (e.g. 200 GeV) need this.
	// Very high energies (e.g. 100 TeV) produce issues with
	// this choice. This is a temp. fix.
	softMu2_ = rootFinder.value(fs, 0.25GeV2, 1.3GeV2);
	softXSec_ = fs.softXSec();
	}catch(GSLBisection::IntervalError){
	throw Exception() <<
	"\n*********************************************************\n" " Inconsistent MPI parameter choice for this beam energy *\n"
	"**********************************************************\n"
	"MPIHandler parameter choice is unable to reproduce\n"
	"the total cross section. Please check arXiv:0806.2949\n"
	"for the allowed parameter space."
	<< Exception::runerror;
	}
	}

	}else{
	//single component model
	TotalXSecBisection fn(this);
	try{
	softXSec_ = rootFinder.value(fn, 0millibarn, 5000millibarn);
	}catch(GSLBisection::IntervalError){
	throw Exception() <<
	"\n**********************************************************\n"
	"* Inconsistent MPI parameter choice for this beam energy *\n"
	"**********************************************************\n"
	"MPIHandler parameter choice is unable to reproduce\n"
	"the total cross section. Please check arXiv:0806.2949\n"
	"for the allowed parameter space."
	<< Exception::runerror;
	}
	}

	//now get the differential cross section at ptmin
	ProHdlPtr qcd = new_ptr(ProcessHandler());
	Energy eps = 0.1*GeV;
	Energy ptminPlus = Ptmin_ + eps;

	Ptr<SimpleKTCut>::pointer ktCut = new_ptr(SimpleKTCut(ptminPlus));
	ktCut->init();
	ktCut->initrun();

	CutsPtr qcdCut = new_ptr(Cuts(2*ptminPlus));
	qcdCut->add(dynamic_ptr_cast<tOneCutPtr>(ktCut));
	qcdCut->init();
	qcdCut->initrun();

	qcd->initialize(subProcesses()[0], qcdCut, eventHandler());
	qcd->initrun();

	// ds/dp_T^2 = 1/2/p_T ds/dp_T
	DiffXSec hardPlus = (hardXSec_-qcd->integratedXSec())/(2Ptmin_eps);

	betaBisection fn2(softXSec_, hardPlus, Ptmin_);
	try{
	beta_ = rootFinder.value(fn2, -10/GeV2, 2/GeV2);
	}catch(GSLBisection::IntervalError){
	try{
	// Very low energies (e.g. 200 GeV) need this.
	// Very high energies (e.g. 100 TeV) produce issues with
	// this choice. This is a temp. fix.
	beta_ = rootFinder.value(fn2, -5/GeV2, 8./GeV2);
	}catch(GSLBisection::IntervalError){
	throw Exception() << "MPIHandler: slope of soft pt spectrum couldn't be "
	<< "determined."
	<< Exception::runerror;
	}
	}
	}

	Probs(UEXSecs);
	//MultDistribution("probs.test");

	UEXSecs.clear();
	}

	void MPIHandler::MultDistribution(string filename) const {
	ofstream file;
	double p(0.0), pold(0.0);
	file.open(filename.c_str());
	//theMultiplicities
	Selector<MPair>::const_iterator it = theMultiplicities.begin();

	while(it != theMultiplicities.end()){
	p = it->first;
	file << it->second.first << " " << it->second.second
	<< " " << p-pold << '\n';
	it++;
	pold = p;
	}

	file << "sum of all probabilities: " << theMultiplicities.sum()
	<< endl;

	file.close();
	}

	void MPIHandler::statistics() const {
	ostream & file = generator()->misc();

	string line = "======================================="
	"=======================================\n";

	for(unsigned int i=0; i<cuts().size(); i++){
	Stat tot;
	if(i == 0)
	file << "Statistics for the UE process: \n";
	else
	file << "Statistics for additional hard Process " << i << ": \n";

	processHandlers()[i]->statistics(file, tot);
	file << "\n";
	}

	if(softInt_){
	file << line
	<< "Eikonalized and soft cross sections:\n\n"
	<< "Model parameters: "
	<< "ptmin: " << Ptmin_/GeV << " GeV"
	<< ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n"
	<< " "
	<< "DL mode: " << DLmode_
	<< ", CMenergy: " << generator()->maximumCMEnergy()/GeV
	<< " GeV" << '\n'
	<< "hard inclusive cross section (mb): "
	<< hardXSec_/millibarn << '\n'
	<< "soft inclusive cross section (mb): "
	<< softXSec_/millibarn << '\n'
	<< "total cross section (mb): "
	<< totalXSecExp()/millibarn << '\n'
	<< "inelastic cross section (mb): "
	<< inelXSec_/millibarn << '\n'
	// TODO: Include diffrative Cross Section in calculation
	// << "diffractive cross section (mb): "
	// << diffratio_*totalXSecExp()/millibarn << '\n'
	<< "soft inv radius (GeV2): "
	<< softMu2_/GeV2 << '\n'
	<< "slope of soft pt spectrum (1/GeV2): "
	<< beta_sqr(1.GeV) << '\n'
	<< "Average hard multiplicity: "
	<< avgNhard_ << '\n'
	<< "Average soft multiplicity: "
	<< avgNsoft_ << '\n' << line << endl;
	}else{
	file << line
	<< "Eikonalized and soft cross sections:\n\n"
	<< "Model parameters: "
	<< "ptmin: " << Ptmin_/GeV << " GeV"
	<< ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n"
	<< " "
	<< ", CMenergy: " << generator()->maximumCMEnergy()/GeV
	<< " GeV" << '\n'
	<< "hard inclusive cross section (mb): "
	<< hardXSec_/millibarn << '\n'
	<< "Average hard multiplicity: "
	<< avgNhard_ << '\n' << line << endl;
	}
	}

	unsigned int MPIHandler::multiplicity(unsigned int sel){
	if(sel==0){//draw from the pretabulated distribution
	MPair m = theMultiplicities.select(UseRandom::rnd());
	softMult_ = m.second;
	return m.first;
	} else{ //fixed multiplicities for the additional hard scatters
	if(additionalMultiplicities_.size() < sel)
	throw Exception() << "MPIHandler::multiplicity: process index "
	<< "is out of range"
	<< Exception::runerror;

	return additionalMultiplicities_[sel-1];
	}
	}


	void MPIHandler::Probs(XSVector UEXSecs) {
	GSLIntegrator integrator;
	unsigned int iH(1), iS(0);
	double P(0.0);
	double P0(0.0);//the probability for i hard and zero soft scatters
	Length bmax(500.0*sqrt(millibarn));
	//only one UE process will be drawn from a probability distribution,
	//so check that.
	assert(UEXSecs.size() == 1);

	for ( XSVector::const_iterator it = UEXSecs.begin();
	it != UEXSecs.end(); ++it ) {

	iH = 0;

	//get the inel xsec
	Eikonalization inelint(this, -1, *it, softXSec_, softMu2_);
	inelXSec_ = integrator.value(inelint, ZERO, bmax);

	avgNhard_ = 0.0;
	avgNsoft_ = 0.0;
	bmax = 10.0*sqrt(millibarn);
	do{//loop over hard ints
	if(Algorithm()>-1 && iH==0){
	iH++;
	continue;
	}
	iS = 0;
	do{//loop over soft ints
	if( ( Algorithm() == -1 && iS==0 && iH==0 ) ){
	iS++;
	continue;
	}

	Eikonalization integrand(this, iH100+iS, it, softXSec_, softMu2_);

	if(Algorithm() > 0){
	P = integrator.value(integrand, ZERO, bmax) / *it;
	}else{
	P = integrator.value(integrand, ZERO, bmax) / inelXSec_;
	}
	//store the probability
	if(Algorithm()>-1){
	theMultiplicities.insert(P, make_pair(iH-1, iS));
	avgNhard_ += P*(iH-1);
	}else{
	theMultiplicities.insert(P, make_pair(iH, iS));
	avgNhard_ += P*(iH);
	}
	avgNsoft_ += P*iS;
	if(iS==0)
	P0 = P;

	iS++;
	} while ( (iS < maxScatters_) && (iS < 5 \|\| P > 1.e-15 ) && softInt_ );
	iH++;
	} while ( (iH < maxScatters_) && (iH < 5 \|\| P0 > 1.e-15) );
	}
	}


	// calculate the integrand
	Length Eikonalization::operator() (Length b) const {
	unsigned int Nhard(0), Nsoft(0);
	//fac is just: d^2b=fac*db despite that large number
	Length fac(Constants::twopi*b);

	double chiTot(( theHandler->OverlapFunction(b)*hardXSec_ +
	theHandler->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0);

	//total cross section wanted
	if(theoption == -2) return 2 * fac * ( 1 - exp(-chiTot) );
	//inelastic cross section
	if(theoption == -1) return fac * ( 1 - exp(- 2.0 * chiTot) );

	if(theoption >= 0){
	//encode multiplicities as: N_hard*100 + N_soft
	Nhard = theoption/100;
	Nsoft = theoption%100;

	if(theHandler->Algorithm() > 0){
	//P_n*sigma_hard: n-1 extra scatters + 1 hard scatterer != hardXSec_
	return fac * Nhard * theHandler->poisson(b, hardXSec_, Nhard) *
	theHandler->poisson(b, softXSec_, Nsoft, softMu2_);
	}else{
	//P_n*sigma_inel: n extra scatters
	return fac * theHandler->poisson(b, hardXSec_, Nhard) *
	theHandler->poisson(b, softXSec_, Nsoft, softMu2_);
	}

	}else{
	throw Exception() << "Parameter theoption in Struct Eikonalization in "
	<< "MPIHandler.cc has not allowed value"
	<< Exception::runerror;
	return 0.0*meter;
	}
	}

	InvEnergy2 slopeBisection::operator() (Energy2 softMu2) const {
	GSLBisection root;
	//single component model
	TotalXSecBisection fn(handler_, softMu2);
	try{
	softXSec_ = root.value(fn, 0millibarn, 5000millibarn);
	}catch(GSLBisection::IntervalError){
	throw Exception() << "MPIHandler 2-Component model didn't work out."
	<< Exception::runerror;
	}
	return handler_->slopeDiff(softXSec_, softMu2);
	}

	LengthDiff slopeInt::operator() (Length b) const {
	//fac is just: d^2b=fac*db
	Length fac(Constants::twopi*b);

	double chiTot(( handler_->OverlapFunction(b)*hardXSec_ +
	handler_->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0);

	InvEnergy2 b2 = sqr(b/hbarc);
	//B*sigma_tot
	return fac * b2 * ( 1 - exp(-chiTot) );
	}

	InvArea MPIHandler::OverlapFunction(Length b, Energy2 mu2) const {
	if(mu2 == ZERO)
	mu2 = invRadius_;

	InvLength mu = sqrt(mu2)/hbarc;
	return (sqr(mu)/96/Constants::pi)pow(mub, 3)(gsl_sf_bessel_Kn(3, mub));
	}

	double MPIHandler::poisson(Length b, CrossSection sigma, unsigned int N, Energy2 mu2) const {
	if(sigma > 0*millibarn){
	return pow(OverlapFunction(b, mu2)*sigma, (double)N)/factorial(N)
	exp(-OverlapFunction(b, mu2)sigma);
	}else{
	return (N==0) ? 1.0 : 0.0;
	}
	}

	CrossSection MPIHandler::totalXSecDiff(CrossSection softXSec,
	Energy2 softMu2) const {
	GSLIntegrator integrator;
	Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2);
	Length bmax = 500.0*sqrt(millibarn);

	CrossSection tot = integrator.value(integrand, ZERO, bmax);
	return (tot-totalXSecExp());
	}

	InvEnergy2 MPIHandler::slopeDiff(CrossSection softXSec,
	Energy2 softMu2) const {
	GSLIntegrator integrator;
	Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2);
	Length bmax = 500.0*sqrt(millibarn);

	CrossSection tot = integrator.value(integrand, ZERO, bmax);

	slopeInt integrand2(this, hardXSec_, softXSec, softMu2);

	return integrator.value(integrand2, ZERO, bmax)/tot - slopeExp();
	}

	CrossSection MPIHandler::totalXSecExp() const {
	if(totalXSecExp_ != 0*millibarn)
	return totalXSecExp_;

	double pom_old = 0.0808;
	CrossSection coef_old = 21.7*millibarn;
	double pom_new_hard = 0.452;
	CrossSection coef_new_hard = 0.0139*millibarn;
	double pom_new_soft = 0.0667;
	CrossSection coef_new_soft = 24.22*millibarn;

	Energy energy(generator()->maximumCMEnergy());

	switch(DLmode_){
	case 1://old DL extrapolation
	return coef_old * pow(energy/GeV, 2*pom_old);
	break;

	case 2://old DL extrapolation fixed to CDF
	return 81.8millibarn pow(energy/1800.0/GeV, 2*pom_old);
	break;

	case 3://new DL extrapolation
	return 0.8coef_new_hard pow(energy/GeV, 2*pom_new_hard) +
	coef_new_soft * pow(energy/GeV, 2*pom_new_soft);
	break;

	default:
	throw Exception() << "MPIHandler::totalXSecExp non-existing mode selected"
	<< Exception::runerror;
	}
	}

	InvEnergy2 MPIHandler::slopeExp() const{
	//Currently return the slope as calculated in the pomeron fit by
	//Donnachie & Landshoff
	Energy energy(generator()->maximumCMEnergy());
	//slope at
	Energy e_0 = 1800*GeV;
	InvEnergy2 b_0 = 17/GeV2;
	return b_0 + log(energy/e_0)/GeV2;
	}

	void MPIHandler::overrideUECuts() {
	if(energyExtrapolation_==1)
	Ptmin_ = EEparamA_ * log(generator()->maximumCMEnergy() / EEparamB_);
	- else if(energyExtrapolation_==2)
	+ else if(energyExtrapolation_==2){
	+ // For small CME <900GeV the old power law does not work (hard cross section diverges)
	Ptmin_ = pT0_*pow(double(generator()->maximumCMEnergy()/refScale_),b_);
	+ }else if(energyExtrapolation_==3){
	+ // New parametrization that works for small CME: pTmin0*((c+sqrt(s))/E0)^b
	+ Ptmin_ = pT0_*pow(double((offset_+generator()->maximumCMEnergy())/refScale_),b_);
	+ cout << Ptmin_/GeV << endl;
	+ }
	else
	assert(false);
	// create a new SimpleKTCut object with the calculated ptmin value
	Ptr<SimpleKTCut>::pointer newUEktCut = new_ptr(SimpleKTCut(Ptmin_));
	newUEktCut->init();
	newUEktCut->initrun();

	// create a new Cuts object with MHatMin = 2 * Ptmin_
	CutsPtr newUEcuts = new_ptr(Cuts(2*Ptmin_));
	newUEcuts->add(dynamic_ptr_cast<tOneCutPtr>(newUEktCut));
	newUEcuts->init();
	newUEcuts->initrun();

	// replace the old Cuts object
	cuts()[0] = newUEcuts;
	}

	void MPIHandler::persistentOutput(PersistentOStream & os) const {
	os << theMultiplicities << theHandler
	<< theSubProcesses << theCuts << theProcessHandlers
	<< additionalMultiplicities_ << identicalToUE_
	<< ounit(PtOfQCDProc_, GeV) << ounit(Ptmin_, GeV)
	<< ounit(hardXSec_, millibarn) << ounit(softXSec_, millibarn)
	<< ounit(beta_, 1/GeV2)
	<< algorithm_ << ounit(invRadius_, GeV2)
	<< numSubProcs_ << colourDisrupt_ << softInt_ << twoComp_
	<< DLmode_ << ounit(totalXSecExp_, millibarn)
	<< energyExtrapolation_ << ounit(EEparamA_, GeV) << ounit(EEparamB_, GeV)
	- << ounit(refScale_,GeV) << ounit(pT0_,GeV) << b_
	+ << ounit(refScale_,GeV) << ounit(pT0_,GeV) << b_ << ounit(offset_,GeV)
	<< avgNhard_ << avgNsoft_ << softMult_
	<< ounit(inelXSec_, millibarn)
	<< ounit(softMu2_, GeV2) << diffratio_;
	}

	void MPIHandler::persistentInput(PersistentIStream & is, int) {
	is >> theMultiplicities >> theHandler
	>> theSubProcesses >> theCuts >> theProcessHandlers
	>> additionalMultiplicities_ >> identicalToUE_
	>> iunit(PtOfQCDProc_, GeV) >> iunit(Ptmin_, GeV)
	>> iunit(hardXSec_, millibarn) >> iunit(softXSec_, millibarn)
	>> iunit(beta_, 1/GeV2)
	>> algorithm_ >> iunit(invRadius_, GeV2)
	>> numSubProcs_ >> colourDisrupt_ >> softInt_ >> twoComp_
	>> DLmode_ >> iunit(totalXSecExp_, millibarn)
	>> energyExtrapolation_ >> iunit(EEparamA_, GeV) >> iunit(EEparamB_, GeV)
	- >> iunit(refScale_,GeV) >> iunit(pT0_,GeV) >> b_
	+ >> iunit(refScale_,GeV) >> iunit(pT0_,GeV) >> b_ >> iunit(offset_,GeV)
	>> avgNhard_ >> avgNsoft_ >> softMult_
	>> iunit(inelXSec_, millibarn)
	>> iunit(softMu2_, GeV2) >> diffratio_;
	currentHandler_ = this;
	}

	void MPIHandler::clean() {
	// ThePEG's event handler's usual event cleanup doesn't reach these
	// XCombs. Need to do it by hand here.
	for ( size_t i = 0; i < theSubProcesses.size(); ++i ) {
	theSubProcesses[i]->pExtractor()->lastXCombPtr()->clean();
	}
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<MPIHandler,Interfaced>
	describeHerwigMPIHandler("Herwig::MPIHandler",
	"JetCuts.so SimpleKTCut.so HwMPI.so");

	void MPIHandler::Init() {

	static ClassDocumentation<MPIHandler> documentation
	("The MPIHandler class is the main administrator of the multiple interaction model",
	"The underlying event was simulated with an eikonal model for multiple partonic interactions."
	"Details can be found in Ref.~\\cite{Bahr:2008dy,Bahr:2009ek}.",
	"%\\cite{Bahr:2008dy}\n"
	"\\bibitem{Bahr:2008dy}\n"
	" M.~Bahr, S.~Gieseke and M.~H.~Seymour,\n"
	" ``Simulation of multiple partonic interactions in Herwig,''\n"
	" JHEP {\\bf 0807}, 076 (2008)\n"
	" [arXiv:0803.3633 [hep-ph]].\n"
	" %%CITATION = JHEPA,0807,076;%%\n"
	"\\bibitem{Bahr:2009ek}\n"
	" M.~Bahr, J.~M.~Butterworth, S.~Gieseke and M.~H.~Seymour,\n"
	" ``Soft interactions in Herwig,''\n"
	" arXiv:0905.4671 [hep-ph].\n"
	" %%CITATION = ARXIV:0905.4671;%%\n"
	);

	static RefVector<MPIHandler,SubProcessHandler> interfaceSubhandlers
	("SubProcessHandlers",
	"The list of sub-process handlers used in this EventHandler. ",
	&MPIHandler::theSubProcesses, -1, false, false, true, false, false);

	static RefVector<MPIHandler,Cuts> interfaceCuts
	("Cuts",
	"List of cuts used for the corresponding list of subprocesses. These cuts "
	"should not be overidden in individual sub-process handlers.",
	&MPIHandler::theCuts, -1, false, false, true, false, false);

	static Parameter<MPIHandler,Energy2> interfaceInvRadius
	("InvRadius",
	"The inverse hadron radius squared used in the overlap function",
	&MPIHandler::invRadius_, GeV2, 2.0GeV2, 0.2GeV2, 4.0*GeV2,
	true, false, Interface::limited);

	static ParVector<MPIHandler,int> interfaceadditionalMultiplicities
	("additionalMultiplicities",
	"specify the multiplicities of secondary hard processes (multiple parton scattering)",
	&MPIHandler::additionalMultiplicities_,
	-1, 0, 0, 3,
	false, false, true);

	static Parameter<MPIHandler,int> interfaceIdenticalToUE
	("IdenticalToUE",
	"Specify which of the hard processes is identical to the UE one (QCD dijets)",
	&MPIHandler::identicalToUE_, -1, 0, 0,
	false, false, Interface::nolimits);

	static Parameter<MPIHandler,Energy> interfacePtOfQCDProc
	("PtOfQCDProc",
	"Specify the value of the pt cutoff for the process that is identical to the UE one",
	&MPIHandler::PtOfQCDProc_, GeV, -1.0*GeV, ZERO, ZERO,
	false, false, Interface::nolimits);

	static Parameter<MPIHandler,double> interfacecolourDisrupt
	("colourDisrupt",
	"Fraction of connections to additional subprocesses, which are colour disrupted.",
	&MPIHandler::colourDisrupt_,
	0.0, 0.0, 1.0,
	false, false, Interface::limited);


	static Parameter<MPIHandler,double> interfaceDiffRatio
	("DiffractiveRatio",
	"Fraction of diffractive cross section in inelastic cross section.",
	&MPIHandler::diffratio_,
	0.2, 0.0, 1.0,
	false, false, Interface::limited);


	static Switch<MPIHandler,bool> interfacesoftInt
	("softInt",
	"Switch to enable soft interactions",
	&MPIHandler::softInt_, true, false, false);

	static SwitchOption interfacesoftIntYes
	(interfacesoftInt,
	"Yes",
	"enable the two component model",
	true);
	static SwitchOption interfacesoftIntNo
	(interfacesoftInt,
	"No",
	"disable the model",
	false);


	static Switch<MPIHandler,unsigned int> interEnergyExtrapolation
	("EnergyExtrapolation",
	"Switch to ignore the cuts object at MPIHandler:Cuts[0]. "
	"Instead, extrapolate the pt cut.",
	&MPIHandler::energyExtrapolation_, 2, false, false);
	static SwitchOption interEnergyExtrapolationLog
	(interEnergyExtrapolation,
	"Log",
	"Use logarithmic dependence, ptmin = A * log (sqrt(s) / B).",
	1);
	static SwitchOption interEnergyExtrapolationPower
	(interEnergyExtrapolation,
	"Power",
	"Use power law, ptmin = pt_0 * (sqrt(s) / E_0)^b.",
	2);
	+ static SwitchOption interEnergyExtrapolationPowerModified
	+ (interEnergyExtrapolation,
	+ "PowerModified",
	+ "Use modified power law with offset to work for small center of mass energies."
	+ "ptmin = pt_0 * ((offset+sqrt(s)) / E_0)^b.",
	+ 3);
	+
	static SwitchOption interEnergyExtrapolationNo
	(interEnergyExtrapolation,
	"No",
	"Use manually set value for the minimal pt, "
	"specified in MPIHandler:Cuts[0]:OneCuts[0]:MinKT.",
	0);

	static Parameter<MPIHandler,Energy> interfaceEEparamA
	("EEparamA",
	"Parameter A in the empirical parametrization "
	"ptmin = A * log (sqrt(s) / B)",
	&MPIHandler::EEparamA_, GeV, 0.6*GeV, ZERO, Constants::MaxEnergy,
	false, false, Interface::limited);

	static Parameter<MPIHandler,Energy> interfaceEEparamB
	("EEparamB",
	"Parameter B in the empirical parametrization "
	"ptmin = A * log (sqrt(s) / B)",
	&MPIHandler::EEparamB_, GeV, 39.0*GeV, ZERO, Constants::MaxEnergy,
	false, false, Interface::limited);

	static Switch<MPIHandler,bool> interfacetwoComp
	("twoComp",
	"switch to enable the model with a different radius for soft interactions",
	&MPIHandler::twoComp_, true, false, false);

	static SwitchOption interfacetwoCompYes
	(interfacetwoComp,
	"Yes",
	"enable the two component model",
	true);
	static SwitchOption interfacetwoCompNo
	(interfacetwoComp,
	"No",
	"disable the model",
	false);


	static Parameter<MPIHandler,CrossSection> interfaceMeasuredTotalXSec
	("MeasuredTotalXSec",
	"Value for the total cross section (assuming rho=0). If non-zero, this "
	"overwrites the Donnachie-Landshoff parametrizations.",
	&MPIHandler::totalXSecExp_, millibarn, 0.0millibarn, 0.0millibarn, 0*millibarn,
	false, false, Interface::lowerlim);


	static Switch<MPIHandler,unsigned int> interfaceDLmode
	("DLmode",
	"Choice of Donnachie-Landshoff parametrization for the total cross section.",
	&MPIHandler::DLmode_, 2, false, false);
	static SwitchOption interfaceDLmodeStandard
	(interfaceDLmode,
	"Standard",
	"Standard parametrization with s**0.08",
	1);
	static SwitchOption interfaceDLmodeCDF
	(interfaceDLmode,
	"CDF",
	"Standard parametrization but normalization fixed to CDF's measured value",
	2);
	static SwitchOption interfaceDLmodeNew
	(interfaceDLmode,
	"New",
	"Parametrization taking hard and soft pomeron contributions into account",
	3);

	static Parameter<MPIHandler,Energy> interfaceReferenceScale
	("ReferenceScale",
	"The reference energy for power law energy extrapolation of pTmin",
	&MPIHandler::refScale_, GeV, 7000.0GeV, 0.0GeV, 20000.*GeV,
	false, false, Interface::limited);

	static Parameter<MPIHandler,Energy> interfacepTmin0
	("pTmin0",
	"The pTmin at the reference scale for power law extrapolation of pTmin.",
	&MPIHandler::pT0_, GeV, 3.11GeV, 0.0GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<MPIHandler,double> interfacePower
	("Power",
	"The power for power law extrapolation of the pTmin cut-off.",
	&MPIHandler::b_, 0.21, 0.0, 10.0,
	false, false, Interface::limited);
	+
	+ static Parameter<MPIHandler,Energy> interfaceOffset
	+ ("Offset",
	+ "The offset used in the modified power law extrapolation of the pTmin cut-off.",
	+ &MPIHandler::offset_, GeV, 622.0GeV, 500.0GeV, 1000.0*GeV,
	+ false, false, Interface::limited);
	+
	}
	+
	diff --git a/UnderlyingEvent/MPIHandler.h b/UnderlyingEvent/MPIHandler.h
	--- a/UnderlyingEvent/MPIHandler.h
	+++ b/UnderlyingEvent/MPIHandler.h
	@@ -1,897 +1,898 @@
	// -- C++ --
	//
	// MPIHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2017 The Herwig Collaboration
	//
	// Herwig is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_MPIHandler_H
	#define HERWIG_MPIHandler_H
	//
	// This is the declaration of the MPIHandler class.
	//
	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Herwig/PDT/StandardMatchers.h"
	#include "Herwig/Utilities/GSLBisection.h"
	//#include "Herwig/Utilities/GSLMultiRoot.h"
	#include "Herwig/Utilities/GSLIntegrator.h"
	#include "Herwig/Shower/UEBase.h"

	#include <cassert>
	#include "ProcessHandler.h"
	#include "MPIHandler.fh"


	namespace Herwig {
	using namespace ThePEG;

	/** \ingroup UnderlyingEvent
	* \class MPIHandler
	* This class is responsible for generating additional
	* semi hard partonic interactions.
	*
	* \author Manuel B\"ahr
	*
	* @see \ref MPIHandlerInterfaces "The interfaces"
	* defined for MPIHandler.
	* @see ProcessHandler
	* @see ShowerHandler
	* @see HwRemDecayer
	*/

	class MPIHandler: public UEBase {

	/**
	* Maximum number of scatters
	*/
	static const unsigned int maxScatters_ = 99;

	/**
	* Class for the integration is a friend to access private members
	*/
	friend struct Eikonalization;
	friend struct TotalXSecBisection;
	friend struct slopeAndTotalXSec;
	friend struct slopeInt;
	friend struct slopeBisection;

	public:

	/** A vector of <code>SubProcessHandler</code>s. */
	typedef vector<SubHdlPtr> SubHandlerList;

	/** A vector of <code>Cut</code>s. */
	typedef vector<CutsPtr> CutsList;

	/** A vector of <code>ProcessHandler</code>s. */
	typedef vector<ProHdlPtr> ProcessHandlerList;

	/** A vector of cross sections. */
	typedef vector<CrossSection> XSVector;

	/** A pair of multiplicities: hard, soft. */
	typedef pair<unsigned int, unsigned int> MPair;

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MPIHandler(): softMult_(0), identicalToUE_(-1),
	PtOfQCDProc_(-1.0GeV), Ptmin_(-1.0GeV),
	hardXSec_(0millibarn), softXSec_(0millibarn),
	totalXSecExp_(0*millibarn),
	softMu2_(ZERO), beta_(100.0/GeV2),
	algorithm_(2), numSubProcs_(0),
	colourDisrupt_(0.0), softInt_(true), twoComp_(true),
	DLmode_(2), avgNhard_(0.0), avgNsoft_(0.0),
	- energyExtrapolation_(2), EEparamA_(0.6*GeV),
	+ energyExtrapolation_(3), EEparamA_(0.6*GeV),
	EEparamB_(37.5GeV), refScale_(7000.GeV),
	- pT0_(3.11*GeV), b_(0.21) {}
	+ pT0_(2.875GeV), b_(0.3101), offset_(622.204GeV) {}

	/**
	* The destructor.
	*/
	virtual ~MPIHandler(){}
	//@}

	public:

	/** @name Methods for the MPI generation. */
	//@{

	/*
	* @return true if for this beam setup MPI can be generated
	*/
	virtual bool beamOK() const;

	/**
	* Return true or false depending on whether soft interactions are enabled.
	*/
	virtual bool softInt() const {return softInt_;}

	/**
	* Get the soft multiplicity from the pretabulated multiplicity
	* distribution. Generated in multiplicity in the first place.
	* @return the number of extra soft events in this collision
	*/
	virtual unsigned int softMultiplicity() const {return softMult_;}

	/**
	* Sample from the pretabulated multiplicity distribution.
	* @return the number of extra events in this collision
	*/
	virtual unsigned int multiplicity(unsigned int sel=0);

	/**
	* Select a StandardXComb according to it's weight
	* @return that StandardXComb Object
	* @param sel is the subprocess that should be returned,
	* if more than one is specified.
	*/
	virtual tStdXCombPtr generate(unsigned int sel=0);
	//@}


	/** @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();

	/**
	* Initialize this Multiple Interaction handler and all related objects needed to
	* generate additional events.
	*/
	virtual void initialize();

	/**
	* Finalize this Multiple Interaction handler and all related objects needed to
	* generate additional events.
	*/
	virtual void finalize();

	/**
	* Clean up the XCombs from our subprocesses after each event.
	* ThePEG cannot see them, so the usual cleaning misses these.
	*/
	virtual void clean();

	/**
	* Write out accumulated statistics about integrated cross sections.
	*/
	void statistics() const;

	/**
	* The level of statistics. Controlls the amount of statistics
	* written out after each run to the <code>EventGenerator</code>s
	* <code>.out</code> file. Simply the EventHandler method is called here.
	*/
	int statLevel() const {return eventHandler()->statLevel();}

	/**
	* Return the hard cross section above ptmin
	*/
	CrossSection hardXSec() const { return hardXSec_; }

	/**
	* Return the soft cross section below ptmin
	*/
	CrossSection softXSec() const { return softXSec_; }

	/**
	* Return the inelastic cross section
	*/
	CrossSection inelasticXSec() const { return inelXSec_; }

	/**
	* Return the non-diffractive cross section assumed by the model.
	* TODO: See comment at diffractiveXSec.
	*/
	CrossSection nonDiffractiveXSec() const {
	return (1.-diffratio_)*inelXSec_;
	}

	/**
	* Return the diffractive cross section assumed by the model.
	* For now the diffractive cross section is seen as a fixed part of the
	* inelastic cross section.
	* TODO: Energy dependence and/or Include diffraction in Eikonalisation.
	*/
	CrossSection diffractiveXSec() const {
	return diffratio_*inelXSec_;
	}


	/** @name Simple access functions. */
	//@{

	/**
	* Return the ThePEG::EventHandler assigned to this handler.
	* This methods shadows ThePEG::StepHandler::eventHandler(), because
	* it is not virtual in ThePEG::StepHandler. This is ok, because this
	* method would give a null-pointer at some stages, whereas this method
	* gives access to the explicitely copied pointer (in initialize())
	* to the ThePEG::EventHandler.
	*/
	tEHPtr eventHandler() const {return theHandler;}

	/**
	* Return the current handler
	*/
	static const MPIHandler * currentHandler() {
	return currentHandler_;
	}

	/**
	* Return theAlgorithm.
	*/
	virtual int Algorithm() const {return algorithm_;}

	/**
	* Return the ptmin parameter of the model
	*/
	virtual Energy Ptmin() const {
	if(Ptmin_ > ZERO)
	return Ptmin_;
	else
	throw Exception() << "MPIHandler::Ptmin called without initialize before"
	<< Exception::runerror;
	}

	/**
	* Return the slope of the soft pt spectrum as calculated.
	*/
	virtual InvEnergy2 beta() const {
	if(beta_ != 100.0/GeV2)
	return beta_;
	else
	throw Exception() << "MPIHandler::beta called without initialization"
	<< Exception::runerror;
	}

	/**
	* Return the pt Cutoff of the Interaction that is identical to the UE
	* one.
	*/
	virtual Energy PtForVeto() const {return PtOfQCDProc_;}

	/**
	* Return the number of additional "hard" processes ( = multiple
	* parton scattering)
	*/
	virtual unsigned int additionalHardProcs() const {return numSubProcs_-1;}

	/**
	* Return the fraction of colour disrupted connections to the
	* suprocesses.
	*/
	virtual double colourDisrupt() const {return colourDisrupt_;}

	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;
	//@}

	private:

	/**
	* Access the list of sub-process handlers.
	*/
	const SubHandlerList & subProcesses()
	const {return theSubProcesses;}

	/**
	* Access the list of sub-process handlers.
	*/
	SubHandlerList & subProcesses() {return theSubProcesses;}

	/**
	* Access the list of cuts.
	*/
	const CutsList & cuts() const {return theCuts;}

	/**
	* Access the list of cuts.
	*/
	CutsList & cuts() {return theCuts;}

	/**
	* Access the list of sub-process handlers.
	*/
	const ProcessHandlerList & processHandlers()
	const {return theProcessHandlers;}

	/**
	* Access the list of sub-process handlers.
	*/
	ProcessHandlerList & processHandlers() {return theProcessHandlers;}


	/**
	* Method to calculate the individual probabilities for N scatters in the event.
	* @param UEXSecs is(are) the inclusiv cross section(s) for the UE process(es).
	*/
	void Probs(XSVector UEXSecs);

	/**
	* Debug method to check the individual probabilities.
	* @param filename is the file the output gets written to
	*/
	void MultDistribution(string filename) const;

	/**
	* Return the value of the Overlap function A(b) for a given impact
	* parameter \a b.
	* @param b impact parameter
	* @param mu2 = inv hadron radius squared. 0 will use the value of
	* invRadius_
	* @return inverse area.
	*/
	InvArea OverlapFunction(Length b, Energy2 mu2=ZERO) const;

	/**
	* Method to calculate the poisson probability for expectation value
	* \f$<n> = A(b)\sigma\f$, and multiplicity N.
	*/
	double poisson(Length b, CrossSection sigma,
	unsigned int N, Energy2 mu2=ZERO) const;

	/**
	* Returns the total cross section for the current CMenergy. The
	* decision which parametrization will be used is steered by a
	* external parameter of this class.
	*/
	CrossSection totalXSecExp() const;

	/**
	* Difference of the calculated total cross section and the
	* experimental one from totalXSecExp.
	* @param softXSec = the soft cross section that is used
	* @param softMu2 = the soft radius, if 0 the hard radius will be used
	*/
	CrossSection totalXSecDiff(CrossSection softXSec,
	Energy2 softMu2=ZERO) const;

	/**
	* Difference of the calculated elastic slope and the
	* experimental one from slopeExp.
	* @param softXSec = the soft cross section that is used
	* @param softMu2 = the soft radius, if 0 the hard radius will be used
	*/
	InvEnergy2 slopeDiff(CrossSection softXSec,
	Energy2 softMu2=ZERO) const;

	/**
	* Returns the value of the elastic slope for the current CMenergy.
	* The decision which parametrization will be used is steered by a
	* external parameter of this class.
	*/
	InvEnergy2 slopeExp() const;


	/**
	* Calculate the minimal transverse momentum from the extrapolation
	*/
	void overrideUECuts();


	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	MPIHandler & operator=(const MPIHandler &) = delete;

	/**
	* A pointer to the EventHandler that calls us. Has to be saved, because the
	* method eventHandler() inherited from ThePEG::StepHandler returns a null-pointer
	* sometimes. Leif changed that in r1053 so that a valid pointer is present, when
	* calling doinitrun().
	*/
	tEHPtr theHandler;

	/**
	* The list of <code>SubProcessHandler</code>s.
	*/
	SubHandlerList theSubProcesses;

	/**
	* The kinematical cuts used for this collision handler.
	*/
	CutsList theCuts;

	/**
	* List of ProcessHandler used to sample different processes independently
	*/
	ProcessHandlerList theProcessHandlers;

	/**
	* A ThePEG::Selector where the individual Probabilities P_N are stored
	* and the actual Multiplicities can be selected.
	*/
	Selector<MPair> theMultiplicities;

	/**
	* Variable to store the soft multiplicity generated for a event. This
	* has to be stored as it is generated at the time of the hard
	* additional interactions but used later on.
	*/
	unsigned int softMult_;

	/**
	* Variable to store the multiplicity of the second hard process
	*/
	vector<int> additionalMultiplicities_;

	/**
	* Variable to store the information, which process is identical to
	* the UE one (QCD dijets).
	* 0 means "real" hard one
	* n>0 means the nth additional hard scatter
	* -1 means no one!
	*/
	int identicalToUE_;

	/**
	* Variable to store the minimal pt of the process that is identical
	* to the UE one. This only has to be set, if it can't be determined
	* automatically (i.e. when reading QCD LesHouches files in).
	*/
	Energy PtOfQCDProc_;

	/**
	* Variable to store the parameter ptmin
	*/
	Energy Ptmin_;

	/**
	* Variable to store the hard cross section above ptmin
	*/
	CrossSection hardXSec_;

	/**
	* Variable to store the final soft cross section below ptmin
	*/
	CrossSection softXSec_;

	/**
	* Variable to store the inelastic cross section
	*/
	CrossSection inelXSec_;

	/**
	* Variable to store the total pp cross section (assuming rho=0!) as
	* measured at LHC. If this variable is set, this value is used in the
	* subsequent run instead of any of the Donnachie-Landshoff
	* parametrizations.
	*/
	CrossSection totalXSecExp_;

	/**
	* Variable to store the soft radius, that is calculated during
	* initialization for the two-component model.
	*/
	Energy2 softMu2_;

	/**
	* slope to the non-perturbative pt spectrum: \f$d\sigma/dp_T^2 = A \exp
	* (- beta p_T^2)\f$. Its value is determined durint initialization.
	*/
	InvEnergy2 beta_;
	/**
	* Switch to be set from outside to determine the algorithm used for
	* UE activity.
	*/
	int algorithm_;

	/**
	* Inverse hadron Radius squared \f$ (\mu^2) \f$. Used inside the overlap function.
	*/
	Energy2 invRadius_;

	/**
	* Member variable to store the actual number of separate SubProcesses
	*/
	unsigned int numSubProcs_;

	/**
	* Variable to store the relative number of colour disrupted
	* connections to additional subprocesses. This variable is used in
	* Herwig::HwRemDecayer but store here, to have access to all
	* parameters through one Object.
	*/
	double colourDisrupt_;

	/**
	* Flag to store whether soft interactions, i.e. pt < ptmin should be
	* simulated.
	*/
	bool softInt_;

	/**
	* Flag to steer wheather the soft part has a different radius, that
	* will be dynamically fixed.
	*/
	bool twoComp_;

	/**
	* Switch to determine which Donnachie & Landshoff parametrization
	* should be used.
	*/
	unsigned int DLmode_;

	/**
	* Variable to store the average hard multiplicity.
	*/
	double avgNhard_;

	/**
	* Variable to store the average soft multiplicity.
	*/
	double avgNsoft_;

	/**
	* The current handler
	*/
	static MPIHandler * currentHandler_;

	/**
	* Flag to store whether to calculate the minimal UE pt according to an
	* extrapolation formula or whether to use MPIHandler:Cuts[0]:OneCuts[0]:MinKT
	*/
	unsigned int energyExtrapolation_;

	/**
	* Parameters for the energy extrapolation formula
	*/
	Energy EEparamA_;
	Energy EEparamB_;
	Energy refScale_;
	Energy pT0_;
	double b_;
	+ Energy offset_;

	/**
	* Parameters to set the fraction of diffractive cross section in the inelastic cross section.
	*/
	double diffratio_=0.2;

	protected:

	/** @cond EXCEPTIONCLASSES */

	/**
	* Exception class used by the MultipleInteractionHandler, when something
	* during initialization went wrong.
	* \todo understand!!!
	*/
	class InitError: public Exception {};

	/** @endcond */

	};

	}

	namespace Herwig {

	/**
	* A struct for the 2D root finding that is necessary to determine the
	* soft cross section and the soft radius that is needed to describe
	* the total cross section correctly.
	* NOT IN USE CURRENTLY
	*/
	struct slopeAndTotalXSec : public GSLHelper<CrossSection, CrossSection> {

	public:

	/**
	* Constructor
	*/
	slopeAndTotalXSec(tcMPIHPtr handler): handler_(handler) {}

	/** second argument type */
	typedef Energy2 ArgType2;

	/** second value type */
	typedef InvEnergy2 ValType2;

	/** first element of the vector like function to find root for
	* @param softXSec soft cross-section
	* @param softMu2 \f$\mu^2\f$
	*/
	CrossSection f1(ArgType softXSec, ArgType2 softMu2) const {
	return handler_->totalXSecDiff(softXSec, softMu2);
	}

	/** second element of the vector like function to find root for
	* @param softXSec soft cross-section
	* @param softMu2 \f$\mu^2\f$
	*/
	InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const {
	return handler_->slopeDiff(softXSec, softMu2);
	}

	/** provide the actual units of use */
	virtual ValType vUnit() const {return 1.0*millibarn;}

	/** otherwise rounding errors may get significant */
	virtual ArgType aUnit() const {return 1.0*millibarn;}

	/** provide the actual units of use */
	ValType2 vUnit2() const {return 1.0/GeV2;}

	/** otherwise rounding errors may get significant */
	ArgType2 aUnit2() const {return GeV2;}

	private:

	/**
	* Pointer to the handler
	*/
	tcMPIHPtr handler_;

	};

	/**
	* A struct for the root finding that is necessary to determine the
	* slope of the soft pt spectrum to match the soft cross section
	*/
	struct betaBisection : public GSLHelper<Energy2, InvEnergy2>{
	public:
	/**
	* Constructor.
	* @param soft = soft cross section, i.e. the integral of the soft
	* pt spectrum f(u=p_T^2) = dsig exp(-beta*u/u_min)
	* @param dsig = dsigma_hard/dp_T^2 at the p_T cutoff
	* @param ptmin = p_T cutoff
	*/
	betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin)
	: softXSec_(soft), dsig_(dsig), ptmin_(ptmin) {}

	/**
	* Operator that is used inside the GSLBisection class
	*/
	virtual Energy2 operator ()(InvEnergy2 beta) const
	{
	if( fabs(beta*GeV2) < 1.E-4 )
	beta = (beta > ZERO) ? 1.E-4/GeV2 : -1.E-4/GeV2;

	return (exp(beta*sqr(ptmin_)) - 1.0)/beta - softXSec_/dsig_;
	}

	/** provide the actual units of use */
	virtual ValType vUnit() const {return 1.0*GeV2;}

	/** provide the actual units of use */
	virtual ArgType aUnit() const {return 1.0/GeV2;}

	private:

	/** soft cross section */
	CrossSection softXSec_;

	/** dsigma/dp_T^2 at ptmin */
	DiffXSec dsig_;

	/** pt cutoff */
	Energy ptmin_;
	};

	/**
	* A struct for the root finding that is necessary to determine the
	* soft cross section and soft mu2 that are needed to describe the
	* total cross section AND elastic slope correctly.
	*/
	struct slopeBisection : public GSLHelper<InvEnergy2, Energy2> {
	public:
	/** Constructor */
	slopeBisection(tcMPIHPtr handler) : handler_(handler) {}

	/**
	* Return the difference of the calculated elastic slope to the
	* experimental one for a given value of the soft mu2. During that,
	* the soft cross section get fixed.
	*/
	InvEnergy2 operator ()(Energy2 arg) const;

	/** Return the soft cross section that has been calculated */
	CrossSection softXSec() const {return softXSec_;}

	private:
	/** const pointer to the MPIHandler to give access to member functions.*/
	tcMPIHPtr handler_;
	/** soft cross section that is determined on the fly.*/
	mutable CrossSection softXSec_;
	};

	/**
	* A struct for the root finding that is necessary to determine the
	* soft cross section that is needed to describe the total cross
	* section correctly.
	*/
	struct TotalXSecBisection : public GSLHelper<CrossSection, CrossSection> {
	public:

	/**
	* Constructor
	* @param handler The handler
	* @param softMu2 \f$\mu^2\f$
	*/
	TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO):
	handler_(handler), softMu2_(softMu2) {}

	/**
	* operator to return the cross section
	* @param argument input cross section
	*/
	CrossSection operator ()(CrossSection argument) const {
	return handler_->totalXSecDiff(argument, softMu2_);
	}

	/** provide the actual units of use */
	virtual ValType vUnit() const {return 1.0*millibarn;}

	/** otherwise rounding errors may get significant */
	virtual ArgType aUnit() const {return 1.0*millibarn;}

	private:

	/**
	* The handler
	*/
	tcMPIHPtr handler_;

	/**
	* \f$\mu^2\f$
	*/
	Energy2 softMu2_;

	};

	/**
	* Typedef for derivative of the length
	*/
	typedef decltype(mm/GeV2) LengthDiff;

	/**
	* A struct for the integrand for the slope
	*/
	struct slopeInt : public GSLHelper<LengthDiff, Length>{

	public:
	/** Constructor
	* @param handler The handler
	* @param hard The hard cross section
	* @param soft The soft cross section
	* @param softMu2 \f$\mu^2\f$
	*/
	slopeInt(tcMPIHPtr handler, CrossSection hard,
	CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
	: handler_(handler), hardXSec_(hard),
	softXSec_(soft), softMu2_(softMu2) {}

	/**
	* Operator to return the answer
	* @param arg The argument
	*/
	ValType operator ()(ArgType arg) const;

	private:

	/**
	* Pointer to the Handler that calls this integrand
	*/
	tcMPIHPtr handler_;

	/**
	* The hard cross section to be eikonalized
	*/
	CrossSection hardXSec_;

	/**
	* The soft cross section to be eikonalized. Default is zero
	*/
	CrossSection softXSec_;

	/**
	* The inv radius^2 of the soft interactions.
	*/
	Energy2 softMu2_;

	};

	/**
	* A struct for the eikonalization of the inclusive cross section.
	*/
	struct Eikonalization : public GSLHelper<Length, Length>{

	/**
	* The constructor
	* @param handler is the pointer to the MPIHandler to get access to
	* MPIHandler::OverlapFunction and member variables of the MPIHandler.
	* @param option is a flag, whether the inelastic or the total
	* @param handler The handler
	* @param hard The hard cross section
	* @param soft The soft cross section
	* @param softMu2 \f$\mu^2\f$
	* cross section should be returned (-2 or -1). For option = N > 0 the integrand
	* is N(A(b)sigma)^N/N! exp(-A(b)sigma) this is the P_Nsigma where
	* P_N is the Probability of having exactly N interaction (including the hard one)
	* This is equation 14 from "Jimmy4: Multiparton Interactions in HERWIG for the LHC"
	*/
	Eikonalization(tcMPIHPtr handler, int option, CrossSection hard,
	CrossSection soft=0*millibarn, Energy2 softMu2=ZERO)
	: theHandler(handler), theoption(option), hardXSec_(hard),
	softXSec_(soft), softMu2_(softMu2) {}

	/**
	* Get the function value
	*/
	Length operator ()(Length argument) const;

	private:
	/**
	* Pointer to the Handler that calls this integrand
	*/
	tcMPIHPtr theHandler;

	/**
	* A flag to switch between the calculation of total and inelastic cross section
	* or calculations for the individual probabilities. See the constructor
	*/
	int theoption;

	/**
	* The hard cross section to be eikonalized
	*/
	CrossSection hardXSec_;

	/**
	* The soft cross section to be eikonalized. Default is zero
	*/
	CrossSection softXSec_;

	/**
	* The inv radius^2 of the soft interactions.
	*/
	Energy2 softMu2_;

	};
	}

	#endif /* HERWIG_MPIHandler_H */
	diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in
	--- a/src/defaults/Hadronization.in
	+++ b/src/defaults/Hadronization.in
	@@ -1,96 +1,96 @@
	# -- ThePEG-repository --

	############################################################
	# Setup of default hadronization
	#
	# There are no user servicable parts inside.
	#
	# Anything that follows below should only be touched if you
	# know what you're doing.
	#############################################################

	cd /Herwig/Particles
	create ThePEG::ParticleData Cluster
	setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1
	create ThePEG::ParticleData Remnant
	setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1

	mkdir /Herwig/Hadronization
	cd /Herwig/Hadronization

	create Herwig::ClusterHadronizationHandler ClusterHadHandler
	create Herwig::PartonSplitter PartonSplitter
	create Herwig::ClusterFinder ClusterFinder
	create Herwig::ColourReconnector ColourReconnector
	create Herwig::ClusterFissioner ClusterFissioner
	create Herwig::LightClusterDecayer LightClusterDecayer
	create Herwig::ClusterDecayer ClusterDecayer
	create Herwig::HwppSelector HadronSelector

	newdef ClusterHadHandler:PartonSplitter PartonSplitter
	newdef ClusterHadHandler:ClusterFinder ClusterFinder
	newdef ClusterHadHandler:ColourReconnector ColourReconnector
	newdef ClusterHadHandler:ClusterFissioner ClusterFissioner
	newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer
	newdef ClusterHadHandler:ClusterDecayer ClusterDecayer
	newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2
	newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter
	newdef ClusterHadHandler:UnderlyingEventHandler NULL

	newdef ClusterFissioner:HadronSelector HadronSelector
	newdef LightClusterDecayer:HadronSelector HadronSelector
	newdef ClusterDecayer:HadronSelector HadronSelector

	newdef ColourReconnector:ColourReconnection Yes
	-newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.57
	-newdef ColourReconnector:ReconnectionProbability 0.89
	+newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7
	+newdef ColourReconnector:ReconnectionProbability 0.95
	newdef ColourReconnector:Algorithm Baryonic
	newdef ColourReconnector:OctetTreatment All

	# Clustering parameters for light quarks
	newdef ClusterFissioner:ClMaxLight 3.649
	newdef ClusterFissioner:ClPowLight 2.780
	newdef ClusterFissioner:PSplitLight 0.899
	newdef ClusterDecayer:ClDirLight 1
	newdef ClusterDecayer:ClSmrLight 0.78

	# Clustering parameters for b-quarks
	newdef ClusterFissioner:ClMaxBottom 3.757
	newdef ClusterFissioner:ClPowBottom 0.547
	newdef ClusterFissioner:PSplitBottom 0.625
	newdef ClusterDecayer:ClDirBottom 1
	newdef ClusterDecayer:ClSmrBottom 0.078
	newdef HadronSelector:SingleHadronLimitBottom 0.000

	# Clustering parameters for c-quarks
	newdef ClusterFissioner:ClMaxCharm 3.950
	newdef ClusterFissioner:ClPowCharm 2.559
	newdef ClusterFissioner:PSplitCharm 0.994
	newdef ClusterDecayer:ClDirCharm 1
	newdef ClusterDecayer:ClSmrCharm 0.163
	newdef HadronSelector:SingleHadronLimitCharm 0.000

	# Clustering parameters for exotic quarks
	# (e.g. hadronizing Susy particles)
	newdef ClusterFissioner:ClMaxExotic 2.7*GeV
	newdef ClusterFissioner:ClPowExotic 1.46
	newdef ClusterFissioner:PSplitExotic 1.00
	newdef ClusterDecayer:ClDirExotic 1
	newdef ClusterDecayer:ClSmrExotic 0.
	newdef HadronSelector:SingleHadronLimitExotic 0.

	#
	newdef PartonSplitter:SplitPwtSquark 0.824135
	newdef PartonSplitter:Split uds

	#
	newdef HadronSelector:PwtDquark 1.0
	newdef HadronSelector:PwtUquark 1.0
	newdef HadronSelector:PwtSquark 0.291717
	newdef HadronSelector:PwtCquark 0.0
	newdef HadronSelector:PwtBquark 0.0
	newdef HadronSelector:PwtDIquark 0.298
	newdef HadronSelector:SngWt 0.74
	newdef HadronSelector:DecWt 0.62
	newdef HadronSelector:Mode 1
	newdef HadronSelector:BelowThreshold All
	diff --git a/src/defaults/UnderlyingEvent.in b/src/defaults/UnderlyingEvent.in
	--- a/src/defaults/UnderlyingEvent.in
	+++ b/src/defaults/UnderlyingEvent.in
	@@ -1,88 +1,89 @@
	# -- ThePEG-repository --

	################################################
	# Set up the handler for the underlying event.
	################################################

	mkdir /Herwig/UnderlyingEvent
	cd /Herwig/UnderlyingEvent

	######## cuts #################################################
	library SimpleKTCut.so

	# cut on pt. Without a specific matcher object, it works on all
	# particles
	create ThePEG::SimpleKTCut KtCut
	newdef KtCut:MinKT 3.06

	# create the cuts object for the Underlying Event
	create ThePEG::Cuts UECuts
	# This should always be 2*MinKT!!
	newdef UECuts:MHatMin 6.12

	insert UECuts:OneCuts 0 KtCut

	######## subprocess ###########################################
	create ThePEG::SubProcessHandler FastQCD

	create Herwig::MEQCD2to2Fast MEQCD2to2Fast HwMEHadronFast.so
	insert FastQCD:MatrixElements 0 MEQCD2to2Fast
	cp /Herwig/Partons/PPExtractor /Herwig/Partons/MPIExtractor
	newdef /Herwig/Partons/MPIExtractor:FlatSHatY 1
	newdef FastQCD:PartonExtractor /Herwig/Partons/MPIExtractor

	# always use LO PDF for MPI and remnants
	newdef /Herwig/Partons/MPIExtractor:FirstPDF /Herwig/Partons/MPIPDF
	newdef /Herwig/Partons/MPIExtractor:SecondPDF /Herwig/Partons/MPIPDF

	######## MPI Handler ##########################################
	library JetCuts.so
	library HwMPI.so
	create Herwig::MPIHandler MPIHandler

	# set the subprocesses and corresponding cuts
	insert MPIHandler:SubProcessHandlers 0 FastQCD
	insert MPIHandler:Cuts 0 UECuts

	# energy extrapolation for the MinKT cut
	-newdef MPIHandler:EnergyExtrapolation Power
	-newdef MPIHandler:Power 0.322
	-newdef MPIHandler:pTmin0 2.55
	+newdef MPIHandler:EnergyExtrapolation PowerModified
	+newdef MPIHandler:Power 0.31
	+newdef MPIHandler:pTmin0 2.87
	+newdef MPIHandler:Offset 622.203*GeV
	newdef MPIHandler:ReferenceScale 7000.*GeV
	newdef MPIHandler:EEparamA 0.585*GeV
	newdef MPIHandler:EEparamB 37.45*GeV

	# MultiPeriph model as default.
	newdef /Herwig/Partons/RemnantDecayer:MultiPeriph Yes

	# Paramters for the ladders for soft MPIs
	-newdef /Herwig/Partons/RemnantDecayer:ladderMult 0.466
	-newdef /Herwig/Partons/RemnantDecayer:ladderbFactor 0.46
	+newdef /Herwig/Partons/RemnantDecayer:ladderMult 0.6
	+newdef /Herwig/Partons/RemnantDecayer:ladderbFactor 0.0
	# Set gaussian width of longitudinal momentum fraction fluctuation
	newdef /Herwig/Partons/RemnantDecayer:gaussWidth 0.03


	# Ratio of diffractive events in MinBias runs
	-newdef MPIHandler:DiffractiveRatio 0.133
	+newdef MPIHandler:DiffractiveRatio 0.25

	# The inverse hadron radius
	-newdef MPIHandler:InvRadius 1.58
	+newdef MPIHandler:InvRadius 1.1

	# Set the details of the soft model

	# Flag to decide whether additional soft interactions (i.e. pt < ptmin)
	# should be simulated at all

	newdef MPIHandler:softInt Yes

	# Flag to decide whether to use the 2 component model. In this model,
	# the hadron radius (or better: mean parton-parton separation) can be
	# different for hard and soft scatters. The radius for the soft scatters
	# is computed during the run startup so that the elastic t-slope, B_el,
	# coincides with current measurements (for CMenergy < 1.8 TeV) or
	# Donnachie-Landshoff pomeron fit (CMenergy > 1.8 TeV). If this model is
	# disabled one has to take care that this observable is correctly
	# described by adjusting InvRadius properly.

	newdef MPIHandler:twoComp Yes
	newdef MPIHandler:DLmode 2

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 3:05 PM (44 m, 59 s)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486836
Default Alt Text: (60 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions