No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	diff --git a/Handlers/StdDependentXComb.cc b/Handlers/StdDependentXComb.cc
	--- a/Handlers/StdDependentXComb.cc
	+++ b/Handlers/StdDependentXComb.cc
	@@ -1,164 +1,165 @@
	// -- C++ --
	//
	// StdDependentXComb.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2010 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the StdDependentXComb class.
	//

	#include "StdDependentXComb.h"
	#include "StdXCombGroup.h"
	#include "ThePEG/MatrixElement/MEGroup.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace ThePEG;

	StdDependentXComb::StdDependentXComb(tStdXCombPtr newHead,
	const PBPair & newPartonBins, tMEPtr newME,
	const DiagramVector & newDiagrams)
	: StandardXComb(newHead->maxEnergy(),newHead->particles(),
	newHead->eventHandlerPtr(),
	const_ptr_cast<tSubHdlPtr>(newHead->subProcessHandler()),
	newHead->pExtractor(),newHead->CKKWHandler(),
	newPartonBins,newHead->cuts(),newME,newDiagrams,newHead->mirror(),
	newHead), resetIncoming(true) {}

	StdDependentXComb::StdDependentXComb()
	: StandardXComb(), resetIncoming(true) {}

	StdDependentXComb::~StdDependentXComb() { }

	tSubProPtr StdDependentXComb::construct() {
	// first get the meMomenta in their CMS, as this may
	// not be the case
	Boost cm = (meMomenta()[0] + meMomenta()[1]).findBoostToCM();
	if ( cm.mag2() > Constants::epsilon ) {
	for ( vector<Lorentz5Momentum>::iterator m = meMomenta().begin();
	m != meMomenta().end(); ++m ) {
	*m = m->boost(cm);
	}
	}
	tSubProPtr sub = StandardXComb::construct();
	sub->head(head()->subProcess());
	sub->groupWeight(lastCrossSection()/head()->lastCrossSection());
	return sub;
	}

	void StdDependentXComb::setProcess() {
	meMomenta().resize(mePartonData().size());
	}

	void StdDependentXComb::setPartonBinInstances(Energy2 scale) {

	PBIPair newBins;

	Direction<0> dir(true);
	newBins.first =
	new_ptr(PartonBinInstance(lastPartons().first,partonBins().first,scale));

	dir.reverse();
	newBins.second =
	new_ptr(PartonBinInstance(lastPartons().second,partonBins().second,scale));

	resetPartonBinInstances(newBins);
	setPartonBinInfo();

	lastPartons().first->scale(partonBinInstances().first->scale());
	lastPartons().second->scale(partonBinInstances().second->scale());

	}

	void StdDependentXComb::setIncomingPartons() {

	if ( !resetIncoming )
	return;

	resetIncoming = false;

	clean();
	createPartonBinInstances();
	setPartonBinInfo();

	lastParticles(head()->lastParticles());

	lastPartons(make_pair(mePartonData()[0]->produceParticle(Lorentz5Momentum()),
	mePartonData()[1]->produceParticle(Lorentz5Momentum())));

	Lorentz5Momentum pFirst = meMomenta()[0];
	Lorentz5Momentum pSecond = meMomenta()[1];

	if ( head()->matrixElement()->wantCMS() ) {
	Boost toLab = (head()->lastPartons().first->momentum() +
	head()->lastPartons().second->momentum()).boostVector();
	if ( toLab.mag2() > Constants::epsilon ) {
	pFirst.boost(toLab);
	pSecond.boost(toLab);
	}
	}

	lastPartons().first->set5Momentum(pFirst);
	lastPartons().second->set5Momentum(pSecond);

	lastS((lastParticles().first->momentum() +
	lastParticles().second->momentum()).m2());
	lastSHat((lastPartons().first->momentum() +
	lastPartons().second->momentum()).m2());
	lastP1P2(make_pair(0.0, 0.0));

	double x1 =
	lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus();
	double x2 =
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus();

	lastX1X2(make_pair(x1,x2));
	lastY(log(lastX1()/lastX2())*0.5);

	}

	CrossSection StdDependentXComb::dSigDR() {

	setProcess();
	setIncomingPartons();

	lastScale(head()->lastScale());
	lastAlphaS(head()->lastAlphaS());
	lastAlphaEM(head()->lastAlphaEM());

	- if ( !willPassCuts() ) {
	+ if ( (!willPassCuts() && !matrixElement()->headCuts()) \|\|
	+ !matrixElement()->apply() ) {
	subProcess(SubProPtr());
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastPDFWeight(head()->lastPDFWeight());

	CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();

	subProcess(SubProPtr());
	lastCrossSection(xsec);

	return xsec;

	}

	void StdDependentXComb::Init() {}

	void StdDependentXComb::persistentOutput(PersistentOStream &) const {
	}

	void StdDependentXComb::persistentInput(PersistentIStream &, int) {
	}

	ClassDescription<StdDependentXComb> StdDependentXComb::initStdDependentXComb;

	diff --git a/Handlers/StdXCombGroup.cc b/Handlers/StdXCombGroup.cc
	--- a/Handlers/StdXCombGroup.cc
	+++ b/Handlers/StdXCombGroup.cc
	@@ -1,366 +1,362 @@
	// -- C++ --
	//
	// StdXCombGroup.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2010 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the StdXCombGroup class.
	//

	#include "StdXCombGroup.h"
	#include "ThePEG/MatrixElement/MEGroup.h"
	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Handlers/SubProcessHandler.h"
	#include "ThePEG/Cuts/Cuts.h"
	#include "ThePEG/PDF/PartonExtractor.h"
	#include "ThePEG/Utilities/Debug.h"
	#include "ThePEG/Utilities/Maths.h"
	#include "ThePEG/PDT/ParticleData.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/SimplePhaseSpace.h"
	#include "ThePEG/Utilities/UtilityBase.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/SubProcessGroup.h"
	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/MatrixElement/MEBase.h"
	#include "ThePEG/MatrixElement/ColourLines.h"
	#include "ThePEG/Handlers/LuminosityFunction.h"
	#include "ThePEG/Handlers/CascadeHandler.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/EventRecord/TmpTransform.h"

	using namespace ThePEG;

	StdXCombGroup::StdXCombGroup(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME,
	const DiagramVector & newDiagrams, bool mir)
	: StandardXComb(newMaxEnergy,inc,newEventHandler,newSubProcessHandler,
	newExtractor, newCKKW, newPartonBins, newCuts,
	newME, newDiagrams, mir),
	theMEGroup(newME), theDependent(), theLastHeadCrossSection(ZERO) {}

	StdXCombGroup::StdXCombGroup()
	: StandardXComb(), theDependent() {}

	void StdXCombGroup::build(const PartonPairVec& allPBins) {
	for ( MEVector::const_iterator me = theMEGroup->dependent().begin();
	me != theMEGroup->dependent().end(); ++me ) {
	StdDependentXCombPtr dep =
	theMEGroup->makeDependentXComb(this,diagrams().front()->partons(),*me,allPBins);
	theDependent.push_back(dep);
	}
	}

	StdXCombGroup::~StdXCombGroup() { }

	CrossSection StdXCombGroup::dSigDR(const pair<double,double> ll, int nr, const double * r) {

	matrixElement()->flushCaches();

	if ( matrixElement()->keepRandomNumbers() ) {
	lastRandomNumbers().resize(nDim());
	copy(r,r+nDim(),lastRandomNumbers().begin());
	}

	pExtractor()->select(this);
	setPartonBinInfo();
	lastP1P2(ll);
	lastS(sqr(maxEnergy())/exp(lastP1() + lastP2()));

	meMomenta().resize(mePartonData().size());
	matrixElement()->setXComb(this);

	PPair partons;

	if ( !matrixElement()->haveX1X2() ) {

	if ( !pExtractor()->generateL(partonBinInstances(),
	r, r + nr - partonDims.second) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	partons = make_pair(partonBinInstances().first->parton(),
	partonBinInstances().second->parton());
	lastSHat(lastS()/exp(partonBinInstances().first->l() +
	partonBinInstances().second->l()));
	meMomenta()[0] = partons.first->momentum();
	meMomenta()[1] = partons.second->momentum();

	} else {
	if ( !matrixElement()->generateKinematics(r + partonDims.first) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
	matrixElement()->setKinematics();

	lastScale(matrixElement()->scale());

	partons.first = mePartonData()[0]->produceParticle(meMomenta()[0]);
	partons.second = mePartonData()[1]->produceParticle(meMomenta()[1]);

	Direction<0> dir(true);
	partonBinInstances().first =
	new_ptr(PartonBinInstance(lastParticles().first,partons.first,
	partonBins().first,lastScale()));
	dir.reverse();
	partonBinInstances().second =
	new_ptr(PartonBinInstance(lastParticles().second,partons.second,
	partonBins().second,lastScale()));
	}

	lastPartons(partons);

	if ( lastSHat() < cuts()->sHatMin() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastY(0.5*(partonBinInstances().second->l() -
	partonBinInstances().first->l()));
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() )
	SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());

	Energy summ = ZERO;
	if ( meMomenta().size() == 3 ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
	} else {
	for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass());
	summ += mePartonData()[i]->massMin();
	}
	if ( sqr(summ) >= lastSHat() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}

	matrixElement()->setXComb(this);
	if ( !matrixElement()->haveX1X2() )
	lastScale(max(lastSHat()/4.0, cuts()->scaleMin()));

	lastSHat(pExtractor()->generateSHat(lastS(), partonBinInstances(),
	r, r + nr - partonDims.second,
	matrixElement()->haveX1X2()));

	if ( !cuts()->sHat(lastSHat()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	r += partonDims.first;

	lastX1X2(make_pair(lastPartons().first->momentum().plus()/
	lastParticles().first->momentum().plus(),
	lastPartons().second->momentum().minus()/
	lastParticles().second->momentum().minus()));

	if ( !cuts()->x1(lastX1()) \|\| !cuts()->x2(lastX2()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	lastY((lastPartons().first->momentum() +
	lastPartons().second->momentum()).rapidity());
	if ( !cuts()->yHat(lastY()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	meMomenta()[0] = lastPartons().first->momentum();
	meMomenta()[1] = lastPartons().second->momentum();
	if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
	if ( matrixElement()->wantCMS() &&
	!matrixElement()->haveX1X2() )
	SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());

	if ( meMomenta().size() == 3 ) {
	if ( !matrixElement()->haveX1X2() )
	meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
	} else {
	if ( sqr(summ) >= lastSHat() ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}

	matrixElement()->setXComb(this);
	if ( !matrixElement()->haveX1X2() ) {
	if ( !matrixElement()->generateKinematics(r) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	}
	lastScale(matrixElement()->scale());
	if ( !cuts()->scale(lastScale()) ) {
	lastCrossSection(ZERO);
	return ZERO;
	}

	pair<bool,bool> evalPDFS =
	make_pair(matrixElement()->havePDFWeight1(),
	matrixElement()->havePDFWeight2());
	if ( mirror() )
	swap(evalPDFS.first,evalPDFS.second);
	lastPDFWeight(pExtractor()->fullFn(partonBinInstances(), lastScale(),
	evalPDFS));
	if ( lastPDFWeight() == 0.0 ) {
	lastCrossSection(ZERO);
	return ZERO;
	}
	matrixElement()->setKinematics();
	CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();

	bool noHeadPass = !willPassCuts() \|\| xsec == ZERO;
	if ( noHeadPass )
	lastCrossSection(ZERO);

	lastAlphaS(matrixElement()->alphaS());
	lastAlphaEM(matrixElement()->alphaEM());

	subProcess(SubProPtr());
	if ( CKKWHandler() && matrixElement()->maxMultCKKW() > 0 &&
	matrixElement()->maxMultCKKW() > matrixElement()->minMultCKKW() ) {
	newSubProcess(theMEGroup->subProcessGroups());
	CKKWHandler()->setXComb(this);
	xsec *= CKKWHandler()->reweightCKKW(matrixElement()->minMultCKKW(),
	matrixElement()->maxMultCKKW());
	}

	if ( matrixElement()->reweighted() ) {
	newSubProcess(theMEGroup->subProcessGroups());
	xsec = matrixElement()->reWeight() matrixElement()->preWeight();
	}

	lastHeadCrossSection(xsec);

	CrossSection depxsec = ZERO;

	if ( !theMEGroup->mcSumDependent() ) {
	for ( StdDepXCVector::const_iterator dep = theDependent.begin();
	dep != theDependent.end(); ++dep ) {
	if ( !(*dep) )
	continue;
	if ( !(**dep).matrixElement()->apply() )
	continue;
	if ( noHeadPass && (**dep).matrixElement()->headCuts() )
	continue;
	depxsec += (**dep).dSigDR();
	}
	} else {
	if ( !noHeadPass \|\| !theMEGroup->lastDependentXComb()->matrixElement()->headCuts() )
	depxsec = theDependent.size()*theMEGroup->lastDependentXComb()->dSigDR();
	}

	if ( xsec != ZERO ) {
	double rw = 1.0 + depxsec/xsec;
	xsec *= rw;
	} else {
	xsec = depxsec;
	}

	lastCrossSection(xsec);

	if ( xsec != ZERO )
	theMEGroup->lastEventStatistics();

	- if ( !theMEGroup->subProcessGroups() ) {
	- lastHeadCrossSection(xsec);
	- }
	-
	return xsec;

	}

	void StdXCombGroup::newSubProcess(bool) {

	// subprocess selection goes here
	// if me group returns an associated
	// selector

	StandardXComb::newSubProcess(theMEGroup->subProcessGroups());

	if ( !theMEGroup->subProcessGroups() )
	return;

	subProcess()->groupWeight(lastHeadCrossSection()/lastCrossSection());

	Ptr<SubProcessGroup>::tptr group =
	dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(subProcess());
	assert(group);

	for ( StdDepXCVector::iterator dep = theDependent.begin();
	dep != theDependent.end(); ++dep ) {
	if ( !(*dep) )
	continue;
	if ( (**dep).lastCrossSection() == ZERO \|\|
	!(**dep).matrixElement()->apply() )
	continue;
	tSubProPtr ds;
	try {
	ds = (**dep).construct();
	} catch(Veto&) {
	continue;
	}
	if ( ds )
	group->add(ds);
	}

	}

	tSubProPtr StdXCombGroup::construct() {

	if ( !theMEGroup->subProcessGroups() )
	return StandardXComb::construct();

	matrixElement()->setXComb(this);

	setPartonBinInfo();
	matrixElement()->setKinematics();

	newSubProcess(true);

	TmpTransform<tSubProPtr>
	tmp(subProcess(), Utilities::getBoostToCM(subProcess()->incoming()));
	if ( !cuts()->passCuts(*subProcess()) ) {
	subProcess()->groupWeight(0.0);
	}

	return subProcess();

	}


	void StdXCombGroup::Init() {}

	void StdXCombGroup::persistentOutput(PersistentOStream & os) const {
	os << theDependent << theMEGroup << ounit(theLastHeadCrossSection,nanobarn);
	}

	void StdXCombGroup::persistentInput(PersistentIStream & is, int) {
	is >> theDependent >> theMEGroup >> iunit(theLastHeadCrossSection,nanobarn);
	}

	ClassDescription<StdXCombGroup> StdXCombGroup::initStdXCombGroup;

	diff --git a/Handlers/StdXCombGroup.h b/Handlers/StdXCombGroup.h
	--- a/Handlers/StdXCombGroup.h
	+++ b/Handlers/StdXCombGroup.h
	@@ -1,187 +1,187 @@
	// -- C++ --
	//
	// StdXCombGroup.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2007 Leif Lonnblad
	// Copyright (C) 2009-2010 Simon Platzer
	//
	// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef ThePEG_StdXCombGroup_H
	#define ThePEG_StdXCombGroup_H
	// This is the declaration of the StdXCombGroup class.

	#include "StandardXComb.h"
	#include "StdDependentXComb.h"
	#include "StdXCombGroup.fh"
	#include "ThePEG/MatrixElement/MEGroup.fh"

	namespace ThePEG {

	/**
	* The StdXCombGroup class represents a 'head' XComb object
	* in association with a group of dependent XComb objects.
	*
	* @see StdDependendXComb
	* @see MEGroup
	*/
	class StdXCombGroup: public StandardXComb {

	/** MEBase needs to be a friend. */
	friend class MEBase;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Standard constructor.
	*/
	StdXCombGroup(Energy newMaxEnergy, const cPDPair & inc,
	tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
	tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
	const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME,
	const DiagramVector & newDiagrams, bool mir);

	/**
	* Default constructor.
	*/
	StdXCombGroup();

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

	public:

	/**
	* Generate a phase space point from a vector \a r of \a nr numbers
	* in the interval ]0,1[ and return the corresponding differential
	* cross section.
	*/
	virtual CrossSection dSigDR(const pair<double,double> ll, int nr, const double * r);

	/**
	+ * Return the cross section calculated from the head matrix element
	+ */
	+ CrossSection lastHeadCrossSection() const { return theLastHeadCrossSection; }
	+
	+ /**
	* Visit the dependent XComb objects
	*/
	const StdDepXCVector& dependent() const { return theDependent; }

	/**
	* Return the matrix element group steered by this
	* XComb group.
	*/
	tcMEGroupPtr meGroup() const { return theMEGroup; }

	/**
	* Initialize this XComb group
	*/
	void build(const PartonPairVec& allPBins);

	/**
	* Construct a sub-process object from the information available.
	*/
	virtual tSubProPtr construct();

	protected:

	/**
	* Construct the corresponding SubProcess object if it hasn't been
	* done before.
	*/
	virtual void newSubProcess(bool);

	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 interface.
	*/
	static void Init();

	protected:

	/**
	* Set the cross section calculated from the head matrix element
	*/
	void lastHeadCrossSection(CrossSection xs) { theLastHeadCrossSection = xs; }

	- /**
	- * Return the cross section calculated from the head matrix element
	- */
	- CrossSection lastHeadCrossSection() const { return theLastHeadCrossSection; }
	-
	private:

	/**
	* The MEGroup object
	*/
	MEGroupPtr theMEGroup;

	/**
	* The dependent XComb objects
	*/
	StdDepXCVector theDependent;

	/**
	* The cross section calculated from the head matrix element
	*/
	CrossSection theLastHeadCrossSection;

	private:

	/**
	* Describe a concrete class with persistent data.
	*/
	static ClassDescription<StdXCombGroup> initStdXCombGroup;

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

	};

	/** @cond TRAITSPECIALIZATIONS */

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

	/**
	* This template specialization informs ThePEG about the name of the
	* StdXCombGroup class.
	*/
	template <>
	struct ClassTraits<StdXCombGroup>:
	public ClassTraitsBase<StdXCombGroup> {
	/** Return the class name. */
	static string className() { return "ThePEG::StdXCombGroup"; }
	};

	/** @endcond */

	}

	#endif /* ThePEG_StdXCombGroup_H */
	diff --git a/MatrixElement/ColourLines.cc b/MatrixElement/ColourLines.cc
	--- a/MatrixElement/ColourLines.cc
	+++ b/MatrixElement/ColourLines.cc
	@@ -1,138 +1,143 @@
	// -- C++ --
	//
	// ColourLines.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 ColourLines class.
	//

	#include "ColourLines.h"
	#include "ColourLines.xh"
	#include "ThePEG/EventRecord/ColourLine.h"
	#include "ThePEG/EventRecord/MultiColour.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Utilities/StringUtils.h"

	using namespace ThePEG;

	ColourLines::ColourLines(string s) {
	+ reset(s);
	+}
	+
	+void ColourLines::reset(string s) {
	+ theLines.clear();
	while ( true ) {
	string line = StringUtils::car(s, ",");
	line = StringUtils::stripws(line);
	Line l;
	while (line!="") {
	string loc = StringUtils::car(line);
	string first = StringUtils::car(loc,":");
	string second = StringUtils::cdr(loc,":");
	if(second!="") {
	int i;
	istringstream is(first);
	is >> i;
	int j;
	istringstream is2(second);
	is2 >> j;
	l.push_back(make_pair(i,j));
	}
	else {
	int i;
	istringstream is(first);
	is >> i;
	l.push_back(make_pair(i,0));
	}
	line = StringUtils::cdr(line);
	};
	if ( l.empty() ) return;
	theLines.push_back(l);
	s = StringUtils::cdr(s, ",");
	}
	}

	void ColourLines::connect(const tPVector & partons) const {
	VertexVector sinks;
	VertexVector sources;
	long np = partons.size();

	// Create each line and connect the specified partons to them. Save
	// all lines coming from a source or ending in a sink.
	for ( LineVector::size_type il = 0; il < theLines.size(); ++il ) {
	const Line & line = theLines[il];
	ColinePtr cline = new_ptr(ColourLine());
	for ( Line::size_type i = 0; i < line.size(); ++i ) {
	if ( line[i].first > np ) {
	// this is a colour source.
	int is = line[i].first - np;
	sources.resize(is);
	sources[is - 1].push_back(cline);
	} else if ( -line[i].first > np ) {
	// this is a colour sink.
	int is = -line[i].first - np;
	sources.resize(is);
	sources[is - 1].push_back(cline);
	} else if ( line[i].first > 0 ) {
	// This is a coloured particle.
	if ( !partons[line[i].first - 1]->hasColour() )
	throw ColourGeometryException(partons, line);
	if(line[i].second==0) {
	cline->addColoured(partons[line[i].first - 1]);
	}
	else {
	Ptr<MultiColour>::pointer colour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(partons[line[i].first - 1]->colourInfo());
	assert(colour);
	colour->colourLine(cline,line[i].second);
	}
	} else {
	if ( !partons[-line[i].first - 1]->hasAntiColour() )
	throw ColourGeometryException(partons, line);
	if(line[i].second==0) {
	cline->addAntiColoured(partons[-line[i].first - 1]);
	}
	else {
	Ptr<MultiColour>::pointer colour =
	dynamic_ptr_cast<Ptr<MultiColour>::pointer>(partons[-line[i].first - 1]->colourInfo());
	assert(colour);
	colour->antiColourLine(cline,line[i].second);
	}
	}
	}
	}

	// Now connect up all lines steming from sources.
	for ( VertexVector::size_type i = 0; i < sources.size(); ++i ) {
	if ( sources[i].empty() ) continue;
	if ( sources[i].size() != 3 ) throw ColourGeometryException(partons,
	vector<pair<int,int> >() );
	sources[i][0]->setSourceNeighbours(sources[i][1], sources[i][2]);
	}

	// Now connect up all lines ending in sinks.
	for ( VertexVector::size_type i = 0; i < sinks.size(); ++i ) {
	if ( sinks[i].empty() ) continue;
	if ( sinks[i].size() != 3 ) throw ColourGeometryException(partons,
	vector<pair<int,int> >());
	sinks[i][0]->setSinkNeighbours(sinks[i][1], sinks[i][2]);
	}
	}

	ColourGeometryException::
	ColourGeometryException(const tPVector & p, const vector<pair<int,int> > & c) {
	if ( c.empty() )
	theMessage << "The number of colour lines steming from one colour source "
	<< "or ending in one colour sink was not equal to three.\n";
	else {
	theMessage << "Cannot connect the following partons:\n";
	for ( unsigned i = 0; i < p.size(); ++i )
	theMessage << " " << p[i]->PDGName();
	theMessage << "\n to the following colour line:\n";
	for ( unsigned i = 0; i < c.size(); ++i ) theMessage << " (" << c[i].first << ","
	<< c[i].second << ") ";
	theMessage << endl;
	}
	severity(maybeabort);
	}

	diff --git a/MatrixElement/ColourLines.h b/MatrixElement/ColourLines.h
	--- a/MatrixElement/ColourLines.h
	+++ b/MatrixElement/ColourLines.h
	@@ -1,78 +1,89 @@
	// -- C++ --
	//
	// ColourLines.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_ColourLines_H
	#define ThePEG_ColourLines_H
	// This is the declaration of the ColourLines class.

	#include "ThePEG/Config/ThePEG.h"

	namespace ThePEG {

	/**
	* The ColourLines class defines the colour flow in a SubProcess. It
	* defines a number of colour lines and specifies which particles are
	* connected to them.
	*
	*/
	class ColourLines: public Base {

	public:

	/** A single colour line */
	typedef vector<pair<int,int> > Line;
	/** A vector of colour lines. */
	typedef vector<Line> LineVector;
	/** A vector of <code>ColourLine</code>. */
	typedef vector<ColinePtr> Vertex;
	/** A vector of vertices. */
	typedef vector<Vertex> VertexVector;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	ColourLines() {}

	/**
	* The standard constructor. The string \a s should contain a
	* comma-separated sequence of integers. Each sequence of numbers
	* indicates a colour line and the integer represents a parton
	* connected to it. If the integer is negative, it means that the
	* line is the corresponding partons anti-colour. Note that the
	* partons are numbered from 1: The first spacelike particle is 1, the second
	* is 2 and the internal time-like and outgoing are numbered after all
	* the spacelike particles.
	*/
	ColourLines(string s);
	//@}

	+ /**
	+ * Reset this ColourLines object. The string \a s should contain a
	+ * comma-separated sequence of integers. Each sequence of numbers
	+ * indicates a colour line and the integer represents a parton
	+ * connected to it. If the integer is negative, it means that the
	+ * line is the corresponding partons anti-colour. Note that the
	+ * partons are numbered from 1: The first incoming is 1, the second
	+ * is 2 and the internal and outgoing are numbered 3 and upwards.
	+ */
	+ void reset(string s);
	+
	public:

	/**
	* Create the corresponding <code>ColourLine</code>s and connect the
	* given \a partons. The partons are assumed to be in the same order
	* as the numbers specified in the constructor.
	*/
	void connect(const tPVector & partons) const;

	private:

	/**
	* The vector of colour lines.
	*/
	LineVector theLines;

	};

	}

	#endif /* ThePEG_ColourLines_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:27 PM (1 d, 10 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805827
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions