No OneTemporary
Actions

Size

80 KB

Subscribers

None

View Options

	diff --git a/EventRecord/SubProcessGroup.h b/EventRecord/SubProcessGroup.h
	--- a/EventRecord/SubProcessGroup.h
	+++ b/EventRecord/SubProcessGroup.h
	@@ -1,165 +1,165 @@
	// -- C++ --
	//
	// SubProcessGroup.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2010 Simon Platzer
	//
	// 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_SubProcessGroup_H
	#define ThePEG_SubProcessGroup_H
	// This is the declaration of the SubProcessGroup class.

	#include "ThePEG/EventRecord/SubProcess.h"

	namespace ThePEG {

	/**
	* A SubProcessGroup object represents a group of SubProcess
	* objects in dependence of a head SubProcess object.
	*
	* @see StdXCombGroup
	* @see MEGroup
	*/
	class SubProcessGroup: public SubProcess {

	public:

	/**
	* Standard constructor.
	* @param newIncoming the two incoming partons.
	* @param newCollision the Collision to which this SubProcessGroup belongs.
	* @param newHandler the MEBase object which generated this SubProcessGroup.
	*/
	SubProcessGroup(const PPair & newIncoming,
	tCollPtr newCollision = tCollPtr(),
	tcEventBasePtr newHandler = tcEventBasePtr());

	/**
	* Destructor.
	*/
	virtual ~SubProcessGroup();

	/**
	* Return a clone of this sub process group.
	*/
	virtual SubProPtr clone() const;

	protected:

	/**
	* Rebind to cloned objects. When a SubProcessGroup is cloned, a shallow
	* copy is done first, then all <code>Particle</code>s etc, are
	* cloned, and finally this method is used to see to that the
	* pointers in the cloned SubProcessGroup points to the cloned
	* <code>Particle</code>s etc.
	*/
	virtual void rebind(const EventTranslationMap & trans);

	public:

	/**
	* Perform a LorentzTransformation of all particles in the sub
	* process.
	*/
	virtual void transform(const LorentzRotation &);

	/**
	* Return the dependent SubProcess objects
	*/
	const SubProcessVector& dependent() const { return theDependent; }

	/**
	* Access the dependent SubProcess objects
	*/
	SubProcessVector& dependent() { return theDependent; }

	/**
	* Add a dependent SubProcess
	*/
	void add(tSubProPtr sub) { dependent().push_back(sub); }

	public:

	/**
	* Standard function for writing to a persistent stream.
	*/
	void persistentOutput(PersistentOStream &) const;

	/**
	* Standard function for reading from a persistent stream.
	*/
	void persistentInput(PersistentIStream &, int);

	/**
	* Standard Init function. @see Base::Init().
	*/
	static void Init();

	private:

	/**
	* The dependent subprocesses
	*/
	SubProcessVector theDependent;

	public:

	/**
	* Put to ostream
	*/
	virtual void printMe(ostream&) const;

	private:

	/**
	* Describe concrete class with persistent data.
	*/
	static ClassDescription<SubProcessGroup> initSubProcessGroup;

	/**
	* Private default constructor must only be used by the
	* PersistentIStream class via the ClassTraits<SubProcessGroup> class .
	*/
	SubProcessGroup() : SubProcess() {}

	/**
	* The ClassTraits<SubProcessGroup> class must be a friend to be able to
	* use the private default constructor.
	*/
	- friend class ClassTraits<SubProcessGroup>;
	+ friend struct ClassTraits<SubProcessGroup>;

	/**
	* Assignment is forbidden.
	*/
	SubProcessGroup & operator=(const SubProcessGroup &);

	};


	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base class of Collision. */
	template <>
	struct BaseClassTrait<SubProcessGroup,1>: public ClassTraitsType {
	/** Typedef of the first base class of SubProcessGroup. */
	typedef EventRecordBase NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the SubProcessGroup class and how to create it. */
	template <>
	struct ClassTraits<SubProcessGroup>: public ClassTraitsBase<SubProcessGroup> {
	/** Return a platform-independent class name */
	static string className() { return "ThePEG::SubProcessGroup"; }
	/** Create a SubProcessGroup object. */
	static TPtr create() { return TPtr::Create(SubProcessGroup()); }
	};

	/** @endcond */

	}

	#endif /* ThePEG_SubProcessGroup_H */
	diff --git a/LesHouches/LesHouchesReader.cc b/LesHouches/LesHouchesReader.cc
	--- a/LesHouches/LesHouchesReader.cc
	+++ b/LesHouches/LesHouchesReader.cc
	@@ -1,1510 +1,1510 @@
	// -- 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),
	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));

	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 && 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"
	<< 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;
	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);
	} else {
	id = idit - lprup.begin();
	}
	++neve;
	++oldeve[id];
	oldsum += hepeup.XWGTUP;
	if ( cacheFile() ) {
	if ( eventWeight() == 0.0 ) {
	++cuteve;
	continue;
	}
	cacheEvent();
	}
	++neweve[id];
	newmax[id] = max(newmax[id], abs(eventWeight()));
	if ( eventWeight() < 0.0 ) negw = true;
	}

	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 ) {
	// 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];
	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 ( 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;

	}

	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);
	if(theMomentumTreatment == 1) mom.rescaleEnergy();
	else if(theMomentumTreatment == 2) mom.rescaleMass();
	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;
	}
	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;
	}
	if( theIncludeSpin && abs(pd->id()) == ParticleID::tauminus &&
	hepeup.SPINUP[i] !=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.-hepeup.SPINUP[i]);
	rho(1,1) = 0.5*(1.+hepeup.SPINUP[i]);
	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]) {
	+ 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]) {
	+ 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) {
	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]) {
	+ 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(iy==0) {
	if(external[iy]->colourInfo()->colourLines().empty()) continue;
	anti2 = external[iy]->colourInfo()->antiColourLines();
	}
	else {
	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]) {
	+ 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 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);
	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);
	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
	<< 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
	>> 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 form 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;
	}

	}

	diff --git a/MatrixElement/MEqq2qq1.h b/MatrixElement/MEqq2qq1.h
	--- a/MatrixElement/MEqq2qq1.h
	+++ b/MatrixElement/MEqq2qq1.h
	@@ -1,161 +1,161 @@
	// -- C++ --
	//
	// MEqq2qq1.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_MEqq2qq_H
	#define ThePEG_MEqq2qq_H
	// This is the declaration of the MEqq2qq class.

	#include "ThePEG/MatrixElement/ME2to2QCD.h"

	namespace ThePEG {

	/**
	* MEqq2qq inherits from the ME2to2QCD and implements the standard
	* \f$q_i\bar{q}_i\rightarrow q_i\bar{q}_i\f$ matrix element.
	*
	* @see \ref MEqq2qqInterfaces "The interfaces"
	* defined for MEqq2qq.
	* @see ME2to2QCD
	*/
	class MEqq2qq: public ME2to2QCD {

	public:

	/** @name Virtual functions required by the MEBase class. */
	//@{
	/**
	* The matrix element for the kinematical configuration
	* previously provided by the last call to setKinematics(), suitably
	* scaled by sHat() to give a dimension-less number.
	* @return the matrix element scaled with sHat() to give a
	* dimensionless number.
	*/
	virtual double me2() const;

	/**
	* Add all possible diagrams with the add() function.
	*/
	virtual void getDiagrams() const;

	/**
	* Return a Selector with possible colour geometries for the selected
	* diagram weighted by their relative probabilities.
	* @param diag the diagram chosen.
	* @return the possible colour geometries weighted by their
	* relative probabilities.
	*/
	virtual Selector<const ColourLines *>
	colourGeometries(tcDiagPtr diag) const;

	/**
	* Get diagram selector. With the information previously supplied with the
	* setKinematics method, a derived class may optionally
	* override this method to weight the given diagrams with their
	* (although certainly not physical) relative probabilities.
	* @param dv the diagrams to be weighted.
	* @return a Selector relating the given diagrams to their weights.
	*/
	virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;
	//@}

	protected:

	/** @name Internal functions returning the matrix element squared
	* for different colour configurations. */
	//@{
	/**
	* Return the matrix element squared (without common pre-factors)
	* for the specific colour configuration.
	*/
	double colA() const
	{
	return (sqr(tHat()) + sqr(uHat()))/sqr(sHat()) +
	- interference()? -double(sqr(uHat())/(3.0sHat()tHat())): 0.0;
	+ (interference()? -double(sqr(uHat())/(3.0sHat()tHat())): 0.0);
	}

	/**
	* Return the matrix element squared (without common pre-factors)
	* for the specific colour configuration.
	*/
	double colB() const
	{
	return (sqr(uHat()) + sqr(sHat()))/sqr(tHat()) +
	- interference()? -double(sqr(uHat())/(3.0sHat()tHat())): 0.0;
	+ (interference()? -double(sqr(uHat())/(3.0sHat()tHat())): 0.0);
	}
	//@}

	public:

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	private:

	/**
	* Describe a concrete class without persistent data.
	*/
	static NoPIOClassDescription<MEqq2qq> initMEqq2qq;

	/**
	* Private and non-existent assignment operator.
	*/
	MEqq2qq & operator=(const MEqq2qq &);

	};

	}


	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of MEqq2qq. */
	template <>
	struct BaseClassTrait<MEqq2qq,1>: public ClassTraitsType {
	/** Typedef of the first base class of MEqq2qq. */
	typedef ME2to2QCD NthBase;
	};

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

	/** @endcond */

	}

	#endif /* ThePEG_MEqq2qq_H */
	diff --git a/Vectors/HepMCConverter.tcc b/Vectors/HepMCConverter.tcc
	--- a/Vectors/HepMCConverter.tcc
	+++ b/Vectors/HepMCConverter.tcc
	@@ -1,315 +1,315 @@
	// -- C++ --
	//
	// HepMCConverter.tcc 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 HepMCConverter class.
	//

	#include "HepMCConverter.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/StandardSelectors.h"
	#include "ThePEG/EventRecord/Collision.h"
	#include "ThePEG/EventRecord/Step.h"
	#include "ThePEG/EventRecord/SubProcess.h"
	#include "ThePEG/Handlers/XComb.h"
	#include "ThePEG/Handlers/EventHandler.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/PDF/PDF.h"
	#include "ThePEG/PDT/StandardMatchers.h"
	#include "ThePEG/Utilities/Throw.h"

	namespace ThePEG {


	template <typename HepMCEventT, typename Traits>
	typename HepMCConverter<HepMCEventT,Traits>::GenEvent *
	HepMCConverter<HepMCEventT,Traits>::
	convert(const Event & ev, bool nocopies, Energy eunit, Length lunit) {
	HepMCConverter<HepMCEventT,Traits> converter(ev, nocopies, eunit, lunit);
	return converter.geneve;
	}

	template <typename HepMCEventT, typename Traits>
	void HepMCConverter<HepMCEventT,Traits>::
	convert(const Event & ev, GenEvent & gev, bool nocopies) {
	HepMCConverter<HepMCEventT,Traits>
	converter(ev, gev, nocopies,
	Traits::momentumUnit(gev), Traits::lengthUnit(gev));
	}

	template <typename HepMCEventT, typename Traits>
	void HepMCConverter<HepMCEventT,Traits>::
	convert(const Event & ev, GenEvent & gev, bool nocopies,
	Energy eunit, Length lunit) {
	HepMCConverter<HepMCEventT,Traits> converter(ev, gev, nocopies, eunit, lunit);
	}

	template <typename HepMCEventT, typename Traits>
	HepMCConverter<HepMCEventT,Traits>::
	HepMCConverter(const Event & ev, bool nocopies, Energy eunit, Length lunit)
	: energyUnit(eunit), lengthUnit(lunit) {

	geneve = Traits::newEvent(ev.number(), ev.weight(), ev.optionalWeights());

	init(ev, nocopies);

	}

	template <typename HepMCEventT, typename Traits>
	HepMCConverter<HepMCEventT,Traits>::
	HepMCConverter(const Event & ev, GenEvent & gev, bool nocopies,
	Energy eunit, Length lunit)
	: energyUnit(eunit), lengthUnit(lunit) {

	geneve = &gev;
	Traits::resetEvent(geneve, ev.number(), ev.weight(), ev.optionalWeights());

	init(ev, nocopies);

	}

	struct ParticleOrderNumberCmp {
	bool operator()(tcPPtr a, tcPPtr b) const {
	return a->number() < b->number();
	}
	};

	template <typename HepMCEventT, typename Traits>
	void HepMCConverter<HepMCEventT,Traits>::init(const Event & ev, bool nocopies) {

	if ( lengthUnit != millimeter && lengthUnit != centimeter )
	Throw<HepMCConverterException>()
	<< "Length unit used for HepMC::GenEvent was not MM nor CM."
	<< Exception::runerror;
	if ( energyUnit != GeV && energyUnit != MeV )
	Throw<HepMCConverterException>()
	<< "Momentum unit used for HepMC::GenEvent was not GEV nor MEV."
	<< Exception::runerror;
	Traits::setUnits(*geneve, energyUnit, lengthUnit);

	tcEHPtr eh;
	if ( ev.primaryCollision() && ( eh =
	dynamic_ptr_cast<tcEHPtr>(ev.primaryCollision()->handler()) ) ) {

	// Get general event info if present.
	Traits::setScaleAndAlphas(*geneve, eh->lastScale(),
	eh->SM().alphaS(eh->lastScale()),
	eh->SM().alphaEM(eh->lastScale()), energyUnit);
	}

	// Extract all particles and order them.
	tcPVector all;
	ev.select(back_inserter(all), SelectAll());
	stable_sort(all.begin(), all.end(), ParticleOrderNumberCmp());
	vertices.reserve(all.size()*2);

	// Create GenParticle's and map them to the ThePEG particles.
	for ( int i = 0, N = all.size(); i < N; ++i ) {
	tcPPtr p = all[i];
	if ( nocopies && p->next() ) continue;
	if ( pmap.find(p) != pmap.end() ) continue;
	GenParticle * gp = pmap[p] = createParticle(p);
	if ( p->hasColourInfo() ) {
	// Check if the particle is connected to colour lines, in which
	// case the lines are mapped to an integer and set in the
	// GenParticle's Flow info.
	tcColinePtr l;
	- if ( l = p->colourLine() ) {
	+ if ( (l = p->colourLine()) ) {
	if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500;
	Traits::setColourLine(*gp, 1, flowmap[l]);
	}
	- if ( l = p->antiColourLine() ) {
	+ if ( (l = p->antiColourLine()) ) {
	if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500;
	Traits::setColourLine(*gp, 2, flowmap[l]);
	}
	}

	if ( !p->children().empty() \|\| p->next() ) {
	// If the particle has children it should have a decay vertex:
	vertices.push_back(Vertex());
	decv[p] = &vertices.back();
	vertices.back().in.insert(p);
	}

	if ( !p->parents().empty() \|\| p->previous() \|\|
	(p->children().empty() && !p->next()) ) {
	// If the particle has parents it should have a production
	// vertex. If neither parents or children it should still have a
	// dummy production vertex.
	vertices.push_back(Vertex());
	prov[p] = &vertices.back();
	vertices.back().out.insert(p);
	}
	}

	// Now go through the the particles again, and join the vertices.
	for ( int i = 0, N = all.size(); i < N; ++i ) {
	tcPPtr p = all[i];
	if ( nocopies ) {
	if ( p->next() ) continue;
	for ( int i = 0, N = p->children().size(); i < N; ++i )
	join(p, p->children()[i]->final());
	tcPPtr pp = p;
	while ( pp->parents().empty() && pp->previous() ) pp = pp->previous();
	for ( int i = 0, N = pp->parents().size(); i < N; ++i )
	join(pp->parents()[i]->final(), p);
	} else {
	for ( int i = 0, N = p->children().size(); i < N; ++i )
	join(p, p->children()[i]);
	if ( p->next() ) join(p, p->next());
	for ( int i = 0, N = p->parents().size(); i < N; ++i )
	join(p->parents()[i], p);
	if ( p->previous() ) join(p->previous(), p);
	}
	}

	// Time to create the GenVertex's
	for ( typename VertexMap::iterator it = prov.begin(); it != prov.end(); ++it )
	if ( !member(vmap, it->second) )
	vmap[it->second] = createVertex(it->second);
	for ( typename VertexMap::iterator it = decv.begin(); it != decv.end(); ++it )
	if ( !member(vmap, it->second) )
	vmap[it->second] = createVertex(it->second);

	// Now find the primary signal process vertex defined to be the
	// decay vertex of the first parton coming into the primary hard
	// sub-collision.
	tSubProPtr sub = ev.primarySubProcess();
	if ( sub && sub->incoming().first ) {
	const Vertex * prim = decv[sub->incoming().first];
	Traits::setSignalProcessVertex(*geneve, vmap[prim]);
	vmap.erase(prim);
	}

	// Then add the rest of the vertices.
	for ( typename GenVertexMap::iterator it = vmap.begin();
	it != vmap.end(); ++it )
	Traits::addVertex(*geneve, it->second);

	// and the incoming beam particles
	Traits::setBeamParticles(*geneve,pmap[ev.incoming().first],
	pmap[ev.incoming().second]);

	// and the PDF info
	setPdfInfo(ev);

	// and the cross section info
	Traits::setCrossSection(*geneve,
	eh->integratedXSec()/picobarn,
	eh->integratedXSecErr()/picobarn);

	}

	template <typename HepMCEventT, typename Traits>
	typename HepMCConverter<HepMCEventT,Traits>::GenParticle *
	HepMCConverter<HepMCEventT,Traits>::createParticle(tcPPtr p) const {
	int status = 1;
	if ( !p->children().empty() \|\| p->next() ) status = 11;
	if ( !p->children().empty() ) {
	long id = p->data().id();
	if ( BaryonMatcher::Check(id) \|\| MesonMatcher::Check(id) \|\|
	id == ParticleID::muminus \|\| id == ParticleID::muplus \|\|
	id == ParticleID::tauminus \|\| id == ParticleID::tauplus )
	if ( p->mass() <= p->data().massMax() &&
	p->mass() >= p->data().massMin() ) status = 2;
	}
	GenParticle * gp =
	Traits::newParticle(p->momentum(), p->id(), status, energyUnit);

	if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) {
	DPair pol = p->spinInfo()->polarization();
	Traits::setPolarization(*gp, pol.first, pol.second);
	}

	return gp;

	}

	template <typename HepMCEventT, typename Traits>
	void HepMCConverter<HepMCEventT,Traits>::join(tcPPtr parent, tcPPtr child) {
	Vertex * dec = decv[parent];
	Vertex * pro = prov[child];
	if ( !pro \|\| !dec ) Throw<HepMCConverterException>()
	<< "Found a reference to a ThePEG::Particle which was not in the Event."
	<< Exception::eventerror;
	if ( pro == dec ) return;
	while ( !pro->in.empty() ) {
	dec->in.insert(*(pro->in.begin()));
	decv[*(pro->in.begin())] = dec;
	pro->in.erase(pro->in.begin());
	}
	while ( !pro->out.empty() ) {
	dec->out.insert(*(pro->out.begin()));
	prov[*(pro->out.begin())] = dec;
	pro->out.erase(pro->out.begin());
	}
	}

	template <typename HepMCEventT, typename Traits>
	typename HepMCConverter<HepMCEventT,Traits>::GenVertex *
	HepMCConverter<HepMCEventT,Traits>::createVertex(Vertex * v) {
	if ( !v ) Throw<HepMCConverterException>()
	<< "Found internal null Vertex." << Exception::abortnow;

	GenVertex * gv = new GenVertex();

	// We assume that the vertex position is the average of the decay
	// vertices of all incoming and the creation vertices of all
	// outgoing particles in the lab. Note that this will probably not
	// be useful information for very small distances.
	LorentzPoint p;
	for ( tcParticleSet::iterator it = v->in.begin();
	it != v->in.end(); ++it ) {
	p += (**it).labDecayVertex();
	Traits::addIncoming(gv, pmap[it]);
	}
	for ( tcParticleSet::iterator it = v->out.begin();
	it != v->out.end(); ++it ) {
	p += (**it).labVertex();
	Traits::addOutgoing(gv, pmap[it]);
	}

	p /= double(v->in.size() + v->out.size());
	Traits::setPosition(*gv, p, lengthUnit);

	return gv;
	}

	template <typename HepMCEventT, typename Traits>
	void HepMCConverter<HepMCEventT,Traits>::setPdfInfo(const Event & e) {
	// ids of the partons going into the primary sub process
	tSubProPtr sub = e.primarySubProcess();
	int id1 = sub->incoming().first ->id();
	int id2 = sub->incoming().second->id();
	// get the event handler
	tcEHPtr eh = dynamic_ptr_cast<tcEHPtr>(e.handler());
	// get the values of x
	double x1 = eh->lastX1();
	double x2 = eh->lastX2();
	// get the pdfs
	pair<PDF,PDF> pdfs;
	pdfs.first = eh->pdf<PDF>(sub->incoming().first );
	pdfs.second = eh->pdf<PDF>(sub->incoming().second);
	// get the scale
	Energy2 scale = eh->lastScale();
	// get the values of the pdfs
	double xf1 = pdfs.first.xfx(sub->incoming().first->dataPtr(), scale, x1);
	double xf2 = pdfs.second.xfx(sub->incoming().second->dataPtr(), scale, x2);

	Traits::setPdfInfo(*geneve, id1, id2, x1, x2, sqrt(scale/GeV2), xf1, xf2);

	}



	}


	diff --git a/configure.ac b/configure.ac
	--- a/configure.ac
	+++ b/configure.ac
	@@ -1,150 +1,153 @@
	dnl Process this file with autoconf to produce a configure script.

	AC_PREREQ([2.59])
	AC_INIT([ThePEG],[1.8.0],[http://www.thep.lu.se/ThePEG/],[ThePEG])
	AC_CONFIG_AUX_DIR([Config])
	AC_CONFIG_MACRO_DIR([m4])

	THEPEG_LIBTOOL_VERSION_INFO(15,0,0)

	AC_CONFIG_SRCDIR([EventRecord/SubProcess.h])
	AC_CONFIG_HEADERS([Config/config.h Config/LWH.h Config/ThePEG_Qty.h:Config/ThePEG_Qty.h.a:Config/ThePEG_Qty.h.in:Config/ThePEG_Qty.h.b])

	AC_CANONICAL_HOST

	case "${host}" in
	-darwin[[0156]].)
	AC_MSG_ERROR([ThePEG requires OS X 10.3 or later])
	;;
	-darwin7.)
	if test "x$MACOSX_DEPLOYMENT_TARGET" != "x10.3"; then
	AC_MSG_ERROR(
	[Please export the MACOSX_DEPLOYMENT_TARGET variable, set to 10.3])
	fi
	;;
	esac

	AC_LANG(C++)
	AM_INIT_AUTOMAKE([1.9 gnu dist-bzip2 -Wall])
	dnl also include std-options once --version and --help exist
	m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
	m4_ifdef([AM_PROG_AR], [AM_PROG_AR])

	dnl Checks for programs.
	AC_PROG_CXX
	AC_PROG_INSTALL
	AC_PROG_MAKE_SET
	AC_PROG_LN_S

	LT_PREREQ([2.2])
	LT_INIT([disable-static dlopen pic-only])

	VL_LIB_READLINE

	THEPEG_UNIT_CHECKING
	THEPEG_CHECK_GSL
	THEPEG_SEARCH_LHAPDF
	THEPEG_CHECK_HEPMC
	THEPEG_CHECK_RIVET

	AX_CHECK_ZLIB

	dnl AX_CHECK_BZ2LIB

	THEPEG_DEFINE_ENVDEFAULT(ThePEG_GZREAD_FILE,GZREAD_FILE,gunzip -c,[The command which, taking the name of a gzipped file as argument, unzips it and prints it to stdout. Default is "gunzip -c"])

	THEPEG_DEFINE_ENVDEFAULT(ThePEG_GZWRITE_FILE,GZWRITE_FILE,[gzip -c > ],[The command which, taking the name of a gzipped file as argument, reads stdin, zips it and writes it to the file. Default is "gzip -c > ".])

	THEPEG_DEFINE_ENVDEFAULT(ThePEG_BZ2READ_FILE,BZ2READ_FILE,bunzip2 -c,[The command which, taking the name of a bzipped file as argument, unzips it and prints it to stdout. Default is "bunzip2 -c".])

	THEPEG_DEFINE_ENVDEFAULT(ThePEG_BZ2WRITE_FILE,BZ2WRITE_FILE,[bzip2 -c > ],[The command which, taking the name of a bzipped file as argument, reads stdin, zips it and writes it to the file. Default is "bzip2 -c > ".])

	THEPEG_CHECK_EXPM1
	THEPEG_CHECK_LOG1P
	THEPEG_CHECK_DLOPEN
	THEPEG_CHECK_FPUCONTROL
	THEPEG_CHECK_FENV
	THEPEG_CHECK_AIDA

	AM_CPPFLAGS="-I\$(top_builddir)/include \$(GSLINCLUDE)"

	AX_COMPILER_VENDOR

	case "${ax_cv_cxx_compiler_vendor}" in
	gnu)
	AM_CXXFLAGS="-ansi -pedantic -Wall -W"
	;;
	+ clang)
	+ AM_CXXFLAGS="-ansi -pedantic -Wall"
	+ ;;
	intel)
	AM_CXXFLAGS="-strict-ansi -Wall -wd13000,1418,981,444,383,1599,1572,2259,980"
	;;
	esac

	THEPEG_CHECK_ABS_BUG

	AC_SUBST(AM_CPPFLAGS)
	AC_SUBST(AM_CXXFLAGS)

	dnl do an actual capability check on ld instead of this workaround
	case "${host}" in
	-darwin)
	;;
	*)
	AM_LDFLAGS="-Wl,--enable-new-dtags"
	;;
	esac

	AC_SUBST(AM_LDFLAGS)

	THEPEG_EMPTY_SUBST

	AC_PATH_PROG(PERL, perl)

	AC_ARG_WITH(javagui,
	[ --with-javagui Compile and install the java-based GUI.])

	if test "x$with_javagui" != "xno"; then
	THEPEG_HAS_JAVA([1.4], [], [with_javagui=no; AC_MSG_NOTICE([Java GUI disabled])])
	fi

	AM_CONDITIONAL([JAVAGUI], [test "x$with_javagui" != "xno"])

	AC_CONFIG_FILES([Helicity/Makefile
	Helicity/WaveFunction/Makefile
	Helicity/Vertex/Makefile
	Helicity/Vertex/Scalar/Makefile
	Helicity/Vertex/Vector/Makefile
	Helicity/Vertex/Tensor/Makefile
	Utilities/Makefile
	include/Makefile
	Interface/Makefile
	LesHouches/Makefile
	Vectors/Makefile
	PDT/Makefile
	PDF/Makefile
	Persistency/Makefile
	Config/Makefile
	Handlers/Makefile
	MatrixElement/Makefile
	Pointer/Makefile
	lib/Makefile
	lib/Makefile.common.install
	src/Makefile
	ACDC/Makefile
	Repository/Makefile
	EventRecord/Makefile
	StandardModel/Makefile
	Cuts/Makefile
	Analysis/Makefile
	Doc/Makefile
	Doc/MakeDocs.in
	Doc/refman.h
	Doc/refman.conf
	java/Makefile
	Makefile])
	AC_CONFIG_FILES([Doc/fixinterfaces.pl],[chmod +x Doc/fixinterfaces.pl])

	THEPEG_OVERVIEW

	AC_CONFIG_COMMANDS([summary],[cat config.thepeg])

	AC_OUTPUT

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 12:46 PM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4017320
Default Alt Text: (80 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions