No OneTemporary
Actions

Size

35 KB

Subscribers

None

View Options

	diff --git a/Handlers/StandardEventHandler.cc b/Handlers/StandardEventHandler.cc
	--- a/Handlers/StandardEventHandler.cc
	+++ b/Handlers/StandardEventHandler.cc
	@@ -1,803 +1,804 @@
	// -- C++ --
	//
	// StandardEventHandler.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 StandardEventHandler class.
	//

	#include "StandardEventHandler.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "ThePEG/Handlers/StdXCombGroup.h"
	#include "ThePEG/Handlers/SubProcessHandler.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Utilities/LoopGuard.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/MatrixElement/MEGroup.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Handlers/LuminosityFunction.h"
	#include "ThePEG/Handlers/SamplerBase.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Config/algorithm.h"
	#include <iomanip>
	#include <sstream>

	using namespace ThePEG;

	StandardEventHandler::StandardEventHandler()
	: EventHandler(false), collisionCuts(true), theBinStrategy(2), weightedEvents(false),
	theLumiDim(0) {
	setupGroups();
	}

	StandardEventHandler::~StandardEventHandler() {}

	void StandardEventHandler::reject(double w) {
	tStdXCombPtr last = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
	if ( !last ) return;
	last->reject(w);
	}

	void StandardEventHandler::reweight(double factor) const {
	tStdXCombPtr last = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
	if ( !currentEvent() \|\| !last )
	return;
	double weight = currentEvent()->weight();
	last->reweight(weight,factor*weight);
	currentEvent()->weight(factor*weight);
	for ( map<string,double>::iterator w =
	currentEvent()->optionalWeights().begin();
	w != currentEvent()->optionalWeights().end(); ++w )
	w->second *= factor;
	}

	void StandardEventHandler::doupdate() {
	EventHandler::doupdate();
	bool redo = touched();
	UpdateChecker::check(theIncomingA, redo);
	UpdateChecker::check(theIncomingB, redo);
	for_each(subProcesses(), UpdateChecker(redo));
	if ( !redo ) return;
	theIncoming.first = theIncomingA;
	theIncoming.second = theIncomingB;
	for ( SubHandlerList::iterator sit = subProcesses().begin();
	sit != subProcesses().end(); ++sit )
	if ( !(**sit).pExtractor()->canHandle(incoming()) )
	throw StandardEventHandlerUpdateException()
	<< "Cannot use the parton extractor '" << (**sit).pExtractor()->name()
	<< "' in the SubProcessHandler '" << (**sit).name() << "' in the "
	<< "StandardEventHandler '" << name() << "' since it cannot handle "
	<< "the requested types of incoming particles ("
	<< theIncomingA->name() << "," << theIncomingB->name() << ").";
	}

	void StandardEventHandler::doinit() {
	EventHandler::doinit();
	if ( !lumiFnPtr() ) throw StandardEventHandlerUpdateException()
	<< "The StandardEventHandler '" << name() << "' does not have any "
	<< "LuminosityFunction object assigned to it, which it needs to be "
	<< "able to generate events." << Exception::warning;
	}

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

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

	void StandardEventHandler::
	addME(Energy maxEnergy, tSubHdlPtr sub, tPExtrPtr extractor, tCutsPtr cuts,
	tCascHdlPtr ckkw, tMEPtr me, const PBPair & pBins,
	const PartonPairVec& allPBins) {
	typedef MEBase::DiagramVector DiagramVector;
	typedef map<string,DiagramVector> DiagramMap;
	cPDPair pin(pBins.first->parton(), pBins.second->parton());
	DiagramVector diag = me->diagrams();
	DiagramMap tdiag;
	DiagramMap tmdiag;
	for ( int i = 0, N = diag.size(); i < N; ++i ) {
	cPDPair din(diag[i]->partons()[0], diag[i]->partons()[1]);
	if (!me->noMirror())
	if ( din.first->id() < din.second->id() ) swap(din.first, din.second);
	if ( din == pin ) tdiag[diag[i]->getTag()].push_back(diag[i]);
	if (!me->noMirror())
	if ( din.first == pin.second && din.second == pin.first )
	tmdiag[diag[i]->getTag()].push_back(diag[i]);
	}

	if ( tdiag.empty() ) tdiag = tmdiag;
	for ( DiagramMap::iterator dit = tdiag.begin(); dit != tdiag.end(); ++dit ) {
	cPDPair din(dit->second.back()->partons()[0],
	dit->second.back()->partons()[1]);
	// check
	assert(me->noMirror() ? din == pin : true);
	StdXCombPtr xcomb = me->makeXComb(maxEnergy, incoming(), this, sub, extractor,
	ckkw, pBins, cuts, dit->second, din != pin,
	allPBins);
	if ( xcomb->checkInit() ) xCombs().push_back(xcomb);
	else generator()->logWarning(
	StandardEventHandlerInitError() << "The matrix element '"
	<< xcomb->matrixElement()->name()
	<< "' cannot generate the diagram '"
	<< dit->first
	<< "' when used together with the parton extractor '"
	<< xcomb->pExtractor()->name()
	<< "'.\nThe corresponding diagram is switched off, "
	<< "check the collision energy and/or the cuts."
	<< Exception::warning);
	}
	}

	void StandardEventHandler::initGroups() {
	tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
	if ( lastXC ) optGroups = lastXC->subProcessHandler()->groups();
	EventHandler::initGroups();
	}

	tCollPtr StandardEventHandler::performCollision() {
	tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
	if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr());
	lastExtractor()->select(lastXC);
	currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this)));
	if ( currentEvent() ) currentEvent()->addCollision(currentCollision());
	currentStep(new_ptr(Step(currentCollision())));
	currentCollision()->addStep(currentStep());

	currentStep()->addSubProcess(lastXC->construct());

	lastExtractor()->construct(lastXC->partonBinInstances(), currentStep());
	if ( collisionCuts )
	if ( !lastCuts().passCuts(*currentCollision()) ) throw Veto();
	initGroups();
	if ( ThePEG_DEBUG_ITEM(1) ) {
	if ( currentEvent() )
	generator()->logfile() << *currentEvent();
	else
	generator()->logfile() << *currentCollision();
	}
	return continueCollision();
	}

	void StandardEventHandler::setScale(Energy2 scale) {
	lastXCombPtr()->lastScale(scale);
	}

	void StandardEventHandler::initialize() {

	theLumiDim = lumiFn().nDim(incoming());
	Energy maxEnergy = lumiFn().maximumCMEnergy();

	xCombs().clear();
	xSecs().clear();

	cuts()->initialize(sqr(maxEnergy), lumiFn().Y());

	for ( SubHandlerList::const_iterator sit = subProcesses().begin();
	sit != subProcesses().end(); ++sit ) {
	CutsPtr kincuts = (sit).cuts()? (sit).cuts(): cuts();
	if ( (**sit).cuts() ) kincuts->initialize(sqr(maxEnergy), lumiFn().Y());
	PExtrPtr pextract = (**sit).pExtractor();

	tCascHdlPtr ckkw = (**sit).CKKWHandler();
	if ( !ckkw ) ckkw = CKKWHandler();

	PartonPairVec vpc = pextract->getPartons(maxEnergy, incoming(), *kincuts);

	// The last parton bin pair was in fact the bins corresponding to
	// the incoming particles, so we remove them, but save them to
	// keep them referenced.
	PBPair orig = vpc.back();
	vpc.pop_back();
	for ( PartonPairVec::iterator ppit = vpc.begin();
	ppit != vpc.end(); ++ppit )
	for ( MEVector::const_iterator meit = (**sit).MEs().begin();
	meit != (**sit).MEs().end(); ++meit ) {
	addME(maxEnergy, sit, pextract, kincuts, ckkw, meit, *ppit,vpc);
	}
	}
	xSecs().resize(xCombs().size());

	theMaxDims.clear();
	switch ( binStrategy() ) {
	case 0: {
	theMaxDims.push_back(0);
	for ( int i = 0, N = xCombs().size(); i < N; ++i )
	theMaxDims[0] = max(theMaxDims[0], xCombs()[i]->nDim());
	break;
	}
	case 1: {
	for ( int i = 0, N = xCombs().size(); i < N; ++i )
	theMEXMap[xCombs()[i]->matrixElement()].push_back(xCombs()[i]);
	MEXMap::const_iterator mei = theMEXMap.begin();
	for ( int i = 0, N = theMEXMap.size(); i < N; ++i, ++mei) {
	theMaxDims.push_back(0);
	for ( int j = 0, M = mei->second.size(); j < M; ++j )
	theMaxDims[i] = max(theMaxDims[i], mei->second[j]->nDim());
	}
	break;
	}
	case 2: {
	for ( int i = 0, N = xCombs().size(); i < N; ++i )
	theMaxDims.push_back(xCombs()[i]->nDim());
	break;
	}
	}

	sampler()->setEventHandler(this);
	sampler()->initialize();
	}

	CrossSection StandardEventHandler::
	dSigDR(const pair<double,double> ll, Energy2 maxS,
	int ibin, int nr, const double * r) {
	PPair inc = make_pair(incoming().first->produceParticle(),
	incoming().second->produceParticle());
	SimplePhaseSpace::CMS(inc, maxS);

	XVector xv;
	switch ( binStrategy() ) {
	case 0:
	xv = xCombs();
	break;
	case 1: {
	MEXMap::iterator mei = theMEXMap.begin();
	for ( int i = 0; i < ibin; ++i) ++mei;
	xv = mei->second;
	break;
	}
	case 2:
	xv = XVector(1, xCombs()[ibin]);
	break;
	}

	xSecs().resize(xv.size());
	for ( int i = 0, N = xv.size(); i < N; ++i ) xv[i]->prepare(inc);
	CrossSection sum = ZERO;
	for ( int i = 0, N = xv.size(); i < N; ++i )
	xSecs()[i] = ( sum += xv[i]->dSigDR(ll, nr, r) );
	return sum;
	}

	tStdXCombPtr StandardEventHandler::select(int bin, double & weight) {

	int i = 0;
	if ( binStrategy() != 2 )
	i = upper_bound(xSecs().begin(), xSecs().end(), rnd()*xSecs().back())
	- xSecs().begin();
	tStdXCombPtr lastXC;
	switch ( binStrategy() ) {
	case 0:
	lastXC = xCombs()[i];
	break;
	case 1: {
	MEXMap::iterator mei = theMEXMap.begin();
	for ( int j = 0; j < bin; ++j) ++mei;
	lastXC = mei->second[i];
	break;
	}
	case 2:
	lastXC = xCombs()[bin];
	break;
	}
	// clean up the old XComb object before switching to a new one
	if ( theLastXComb && theLastXComb != lastXC ) theLastXComb->clean();
	theLastXComb = lastXC;
	lastXC->matrixElement()->setXComb(lastXC);
	weight /= lastXC->matrixElement()->preWeight();
	lastXC->select(weight);
	lastXC->accept();
	return lastXC;
	}

	int StandardEventHandler::nBins() const {
	switch ( binStrategy() ) {
	case 0: return 1;
	case 1: return theMEXMap.size();
	case 2: return xCombs().size();
	}
	return -1;
	}

	struct Stat {

	Stat() : attempted(0), accepted(0), sumw(0.0), sumw2(),
	maxXSec(CrossSection()), totsum(0.0), totrerr(0.0) {}
	Stat(long att, long acc, double w, double w2, CrossSection x,
	double sumw, double genrerr)
	: attempted(att), accepted(acc), sumw(w), sumw2(w2), maxXSec(x),
	totsum(sumw), totrerr(genrerr) {}

	inline CrossSection xSec() const {
	return totsum != 0.0? maxXSec*sumw/totsum: maxXSec;
	}

	inline CrossSection xSecErr() const {
	if ( totsum == 0.0 ) return maxXSec;
	if ( sumw == 0.0 ) return xSec();
	return abs(xSec())*sqrt(sqr(totrerr) + sumw2/sqr(sumw));
	}

	long attempted;
	long accepted;
	double sumw;
	double sumw2;
	CrossSection maxXSec;
	double totsum;
	double totrerr;

	const Stat & operator+= (const Stat & s) {
	attempted += s.attempted;
	accepted += s.accepted;
	sumw += s.sumw;
	sumw2 += s.sumw2;

	// replaced by taking whatever is in the one to be added, see
	// comments below; SP 2013/03/14

	/*
	totsum = max(totsum, s.totsum);
	if ( totsum != 0.0 )
	maxXSec = max(maxXSec, s.maxXSec);
	else
	maxXSec += s.maxXSec;
	*/

	// whereever Stat is used nothing happens with the above code,
	// because those quantities are always the same for all Stat
	// objects involved, so long as cross sections and sum of weights
	// are positive; but it fails, if we have negative cross sections
	// (which can happen when running real-subtraction of a NLO
	// without the Born and virtual contributions). To this extent
	// it's commented out and replaced by just taking whatever is in
	// the statistics to be added. There are more issues in the
	// statistics treatment going on here which urgently need to be
	// looked at, and the Stat object itself needs to be made much
	// more transparent and cleaned up from misleading constructions
	// like this and the old one below the next two lines.

	totsum = s.totsum;
	maxXSec = s.maxXSec;

	totrerr = s.totrerr;
	return *this;
	}
	};

	void StandardEventHandler::statistics(ostream & os) const {
	if ( statLevel() == 0 ) return;
	map<cPDPair, Stat> partonMap;
	map<MEPtr, Stat> meMap;
	map<PExtrPtr, Stat> extractMap;
	Stat tot;
	double genrerr = sampler()->integratedXSecErr()/sampler()->integratedXSec();

	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	Stat s(x.stats().attempts(), x.stats().accepted(),
	x.stats().sumWeights(), x.stats().sumWeights2(),
	sampler()->integratedXSec(), sampler()->sumWeights(), genrerr);
	partonMap[x.partons()] += s;
	meMap[x.matrixElement()] += s;
	extractMap[x.pExtractor()] += s;
	tot += s;
	}

	string line = string(78, '=') + "\n";

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

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

	os << line << "Total (from weighted events): including vetoed events" << setw(23)
	<< ouniterr(sampler()->integratedXSec(),
	sampler()->integratedXSecErr(), nanobarn)
	<< endl;
	- os << line << "Total (from unweighted events):" << setw(17) << tot.accepted << setw(13)
	+ os << line << "Total (from "
	+ << (weighted() ? " weighted" : "unweighted") << " events):"
	+ << setw(17) << tot.accepted << setw(13)
	<< tot.attempted << setw(17)
	<< ouniterr(tot.xSec(),tot.xSecErr() , nanobarn)
	<< endl << line;

	if ( statLevel() == 1 ) return;

	os << "Per matrix element breakdown:\n";
	for ( map<MEPtr, Stat>::iterator i = meMap.begin();
	i != meMap.end(); ++i ) {
	string n = i->first->name();
	n.resize(37, ' ');
	os << n << setw(11) << i->second.accepted << setw(13)
	<< i->second.attempted << setw(17)
	<< ouniterr(i->second.xSec(), i->second.xSecErr(), nanobarn)
	<< endl;
	}
	os << line;

	if ( statLevel() == 2 ) return;

	os << "Per parton extractor breakdown:\n";
	for ( map<PExtrPtr, Stat>::iterator i = extractMap.begin();
	i != extractMap.end(); ++i ) {
	string n = i->first->name();
	n.resize(37, ' ');
	os << n << setw(11) << i->second.accepted << setw(13)
	<< i->second.attempted << setw(17)
	<< ouniterr(i->second.xSec(), i->second.xSecErr(), nanobarn)
	<< endl;
	}
	os << line;

	os << "Per incoming partons breakdown:\n";
	for ( map<cPDPair, Stat>::iterator i = partonMap.begin();
	i != partonMap.end(); ++i ) {
	string n = i->first.first->PDGName() + " " + i->first.second->PDGName();
	n.resize(37, ' ');
	os << n << setw(11) << i->second.accepted << setw(13)
	<< i->second.attempted << setw(17)
	<< ouniterr(i->second.xSec(), i->second.xSecErr(), nanobarn)
	<< endl;
	}
	os << line;

	if ( statLevel() == 3 ) return;

	os << "Detailed breakdown:\n";
	CrossSection xsectot = sampler()->integratedXSec()/sampler()->sumWeights();

	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	os << "(" << x.pExtractor()->name() << ") "
	<< x.partons().first->PDGName() << " "
	<< x.partons().second->PDGName()

	<< " (" << x.matrixElement()->name() << " "
	<< x.lastDiagram()->getTag() << ") " << endl
	<< setw(48) << x.stats().accepted() << setw(13) << x.stats().attempts()
	<< setw(17)
	<< ouniterr(x.stats().sumWeights()*xsectot,
	sqrt(x.stats().sumWeights2())*xsectot, nanobarn) << endl;
	}

	os << line;

	}

	CrossSection StandardEventHandler::histogramScale() const {
	Stat tot;
	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	Stat s;
	s = Stat(x.stats().attempts(), x.stats().accepted(),
	x.stats().sumWeights(), x.stats().sumWeights2(),
	sampler()->integratedXSec(), sampler()->sumWeights(), 1.0);
	tot += s;
	}

	return tot.xSec()/tot.sumw;
	}

	CrossSection StandardEventHandler::integratedXSec() const {
	if ( sampler()->integratedXSec() == ZERO )
	return sampler()->maxXSec();

	Stat tot;
	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	Stat s;
	s = Stat(x.stats().attempts(), x.stats().accepted(),
	x.stats().sumWeights(), x.stats().sumWeights2(),
	sampler()->integratedXSec(), sampler()->sumWeights(), 1.0);
	tot += s;
	}

	return tot.xSec();
	-
	}

	CrossSection StandardEventHandler::integratedXSecErr() const {
	if ( sampler()->integratedXSec() == ZERO )
	return sampler()->maxXSec();

	Stat tot;
	double genrerr = sampler()->integratedXSecErr()/sampler()->integratedXSec();
	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	Stat s;
	s = Stat(x.stats().attempts(), x.stats().accepted(),
	x.stats().sumWeights(), x.stats().sumWeights2(),
	sampler()->integratedXSec(), sampler()->sumWeights(), genrerr);
	tot += s;
	}

	return tot.xSecErr();

	}

	CrossSection StandardEventHandler::integratedXSecNoReweight() const {
	if ( sampler()->integratedXSec() == ZERO )
	return sampler()->maxXSec();

	Stat tot;
	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	Stat s;
	s = Stat(x.stats().attempts(), x.stats().accepted(),
	x.stats().sumWeightsNoReweight(), x.stats().sumWeights2NoReweight(),
	sampler()->integratedXSec(), sampler()->sumWeights(), 1.0);
	tot += s;
	}

	return tot.xSec();
	}

	CrossSection StandardEventHandler::integratedXSecErrNoReweight() const {
	if ( sampler()->integratedXSec() == ZERO )
	return sampler()->maxXSec();

	Stat tot;
	double genrerr = sampler()->integratedXSecErr()/sampler()->integratedXSec();
	for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
	const StandardXComb & x = *xCombs()[i];
	Stat s;
	s = Stat(x.stats().attempts(), x.stats().accepted(),
	x.stats().sumWeightsNoReweight(), x.stats().sumWeights2NoReweight(),
	sampler()->integratedXSec(), sampler()->sumWeights(), genrerr);
	tot += s;
	}

	return tot.xSecErr();

	}


	void StandardEventHandler::doinitrun() {
	EventHandler::doinitrun();
	for ( SubHandlerList::iterator sit = subProcesses().begin();
	sit != subProcesses().end(); ++sit )
	(**sit).initrun();
	sampler()->initrun();
	for ( int i = 0, N = xCombs().size(); i < N; ++i )
	xCombs()[i]->reset();
	}

	CrossSection StandardEventHandler::dSigDR(const vector<double> & r) {
	double jac = 1.0;
	pair<double,double> ll = lumiFn().generateLL(&r[0], jac);
	Energy2 maxS = sqr(lumiFn().maximumCMEnergy())/exp(ll.first + ll.second);
	int bin = sampler()->lastBin();
	CrossSection x = jac*lumiFn().value(incoming(), ll.first, ll.second)
	*dSigDR(ll, maxS, bin, nDim(bin) - lumiDim(), &r[lumiDim()]);
	return x;
	}

	EventPtr StandardEventHandler::generateEvent() {

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

	while (1) {
	loopGuard();

	EventHandler::clean();

	double weight = sampler()->generate();

	tStdXCombPtr lastXC = select(sampler()->lastBin(), weight);

	try {

	lumiFn().select(lastXC);

	currentEventBoost() = lumiFn().getBoost();

	currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
	generator()->currentEventNumber(), weight)));

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

	currentEvent()->transform(currentEventBoost());

	return currentEvent();

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

	EventPtr StandardEventHandler::continueEvent() {
	if ( !generator() ) throw StandardEventHandlerInitError()
	<< "The event handler '" << name() << "' had not been isolated "
	<< "from the setup phase before it was used." << Exception::maybeabort;
	try {
	continueCollision();
	}
	catch (Veto) {
	reject(currentEvent()->weight());
	}
	catch (Stop) {
	}
	catch (Exception &) {
	reject(currentEvent()->weight());
	throw;
	}
	return currentEvent();
	}

	void StandardEventHandler::select(tXCombPtr newXComb) {
	EventHandler::select(newXComb);
	lastExtractor()->select(newXComb);
	}

	void StandardEventHandler::clean() {
	if ( theLastXComb ) theLastXComb->clean();
	for (size_t i=0; i < theXCombs.size(); ++i )
	if ( theXCombs[i] ) theXCombs[i]->clean();
	EventHandler::clean();
	}

	void StandardEventHandler::dofinish() {
	clean();
	EventHandler::dofinish();
	}

	ClassDescription<StandardEventHandler>
	StandardEventHandler::initStandardEventHandler; //

	void StandardEventHandler::setIncomingA(PDPtr p) {
	theIncomingA = p;
	theIncoming.first = p;
	}

	void StandardEventHandler::setIncomingB(PDPtr p) {
	theIncomingB = p;
	theIncoming.second = p;
	}

	void StandardEventHandler::Init() {

	static ClassDocumentation<StandardEventHandler> documentation
	("This is the standard event handler to generate hard sub-processes "
	"within ThePEG. It must specify a pair of incoming particle beams "
	"in <interface>BeamA</interface> and <interface>BeamB</interface> "
	"and a suiteable <interface>LuminosityFunction</interface>. In "
	"addition at least one object describing the sub-processes to be "
	"generated must be specified in "
	"<interface>SubProcessHandlers</interface>.");

	static Reference<StandardEventHandler,ParticleData> interfaceIncomingA
	("BeamA",
	"The type of particles in first beam",
	&StandardEventHandler::theIncomingA, false, false, true, false,
	&StandardEventHandler::setIncomingA);

	static Reference<StandardEventHandler,ParticleData> interfaceIncomingB
	("BeamB",
	"The type of particles in second beam",
	&StandardEventHandler::theIncomingB, false, false, true, false,
	&StandardEventHandler::setIncomingB);

	static RefVector<StandardEventHandler,SubProcessHandler> interfaceSubhandlers
	("SubProcessHandlers",
	"The list of sub-process handlers used in this StandardEventHandler. ",
	&StandardEventHandler::theSubProcesses, 0, false, false, true, false);

	static Reference<StandardEventHandler,Cuts> interfaceCuts
	("Cuts",
	"Common kinematical cuts for this StandardEventHandler. These cuts "
	"may be overidden in individual sub-process handlers.",
	&StandardEventHandler::theCuts, false, false, true, false);

	static Switch<StandardEventHandler,bool> interfaceCollisionCuts
	("CollisionCuts",
	"Switch on or off cuts on collision objects",
	&StandardEventHandler::collisionCuts, true, false, false);
	static SwitchOption interfaceCollisionCutsOn
	(interfaceCollisionCuts,
	"On",
	"Switch on cuts on collision objects",
	true);
	static SwitchOption interfaceCollisionCutsOff
	(interfaceCollisionCuts,
	"Off",
	"Switch off cuts on collision cuts",
	false);

	static Switch<StandardEventHandler,int> interfaceBinStrategy
	("BinStrategy",
	"The strategy to be used when sampling different ThePEG::XComb "
	"objects. An ThePEG::XComb objet represents a pair of incoming "
	"parton types as defined by a THePEG::PartonExtractor and a "
	"matrix element.",
	&StandardEventHandler::theBinStrategy, 2, false, false);

	static SwitchOption interfaceBinStrategy0
	(interfaceBinStrategy,
	"AllAtOnce",
	"All bins are sampled together.",
	0);

	static SwitchOption interfaceBinStrategy1
	(interfaceBinStrategy,
	"PerME",
	"All bins which have the same matrix element object are sampled together.",
	1);

	static SwitchOption interfaceBinStrategy2
	(interfaceBinStrategy,
	"Individual",
	"All bins are sampled individually.",
	2);

	static Switch<StandardEventHandler,bool> interfaceWeighted
	("Weighted",
	"If switched on, this event Handler will produce weighted events",
	&StandardEventHandler::weightedEvents, false, false);
	static SwitchOption interfaceWeightedTrue
	(interfaceWeighted, "On",
	"This EventHandler produces weighted events.", true);
	static SwitchOption interfaceWeightedFalse
	(interfaceWeighted, "Off",
	"This EventHandler produces unweighted events.", false);

	static Reference<StandardEventHandler,SamplerBase> interfaceSampler
	("Sampler",
	"The phase space sampler responsible for generating phase space"
	"points according to the cross section given by this event handler",
	&StandardEventHandler::theSampler, false, false, true, true);

	interfaceSubhandlers.rank(11);
	interfaceIncomingA.rank(3);
	interfaceIncomingB.rank(2);

	}

	void StandardEventHandler::persistentOutput(PersistentOStream & os) const {
	os << theIncomingA << theIncomingB << theSubProcesses << theCuts << collisionCuts
	<< theXCombs << ounit(theXSecs, nanobarn)
	<< theBinStrategy << theMaxDims << theMEXMap
	<< weightedEvents
	<< theSampler << theLumiDim;
	}

	void StandardEventHandler::persistentInput(PersistentIStream & is, int) {
	is >> theIncomingA >> theIncomingB >> theSubProcesses >> theCuts >> collisionCuts
	>> theXCombs >> iunit(theXSecs, nanobarn)
	>> theBinStrategy >> theMaxDims>> theMEXMap
	>> weightedEvents
	>> theSampler >> theLumiDim;
	}

	diff --git a/Utilities/UnitIO.h b/Utilities/UnitIO.h
	--- a/Utilities/UnitIO.h
	+++ b/Utilities/UnitIO.h
	@@ -1,297 +1,305 @@
	// -- C++ --
	//
	// UnitIO.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_UnitIO_H
	#define ThePEG_UnitIO_H
	// This is the declaration of the IUnit and OUnit classes and
	// associated templated functions.

	#include <complex>
	#include <iomanip>
	#include <sstream>
	#include <cstdlib>
	#include <cmath>

	// Workarounds for OS X
	#if defined __APPLE__ && defined __MACH__
	extern "C" int isnan(double) throw();
	extern "C" int isinf(double) throw();
	#endif

	namespace ThePEG {

	/**
	* The OUnit< class is used to
	* facilitate output of unitful numbers to a
	* persistent stream. An Energy can hence be written like
	* this:<BR> <code>os
	* << ounit(x, GeV);</code><BR> Also containers of unitful
	* numbers may be written like this, as well as LorentzVector and
	* ThreeVector.
	*
	* @see PersistentOStream
	* @see PersistentIStream
	*
	*/
	template <typename T, typename UT>
	struct OUnit {

	/** Constructor given an object to be written assuming the given
	* unit. */
	OUnit(const T & t, const UT & u): theX(t), theUnit(u) {}

	/** Copy constructor */
	OUnit(const OUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}

	/** Reference to the object to be written. */
	const T & theX;

	/** The unit assumed when writing the object. */
	const UT & theUnit;
	};

	/**
	* The IUnit class is used to facilitate input of unitful numbers from
	* and to a persistent stream. An Energy can hence be read like
	* this:<BR> <code>is >> iunit(x, GeV);</code><BR> Also containers of
	* unitful numbers may be read like this, as well as LorentzVector and
	* ThreeVector.
	*
	* @see PersistentOStream
	* @see PersistentIStream
	*
	*/
	template <typename T, typename UT>
	struct IUnit {

	/** Constructor given an object to be read assuming the given
	* unit. */
	IUnit(T & t, const UT & u): theX(t), theUnit(u) {}

	/** Copy constructor */
	IUnit(const IUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}

	/** Reference to the object to be read. */
	T & theX;

	/** The unit assumed when reading the object. */
	const UT & theUnit;

	};

	/** Helper function creating a OUnit object given an object and a
	* unit. */
	template <typename T, typename UT>
	inline OUnit<T,UT> ounit(const T & t, const UT & ut) {
	return OUnit<T,UT>(t, ut);
	}

	/** Helper function creating a IUnit object given an object and a
	* unit. */
	template <typename T, typename UT>
	inline IUnit<T,UT> iunit(T & t, const UT & ut) {
	return IUnit<T,UT>(t, ut);
	}

	/** Helper function writing out an object with a given unit to an
	* output stream. */
	template <typename OStream, typename T, typename UT>
	void ounitstream(OStream & os, const T & t, UT & u) {
	os << t/u;
	}

	/** Helper function reading an object with a given unit from an
	* input stream. */
	template <typename IStream, typename T, typename UT>
	void iunitstream(IStream & is, T & t, UT & u) {
	double d;
	is >> d;
	t = d*u;;
	}

	/** Helper function reading a complex object with a given unit from an
	* input stream. */
	template <typename IStream, typename T, typename UT>
	void iunitstream(IStream & is, std::complex<T> & t, UT & u) {
	std::complex<double> d;
	is >> d;
	t = d*u;;
	}

	/** Output an OUnit object to a stream. */
	template <typename OStream, typename T, typename UT>
	OStream & operator<<(OStream & os, const OUnit<T,UT> & u) {
	ounitstream(os, u.theX, u.theUnit);
	return os;
	}

	/** Input an IUnit object from a stream. */
	template <typename IStream, typename T, typename UT>
	IStream & operator>>(IStream & is, const IUnit<T,UT> & u) {
	iunitstream(is, u.theX, u.theUnit);
	return is;
	}

	/**
	* OUnitErr is used to write out unitful numbers with an error
	* estimate on a standard ostream. using the helper function ouniterr
	* an energy <code>e</code> with an error estimate <code>de</code> can
	* be written out as eg. <code>cout << ouniterr(e, de,
	* GeV);</code>. The result will be presented in scientific format
	* (with the exponent divisible by three) with the relevant number of
	* significant digits with a single digit in parenthesis indicating
	* the error in the least significant digit,
	* eg. <code>1.23(2)e+03</code>.
	*/
	template <typename T, typename UT>
	struct OUnitErr {

	/** Constructor given an object to be written assuming the given
	* unit. */
	OUnitErr(const T & t, const T & dt, const UT & u): x(t/u), dx(dt/u) {}

	/** The number to be written. */
	double x;

	/** The estimated error of the number to be written. */
	double dx;

	};

	/** Helper function creating a OUnitErr object. */
	template <typename T, typename UT>
	inline OUnitErr<T,UT> ouniterr(const T & t, const T & dt, const UT & ut) {
	return OUnitErr<T,UT>(t, dt, ut);
	}

	/** Helper function creating a OUnitErr object. */
	inline OUnitErr<double,double> ouniterr(double t, double dt) {
	return OUnitErr<double,double>(t, dt, 1.0);
	}

	/** Output an OUnitErr object to a stream. */
	template <typename OStream, typename T, typename UT>
	OStream & operator<<(OStream & os, const OUnitErr<T,UT> & u) {
	if ( isnan(u.x) \|\| isinf(u.x) ) return os << u.x;
	- if ( isnan(u.dx) \|\| isinf(u.dx) ) return os << u.x << '(' << u.dx << ')';
	+ if ( isnan(u.dx) \|\| isinf(u.dx) ) {
	+ ostringstream out;
	+ out << u.x << '(' << u.dx << ')';
	+ return os << out.str();
	+ }
	double dx = min(u.dx, abs(u.x));
	if ( dx <= 0.0 ) return os << u.x;
	+ double x = abs(u.x);
	ostringstream osse;
	osse << std::scientific << setprecision(0) << dx;
	string sse = osse.str();
	string::size_type ee = sse.find('e');
	- long m = static_cast<long>(round(abs(u.x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
	+ long m = static_cast<long>(round(abs(x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
	int powx = m <= 0? os.precision(): int(log10(double(m)));
	if ( m <= 0 \|\| powx > os.precision() ) sse[0]='0';
	ostringstream oss;
	- oss << std::scientific << setprecision(powx) << u.x;
	+ oss << std::scientific << setprecision(powx) << x;
	string ss = oss.str();
	string::size_type e = ss.find('e');
	ostringstream out;
	int pp = std::atoi(ss.substr(e + 1).c_str());
	if ( pp%3 == 0 )
	out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
	else if ( (pp - 1)%3 == 0 ) {
	ostringstream oss;
	- oss << std::scientific << setprecision(powx) << u.x/10.0;
	+ oss << std::scientific << setprecision(powx) << x/10.0;
	string ss = oss.str();
	string::size_type e = ss.find('e');
	if ( powx == 0 )
	out << ss.substr(0, e) << "0(" << sse[0] << "0)" << ss.substr(e);
	else if ( powx == 1 )
	out << ss.substr(0, ss.find('.'))
	<< ss.substr(ss.find('.') + 1, e - ss.find('.') - 1)
	<< "(" << sse[0] << ")" << ss.substr(e);
	else {
	swap(ss[ss.find('.')], ss[ss.find('.') + 1]);
	out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
	}
	}
	else {
	ostringstream oss;
	- oss << std::scientific << setprecision(powx) << u.x*10.0;
	+ oss << std::scientific << setprecision(powx) << x*10.0;
	string ss = oss.str();
	string::size_type e = ss.find('e');
	if ( powx == 0 )
	out << "0." << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
	else {
	swap(ss[ss.find('.')], ss[ss.find('.') - 1]);
	out << ss.substr(0, ss.find('.')) << "0" << ss.substr(ss.find('.'), e)
	<< "(" << sse[0] << ")" << ss.substr(e);
	}
	}
	- return os << out.str();
	+ string res = out.str();
	+ if ( u.x < 0.0 )
	+ res = "-" + res;
	+ return os << res;
	}

	/**
	* The IUnitErr class is used to facilitate input of unitful numbers
	* with error estimates written out using the OUnitErr class.
	*
	*/
	template <typename T, typename UT>
	struct IUnitErr {

	/** Constructor given an object to be read assuming the given
	* unit. */
	IUnitErr(T & t, T & dt, const UT & u): x(t), dx(dt), ut(u) {}

	/** Reference to the object to be read. */
	T & x;

	/** The estimated error of the number to be read. */
	T & dx;

	/** The unit assumed when reading the object. */
	UT ut;

	};

	/** Helper function creating a IUnitErr object. */
	template <typename T, typename UT>
	inline IUnitErr<T,UT> iuniterr(T & t, T & dt, const UT & ut) {
	return IUnitErr<T,UT>(t, dt, ut);
	}

	/** Helper function creating a OUnitErr object. */
	inline IUnitErr<double,double> iuniterr(double & t, double & dt) {
	return IUnitErr<double,double>(t, dt, 1.0);
	}

	/** Input an IUnit object from a stream. */
	template <typename IStream, typename T, typename UT>
	IStream & operator>>(IStream & is, const IUnitErr<T,UT> & u) {
	string s;
	double x = 0.0;
	double dx = 0.0;
	double ex = 1.0;
	is >> s;
	string::size_type open = s.find('(');
	string::size_type close = s.find(')');
	string se = "0";
	string sp = "1";
	double pe = 1.0;
	if ( open != string::npos && close != string::npos ) {
	se = s.substr(open + 1);
	sp += s.substr(close + 1);
	string::size_type dot = s.find('.');
	if ( dot != string::npos && dot < open ) pe = std::pow(10.0, 1.0 - (open - dot));
	}

	istringstream(s) >> x;
	istringstream(se) >> dx;
	istringstream(sp) >> ex;

	u.x = xexu.ut;
	u.dx = dxexpe*u.ut;

	return is;
	}

	}

	#endif /* ThePEG_UnitIO_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:02 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805956
Default Alt Text: (35 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions