No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	diff --git a/UnderlyingEvent/MPIHandler.cc b/UnderlyingEvent/MPIHandler.cc
	--- a/UnderlyingEvent/MPIHandler.cc
	+++ b/UnderlyingEvent/MPIHandler.cc
	@@ -1,849 +1,848 @@
	// -- 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){
	// 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(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(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);

	}

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:41 AM (10 h, 53 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111489
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions