No OneTemporary
Actions

Size

156 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/FxFx/FxFxAnalysis.cc b/MatrixElement/FxFx/FxFxAnalysis.cc
	--- a/MatrixElement/FxFx/FxFxAnalysis.cc
	+++ b/MatrixElement/FxFx/FxFxAnalysis.cc
	@@ -1,389 +1,389 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FxFxAnalysis class.
	//

	#include "FxFxAnalysis.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Vectors/HepMCConverter.h"
	#include "ThePEG/Config/HepMCHelper.h"
	#include "ThePEG/EventRecord/Event.h"
	#include "ThePEG/EventRecord/SubProcess.h"
	#include "ThePEG/EventRecord/SubProcessGroup.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "HepMC/GenEvent.h"
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/Tools/Logging.hh"
	#include "HepMC/IO_AsciiParticles.h"
	#include "HepMC/IO_GenEvent.h"
	#include "ThePEG/Utilities/XSecStat.h"

	using namespace ThePEG;

	FxFxAnalysis::FxFxAnalysis()
	: _remnantId(82), _format(1),_unitchoice(),
	- _geneventPrecision(16), debug(false), _rivet(), _nevent(0), useoptweights(false), normoptweights(false) {}
	+ _geneventPrecision(16), debug(false),useoptweights(false), normoptweights(false), _rivet(), _nevent(0) {}

	HepMC::GenEvent * FxFxAnalysis::makeEvent(tEventPtr event, tSubProPtr sub, long no,
	Energy eUnit, Length lUnit,
	CrossSection xsec, CrossSection xsecErr) const {
	//convert the event from the Herwig format to the HepMC format and write it to the common block
	HepMC::GenEvent * ev = HepMCConverter<HepMC::GenEvent>::convert(*event, false,eUnit, lUnit);

	//reset the event
	HepMCTraits<HepMC::GenEvent>::resetEvent(ev, no, event->weight()*sub->groupWeight(), event->optionalWeights());

	//set the cross section
	HepMCTraits<HepMC::GenEvent>::setCrossSection(*ev,xsec/picobarn,
	xsecErr/picobarn);
	return ev;
	}

	-HepMC::GenEvent * FxFxAnalysis::makeEventW(tEventPtr event, tSubProPtr sub, long no,
	+HepMC::GenEvent * FxFxAnalysis::makeEventW(tEventPtr event, long no,
	Energy eUnit, Length lUnit,
	CrossSection xsec, CrossSection xsecErr, double evoptweight, double centralweight) const {

	//convert the event from the Herwig format to the HepMC format and write it to the common block
	HepMC::GenEvent * ev = HepMCConverter<HepMC::GenEvent>::convert(*event, false,eUnit, lUnit);

	if(normoptweights) { evoptweight /= centralweight; }

	//reset the event
	HepMCTraits<HepMC::GenEvent>::resetEvent(ev, no, evoptweight, event->optionalWeights());
	//set the cross section
	HepMCTraits<HepMC::GenEvent>::setCrossSection(*ev,xsec/picobarn,
	xsecErr/picobarn);
	return ev;
	}

	void FxFxAnalysis::analyze(ThePEG::tEventPtr event, long ieve, int loop, int state) {
	Energy eUnit;
	Length lUnit;
	switch (_unitchoice) {
	default: eUnit = GeV; lUnit = millimeter; break;
	case 1: eUnit = MeV; lUnit = millimeter; break;
	case 2: eUnit = GeV; lUnit = centimeter; break;
	case 3: eUnit = MeV; lUnit = centimeter; break;
	}

	tcFxFxEventHandlerPtr eh = dynamic_ptr_cast<tcFxFxEventHandlerPtr>(event->primaryCollision()->handler());
	assert(eh);

	CrossSection xsec = eh->integratedXSec();
	CrossSection xsecErr = eh->integratedXSecErr();

	optxsec = eh->optintegratedXSecMap();


	int ii = 0;
	if(useoptweights) {
	for (map<string,CrossSection>::const_iterator it= optxsec.begin(); it!=optxsec.end(); ++it){
	OptXS[ii] = it->second/picobarn;
	ii++;
	}
	}

	tSubProPtr sub = event->primarySubProcess();
	Ptr<SubProcessGroup>::tptr grp =
	dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);

	AnalysisHandler::analyze(event, ieve, loop, state);
	// Rotate to CMS, extract final state particles and call analyze(particles).
	// convert to hepmc
	HepMC::GenEvent * hepmc =
	makeEvent(event,sub,_nevent,eUnit,lUnit,xsec,xsecErr);

	//count weights here
	std::vector<std::string> strs;
	if(_numweights == -999) {
	_numweights = 0;
	for (map<string,double>::const_iterator it= event->optionalWeights().begin(); it!=event->optionalWeights().end(); ++it){
	string first_piece = it->first;
	string word;
	istringstream iss(first_piece, istringstream::in);
	while( iss >> word ) strs.push_back(word);
	if(strs[0] != "np") { _numweights++; }
	strs.clear();
	}
	}
	double CentralWeight = 1.;
	if(useoptweights) {
	//find the weights
	_i = 0;//counter
	OptWeights.clear();
	for (map<string,double>::const_iterator it= event->optionalWeights().begin(); it!=event->optionalWeights().end(); ++it){
	// std::cout << it->first << " => " << it->second << '\n';
	string first_piece = it->first;
	string word;
	istringstream iss(first_piece, istringstream::in);
	while( iss >> word ) {
	strs.push_back(word);
	}
	std::pair<int,double> OptWeightsTemp;
	/* if(strs[0] == "PDF") {
	OptWeightsTemp.first = atoi(strs[2].c_str());
	OptWeightsTemp.second= it->second;
	OptWeights.push_back(OptWeightsTemp);
	_i++;
	}
	if(strs[0] == "SC") {
	OptWeightsTemp.first = atoi(strs[3].c_str());
	OptWeightsTemp.second = it->second;
	// cout << "OptWeightsTemp.first = " << OptWeightsTemp.first << " OptWeightsTemp.second = " << OptWeightsTemp.second << endl;
	if(OptWeightsTemp.first == 1001) { cout << "OptWeightsTemp.second = " << OptWeightsTemp.second << endl; CentralWeight = OptWeightsTemp.second; }
	OptWeights.push_back(OptWeightsTemp);
	_i++;
	}*/
	if(strs[0] != "np") {
	OptWeightsTemp.first = atoi((it->first).c_str());
	OptWeightsTemp.second= it->second;
	if(OptWeightsTemp.first == 1001) { /cout << "OptWeightsTemp.second = " << OptWeightsTemp.second << endl;/ CentralWeight = OptWeightsTemp.second; }
	OptWeights.push_back(OptWeightsTemp);
	_i++;
	}
	strs.clear();
	}
	}

	/* multiple hepmcs for scale/pdf variations
	*/
	vector<HepMC::GenEvent*> hepmcMULTI;// = new HepMC::GenEvent[_numweights];
	HepMC::GenEvent * hepmcMULTIi;
	if(useoptweights) {
	for(int rr = 0; rr < _numweights; rr++) {
	double xsrr = optxsec[std::to_string(OptWeights[rr].first)]/picobarn;
	/* cout << "xsec = " << xsec/picobarn << endl;
	cout << "OptWeights[rr].second = " << OptWeights[rr].second << endl;
	cout << "xsrr = " << xsrr << endl;*/
	- hepmcMULTIi = makeEventW(event,sub,_nevent,eUnit,lUnit,xsrr*picobarn,xsecErr,OptWeights[rr].second, CentralWeight);
	+ hepmcMULTIi = makeEventW(event,_nevent,eUnit,lUnit,xsrr*picobarn,xsecErr,OptWeights[rr].second, CentralWeight);
	hepmcMULTI.push_back(hepmcMULTIi);
	}
	}

	CurrentGenerator::Redirect stdout(cout);

	if ( _rivet ) _rivet->analyze(*hepmc);

	if(useoptweights) {
	for(int rr = 0; rr < _numweights; rr++) {
	if ( _rivetMULTI[rr] ) _rivetMULTI[rr]->analyze(*hepmcMULTI[rr]);
	}
	}
	// delete hepmc events
	delete hepmc;
	if(useoptweights) {
	for(int rr = 0; rr < _numweights; rr++) {
	delete hepmcMULTI[rr];
	}
	}

	if ( grp ) {
	for ( SubProcessVector::const_iterator s = grp->dependent().begin();
	s != grp->dependent().end(); ++s ) {

	hepmc = makeEvent(event,*s,_nevent,eUnit,lUnit,xsec,xsecErr);

	if ( _rivet ) _rivet->analyze(*hepmc);
	// delete hepmc event
	delete hepmc;


	}

	}

	++_nevent;
	}


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

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

	void FxFxAnalysis::persistentOutput(ThePEG::PersistentOStream & os) const {
	os << _analyses << filename << debug << useoptweights << normoptweights;
	}

	void FxFxAnalysis::persistentInput(ThePEG::PersistentIStream & is, int) {
	is >> _analyses >> filename >> debug >> useoptweights >> normoptweights;
	}

	ThePEG::ClassDescription<FxFxAnalysis> FxFxAnalysis::initFxFxAnalysis;
	// Definition of the static class description member.

	void FxFxAnalysis::Init() {
	static ThePEG::ClassDocumentation<FxFxAnalysis> documentation
	("The FxFxAnalysis class is a simple class to allow analyses"
	" from the Rivet library to be called from ThePEG");

	static ThePEG::ParVector<FxFxAnalysis,string> interfaceAnalyses
	("Analyses",
	"The names of the Rivet analyses to use",
	&FxFxAnalysis::_analyses, -1, "", "","" "",
	false, false, ThePEG::Interface::nolimits);

	static Parameter<FxFxAnalysis,long> interfaceRemnantId
	("RemnantId",
	"Set the PDG id to be used for remnants.",
	&FxFxAnalysis::_remnantId, 82, 0, 0,
	false, false, Interface::nolimits);

	static Parameter<FxFxAnalysis,string> interfaceFilename
	("Filename",
	"The name of the file where the YODA histograms are put. If empty, "
	"the run name will be used instead. '.yoda' will in any case be "
	"appended to the file name.",
	&FxFxAnalysis::filename, "", true, false);


	static Switch<FxFxAnalysis,bool> interfaceDebug
	("Debug",
	"Enable debug information from Rivet",
	&FxFxAnalysis::debug, false, true, false);
	static SwitchOption interfaceDebugNo
	(interfaceDebug,
	"No",
	"Disable debug information.",
	false);
	static SwitchOption interfaceDebugYes
	(interfaceDebug,
	"Yes",
	"Enable debug information from Rivet.",
	true);


	static Switch<FxFxAnalysis,bool> interfaceUseOptWeights
	("UseOptWeights",
	"Enable debug information from Rivet",
	&FxFxAnalysis::useoptweights, false, true, false);
	static SwitchOption interfaceUseOptWeightsNo
	(interfaceUseOptWeights,
	"No",
	"Disable optional weights",
	false);
	static SwitchOption interfaceUseOptWeightsYes
	(interfaceUseOptWeights,
	"Yes",
	"Enable debug information from Rivet.",
	true);

	static Switch<FxFxAnalysis,bool> interfaceNormOptWeights
	("NormOptWeights",
	"Enable debug information from Rivet",
	&FxFxAnalysis::normoptweights, false, true, false);
	static SwitchOption interfaceNormOptWeightsNo
	(interfaceNormOptWeights,
	"No",
	"Do not normalize optional weights",
	false);
	static SwitchOption interfacedNormOptWeightsYes
	(interfaceNormOptWeights,
	"Yes",
	"Normalize optional weights.",
	true);




	interfaceAnalyses.rank(10);

	}

	void FxFxAnalysis::dofinish() {
	AnalysisHandler::dofinish();
	if( _nevent > 0 && _rivet ) {
	CurrentGenerator::Redirect stdout(cout);
	_rivet->setCrossSection(generator()->integratedXSec()/picobarn);
	_rivet->finalize();

	string fname = filename;
	fname = generator()->runName() + ".yoda";
	_rivet->writeData(fname);
	}

	delete _rivet;
	_rivet = 0;

	if(useoptweights) {
	for(int rr = 0; rr < _numweights; rr++) {
	cout << (OptWeights[rr].first) << ", cross section = " << OptXS[rr] << endl;
	if( _nevent > 0 && _rivetMULTI[rr] ) {
	double xsrr = optxsec[std::to_string(OptWeights[rr].first)]/picobarn;
	// _rivetMULTI[rr]->setCrossSection(OptXS[rr]);
	_rivetMULTI[rr]->setCrossSection(xsrr);
	_rivetMULTI[rr]->finalize();

	string fname = filename;
	fname = generator()->runName() + "_" + std::to_string(OptWeights[rr].first) + ".yoda";

	_rivetMULTI[rr]->writeData(fname);
	}
	delete _rivetMULTI[rr];
	_rivetMULTI[rr] = 0;
	}
	}
	}

	void FxFxAnalysis::doinit() {
	_numweights = -999;

	for(int rr = 0; rr < 120; rr++) {
	OptXS.push_back(0.);
	}
	AnalysisHandler::doinit();
	if(_analyses.empty())
	throw ThePEG::Exception() << "Must have at least one analysis loaded in "
	<< "FxFxAnalysis::doinitrun()"
	<< ThePEG::Exception::runerror;

	// check that analysis list is available
	_rivet = new Rivet::AnalysisHandler; //(fname);
	_rivet->addAnalyses(_analyses);
	if ( _rivet->analysisNames().size() != _analyses.size() ) {
	throw ThePEG::Exception()
	<< "Rivet could not find all requested analyses.\n"
	<< "Use 'rivet --list-analyses' to check availability.\n"
	<< ThePEG::Exception::runerror;
	}

	delete _rivet;
	_rivet = 0;
	}

	void FxFxAnalysis::doinitrun() {
	_numweights = -999;
	AnalysisHandler::doinitrun();
	// create FxFx analysis handler
	if(useoptweights) { cout << "Warning: Using optional weights launches multiple rivet analyses. This may slow down your run substantially!" << endl; }
	CurrentGenerator::Redirect stdout(cout);
	_rivet = new Rivet::AnalysisHandler; //(fname);
	_rivet->addAnalyses(_analyses);
	if(useoptweights) {
	for(int rr = 0; rr < 110; rr++) {
	OptXS.push_back(0.);
	_rivetMULTI[rr] = new Rivet::AnalysisHandler; //(fname);
	_rivetMULTI[rr]->addAnalyses(_analyses);
	}
	}

	// check that analysis list is still available
	if ( _rivet->analysisNames().size() != _analyses.size() ) {
	throw ThePEG::Exception()
	<< "Rivet could not find all requested analyses.\n"
	<< "Use 'rivet --list-analyses' to check availability.\n"
	<< ThePEG::Exception::runerror;
	}
	if ( debug )
	Rivet::Log::setLevel("Rivet",Rivet::Log::DEBUG);
	}
	diff --git a/MatrixElement/FxFx/FxFxAnalysis.h b/MatrixElement/FxFx/FxFxAnalysis.h
	--- a/MatrixElement/FxFx/FxFxAnalysis.h
	+++ b/MatrixElement/FxFx/FxFxAnalysis.h
	@@ -1,262 +1,262 @@
	// -- C++ --
	#ifndef THEPEG_FxFxAnalysis_H
	#define THEPEG_FxFxAnalysis_H
	//
	// This is the declaration of the FxFxAnalysis class.
	//

	#include "ThePEG/Handlers/AnalysisHandler.h"
	#include "Rivet/AnalysisHandler.hh"
	#include "FxFxReader.h"
	#include "FxFxEventHandler.h"

	namespace ThePEG {

	/**
	* Here is the documentation of the FxFxAnalysis class.
	*
	* @see \ref FxFxAnalysisInterfaces "The interfaces"
	* defined for FxFxAnalysis.
	*/
	class FxFxAnalysis: public ThePEG::AnalysisHandler {

	public:

	/**
	* The default constructor.
	*/
	FxFxAnalysis();

	public:

	/** @name Virtual functions required by the AnalysisHandler class. */
	//@{
	/**
	* Analyze a given Event. Note that a fully generated event
	* may be presented several times, if it has been manipulated in
	* between. The default version of this function will call transform
	* to make a lorentz transformation of the whole event, then extract
	* all final state particles and call analyze(tPVector) of this
	* analysis object and those of all associated analysis objects. The
	* default version will not, however, do anything on events which
	* have not been fully generated, or have been manipulated in any
	* way.
	* @param event pointer to the Event to be analyzed.
	* @param ieve the event number.
	* @param loop the number of times this event has been presented.
	* If negative the event is now fully generated.
	* @param state a number different from zero if the event has been
	* manipulated in some way since it was last presented.
	*/
	virtual void analyze(ThePEG::tEventPtr event, long ieve, int loop, int state);

	/**
	* Produca a HepMC event for the given subprocess
	*/
	HepMC::GenEvent * makeEvent(tEventPtr event, tSubProPtr sub, long no,
	Energy eUnit, Length lUnit,
	CrossSection xsec, CrossSection xsecErr) const;

	/**
	* Produca a HepMC event for the given subprocess
	*/
	- HepMC::GenEvent * makeEventW(tEventPtr event, tSubProPtr sub, long no,
	+ HepMC::GenEvent * makeEventW(tEventPtr event, long no,
	Energy eUnit, Length lUnit,
	CrossSection xsec, CrossSection xsecErr, double evoptweight, double centralweight) const;

	int _i;

	int _numweights;

	//@}

	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(ThePEG::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(ThePEG::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();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual ThePEG::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 ThePEG::IBPtr fullclone() const;
	//@}

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object. Called in the read phase.
	*/
	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();
	//@}

	private:

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

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

	private:
	/**
	* The PDG ID to be used for remnants
	*/
	long _remnantId;

	/**
	* The HepMC format
	*/
	int _format;

	/**
	* Selector for the choice of units
	*/
	int _unitchoice;

	/**
	* Choice of output precision in GenEvent format
	*/
	unsigned int _geneventPrecision;

	/**
	* The Analyses to use
	*/
	vector<string> _analyses;

	/**
	* The base name of the output file.
	*/
	string filename;

	/**
	* Enable debugging information from FxFx
	*/
	bool debug;

	/**
	* Enable use of optional weights in analysis
	*/
	bool useoptweights = false;

	/**
	* normalize optional weights to the central weight
	*/

	bool normoptweights = false;

	/**
	* The FxFxAnalysisHandler
	*/
	Rivet::AnalysisHandler * _rivet;

	/**
	* The FxFxAnalysisHandlers for multiple weights
	*/
	Rivet::AnalysisHandler * _rivetMULTI[120];


	/**
	* holders of weights and cross section
	*/

	std::vector< std::pair<int,double> > OptWeights;
	std::vector<double> OptXS;
	map<string,CrossSection> optxsec;


	/**
	* Event count
	*/
	unsigned long _nevent;
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FxFxAnalysis class and the shared object where it is defined. */
	template <>
	struct ClassTraits<FxFxAnalysis>
	: public ClassTraitsBase<FxFxAnalysis> {
	/** Return a platform-independent class name */
	static string className() { return "ThePEG::FxFxAnalysis"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FxFxAnalysis is implemented. It may also include several, space-separated,
	* libraries if the class FxFxAnalysis depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "FxFxAnalysis.so"; }
	};

	/** @endcond */

	}

	#endif /* THEPEG_FxFxAnalysis_H */
	diff --git a/MatrixElement/FxFx/FxFxEventHandler.cc b/MatrixElement/FxFx/FxFxEventHandler.cc
	--- a/MatrixElement/FxFx/FxFxEventHandler.cc
	+++ b/MatrixElement/FxFx/FxFxEventHandler.cc
	@@ -1,652 +1,652 @@
	// -- C++ --
	//
	// Based on:
	// FxFxEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the FxFxEventHandler class.
	//

	#include "FxFxEventHandler.h"
	#include "FxFxReader.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;

	FxFxEventHandler::~FxFxEventHandler() {}

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

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

	void FxFxEventHandler::doinit() {
	EventHandler::doinit();

	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	readers()[i]->init();
	//FxFxReader & reader = *readers()[i];
	//reader.initialize(*this);
	//weightnames = reader.optWeightsNamesFunc();
	}
	/*
	XSecStat* initxsecs = new XSecStat[weightnames.size()];

	for(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;

	}

	void FxFxEventHandler::initialize() {

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

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

	// Go through all the readers and collect information about cross
	// sections and processes.
	typedef map<int,tFxFxReaderPtr> ProcessMap;
	ProcessMap processes;
	PDPair incoming;
	Energy MaxEA = ZERO;
	Energy MaxEB = ZERO;
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	FxFxReader & reader = *readers()[i];
	reader.initialize(*this);
	weightnames = reader.optWeightsNamesFunc();
	// Check that the incoming particles are consistent between the
	// readers.
	if ( !incoming.first ) {
	incoming.first = getParticleData(reader.heprup.IDBMUP.first);
	if ( !incoming.first )
	Throw<FxFxInitError>()
	<< "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<FxFxInitError>()
	<< "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 )
	Repository::clog()
	<< "The different FxFxReader objects in the "
	<< "FxFxEventHandler '" << 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 FxFxInitError()
	<< "The reader '" << reader.name()
	<< "' contains negatively weighted events, "
	<< "which is not allowed for the FxFxEventHandler '"
	<< 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<FxFxPNumException>()
	<< "In the FxFxEventHandler '"
	<< 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);
	}
	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());
	}
	// Check that we have any cross section at all.
	if ( stats.maxXSec() <= ZERO )
	throw FxFxInitError()
	<< "The sum of the cross sections of the readers in the "
	<< "FxFxEventHandler '" << 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 FxFxEventHandler::doinitrun() {
	EventHandler::doinitrun();
	stats.reset();
	histStats.reset();

	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	readers()[i]->init();
	FxFxReader & reader = *readers()[i];
	reader.initialize(*this);
	weightnames = reader.optWeightsNamesFunc();
	}

	// XSecStat initxsecs;
	XSecStat* initxsecs = new XSecStat[weightnames.size()];
	- for(int ww = 0; ww < weightnames.size(); ww++){
	+ for(size_t 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));
	optstats.insert(std::make_pair(weightnames[ww], initxsecs[ww]));
	opthistStats.insert(std::make_pair(weightnames[ww], initxsecs[ww]));
	optxs.insert(std::make_pair(weightnames[ww], initxs));
	}
	ntries = 0;

	}

	EventPtr FxFxEventHandler::generateEvent() {

	LoopGuard<EventLoopException,FxFxEventHandler>
	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);
	}

	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);
	}
	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;


	// fact for weight normalization
	double fact = theNormWeight ? double(selector().sum()/picobarn) : 1.;

	try {

	theLastXComb = currentReader()->getXComb();

	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
	generator()->currentEventNumber(), weight*fact)));
	currentEvent()->optionalWeights() = currentReader()->optionalEventWeights();
	// normalize the optional weights
	for(map<string,double>::iterator it = currentEvent()->optionalWeights().begin();
	it!=currentEvent()->optionalWeights().end();++it) {
	if(it->first!="ecom"&& it->second!=-999 && it->second!=-111 && it->second!=-222 && it->second!=-333) { it->second *= fact; }
	}


	//print optional weights here
	// cout << "event weight = " << weight << " fact = " << fact << endl;
	/*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
	// 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 FxFxEventHandler::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 FxFxEventHandler::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);
	}
	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 FxFxEventHandler::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());

	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 FxFxEventHandler::continueEvent() {
	try {
	continueCollision();
	}
	catch (Veto) {
	const double fact =
	theNormWeight ?
	double(selector().sum()/picobarn) : 1.;
	reject(currentEvent()->weight()/fact);
	}
	catch (Stop) {
	}
	catch (Exception &) {
	const double fact =
	theNormWeight ?
	double(selector().sum()/picobarn) : 1.;
	reject(currentEvent()->weight()/fact);
	throw;
	}
	return currentEvent();
	}

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

	void FxFxEventHandler::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;
	}

	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 ) {
	FxFxReader & 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;
	}
	os << line;
	} else {

	os << "Per Les Houches Reader (and process #) breakdown:\n";

	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	FxFxReader & 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;
	CrossSection xsectot = reader.stats.xSec();
	if ( xsectot != ZERO ) xsectot /= reader.stats.sumWeights();
	typedef FxFxReader::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;
	}
	os << line;
	}
	}

	string warn = "Warning: Result may be statistically incorrect since\n"
	" the following FxFxReaders were oversampled:\n";
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	FxFxReader & 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;
	}
	}
	}


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

	void FxFxEventHandler::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();
	}
	}

	void FxFxEventHandler::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++;
	}
	}



	map<string,CrossSection> FxFxEventHandler::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;
	}

	return result;
	}


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

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

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

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

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

	}

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

	ClassDescription<FxFxEventHandler>
	FxFxEventHandler::initFxFxEventHandler;
	// Definition of the static class description member.

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

	void FxFxEventHandler::Init() {

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


	static RefVector<FxFxEventHandler,FxFxReader>
	interfaceFxFxReaders
	("FxFxReaders",
	"Objects capable of reading events from an event file or an "
	"external matrix element generator.",
	&FxFxEventHandler::theReaders, -1, false, false, true, false, false);

	static Switch<FxFxEventHandler,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.",
	&FxFxEventHandler::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<FxFxEventHandler,bool> interfaceWarnPNum
	("WarnPNum",
	"Warn if the same process number is used in more than one "
	"FxFxReader.",
	&FxFxEventHandler::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<FxFxEventHandler,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.",
	&FxFxEventHandler::theUnitTolerance, 1.0e-6, 0.0, 0,
	true, false, Interface::lowerlim,
	&FxFxEventHandler::setUnitTolerance,
	(double(FxFxEventHandler::*)()const)(0),
	(double(FxFxEventHandler::*)()const)(0),
	(double(FxFxEventHandler::*)()const)(0),
	(double(FxFxEventHandler::*)()const)(0));

	static Switch<FxFxEventHandler,unsigned int> interfaceWeightNormalization
	("WeightNormalization",
	"How to normalize the output weights",
	&FxFxEventHandler::theNormWeight, 0, false, false);
	static SwitchOption interfaceWeightNormalizationUnit
	(interfaceWeightNormalization,
	"Normalized",
	"Standard normalization, i.e. +/- for unweighted events",
	0);
	static SwitchOption interfaceWeightNormalizationCrossSection
	(interfaceWeightNormalization,
	"CrossSection",
	"Normalize the weights to the max cross section in pb",
	1);


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

	}

	diff --git a/MatrixElement/FxFx/FxFxHandler.cc b/MatrixElement/FxFx/FxFxHandler.cc
	--- a/MatrixElement/FxFx/FxFxHandler.cc
	+++ b/MatrixElement/FxFx/FxFxHandler.cc
	@@ -1,1656 +1,1655 @@
	#include "FxFxHandler.h"
	#include "FxFxReader.h"
	#include "FxFxReader.fh"
	#include "FxFxEventHandler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
	#include "Herwig/PDF/HwRemDecayer.h"
	#include <queue>
	#include "ThePEG/Utilities/Throw.h"
	#include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h"
	#include "fastjet/PseudoJet.hh"
	#include "fastjet/ClusterSequence.hh"
	#include "gsl/gsl_rng.h"
	#include "gsl/gsl_randist.h"
	#include <boost/algorithm/string.hpp>


	using namespace Herwig;
	using namespace ThePEG;

	bool recordEntry(PPtr i,PPtr j) {
	return (i->number()<j->number());
	}
	bool pTsortFunction(PPtr i,PPtr j) {
	return (i->momentum().perp2()>j->momentum().perp2());
	}
	bool ETsortFunction(pair<Energy, Lorentz5Momentum> i,pair<Energy, Lorentz5Momentum> j) {
	return (i.first>j.first);
	}
	bool isMomLessThanEpsilon(Lorentz5Momentum p,Energy epsilon) {
	return (abs(p.x())<epsilon&&abs(p.y())<epsilon&&
	abs(p.z())<epsilon&&abs(p.t())<epsilon);
	}

	FxFxHandler::FxFxHandler()
	: ncy_(100),ncphi_(60),ihvy_(-999),nph_(-999),nh_(-999),
	- etclusmean_(20*GeV),rclus_(0.4),etaclmax_(5.0),rclusfactor_(1.5),
	- ihrd_(-999),njets_(-999),drjmin_(-999), highestMultiplicity_(false),
	- ycmax_(5.4),ycmin_(-5.4),jetAlgorithm_(1),vetoIsTurnedOff_(false),vetoSoftThanMatched_(false), etclusfixed_(true),epsetclus_(2.5*GeV), mergemode_(0), vetoHeavyQ_(true), hpdetect_(true), vetoHeavyFlavour_(false)
	+ mergemode_(0), vetoHeavyQ_(true), vetoHeavyFlavour_(false), etclusmean_(20GeV), jetAlgorithm_(1),vetoIsTurnedOff_(false),vetoSoftThanMatched_(false),etclusfixed_(true),epsetclus_(2.5GeV),rclus_(0.4),etaclmax_(5.0),rclusfactor_(1.5), hpdetect_(true), ihrd_(-999),njets_(-999),drjmin_(-999), highestMultiplicity_(false), ycmax_(5.4),ycmin_(-5.4)
	{}

	void FxFxHandler::doinitrun() {
	QTildeShowerHandler::doinitrun();
	// et_ holds the ET deposited in the (ncy_ x ncphi_) calorimeter cells.
	et_.resize(ncy_);
	for(unsigned int ixx=0; ixx<et_.size(); ixx++) et_[ixx].resize(ncphi_);
	// jetIdx_ for a given calorimeter cell this holds the index of the jet
	// that the cell was clustered into.
	jetIdx_.resize(ncy_);
	for(unsigned int ixx=0; ixx<jetIdx_.size(); ixx++) jetIdx_[ixx].resize(ncphi_);
	if(mergemode_ == 0) {
	cout << "Merging mode is FxFx." << endl;
	} else if (mergemode_ == 1) {
	cout << "Merging mode is Tree level." << endl;
	} else if (mergemode_ == 2) {
	cout << "Merging mode is Tree level, using MadGraph Les Houches information." << endl;
	}
	if(hpdetect_ && mergemode_ != 2) {
	cout << "Automatic detection of hard process enabled." << endl;
	ihrd_ = -999;
	}
	}

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

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

	void FxFxHandler::persistentOutput(PersistentOStream & os) const {
	os << alphaS_
	<< ncy_ << ncphi_ << ihvy_ << nph_ << nh_
	<< ounit(etclusmean_,GeV) << rclus_ << etaclmax_ << rclusfactor_
	<< ihrd_ << njets_ << drjmin_ << highestMultiplicity_
	<< ycmax_ << ycmin_ << jetAlgorithm_ << vetoIsTurnedOff_ << vetoSoftThanMatched_ << etclusfixed_
	<< cphcal_ << sphcal_ << cthcal_ << sthcal_ << ounit(epsetclus_,GeV) << vetoHeavyQ_ << mergemode_ << hpdetect_ << vetoHeavyFlavour_;
	}

	void FxFxHandler::persistentInput(PersistentIStream & is, int) {
	is >> alphaS_
	>> ncy_ >> ncphi_ >> ihvy_ >> nph_ >> nh_
	>> iunit(etclusmean_,GeV) >> rclus_ >> etaclmax_ >> rclusfactor_
	>> ihrd_ >> njets_ >> drjmin_ >> highestMultiplicity_
	>> ycmax_ >> ycmin_ >> jetAlgorithm_ >> vetoIsTurnedOff_ >> vetoSoftThanMatched_ >> etclusfixed_
	>> cphcal_ >> sphcal_ >> cthcal_ >> sthcal_ >> iunit(epsetclus_,GeV) >> vetoHeavyQ_ >> mergemode_ >> hpdetect_ >> vetoHeavyFlavour_;
	}

	ClassDescription<FxFxHandler> FxFxHandler::initFxFxHandler;
	// Definition of the static class description member.

	void FxFxHandler::Init() {

	static ClassDocumentation<FxFxHandler> documentation
	("The FxFxHandler class performs MEPS merging "
	"using the MLM procedure.");

	static Reference<FxFxHandler,ShowerAlpha> interfaceShowerAlpha
	("ShowerAlpha",
	"The object calculating the strong coupling constant",
	&FxFxHandler::alphaS_, false, false, true, false, false);

	static Parameter<FxFxHandler,int> interfaceihvy
	("ihvy",
	"heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t)",
	&FxFxHandler::ihvy_, -999, -999, 7,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,int> interfacenph
	("nph",
	"Number of photons in the AlpGen process",
	&FxFxHandler::nph_, -999, -999, 7,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,int> interfacenh
	("nh",
	"Number of higgses in the AlpGen process",
	&FxFxHandler::nph_, -999, -999, 7,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,Energy> interfaceETClus
	("ETClus",
	"The ET threshold defining a jet in the merging procedure",
	&FxFxHandler::etclusmean_, GeV, 20GeV, 0GeV, 14000*GeV,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,double> interfaceRClus
	("RClus",
	"The cone size used to define a jet in the merging procedure",
	&FxFxHandler::rclus_, 0.4, 0.0, 4.0,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,double> interfaceEtaClusMax
	("EtaClusMax",
	"The maximum \|eta\| used to define a jet in the merging procedure",
	&FxFxHandler::etaclmax_, 5.0, 0.0, 15.0,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,double> interfaceRClusFactor
	("RClusFactor",
	"The prefactor for RClus used to define the jet-parton matching "
	"distance",
	&FxFxHandler::rclusfactor_, 1.5, 0.0, 4.0,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,int> interfaceihrd
	("ihrd",
	"The hard process code",
	&FxFxHandler::ihrd_, -999, 0, 10000,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,int> interfacenjetsmax
	("njetsmax",
	"The number of light jets in the maximum-multiplicity process",
	&FxFxHandler::njets_, -999, 0, 10000,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,double> interfacedrjmin
	("drjmin",
	"Mimimum parton-parton R-sep used for generation.",
	&FxFxHandler::drjmin_, 0.7, 0.0, 4.0,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,bool> interfacehighestMultiplicity
	("highestMultiplicity",
	"If true it indicates that this is the highest multiplicity input "
	"ME-level configuration to be processed.",
	&FxFxHandler::highestMultiplicity_, 0, 0, 1,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,bool> interfaceETClusFixed
	("ETClusFixed",
	"If false, indicates that the jet merging scale, etclus_ is allowed to vary"
	"according to epsetclus_",
	&FxFxHandler::etclusfixed_, 1, 0, 1,
	false, false, Interface::limited);

	static Parameter<FxFxHandler,Energy> interfaceEpsilonETClus
	("EpsilonETClus",
	"The ET threshold defining a jet in the merging procedure",
	&FxFxHandler::epsetclus_, GeV, 2.5GeV, 0GeV, 100.0*GeV,
	false, false, Interface::limited);

	static Switch<FxFxHandler,int> interfaceJetAlgorithm
	("JetAlgorithm",
	"Determines the jet algorithm for finding jets in parton-jet "
	"matching in the MLM procedure.",
	&FxFxHandler::jetAlgorithm_, 1, false, false);
	static SwitchOption AntiKt
	(interfaceJetAlgorithm,
	"AntiKt",
	"The anti-kt jet algorithm.",
	-1);
	static SwitchOption CambridgeAachen
	(interfaceJetAlgorithm,
	"CambridgeAachen",
	"The Cambridge-Aachen jet algorithm.",
	0);
	static SwitchOption Kt
	(interfaceJetAlgorithm,
	"Kt",
	"The Kt jet algorithm.",
	1);

	static Switch<FxFxHandler,int> interfaceMergeMode
	("MergeMode",
	"The choice of merging mode",
	&FxFxHandler::mergemode_, 0, false, false);
	static SwitchOption FxFx
	(interfaceMergeMode,
	"FxFx",
	"FxFx merging.",
	0);
	static SwitchOption Tree
	(interfaceMergeMode,
	"Tree",
	"Tree-level merging.",
	1);
	static SwitchOption TreeMG
	(interfaceMergeMode,
	"TreeMG5",
	"Tree-level merging using the MadGraph pt clustering information.",
	2);


	static Switch<FxFxHandler,bool> interfaceHardProcessDetection
	("HardProcessDetection",
	"The choice of merging mode",
	&FxFxHandler::hpdetect_, true, false, false);
	static SwitchOption Automatic
	(interfaceHardProcessDetection,
	"Automatic",
	"Automatically determine which particles to include in the merging.",
	true);
	static SwitchOption Manual
	(interfaceHardProcessDetection,
	"Manual",
	"Use the ihrd code to determine which particles to include in the merging.",
	false);

	static Switch<FxFxHandler,bool> interfaceVetoIsTurnedOff
	("VetoIsTurnedOff",
	"Allows the vetoing mechanism to be switched off.",
	&FxFxHandler::vetoIsTurnedOff_, false, false, false);
	static SwitchOption VetoingIsOn
	(interfaceVetoIsTurnedOff,
	"VetoingIsOn",
	"The MLM merging veto mechanism is switched ON.",
	false);
	static SwitchOption VetoingIsOff
	(interfaceVetoIsTurnedOff,
	"VetoingIsOff",
	"The MLM merging veto mechanism is switched OFF.",
	true);


	static Switch<FxFxHandler,bool> interfaceVetoHeavyFlavour
	("VetoHeavyFlavour",
	"Allows the heavy flavour vetoing mechanism to be switched off.",
	&FxFxHandler::vetoHeavyFlavour_, false, false, false);
	static SwitchOption HeavyFVetoingIsOn
	(interfaceVetoHeavyFlavour,
	"Yes",
	"The MLM merging veto mechanism for heavy flavour is switched ON.",
	true);
	static SwitchOption HeavyFVetoingIsOff
	(interfaceVetoHeavyFlavour,
	"No",
	"The MLM merging veto mechanism for heavy flavour is switched OFF.",
	false);

	static Switch<FxFxHandler,bool> interfaceHeavyQVeto
	("HeavyQVeto",
	"Allows the vetoing mechanism on the heavy quark products to be switched off.",
	&FxFxHandler::vetoHeavyQ_, false, false, false);
	static SwitchOption HQVetoingIsOn
	(interfaceHeavyQVeto,
	"Yes",
	"The MLM merging veto on Heavy quark decay produts mechanism is switched ON.",
	true);
	static SwitchOption HQVetoingIsOff
	(interfaceHeavyQVeto,
	"No",
	"The MLM merging veto on Heavy quark decay products mechanism is switched OFF.",
	false);

	static Switch<FxFxHandler,bool> interfaceVetoSoftThanMatched
	("VetoSoftThanMatched",
	"Allows the vetoing mechanism to be switched off.",
	&FxFxHandler::vetoSoftThanMatched_, false, false, false);
	static SwitchOption VetoSoftIsOn
	(interfaceVetoSoftThanMatched,
	"VetoSoftIsOn",
	"The vetoing of highest-mult. events with jets softer than matched ones is ON",
	true);
	static SwitchOption VetoSoftIsOff
	(interfaceVetoSoftThanMatched,
	"VetoSoftIsOff",
	"The vetoing of highest-mult. events with jets softer than matched ones is OFF.",
	false);


	}

	void FxFxHandler::dofinish() {
	QTildeShowerHandler::dofinish();
	}

	void FxFxHandler::doinit() {

	//print error if HardProcID is not set in input file
	if(ihrd_ == -999 && !hpdetect_) { cout << "Error: FxFxHandler:ihrd not set and FxFx:HardProcessDetection set to Manual!" << endl; exit(1); }
	QTildeShowerHandler::doinit();

	// Compute calorimeter edges in rapidity for GetJet algorithm.
	ycmax_=etaclmax_+rclus_;
	ycmin_=-ycmax_;
	}

	// Throws a veto according to MLM strategy ... when we finish writing it.
	bool FxFxHandler::showerHardProcessVeto() const {
	int debug_mode = 0;
	if(vetoIsTurnedOff_) {
	// cout << "Vetoing is turned OFF." << endl;
	return false;
	}

	//if(debug_mode) { cout << "debug_mode = " << 5 << endl; }

	// Skip veto for processes in which merging is not implemented:
	if(ihrd_==7\|\|ihrd_==8\|\|ihrd_==13) {
	ostringstream wstring;
	wstring << "FxFxHandler::showerHardProcessVeto() - warning."
	<< "MLM merging not implemented "
	<< "processes 4Q (ihrd=7), QQh (ihrd=8), "
	<< "(single) top (ihrd=13) \n";
	generator()->logWarning( Exception(wstring.str(),
	Exception::warning) );
	return false;
	}

	// Fill preshowerISPs_ pair and preshowerFSPs_ particle pointer vector.
	getPreshowerParticles();

	// Fill showeredISHs_, showeredISPs and showeredRems pairs, as well as
	// showeredFSPs_ particle pointer vector.
	getShoweredParticles();

	// Turn on some screen output debugging: 0 = none ---> 5 = very verbose.
	doSanityChecks(debug_mode);


	// Dimensions of each calorimter cell in y and phi.
	dely_ = (ycmax_-ycmin_)/double(ncy_);
	delphi_ = 2*M_PI/double(ncphi_);

	// Fill partonsToMatch_ with only those pre-shower partons intended to
	// used in jet-parton matching and fill particlesToCluster_ using only
	// those final state particles (post-shower) which are supposed to go
	// in the jet clustering used to do merging.
	partonsToMatch_ = preshowerFSPs_;
	particlesToCluster_ = showeredFSPs_ ;

	// Filter out all but the 'extra' light-parton progenitors and their
	// associated final state particles.
	if(mergemode_ == 0 \|\| mergemode_ == 1) { caldel_m(); }
	else if(mergemode_ == 2) { caldel_mg(); }

	double prob(1);
	//if etclusfixed_ then set the etclus_ to the fixed chosen value
	if(etclusfixed_) {
	etclus_ = etclusmean_;
	} else {
	//else, if we wish to vary etclus_, we use the probability distribution
	//choose a probability between 0 and 1
	prob = rnd();
	etclus_ = etclusran_(prob);
	}

	// Cluster particlesToCluster_ into jets with FastJet.
	getFastJets(rclus_,etclus_,etaclmax_);

	if(mergemode_ == 0) {

	// Get npLO_ and npNLO_ for FxFx matching
	getnpFxFx();

	// print the npXLO_ values obtained
	// cout << "HANDLER:\t\t\t\t" << npLO_ << "\t\t" << npNLO_ << endl;

	//FxFx modifications start here.

	// Sort partonsToMatch_ from high to low pT.
	sort(partonsToMatch_.begin(),partonsToMatch_.end(),pTsortFunction);

	// Count the number of jets.
	int njets_found(pjet_.size());

	// If the number of jets found is not equal to the number of partons in the Born
	// (i.e., the number of partons in the S-event, or one less than the number of partons in an H-event),
	// the jets cannot be matched and the event has to be rejected. The number of partons in the Born is written in the event file with a name “npNLO”
	// if there are no jets to match and no jets have been found, do not veto.
	if(njets_found == 0 && npNLO_ == 0) { /cout << "njets_found = " << njets_found << " and npNLO = " << npNLO_ << ", accepting" << endl;/ return false; }

	//if the number of jets is smaller than npNLO -> reject the event.
	if(njets_found < npNLO_) { /cout << "njets_found = " << njets_found << " and npNLO = " << npNLO_ << ", rejecting" << endl;/ return true; }
	// For the maximum-multiplicity sample, the number of jets obtained does not have to be exactly equal to npNLO, it may also be larger;
	if(njets_found > npNLO_ && npNLO_ != njets_) { /cout << "njets_found = " << njets_found << " and npNLO = " << npNLO_ << ", rejecting" << endl;/ return true; }

	// Create the matrix-element jets.
	// Cluster also the partons at the hard-matrix element level into jets with the same algorithm as above,
	// but without the requirement of a minimal pT on the jets (or set it very small).
	// By construction, for S-events you should find exactly npNLO jets, while for the H-events it is either npNLO or npNLO+1.
	// Cluster partonsToMatch_ into jets with FastJet.
	- getFastJetsToMatch(rclus_,0*GeV,etaclmax_);
	+ getFastJetsToMatch(rclus_);

	- int me_njets_found(pjetME_.size());
	+ //int me_njets_found(pjetME_.size());

	// cout << "number of ME jets found = " << me_njets_found << "partons to match: " << partonsToMatch_.size() << endl;

	// Match light progenitors to jets.
	vector<int> jetToPartonMap(pjetME_.size(),-999);
	Energy etmin(777e100*GeV);

	// Match the jets.
	// Try to match the “npNLO” hardest jets created post-shower with any of the jets pre-shower. Two jets are matched if the distance between them is smaller than 1.5*DeltaR.
	// If not all the npNLO hardest shower jets are matched the event has to be rejected.
	// Note that if the current event does not belong to the maximum multiplicity sample, this means that all the shower jets need to be matched, because the requirement above already rejects
	// events that do not have npNLO shower jets.
	// For those events, at the level of the matrix elements there can either be npNLO or npNLO+1 matrix-element jets, depending on S- or H-events and the kinematics of those partons.
	// Still only the shower jets need to be matched, so an event should not be rejected if a matrix-element jet cannot be matched.
	// For each parton, starting with the hardest one ...
	- for(unsigned int ixx=0; ixx<npNLO_; ixx++) {
	+ for(int ixx=0; ixx<npNLO_; ixx++) {
	// ... loop over all jets not already matched.
	double DRmin(777e100);
	int jetIndexForDRmin(-999);
	for(unsigned int jxx=0; jxx<pjetME_.size(); jxx++) {
	// ... and tag closest of the remaining ones
	double DRpartonJet(partonJetDeltaR(pjetME_[ixx],pjet_[jxx]));
	if(jetToPartonMap[jxx]<0&&DRpartonJet<DRmin) {
	DRmin=DRpartonJet;
	jetIndexForDRmin=jxx;
	}
	}
	// If the parton-jet distance is less than the matching
	// distance, the parton and jet match.
	if(DRmin<rclus_*rclusfactor_&&jetIndexForDRmin>=0) {
	jetToPartonMap[jetIndexForDRmin]=ixx;
	if(ixx==0\|\|etjet_[jetIndexForDRmin]<etmin)
	etmin=etjet_[jetIndexForDRmin];
	// Otherwise this parton is not matched so veto the event.
	} else return true;
	}

	// Veto events where matched jets are softer than non-matched ones,
	// in the inclusive (highestMultiplicity_ = true) mode, unless we
	// are dealing with NLO input events.
	if(npNLO_ == njets_ && vetoSoftThanMatched_) {
	//cout << "highest mult. event being tested for softer jets than matched..." << endl;
	for(unsigned int iyy=0; iyy<pjet_.size(); iyy++) {
	//cout << "etjet_[iyy] = " << etjet_[iyy]/GeV << endl;
	if(jetToPartonMap[iyy]<0&&etmin<etjet_[iyy]) { /cout << "VETO!" << endl;/ return true; }
	}
	}


	if(!vetoHeavyQ_) {
	//cout << "no heavy quark decay product veto!" << endl;
	return false;
	}

	//end of FxFx part
	// **************************************************************** //
	// * Now look to the non-light partons for heavy quark processes. * //
	// **************************************************************** //

	if( (ihrd_<=2\|\|ihrd_==6\|\|ihrd_==10\|\|ihrd_==15\|\|ihrd_==16) \|\| (hpdetect_==true && hvqfound==true) ) {
	// Extract heavy quark progenitors and the radiation they
	// produce and put it in the calorimeter.
	caldel_hvq();

	// Cluster particlesToCluster_ into jets with FastJet.
	getFastJets(rclus_,etclus_,etaclmax_);

	// If the radiation from the heavy quarks does not give rise
	// to any jets we accept event.
	if(pjet_.size() == 0) return false;

	// If extra jets emerge from the jet clustering we only
	// accept events where the jets formed by radiation from
	// b and c quarks lies within drjmin_ of the heavy quark
	// progenitor.
	int nmjet(pjet_.size());
	for(unsigned int ixx=0; ixx<pjet_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
	if(!(abs(partonsToMatch_[jxx]->id())==4\|\|abs(partonsToMatch_[jxx]->id())==5)) continue;
	if(partonJetDeltaR(partonsToMatch_[jxx],pjet_[ixx])<drjmin_) {
	nmjet--; // Decrease the number of unmatched jets.
	etjet_[ixx]=0*GeV; // Set jet ET to zero to indicate it is 'matched'.
	}
	}
	}

	// If every jet matched to _at_least_one_ progenitor accept the event.
	if(nmjet<=0) return false;

	}
	}

	/*****
	*****
	***** END of mergemode_ == 0 (i.e. FxFx MERGING BLOCK)
	*****/



	if(mergemode_ == 1 \|\| mergemode_ == 2) {

	// determine whether event is of highest multiplicity or not
	- if(partonsToMatch_.size()==njets_) { highestMultiplicity_ = true; }
	+ unsigned long int njets_long = njets_;
	+ if(partonsToMatch_.size()==njets_long) { highestMultiplicity_ = true; }

	// cout << "jet finding gives pjet_.size() = " << pjet_.size() << endl;
	// If there are less jets than partons then parton-jet matching is
	// bound to fail: reject the event already. Also, if the input is
	// an NLO event file it will 99.5% of the time contain a number of
	// light partons in the F.S. equal to that in the real emission
	// process in the NLO calculation, moreover, it has already
	// effectively merged njets_-1 and njets jet events. So in that
	// case we do not reject events on the grounds they have jet
	// multiplicity less than partonsToMatch_.size() but rather less
	// jets than partonsToMatch.size()-1; such events are better
	// described by the lower-by-one-unit (partonsToMatch_.size()-1)
	// of multiplicity NLO event file, or the lower-by-two-units
	// (partonsToMatch_.size()-2) of multiplicity LO event file.

	// If it is not jet production apply rejection criterion as above.
	if(ihrd_!=9) {
	if(pjet_.size() < partonsToMatch_.size()) return true;
	}
	// Otherwise, in the case of jet production allow the lowest
	// contributing multiplicity process (just at NLO), namely,
	// dijet production, to give rise to 1-jet and even 0-jet
	// events, since these can contribute to, for example, the
	// inclusive jet cross section i.e. in this case the rejection
	// is only applied in the case of the next-to-lowest multiplicity
	// processes (>2 parton events at LO and >3 parton events at NLO).
	else {
	// KH - March 5th
	// Removed the following line giving special treatment
	// also to the LO events, to maintain consistency with
	// the fortran algorithm, at least for now. So now jet
	// production at LO is being treated the same as all
	// other processes.
	// if(partonsToMatch_.size()==2 && pjet_.size()<2) return false;
	if(pjet_.size() < partonsToMatch_.size()) return true;
	}

	// Sort partonsToMatch_ from high to low pT.
	sort(partonsToMatch_.begin(),partonsToMatch_.end(),pTsortFunction);

	// Match light progenitors to jets.
	vector<int> jetToPartonMap(pjet_.size(),-999);
	Energy etmin(777e100*GeV);

	// For each parton, starting with the hardest one ...
	for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
	// ... loop over all jets not already matched.
	double DRmin(777e100);
	int jetIndexForDRmin(-999);
	for(unsigned int jxx=0; jxx<pjet_.size(); jxx++) {
	// ... and tag closest of the remaining ones
	double DRpartonJet(partonJetDeltaR(partonsToMatch_[ixx],pjet_[jxx]));
	if(jetToPartonMap[jxx]<0&&DRpartonJet<DRmin) {
	DRmin=DRpartonJet;
	jetIndexForDRmin=jxx;
	}
	}
	// If the parton-jet distance is less than the matching
	// distance, the parton and jet match.
	if(DRmin<rclus_*rclusfactor_&&jetIndexForDRmin>=0) {
	jetToPartonMap[jetIndexForDRmin]=ixx;
	if(ixx==0\|\|etjet_[jetIndexForDRmin]<etmin)
	etmin=etjet_[jetIndexForDRmin];
	// Otherwise this parton is not matched so veto the event.
	} else return true;
	}


	// Veto events with larger jet multiplicity from exclusive sample.
	if(!highestMultiplicity_&&pjet_.size()>partonsToMatch_.size()) return true;

	// Veto events where matched jets are softer than non-matched ones,
	// in the inclusive (highestMultiplicity_ = true) mode, unless we
	// are dealing with NLO input events.
	if(highestMultiplicity_) {
	for(unsigned int ixx=0; ixx<pjet_.size(); ixx++)
	if(jetToPartonMap[ixx]<0&&etmin<etjet_[ixx]) return true;
	}

	if(!vetoHeavyQ_) {
	//cout << "no heavy quark decay product veto!" << endl;
	return false;
	}


	// **************************************************************** //
	// * Now look to the non-light partons for heavy quark processes. * //
	// **************************************************************** //

	if( (ihrd_<=2\|\|ihrd_==6\|\|ihrd_==10\|\|ihrd_==15\|\|ihrd_==16) \|\| (hpdetect_==true && hvqfound==true)) {

	// Extract heavy quark progenitors and the radiation they
	// produce and put it in the calorimeter.
	caldel_hvq();

	// Cluster particlesToCluster_ into jets with FastJet.
	getFastJets(rclus_,etclus_,etaclmax_);

	// If the radiation from the heavy quarks does not give rise
	// to any jets we accept event.
	if(pjet_.size() == 0) return false;

	// If extra jets emerge from the jet clustering we only
	// accept events where the jets formed by radiation from
	// b and c quarks lies within drjmin_ of the heavy quark
	// progenitor.
	int nmjet(pjet_.size());
	for(unsigned int ixx=0; ixx<pjet_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
	if(!(abs(partonsToMatch_[jxx]->id())==4\|\|abs(partonsToMatch_[jxx]->id())==5)) continue;
	if(partonJetDeltaR(partonsToMatch_[jxx],pjet_[ixx])<drjmin_) {
	nmjet--; // Decrease the number of unmatched jets.
	etjet_[ixx]=0*GeV; // Set jet ET to zero to indicate it is 'matched'.
	}
	}
	}

	// If every jet matched to _at_least_one_ progenitor accept the event.
	if(nmjet<=0) return false;
	else {
	// If unmatched jets remain, reject the event if highestMultiplicity_!=1
	if(!highestMultiplicity_) return true;
	else {
	// If unmatched jets remain and highestMultiplicity is true then check
	// that these are softer than all the matched ones (from the light-parton
	// matching round).
	Energy etmax(0.*GeV);
	for(unsigned int ixx=0; ixx<pjet_.size(); ixx++) etmax=max(etjet_[ixx],etmax);
	if(etmax>etmin) return true;
	}
	}
	}

	}


	// Otherwise we accept the event ...
	return false;

	}

	/* Function that returns the R distance
	between a particle and a jet. */
	double FxFxHandler::partonJetDeltaR(ThePEG::tPPtr partonptr, LorentzMomentum jetmom) const {
	LorentzMomentum partonmom(partonptr->momentum());
	// Calculate DY, DPhi and then DR
	double DY(partonmom.eta()-jetmom.eta());
	double DPhi(partonmom.phi()-jetmom.phi());
	if(DPhi>M_PI) DPhi=2*M_PI-DPhi;
	double DR(sqrt(sqr(DY)+sqr(DPhi)));
	return DR;
	}

	double FxFxHandler::partonJetDeltaR(LorentzMomentum jetmom1, LorentzMomentum jetmom2) const {
	// Calculate DY, DPhi and then DR
	double DY(jetmom1.eta()-jetmom2.eta());
	double DPhi(jetmom1.phi()-jetmom2.phi());
	if(DPhi>M_PI) DPhi=2*M_PI-DPhi;
	double DR(sqrt(sqr(DY)+sqr(DPhi)));
	return DR;
	}

	// Get FastJets
	void FxFxHandler::getFastJets(double rjet, Energy ejcut, double etajcut) const {

	vector<fastjet::PseudoJet> particlesToCluster;
	for(unsigned int ipar=0; ipar<particlesToCluster_.size(); ipar++) {
	double y(particlesToCluster_[ipar]->momentum().eta());
	if(y>=ycmin_&&y<=ycmax_) {
	int absId(abs(particlesToCluster_[ipar]->id()));
	// If it's not a lepton / top / photon it may go in the jet finder.
	if(!(absId>=11&&absId<=16) && absId!=6 && absId!=22) {
	// input particles into fastjet pseudojet
	fastjet::PseudoJet p(particlesToCluster_[ipar]->momentum().x()/GeV,
	particlesToCluster_[ipar]->momentum().y()/GeV,
	particlesToCluster_[ipar]->momentum().z()/GeV,
	particlesToCluster_[ipar]->momentum().e()/GeV);
	p.set_user_index(ipar);
	particlesToCluster.push_back(p);
	}
	}
	}

	fastjet::RecombinationScheme recombinationScheme = fastjet::E_scheme;
	fastjet::Strategy strategy = fastjet::Best;
	double R(rjet);
	fastjet::JetDefinition theJetDefinition;
	switch (jetAlgorithm_) {
	case -1: theJetDefinition=fastjet::JetDefinition(fastjet::antikt_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	case 0: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	case 1: theJetDefinition=fastjet::JetDefinition(fastjet::kt_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	default: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	}
	fastjet::ClusterSequence fastjetEvent(particlesToCluster,theJetDefinition);
	vector<fastjet::PseudoJet> inclusiveJets = fastjetEvent.inclusive_jets();
	inclusiveJets = fastjet::sorted_by_pt(inclusiveJets);

	// Fill the array of jet momenta for the rest of the veto procedure.
	pjet_.clear();
	pjet_.resize(inclusiveJets.size());
	etjet_.clear();
	etjet_.resize(inclusiveJets.size());
	for(unsigned int ffj=0; ffj<pjet_.size();ffj++) {
	pjet_[ffj] = Lorentz5Momentum(inclusiveJets[ffj].px()*GeV,
	inclusiveJets[ffj].py()*GeV,
	inclusiveJets[ffj].pz()*GeV,
	inclusiveJets[ffj].e()*GeV);
	pjet_[ffj].rescaleMass();
	etjet_[ffj] = pjet_[ffj].et();
	}

	// Throw the jet away if it's outside required eta region or
	// has transverse energy below ejcut.
	for(unsigned int fj=0; fj<pjet_.size(); fj++)
	if(etjet_[fj]<ejcut\|\|fabs(pjet_[fj].eta())>etajcut) {
	pjet_.erase(pjet_.begin()+fj);
	etjet_.erase(etjet_.begin()+fj);
	fj--;
	}

	// Sort jets from high to low ET.
	vector<pair<Energy, Lorentz5Momentum> > etjet_pjet;
	for(unsigned int ixx=0; ixx<etjet_.size(); ixx++)
	etjet_pjet.push_back(make_pair(etjet_[ixx],pjet_[ixx]));
	sort(etjet_pjet.begin(),etjet_pjet.end(),ETsortFunction);
	for(unsigned int ixx=0; ixx<etjet_.size(); ixx++) {
	etjet_[ixx]=etjet_pjet[ixx].first;
	pjet_[ixx]=etjet_pjet[ixx].second;
	}

	return;
	}

	// Get FastJets from partonsToMatch_
	-void FxFxHandler::getFastJetsToMatch(double rjet, Energy ejcut, double etajcut) const {
	+void FxFxHandler::getFastJetsToMatch(double rjet) const {

	vector<fastjet::PseudoJet> particlesToCluster;
	for(unsigned int ipar=0; ipar<partonsToMatch_.size(); ipar++) {
	double y(partonsToMatch_[ipar]->momentum().eta());
	if(y>=ycmin_&&y<=ycmax_) {
	int absId(abs(partonsToMatch_[ipar]->id()));
	// If it's not a lepton / top / photon it may go in the jet finder.
	if(!(absId>=11&&absId<=16) && absId!=6 && absId!=22) {
	// input particles into fastjet pseudojet
	fastjet::PseudoJet p(partonsToMatch_[ipar]->momentum().x()/GeV,
	partonsToMatch_[ipar]->momentum().y()/GeV,
	partonsToMatch_[ipar]->momentum().z()/GeV,
	partonsToMatch_[ipar]->momentum().e()/GeV);
	p.set_user_index(ipar);
	particlesToCluster.push_back(p);
	}
	}
	}

	fastjet::RecombinationScheme recombinationScheme = fastjet::E_scheme;
	fastjet::Strategy strategy = fastjet::Best;
	double R(rjet);
	fastjet::JetDefinition theJetDefinition;
	switch (jetAlgorithm_) {
	case -1: theJetDefinition=fastjet::JetDefinition(fastjet::antikt_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	case 0: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	case 1: theJetDefinition=fastjet::JetDefinition(fastjet::kt_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	default: theJetDefinition=fastjet::JetDefinition(fastjet::cambridge_algorithm,
	R,
	recombinationScheme,
	strategy); break;
	}
	fastjet::ClusterSequence fastjetEvent(particlesToCluster,theJetDefinition);
	vector<fastjet::PseudoJet> inclusiveJets = fastjetEvent.inclusive_jets();
	inclusiveJets = fastjet::sorted_by_pt(inclusiveJets);

	// Fill the array of jet momenta for the rest of the veto procedure.
	pjetME_.clear();
	pjetME_.resize(inclusiveJets.size());
	for(unsigned int ffj=0; ffj<pjetME_.size();ffj++) {
	pjetME_[ffj] = Lorentz5Momentum(inclusiveJets[ffj].px()*GeV,
	inclusiveJets[ffj].py()*GeV,
	inclusiveJets[ffj].pz()*GeV,
	inclusiveJets[ffj].e()*GeV);
	pjetME_[ffj].rescaleMass();
	}

	return;
	}

	// Deletes particles from partonsToMatch_ and particlesToCluster_
	// vectors so that these contain only the partons to match to the
	// jets and the particles used to build jets respectively. By and
	// large the candidates for deletion are: vector bosons and their
	// decay products, Higgs bosons, photons as well as _primary_, i.e.
	// present in the lowest multiplicity process, heavy quarks and
	// any related decay products.
	void FxFxHandler::getDescendents(PPtr theParticle) const {
	ParticleVector theChildren(theParticle->children());
	for (unsigned int ixx=0; ixx<theChildren.size(); ixx++)
	if(theChildren[ixx]->children().size()==0)
	tmpList_.push_back(theChildren[ixx]);
	else
	getDescendents(theChildren[ixx]);
	return;
	}

	void FxFxHandler::caldel_m() const {

	preshowerFSPsToDelete_.clear();
	showeredFSPsToDelete_.clear();

	hvqfound = false;

	if(hpdetect_ && mergemode_!=2) {
	for(unsigned int ixx=0; ixx<preshowerFSPs_.size(); ixx++) {
	tmpList_.clear();
	/* Exclude the top quarks and any children
	they may have produced from the jet parton matching. */


	if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6) {
	hvqfound = true;
	// cout << "preshowerFSPs_[ixx]->id() = " << preshowerFSPs_[ixx]->id() << " " << preshowerFSPs_[ixx]->momentum().perp2() << endl;
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]->parents()[0]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++) {
	// cout << "tmpList_[jxx]->id() = " << tmpList_[jxx]->id() << " " << tmpList_[jxx]->momentum().perp2() << endl;
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	continue;
	}

	/* Exclude the v.bosons and any children they may
	have produced from the jet parton matching. */
	if( (abs(preshowerFSPs_[ixx]->parents()[0]->id())==23\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24\|\|
	abs(preshowerFSPs_[ixx]->id())==22\|\|
	abs(preshowerFSPs_[ixx]->id())==25\|\| (abs(preshowerFSPs_[ixx]->id()) < 17 && abs(preshowerFSPs_[ixx]->id()) > 10)) && (abs(preshowerFSPs_[ixx]->parents()[0]->id())!=6)) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	continue;
	}
	/* Exclude the bottom quarks and any children
	they may have produced from the jet parton matching. */
	if(abs(preshowerFSPs_[ixx]->id())==5&&ixx<2) {
	hvqfound = true;
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	continue;
	}
	}
	}
	if(!hpdetect_) {
	for(unsigned int ixx=0; ixx<preshowerFSPs_.size(); ixx++) {
	tmpList_.clear();
	if(ihrd_<=2) {
	/* wqq... , zqq... */
	/* Exclude the heavy quarks and any children they may
	have produced as well as the v.boson and any children
	it may have produced from the jet parton matching. */
	if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==23\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	if(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=4) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "wqq / zqq process should have 4 particles to omit from"
	<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
	}
	} else if(ihrd_<=4) {
	/* zjet... */
	/* Exclude the v.boson and any children
	it may have produced from the jet parton matching. */
	if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==23\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	/*if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=2) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "zjet process should have 2 particles to omit from"
	<< "jet-parton matching." << Exception::eventerror;
	}*/
	} else if(ihrd_==5) {
	/* vbjet... */
	/* Exclude the v.bosons and any children they may
	have produced from the jet parton matching. */
	if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==23\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24\|\|
	abs(preshowerFSPs_[ixx]->id())==22\|\|
	abs(preshowerFSPs_[ixx]->id())==25) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	} else if(ihrd_==6) {
	/* 2Q... */
	/* Exclude the heavy quarks and any children
	they may have produced from the jet parton matching. */
	if(ihvy_==6) {
	if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	// cout << "preshowerFSPs_[ixx]->id() = " << preshowerFSPs_[ixx]->id() << " " << preshowerFSPs_[ixx]->momentum().perp2() << endl;
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	}
	if(abs(preshowerFSPs_[ixx]->parents()[0]->id())==6) {
	getDescendents(preshowerFSPs_[ixx]->parents()[0]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++) {
	// cout << "tmpList_[jxx]->id() = " << tmpList_[jxx]->id() << " " << tmpList_[jxx]->momentum().perp2() << endl;
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	}
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=6) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "2Q process should have 6 particles to omit from"
	<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
	}
	} else {
	if(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=2) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "2Q process should have 2 particles to omit from"
	<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
	}
	}
	} else if(ihrd_==7) {
	/* 4Q... */
	/* There are no light jets for this process, so nothing to match. */
	} else if(ihrd_==9) {
	/* Njet... */
	} else if(ihrd_==10) {
	/* wcjet... */
	/* Exclude the charm quark and any children it may
	have produced as well as the v.boson and any children
	it may have produced from the jet parton matching. */
	if((abs(preshowerFSPs_[ixx]->id())==4&&ixx<1)\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=3) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "wcjet process should have 3 particles to omit from"
	<< "jet-parton matching." << Exception::eventerror;
	}
	} else if(ihrd_==11) {
	/* phjet... */
	/* Exclude the hard photons from the jet parton matching. */
	if(abs(preshowerFSPs_[ixx]->id())==22) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	unsigned int tmpUnsignedInt(nph_);
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "phjet process should have " << nph_ << " particles to omit from"
	<< "jet-parton matching." << Exception::eventerror;
	}
	} else if(ihrd_==12) {
	/* hjet... */
	/* Exclude the higgs and any children it may have
	produced from the jet parton matching. */
	if(abs(preshowerFSPs_[ixx]->id())==25) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=1) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "hjet process should have 1 particle to omit from"
	<< "jet-parton matching." << Exception::eventerror;
	}
	} else if(ihrd_==14) {
	/* wphjet... */
	/* Exclude the v.boson and any children it may have
	produced from the jet parton matching. */
	// AND WHAT ABOUT THE PHOTON? <--- CHECK THIS WITH AlpGen GUYs.
	if(abs(preshowerFSPs_[ixx]->id())==22\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	unsigned int tmpUnsignedInt(2+nph_);
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "wphjet process should have " << 2+nph_ << " particles to omit from"
	<< "jet-parton matching." << Exception::eventerror;
	}
	} else if(ihrd_==15) {
	/* wphqq... <--- N.B. if q = top, it is not decayed. */
	/* Exclude the heavy quarks and any children they may
	have produced as well as the v.boson and any children
	it may have produced from the jet parton matching. */
	// AND WHAT ABOUT THE PHOTON? <--- CHECK THIS WITH AlpGen GUYs.
	if(abs(preshowerFSPs_[ixx]->id())==22\|\|
	(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx==(preshowerFSPs_.size()-(2+nph_+1)))\|\|
	(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx==(preshowerFSPs_.size()-(2+nph_+2)))\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	unsigned int tmpUnsignedInt(4+nph_);
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "wphjet process should have " << 4+nph_ << " particles to omit from"
	<< "jet-parton matching." << Exception::eventerror;
	}
	} else if(ihrd_==16) {
	/* 2Qph... <--- if Q is a top it will decay. */
	/* Exclude the hard photons and any children they
	may have produced from the jet parton matching
	as well as the heavy quarks and any children it
	may have produced. */
	// AND WHAT ABOUT THE PHOTON?
	if(ihvy_==6) {
	if(abs(preshowerFSPs_[ixx]->id())==22\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==6\|\|
	abs(preshowerFSPs_[ixx]->parents()[0]->id())==24) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	unsigned int tmpUnsignedInt(6+nph_);
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "wphjet process should have " << 6+nph_ << " particles to omit from"
	<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
	}
	} else {
	if(abs(preshowerFSPs_[ixx]->id())==22\|\|
	(abs(preshowerFSPs_[ixx]->id())==ihvy_&&ixx<2)) {
	preshowerFSPsToDelete_.push_back(preshowerFSPs_[ixx]);
	getDescendents(preshowerFSPs_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	unsigned int tmpUnsignedInt(2+nph_);
	if(ixx==preshowerFSPs_.size()-1&&preshowerFSPsToDelete_.size()!=tmpUnsignedInt) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "wphjet process should have " << 2+nph_ << " particles to omit from"
	<< "jet-parton matching for ihvy=" << ihvy_ << "." << Exception::eventerror;
	}
	}
	}
	}
	}
	// cout << "partonsToMatch_.size()= " << partonsToMatch_.size() << " preshowerFSPsToDelete_.size() = " << preshowerFSPsToDelete_.size() << endl;
	// cout << "deleting preshowerFSPs" << endl;
	for(unsigned int ixx=0; ixx<preshowerFSPsToDelete_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
	if(preshowerFSPsToDelete_[ixx]==partonsToMatch_[jxx]) {
	// cout << "deleting " << preshowerFSPsToDelete_[ixx]->id() << " with parent " << preshowerFSPsToDelete_[ixx]->parents()[0]->id() << endl; // " energy = " << preshowerFSPsToDelete_[ixx]->momentum().e()*GeV << endl;
	partonsToMatch_.erase(partonsToMatch_.begin()+jxx);
	break;
	}
	}
	}
	//cout << "partonsToMatch_.size() (AFTER) = " << partonsToMatch_.size() << endl;
	//cout << "deleting showeredFSPs" << endl;
	for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
	if(showeredFSPsToDelete_[ixx]==particlesToCluster_[jxx]) {
	// cout << "deleting " << showeredFSPsToDelete_[ixx]->id() << " with parent " << showeredFSPsToDelete_[ixx]->parents()[0]->id() << endl; //" energy = " << preshowerFSPsToDelete_[ixx]->momentum().e()*GeV << endl;
	particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
	break;
	}
	}
	}

	// Sanity check!
	if(partonsToMatch_.size()>particlesToCluster_.size()) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "No. of ME level partons to be matched to jets = "
	<< partonsToMatch_.size() << "\n"
	<< "No. of showered particles to build jets from = "
	<< particlesToCluster_.size() << "\n"
	<< "There should be at least as many partons to\n"
	<< "cluster as there are partons to match to.\n"
	<< Exception::eventerror;
	}
	// cout << "partonsToMatch_.size() (AFTER2) = " << partonsToMatch_.size() << endl;


	// Acid test.
	/* unsigned int tmpUnsignedInt(njets_);
	if(!inputIsNLO_&&partonsToMatch_.size()!=tmpUnsignedInt) {
	for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
	if(abs(partonsToMatch_[ixx]->id())>=6&&
	abs(partonsToMatch_[ixx]->id())!=21)
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "Found a parton to match to which is not a quark or gluon!"
	<< *partonsToMatch_[ixx] << "\n"
	<< Exception::eventerror;
	}
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "No. of ME level partons to be matched to jets = "
	<< partonsToMatch_.size() << "\n"
	<< "No. of light jets (njets) in AlpGen process = "
	<< njets_ << "\n"
	<< "These should be equal." << "\n"
	<< Exception::eventerror;
	} */
	//cout << "partonsToMatch_.size() (AFTER3) = " << partonsToMatch_.size() << endl;

	return;
	}

	// This looks for all descendents of a top up to but not including
	// the W and b children.
	void FxFxHandler::getTopRadiation(PPtr theParticle) const {
	ParticleVector theChildren(theParticle->children());
	for (unsigned int ixx=0; ixx<theChildren.size(); ixx++)
	if(theChildren[ixx]->children().size()==0)
	tmpList_.push_back(theChildren[ixx]);
	else if(abs(theChildren[ixx]->id())==5\|\|abs(theChildren[ixx]->id())==24)
	return;
	else
	getTopRadiation(theChildren[ixx]);
	return;
	}

	void FxFxHandler::caldel_mg() const {

	preshowerFSPsToDelete_.clear();
	showeredFSPsToDelete_.clear();
	/*
	* Get the MadGraph clustering information
	* from the Les Houches event tags
	*/
	ptclust_.clear();
	getptclust();

	//for(unsigned int izz = 0; izz<ptclust_.size(); izz++) cout << "ptclust_[" << izz << "] = " << ptclust_[izz] << endl;

	/*
	* Get the COM energy of the colliding hadrons
	*
	*/
	getECOM();
	// cout << "ECOM_ = " << ECOM_ << endl;

	for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
	if(ptclust_[ixx] == ECOM_){
	preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
	tmpList_.clear();
	getDescendents(partonsToMatch_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	}

	for(unsigned int ixx=0; ixx<preshowerFSPsToDelete_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
	if(preshowerFSPsToDelete_[ixx]==partonsToMatch_[jxx]) {
	//cout << "deleting " << preshowerFSPsToDelete_[ixx]->id() << endl;
	partonsToMatch_.erase(partonsToMatch_.begin()+jxx);
	break;
	}
	}
	}

	for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
	if(showeredFSPsToDelete_[ixx]==particlesToCluster_[jxx]) {
	// cout << "deleting " << showeredFSPsToDelete_[ixx]->id() << " with parent " << showeredFSPsToDelete_[ixx]->parents()[0]->id() << endl; //" energy = " << preshowerFSPsToDelete_[ixx]->momentum().e()*GeV << endl;
	particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
	break;
	}
	}
	}

	// Sanity check!
	if(partonsToMatch_.size()>particlesToCluster_.size()) {
	throw Exception()
	<< "FxFxHandler::caldel_m() - ERROR!\n"
	<< "No. of ME level partons to be matched to jets = "
	<< partonsToMatch_.size() << "\n"
	<< "No. of showered particles to build jets from = "
	<< particlesToCluster_.size() << "\n"
	<< "There should be at least as many partons to\n"
	<< "cluster as there are partons to match to.\n"
	<< Exception::eventerror;
	}



	}


	void FxFxHandler::caldel_hvq() const {

	// Fill partonsToMatch_ with only those pre-shower partons intended to
	// be used in heavy-quark-jet matching and fill particlesToCluster_ using
	// only those final state particles (post-shower) which are supposed
	// in the heavy-quark-jet clustering used to do merging. To begin with
	// these are made from the corresponding sets of particles that were
	// omitted from the initial jet-parton matching run.
	partonsToMatch_ = preshowerFSPsToDelete_;
	particlesToCluster_.resize(showeredFSPsToDelete_.size());
	for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++)
	particlesToCluster_[ixx] = showeredFSPsToDelete_[ixx];

	// Reset the arrays of particles to delete so that the partonsToMatch_
	// and particlesToCluster_ vectors can undergo further filtering.
	preshowerFSPsToDelete_.clear();
	showeredFSPsToDelete_.clear();

	// Determine further particles in partonsToMatch_ and particlesToCluster_
	// for deletion.
	for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
	// If the progenitor particle is not a heavy
	// quark we delete it and its descendents.
	if(abs(partonsToMatch_[ixx]->id())<4\|\|abs(partonsToMatch_[ixx]->id())>6) {
	preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
	tmpList_.clear();
	getDescendents(partonsToMatch_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	// If the progenitor is a b quark from a top decay drop
	// it & it's descendents too!
	} else if(abs(partonsToMatch_[ixx]->id())==5&&
	partonsToMatch_[ixx]->parents().size()>0&&
	abs(partonsToMatch_[ixx]->parents()[0]->id())==6) {
	preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
	tmpList_.clear();
	getDescendents(partonsToMatch_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	// If (it's a hvy quark not from a top decay and) it has a W/Z/H
	// as a parent [ditto].
	} else if(partonsToMatch_[ixx]->parents().size()>0&&
	(abs(partonsToMatch_[ixx]->parents()[0]->id())==23\|\|
	abs(partonsToMatch_[ixx]->parents()[0]->id())==24\|\|
	abs(partonsToMatch_[ixx]->parents()[0]->id())==25)) {
	preshowerFSPsToDelete_.push_back(partonsToMatch_[ixx]);
	tmpList_.clear();
	getDescendents(partonsToMatch_[ixx]);
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	showeredFSPsToDelete_.push_back(tmpList_[jxx]);
	}
	}

	// Now do the necessary deleting from partonsToMatch_ and
	// particlesToCluster_.
	for(unsigned int ixx=0; ixx<preshowerFSPsToDelete_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<partonsToMatch_.size(); jxx++) {
	if(preshowerFSPsToDelete_[ixx]==partonsToMatch_[jxx]) {
	partonsToMatch_.erase(partonsToMatch_.begin()+jxx);
	break;
	}
	}
	}
	for(unsigned int ixx=0; ixx<showeredFSPsToDelete_.size(); ixx++) {
	for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
	if(showeredFSPsToDelete_[ixx]==particlesToCluster_[jxx]) {
	particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
	break;
	}
	}
	}

	// Now we return to get the decaying top quarks and any
	// radiation they produced.
	ParticleVector intermediates(lastXCombPtr()->subProcess()->intermediates());
	for(unsigned int ixx=0; ixx<intermediates.size(); ixx++) {
	if(abs(intermediates[ixx]->id())==6) {
	partonsToMatch_.push_back(intermediates[ixx]);
	tmpList_.clear();
	getTopRadiation(partonsToMatch_.back());
	for(unsigned int jxx=0; jxx<tmpList_.size(); jxx++)
	particlesToCluster_.push_back(tmpList_[jxx]);
	}
	}

	// If there are any heavy quark progenitors we have to remove
	// the final (showered) instance of them from particlesToCluster.
	ParticleVector evolvedHeavyQuarks;
	for(unsigned int ixx=0; ixx<partonsToMatch_.size(); ixx++) {
	if(abs(partonsToMatch_[ixx]->id())>=4&&abs(partonsToMatch_[ixx]->id())<=6) {
	theProgenitor = partonsToMatch_[ixx];
	// Follow the heavy quark line down to where it stops branching.
	while(theProgenitor->children().size()>0) {
	theLastProgenitor = theProgenitor;
	for(unsigned int jxx=0; jxx<theProgenitor->children().size(); jxx++) {
	if(theProgenitor->children()[jxx]->id()==theProgenitor->id())
	theProgenitor=theProgenitor->children()[jxx];
	}
	// If the progenitor had children but none of them had
	// the same particle id as it, then it must have undergone
	// a decay rather than a branching, i.e. it is the end of
	// the evolution line, so,
	if(theProgenitor==theLastProgenitor) break;
	}
	evolvedHeavyQuarks.push_back(theProgenitor);
	}
	}
	// Now delete the evolved heavy quark from the particlesToCluster.
	for(unsigned int ixx=0; ixx<evolvedHeavyQuarks.size(); ixx++) {
	for(unsigned int jxx=0; jxx<particlesToCluster_.size(); jxx++) {
	if(evolvedHeavyQuarks[ixx]==particlesToCluster_[jxx]) {
	particlesToCluster_.erase(particlesToCluster_.begin()+jxx);
	break;
	}
	}
	}

	return;
	}

	// get npLO_ and npNLO_
	void FxFxHandler::getnpFxFx() const {
	split_vector_type SplitVec;
	// pull the optional weights from the current event
	map<string,double> optionalEventWeights = eventHandler()->currentEvent()->optionalWeights();
	// loop over the optional weights and find np values
	for (map<string,double>::const_iterator it=optionalEventWeights.begin(); it!=optionalEventWeights.end(); ++it){
	// split the line
	boost::split( SplitVec, it->first, boost::is_any_of(" ") );
	// if np is found, store the information
	if(SplitVec[0] == "np") {
	npLO_ = atof(SplitVec[1].c_str());
	npNLO_ = atof(SplitVec[2].c_str());
	}
	}
	return;
	}

	// get hadron COM energy
	void FxFxHandler::getECOM() const {
	split_vector_type SplitVec;
	// pull the optional weights from the current event
	map<string,double> optionalEventWeights = eventHandler()->currentEvent()->optionalWeights();
	// loop over the optional weights and find np values
	for (map<string,double>::const_iterator it=optionalEventWeights.begin(); it!=optionalEventWeights.end(); ++it){
	// split the line
	boost::split( SplitVec, it->first, boost::is_any_of(" ") );
	// if np is found, store the information
	if(SplitVec[0] == "ecom") {
	ECOM_ = it->second;
	}
	}
	return;
	}


	// get pt_clust information
	void FxFxHandler::getptclust() const {
	split_vector_type SplitVec;
	// pull the optional weights from the current event
	map<string,double> optionalEventWeights = eventHandler()->currentEvent()->optionalWeights();
	// loop over the optional weights and find np values
	string str_eq = "=";
	string str_quote = "\"";
	int countptc_(0);
	for (map<string,double>::const_iterator it=optionalEventWeights.begin(); it!=optionalEventWeights.end(); ++it){
	// split the line
	double wgtid = it->second;
	//cout << "wgtdid = " << wgtid << endl;
	if(wgtid == -333) {
	// cout << "it->first for -333 = " << it->first << endl;
	boost::split( SplitVec, it->first, boost::is_any_of(" ") );
	// if np is found, store the information
	double ptclust(0.);
	string stringtohandle("");
	for(unsigned int i_ = 0; i_ < SplitVec.size(); ++i_) {
	if(SplitVec[i_].find("pt_clust_") != std::string::npos) {
	//cout << "pt_clust_ found in " << SplitVec[i_] << endl;
	countptc_++;
	// cout << "SplitVec[i_] = " << SplitVec[i_] << endl;
	stringtohandle = SplitVec[i_];
	stringtohandle.erase(0,10);
	erase_substr(stringtohandle,str_eq);
	erase_substr(stringtohandle,str_quote);
	// cout << "stringtohandle = " << stringtohandle << endl;
	ptclust = atof(stringtohandle.c_str());
	ptclust_.push_back(ptclust);
	}
	}
	}
	}
	return;
	}

	void FxFxHandler::getPreshowerParticles() const {
	// LH file initial-state partons:
	preshowerISPs_ = lastXCombPtr()->subProcess()->incoming();
	// LH file final-state partICLEs:
	preshowerFSPs_ = lastXCombPtr()->subProcess()->outgoing();
	return;
	}

	void FxFxHandler::getShoweredParticles() const {
	// Post-shower initial-state hadrons:
	showeredISHs_ = eventHandler()->currentEvent()->incoming();

	// Post-shower initial-state partons:
	for(unsigned int ixx=0; ixx<(showeredISHs_.first)->children().size(); ixx++)
	if(((showeredISHs_.first)->children()[ixx]->id())<6\|\|
	((showeredISHs_.first)->children()[ixx]->id())==21)
	showeredISPs_.first=(showeredISHs_.first)->children()[ixx];
	for(unsigned int ixx=0; ixx<(showeredISHs_.second)->children().size(); ixx++)
	if(((showeredISHs_.second)->children()[ixx]->id())<6\|\|
	((showeredISHs_.second)->children()[ixx]->id())==21)
	showeredISPs_.second=(showeredISHs_.second)->children()[ixx];

	// Post-shower final-state partICLEs plus remnants (to be removed later):
	showeredFSPs_ = eventHandler()->currentEvent()->getFinalState();

	// Post-shower final-state remnants:
	for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++) {
	if(showeredFSPs_[ixx]->PDGName()=="Rem:p+"\|\|
	showeredFSPs_[ixx]->PDGName()=="Rem:pbar-") {
	if(showeredFSPs_[ixx]->parents()[0]->parents()[0]==
	showeredISHs_.first)
	showeredRems_.first=showeredFSPs_[ixx];
	else if(showeredFSPs_[ixx]->parents()[0]->parents()[0]==
	showeredISHs_.second)
	showeredRems_.second=showeredFSPs_[ixx];
	}
	}

	// Now delete found remnants from the showeredFSPs vector for consistency.
	for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
	if(showeredFSPs_[ixx]->PDGName()=="Rem:p+")
	showeredFSPs_.erase(showeredFSPs_.begin()+ixx);
	for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
	if(showeredFSPs_[ixx]->PDGName()=="Rem:pbar-")
	showeredFSPs_.erase(showeredFSPs_.begin()+ixx);
	sort(showeredFSPs_.begin(),showeredFSPs_.end(),recordEntry);

	return;
	}

	void FxFxHandler::doSanityChecks(int debugLevel) const {

	// When checking momentum conservation in the form
	// p_in - p_out, any momentum component bigger / less
	// than + / - epsilon will result in the p_in - p_out
	// vector being flagged as "non-null", triggering a
	// warning that momentum conservation is violated.
	Energy epsilon(0.5*GeV);
	if(debugLevel>=5) epsilon=1e-9*GeV;

	// Print out what was found for the incoming and outgoing
	// partons in the lastXCombPtr regardless.
	if(debugLevel>=5) {
	cout << "\n\n\n\n";
	cout << "****************************************************" << endl;
	cout << " The following are the hard subprocess momenta from " << "\n"
	<< " lastXCombPtr and should be basically identical to " << "\n"
	<< " the input LH file momenta." << "\n\n";
	cout << " Incoming particles:" << "\n"
	<< *(preshowerISPs_.first) << "\n"
	<< *(preshowerISPs_.second) << endl;
	cout << " Outgoing particles:" << endl;
	for(unsigned int ixx=0; ixx<preshowerFSPs_.size(); ixx++)
	cout << *(preshowerFSPs_[ixx]) << endl;
	}
	// Print out what was found for the incoming and outgoing
	// partons after the shower.
	if(debugLevel>=5) {
	cout << "\n\n";
	cout << "****************************************************" << endl;
	cout << " The following are the particles left at the end of" << "\n"
	<< " the showering step." << "\n\n";
	cout << " Incoming hadrons:" << "\n"
	<< *(showeredISHs_.first) << "\n"
	<< *(showeredISHs_.second) << endl;
	cout << " Incoming partons:" << "\n"
	<< *(showeredISPs_.first) << "\n"
	<< *(showeredISPs_.second) << endl;
	cout << " Outgoing partons:" << endl;
	for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
	cout << *(showeredFSPs_[ixx]) << endl;
	cout << " Outgoing remnants:" << "\n"
	<< *(showeredRems_.first) << "\n"
	<< *(showeredRems_.second) << endl;
	}

	// Check if we correctly identified all initial and final-state
	// particles by testing momentum is conserved.
	if(debugLevel>=4) {
	Lorentz5Momentum tmpMom;
	tmpMom += showeredISPs_.first->momentum();
	tmpMom += showeredISPs_.second->momentum();
	for(unsigned int ixx=0; ixx<showeredFSPs_.size(); ixx++)
	tmpMom -= showeredFSPs_[ixx]->momentum();
	if(!isMomLessThanEpsilon(tmpMom,epsilon))
	cout << "Total parton mom.in - total parton mom.out = "
	<< tmpMom/GeV << endl;
	tmpMom = showeredISHs_.first->momentum()
	- showeredRems_.first->momentum() -showeredISPs_.first->momentum();
	if(!isMomLessThanEpsilon(tmpMom,epsilon))
	cout << "First p_hadron-p_remnant-p_incoming " << tmpMom/GeV << endl;
	tmpMom = showeredISHs_.second->momentum()
	- showeredRems_.second->momentum()-showeredISPs_.second->momentum();
	if(!isMomLessThanEpsilon(tmpMom,epsilon))
	cout << "Second p_hadron-p_remnant-p_incoming " << tmpMom/GeV << endl;
	}

	// Check if what we found to be the remnant is consistent with
	// what we identified as the parent incoming hadron i.e. p+
	// goes with Rem:p+ and pbar- goes with Rem:pbar-.
	if(debugLevel>=0) {
	string tmpString;
	tmpString=showeredRems_.first->PDGName();
	tmpString=tmpString.substr(tmpString.find_first_of(":")+1,
	string::npos);
	if(showeredISHs_.first->PDGName()!=tmpString) {
	cout << "FxFxHandler::showerHardProcessVeto" << "\n"
	<< "Fatal error in pairing of remnant and parent hadron." << "\n"
	<< "Remnant = " << *(showeredRems_.first) << "\n"
	<< "Parent hadron = " << *(showeredISHs_.first)
	<< endl;
	cout << showeredISHs_.first->PDGName() << endl;
	cout << tmpString << endl;
	}
	tmpString=showeredRems_.second->PDGName();
	tmpString=tmpString.substr(tmpString.find_first_of(":")+1,
	string::npos);
	if(showeredISHs_.second->PDGName()!=tmpString) {
	cout << "FxFxHandler::showerHardProcessVeto" << "\n"
	<< "Fatal error in pairing of remnant and parent hadron." << "\n"
	<< "Remnant = " << *(showeredRems_.second) << "\n"
	<< "Parent hadron = " << *(showeredISHs_.second)
	<< endl;
	cout << showeredISHs_.second->PDGName() << endl;
	cout << tmpString << endl;
	}
	}
	return;
	}

	void FxFxHandler::printMomVec(vector<Lorentz5Momentum> momVec) {
	cout << "\n\n";

	// Label columns.
	printf("%5s %9s %9s %9s %9s %9s %9s %9s %9s %9s\n",
	"jet #",
	"px","py","pz","E",
	"eta","phi","pt","et","mass");

	// Print out the details for each jet
	for (unsigned int ixx=0; ixx<momVec.size(); ixx++) {
	printf("%5u %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f\n",
	ixx,
	double(momVec[ixx].x()/GeV),double(momVec[ixx].y()/GeV),
	double(momVec[ixx].z()/GeV),double(momVec[ixx].t()/GeV),
	double(momVec[ixx].eta()) ,double(momVec[ixx].phi()),
	double(momVec[ixx].perp()/GeV),
	double(momVec[ixx].et()/GeV),
	double(momVec[ixx].m()/GeV));
	}

	}

	Energy FxFxHandler::etclusran_(double petc) const {
	return (((2 * epsetclus_)/M_PI) * asin(2 * petc - 1) + etclusmean_);
	}


	void FxFxHandler::erase_substr(std::string& subject, const std::string& search) const {
	size_t pos = 0;
	while((pos = subject.find(search, pos)) != std::string::npos) {
	subject.erase( pos, search.length() );
	}
	}
	diff --git a/MatrixElement/FxFx/FxFxHandler.h b/MatrixElement/FxFx/FxFxHandler.h
	--- a/MatrixElement/FxFx/FxFxHandler.h
	+++ b/MatrixElement/FxFx/FxFxHandler.h
	@@ -1,640 +1,658 @@
	// -- C++ --
	#ifndef HERWIG_FxFxHandler_H
	#define HERWIG_FxFxHandler_H
	//
	// This is the declaration of the FxFxHandler class.
	//

	#include "Herwig/Shower/QTilde/QTildeShowerHandler.h"
	#include "Herwig/Shower/ShowerHandler.h"
	#include "ThePEG/Config/Pointers.h"
	#include "Herwig/Shower/ShowerAlpha.h"
	#include "fastjet/PseudoJet.hh"
	#include "fastjet/ClusterSequence.hh"
	#include "ThePEG/Utilities/CompSelector.h"
	#include "ThePEG/Utilities/XSecStat.h"

	namespace Herwig {
	class FxFxHandler;
	}

	//declaration of thepeg ptr
	namespace ThePEG {
	ThePEG_DECLARE_POINTERS(Herwig::FxFxHandler,FxFxHandlerPtr);
	}

	namespace Herwig {

	using namespace ThePEG;

	typedef vector< string > split_vector_type;

	/**
	* Here is the documentation of the FxFxHandler class.
	*
	* @see \ref FxFxHandlerInterfaces "The interfaces"
	* defined for FxFxHandler.
	*/
	class FxFxHandler: public QTildeShowerHandler {

	/**
	* FxFxHandler should have access to our private parts.
	*/
	friend class FxFxEventHandler;

	friend class FxFxReader;


	public:

	/**
	* The default constructor.
	*/
	FxFxHandler();

	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();

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Finalize the object
	*/
	virtual void dofinish();

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

	public:
	/**
	* Hook to allow vetoing of event after showering hard sub-process
	* as in e.g. MLM merging.
	*/
	virtual bool showerHardProcessVeto() const;

	/**
	* information for FxFx merging
	*/
	mutable int npLO_;
	mutable int npNLO_;

	/**
	* information for tree-level merging
	*/
	mutable vector<double> ptclust_;

	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:

	/*
	* whether a heavy quark has been found in the merging
	*/
	mutable bool hvqfound = false;

	/*
	* Run MLM jet-parton matching on the 'extra' jets.
	*/
	bool lightJetPartonVeto();

	/*
	* Function that calculates deltaR between a parton and a jet
	*/
	double partonJetDeltaR(ThePEG::tPPtr partonptr, LorentzMomentum jetmom) const;


	/*
	* Function that calculates deltaR between two jets
	*/
	double partonJetDeltaR(LorentzMomentum jetmom1, LorentzMomentum jetmom2) const;

	/**
	* Find jets using the FastJet package on particlesToCluster_.
	*/
	void getFastJets(double rjet, Energy ejcut, double etajcut) const;

	/**
	* Find jets using the FastJet package on particlesToCluster_.
	*/
	- void getFastJetsToMatch(double rjet, Energy ejcut, double etajcut) const;
	+ void getFastJetsToMatch(double rjet) const;

	/**
	* Deletes particles from partonsToMatch_ and particlesToCluster_
	* vectors so that these contain only the partons to match to the
	* jets and the particles used to build jets respectively. By and
	* large the candidates for deletion are: vector bosons and their
	* decay products, Higgs bosons, photons as well as _primary_, i.e.
	* present in the lowest multiplicity process, heavy quarks and
	* any related decay products.
	*/
	void caldel_m() const;

	/**
	* Deletes particles from partonsToMatch_ and particlesToCluster_
	* vectors so that these contain only the partons to match to the
	* jets and the particles used to build jets respectively. The candidates
	* are chosen according to the information passed from madgraph.
	*/
	void caldel_mg() const;

	/**
	* c++ translation of subroutine of same name from alpsho.f.
	* Label all particles with status between ISTLO and ISTHI
	* (until a particle with status ISTOP is found) as final-state,
	* call calsim_m and then put labels back to normal. This
	* version keeps only all IST=1 particles rejected by caldel as
	* daughters of vetoed heavy-quark mothers: jets complementary
	* to those reconstructed by caldel.
	*/
	void caldel_hvq() const;

	/**
	* get the MG5_aMC information required for FxFx merging
	*/
	void getnpFxFx() const;

	/**
	* get the MG5_aMC information required for FxFx merging
	*/
	void getECOM() const;


	/**
	* get the MG5_aMC information required for tree-level merging
	*/
	void getptclust() const;

	/**
	* Erases all occurences of a substring from a string
	*/

	void erase_substr(std::string& subject, const std::string& search) const;


	/**
	* Get the particles from lastXCombPtr filling the pair
	* preshowerISPs_ and particle pointer vector preshowerFSPs_.
	*/
	void getPreshowerParticles() const;

	/**
	* Get the particles from eventHandler()->currentEvent()->...
	* filling the particle pairs showeredISHs_, showeredISPs_,
	* showeredRems_ and the particle pointer vector showeredFSPs_.
	*/
	void getShoweredParticles() const;

	/**
	* Allows printing of debug output and sanity checks like
	* total momentum consrvation to be carried out.
	* debugLevel = -1, 0, ...5
	* = no debugging, minimal debugging, ... verbose.
	*/
	void doSanityChecks(int debugLevel) const;

	/**
	* Given a pointer to a particle this finds all its final state
	* descendents.
	*/
	void getDescendents(PPtr theParticle) const;

	/**
	* Accumulates all descendents of tops down to the b and W
	* but not including them.
	*/
	void getTopRadiation(PPtr theParticle) const;

	/**
	* Sorts a given vector of particles by descending pT or ETJET
	*/

	ParticleVector pTsort(ParticleVector unsortedVec);
	pair< vector<Energy>, vector<Lorentz5Momentum> > ETsort(vector<Energy> unsortedetjet, vector<Lorentz5Momentum> unsortedVec);

	/*
	* A function that prints a vector of Lorentz5Momenta in a fancy way
	*/
	void printMomVec(vector<Lorentz5Momentum> momVec);


	/*
	* A probability function for varying etclus_ about the mean value
	*/
	Energy etclusran_(double petc) const;

	private:

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

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

	private:


	/**
	* Initial-state incoming partons prior to showering
	* (i.e. from lastXCombPtr).
	*/
	mutable PPair preshowerISPs_;

	/**
	* Final-state outgoing partICLEs prior to showering
	* (i.e. from lastXCombPtr).
	*/
	mutable ParticleVector preshowerFSPs_;

	/**
	* Final-state outgoing partICLEs prior to showering _to_be_removed_
	* from preShowerFSPs_ prior to the light-parton-light-jet matching
	* step. This same list is the starting point for determining
	* partonsToMatch_ for the case of merging in heavy quark production.
	*/
	mutable ParticleVector preshowerFSPsToDelete_;

	/**
	* Initial-state incoming hadrons after shower of hard process
	* (eventHandler()->currentEvent()->incoming()).
	*/
	mutable PPair showeredISHs_;

	/**
	* Initial-state incoming partons after shower of hard process
	* (look for partonic children of showeredISHs_).
	*/
	mutable PPair showeredISPs_;

	/**
	* Final-state outgoing partICLEs after shower of hard process
	* (eventHandler()->currentEvent()->getFinalState()).
	*/
	mutable tPVector showeredFSPs_;

	/**
	* Final-state outgoing partICLEs after shower of hard process
	* _to_be_removed_ from showeredFSPs_ prior to the
	* light-parton-light-jet matching step. This same list is the
	* starting point for determining particlesToCluster_ for the
	* case of merging in heavy quark production.
	*/
	mutable ParticleVector showeredFSPsToDelete_;

	/**
	* ONLY the final-state partons from preshowerFSPs_ that are
	* supposed to enter the jet-parton matching.
	*/
	mutable ParticleVector partonsToMatch_;

	/*
	* The shower progenitors
	*/

	mutable PPtr theProgenitor;
	mutable PPtr theLastProgenitor;

	/**
	* ONLY the final-state particles from showeredFSPs_ (and maybe
	* also showeredRems_) that are supposed to go for jet clustering.
	*/
	mutable tPVector particlesToCluster_;

	/**
	* Final-state remnants after shower of hard process
	* (look for remnants initially in showeredFSPs_).
	*/
	mutable PPair showeredRems_;

	/**
	* the COM of the incoming hadrons
	*/

	mutable double ECOM_;

	/**
	* Pointer to the object calculating the strong coupling
	*/
	ShowerAlphaPtr alphaS_;

	/**
	* Information extracted from the XComb object
	*/
	//@{
	/**
	* The fixed factorization scale used in the MEs.
	*/
	Energy pdfScale_;

	/**
	* Centre of mass energy
	*/
	Energy2 sHat_;

	/**
	* Constant alphaS used to generate LH events - if not already
	* using CKKW scale (ickkw = 1 in AlpGen for example).
	*/
	double alphaSME_;
	//@}

	/*
	* Number of rapidity segments of the calorimeter.
	*/
	unsigned int ncy_;

	/*
	* Number of phi segments of the calorimeter.
	*/
	unsigned int ncphi_;

	/*
	* Heavy flavour in WQQ,ZQQ,2Q etc (4=c, 5=b, 6=t).
	*/
	int ihvy_;

	/*
	* Number of photons in the AlpGen process.
	*/
	int nph_;

	/*
	* Number of higgses in the AlpGen process.
	*/
	int nh_;

	/*
	* Jet ET cut to apply in jet clustering (in merging).
	*/
	mutable Energy etclus_;

	+
	+
	+ /*
	+ * The merging mode (FxFx vs tree-level) used.
	+ */
	+ int mergemode_;
	+
	+
	+
	+ /*
	+ * Allows the vetoing on heavy quark decay products to be turned off.
	+ */
	+ bool vetoHeavyQ_;
	+
	+ /*
	+ * Allows vetoing of heavy flavour
	+ */
	+
	+ bool vetoHeavyFlavour_;
	+
	+
	/*
	* Mean Jet ET cut to apply in jet clustering (in merging).
	*/
	Energy etclusmean_;

	+
	+ /*
	+ * The jet algorithm used for parton-jet matching in the MLM procedure.
	+ */
	+ int jetAlgorithm_;
	+
	+ /*
	+ * Allows the vetoing to be turned off completely - just for convenience.
	+ */
	+ bool vetoIsTurnedOff_;
	+
	+ /*
	+ * Veto if there exist softer unmatched jets than matched
	+ */
	+
	+ bool vetoSoftThanMatched_;
	+
	+ /*
	+ * This flags whether the etclus_ (merging scale) should be fixed or variable according to a prob. distribution around the mean
	+ */
	+ bool etclusfixed_;
	+
	+
	+
	/*
	* maximum deviation from mean Jet ET cut to apply in jet clustering (in merging).
	*/
	Energy epsetclus_;



	/*
	* Cone size used in jet clustering (in merging).
	*/
	double rclus_;

	/*
	* Max \|eta\| for jets in clustering (in merging).
	*/
	double etaclmax_;

	/*
	* Default 1.5 factor used to decide if a jet matches a parton
	* in merging: if DR(parton,jet)<rclusfactor*rclus the parton
	* and jet are said to have been matched.
	*/
	double rclusfactor_;

	/*
	* Determines whether to detect the hard process or to manually determine which particles
	* to include in the merging. If False, then the ihrd code below is used.
	*/
	bool hpdetect_;

	/*
	* The AlpGen hard process code. Relation to the AlpGen process names:
	* 1: wqq, 2: zqq, 3: wjet, 4: zjet, 5: vbjet, 6: 2Q, 8: QQh, 9: Njet,
	* 10: wcjet, 11: phjet, 12: hjet, 13: top, 14: wphjet, 15: wphqq,
	* 16: 2Qph.
	*/
	int ihrd_;

	/*
	* The number of light jets in the AlpGen process (i.e. the 'extra' ones).
	*/
	int njets_;

	/*
	* Mimimum parton-parton R-sep used for generation (used for hvq merging).
	*/
	double drjmin_;

	/*
	* This flags that the highest multiplicity ME-level process is
	* being processed.
	*/
	mutable bool highestMultiplicity_;

	- /*
	- * This flags whether the etclus_ (merging scale) should be fixed or variable according to a prob. distribution around the mean
	- */
	- bool etclusfixed_;
	-
	- /*
	+ /*
	* The forwards rapidity span of the calorimeter.
	*/
	double ycmax_;

	- /*
	+
	+ /*
	* The backwards rapidity span of the calorimeter.
	*/
	double ycmin_;

	- /*
	- * The jet algorithm used for parton-jet matching in the MLM procedure.
	- */
	- int jetAlgorithm_;

	- /*
	- * The merging mode (FxFx vs tree-level) used.
	- */
	- int mergemode_;

	- /*
	- * Allows the vetoing to be turned off completely - just for convenience.
	- */
	- bool vetoIsTurnedOff_;


	- /*
	- * Allows the vetoing on heavy quark decay products to be turned off.
	- */
	- bool vetoHeavyQ_;
	-
	- /*
	- * Allows vetoing of heavy flavour
	- */
	-
	- bool vetoHeavyFlavour_;


	- /*
	- * Veto if there exist softer unmatched jets than matched
	- */

	- bool vetoSoftThanMatched_;
	+
	+

	+
	/*
	* Cosine of phi values of calorimeter cell centres.
	* Goes phi~=0 to phi~=2*pi (index = 0 ---> ncphi).
	* ==> Cosine goes from +1 ---> +1 (index = 0 ---> ncphi).
	*/
	vector<double> cphcal_;

	/*
	* Sine of phi values of calorimeter cell centres.
	* Goes phi~=0 to phi~=2*pi (index = 0 ---> ncphi).
	* ==> Sine goes 0 -> 1 -> 0 -> -1 -> 0 (index = 0 ---> ncphi).
	*/
	vector<double> sphcal_;

	/*
	* Cosine of theta values of calorimeter cell centres in Y.
	* Goes bwds th~=pi to fwds th~=0 (index = 0 ---> ncy).
	* ==> Cosine goes from -1 ---> +1 (index = 0 ---> ncy).
	*/
	vector<double> cthcal_;

	/*
	* Sine of theta values of calorimeter cell centres in Y.
	* Goes bwds th~=pi to fwds th~=0 (index = 0 ---> ncy).
	* ==> Sine goes from 0 ---> +1 ---> 0 (index = 0 ---> ncy).
	*/
	vector<double> sthcal_;

	/*
	* Transverse energy deposit in a given calorimeter cell.
	* First array index corresponds to rapidity index of cell,
	* second array index corresponds to phi cell index.
	*/
	vector<vector<Energy> > et_;

	/*
	* For a given calorimeter cell this holds the index of the jet
	* that the cell was clustered into.
	*/
	vector<vector<int> > jetIdx_;

	/*
	* Vector holding the Lorentz 5 momenta of each jet.
	*/
	mutable vector<Lorentz5Momentum> pjet_;

	/*
	* Vector holding the Lorentz 5 momenta of each jet from ME partons
	*/
	mutable vector<Lorentz5Momentum> pjetME_;


	/*
	* Vector holding the list of FS particles resulting from
	* the particle input to getDescendents.
	*/
	mutable ParticleVector tmpList_;

	/*
	* Variables for the C++ translation of the calini_m(), calsim_m(),
	* getjet_m(...) and caldel_m() functions
	*/
	mutable vector<Energy> etjet_;
	vector<Energy> etjetME_;
	mutable double dely_, delphi_;

	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the FxFxHandler class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::FxFxHandler>
	: public ClassTraitsBase<Herwig::FxFxHandler> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::FxFxHandler"; }
	/**
	* The name of a file containing the dynamic library where the class
	* FxFxHandler is implemented. It may also include several, space-separated,
	* libraries if the class FxFxHandler depends on other classes (base classes
	* excepted). In this case the listed libraries will be dynamically
	* linked in the order they are specified.
	*/
	static string library() { return "FxFxHandler.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_FxFxHandler_H */
	diff --git a/MatrixElement/FxFx/FxFxReader.h b/MatrixElement/FxFx/FxFxReader.h
	--- a/MatrixElement/FxFx/FxFxReader.h
	+++ b/MatrixElement/FxFx/FxFxReader.h
	@@ -1,1006 +1,1006 @@
	// -- C++ --
	//
	// FxFxReader.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef THEPEG_FxFxReader_H
	#define THEPEG_FxFxReader_H
	// This is the declaration of the FxFxReader class.

	#include "FxFx.h"
	#include "ThePEG/Handlers/HandlerBase.h"
	#include "ThePEG/Utilities/ObjectIndexer.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "ThePEG/Utilities/XSecStat.h"
	#include "ThePEG/PDF/PartonBinInstance.h"
	#include "ThePEG/PDF/PartonBin.fh"
	#include "ThePEG/MatrixElement/ReweightBase.h"
	#include "FxFxEventHandler.fh"
	#include "FxFxReader.fh"
	#include "ThePEG/Utilities/CFile.h"
	#include <cstdio>
	#include <cstring>

	namespace ThePEG {

	/**
	* FxFxReader is an abstract base class to be used for objects
	* which reads event files or streams from matrix element
	* generators. Derived classes must at least implement the open() and
	* doReadEvent() methods to read in information about the whole run into
	* the HEPRUP variable and next event into the HEPEUP variable
	* respectively. Also the close() function to close the file or stream
	* read must be implemented. Although these functions are named as if
	* we are reading from event files, they could just as well implement
	* the actual generation of events.
	*
	* After filling the HEPRUP and HEPEUP variables, which are protected
	* and easily accesible from the sub-class, this base class will then
	* be responsible for transforming this data to the ThePEG Event
	* record in the getEvent() method. <code>FxFxReader</code>s can
	* only be used inside FxFxEventHandler objects.
	*
	* In the initialization the virtual open() and scan() functions are
	* called. Here the derived class must provide the information about
	* the processes in the variables corresponding to the HEPRUP common
	* block. Note that the IDWTUP is required to be +/- 1, and sub
	* classes are required to change the information accordingly to
	* ensure the correct corss section sampling. Note also that the
	* controlling FxFxEventHandler may choose to generate weighted
	* events even if IDWTUP is 1.
	*
	* Note that the information given per process in e.g. the XSECUP and
	* XMAXUP vectors is not used by the FxFxEventHandler and by
	* default the FxFxReader is not assumed to be able to actively
	* choose between the sub-processes. Instead, the
	* FxFxEventHandler can handle several FxFxReader objects
	* and choose between them. However, a sub-class of FxFxReader
	* may set the flag isActive, in which case it is assumed to be able
	* to select between its sub-processes itself.
	*
	* The FxFxReader may be assigned a number ReweightBase objects
	* which either completely reweights the events produced (in the
	* reweights vector), or only biases the selection without influencing
	* the cross section (in the preweights vector). Note that it is the
	* responsibility of a sub-class to call the reweight() function and
	* multiply the weight according to its return value (typically done
	* in the readEvent() function).
	*
	* @see \ref FxFxReaderInterfaces "The interfaces"
	* defined for FxFxReader.
	* @see Event
	* @see FxFxEventHandler
	*/
	class FxFxReader: public HandlerBase, public LastXCombInfo<> {

	/**
	* FxFxEventHandler should have access to our private parts.
	*/
	friend class FxFxEventHandler;

	/**
	* Map for accumulating statistics of cross sections per process
	* number.
	*/
	typedef map<int,XSecStat> StatMap;

	/**
	* Map of XComb objects describing the incoming partons indexed by
	* the corresponding PartonBin pair.
	*/
	typedef map<tcPBPair,XCombPtr> XCombMap;

	/**
	* A vector of pointers to ReweightBase objects.
	*/
	typedef vector<ReweightPtr> ReweightVector;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor. If the optional argument is true, the reader
	* is assumed to be able to produce events on demand for a given
	* process.
	*/
	FxFxReader(bool active = false);

	/**
	* Copy-constructor.
	*/
	FxFxReader(const FxFxReader &);

	/**
	* Destructor.
	*/
	virtual ~FxFxReader();
	//@}

	public:

	/** @name Main virtual fuctions to be overridden in
	* sub-classes. They are named as if we are reading from event
	* files, but could equally well implement the actual generation of
	* events. */
	//@{
	/**
	* Open a file or stream with events and read in the run information
	* into the heprup variable.
	*/
	virtual void open() = 0;

	/**
	* Read the next event from the file or stream into the
	* corresponding protected variables. Return false if there is no
	* more events.
	*/
	virtual bool doReadEvent() = 0;

	/**
	* Close the file or stream from which events have been read.
	*/
	virtual void close() = 0;

	/**
	* return the weight names
	*/
	// virtual vector<string> optWeightsNamesFunc();
	virtual vector<string> optWeightsNamesFunc() = 0;
	//virtual vector<string*> optWeightNamesFunc() = 0;
	vector<string> optionalWeightsNames;

	/**
	* The ID (e.g. 100x, 2001) for the weight
	*/

	// vector<string> optionalWeightsNames;


	//@}

	/** @name Other important function which may be overridden in
	* sub-classes which wants to bypass the basic HEPRUP or HEPEUP
	* variables or otherwise facilitate the conversion to ThePEG
	* objects. */
	//@{
	/**
	* Initialize. This function is called by the FxFxEventHandler
	* to which this object is assigned.
	*/
	virtual void initialize(FxFxEventHandler & eh);

	/**
	* Calls readEvent() or uncacheEvent() to read information into the
	* FxFx common block variables. This function is called by the
	* FxFxEventHandler if this reader has been selectod to
	* produce an event.
	*
	* @return the weight asociated with this event. If negative weights
	* are allowed it should be between -1 and 1, otherwise between 0
	* and 1. If outside these limits the previously estimated maximum
	* is violated. Note that the estimated maximum then should be
	* updated from the outside.
	*/
	virtual double getEvent();

	/**
	* Calls doReadEvent() and performs pre-defined reweightings. A
	* sub-class overrides this function it must make sure that the
	* corresponding reweightings are done.
	*/
	virtual bool readEvent();

	/**
	* Skip \a n events. Used by FxFxEventHandler to make sure
	* that a file is scanned an even number of times in case the events
	* are not ramdomly distributed in the file.
	*/
	virtual void skip(long n);

	/**
	* Get an XComb object. Converts the information in the Les Houches
	* common block variables to an XComb object describing the sub
	* process. This is the way information is conveyed from the reader
	* to the controlling FxFxEventHandler.
	*/
	tXCombPtr getXComb();

	/**
	* Get a SubProcess object corresponding to the information in the
	* Les Houches common block variables.
	*/
	tSubProPtr getSubProcess();

	/**
	* Scan the file or stream to obtain information about cross section
	* weights and particles etc. This function should fill the
	* variables corresponding to the /HEPRUP/ common block. The
	* function returns the number of events scanned.
	*/
	virtual long scan();

	/**
	* Take the information corresponding to the HEPRUP common block and
	* initialize the statistics for this reader.
	*/
	virtual void initStat();

	/**
	* Reweights the current event using the reweights and preweights
	* vectors. It is the responsibility of the sub-class to call this
	* function after the HEPEUP information has been retrieved.
	*/
	double reweight();

	/**
	* Converts the information in the Les Houches common block
	* variables into a Particle objects.
	*/
	virtual void fillEvent();

	/**
	* Removes the particles created in the last generated event,
	* preparing to produce a new one.
	*/
	void reset();

	/**
	* Possibility for subclasses to recover from non-conformant
	* settings of XMAXUP when an event file has been scanned with \a
	* neve events. Should set weightScale so that the average XMAXUP
	* times weightScale gives the cross section for a process. (This is
	* needed for MadEvent).
	*/
	virtual void setWeightScale(long neve);

	//@}

	/** @name Access information about the current event. */
	//@{

	/**
	* Return the size of this event in bytes. To be used for the cache
	* file. \a npart is the number of particles. If \a npart is 0, the
	* number is taken from NUP.
	*/
	static size_t eventSize(int N) {
	return (N + 1)*sizeof(int) + // IDPRUP, ISTUP
	(7N + 4)sizeof(double) + // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP,
	// VTIMUP, SPINUP
	N*sizeof(long) + // IDUP
	2Nsizeof(pair<int,int>) + // MOTHUP, ICOLUP
	sizeof(pair<double,double>) + // XPDWUP.
	2*sizeof(double); // lastweight and preweight
	}

	/**
	* The current event weight given by XWGTUP times possible
	* reweighting. Note that this is not necessarily the same as what
	* is returned by getEvent(), which is scaled with the maximum
	* weight.
	*/
	double eventWeight() const { return hepeup.XWGTUP*lastweight; }

	/**
	* Return the optional named weights associated to the current event.
	*/
	const map<string,double>& optionalEventWeights() const { return optionalWeights; }

	/**
	* Return the optional npLO and npNLO
	*/
	const int& optionalEventnpLO() const { return optionalnpLO; }
	const int& optionalEventnpNLO() const { return optionalnpNLO; }

	/**
	* The pair of PartonBinInstance objects describing the current
	* incoming partons in the event.
	*/
	const PBIPair & partonBinInstances() const { return thePartonBinInstances; }
	/**
	* Return the instances of the beam particles for the current event.
	*/
	const PPair & beams() const { return theBeams; }
	/**
	* Return the instances of the incoming particles to the sub process
	* for the current event.
	*/
	const PPair & incoming() const { return theIncoming; }
	/**
	* Return the instances of the outgoing particles from the sub process
	* for the current event.
	*/
	const PVector & outgoing() const { return theOutgoing; }
	/**
	* Return the instances of the intermediate particles in the sub
	* process for the current event.
	*/
	const PVector & intermediates() const { return theIntermediates; }
	/**
	* 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.
	*/
	int maxMultCKKW() const { return theMaxMultCKKW; }
	/**
	* 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.
	*/
	int minMultCKKW() const { return theMinMultCKKW; } //@}

	/** @name Other inlined access functions. */
	//@{
	/**
	* The number of events found in this reader. If less than zero the
	* number of events are unlimited.
	*/
	long NEvents() const { return theNEvents; }

	/**
	* The number of events produced so far. Is reset to zero if an
	* event file is reopened.
	*/
	long currentPosition() const { return position; }

	/**
	* The maximum number of events to scan to collect information about
	* processes and cross sections. If less than 0, all events will be
	* scanned.
	*/
	long maxScan() const { return theMaxScan; }

	/**
	* Return true if this reader is active.
	*/
	bool active() const { return isActive; }

	/**
	* True if negative weights may be produced.
	*/
	bool negativeWeights() const { return heprup.IDWTUP < 0; }

	/**
	* The collected cross section statistics for this reader.
	*/
	const XSecStat & xSecStats() const { return stats; }

	/**
	* Collected statistics about the individual processes.
	*/
	const StatMap & processStats() const { return statmap; }

	/**
	* Select the current event. It will later be rejected with a
	* probability given by \a weight.
	*/
	void select(double weight) {
	stats.select(weight);
	statmap[hepeup.IDPRUP].select(weight);
	}

	/**
	* Accept the current event assuming it was previously selcted.
	*/
	void accept() {
	stats.accept();
	statmap[hepeup.IDPRUP].accept();
	}

	/**
	* Reject the current event assuming it was previously accepted.
	*/
	void reject(double w) {
	stats.reject(w);
	statmap[hepeup.IDPRUP].reject(w);
	}

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

	/**
	* The PartonExtractor object used to construct remnants.
	*/
	tPExtrPtr partonExtractor() const { return thePartonExtractor; }

	/**
	* Return a possibly null pointer to a CascadeHandler to be used for
	* CKKW-reweighting.
	*/
	tCascHdlPtr CKKWHandler() const { return theCKKW; }

	/**
	* The pairs of PartonBin objects describing the partons which can
	* be extracted by the PartonExtractor object.
	*/
	const PartonPairVec & partonBins() const { return thePartonBins; }

	/**
	* The map of XComb objects indexed by the corresponding PartonBin
	* pair.
	*/
	const XCombMap & xCombs() const { return theXCombs; }

	/**
	* The Cuts object to be used for this reader.
	*/
	const Cuts & cuts() const { return *theCuts; }

	//@}

	protected:

	/** @name Functions for manipulating cache files. */
	//@{

	/**
	* Name of file used to cache the events form the reader in a
	* fast-readable form. If empty, no cache file will be generated.
	*/
	string cacheFileName() const { return theCacheFileName; }

	/**
	* Determines whether to apply cuts to events converting them to
	* ThePEG format.
	*/
	bool cutEarly() const { return doCutEarly; }

	/**
	* File stream for the cache.
	*/
	CFile cacheFile() const { return theCacheFile;}

	/**
	* Open the cache file for reading.
	*/
	void openReadCacheFile();

	/**
	* Open the cache file for writing.
	*/
	void openWriteCacheFile();

	/**
	* Close the cache file;
	*/
	void closeCacheFile();

	/**
	* Write the current event to the cache file.
	*/
	void cacheEvent() const;

	/**
	* Read an event from the cache file. Return false if something went wrong.
	*/
	bool uncacheEvent();

	/**
	* Reopen a reader. If we have reached the end of an event file,
	* reopen it and issue a warning if we have used up a large fraction
	* of it.
	*/
	void reopen();

	/**
	* Helper function to write a variable to a memory location
	*/
	template <typename T>
	static char * mwrite(char * pos, const T & t, size_t n = 1) {
	std::memcpy(pos, &t, n*sizeof(T));
	return pos + n*sizeof(T);
	}

	/**
	* Helper function to read a variable from a memory location
	*/
	template <typename T>
	static const char * mread(const char * pos, T & t, size_t n = 1) {
	std::memcpy(&t, pos, n*sizeof(T));
	return pos + n*sizeof(T);
	}

	//@}

	/** @name Auxilliary virtual methods which may be verridden by sub-classes. */
	//@{
	/**
	* Check the existence of a pair of PartonBin objects corresponding
	* to the current event.
	*
	* @return false if no pair of suitable PartonBin objects was found.
	*/
	virtual bool checkPartonBin();

	/**
	* Create instances of all particles in the event and store them
	* in particleIndex.
	*/
	virtual void createParticles();

	/**
	* Using the already created particles create a pair of
	* PartonBinInstance objects corresponding to the incoming
	* partons. Return the corresponding PartonBin objects.
	*/
	virtual tcPBPair createPartonBinInstances();

	/**
	* Create instances of the incoming beams in the event and store
	* them in particleIndex. If no beam particles are included in the
	* event they are created from the run info.
	*/
	virtual void createBeams();

	/**
	* Go through the mother indices and connect up the Particles.
	*/
	virtual void connectMothers();
	//@}

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

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

	protected:

	/** @name Set functions for some variables not in the Les Houches accord. */
	//@{
	/**
	* The number of events in this reader. If less than zero the number
	* of events is unlimited.
	*/
	void NEvents(long x) { theNEvents = x; }

	/**
	* The map of XComb objects indexed by the corresponding PartonBin
	* pair.
	*/
	XCombMap & xCombs() { return theXCombs; }
	//@}

	/** @name Standard (and non-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() {
	close();
	HandlerBase::dofinish();
	}

	/**
	* Return true if this object needs to be initialized before all
	* other objects because it needs to extract PDFs from the event file.
	*/
	virtual bool preInitialize() const;

	/**
	* Called from doinit() to extract PDFs from the event file and add
	* the corresponding objects to the current EventGenerator.
	*/
	virtual void initPDFs();
	//@}

	protected:

	/**
	* The HEPRUP common block.
	*/
	HEPRUP heprup;

	/**
	* The HEPEUP common block.
	*/
	HEPEUP hepeup;

	/**
	* The ParticleData objects corresponding to the incoming particles.
	*/
	tcPDPair inData;

	/**
	* The PDFBase objects which has been used for the beam particle
	* when generating the events being read. Specified in the interface
	* or derived from PDFGUP and PDFSUP.
	*/
	pair<PDFPtr,PDFPtr> inPDF;

	/**
	* The PDFBase object to be used in the subsequent generation.
	*/
	pair<cPDFPtr,cPDFPtr> outPDF;

	/**
	* The PartonExtractor object used to construct remnants.
	*/
	PExtrPtr thePartonExtractor;

	/**
	* A pointer to a CascadeHandler to be used for CKKW-reweighting.
	*/
	tCascHdlPtr theCKKW;

	/**
	* The pairs of PartonBin objects describing the partons which can
	* be extracted by the PartonExtractor object.
	*/
	PartonPairVec thePartonBins;

	/**
	* The map of XComb objects indexed by the corresponding PartonBin
	* pair.
	*/
	XCombMap theXCombs;

	/**
	* The Cuts object to be used for this reader.
	*/
	CutsPtr theCuts;

	/**
	* The number of events in this reader. If less than zero the number
	* of events is unlimited.
	*/
	long theNEvents;

	/**
	* The number of events produced by this reader so far. Is reset
	* every time an event file is reopened.
	*/
	long position;

	/**
	* The number of times this reader has been reopened.
	*/
	int reopened;

	/**
	* The maximum number of events to scan to collect information about
	* processes and cross sections. If less than 0, all events will be
	* scanned.
	*/
	long theMaxScan;

	/**
	* Flag to tell whether we are in the process of scanning.
	*/
	bool scanning;

	/**
	* True if this is an active reader.
	*/
	bool isActive;

	/**
	* Name of file used to cache the events form the reader in a
	* fast-readable form. If empty, no cache file will be generated.
	*/
	string theCacheFileName;

	/**
	* Determines whether to apply cuts to events before converting them
	* to ThePEG format.
	*/
	bool doCutEarly;

	/**
	* Collect statistics for this reader.
	*/
	XSecStat stats;

	/**
	* Collect statistics for each individual process.
	*/
	StatMap statmap;

	/**
	* The pair of PartonBinInstance objects describing the current
	* incoming partons in the event.
	*/
	PBIPair thePartonBinInstances;

	/**
	* Association between ColourLines and colour indices in the current
	* translation.
	*/
	ObjectIndexer<long,ColourLine> colourIndex;

	/**
	* Association between Particles and indices in the current
	* translation.
	*/
	ObjectIndexer<long,Particle> particleIndex;

	/**
	* The instances of the beam particles for the current event.
	*/
	PPair theBeams;

	/**
	* The instances of the incoming particles to the sub process for
	* the current event.
	*/
	PPair theIncoming;

	/**
	* The instances of the outgoing particles from the sub process for
	* the current event.
	*/
	PVector theOutgoing;

	/**
	* The instances of the intermediate particles in the sub process for
	* the current event.
	*/
	PVector theIntermediates;

	/**
	* File stream for the cache.
	*/
	CFile theCacheFile;

	/**
	* The reweight objects modifying the weights of this reader.
	*/
	ReweightVector reweights;

	/**
	* The preweight objects modifying the weights of this reader.
	*/
	ReweightVector preweights;

	/**
	* The factor with which this reader was last pre-weighted.
	*/
	double preweight;

	/**
	* Should the event be reweighted by PDFs used by the PartonExtractor?
	*/
	bool reweightPDF;

	/**
	* Should PDFBase objects be constructed from the information in the
	* event file in the initialization?
	*/
	bool doInitPDFs;

	/**
	* 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.
	*/
	int theMaxMultCKKW;

	/**
	* 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.
	*/
	int theMinMultCKKW;

	/**
	* The weight multiplying the last read event due to PDF
	* reweighting, CKKW reweighting or assigned reweight and preweight
	* objects.
	*/
	double lastweight;

	/**
	* The optional weights associated to the last read events.
	*/
	map<string,double> optionalWeights;

	+ /**
	+ * If the maximum cross section of this reader has been increased
	+ * with increaseMaxXSec(), this is the total factor with which it
	+ * has been increased.
	+ */
	+ double maxFactor;
	+
	/**
	* npLO for FxFx merging
	*/
	int optionalnpLO;

	/**
	* npNLO for FxFx merging
	*/
	int optionalnpNLO;

	/**
	- * If the maximum cross section of this reader has been increased
	- * with increaseMaxXSec(), this is the total factor with which it
	- * has been increased.
	- */
	- double maxFactor;
	-
	- /**
	* The (reweighted) XWGTUP value should be scaled with this cross
	* section when compared to the overestimated cross section.
	*/
	CrossSection weightScale;

	/**
	* Individual scales for different sub-processes if reweighted.
	*/
	vector<double> xSecWeights;

	/**
	* Individual maximum weights for individual (possibly reweighted)
	* processes.
	*/
	map<int,double> maxWeights;

	/**
	* Is set to true when getEvent() is called from skip(int).
	*/
	bool skipping;

	/**
	* Option for the treatment of the momenta supplied
	*/
	unsigned int theMomentumTreatment;

	/**
	* Set to true if warnings about possible weight incompatibilities
	* should be issued.
	*/
	bool useWeightWarnings;

	/**
	* Option to allow reopening of the file
	*/
	bool theReOpenAllowed;

	/**
	* Use the spin information
	*/
	bool theIncludeSpin;

	private:

	/** Access function for the interface. */
	void setBeamA(long id);
	/** Access function for the interface. */
	long getBeamA() const;
	/** Access function for the interface. */
	void setBeamB(long id);
	/** Access function for the interface. */
	long getBeamB() const;
	/** Access function for the interface. */
	void setEBeamA(Energy e);
	/** Access function for the interface. */
	Energy getEBeamA() const;
	/** Access function for the interface. */
	void setEBeamB(Energy e);
	/** Access function for the interface. */
	Energy getEBeamB() const;
	/** Access function for the interface. */
	void setPDFA(PDFPtr);
	/** Access function for the interface. */
	PDFPtr getPDFA() const;
	/** Access function for the interface. */
	void setPDFB(PDFPtr);
	/** Access function for the interface. */
	PDFPtr getPDFB() const;

	private:

	/**
	* Describe an abstract base class with persistent data.
	*/
	static AbstractClassDescription<FxFxReader> initFxFxReader;

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

	public:

	/** @cond EXCEPTIONCLASSES */
	/** Exception class used by FxFxReader in case inconsistencies
	* are encountered. */
	class FxFxInconsistencyError: public Exception {};

	/** Exception class used by FxFxReader in case more events
	than available are requested. */
	class FxFxReopenWarning: public Exception {};

	/** Exception class used by FxFxReader in case reopening an
	event file fails. */
	class FxFxReopenError: public Exception {};

	/** Exception class used by FxFxReader in case there is
	information missing in the initialization phase. */
	class FxFxInitError: public InitException {};
	/** @endcond */

	};

	/// Stream output for HEPEUP
	ostream & operator<<(ostream & os, const HEPEUP & h);

	}


	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	};

	/** @endcond */

	}

	#endif /* THEPEG_FxFxReader_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 6:33 PM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805614
Default Alt Text: (156 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions