No OneTemporary
Actions

Size

39 KB

Subscribers

None

View Options

	diff --git a/LesHouches/LesHouchesEventHandler.cc b/LesHouches/LesHouchesEventHandler.cc
	--- a/LesHouches/LesHouchesEventHandler.cc
	+++ b/LesHouches/LesHouchesEventHandler.cc
	@@ -1,623 +1,584 @@
	// -- C++ --
	//
	// LesHouchesEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the LesHouchesEventHandler class.
	//

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

	using namespace ThePEG;

	LesHouchesEventHandler::~LesHouchesEventHandler() {}

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

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

	void LesHouchesEventHandler::doinit() {

	EventHandler::doinit();

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

	ntries = 0;
	-
	}

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

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

	// Go through all the readers and collect information about cross
	// sections and processes.
	typedef map<int,tLesHouchesReaderPtr> ProcessMap;
	ProcessMap processes;
	PDPair incoming;
	Energy MaxEA = ZERO;
	Energy MaxEB = ZERO;
	for ( int i = 0, N = readers().size(); i < N; ++i ) {
	LesHouchesReader & reader = *readers()[i];
	reader.initialize(*this);
	+
	+
	+ // TODO check if loop similarity check is needed
	+
	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<LesHouchesInitError>()
	- << "Unknown beam PID " << reader.heprup.IDBMUP.first
	- << ". Have you created a matching BeamParticle object?"
	- << Exception::runerror;
	+ Throw<LesHouchesInitError>()
	+ << "Unknown beam PID " << reader.heprup.IDBMUP.first
	+ << ". Have you created a matching BeamParticle object?"
	+ << Exception::runerror;
	}
	if ( !incoming.second ) {
	incoming.second = getParticleData(reader.heprup.IDBMUP.second);
	if ( !incoming.second )
	- Throw<LesHouchesInitError>()
	- << "Unknown beam PID " << reader.heprup.IDBMUP.first
	- << ". Have you created a matching BeamParticle object?"
	- << Exception::runerror;
	+ Throw<LesHouchesInitError>()
	+ << "Unknown beam PID " << reader.heprup.IDBMUP.first
	+ << ". Have you created a matching BeamParticle object?"
	+ << Exception::runerror;
	}
	- if ( incoming.first->id() != reader.heprup.IDBMUP.first \|\|
	- incoming.second->id() != reader.heprup.IDBMUP.second )
	+ if ( incoming.first->id() != reader.heprup.IDBMUP.first \|\|
	+ incoming.second->id() != reader.heprup.IDBMUP.second )
	Repository::clog()
	- << "The different LesHouchesReader objects in the "
	- << "LesHouchesEventHandler '" << name() << "' have different "
	- << "types of colliding particles." << Exception::warning;
	+ << "The different LesHouchesReader objects in the "
	+ << "LesHouchesEventHandler '" << name() << "' have different "
	+ << "types of colliding particles." << Exception::warning;
	MaxEA = max(MaxEA, reader.heprup.EBMUP.first*GeV);
	MaxEB = max(MaxEB, reader.heprup.EBMUP.second*GeV);

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

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

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

	}

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

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

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

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

	}

	EventPtr LesHouchesEventHandler::generateEvent() {

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

	while ( true ) {
	loopGuard();

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

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

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

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

	select(weight/currentReader()->preweight);
	histStats.select(weight);

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

	accept();

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

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

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

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

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

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

	void LesHouchesEventHandler::skipEvents() {

	// Don't do this for readers which seem to generate events on the fly.
	if ( currentReader()->active() \|\| currentReader()->NEvents() <= 0 ) return;

	// Estimate the fration of the total events available from
	// currentReader() which will be requested.
	double frac = currentReader()->stats.maxXSec()/stats.maxXSec();
	if ( stats.accepted() > 0 )
	frac *= double(stats.attempts())/double(stats.accepted());
	else
	frac *= double(stats.attempts() + 1);
	double xscan = generator()->N()*frac/currentReader()->NEvents();

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

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

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

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

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

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

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

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

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

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

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

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

	if ( statLevel() == 1 ) return;

	if ( statLevel() == 2 ) {

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

	void LesHouchesEventHandler::Init() {

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


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

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

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

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

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

	}

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

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

	namespace ThePEG {


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

	public:

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

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

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

	- friend class LesHouchesHandler;
	-
	-
	public:

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

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

	public:

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

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

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

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

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

	virtual int ntriesinternal() const;

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

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

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

	//@}

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

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

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

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

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

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

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

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

	//@}

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

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

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

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

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


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

	//@}

	public:

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

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

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

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

	protected:

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

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

	protected:

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

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

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

	protected:

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

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

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

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

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

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

	vector<string> weightnames;

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

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

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

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

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


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



	public:

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

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

	private:

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

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

	};

	}

	// CLASSDOC OFF

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#endif /* HERWIG_LesHouchesEventHandler_H */

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions