No OneTemporary
Actions

Size

99 KB

Subscribers

None

View Options

	diff --git a/LesHouches/LesHouchesEventHandler.cc b/LesHouches/LesHouchesEventHandler.cc
	--- a/LesHouches/LesHouchesEventHandler.cc
	+++ b/LesHouches/LesHouchesEventHandler.cc
	@@ -1,650 +1,586 @@
	// -- C++ --
	//
	// LesHouchesEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG 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 LesHouchesEventHandler class.
	//

	#include "LesHouchesEventHandler.h"
	#include "LesHouchesReader.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Repository/Repository.h"
	#include "ThePEG/Handlers/LuminosityFunction.h"
	#include "ThePEG/Handlers/XComb.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/Utilities/LoopGuard.h"
	#include "ThePEG/EventRecord/Event.h"
	#include "ThePEG/EventRecord/Collision.h"
	#include "ThePEG/EventRecord/Step.h"
	#include "ThePEG/EventRecord/SubProcess.h"
	#include "ThePEG/Utilities/EnumIO.h"
	#include "ThePEG/Utilities/Maths.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace ThePEG;

	LesHouchesEventHandler::~LesHouchesEventHandler() {}

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

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

	void LesHouchesEventHandler::doinit() {

	EventHandler::doinit();

	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	readers()[i]->init();
	}

	ntries = 0;

	}

	void LesHouchesEventHandler::initialize() {
	- if ( lumiFnPtr() ) Repository::clog()
	- << "The LuminosityFunction '" << lumiFnPtr()->name()
	- << "' assigned to the LesHouchesEventHandler '" << name()
	- << "' will not be active in this run. Instead the incoming "
	- << "particles will be determined by the used LesHouchesReader objects.\n"
	- << Exception::warning;
	+ if ( lumiFnPtr() )
	+ Repository::clog()
	+ << "The LuminosityFunction '" << lumiFnPtr()->name()
	+ << "' assigned to the LesHouchesEventHandler '" << name()
	+ << "' will not be active in this run. Instead the incoming "
	+ << "particles will be determined by the used LesHouchesReader objects.\n"
	+ << Exception::warning;

	if ( readers().empty() )
	throw LesHouchesInitError()
	<< "No readers were defined for the LesHouchesEventHandler '"
	<< name() << "'" << Exception::warning;

	// Go through all the readers and collect information about cross
	// sections and processes.
	typedef map<int,tLesHouchesReaderPtr> ProcessMap;
	ProcessMap processes;
	PDPair incoming;
	Energy MaxEA = ZERO;
	Energy MaxEB = ZERO;
	- //temporary holder for weightnames: to check that they are the same between all readers
	- vector<vector<string> > weightnames_tmp;
	+ weightnames.clear();
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	LesHouchesReader & reader = *readers()[i];
	reader.initialize(*this);
	- weightnames_tmp.push_back(reader.optWeightsNamesFunc());
	+
	+ if ( i == 0 )
	+ weightnames = reader.optWeightsNamesFunc();
	+ else if ( reader.optWeightsNamesFunc() != weightnames )
	+ throw LesHouchesInitError()
	+ << "the optional weights names for the LesHouchesEventHandler do not match '"
	+ << name() << "'" << Exception::warning;
	+
	// Check that the incoming particles are consistent between the
	// readers.
	if ( !incoming.first ) {
	incoming.first = getParticleData(reader.heprup.IDBMUP.first);
	if ( !incoming.first )
	- Throw<LesHouchesInitError>()
	- << "Unknown beam PID " << reader.heprup.IDBMUP.first
	- << ". Have you created a matching BeamParticle object?"
	- << Exception::runerror;
	+ Throw<LesHouchesInitError>()
	+ << "Unknown beam PID " << reader.heprup.IDBMUP.first
	+ << ". Have you created a matching BeamParticle object?"
	+ << Exception::runerror;
	}
	if ( !incoming.second ) {
	incoming.second = getParticleData(reader.heprup.IDBMUP.second);
	if ( !incoming.second )
	- Throw<LesHouchesInitError>()
	- << "Unknown beam PID " << reader.heprup.IDBMUP.first
	- << ". Have you created a matching BeamParticle object?"
	- << Exception::runerror;
	+ Throw<LesHouchesInitError>()
	+ << "Unknown beam PID " << reader.heprup.IDBMUP.first
	+ << ". Have you created a matching BeamParticle object?"
	+ << Exception::runerror;
	}
	- if ( incoming.first->id() != reader.heprup.IDBMUP.first \|\|
	- incoming.second->id() != reader.heprup.IDBMUP.second )
	+ if ( incoming.first->id() != reader.heprup.IDBMUP.first \|\|
	+ incoming.second->id() != reader.heprup.IDBMUP.second )
	Repository::clog()
	- << "The different LesHouchesReader objects in the "
	- << "LesHouchesEventHandler '" << name() << "' have different "
	- << "types of colliding particles." << Exception::warning;
	+ << "The different LesHouchesReader objects in the "
	+ << "LesHouchesEventHandler '" << name() << "' have different "
	+ << "types of colliding particles." << Exception::warning;
	MaxEA = max(MaxEA, reader.heprup.EBMUP.first*GeV);
	MaxEB = max(MaxEB, reader.heprup.EBMUP.second*GeV);

	// Check that the weighting of the events in the different readers
	// is consistent with the ones requested for this event
	// handler. Also collect the sum of the maximum weights.
	if ( reader.negativeWeights() && weightOption() > 0 )
	throw LesHouchesInitError()
	- << "The reader '" << reader.name()
	- << "' contains negatively weighted events, "
	- << "which is not allowed for the LesHouchesEventHandler '"
	- << name() << "'." << Exception::warning;
	+ << "The reader '" << reader.name()
	+ << "' contains negatively weighted events, "
	+ << "which is not allowed for the LesHouchesEventHandler '"
	+ << name() << "'." << Exception::warning;

	// Check that we do not have the same process numbers in different
	// readers.
	for ( int ip = 0; ip < reader.heprup.NPRUP; ++ip ) {
	if ( reader.heprup.LPRUP[ip] ) {
	- ProcessMap::iterator pit = processes.find(reader.heprup.LPRUP[ip]);
	- if ( pit == processes.end() )
	- processes[reader.heprup.LPRUP[ip]] = readers()[i];
	- else if ( warnPNum ) {
	- Throw<LesHouchesPNumException>()
	- << "In the LesHouchesEventHandler '"
	- << name() << "', both the '" << pit->second->name() << "' and '"
	- << reader.name() << "' contains sub-process number " << pit->first
	- << ". This process may be double-counted in this run."
	- << Exception::warning;
	- }
	+ ProcessMap::iterator pit = processes.find(reader.heprup.LPRUP[ip]);
	+ if ( pit == processes.end() )
	+ processes[reader.heprup.LPRUP[ip]] = readers()[i];
	+ else if ( warnPNum ) {
	+ Throw<LesHouchesPNumException>()
	+ << "In the LesHouchesEventHandler '"
	+ << name() << "', both the '" << pit->second->name() << "' and '"
	+ << reader.name() << "' contains sub-process number " << pit->first
	+ << ". This process may be double-counted in this run."
	+ << Exception::warning;
	+ }
	}
	}
	selector().insert(reader.stats.maxXSec(), i);
	}
	- // Check that all the weightnames are equal amongst readers
	- if(readers().size() > 1) {
	- for ( int i = 1, N = readers().size(); i < N; ++i ) {
	- // already throw exception if the weightnames are not the same size
	- if(weightnames_tmp[0].size() != weightnames_tmp[i].size()) {
	- throw LesHouchesInitError()
	- << "the optional weights names for the LesHouchesEventHandler do not match '"
	- << name() << "'" << Exception::warning;
	- }
	- for(int j = 0; j < weightnames_tmp[0].size(); j++) {
	- if(weightnames_tmp[0][j] != weightnames_tmp[i][j]) {
	- throw LesHouchesInitError()
	- << "the optional weights names for the LesHouchesEventHandler do not match '"
	- << name() << "'" << Exception::warning;
	- }
	- }
	- }
	- }
	-
	- // Now set the weight names to that of the first reader
	- LesHouchesReader & reader = *readers()[0];
	- weightnames = reader.optWeightsNamesFunc();
	-
	-
	-
	stats.maxXSec(selector().sum());
	histStats.maxXSec(selector().sum());
	- for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
	- (it->second).maxXSec(selector().sum());
	+ for (map<string,OptWeight>::iterator it= opt.begin(); it!=opt.end(); ++it) {
	+ (it->second).stats.maxXSec(selector().sum());
	+ (it->second).histStats.maxXSec(selector().sum());
	}
	// Check that we have any cross section at all.
	if ( stats.maxXSec() <= ZERO )
	throw LesHouchesInitError()
	<< "The sum of the cross sections of the readers in the "
	<< "LesHouchesEventHandler '" << name()
	<< "' was zero." << Exception::warning;

	// We now create a LuminosityFunction object to inform others about
	// the energy of the beam.
	theIncoming = incoming;
	lumiFn(new_ptr(LuminosityFunction(MaxEA, MaxEB)));

	}

	void LesHouchesEventHandler::doinitrun() {
	EventHandler::doinitrun();
	stats.reset();
	histStats.reset();

	- vector<vector<string> > weightnames_tmp;
	+ weightnames.clear();

	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	+
	readers()[i]->initrun();
	LesHouchesReader & reader = *readers()[i];
	reader.initialize(*this);
	- weightnames_tmp.push_back(reader.optWeightsNamesFunc());
	+
	+ if ( i == 0 )
	+ weightnames = reader.optWeightsNamesFunc();
	+ else if ( reader.optWeightsNamesFunc() != weightnames )
	+ throw LesHouchesInitError()
	+ << "the optional weights names for the LesHouchesEventHandler do not match '"
	+ << name() << "'" << Exception::warning;
	+ }
	+
	+ for(unsigned int ww = 0; ww < weightnames.size(); ++ww) {
	+ const OptWeight optweight = { XSecStat(), ZERO, XSecStat() };
	+ opt.insert(make_pair(weightnames[ww], optweight));
	}

	- // Check that all the weightnames are equal amongst readers
	- if(readers().size() > 1) {
	- for ( int i = 1, N = readers().size(); i < N; ++i ) {
	- // already throw exception if the weightnames are not the same size
	- if(weightnames_tmp[0].size() != weightnames_tmp[i].size()) {
	- throw LesHouchesInitError()
	- << "the optional weights names for the LesHouchesEventHandler do not match '"
	- << name() << "'" << Exception::warning;
	- }
	- for(int j = 0; j < weightnames_tmp[0].size(); j++) {
	- if(weightnames_tmp[0][j] != weightnames_tmp[i][j]) {
	- throw LesHouchesInitError()
	- << "the optional weights names for the LesHouchesEventHandler do not match '"
	- << name() << "'" << Exception::warning;
	- }
	- }
	- }
	- }
	-
	- // Now set the weight names to that of the first reader
	- LesHouchesReader & reader = *readers()[0];
	- weightnames = reader.optWeightsNamesFunc();
	-
	- XSecStat* initxsecs = new XSecStat[weightnames.size()];
	- for(unsigned int ww = 0; ww < weightnames.size(); ww++){
	- initxsecs[ww].reset();
	- optstats.insert(std::make_pair<string,XSecStat>(weightnames[ww], initxsecs[ww]));
	- opthistStats.insert(std::make_pair<string,XSecStat>(weightnames[ww], initxsecs[ww]));
	- CrossSection initxs = 0.*picobarn;
	- optxs.insert(std::make_pair<string,CrossSection>(weightnames[ww], initxs));
	- }
	ntries = 0;

	}

	EventPtr LesHouchesEventHandler::generateEvent() {

	LoopGuard<EventLoopException,LesHouchesEventHandler>
	loopGuard(*this, maxLoop());

	while ( true ) {
	loopGuard();

	currentReader(readers()[selector().select(UseRandom::current())]);

	skipEvents();
	currentReader()->reset();

	double weight = currentReader()->getEvent();
	if ( weightOption() == unitweight && weight < 0.0 ) weight = 0.0;

	if ( weightOption() == unitweight \|\| weightOption() == unitnegweight ) {
	CrossSection newmax = selector().reweight(weight);
	if ( newmax > CrossSection() )
	- increaseMaxXSec(newmax);
	+ increaseMaxXSec(newmax);
	}
	-
	select(weight/currentReader()->preweight);
	histStats.select(weight);

	if ( !weighted() ) {
	if ( weightOption() == unitweight \|\| weightOption() == unitnegweight ) {
	- if ( !rndbool(abs(weight)) ) continue;
	- weight = Math::sign(1.0, weight);
	+ if ( !rndbool(abs(weight)) ) continue;
	+ weight = Math::sign(1.0, weight);
	}
	else if ( weight == 0.0 ) continue;
	} else if ( weight == 0.0 ) continue;

	accept();

	// Divide by the bias introduced by the preweights in the reader.
	-
	-
	weight /= currentReader()->preweight;

	-
	try {
	-
	theLastXComb = currentReader()->getXComb();

	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
	- generator()->currentEventNumber(), weight)));
	+ generator()->currentEventNumber(), weight)));
	currentEvent()->optionalWeights() = currentReader()->optionalEventWeights();

	-
	- //For debugging purposes: print optional weights here
	- /* for (map<string,double>::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){
	- std::cout << it->first << " => " << it->second << '\n';
	- }
	- cout << endl;*/
	-
	- //print npLO and npNLO (FxFx merging)
	- // cout << currentReader()->optionalEventnpLO() << "\t" << currentReader()->optionalEventnpNLO() << endl;
	-
	performCollision();

	if ( !currentCollision() ) throw Veto();

	return currentEvent();
	}
	catch (Veto) {
	reject(weight);
	}
	catch (Stop) {
	}
	catch (Exception &) {
	reject(weight);
	throw;
	}
	}
	}

	void LesHouchesEventHandler::skipEvents() {
	-
	+
	+ if ( weightOption() == 2 \|\| weightOption() == -2 ) return; //does it make sense to skip events if we are using varying weights?
	+
	// Don't do this for readers which seem to generate events on the fly.
	if ( currentReader()->active() \|\| currentReader()->NEvents() <= 0 ) return;
	-
	+
	// Estimate the fration of the total events available from
	// currentReader() which will be requested.
	double frac = currentReader()->stats.maxXSec()/stats.maxXSec();
	if ( stats.accepted() > 0 )
	frac *= double(stats.attempts())/double(stats.accepted());
	else
	frac *= double(stats.attempts() + 1);
	double xscan = generator()->N()*frac/currentReader()->NEvents();

	// Estimate the number of times we need to go through the events for
	// the currentReader(), and how many events on average we need to
	// skip for each attempted event to go through the file an integer
	// number of times.
	double nscan = ceil(xscan);
	double meanskip = nscan/xscan - 1.0;
	// Skip an average numer of steps with a Poissonian distribution.
	currentReader()->
	skip(UseRandom::rndPoisson(meanskip)%currentReader()->NEvents());
	}

	void LesHouchesEventHandler::select(double weight) {
	stats.select(weight);
	currentReader()->select(weight);
	- vector<double> w;
	- for (map<string,double>::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){
	- w.push_back(it->second);
	+ for (map<string,double>::const_iterator it = currentReader()->optionalEventWeights().begin();
	+ it != currentReader()->optionalEventWeights().end();
	+ ++it) {
	+ const double w = it->second;
	+ OptWeight & o = opt[it->first];
	+ o.histStats.select(w);
	+ o.stats.select(w);
	}
	- int ii = 0;
	- for (map<string,XSecStat>::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){
	- (it->second).select(w[ii]);
	- ii++;
	- }
	- ii = 0;
	- for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
	- (it->second).select(w[ii]);
	- ii++;
	- }
	-
	}

	tCollPtr LesHouchesEventHandler::performCollision() {
	lastExtractor()->select(lastXCombPtr());
	if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr());
	currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this)));
	if ( currentEvent() ) currentEvent()->addCollision(currentCollision());
	currentStep(new_ptr(Step(currentCollision(), this)));
	currentCollision()->addStep(currentStep());
	currentStep()->addSubProcess(currentReader()->getSubProcess());
	lastExtractor()->constructRemnants(lastXCombPtr()->partonBinInstances(),
	- subProcess(), currentStep());
	+ subProcess(), currentStep());

	if ( !currentReader()->cuts().passCuts(*currentCollision()) ) throw Veto();

	initGroups();
	if ( ThePEG_DEBUG_ITEM(1) ) {
	if ( currentEvent() )
	generator()->logfile() << *currentEvent();
	else
	generator()->logfile() << *currentCollision();
	}
	return continueCollision();
	}

	EventPtr LesHouchesEventHandler::continueEvent() {
	try {
	continueCollision();
	}
	catch (Veto) {
	reject(currentEvent()->weight());
	}
	catch (Stop) {
	}
	catch (Exception &) {
	reject(currentEvent()->weight());
	throw;
	}
	return currentEvent();
	}

	void LesHouchesEventHandler::dofinish() {
	EventHandler::dofinish();
	- if ( selector().compensating() ) generator()->log()
	- << "Warning: The run was ended while the LesHouchesEventHandler '"
	- << name() << "' was still trying to compensate for weights larger than 1. "
	- << "The cross section estimates may therefore be statistically "
	- << "inaccurate." << endl;
	+ if ( selector().compensating() )
	+ generator()->log()
	+ << "Warning: The run was ended while the LesHouchesEventHandler '"
	+ << name() << "' was still trying to compensate for weights larger than 1. "
	+ << "The cross section estimates may therefore be statistically "
	+ << "inaccurate." << endl;
	}

	void LesHouchesEventHandler::statistics(ostream & os) const {
	if ( statLevel() == 0 ) return;

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

	if ( stats.accepted() <= 0 ) {
	os << line << "No events generated by event handler '" << name() << "'."
	<< endl;
	- return;
	+ return;
	}

	os << line << "Statistics for Les Houches event handler \'" << name() << "\':\n"
	<< " "
	<< "generated number of Cross-section\n"
	<< " "
	<< " events attempts (nb)\n";

	os << line << "Total:" << setw(42) << stats.accepted() << setw(13)
	<< stats.attempts() << setw(17)
	<< ouniterr(stats.xSec(), stats.xSecErr(), nanobarn) << endl
	<< line;

	if ( statLevel() == 1 ) return;

	if ( statLevel() == 2 ) {

	os << "Per Les Houches Reader breakdown:\n";
	-
	+
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	LesHouchesReader & reader = *readers()[i];
	string n = reader.name();
	n.resize(37, ' ');
	os << n << setw(11) << reader.stats.accepted() << setw(13)
	- << reader.stats.attempts() << setw(17)
	- << ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn)
	- << endl;
	+ << reader.stats.attempts() << setw(17)
	+ << ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn)
	+ << endl;
	}
	os << line;
	} else {

	os << "Per Les Houches Reader (and process #) breakdown:\n";
	-
	+
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	LesHouchesReader & reader = *readers()[i];
	string n = reader.name() + " (all)";
	n.resize(37, ' ');
	os << n << setw(11) << reader.stats.accepted() << setw(13)
	- << reader.stats.attempts() << setw(17)
	- << ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn)
	- << endl;
	+ << reader.stats.attempts() << setw(17)
	+ << ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn)
	+ << endl;
	CrossSection xsectot = reader.stats.xSec();
	if ( xsectot != ZERO ) xsectot /= reader.stats.sumWeights();
	typedef LesHouchesReader::StatMap::const_iterator const_iterator;
	for ( const_iterator i = reader.statmap.begin();
	- i != reader.statmap.end(); ++i ) {
	- ostringstream ss;
	- ss << reader.name() << " (" << i->first << ")";
	- string n = ss.str();
	- n.resize(37, ' ');
	- os << n << setw(11) << i->second.accepted() << setw(13)
	- << i->second.attempts() << setw(17)
	- << ouniterr(i->second.sumWeights()*xsectot,
	- sqrt(i->second.sumWeights2())*xsectot, nanobarn) << endl;
	+ i != reader.statmap.end(); ++i ) {
	+ ostringstream ss;
	+ ss << reader.name() << " (" << i->first << ")";
	+ string n = ss.str();
	+ n.resize(37, ' ');
	+ os << n << setw(11) << i->second.accepted() << setw(13)
	+ << i->second.attempts() << setw(17)
	+ << ouniterr(i->second.sumWeights()*xsectot,
	+ sqrt(i->second.sumWeights2())*xsectot, nanobarn) << endl;
	}
	os << line;
	}
	}

	string warn = "Warning: Result may be statistically incorrect since\n"
	" the following LesHouchesReaders were oversampled:\n";
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	LesHouchesReader & reader = *readers()[i];
	if ( reader.NEvents() > 0 && reader.stats.attempts() > reader.NEvents() ) {
	os << warn;
	warn = "";
	os << "'" << reader.name() << "' (by a factor "
	- << double(reader.stats.attempts())/double(reader.NEvents())
	- << ")" << endl;
	+ << double(reader.stats.attempts())/double(reader.NEvents())
	+ << ")" << endl;
	}
	}
	}


	void LesHouchesEventHandler::increaseMaxXSec(CrossSection maxxsec) {
	stats.maxXSec(selector().sum());
	histStats.maxXSec(selector().sum());
	currentReader()->increaseMaxXSec(maxxsec);
	}

	void LesHouchesEventHandler::accept() {
	ntries++;
	stats.accept();
	histStats.accept();
	currentReader()->accept();
	- for (map<string,XSecStat>::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){
	- (it->second).accept();
	- }
	- for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
	- (it->second).accept();
	+ for (map<string,OptWeight>::iterator it = opt.begin(); it!=opt.end(); ++it) {
	+ OptWeight & tmp = it->second;
	+ tmp.histStats.accept();
	+ tmp.stats.accept();
	}
	}

	void LesHouchesEventHandler::reject(double w) {
	ntries++;
	stats.reject(w);
	histStats.reject(w);
	currentReader()->reject(w);
	- vector<double> wv;
	- for (map<string,double>::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){
	- wv.push_back(it->second);
	- }
	- int ii = 0;
	- for (map<string,XSecStat>::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){
	- (it->second).reject(wv[ii]);
	- ii++;
	- }
	- ii = 0;
	- for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
	- (it->second).reject(wv[ii]);
	- ii++;
	+
	+ for (map<string,double>::const_iterator it = currentReader()->optionalEventWeights().begin();
	+ it != currentReader()->optionalEventWeights().end();
	+ ++it) {
	+ const double w = it->second;
	+ OptWeight & o = opt[it->first];
	+ o.histStats.reject(w);
	+ o.stats.reject(w);
	}
	}

	-
	-
	-map<string,CrossSection> LesHouchesEventHandler::optintegratedXSecMap() const {
	- map<string,CrossSection> result;
	- for (map<string,XSecStat>::const_iterator it= optstats.begin(); it!=optstats.end(); ++it){
	- result[it->first] = (it->second.sumWeights()/it->second.attempts()) * picobarn;
	+const map<string,CrossSection> & LesHouchesEventHandler::optintegratedXSecMap() const {
	+ static map<string,CrossSection> result;
	+ result.clear();
	+ for ( map<string,OptWeight>::const_iterator it= opt.begin(); it!=opt.end(); ++it ) {
	+ result[it->first] = ( it->second.stats.sumWeights() / it->second.stats.attempts() ) * picobarn;
	}
	-
	return result;
	}

	-
	CrossSection LesHouchesEventHandler::histogramScale() const {
	return histStats.xSec()/histStats.sumWeights();
	}

	CrossSection LesHouchesEventHandler::integratedXSec() const {
	return histStats.xSec();
	}

	CrossSection LesHouchesEventHandler::integratedXSecErr() const {
	return histStats.xSecErr();
	}

	int LesHouchesEventHandler::ntriesinternal() const {
	return stats.attempts();
	}

	void LesHouchesEventHandler::persistentOutput(PersistentOStream & os) const {
	os << stats << histStats << theReaders << theSelector
	<< oenum(theWeightOption) << theUnitTolerance << theCurrentReader << warnPNum;
	}

	void LesHouchesEventHandler::persistentInput(PersistentIStream & is, int) {
	is >> stats >> histStats >> theReaders >> theSelector
	>> ienum(theWeightOption) >> theUnitTolerance >> theCurrentReader >> warnPNum;
	}

	ClassDescription<LesHouchesEventHandler>
	LesHouchesEventHandler::initLesHouchesEventHandler;
	// Definition of the static class description member.

	void LesHouchesEventHandler::setUnitTolerance(double x) {
	theUnitTolerance = x;
	selector().tolerance(unitTolerance());
	}

	void LesHouchesEventHandler::Init() {

	static ClassDocumentation<LesHouchesEventHandler> documentation
	("This is the main class administrating the selection of hard "
	"subprocesses from a set of ThePEG::LesHouchesReader objects.");


	static RefVector<LesHouchesEventHandler,LesHouchesReader>
	interfaceLesHouchesReaders
	("LesHouchesReaders",
	"Objects capable of reading events from an event file or an "
	"external matrix element generator.",
	&LesHouchesEventHandler::theReaders, -1, false, false, true, false, false);

	static Switch<LesHouchesEventHandler,WeightOpt> interfaceWeightOption
	("WeightOption",
	"The different ways to weight events in the Les Houches event handler. "
	"Whether weighted or not and whether or not negative weights are allowed.",
	&LesHouchesEventHandler::theWeightOption, unitweight, true, false);
	static SwitchOption interfaceWeightOptionUnitWeight
	(interfaceWeightOption,
	"UnitWeight",
	"All events have unit weight.",
	unitweight);
	static SwitchOption interfaceWeightOptionNegUnitWeight
	(interfaceWeightOption,
	"NegUnitWeight",
	"All events have weight +1 or maybe -1.",
	unitnegweight);
	static SwitchOption interfaceWeightOptionVarWeight
	(interfaceWeightOption,
	"VarWeight",
	"Events may have varying but positive weights.",
	varweight);
	static SwitchOption interfaceWeightOptionVarNegWeight
	(interfaceWeightOption,
	"VarNegWeight",
	"Events may have varying weights, both positive and negative.",
	varnegweight);

	static Switch<LesHouchesEventHandler,bool> interfaceWarnPNum
	("WarnPNum",
	"Warn if the same process number is used in more than one "
	"LesHouchesReader.",
	&LesHouchesEventHandler::warnPNum, true, true, false);
	static SwitchOption interfaceWarnPNumWarning
	(interfaceWarnPNum,
	"Warning",
	"Give a warning message.",
	true);
	static SwitchOption interfaceWarnPNumNoWarning
	(interfaceWarnPNum,
	"NoWarning",
	"Don't give a warning message.",
	false);

	static Parameter<LesHouchesEventHandler,double> interfaceUnitTolerance
	("UnitTolerance",
	"If the <interface>WeightOption</interface> is set to unit weight, do not start compensating unless the a weight is found to be this much larger than unity.",
	&LesHouchesEventHandler::theUnitTolerance, 1.0e-6, 0.0, 0,
	true, false, Interface::lowerlim,
	&LesHouchesEventHandler::setUnitTolerance,
	(double(LesHouchesEventHandler::*)()const)(0),
	(double(LesHouchesEventHandler::*)()const)(0),
	(double(LesHouchesEventHandler::*)()const)(0),
	(double(LesHouchesEventHandler::*)()const)(0));

	interfaceLesHouchesReaders.rank(10);
	interfaceWeightOption.rank(9);

	}

	diff --git a/LesHouches/LesHouchesEventHandler.h b/LesHouches/LesHouchesEventHandler.h
	--- a/LesHouches/LesHouchesEventHandler.h
	+++ b/LesHouches/LesHouchesEventHandler.h
	@@ -1,474 +1,470 @@
	// -- C++ --
	//
	// LesHouchesEventHandler.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef THEPEG_LesHouchesEventHandler_H
	#define THEPEG_LesHouchesEventHandler_H
	//
	// This is the declaration of the LesHouchesEventHandler class.
	//

	#include "ThePEG/Handlers/EventHandler.h"
	#include "LesHouchesEventHandler.fh"
	#include "LesHouchesReader.fh"
	-//#include "LesHouchesAnalysis.h"
	#include "ThePEG/Utilities/CompSelector.h"
	#include "ThePEG/Utilities/XSecStat.h"

	namespace ThePEG {


	/**
	* The LesHouchesEventHandler inherits from the general EventHandler
	* class and administers the reading of events generated by external
	* matrix element generator programs according to the Les Houches
	* accord.
	*
	* The class has a list of <code>LesHouchesReader</code>s which
	* typically are connected to files with event data produced by
	* external matrix element generator programs. When an event is
	* requested by LesHouchesEventHandler, one of the readers are chosen,
	* an event is read in and then passed to the different
	* <code>StepHandler</code> defined in the underlying
	* EventHandler class.
	*
	* @see \ref LesHouchesEventHandlerInterfaces "The interfaces"
	* defined for LesHouchesEventHandler.
	*/
	class LesHouchesEventHandler: public EventHandler {

	public:

	/**
	* A vector of LesHouchesReader objects.
	*/
	typedef vector<LesHouchesReaderPtr> ReaderVector;

	/**
	* A selector of readers.
	*/
	typedef CompSelector<int,CrossSection> ReaderSelector;

	/**
	* Enumerate the weighting options.
	*/
	enum WeightOpt {
	unitweight = 1, /*< All events have unit weight. /
	unitnegweight = -1, /*< All events have wight +/- 1. /
	varweight = 2, /*< Varying positive weights. /
	varnegweight = -2 /*< Varying positive or negative weights. /
	};

	- friend class LesHouchesHandler;
	-
	-
	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	LesHouchesEventHandler()
	: theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true)
	{
	selector().tolerance(unitTolerance());
	}

	/**
	* The destructor.
	*/
	virtual ~LesHouchesEventHandler();
	//@}

	public:

	/** @name Initialization and finalization functions. */
	//@{
	/**
	* Initialize this event handler and all related objects needed to
	* generate events.
	*/
	virtual void initialize();

	/**
	* Write out accumulated statistics about intergrated cross sections
	* and stuff.
	*/
	virtual void statistics(ostream &) const;

	/**
	* Histogram scale. A histogram bin which has been filled with the
	* weights associated with the Event objects should be scaled by
	* this factor to give the correct cross section.
	*/
	virtual CrossSection histogramScale() const;

	/**
	* The estimated total integrated cross section of the processes
	* generated in this run.
	* @return 0 if no integrated cross section could be estimated.
	*/
	virtual CrossSection integratedXSec() const;

	/**
	* The number of attempts inside the statistics object
	*/

	virtual int ntriesinternal() const;

	/**
	* The estimated error in the total integrated cross section of the
	* processes generated in this run.
	* @return 0 if no integrated cross section error could be estimated.
	*/
	virtual CrossSection integratedXSecErr() const;

	/**
	* Map to aid the calculation of the optional weights' integrated cross section
	*/

	- virtual map<string,CrossSection> optintegratedXSecMap() const;
	+ virtual const map<string,CrossSection> & optintegratedXSecMap() const;

	//@}

	/** @name Functions used for the actual generation */
	//@{
	/**
	* Generate an event.
	*/
	virtual EventPtr generateEvent();

	/**
	* Create the Event and Collision objects. Used by the
	* generateEvent() function.
	*/
	virtual tCollPtr performCollision();

	/**
	* Continue generating an event if the generation has been stopped
	* before finishing.
	*/
	virtual EventPtr continueEvent();
	//@}

	/** @name Functions to manipulate statistics. */
	//@{
	/**
	* An event has been selected. Signal that an event has been
	* selected with the given \a weight. If unit weights are requested,
	* the event will be accepted with that weight. This also takes care
	* of the statistics collection of the selected reader object.
	*/
	void select(double weight);

	/**
	* Accept the current event, taking care of the statistics
	* collection of the corresponding reader objects.
	*/
	void accept();

	/**
	* Reject the current event, taking care of the statistics
	* collection of the corresponding reader objects.
	*/
	void reject(double weight);

	/**
	* Increase the overestimated cross section for the selected reader.
	*/
	void increaseMaxXSec(CrossSection maxxsec);

	/**
	* Skip some events. To ensure a reader file is scanned an even
	* number of times, skip a number of events for the selected reader.
	*/
	void skipEvents();

	//@}

	/** @name Simple access functions. */
	//@{
	/**
	* The way weights are to be treated.
	*/
	WeightOpt weightOption() const { return theWeightOption; }

	/**
	* If the weight option is set to unit weight, do not start
	* compensating unless the weight is this much larger than unity.
	*/
	double unitTolerance() const { return theUnitTolerance; }

	/**
	* Access the list of readers.
	*/
	const ReaderVector & readers() const { return theReaders; }

	/**
	* The selector to choose readers according to their overestimated
	* cross section.
	*/
	const ReaderSelector & selector() const { return theSelector; }

	/**
	* The currently selected reader object.
	*/
	tLesHouchesReaderPtr currentReader() const { return theCurrentReader; }


	/**
	* Set the currently selected reader object.
	*/
	void currentReader(tLesHouchesReaderPtr x) { theCurrentReader = x; }

	//@}

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	/**
	* The currently selected reader object.
	*/
	tLesHouchesReaderPtr theCurrentReader;

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

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Initialize this object. Called in the run phase just before
	* a run begins.
	*/
	virtual void doinitrun();

	/**
	* Finalize this object. Called in the run phase just after a
	* run has ended. Used eg. to write out statistics.
	*/
	virtual void dofinish();
	//@}

	protected:

	/**
	* Access the list of readers.
	*/
	ReaderVector & readers() { return theReaders; }

	/**
	* The selector to choose readers according to their overestimated
	* cross section.
	*/
	ReaderSelector & selector() { return theSelector; }

	/**
	* Helper function for the interface;
	*/
	void setUnitTolerance(double);

	/**
	* Collect statistics for this event handler.
	*/
	XSecStat stats;

	- /*
	- * Collect statistics for the optional weights
	- */
	- map<string,XSecStat> optstats;
	-
	- /*
	- * Calculate the cross section for the optional weights
	- */
	-
	- map<string,CrossSection> optxs;
	-
	- /*
	- * Counter for the number of tries for the purpose of statistics
	- */
	-
	- int ntries;
	-
	- /*
	- * Return the optional weights' statistics
	- */
	-
	- map<string,XSecStat> OptStatsFunc() { return optstats; }
	-
	- /*
	- * Return the optional weights' cross sections
	- */
	-
	- map<string,CrossSection> OptXsFunc() { return optxs; }
	-
	-
	/**
	* Collect statistics for this event handler. To be used for
	* histogram scaling.
	*/
	XSecStat histStats;

	- /**
	- * Collect statistics for this event handler's optional weights. To be used for
	- * histogram scaling.
	- */
	-
	- map<string,XSecStat> opthistStats;
	-
	- /*
	+ /**
	* The weight identifiers for the events
	*/

	vector<string> weightnames;

	+ /**
	+ * Collect statistics for this event handler's optional weights.
	+ */
	+ struct OptWeight {
	+
	+ /**
	+ * Collect statistics for the optional weights
	+ */
	+ XSecStat stats;
	+
	+ /**
	+ * Calculate the cross section for the optional weights
	+ */
	+ CrossSection xs;
	+
	+ /**
	+ * Collect statistics for this event handler's optional weights. To be used for
	+ * histogram scaling.
	+ */
	+ XSecStat histStats;
	+ };
	+
	+ /** Map statistics to weight name strings */
	+ map <string,OptWeight> opt;
	+
	+ /**
	+ * Counter for the number of tries for the purpose of statistics
	+ */
	+
	+ int ntries;
	+
	+ /**
	+ * Return the optional weights' statistics
	+ */
	+
	+ const map<string,OptWeight> & optWeights() const { return opt; }
	+
	private:

	/**
	* The list of readers.
	*/
	ReaderVector theReaders;

	/**
	* The selector to choose readers according to their overestimated
	* cross section.
	*/
	ReaderSelector theSelector;

	/**
	* The way weights are to be treated.
	*/
	WeightOpt theWeightOption;

	/**
	* If the weight option is set to unit weight, do not start
	* compensating unless the weight is this much larger than unity.
	*/
	double theUnitTolerance;


	/**
	* Warn if the same process number is used in more than one
	* LesHouchesReader.
	*/
	bool warnPNum;



	public:

	/** @cond EXCEPTIONCLASSES */
	/**
	* Exception class used if no readers were assigned.
	*/
	class LesHouchesInitError: public InitException {};

	/**
	* Exception class used if the same process number is used by more
	* than ne reader.
	*/
	class LesHouchesPNumException: public InitException {};
	/** @endcond */

	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<LesHouchesEventHandler> initLesHouchesEventHandler;

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

	};

	}

	// CLASSDOC OFF

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of LesHouchesEventHandler. */
	template <>
	struct BaseClassTrait<LesHouchesEventHandler,1> {
	/** Typedef of the first base class of LesHouchesEventHandler. */
	typedef EventHandler NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the LesHouchesEventHandler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<LesHouchesEventHandler>
	: public ClassTraitsBase<LesHouchesEventHandler> {
	/** Return a platform-independent class name */
	static string className() { return "ThePEG::LesHouchesEventHandler"; }
	/** Return the name of the shared library be loaded to get access to
	* the LesHouchesEventHandler class and every other class it uses
	* (except the base class). */
	static string library() { return "LesHouches.so"; }
	};

	/** @endcond */

	}

	-#endif /* HERWIG_LesHouchesEventHandler_H */
	+#endif /* THEPEG_LesHouchesEventHandler_H */
	diff --git a/LesHouches/LesHouchesReader.cc b/LesHouches/LesHouchesReader.cc
	--- a/LesHouches/LesHouchesReader.cc
	+++ b/LesHouches/LesHouchesReader.cc
	@@ -1,1566 +1,1583 @@
	// -- C++ --
	//
	// LesHouchesReader.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG 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 LesHouchesReader class.
	//

	#include "LesHouchesReader.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Command.h"
	#include "config.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/PDF/NoPDF.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/EventRecord/TmpTransform.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/Handlers/XComb.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "LesHouchesEventHandler.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/HoldFlag.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"

	using namespace ThePEG;

	LesHouchesReader::LesHouchesReader(bool active)
	: theNEvents(0), position(0), reopened(0), theMaxScan(-1), scanning(false),
	isActive(active), theCacheFileName(""), doCutEarly(true),
	preweight(1.0), reweightPDF(false), doInitPDFs(false),
	theMaxMultCKKW(0), theMinMultCKKW(0), lastweight(1.0), maxFactor(1.0), optionalnpLO(0), optionalnpNLO(0),
	weightScale(1.0*picobarn), skipping(false), theMomentumTreatment(0),
	useWeightWarnings(true),theReOpenAllowed(true), theIncludeSpin(true) {}

	LesHouchesReader::LesHouchesReader(const LesHouchesReader & x)
	: HandlerBase(x), LastXCombInfo<>(x), heprup(x.heprup), hepeup(x.hepeup),
	inData(x.inData), inPDF(x.inPDF), outPDF(x.outPDF),
	thePartonExtractor(x.thePartonExtractor), thePartonBins(x.thePartonBins),
	theXCombs(x.theXCombs), theCuts(x.theCuts),
	theNEvents(x.theNEvents), position(x.position), reopened(x.reopened),
	theMaxScan(x.theMaxScan), scanning(false),
	isActive(x.isActive),
	theCacheFileName(x.theCacheFileName), doCutEarly(x.doCutEarly),
	stats(x.stats), statmap(x.statmap),
	thePartonBinInstances(x.thePartonBinInstances),
	reweights(x.reweights), preweights(x.preweights),
	preweight(x.preweight), reweightPDF(x.reweightPDF),
	doInitPDFs(x.doInitPDFs),
	theMaxMultCKKW(x.theMaxMultCKKW), theMinMultCKKW(x.theMinMultCKKW),
	lastweight(x.lastweight), maxFactor(x.maxFactor),
	weightScale(x.weightScale), xSecWeights(x.xSecWeights),
	maxWeights(x.maxWeights), skipping(x.skipping),
	theMomentumTreatment(x.theMomentumTreatment),
	useWeightWarnings(x.useWeightWarnings),
	theReOpenAllowed(x.theReOpenAllowed),
	theIncludeSpin(x.theIncludeSpin) {}

	LesHouchesReader::~LesHouchesReader() {}

	void LesHouchesReader::doinitrun() {
	HandlerBase::doinitrun();
	stats.reset();
	for ( StatMap::iterator i = statmap.begin(); i != statmap.end(); ++i )
	i->second.reset();
	open();
	if ( cacheFileName().length() ) openReadCacheFile();
	position = 0;
	reopened = 0;
	}

	bool LesHouchesReader::preInitialize() const {
	if ( HandlerBase::preInitialize() ) return true;
	if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) return true;
	return false;
	}

	void LesHouchesReader::doinit() {
	HandlerBase::doinit();
	open();
	close();
	if ( !heprup.IDBMUP.first \|\| !heprup.IDBMUP.second )
	Throw<LesHouchesInitError>()
	<< "No information about incoming particles were found in "
	<< "LesHouchesReader '" << name() << "'." << Exception::warning;
	inData = make_pair(getParticleData(heprup.IDBMUP.first),
	getParticleData(heprup.IDBMUP.second));
	if ( heprup.EBMUP.first <= 0.0 \|\| heprup.EBMUP.second <= 0.0 )
	Throw<LesHouchesInitError>()
	<< "No information about the energy of incoming particles were found in "
	<< "LesHouchesReader '" << name() << "'." << Exception::warning;

	if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) {
	initPDFs();
	if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
	<< "LesHouchesReader '" << name()
	<< "' could not create PDFBase objects in pre-initialization."
	<< Exception::warning;
	}
	else if ( !inPDF.first \|\| !inPDF.second ) Throw<LesHouchesInitError>()
	<< "No information about the PDFs of incoming particles were found in "
	<< "LesHouchesReader '" << name() << "'." << Exception::warning;
	}

	void LesHouchesReader::initPDFs() {
	if ( inPDF.first && inPDF.second ) return;

	string remhname;
	if ( heprup.PDFSUP.first && !inPDF.first) {
	inPDF.first = dynamic_ptr_cast<PDFPtr>
	(generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFA",
	"ThePEGLHAPDF.so"));
	if ( !inPDF.first ) {
	Throw<InitException>()
	<< "LesHouchesReader '" << name() << "' could not use information "
	<< "about the PDFs used because the LHAPDF library was not properly "
	"defined." << Exception::warning;
	return;
	}
	remhname = fullName() + "/DummyRemH";
	generator()->preinitCreate("ThePEG::NoRemnants", remhname);
	generator()->preinitInterface(inPDF.first, "RemnantHandler",
	"set", remhname);
	if ( heprup.PDFGUP.first > 0 && heprup.PDFGUP.first < 10 ) {
	ostringstream os;
	os << heprup.PDFGUP.first << " " << heprup.PDFSUP.first;
	generator()->preinitInterface(inPDF.first, "PDFLIBNumbers",
	"set", os.str());
	} else {
	ostringstream os;
	os << heprup.PDFGUP.first*1000 + heprup.PDFSUP.first;
	generator()->preinitInterface(inPDF.first, "PDFNumber",
	"set", os.str());
	}
	generator()->preinitInterface(inPDF.first, "RangeException",
	"newdef", "Freeze");
	}

	if ( heprup.PDFSUP.second && !inPDF.second) {
	inPDF.second = dynamic_ptr_cast<PDFPtr>
	(generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFB",
	"ThePEGLHAPDF.so"));
	if ( !inPDF.second ) {
	Throw<InitException>()
	<< "LesHouchesReader '" << name() << "' could not use information "
	<< "about the PDFs used because the LHAPDF library was not properly "
	"defined." << Exception::warning;
	return;
	}
	if ( remhname == "" ) {
	remhname = fullName() + "/DummyRemH";
	generator()->preinitCreate("ThePEG::NoRemnants", remhname);
	}
	generator()->preinitInterface(inPDF.second, "RemnantHandler",
	"set", remhname);
	if ( heprup.PDFGUP.second > 0 && heprup.PDFGUP.second < 10 ) {
	ostringstream os;
	os << heprup.PDFGUP.second << " " << heprup.PDFSUP.second;
	generator()->preinitInterface(inPDF.second, "PDFLIBNumbers",
	"set", os.str());
	} else {
	ostringstream os;
	os << heprup.PDFGUP.second*1000 + heprup.PDFSUP.second;
	generator()->preinitInterface(inPDF.second, "PDFNumber",
	"set", os.str());
	}
	generator()->preinitInterface(inPDF.second, "RangeException",
	"newdef", "Freeze");
	}

	if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
	<< "LesHouchesReader '" << name()
	<< "' could not find information about the PDFs used."
	<< Exception::warning;
	}

	void LesHouchesReader::initialize(LesHouchesEventHandler & eh) {
	Energy2 Smax = ZERO;
	double Y = 0.0;
	if ( !theCuts ) {
	theCuts = eh.cuts();
	if ( !theCuts ) Throw<LesHouchesInitError>()
	<< "No Cuts object was assigned to the LesHouchesReader '"
	<< name() << "' nor was one\nassigned to the controlling "
	<< "LesHouchesEventHandler '" << eh.name() << "'.\nAt least one of them "
	<< "needs to have a Cuts object." << Exception::runerror;

	Smax = cuts().SMax();
	Y = cuts().Y();
	}

	theCKKW = eh.CKKWHandler();

	if ( !partonExtractor() ) {
	thePartonExtractor = eh.partonExtractor();
	if ( !partonExtractor() ) Throw<LesHouchesInitError>()
	<< "No PartonExtractor object was assigned to the LesHouchesReader '"
	<< name() << "' nor was one\nassigned to the controlling "
	<< "LesHouchesEventHandler '" << eh.name() << "'.\nAt least one of them "
	<< "needs to have a PartonExtractor object." << Exception::runerror;
	}
	open();

	Energy emax = 2.0sqrt(heprup.EBMUP.firstheprup.EBMUP.second)*GeV;
	theCuts->initialize(sqr(emax),
	0.5*log(heprup.EBMUP.first/heprup.EBMUP.second));
	if ( Smax > ZERO && ( Smax != cuts().SMax() \|\| Y != cuts().Y() ) )
	Throw<LesHouchesInitError>()
	<< "The LesHouchesReader '" << name() << "' uses the same Cuts object "
	<< "as another LesHouchesReader which has not got the same energies of "
	<< "the colliding particles. For the generation to work properly "
	<< "different LesHouchesReader object with different colliding particles "
	<< "must be assigned different (although possibly identical) Cuts "
	<< "objects." << Exception::warning;

	thePartonBins = partonExtractor()->getPartons(emax, inData, cuts());
	for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
	theXCombs[partonBins()[i]] =
	new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(),
	partonBins()[i], theCuts));
	partonExtractor()->nDims(partonBins()[i]);
	}
	outPDF = make_pair(partonExtractor()->getPDF(inData.first),
	partonExtractor()->getPDF(inData.second));

	// SP: re-interpret 3/4 -> 1/2; see discussion w/ Leif and Keith at LH 2013
	-
	- if ( abs(heprup.IDWTUP) == 3 )
	+
	+ /* if ( abs(heprup.IDWTUP) == 3 )
	heprup.IDWTUP = heprup.IDWTUP < 0 ? -1 : 1;

	if ( abs(heprup.IDWTUP) == 4 )
	- heprup.IDWTUP = heprup.IDWTUP < 0 ? -2 : 2;
	+ heprup.IDWTUP = heprup.IDWTUP < 0 ? -2 : 2;*/

	close();

	if ( !heprup.IDWTUP && useWeightWarnings )
	Throw<LesHouchesInitError>()
	<< "No information about the weighting scheme was found. The events "
	<< "produced by LesHouchesReader " << name()
	<< " may not be sampled correctly." << Exception::warning;

	if ( abs(heprup.IDWTUP) > 1 && !eh.weighted() && useWeightWarnings )
	Throw<LesHouchesInitError>()
	<< "LesHouchesReader " << name() << " has the IDWTUP flag set to "
	<< heprup.IDWTUP << " which is not supported by this reader, the "
	<< "produced events may not be sampled correctly. It is up to "
	<< "sub-classes of LesHouchesReader to correctly convert to match IDWTUP "
	<< "+/- 1. Will try to make intelligent guesses to get "
	<< "correct statistics.\nIn most cases this should be sufficient. "
	<< "Unset <interface>WeightWarnings</interface> to avoid this message,"
	<< "or set <interface>Weighted</interface> to on."
	<< Exception::warning;

	if ( heprup.IDWTUP != eh.weightOption() && abs(heprup.IDWTUP) < 3 &&
	useWeightWarnings )
	Throw<LesHouchesInitError>()
	<< "LesHouchesReader " << name() << " has the IDWTUP flag set to "
	<< heprup.IDWTUP
	<< ", which does not correspond\nto the weight option "
	<< eh.weightOption() << " set in "
	<< "the LesHouchesEventHandler " << eh.name() << ".\n\n"
	<< "Use the following handler setting instead:\n"
	<< " set " << eh.name() << ":WeightOption " << heprup.IDWTUP
	<< "\nWill try to make intelligent guesses to get "
	<< "correct statistics. In most cases this should be sufficient. "
	<< "Unset <interface>WeightWarnings</interface> to avoid this message"
	<< Exception::warning;

	scan();
	initStat();
	}


	long LesHouchesReader::scan() {

	open();

	// Shall we write the events to a cache file for fast reading? If so
	// we write to a temporary file if the caches events should be
	// randomized.
	if ( cacheFileName().length() ) openWriteCacheFile();

	// Keep track of the number of events scanned.
	long neve = 0;
	long cuteve = 0;
	bool negw = false;

	// If the open() has not already gotten information about subprocesses
	// and cross sections we have to scan through the events.
	if ( !heprup.NPRUP \|\| cacheFile() \|\| abs(heprup.IDWTUP) != 1 ) { // why scan if IDWTUP != 1?

	HoldFlag<> isScanning(scanning);

	double oldsum = 0.0;
	vector<int> lprup;
	vector<double> newmax;
	vector<long> oldeve;
	vector<long> neweve;
	+ vector<double> sumlprup;
	+ vector<double> sumsqlprup;
	+ vector<long> nscanned;
	for ( int i = 0; ( maxScan() < 0 \|\| i < maxScan() ) && readEvent(); ++i ) {
	if ( !checkPartonBin() ) Throw<LesHouchesInitError>()
	<< "Found event in LesHouchesReader '" << name()
	<< "' which cannot be handeled by the assigned PartonExtractor '"
	<< partonExtractor()->name() << "'." << Exception::runerror;
	vector<int>::iterator idit =
	find(lprup.begin(), lprup.end(), hepeup.IDPRUP);
	int id = lprup.size();
	if ( idit == lprup.end() ) {
	lprup.push_back(hepeup.IDPRUP);
	newmax.push_back(0.0);
	neweve.push_back(0);
	oldeve.push_back(0);
	+ sumlprup.push_back(0.);
	+ sumsqlprup.push_back(0.);
	+ nscanned.push_back(0);
	} else {
	id = idit - lprup.begin();
	}
	++neve;
	++oldeve[id];
	oldsum += hepeup.XWGTUP;
	+ sumlprup[id] += hepeup.XWGTUP;
	+ sumsqlprup[id] += sqr(hepeup.XWGTUP);
	+ ++nscanned[id];
	if ( cacheFile() ) {
	if ( eventWeight() == 0.0 ) {
	++cuteve;
	continue;
	}
	cacheEvent();
	}
	++neweve[id];
	newmax[id] = max(newmax[id], abs(eventWeight()));
	if ( eventWeight() < 0.0 ) negw = true;
	- }
	- // std::cout << "eventWeight() = " << eventWeight() << endl;
	+ } //end of scanning events
	+ /*cout << "cross section estimates for different subproc:" << endl;
	+ for(int jj = 0; jj < lprup.size(); jj++) {
	+ cout << jj << " " << sumlprup[jj]/nscanned[jj] << endl;
	+ }*/
	xSecWeights.resize(oldeve.size(), 1.0);
	for ( int i = 0, N = oldeve.size(); i < N; ++i )
	if ( oldeve[i] ) xSecWeights[i] = double(neweve[i])/double(oldeve[i]);

	if ( maxScan() < 0 \|\| neve > NEvents() ) NEvents(neve - cuteve);

	if ( lprup.size() == heprup.LPRUP.size() ) {
	for ( int id = 0, N = lprup.size(); id < N; ++id ) {
	vector<int>::iterator idit =
	find(heprup.LPRUP.begin(), heprup.LPRUP.end(), hepeup.IDPRUP);
	if ( idit == heprup.LPRUP.end() ) {
	Throw<LesHouchesInitError>()
	<< "When scanning events, the LesHouschesReader '" << name()
	<< "' found undeclared processes." << Exception::warning;
	heprup.NPRUP = 0;
	break;
	}
	int idh = idit - heprup.LPRUP.begin();
	heprup.XMAXUP[idh] = newmax[id];
	}
	}
	if ( heprup.NPRUP == 0 ) {
	// No heprup block was supplied or something went wrong.
	heprup.NPRUP = lprup.size();
	heprup.LPRUP.resize(lprup.size());
	heprup.XMAXUP.resize(lprup.size());
	for ( int id = 0, N = lprup.size(); id < N; ++id ) {
	heprup.LPRUP[id] = lprup[id];
	heprup.XMAXUP[id] = newmax[id];
	}
	- } else if ( abs(heprup.IDWTUP) != 1 ) {
	+ }
	+ if ( abs(heprup.IDWTUP) != 1 ) {
	// Try to fix things if abs(heprup.IDWTUP) != 1.
	double sumxsec = 0.0;
	- for ( int id = 0; id < heprup.NPRUP; ++id ) sumxsec += heprup.XSECUP[id];
	+ for ( int id = 0; id < heprup.NPRUP; ++id ) {
	+ //set the cross section directly from the event weights read
	+ heprup.XSECUP[id] = sumlprup[id]/nscanned[id];
	+ heprup.XERRUP[id] = sqrt( (sumsqlprup[id]/nscanned[id] - sqr(sumlprup[id]/nscanned[id])) / nscanned[id] );
	+ heprup.XMAXUP[id] = newmax[id];
	+ sumxsec += heprup.XSECUP[id];
	+ }
	weightScale = picobarnnevesumxsec/oldsum;
	}
	}

	if ( cacheFile() ) closeCacheFile();

	if ( negw ) heprup.IDWTUP = min(-abs(heprup.IDWTUP), -1);
	-
	+
	return neve;

	}

	void LesHouchesReader::setWeightScale(long) {}

	void LesHouchesReader::initStat() {

	stats.reset();
	statmap.clear();
	if ( heprup.NPRUP <= 0 ) return;

	double sumx = 0.0;
	xSecWeights.resize(heprup.NPRUP, 1.0);
	maxWeights.clear();
	for ( int ip = 0; ip < heprup.NPRUP; ++ip ) {
	sumx = max(heprup.XMAXUP[ip]*xSecWeights[ip], sumx);
	statmap[heprup.LPRUP[ip]] =
	XSecStat(heprup.XMAXUP[ip]weightScalexSecWeights[ip]);
	maxWeights[heprup.LPRUP[ip]] = heprup.XMAXUP[ip];
	}
	stats.maxXSec(sumx*weightScale);
	maxFactor = 1.0;
	}

	void LesHouchesReader::increaseMaxXSec(CrossSection maxxsec) {
	for ( int i = 0; i < heprup.NPRUP; ++i )
	statmap[heprup.LPRUP[i]].maxXSec(statmap[heprup.LPRUP[i]].maxXSec()*
	maxxsec/stats.maxXSec());
	maxFactor *= maxxsec/stats.maxXSec();
	stats.maxXSec(maxxsec);
	}

	tXCombPtr LesHouchesReader::getXComb() {
	if ( lastXCombPtr() ) return lastXCombPtr();
	fillEvent();
	connectMothers();
	tcPBPair sel = createPartonBinInstances();
	tXCombPtr lastXC = xCombs()[sel];
	// clean up the old XComb object before switching to a new one
	if ( theLastXComb && theLastXComb != lastXC )
	theLastXComb->clean();
	theLastXComb = lastXC;
	lastXCombPtr()->subProcess(SubProPtr());
	lastXCombPtr()->setPartonBinInstances(partonBinInstances(),
	sqr(hepeup.SCALUP)*GeV2);
	lastXCombPtr()->lastAlphaS(hepeup.AQCDUP);
	lastXCombPtr()->lastAlphaEM(hepeup.AQEDUP);
	return lastXCombPtr();
	}

	tSubProPtr LesHouchesReader::getSubProcess() {
	getXComb();
	if ( subProcess() ) return subProcess();
	lastXCombPtr()->subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), this)));
	subProcess()->setOutgoing(outgoing().begin(), outgoing().end());
	subProcess()->setIntermediates(intermediates().begin(),
	intermediates().end());
	return subProcess();
	}

	void LesHouchesReader::fillEvent() {
	if ( !particleIndex.empty() ) return;
	particleIndex.clear();
	colourIndex.clear();
	colourIndex(0, tColinePtr());
	createParticles();
	createBeams();

	}

	void LesHouchesReader::reopen() {
	// If we didn't know how many events there were, we know now.
	if ( NEvents() <= 0 ) NEvents(position);
	++reopened;
	// How large fraction of the events have we actually used? And how
	// large will we have used if we go through the file again?
	double frac = double(stats.attempts())/double(NEvents());
	if ( frac*double(reopened + 1)/double(reopened) > 1.0 &&
	NEvents() - stats.attempts() <
	generator()->N() - generator()->currentEventNumber() ) {
	if(theReOpenAllowed)
	generator()->logWarning(LesHouchesReopenWarning()
	<< "Reopening LesHouchesReader '" << name()
	<< "' after accessing " << stats.attempts()
	<< " events out of "
	<< NEvents() << Exception::warning);
	else
	throw LesHouchesReopenWarning()
	<< "More events requested than available in LesHouchesReader "
	<< name() << Exception::runerror;
	}
	if ( cacheFile() ) {
	closeCacheFile();
	openReadCacheFile();
	if ( !uncacheEvent() ) Throw<LesHouchesReopenError>()
	<< "Could not reopen LesHouchesReader '" << name()
	<< "'." << Exception::runerror;
	} else {
	close();
	open();
	if ( !readEvent() ) Throw<LesHouchesReopenError>()
	<< "Could not reopen LesHouchesReader '" << name()
	<< "'." << Exception::runerror;
	}
	}

	void LesHouchesReader::reset() {
	particleIndex.clear();
	colourIndex.clear();
	if ( theLastXComb ) theLastXComb->clean();
	theLastXComb = tXCombPtr();
	}

	bool LesHouchesReader::readEvent() {

	reset();

	if ( !doReadEvent() ) return false;

	// If we are just skipping event we do not need to reweight or do
	// anything fancy.
	- if ( skipping ) { return true; }
	+ if ( skipping ) { return true; }

	if ( cacheFile() && !scanning ) { return true; }

	// Reweight according to the re- and pre-weights objects in the
	// LesHouchesReader base class.
	lastweight = reweight();
	-
	+
	if ( !reweightPDF && !cutEarly() ) return true;
	// We should try to reweight the PDFs or make early cuts here.

	fillEvent();

	double x1 = incoming().first->momentum().plus()/
	beams().first->momentum().plus();

	if ( reweightPDF &&
	inPDF.first && outPDF.first && inPDF.first != outPDF.first ) {
	if ( hepeup.XPDWUP.first <= 0.0 )
	hepeup.XPDWUP.first =
	inPDF.first->xfx(inData.first, incoming().first->dataPtr(),
	sqr(hepeup.SCALUP*GeV), x1);
	double xf = outPDF.first->xfx(inData.first, incoming().first->dataPtr(),
	sqr(hepeup.SCALUP*GeV), x1);
	lastweight *= xf/hepeup.XPDWUP.first;
	hepeup.XPDWUP.first = xf;
	}

	double x2 = incoming().second->momentum().minus()/
	beams().second->momentum().minus();

	if ( reweightPDF &&
	inPDF.second && outPDF.second && inPDF.second != outPDF.second ) {
	if ( hepeup.XPDWUP.second <= 0.0 )
	hepeup.XPDWUP.second =
	inPDF.second->xfx(inData.second, incoming().second->dataPtr(),
	sqr(hepeup.SCALUP*GeV), x2);
	double xf =
	outPDF.second->xfx(inData.second, incoming().second->dataPtr(),
	sqr(hepeup.SCALUP*GeV), x2);
	lastweight *= xf/hepeup.XPDWUP.second;
	hepeup.XPDWUP.second = xf;
	}

	if ( cutEarly() ) {
	if ( !cuts().initSubProcess((incoming().first->momentum() +
	incoming().second->momentum()).m2(),
	0.5*log(x1/x2)) ) lastweight = 0.0;
	tSubProPtr sub = getSubProcess();
	TmpTransform<tSubProPtr> tmp(sub, Utilities::getBoostToCM(sub->incoming()));
	if ( !cuts().passCuts(*sub) ) lastweight = 0.0;
	}

	return true;
	}

	double LesHouchesReader::getEvent() {
	if ( cacheFile() ) {
	if ( !uncacheEvent() ) reopen();
	} else {
	if ( !readEvent() ) reopen();
	}
	++position;

	double max = maxWeights[hepeup.IDPRUP]*maxFactor;
	-
	- return max != 0.0? eventWeight()/max: 0.0;
	-
	+ return max != 0.0? eventWeight()/max: 0.0;
	}

	void LesHouchesReader::skip(long n) {
	HoldFlag<> skipflag(skipping);
	while ( n-- ) getEvent();
	}

	double LesHouchesReader::reweight() {
	preweight = 1.0;
	if ( reweights.empty() && preweights.empty() &&
	!( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) )
	return 1.0;
	fillEvent();
	getSubProcess();
	for ( int i = 0, N = preweights.size(); i < N; ++i ) {
	preweights[i]->setXComb(lastXCombPtr());
	preweight *= preweights[i]->weight();
	}
	double weight = preweight;
	for ( int i = 0, N = reweights.size(); i < N; ++i ) {
	reweights[i]->setXComb(lastXCombPtr());
	weight *= reweights[i]->weight();
	}

	// If we are caching events we do not want to do CKKW reweighting.
	if ( cacheFile() ) return weight;

	if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) {
	CKKWHandler()->setXComb(lastXCombPtr());
	weight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW());
	}
	return weight;
	}

	bool LesHouchesReader::checkPartonBin() {

	// First find the positions of the incoming partons.
	pair< vector<int>, vector<int> > inc;
	for ( int i = 0; i < hepeup.NUP; ++i ) {
	if ( hepeup.ISTUP[i] == -9 ) {
	if ( inc.first.empty() ) inc.first.push_back(i);
	else if ( inc.second.empty() ) inc.second.push_back(i);
	}
	else if ( hepeup.ISTUP[i] == -1 ) {
	if ( inc.first.size() &&
	hepeup.MOTHUP[i].first == inc.first.back() + 1 )
	inc.first.push_back(i);
	else if ( inc.second.size() &&
	hepeup.MOTHUP[i].first == inc.second.back() + 1 )
	inc.second.push_back(i);
	else if ( inc.first.empty() ) {
	inc.first.push_back(-1);
	inc.first.push_back(i);
	}
	else if ( inc.second.empty() ) {
	inc.second.push_back(-1);
	inc.second.push_back(i);
	}
	else if ( inc.first.size() <= inc.second.size() )
	inc.first.push_back(i);
	else
	inc.second.push_back(i);
	}
	}

	// Now store the corresponding id numbers
	pair< vector<long>, vector<long> > ids;
	ids.first.push_back(inc.first[0] < 0? heprup.IDBMUP.first:
	hepeup.IDUP[inc.first[0]]);
	for ( int i = 1, N = inc.first.size(); i < N; ++i )
	ids.first.push_back(hepeup.IDUP[inc.first[i]]);
	ids.second.push_back(inc.second[0] < 0? heprup.IDBMUP.second:
	hepeup.IDUP[inc.second[0]]);
	for ( int i = 1, N = inc.second.size(); i < N; ++i )
	ids.second.push_back(hepeup.IDUP[inc.second[i]]);

	// Find the correct pair of parton bins.
	PBPair pbp;
	for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
	tcPBPtr curr = partonBins()[i].first;
	int icurr = inc.first.size() - 1;
	while ( curr && icurr >= 0 ) {
	if ( curr->parton()->id () != ids.first[icurr] ) break;
	curr = curr->incoming();
	--icurr;
	}
	if(!(!partonBins()[i].first->incoming() &&
	!partonBins()[i].first->particle() &&
	partonBins()[i].first->parton()->id () == ids.first[0] &&
	( inc.first.size()==1 \|\|
	(inc.first.size()==2 && ids.first[0]==ids.first[1]))) &&
	( curr \|\| icurr >= 0 ) ) continue;

	curr = partonBins()[i].second;
	icurr = inc.second.size() - 1;
	while ( curr && icurr >= 0 ) {
	if ( curr->parton()->id () != ids.second[icurr] ) break;
	curr = curr->incoming();
	--icurr;
	}
	if(!(!partonBins()[i].second->incoming() &&
	!partonBins()[i].second->particle() &&
	partonBins()[i].second->parton()->id () == ids.second[0] &&
	( inc.second.size()==1 \|\|
	(inc.second.size()==2 && ids.second[0]==ids.second[1]))) &&
	( curr \|\| icurr >= 0 ) ) continue;

	pbp = partonBins()[i];
	}

	// If we are only checking we return here.
	return ( pbp.first && pbp.second );

	}

	namespace {
	bool recursionNotNull(tcPBPtr bin, tcPPtr p) {
	while ( bin && p ) {
	if ( p->dataPtr() != bin->parton() ) break;
	bin = bin->incoming();
	p = p->parents().size()? p->parents()[0]: tPPtr();
	}
	return bin \|\| p;
	}
	}


	tcPBPair LesHouchesReader::createPartonBinInstances() {
	tcPBPair sel;
	for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
	tcPBPtr bin = partonBins()[i].first;
	tcPPtr p = incoming().first;
	if ( recursionNotNull(bin,p) ) continue;
	bin = partonBins()[i].second;
	p = incoming().second;
	if ( recursionNotNull(bin,p) ) continue;
	sel = partonBins()[i];
	break;
	}
	if ( !sel.first \|\| !sel.second ) Throw<LesHouchesInconsistencyError>()
	<< "Could not find appropriate PartonBin objects for event produced by "
	<< "LesHouchesReader '" << name() << "'." << Exception::runerror;

	Direction<0> dir(true);
	thePartonBinInstances.first =
	new_ptr(PartonBinInstance(incoming().first, sel.first,
	-sqr(hepeup.SCALUP*GeV)));
	if ( thePartonBinInstances.first->xi() > 1.00001 ) {
	Throw<LesHouchesInconsistencyError>()
	<< "Found an event with momentum fraction larger than unity (x1="
	<< thePartonBinInstances.first->xi()
	<< "). The event will be skipped." << Exception::warning;
	throw Veto();
	}
	dir.reverse();
	thePartonBinInstances.second =
	new_ptr(PartonBinInstance(incoming().second, sel.second,
	-sqr(hepeup.SCALUP*GeV)));

	if ( thePartonBinInstances.second->xi() > 1.00001 ) {
	Throw<LesHouchesInconsistencyError>()
	<< "Found an event with momentum fraction larger than unity (x2="
	<< thePartonBinInstances.second->xi()
	<< "). The event will be skipped." << Exception::warning;
	throw Veto();
	}
	return sel;

	}


	void LesHouchesReader::createParticles() {
	theBeams = PPair();
	theIncoming = PPair();
	theOutgoing = PVector();
	theIntermediates = PVector();
	for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
	if ( !hepeup.IDUP[i] ) continue;
	Lorentz5Momentum mom(hepeup.PUP[i][0]GeV, hepeup.PUP[i][1]GeV,
	hepeup.PUP[i][2]GeV, hepeup.PUP[i][3]GeV,
	hepeup.PUP[i][4]*GeV);
	// cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
	if(theMomentumTreatment == 1) mom.rescaleEnergy();
	else if(theMomentumTreatment == 2) mom.rescaleMass();
	// cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
	PDPtr pd = getParticleData(hepeup.IDUP[i]);
	if (!pd) {
	Throw<LesHouchesInitError>()
	<< "LesHouchesReader '" << name() << "' found unknown particle ID "
	<< hepeup.IDUP[i]
	<< " in Les Houches common block structure.\n"
	<< "You need to define the new particle in an input file.\n"
	<< Exception::runerror;
	}
	if ( ! pd->coloured()
	&& ( hepeup.ICOLUP[i].first != 0 \|\| hepeup.ICOLUP[i].second != 0 ) ) {
	Throw<LesHouchesInconsistencyError>()
	<< "LesHouchesReader " << name() << ": " << pd->PDGName()
	<< " is not a coloured particle.\nIt should not have "
	<< "(anti-)colour lines " << hepeup.ICOLUP[i].first
	<< ' ' << hepeup.ICOLUP[i].second
	<< " set; the event file needs to be fixed."
	<< Exception::runerror;
	}
	PPtr p = pd->produceParticle(mom);
	if(hepeup.ICOLUP[i].first>=0 && hepeup.ICOLUP[i].second >=0) {
	tColinePtr c = colourIndex(hepeup.ICOLUP[i].first);
	if ( c ) {
	c->addColoured(p);
	}
	c = colourIndex(hepeup.ICOLUP[i].second);
	if ( c ) c->addAntiColoured(p);
	}
	else {
	tColinePtr c1 = colourIndex(abs(hepeup.ICOLUP[i].first ));
	tColinePtr c2 = colourIndex(abs(hepeup.ICOLUP[i].second));
	if(pd->hasColour()) {
	c1->addColouredIndexed(p,1);
	c2->addColouredIndexed(p,2);
	}
	else {
	c1->addAntiColouredIndexed(p,1);
	c2->addAntiColouredIndexed(p,2);
	}
	}
	particleIndex(i + 1, p);
	switch ( hepeup.ISTUP[i] ) {
	case -9:
	if ( !theBeams.first ) theBeams.first = p;
	else if ( !theBeams.second ) theBeams.second = p;
	else Throw<LesHouchesInconsistencyError>()
	<< "To many incoming beam particles in the LesHouchesReader '"
	<< name() << "'." << Exception::runerror;
	break;
	case -1:
	if ( !theIncoming.first ) theIncoming.first = p;
	else if ( !theIncoming.second ) theIncoming.second = p;
	else if ( particleIndex(theIncoming.first) == hepeup.MOTHUP[i].first )
	theIncoming.first = p;
	else if ( particleIndex(theIncoming.second) == hepeup.MOTHUP[i].first )
	theIncoming.second = p;
	else Throw<LesHouchesInconsistencyError>()
	<< "To many incoming particles to hard subprocess in the "
	<< "LesHouchesReader '" << name() << "'." << Exception::runerror;
	p->scale(sqr(hepeup.SCALUP*GeV));
	break;
	case 1:
	theOutgoing.push_back(p);
	p->scale(sqr(hepeup.SCALUP*GeV));
	break;
	case -2:
	case 2:
	case 3:
	theIntermediates.push_back(p);
	break;
	default:
	Throw<LesHouchesInconsistencyError>()
	<< "Unknown status code (" << hepeup.ISTUP[i]
	<< ") in the LesHouchesReader '" << name() << "'."
	<< Exception::runerror;
	}

	// value 9 is defined as "Unknown or unpolarized particles"
	double spinup = hepeup.SPINUP[i];
	if ( abs(spinup - 9) < 1.0e-3 )
	spinup = 0.;
	if ( spinup < -1. \|\| spinup > 1. ) {
	Throw<LesHouchesInconsistencyError>()
	<< "Polarization must be between -1 and 1, not "
	<< spinup << " as found in the "
	<< "LesHouches event file.\nThe event file needs to be fixed."
	<< Exception::runerror;
	}
	if( theIncludeSpin
	&& abs(pd->id()) == ParticleID::tauminus
	&& spinup !=0) {
	if(pd->iSpin() == PDT::Spin1Half ) {
	vector<Helicity::SpinorWaveFunction> wave;
	Helicity::SpinorWaveFunction(wave,p,Helicity::outgoing,true);
	RhoDMatrix rho(pd->iSpin(),true);
	rho(0,0) = 0.5*(1.-spinup);
	rho(1,1) = 0.5*(1.+spinup);
	p->spinInfo()->rhoMatrix() = rho;
	p->spinInfo()-> DMatrix() = rho;
	}
	}
	}

	// check the colour flows, and if necessary create any sources/sinks
	// hard process
	// get the particles in the hard process
	PVector external;
	for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
	unsigned int moth;
	switch ( hepeup.ISTUP[i] ) {
	case -1:
	external.push_back(particleIndex.find(i+1));
	break;
	case 1: case 2: case 3:
	moth = hepeup.MOTHUP[i].first;
	if(moth!=0 && (hepeup.ISTUP[moth]==-1\|\|hepeup.ISTUP[moth]==-2\|\|
	hepeup.ISTUP[moth]==-9))
	external.push_back(particleIndex.find(i+1));
	moth = hepeup.MOTHUP[i].second;
	if(moth!=0 && (hepeup.ISTUP[moth]==-1\|\|hepeup.ISTUP[moth]==-2\|\|
	hepeup.ISTUP[moth]==-9))
	external.push_back(particleIndex.find(i+1));
	break;
	case -2: case -9: default:
	break;
	}
	}
	// check the incoming/outgoing lines match
	vector<tColinePtr> unMatchedColour,unMatchedAntiColour;
	for(unsigned int ix=0;ix<external.size();++ix) {
	vector<tcColinePtr>
	col = external[ix]->colourInfo()-> colourLines();
	vector<tcColinePtr>
	anti = external[ix]->colourInfo()->antiColourLines();
	if(hepeup.ISTUP[particleIndex(external[ix])-1]<0)
	swap(col,anti);
	if(!col.empty()) {
	for(unsigned int ic1=0;ic1<col.size();++ic1) {
	bool matched=false;
	for(unsigned int iy=0;iy<external.size();++iy) {
	vector<tcColinePtr> col2;
	if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
	if(external[iy]->colourInfo()->colourLines().empty()) continue;
	col2 = external[iy]->colourInfo()->colourLines();
	}
	else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
	if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
	col2 = external[iy]->colourInfo()->antiColourLines();
	}
	for(unsigned int ic2=0;ic2<col2.size();++ic2) {
	if(col[ic1]==col2[ic2]) {
	matched=true;
	break;
	}
	}
	if(matched) break;
	}
	if(!matched) unMatchedColour.push_back(const_ptr_cast<tColinePtr>(col[ic1]));
	}
	}
	if(!anti.empty()) {
	for(unsigned int ic1=0;ic1<col.size();++ic1) {
	bool matched=false;
	for(unsigned int iy=0;iy<external.size();++iy) {
	vector<tcColinePtr> anti2;
	if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
	if(external[iy]->colourInfo()->colourLines().empty()) continue;
	anti2 = external[iy]->colourInfo()->antiColourLines();
	}
	else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
	if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
	anti2 = external[iy]->colourInfo()->colourLines();
	}
	for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
	if(col[ic1]==anti2[ic2]) {
	matched=true;
	break;
	}
	}
	if(matched) break;
	}
	if(!matched) unMatchedAntiColour.push_back(const_ptr_cast<tColinePtr>(anti[ic1]));
	}
	}
	}

	// might have source/sink
	if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) {
	if(unMatchedColour.size() == 3 ) {
	unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1],
	unMatchedColour[2]);
	}
	else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
	Throw<LesHouchesInconsistencyError>()
	<< "LesHouchesReader '" << name() << "' found inconsistent colour "
	<< "flow in Les Houches common block structure for hard process.\n"
	<< hepeup << Exception::runerror;
	}
	if(unMatchedAntiColour.size() == 3 ) {
	unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1],
	unMatchedAntiColour[2]);
	}
	else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
	Throw<LesHouchesInconsistencyError>()
	<< "LesHouchesReader '" << name() << "' found inconsistent colour "
	<< "flow in Les Houches common block structure for hard process.\n"
	<< hepeup << Exception::runerror;
	}
	}

	// any subsequent decays
	for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
	if(hepeup.ISTUP[i] !=2 && hepeup.ISTUP[i] !=3) continue;
	PVector external;
	external.push_back(particleIndex.find(i+1));
	for ( int j = 0; j < N; ++j ) {
	if(hepeup.MOTHUP[j].first==i+1\|\| hepeup.MOTHUP[j].second==i+1)
	external.push_back(particleIndex.find(j+1));
	}
	// check the incoming/outgoing lines match
	vector<tColinePtr> unMatchedColour,unMatchedAntiColour;
	for(unsigned int ix=0;ix<external.size();++ix) {
	vector<tcColinePtr>
	col = external[ix]->colourInfo()-> colourLines();
	vector<tcColinePtr>
	anti = external[ix]->colourInfo()->antiColourLines();
	if(ix==0) swap(col,anti);
	if(!col.empty()) {
	for(unsigned int ic1=0;ic1<col.size();++ic1) {
	bool matched=false;
	for(unsigned int iy=0;iy<external.size();++iy) {
	if(iy==ix) continue;
	vector<tcColinePtr> col2;
	if(iy==0) {
	if(external[iy]->colourInfo()->colourLines().empty()) continue;
	col2 = external[iy]->colourInfo()->colourLines();
	}
	else {
	if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
	col2 = external[iy]->colourInfo()->antiColourLines();
	}
	for(unsigned int ic2=0;ic2<col2.size();++ic2) {
	if(col[ic1]==col2[ic2]) {
	matched=true;
	break;
	}
	}
	if(matched) break;
	}
	if(!matched) unMatchedColour.push_back(const_ptr_cast<tColinePtr>(col[ic1]));
	}
	}
	if(!anti.empty()) {
	for(unsigned int ic1=0;ic1<anti.size();++ic1) {
	bool matched=false;
	for(unsigned int iy=0;iy<external.size();++iy) {
	if(iy==ix) continue;
	vector<tcColinePtr> anti2;
	if(iy==0) {
	if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
	anti2 = external[iy]->colourInfo()->antiColourLines();
	}
	else {
	if(external[iy]->colourInfo()->colourLines().empty()) continue;
	anti2 = external[iy]->colourInfo()->colourLines();
	}
	for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
	if(anti[ic1]==anti2[ic2]) {
	matched=true;
	break;
	}
	}
	if(matched) break;
	}
	if(!matched) unMatchedAntiColour.push_back(const_ptr_cast<tColinePtr>(anti[ic1]));
	}
	}
	}
	// might have source/sink
	if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) {
	if(unMatchedColour.size() == 3 ) {
	unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1],
	unMatchedColour[2]);
	}
	else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
	Throw<LesHouchesInconsistencyError>()
	<< "LesHouchesReader '" << name() << "' found inconsistent colour "
	<< "flow in Les Houches common block structure for decay of \n"
	<< *external[0] << "\n"
	<< hepeup << Exception::runerror;
	}
	if(unMatchedAntiColour.size() == 3 ) {
	unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1],
	unMatchedAntiColour[2]);
	}
	else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
	Throw<LesHouchesInconsistencyError>()
	<< "LesHouchesReader '" << name() << "' found inconsistent colour "
	<< "flow in Les Houches common block structure for decay of\n"
	<< *external[0] << "\n"
	<< hepeup << Exception::runerror;
	}
	}
	}

	}

	void LesHouchesReader::createBeams() {

	if ( !theBeams.first && dynamic_ptr_cast<Ptr<NoPDF>::tcp>(inPDF.first) ) {
	theBeams.first = theIncoming.first;
	}
	else if ( !theBeams.first ) {
	theBeams.first = getParticleData(heprup.IDBMUP.first)->produceParticle();
	double m = theBeams.first->mass()/GeV;
	theBeams.first->set5Momentum
	(Lorentz5Momentum(ZERO, ZERO,
	sqrt(sqr(heprup.EBMUP.first) - sqr(m))*GeV,
	heprup.EBMUP.firstGeV, mGeV));
	hepeup.IDUP.push_back(heprup.IDBMUP.first);
	hepeup.ISTUP.push_back(-9);
	hepeup.MOTHUP.push_back(make_pair(0, 0));
	hepeup.ICOLUP.push_back(make_pair(0, 0));
	hepeup.VTIMUP.push_back(0.0);
	hepeup.SPINUP.push_back(0.0);
	particleIndex(hepeup.IDUP.size(), theBeams.first);
	hepeup.MOTHUP[particleIndex(theIncoming.first) - 1].first =
	hepeup.IDUP.size();
	}
	if ( !theBeams.second && dynamic_ptr_cast<Ptr<NoPDF>::tcp>(inPDF.second) ) {
	theBeams.second = theIncoming.second;
	}
	else if ( !theBeams.second ) {
	theBeams.second = getParticleData(heprup.IDBMUP.second)->produceParticle();
	double m = theBeams.second->mass()/GeV;
	theBeams.second->set5Momentum
	(Lorentz5Momentum(ZERO, ZERO,
	-sqrt(sqr(heprup.EBMUP.second) - sqr(m))*GeV,
	heprup.EBMUP.secondGeV, mGeV));
	hepeup.IDUP.push_back(heprup.IDBMUP.second);
	hepeup.ISTUP.push_back(-9);
	hepeup.MOTHUP.push_back(make_pair(0, 0));
	hepeup.ICOLUP.push_back(make_pair(0, 0));
	hepeup.VTIMUP.push_back(0.0);
	hepeup.SPINUP.push_back(0.0);
	particleIndex(hepeup.IDUP.size(), theBeams.second);
	hepeup.MOTHUP[particleIndex(theIncoming.second) - 1].first =
	hepeup.IDUP.size();
	}
	}

	void LesHouchesReader::connectMothers() {
	const ObjectIndexer<long,Particle> & pi = particleIndex;
	for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
	if ( pi(hepeup.MOTHUP[i].first) )
	pi(hepeup.MOTHUP[i].first)->addChild(pi(i + 1));
	if ( pi(hepeup.MOTHUP[i].second)
	&& hepeup.MOTHUP[i].second != hepeup.MOTHUP[i].first )
	pi(hepeup.MOTHUP[i].second)->addChild(pi(i + 1));
	}
	}

	void LesHouchesReader::openReadCacheFile() {
	if ( cacheFile() ) closeCacheFile();
	cacheFile().open(cacheFileName(), "r");
	position = 0;
	}

	void LesHouchesReader::openWriteCacheFile() {
	if ( cacheFile() ) closeCacheFile();
	cacheFile().open(cacheFileName(), "w");
	}

	void LesHouchesReader::closeCacheFile() {
	cacheFile().close();
	}

	void LesHouchesReader::cacheEvent() const {
	static vector<char> buff;
	cacheFile().write(&hepeup.NUP, sizeof(hepeup.NUP));
	buff.resize(eventSize(hepeup.NUP));
	char * pos = &buff[0];
	pos = mwrite(pos, hepeup.IDPRUP);
	pos = mwrite(pos, hepeup.XWGTUP);
	pos = mwrite(pos, hepeup.XPDWUP);
	pos = mwrite(pos, hepeup.SCALUP);
	pos = mwrite(pos, hepeup.AQEDUP);
	pos = mwrite(pos, hepeup.AQCDUP);
	pos = mwrite(pos, hepeup.IDUP[0], hepeup.NUP);
	pos = mwrite(pos, hepeup.ISTUP[0], hepeup.NUP);
	pos = mwrite(pos, hepeup.MOTHUP[0], hepeup.NUP);
	pos = mwrite(pos, hepeup.ICOLUP[0], hepeup.NUP);
	for ( int i = 0; i < hepeup.NUP; ++i )
	pos = mwrite(pos, hepeup.PUP[i][0], 5);
	pos = mwrite(pos, hepeup.VTIMUP[0], hepeup.NUP);
	pos = mwrite(pos, hepeup.SPINUP[0], hepeup.NUP);
	pos = mwrite(pos, lastweight);
	pos = mwrite(pos, optionalWeights);
	for(int ff = 0; ff < optionalWeightsNames.size(); ff++) {
	pos = mwrite(pos, optionalWeightsNames[ff]);
	}
	pos = mwrite(pos, optionalnpLO);
	pos = mwrite(pos, optionalnpNLO);
	pos = mwrite(pos, preweight);
	cacheFile().write(&buff[0], buff.size(), 1);
	}

	bool LesHouchesReader::uncacheEvent() {
	reset();
	static vector<char> buff;
	if ( cacheFile().read(&hepeup.NUP, sizeof(hepeup.NUP)) != 1 )
	return false;
	buff.resize(eventSize(hepeup.NUP));
	if ( cacheFile().read(&buff[0], buff.size()) != 1 ) return false;
	const char * pos = &buff[0];
	pos = mread(pos, hepeup.IDPRUP);
	pos = mread(pos, hepeup.XWGTUP);
	pos = mread(pos, hepeup.XPDWUP);
	pos = mread(pos, hepeup.SCALUP);
	pos = mread(pos, hepeup.AQEDUP);
	pos = mread(pos, hepeup.AQCDUP);
	hepeup.IDUP.resize(hepeup.NUP);
	pos = mread(pos, hepeup.IDUP[0], hepeup.NUP);
	hepeup.ISTUP.resize(hepeup.NUP);
	pos = mread(pos, hepeup.ISTUP[0], hepeup.NUP);
	hepeup.MOTHUP.resize(hepeup.NUP);
	pos = mread(pos, hepeup.MOTHUP[0], hepeup.NUP);
	hepeup.ICOLUP.resize(hepeup.NUP);
	pos = mread(pos, hepeup.ICOLUP[0], hepeup.NUP);
	hepeup.PUP.resize(hepeup.NUP, vector<double>(5));
	for ( int i = 0; i < hepeup.NUP; ++i )
	pos = mread(pos, hepeup.PUP[i][0], 5);
	hepeup.VTIMUP.resize(hepeup.NUP);
	pos = mread(pos, hepeup.VTIMUP[0], hepeup.NUP);
	hepeup.SPINUP.resize(hepeup.NUP);
	pos = mread(pos, hepeup.SPINUP[0], hepeup.NUP);
	pos = mread(pos, lastweight);
	pos = mread(pos, optionalWeights);
	for(int ff = 0; ff < optionalWeightsNames.size(); ff++) {
	pos = mread(pos, optionalWeightsNames[ff]);
	}
	pos = mread(pos, optionalnpLO);
	pos = mread(pos, optionalnpNLO);
	pos = mread(pos, preweight);

	// If we are skipping, we do not have to do anything else.
	if ( skipping ) return true;

	if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) {
	// The cached event has not been submitted to CKKW reweighting, so
	// we do that now.
	fillEvent();
	getSubProcess();
	CKKWHandler()->setXComb(lastXCombPtr());
	lastweight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW());
	}
	return true;
	}

	void LesHouchesReader::persistentOutput(PersistentOStream & os) const {
	os << heprup.IDBMUP << heprup.EBMUP << heprup.PDFGUP << heprup.PDFSUP
	<< heprup.IDWTUP << heprup.NPRUP << heprup.XSECUP << heprup.XERRUP
	<< heprup.XMAXUP << heprup.LPRUP << hepeup.NUP << hepeup.IDPRUP
	<< hepeup.XWGTUP << hepeup.XPDWUP << hepeup.SCALUP << hepeup.AQEDUP
	<< hepeup.AQCDUP << hepeup.IDUP << hepeup.ISTUP << hepeup.MOTHUP
	<< hepeup.ICOLUP << hepeup.PUP << hepeup.VTIMUP << hepeup.SPINUP
	<< inData << inPDF << outPDF << thePartonExtractor << theCKKW
	<< thePartonBins << theXCombs << theCuts << theNEvents << position
	<< reopened << theMaxScan << isActive
	<< theCacheFileName << doCutEarly << stats << statmap
	<< thePartonBinInstances
	<< theBeams << theIncoming << theOutgoing << theIntermediates
	<< reweights << preweights << preweight << reweightPDF << doInitPDFs
	<< theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO
	<< maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights
	<< theMomentumTreatment << useWeightWarnings << theReOpenAllowed
	<< theIncludeSpin;
	}

	void LesHouchesReader::persistentInput(PersistentIStream & is, int) {
	if ( cacheFile() ) closeCacheFile();
	is >> heprup.IDBMUP >> heprup.EBMUP >> heprup.PDFGUP >> heprup.PDFSUP
	>> heprup.IDWTUP >> heprup.NPRUP >> heprup.XSECUP >> heprup.XERRUP
	>> heprup.XMAXUP >> heprup.LPRUP >> hepeup.NUP >> hepeup.IDPRUP
	>> hepeup.XWGTUP >> hepeup.XPDWUP >> hepeup.SCALUP >> hepeup.AQEDUP
	>> hepeup.AQCDUP >> hepeup.IDUP >> hepeup.ISTUP >> hepeup.MOTHUP
	>> hepeup.ICOLUP >> hepeup.PUP >> hepeup.VTIMUP >> hepeup.SPINUP
	>> inData >> inPDF >> outPDF >> thePartonExtractor >> theCKKW
	>> thePartonBins >> theXCombs >> theCuts >> theNEvents >> position
	>> reopened >> theMaxScan >> isActive
	>> theCacheFileName >> doCutEarly >> stats >> statmap
	>> thePartonBinInstances
	>> theBeams >> theIncoming >> theOutgoing >> theIntermediates
	>> reweights >> preweights >> preweight >> reweightPDF >> doInitPDFs
	>> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO
	>> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights
	>> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed
	>> theIncludeSpin;
	}

	AbstractClassDescription<LesHouchesReader>
	LesHouchesReader::initLesHouchesReader;
	// Definition of the static class description member.

	void LesHouchesReader::setBeamA(long id) { heprup.IDBMUP.first = id; }
	long LesHouchesReader::getBeamA() const { return heprup.IDBMUP.first; }
	void LesHouchesReader::setBeamB(long id) { heprup.IDBMUP.second = id; }
	long LesHouchesReader::getBeamB() const { return heprup.IDBMUP.second; }
	void LesHouchesReader::setEBeamA(Energy e) { heprup.EBMUP.first = e/GeV; }
	Energy LesHouchesReader::getEBeamA() const { return heprup.EBMUP.first*GeV; }
	void LesHouchesReader::setEBeamB(Energy e) { heprup.EBMUP.second = e/GeV; }
	Energy LesHouchesReader::getEBeamB() const { return heprup.EBMUP.second*GeV; }
	void LesHouchesReader::setPDFA(PDFPtr pdf) { inPDF.first = pdf; }
	PDFPtr LesHouchesReader::getPDFA() const { return inPDF.first; }
	void LesHouchesReader::setPDFB(PDFPtr pdf) { inPDF.second = pdf; }
	PDFPtr LesHouchesReader::getPDFB() const { return inPDF.second; }

	void LesHouchesReader::Init() {

	static ClassDocumentation<LesHouchesReader> documentation
	("ThePEG::LesHouchesReader is an abstract base class to be used "
	"for objects which reads event files or streams from matrix element "
	"generators.");

	static Parameter<LesHouchesReader,long> interfaceBeamA
	("BeamA",
	"The PDG id of the incoming particle along the positive z-axis. "
	"If zero the corresponding information is to be deduced from the "
	"event stream/file.",
	0, 0, 0, 0,
	true, false, false,
	&LesHouchesReader::setBeamA,
	&LesHouchesReader::getBeamA, 0, 0, 0);

	static Parameter<LesHouchesReader,long> interfaceBeamB
	("BeamB",
	"The PDG id of the incoming particle along the negative z-axis. "
	"If zero the corresponding information is to be deduced from the "
	"event stream/file.",
	0, 0, 0, 0,
	true, false, false,
	&LesHouchesReader::setBeamB,
	&LesHouchesReader::getBeamB, 0, 0, 0);

	static Parameter<LesHouchesReader,Energy> interfaceEBeamA
	("EBeamA",
	"The energy of the incoming particle along the positive z-axis. "
	"If zero the corresponding information is to be deduced from the "
	"event stream/file.",
	0, GeV, ZERO, ZERO, 1000000000.0*GeV,
	true, false, true,
	&LesHouchesReader::setEBeamA, &LesHouchesReader::getEBeamA, 0, 0, 0);

	static Parameter<LesHouchesReader,Energy> interfaceEBeamB
	("EBeamB",
	"The energy of the incoming particle along the negative z-axis. "
	"If zero the corresponding information is to be deduced from the "
	"event stream/file.",
	0, GeV, ZERO, ZERO, 1000000000.0*GeV,
	true, false, true,
	&LesHouchesReader::setEBeamB, &LesHouchesReader::getEBeamB, 0, 0, 0);

	static Reference<LesHouchesReader,PDFBase> interfacePDFA
	("PDFA",
	"The PDF used for incoming particle along the positive z-axis. "
	"If null the corresponding information is to be deduced from the "
	"event stream/file.",
	0, true, false, true, true, false,
	&LesHouchesReader::setPDFA, &LesHouchesReader::getPDFA, 0);

	static Reference<LesHouchesReader,PDFBase> interfacePDFB
	("PDFB",
	"The PDF used for incoming particle along the negative z-axis. "
	"If null the corresponding information is to be deduced from the "
	"event stream/file.",
	0, true, false, true, true, false,
	&LesHouchesReader::setPDFB, &LesHouchesReader::getPDFB, 0);

	static Parameter<LesHouchesReader,long> interfaceMaxScan
	("MaxScan",
	"The maximum number of events to scan to obtain information about "
	"processes and cross section in the intialization.",
	&LesHouchesReader::theMaxScan, -1, 0, 0,
	true, false, false);

	static Parameter<LesHouchesReader,string> interfaceCacheFileName
	("CacheFileName",
	"Name of file used to cache the events from the reader in a fast-readable "
	"form. If empty, no cache file will be generated.",
	&LesHouchesReader::theCacheFileName, "",
	true, false);
	interfaceCacheFileName.fileType();

	static Switch<LesHouchesReader,bool> interfaceCutEarly
	("CutEarly",
	"Determines whether to apply cuts to events before converting to "
	"ThePEG format.",
	&LesHouchesReader::doCutEarly, true, true, false);
	static SwitchOption interfaceCutEarlyYes
	(interfaceCutEarly,
	"Yes",
	"Event are cut before converted.",
	true);
	static SwitchOption interfaceCutEarlyNo
	(interfaceCutEarly,
	"No",
	"Events are not cut before converted.",
	false);

	static Reference<LesHouchesReader,PartonExtractor> interfacePartonExtractor
	("PartonExtractor",
	"The PartonExtractor object used to construct remnants. If no object is "
	"provided the LesHouchesEventHandler object must provide one instead.",
	&LesHouchesReader::thePartonExtractor, true, false, true, true, false);


	static Reference<LesHouchesReader,Cuts> interfaceCuts
	("Cuts",
	"The Cuts object to be used for this reader. Note that these "
	"must not be looser cuts than those used in the actual generation. "
	"If no object is provided the LesHouchesEventHandler object must "
	"provide one instead.",
	&LesHouchesReader::theCuts, true, false, true, true, false);

	static RefVector<LesHouchesReader,ReweightBase> interfaceReweights
	("Reweights",
	"A list of ThePEG::ReweightBase objects to modify this the weight of "
	"this reader.",
	&LesHouchesReader::reweights, 0, false, false, true, false);

	static RefVector<LesHouchesReader,ReweightBase> interfacePreweights
	("Preweights",
	"A list of ThePEG::ReweightBase objects to bias the phase space for this "
	"reader without influencing the actual cross section.",
	&LesHouchesReader::preweights, 0, false, false, true, false);

	static Switch<LesHouchesReader,bool> interfaceReweightPDF
	("ReweightPDF",
	"If the PDFs used in the generation for this reader is different "
	"from the ones assumed by the associated PartonExtractor object, "
	"should the events be reweighted to fit the latter?",
	&LesHouchesReader::reweightPDF, false, true, false);
	static SwitchOption interfaceReweightPDFNo
	(interfaceReweightPDF, "No", "The event weights are kept as they are.",
	false);
	static SwitchOption interfaceReweightPDFYes
	(interfaceReweightPDF,
	"Yes", "The events are reweighted.", true);

	static Switch<LesHouchesReader,bool> interfaceInitPDFs
	("InitPDFs",
	"If no PDFs were specified in <interface>PDFA</interface> or "
	"<interface>PDFB</interface>for this reader, try to extract the "
	"information from the event file and assign the relevant PDFBase"
	"objects when the reader is initialized.",
	&LesHouchesReader::doInitPDFs, false, true, false);
	static SwitchOption interfaceInitPDFsYes
	(interfaceInitPDFs,
	"Yes",
	"Extract PDFs during initialization.",
	true);
	static SwitchOption interfaceInitPDFsNo
	(interfaceInitPDFs,
	"No",
	"Do not extract PDFs during initialization.",
	false);

	static Parameter<LesHouchesReader,int> interfaceMaxMultCKKW
	("MaxMultCKKW",
	"If this reader is to be used (possibly together with others) for CKKW-"
	"reweighting and veto, this should give the multiplicity of outgoing "
	"particles in the highest multiplicity matrix element in the group. "
	"If set to zero, no CKKW procedure should be applied.",
	&LesHouchesReader::theMaxMultCKKW, 0, 0, 0,
	true, false, Interface::lowerlim);

	static Parameter<LesHouchesReader,int> interfaceMinMultCKKW
	("MinMultCKKW",
	"If this reader is to be used (possibly together with others) for CKKW-"
	"reweighting and veto, this should give the multiplicity of outgoing "
	"particles in the lowest multiplicity matrix element in the group. If "
	"larger or equal to <interface>MaxMultCKKW</interface>, no CKKW "
	"procedure should be applied.",
	&LesHouchesReader::theMinMultCKKW, 0, 0, 0,
	true, false, Interface::lowerlim);


	static Switch<LesHouchesReader,unsigned int> interfaceMomentumTreatment
	("MomentumTreatment",
	"Treatment of the momenta supplied by the interface",
	&LesHouchesReader::theMomentumTreatment, 0, false, false);
	static SwitchOption interfaceMomentumTreatmentAccept
	(interfaceMomentumTreatment,
	"Accept",
	"Just accept the momenta given",
	0);
	static SwitchOption interfaceMomentumTreatmentRescaleEnergy
	(interfaceMomentumTreatment,
	"RescaleEnergy",
	"Rescale the energy supplied so it is consistent with the mass",
	1);
	static SwitchOption interfaceMomentumTreatmentRescaleMass
	(interfaceMomentumTreatment,
	"RescaleMass",
	"Rescale the mass supplied so it is consistent with the"
	" energy and momentum",
	2);


	static Switch<LesHouchesReader,bool> interfaceWeightWarnings
	("WeightWarnings",
	"Determines if warnings about possible weight incompatibilities should "
	"be issued when this reader is initialized.",
	&LesHouchesReader::useWeightWarnings, true, true, false);
	static SwitchOption interfaceWeightWarningsWarnAboutWeights
	(interfaceWeightWarnings,
	"WarnAboutWeights",
	"Warn about possible incompatibilities with the weight option in the "
	"Les Houches common block and the requested weight treatment.",
	true);
	static SwitchOption interfaceWeightWarningsDontWarnAboutWeights
	(interfaceWeightWarnings,
	"DontWarnAboutWeights",
	"Do not warn about possible incompatibilities with the weight option "
	"in the Les Houches common block and the requested weight treatment.",
	false);

	static Switch<LesHouchesReader,bool> interfaceAllowedTopReOpen
	("AllowedToReOpen",
	"Can the file be reopened if more events are requested than the file contains?",
	&LesHouchesReader::theReOpenAllowed, true, false, false);
	static SwitchOption interfaceAllowedTopReOpenYes
	(interfaceAllowedTopReOpen,
	"Yes",
	"Allowed to reopen the file",
	true);
	static SwitchOption interfaceAllowedTopReOpenNo
	(interfaceAllowedTopReOpen,
	"No",
	"Not allowed to reopen the file",
	false);

	static Switch<LesHouchesReader,bool> interfaceIncludeSpin
	("IncludeSpin",
	"Use the spin information present in the event file, for tau leptons"
	" only as this is the only case which makes any sense",
	&LesHouchesReader::theIncludeSpin, true, false, false);
	static SwitchOption interfaceIncludeSpinYes
	(interfaceIncludeSpin,
	"Yes",
	"Use the spin information",
	true);
	static SwitchOption interfaceIncludeSpinNo
	(interfaceIncludeSpin,
	"No",
	"Don't use the spin information",
	false);



	interfaceCuts.rank(8);
	interfacePartonExtractor.rank(7);
	interfaceBeamA.rank(5);
	interfaceBeamB.rank(4);
	interfaceEBeamA.rank(3);
	interfaceEBeamB.rank(2);
	interfaceMaxMultCKKW.rank(1.5);
	interfaceMinMultCKKW.rank(1.0);

	interfaceBeamA.setHasDefault(false);
	interfaceBeamB.setHasDefault(false);
	interfaceEBeamA.setHasDefault(false);
	interfaceEBeamB.setHasDefault(false);
	interfaceMaxMultCKKW.setHasDefault(false);
	interfaceMinMultCKKW.setHasDefault(false);

	}

	namespace ThePEG {

	ostream & operator<<(ostream & os, const HEPEUP & h) {
	os << "<event>\n"
	<< " " << setw(4) << h.NUP
	<< " " << setw(6) << h.IDPRUP
	<< " " << setw(14) << h.XWGTUP
	<< " " << setw(14) << h.SCALUP
	<< " " << setw(14) << h.AQEDUP
	<< " " << setw(14) << h.AQCDUP << "\n";

	for ( int i = 0; i < h.NUP; ++i )
	os << " " << setw(8) << h.IDUP[i]
	<< " " << setw(2) << h.ISTUP[i]
	<< " " << setw(4) << h.MOTHUP[i].first
	<< " " << setw(4) << h.MOTHUP[i].second
	<< " " << setw(4) << h.ICOLUP[i].first
	<< " " << setw(4) << h.ICOLUP[i].second
	<< " " << setw(14) << h.PUP[i][0]
	<< " " << setw(14) << h.PUP[i][1]
	<< " " << setw(14) << h.PUP[i][2]
	<< " " << setw(14) << h.PUP[i][3]
	<< " " << setw(14) << h.PUP[i][4]
	<< " " << setw(1) << h.VTIMUP[i]
	<< " " << setw(1) << h.SPINUP[i] << std::endl;
	os << "</event>" << std::endl;
	return os;
	}

	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:57 PM (1 d, 2 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806183
Default Alt Text: (99 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions