No OneTemporary
Actions

Size

59 KB

Subscribers

None

View Options

	diff --git a/EventRecord/Particle.cc b/EventRecord/Particle.cc
	--- a/EventRecord/Particle.cc
	+++ b/EventRecord/Particle.cc
	@@ -1,524 +1,525 @@
	// -- C++ --
	//
	// Particle.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2017 Leif Lonnblad
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined member functions of
	// the Particle class.
	//
	#include "Particle.h"
	#include "ThePEG/EventRecord/Step.h"
	#include "ThePEG/EventRecord/Event.h"
	#include "ThePEG/EventRecord/ColourLine.h"
	#include "ThePEG/Utilities/Rebinder.h"
	#include "ThePEG/Config/algorithm.h"
	#include "ThePEG/EventRecord/ParticleTraits.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include <iostream>
	#include <iomanip>
	#include <cctype>

	#ifdef ThePEG_TEMPLATES_IN_CC_FILE
	#include "Particle.tcc"
	#endif

	using namespace ThePEG;

	Particle::ParticleRep::ParticleRep(const ParticleRep & p)
	: theParents(p.theParents), theChildren(p.theChildren),
	thePrevious(p.thePrevious), theNext(p.theNext),
	theBirthStep(p.theBirthStep), theVertex(p.theVertex),
	theLifeLength(p.theLifeLength), theScale(p.theScale),
	theVetoScale(p.theVetoScale), theNumber(p.theNumber),
	theExtraInfo(p.theExtraInfo.size()) {
	if ( p.theColourInfo )
	theColourInfo = dynamic_ptr_cast<CBPtr>(p.theColourInfo->clone());
	if ( p.theSpinInfo )
	theSpinInfo = dynamic_ptr_cast<SpinPtr>(p.theSpinInfo->clone());
	for ( int i = 0, N = p.theExtraInfo.size(); i < N; ++i )
	theExtraInfo[i] = p.theExtraInfo[i]->clone();
	}

	Particle::Particle(const Particle & p)
	- : Base(p), theData(p.theData), theMomentum(p.theMomentum), theRep(p.theRep) {
	+ : Base(p), theData(p.theData), theMomentum(p.theMomentum),
	+ theRep(p.theRep), theStatus(p.theStatus) {
	if ( p.theRep ) {
	theRep = new ParticleRep(*p.theRep);
	theRep->theParents.clear();
	}
	}

	Particle::~Particle() {
	if ( theRep ) {
	if ( colourLine() ) colourLine()->removeColoured(this);
	if ( antiColourLine() ) antiColourLine()->removeAntiColoured(this);
	delete theRep;
	}
	theRep = 0;
	theData = cEventPDPtr();
	}

	void Particle::initFull() {
	if ( theRep ) return;
	theRep = new ParticleRep;

	Energy width = data().generateWidth(mass());

	if ( width > ZERO ) {
	Time lifetime = data().generateLifeTime(mass(), width);
	theRep->theLifeLength.setTau(lifetime);
	theRep->theLifeLength.
	setVect((momentum().vect()*(lifetime /
	max(mass(), Constants::epsilon*GeV))));
	theRep->theLifeLength.rescaleEnergy();
	}
	}

	PPtr Particle::clone() const {
	return ptr_new<PPtr>(*this);
	}

	void Particle::rebind(const EventTranslationMap & trans) {
	for ( ParticleVector::iterator pit = rep().theChildren.begin();
	pit != rep().theChildren.end(); ++pit ) pit = trans.translate(pit);
	for ( tParticleVector::iterator pit = rep().theParents.begin();
	pit != rep().theParents.end(); ++pit ) pit = trans.translate(pit);
	rep().thePrevious = trans.translate(rep().thePrevious);
	rep().theNext = trans.translate(rep().theNext);
	if ( hasColourInfo() ) colourInfo()->rebind(trans);
	if ( spinInfo() ) spinInfo()->rebind(trans);
	rep().theBirthStep = trans.translate(rep().theBirthStep);
	for ( EIVector::const_iterator ie = rep().theExtraInfo.begin();
	ie != rep().theExtraInfo.end(); ++ie ) (**ie).rebind(trans);
	}

	tParticleSet Particle::siblings() const {
	tParticleSet theSiblings;
	for ( tParticleVector::const_iterator pit = parents().begin();
	pit != parents().end(); ++pit )
	theSiblings.insert((pit)->children().begin(), (pit)->children().end());
	theSiblings.erase(const_cast<Particle *>(this));
	return theSiblings;
	}

	void Particle::colourNeighbour(tPPtr p, bool anti) {
	tColinePtr line = colourLine(!anti);
	if ( !line ) line = ColourLine::create(this, !anti);
	line->addColoured(p, anti);
	}

	void Particle::outgoingColour(tPPtr p, bool anti) {
	tColinePtr line = colourLine(anti);
	if ( !line ) line = ColourLine::create(this, anti);
	line->addColoured(p, anti);
	}

	tPPtr Particle::incomingColour(bool anti) const {
	if ( !hasColourInfo() ) return tPPtr();
	tColinePtr line = colourLine(anti);
	if ( !line ) return tPPtr();
	for ( int i = 0, N = parents().size(); i < N; ++i )
	if ( parents()[i]->hasColourLine(line, anti) ) return parents()[i];
	return tPPtr();
	}

	tPPtr Particle::outgoingColour(bool anti) const {
	if ( !hasColourInfo() ) return tPPtr();
	tColinePtr line = colourLine(anti);
	if ( !line ) return tPPtr();
	for ( int i = 0, N = children().size(); i < N; ++i )
	if ( children()[i]->hasColourLine(line, anti) ) return children()[i];
	return tPPtr();
	}

	LorentzPoint Particle::labVertex() const {
	LorentzPoint r(rep().theBirthStep && rep().theBirthStep->collision()?
	vertex() + rep().theBirthStep->collision()->vertex():
	vertex());
	return r;
	}

	void Particle::setLabVertex(const LorentzPoint & p) {
	rep().theVertex = ( rep().theBirthStep && rep().theBirthStep->collision()?
	p - rep().theBirthStep->collision()->vertex() : p );
	}

	void Particle::transform(const LorentzRotation & r) {
	if ( hasRep() && spinInfo() ) spinInfo()->transform(momentum(), r);
	theMomentum.transform(r);
	if ( !hasRep() ) return;
	rep().theVertex.transform(r);
	rep().theLifeLength.transform(r);
	}

	void Particle::deepTransform(const LorentzRotation & r) {
	transform(r);
	if ( !theRep ) return;
	for ( int i = 0, N = children().size(); i < N; ++i )
	rep().theChildren[i]->deepTransform(r);
	if ( rep().theNext ) rep().theNext->deepTransform(r);
	}

	void Particle::rotateX(double a) {
	LorentzRotation r;
	r.rotateX(a);
	transform(r);
	}

	void Particle::deepRotateX(double a) {
	LorentzRotation r;
	r.rotateX(a);
	deepTransform(r);
	}

	void Particle::rotateY(double a) {
	LorentzRotation r;
	r.rotateY(a);
	transform(r);
	}

	void Particle::deepRotateY(double a) {
	LorentzRotation r;
	r.rotateY(a);
	deepTransform(r);
	}

	void Particle::rotateZ(double a) {
	LorentzRotation r;
	r.rotateZ(a);
	transform(r);
	}

	void Particle::deepRotateZ(double a) {
	LorentzRotation r;
	r.rotateZ(a);
	deepTransform(r);
	}

	void Particle::rotate(double a, const Axis & axis) {
	LorentzRotation r;
	r.rotate(a, axis);
	transform(r);
	}

	void Particle::deepRotate(double a, const Axis & axis) {
	LorentzRotation r;
	r.rotate(a, axis);
	deepTransform(r);
	}

	string Particle::outputFormat =
	"%n3%s10 %i7 %p[,]0 %c(,) %^^0%vv0 %>>0%<>0 %l{,}0\n"
	" %x10.3%y10.3%z10.3%e10.3%m10.3\n";

	int getNumber(string::const_iterator & pos, int def) {
	if ( !isdigit(*pos) ) return def;
	def = *pos++ - '0';
	while ( isdigit(pos) ) def = 10def + *pos++ - '0';
	return def;
	}

	void writePrecision(ostream & os, string::const_iterator & pos,
	int defw, int defp, double x) {
	defw = getNumber(pos, defw);
	if ( *pos == '.' ) defp = getNumber(++pos, defp);
	int oldp = os.precision();
	os << setprecision(defp) << setw(defw) << x << setprecision(oldp);
	}

	void writeStringAdjusted(ostream & os, bool left, int w, string str) {
	while ( !left && w-- > int(str.size()) ) os << ' ';
	os << str;
	while ( left && w-- > int(str.size()) ) os << ' ';
	}

	template <typename Container>
	void writeParticleRanges(ostream & os, const Container & co, char sep, int w) {
	set<int> cnum;
	for ( typename Container::const_iterator it = co.begin();
	it != co.end(); ++it) cnum.insert((**it).number());

	bool elipsis = false;
	int last = -10;
	for ( set<int>::iterator it = cnum.begin(); it != cnum.end(); ++it) {
	int n = *it;
	int next = 0;
	set<int>::iterator itn = it;
	if ( ++itn != cnum.end() ) next = *itn;
	bool writeit = true;
	bool writesep = false;
	if ( elipsis && ( n != last + 1 \|\| n != next - 1 ) )
	elipsis = false;
	else if ( !elipsis && n == last + 1 && n == next -1 ) {
	os << "..";
	elipsis = true;
	writeit = false;
	}
	else if ( elipsis && n == last + 1 && n == next -1 )
	writeit = false;
	else if ( it != cnum.begin() )
	writesep = true;
	if ( writeit ) {
	if ( writesep ) os << sep;
	os << setw(w) << n;
	}
	last = n;
	}
	}

	ostream & ThePEG::operator<<(ostream & os, const Particle & p) {
	return p.print(os, p.birthStep());
	}

	ostream & Particle::print(ostream & os, tcStepPtr step) const {
	if ( !step ) step = birthStep();
	tCollPtr coll = step? step->collision(): tCollPtr();
	tEventPtr event = coll? coll->event(): tEventPtr();
	string::const_iterator pos = Particle::outputFormat.begin();
	ios::fmtflags saveflags = os.setf(ios::fixed, ios::floatfield);
	while ( pos != Particle::outputFormat.end() ) {
	if ( *pos == '%' && ++pos != Particle::outputFormat.end() ) {
	bool left = false;
	if ( *pos == '-' ) {
	left = true;
	os.setf(ios::left, ios::adjustfield);
	++pos;
	} else {
	os.setf(ios::right, ios::adjustfield);
	}
	char mark;
	char open;
	char close;
	char sep;
	int w;
	string str;
	string fill;
	if ( pos == Particle::outputFormat.end() ) break;
	bool fullColour = false;
	switch ( *pos ) {
	case 'n':
	os << setw(getNumber(++pos, 3)) << number();
	break;
	case 'i':
	os << setw(getNumber(++pos, 8)) << id();
	break;
	case 's':
	writeStringAdjusted(os, left, getNumber(++pos, 8), PDGName());
	break;
	case 'x':
	writePrecision(os, ++pos, 10, 3, momentum().x()/GeV);
	break;
	case 'y':
	writePrecision(os, ++pos, 10, 3, momentum().y()/GeV);
	break;
	case 'z':
	writePrecision(os, ++pos, 10, 3, momentum().z()/GeV);
	break;
	case 'e':
	writePrecision(os, ++pos, 10, 3, momentum().e()/GeV);
	break;
	case 'm':
	writePrecision(os, ++pos, 10, 3, momentum().mass()/GeV);
	break;
	case 'P':
	fullColour = true;
	[[fallthrough]];
	case 'p':
	open = *++pos;
	sep = *++pos;
	close = *++pos;
	w = getNumber(++pos, 0);
	if ( parents().empty() ) break;
	if ( open ) os << open;
	writeParticleRanges(os, parents(), sep, w);
	if ( fullColour && hasColourInfo() &&
	( incomingColour() \|\| incomingAntiColour() ) ) {
	if ( close ) os << open;
	if ( incomingColour() )
	os << "+" << incomingColour()->number();
	if ( incomingAntiColour() )
	os << "-" << incomingAntiColour()->number();
	if ( close ) os << close;
	}
	if ( close ) os << close;
	break;
	case 'l':
	open = *++pos;
	sep = *++pos;
	close = *++pos;
	w = getNumber(++pos, 0);
	if ( hasColourInfo() &&
	( colourLine() \|\| antiColourLine() ) && event) {
	if ( open ) os << open;
	vector<tcColinePtr> clines = colourInfo()->colourLines();
	for ( int i = 0, N = clines.size(); i < N; ++i ) {
	if ( i > 0 && sep ) os << sep;
	clines[i]->write(os, event, false);
	}
	vector<tcColinePtr> aclines = colourInfo()->antiColourLines();
	for ( int i = 0, N = aclines.size(); i < N; ++i ) {
	if ( ( i > 0 \|\| clines.size() ) && sep ) os << sep;
	aclines[i]->write(os, event, true);
	}
	if ( close ) os << close;
	}
	break;
	case 'C':
	fullColour = true;
	[[fallthrough]];
	case 'c':
	open = *++pos;
	sep = *++pos;
	close = *++pos;
	w = getNumber(++pos, 0);
	if ( children().empty() ) break;
	if ( open ) os << open;
	writeParticleRanges(os, children(), sep, w);
	if ( fullColour && hasColourInfo() &&
	( outgoingColour() \|\| outgoingAntiColour() ) ) {
	if ( close ) os << open;
	if ( outgoingColour() )
	os << "+" << outgoingColour()->number();
	if ( outgoingAntiColour() )
	os << "-" << outgoingAntiColour()->number();
	if ( close ) os << close;
	}
	if ( close ) os << close;
	break;
	case '>':
	mark = *++pos;
	w = getNumber(++pos, 0);
	if ( hasColourInfo() && step && step->colourNeighbour(this) ) {
	os << setw(w-1) << step->colourNeighbour(this)->number() << mark;
	}
	break;
	case '<':
	mark = *++pos;
	w = getNumber(++pos, 0);
	if ( hasColourInfo() && step && step->antiColourNeighbour(this) ) {
	int n = step->antiColourNeighbour(this)->number();
	ostringstream oss;
	oss << mark << n;
	writeStringAdjusted(os, left, w, oss.str());
	}
	break;
	case 'v':
	mark = *++pos;
	w = getNumber(++pos, 0);
	if ( next() ) {
	if ( left && mark ) os << mark;
	os << setw(w) << next()->number();
	if ( !left && mark ) os << mark;
	}
	break;
	case '^':
	mark = *++pos;
	w = getNumber(++pos, 0);
	if ( previous() ) {
	if ( left && mark ) os << mark;
	os << setw(w) << previous()->number();
	if ( !left && mark ) os << mark;
	}
	break;
	case 'd':
	switch ( *++pos ) {
	case 'x':
	writePrecision(os, ++pos, 10, 3, lifeLength().x()/mm);
	break;
	case 'y':
	writePrecision(os, ++pos, 10, 3, lifeLength().y()/mm);
	break;
	case 'z':
	writePrecision(os, ++pos, 10, 3, lifeLength().z()/mm);
	break;
	case 't':
	writePrecision(os, ++pos, 10, 3, lifeLength().e()/mm);
	break;
	case 'T':
	writePrecision(os, ++pos, 10, 3, lifeLength().tau()/mm);
	break;
	}
	break;
	case 'V':
	switch ( *++pos ) {
	case 'x':
	writePrecision(os, ++pos, 10, 3, vertex().x()/mm);
	break;
	case 'y':
	writePrecision(os, ++pos, 10, 3, vertex().y()/mm);
	break;
	case 'z':
	writePrecision(os, ++pos, 10, 3, vertex().z()/mm);
	break;
	case 't':
	writePrecision(os, ++pos, 10, 3, vertex().e()/mm);
	break;
	}
	break;
	case 'L':
	switch ( *++pos ) {
	case 'x':
	writePrecision(os, ++pos, 10, 3, labVertex().x()/mm);
	break;
	case 'y':
	writePrecision(os, ++pos, 10, 3, labVertex().y()/mm);
	break;
	case 'z':
	writePrecision(os, ++pos, 10, 3, labVertex().z()/mm);
	break;
	case 't':
	writePrecision(os, ++pos, 10, 3, labVertex().e()/mm);
	break;
	}
	break;
	default:
	os << *pos++;
	}
	} else {
	if ( pos != Particle::outputFormat.end() ) os << *pos++;
	}
	}
	os.flags(saveflags);
	return os;
	}

	void Particle::debugme() const {
	cerr << *this;
	EventRecordBase::debugme();
	}

	void Particle::persistentOutput(PersistentOStream & os) const {
	EventConfig::putParticleData(os, theData);
	- os << ounit(theMomentum, GeV) << bool( theRep != 0 );
	+ os << ounit(theMomentum, GeV) << theStatus << bool( theRep != 0 );
	if ( !theRep ) return;
	os << rep().theParents << rep().theChildren
	<< rep().thePrevious << rep().theNext << rep().theBirthStep
	<< ounit(rep().theVertex, mm) << ounit(rep().theLifeLength, mm)
	<< ounit(rep().theScale, GeV2) << ounit(rep().theVetoScale, GeV2)
	<< rep().theNumber << rep().theDecayMode
	<< rep().theColourInfo << rep().theSpinInfo << rep().theExtraInfo;
	}

	void Particle::persistentInput(PersistentIStream & is, int) {
	bool readRep;
	EventConfig::getParticleData(is, theData);
	- is >> iunit(theMomentum, GeV) >> readRep;
	+ is >> iunit(theMomentum, GeV) >> theStatus >> readRep;
	if ( !readRep ) return;
	if ( !hasRep() ) theRep = new ParticleRep;

	is >> rep().theParents >> rep().theChildren
	>> rep().thePrevious >> rep().theNext >> rep().theBirthStep
	>> iunit(rep().theVertex, mm) >> iunit(rep().theLifeLength, mm)
	>> iunit(rep().theScale, GeV2) >> iunit(rep().theVetoScale, GeV2)
	>> rep().theNumber >> rep().theDecayMode
	>> rep().theColourInfo >> rep().theSpinInfo >> rep().theExtraInfo;
	}

	ClassDescription<Particle> Particle::initParticle;

	void Particle::Init() {}
	diff --git a/EventRecord/Particle.h b/EventRecord/Particle.h
	--- a/EventRecord/Particle.h
	+++ b/EventRecord/Particle.h
	@@ -1,1180 +1,1195 @@
	// -- C++ --
	//
	// Particle.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2017 Leif Lonnblad
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef ThePEG_Particle_H
	#define ThePEG_Particle_H
	// This is the decalaration of the Particle class.

	#include "EventConfig.h"
	#include "ThePEG/Vectors/Lorentz5Vector.h"
	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/Utilities/ClassDescription.h"
	#include "ThePEG/EventRecord/MultiColour.h"
	#include "ThePEG/EventRecord/SpinInfo.h"
	#include "ThePEG/PDT/ParticleData.h"

	namespace ThePEG {

	/**
	* The Particle class is used to describe an instance of a
	* particle. Properties of the corresponding particle type can be
	* accessed through a pointer to a ParticleData object.
	*
	* A Particle object contains pointers to other particles, such as a
	* list of parents and a list of children. It may also contain a
	* pointer to the previous or next instance of the same physical
	* particle if the properties of a given particle has been changed
	* during the generation. Coloured particles contains pointers to
	* ColourLine defining the colour connections to other particles.
	*
	* The Particle also has a pointer to the Step object where it was
	* first introduced in the Event.
	*
	* When printing a particle the format of the output is governed by
	* the static <code>outputFormat</code> string. When a particle is sent
	* to an <code>ostream</code>, the format string is written but with
	* keys prefixed by the <code>\%</code> character replaced with
	* infromation about the particle as follows:<BR> <code>\%\%</code> is
	* replaced by a singel <code>\%</code><BR> <code>\%C</code> sets a flag so
	* that subsequent output of children and parents etc. will contain
	* colour information.<BR> <code>\%n</code> is replaced by the particles
	* number in a fied ow width<BR> <code>\%s</code> is replaced by the name
	* of the particle<BR> <code>\%i</code> is replaced by the id number of
	* the particle type<BR><code>\%x, \%y, \%z, \%e, \%m</code> is replaced by
	* the x-, y-, z-, energy- and mass-component of the particles
	* momentum respectively<BR><code>\%dx, \%dy, \%dz, \%dt, \%dT</code> is
	* replaced by the x-, y-, z-, time- and invariant time-component of
	* the particles lifeLength respectively<BR><code>\%Vx, \%Vy, \%Vz,
	* \%Vt</code> is replaced by the x-, y-, z- and time-component of the
	* creation point relative to the vertex of the
	* collision.<BR><code>\%Lx, \%Ly, \%Lz, \%Lt</code> is replaced by the x-,
	* y-, z- and time-component of the creation point in the current
	* lab-system.<BR> <code>\%p[,]</code> is replaced by a list (of numbers)
	* of the particles parents enclosed by <code>[</code> and <code>]</code>
	* and separated by <code>,</code>, The parent from which the particle
	* inherits its (anti-) colour is marked by a (-)+<BR>
	* <code>\%c(,)</code> is replaced by a list (of numbers) of the
	* particles children enclosed by <code>(</code> and <code>)</code> and
	* separated by <code>,</code>, The child which inherits the particles
	* (anti-) colour is marked by a (-)+<BR> <code>\%></code> is replaced
	* by the number of the colour neighbor<BR> <code>\%<</code> is
	* replaced by the number of the anti-colour neighbor<BR>
	* <code>\%^</code> is replaced by the number of the previous instance of
	* the same physical particle<BR> <code>\%v</code> is replaced by the
	* number of the next instance of the same physical particle.<BR>
	* <code>\%l{,}</code> is replaced by the indices of the colour lines to
	* which this particle is connected enclosed by <code>{</code> and
	* <code>}</code> and separated by <code>,</code>, The line corresponding
	* to the (anti-) colour of the particle is prefixed by a (-)+
	*
	* @see Event
	* @see Collision
	* @see Step
	* @see SubProcess
	* @see Lorentz5Vector
	* @see ColourLine
	* @see ColourBase
	*/
	class Particle: public EventRecordBase {

	public:

	/** Most of the Event classes are friends with each other. */
	friend class Event;
	/** Most of the Event classes are friends with each other. */
	friend class Collision;
	/** Most of the Event classes are friends with each other. */
	friend class Step;
	/** Most of the Event classes are friends with each other. */
	friend class SubProcess;
	/** ParticleData needs to be a friend. */
	friend class ParticleData;

	struct ParticleRep;

	public:

	/**
	* Standard Constructor. Note that the default constructor is
	* private - there is no particle without a pointer to a
	* ParticleData object.
	*/
	- Particle(tcEventPDPtr newData) : theData(newData), theRep(0) {}
	+ Particle(tcEventPDPtr newData) : theData(newData), theRep(0), theStatus(0) {}

	/**
	* Copy constructor.
	*/
	Particle(const Particle &);

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

	/** @name Functions relating to ancestry of particles. */
	//@{
	/**
	* Returns true if and only if this particle has decayed.
	*/
	bool decayed() const {
	return hasRep() && !rep().theChildren.empty();
	}

	/**
	* The list of decay products.
	*/
	const ParticleVector & children() const {
	static const ParticleVector null;
	return hasRep() ? rep().theChildren : null;
	}

	/**
	* Add a child (the childs parent pointer will be set accordingly).
	*/
	void addChild(tPPtr c) {
	rep().theChildren.push_back(c);
	(c->rep()).theParents.push_back(this);
	}

	/**
	* Remove the given child from the list of children of this particle
	* (the corresponding parent pointer of the child will also be
	* removed).
	*/
	void abandonChild(tPPtr child) {
	removeChild(child);
	child->removeParent(this);
	}

	/**
	* The list of parent particles.
	*/
	const tParticleVector & parents() const {
	static const tParticleVector null;
	return hasRep() ? rep().theParents : null;
	}

	/**
	* Return a set of neighboring particles coming from the same decay
	* as this one. The return value is a newly recalculated set
	* every time. It must be stored to be used further, do not directly call
	* e.g. siblings().begin() or siblings().end()!
	*/
	tParticleSet siblings() const;

	/**
	* Undo the decay of this particle, removing all children (and
	* grand children ...) from the event record
	*/
	void undecay() {
	if ( hasRep() ) {
	rep().theChildren.clear();
	rep().theNext = tPPtr();
	}
	}

	/**
	* If this particle has decayed set the corresponding decay mode.
	*/
	void decayMode(tDMPtr dm) { rep().theDecayMode = dm; }

	/**
	* If this particle has decayed get the corresponding decay mode.
	*/
	tDMPtr decayMode() const {
	return hasRep() ? rep().theDecayMode : tDMPtr();
	}

	/**
	* Next instance. Pointer to another instance of the same
	* physical particle in later steps.
	*/
	tPPtr next() const {
	return hasRep() ? rep().theNext : PPtr();
	}

	/**
	* Previous instance. Pointer to another instance of the same
	* physical particle in earlier steps.
	*/
	tPPtr previous() const {
	return hasRep() ? rep().thePrevious : tPPtr();
	}

	/**
	* Original instance. If there exists another previous instance of
	* this particle return this instance (recursively).
	*/
	tcPPtr original() const {
	return previous() ? tcPPtr(previous()->original()) : tcPPtr(this);
	}

	/**
	* Original instance. If there exists another previous instance of
	* this particle return this instance (recursively).
	*/
	tPPtr original() {
	return previous() ? previous()->original() : tPPtr(this);
	}


	/**
	* Final instance. If there exists another subsequent instance of
	* this particle return this instance (recursively).
	*/
	tcPPtr final() const {
	return next() ? tcPPtr(next()->final()) : tcPPtr(this);
	}

	/**
	* Final instance. If there exists another subsequent instance of
	* this particle return this instance (recursively).
	*/
	tPPtr final() {
	return next() ? next()->final() : tPPtr(this);
	}

	//@}

	/** @name Relations to the Event and Step. */
	//@{
	/**
	* Get the first Step object where this particle occurred.
	*/
	tStepPtr birthStep() const {
	return hasRep() ? rep().theBirthStep : tStepPtr();
	}

	/**
	* Get the order-number for this particle in the current event.
	*/
	int number() const {
	return hasRep() ? rep().theNumber : 0;
	}
	+
	+ /**
	+ * Get the status code of the particle
	+ */
	+ int status() const { return theStatus; }
	+
	+ /**
	+ * Set the status code of the particle
	+ */
	+ void status(int n) { theStatus = n; }
	//@}

	/** @name Access the underlying ParticleData object. */
	//@{
	/**
	* Access the ParticleData object of this particle type
	*/
	const ParticleDataClass & data() const { return *theData; }

	/**
	* Access the ParticleData object of this particle type
	*/
	tcEventPDPtr dataPtr() const { return theData; }

	/**
	* Return the PDG name of this particle.
	*/
	const string & PDGName() const { return data().PDGName(); }

	/**
	* Return the PDG id number of this particle.
	*/
	long id() const { return data().id(); }
	//@}

	/** @name Functions to access the momentum. */
	//@{
	/**
	* Return the momentum of this particle.
	*/
	const Lorentz5Momentum & momentum() const { return theMomentum; }

	/**
	* Set the 3-momentum of this particle. The energy is set to be
	* consistent with the mass.
	*/
	void set3Momentum(const Momentum3 & p) {
	theMomentum.setVect(p);
	theMomentum.rescaleEnergy();
	}

	/**
	* Set the momentum of this particle. Afterwards, the underlying
	* Lorentz5Momentum may have inconsistent mass.
	*/
	void setMomentum(const LorentzMomentum & p) {
	theMomentum = p;
	}

	/**
	* Set the momentum and mass.
	*/
	void set5Momentum(const Lorentz5Momentum & p) {
	theMomentum = p;
	}

	/**
	* Acces the mass of this particle.
	*/
	Energy mass() const { return momentum().mass(); }

	/**
	* Acces the mass of this particle type.
	*/
	Energy nominalMass() const { return data().mass(); }

	/**
	* Get the scale at which this particle is considered resolved.
	*/
	Energy2 scale() const {
	return hasRep() ? rep().theScale : -1.0*GeV2;
	}

	/**
	* Set the scale at which this particle is considered resolved.
	*/
	void scale(Energy2 q2) { rep().theScale = q2; }

	/**
	* Get the scale above which this particle should
	* not radiate.
	*/
	Energy2 vetoScale() const {
	return hasRep() ? rep().theVetoScale : -1.0*GeV2;
	}

	/**
	* Set the scale above which this particle should
	* not radiate.
	*/
	void vetoScale(Energy2 q2) { rep().theVetoScale = q2; }

	/**
	* Return the transverse mass (squared), calculated from the energy
	* and the longitudinal momentum.
	*/
	Energy2 mt2() const { return sqr(momentum().t()) - sqr(momentum().z()); }

	/**
	* Return the transverse mass (squared), calculated from the energy
	* and the longitudinal momentum.
	*/
	Energy mt() const { return sqrt(mt2()); }

	/**
	* Return the transverse mass (squared), calculated from the mass
	* and the transverse momentum.
	*/
	Energy2 perpmass2() const { return momentum().perp2() + momentum().mass2(); }

	/**
	* Return the transverse mass (squared), calculated from the mass
	* and the transverse momentum.
	*/
	Energy perpmass() const { return sqrt(perpmass2()); }

	/**
	* Return the (pseudo) rapidity.
	*/
	double rapidity() const {
	return ( Pplus() > ZERO && Pminus() > ZERO )?
	0.5*log(Pplus()/Pminus()) : Constants::MaxFloat;
	}

	/**
	* Return the (pseudo) rapidity.
	*/
	double eta() const {
	Energy rho = momentum().rho();
	return rho > abs(momentum().z())?
	0.5*log((rho+momentum().z())/(rho-momentum().z())) : Constants::MaxFloat;
	}

	/**
	* Return the positive and negative light-cone momenta.
	*/
	Energy Pplus() const { return momentum().plus(); }
	/**
	* Return the positive and negative light-cone momenta.
	*/
	Energy Pminus() const { return momentum().minus(); }
	//@}

	/** @name Functions to access the position. */
	//@{
	/**
	* The creation vertex of this particle. The point is given
	* relative to the collision vertex.
	*/
	const LorentzPoint & vertex() const {
	static const LorentzPoint null;
	return hasRep() ? rep().theVertex : null;
	}

	/**
	* The creation vertex of this particle. The absolute
	* position in the lab is given.
	*/
	LorentzPoint labVertex() const;

	/**
	* The decay vertex of this particle. The point is given
	* relative to the collision vertex.
	*/
	LorentzPoint decayVertex() const {
	return vertex() + lifeLength();
	}

	/**
	* The decay vertex of this particle. The absolute
	* position in the lab is given.
	*/
	LorentzPoint labDecayVertex() const {
	return labVertex() + lifeLength();
	}

	/**
	* The life time/length. Return the Lorentz vector connecting the
	* creation to the decay vertes.
	*/
	const Lorentz5Distance & lifeLength() const {
	static const Lorentz5Distance null;
	return hasRep() ? rep().theLifeLength : null;
	}

	/**
	* Set the creation vertex relative to the collision vertex.
	*/
	void setVertex(const LorentzPoint & p) {
	rep().theVertex = p;
	}

	/**
	* Set the creation vertex in the lab frame of this particle.
	*/
	void setLabVertex(const LorentzPoint &);

	/**
	* Set the life length of this particle. The life time will be
	* automatically rescaled to be consistent with the invariant
	* distance.
	*/
	void setLifeLength(const Distance & d) {
	rep().theLifeLength.setVect(d);
	rep().theLifeLength.rescaleEnergy();
	}

	/**
	* Set the life time/length of a particle. The invariant distance
	* may become inconsistent.
	*/
	void setLifeLength(const LorentzDistance & d) {
	rep().theLifeLength = d;
	}

	/**
	* Set the life time/length of a particle.
	*/
	void setLifeLength(const Lorentz5Distance & d) {
	rep().theLifeLength = d;
	}

	/**
	* The invariant life time of this particle.
	*/
	Time lifeTime() const { return lifeLength().m(); }

	//@}

	/** @name Functions for (Lorentz) transformations. */
	//@{
	/**
	* Do Lorentz transformations on this particle.
	*/
	void transform(const LorentzRotation & r);

	/**
	* Do Lorentz transformations on this particle. \a bx, \a by and \a
	* bz are the boost vector components.
	*/
	void boost(double bx, double by, double bz) {
	transform(LorentzRotation(Boost(bx, by, bz)));
	}

	/**
	* Do Lorentz transformations on this particle. \a b is the boost
	* vector.
	*/
	void boost(const Boost & b) { transform(LorentzRotation(b)); }

	/**
	* Rotate around the x-axis.
	*/
	void rotateX(double a);

	/**
	* Rotate around the y-axis.
	*/
	void rotateY(double a);

	/**
	* Rotate around the z-axis.
	*/
	void rotateZ(double a);

	/**
	* Rotate around the given \a axis.
	*/
	void rotate(double a, const Axis & axis);

	/**
	* Mirror in the xy-plane.
	*/
	void mirror() { theMomentum.setZ(-theMomentum.z()); }

	/**
	* Do Lorentz transformations on this particle and its decendants.
	*/
	void deepTransform(const LorentzRotation & r);

	/**
	* Do Lorentz transformations on this particle and its
	* decendants. \a bx, \a by and \a bz are the boost vector
	* components.
	*/
	void deepBoost(double bx, double by, double bz) {
	deepTransform(LorentzRotation(Boost(bx, by, bz)));
	}

	/**
	* Do Lorentz transformations on this particle and its
	* decendants. \a b is the boost vector.
	*/
	void deepBoost(const Boost & b) { deepTransform(LorentzRotation(b)); }

	/**
	* Rotate this particle and its decendants around the x-axis.
	*/
	void deepRotateX(double a);

	/**
	* Rotate this particle and its decendants around the y-axis.
	*/
	void deepRotateY(double a);

	/**
	* Rotate this particle and its decendants around the z-axis.
	*/
	void deepRotateZ(double a);

	/**
	* Rotate this particle and its decendants around the given \a axis.
	*/
	void deepRotate(double a, const Axis & axis);

	//@}

	/** @name Functions controlling possible mass/momentum inconsistencies. */
	//@{
	/**
	* Return the relative inconsistency in the mass component.
	*/
	double massError() const { return theMomentum.massError(); }

	/**
	* Return the relative inconsistency in the energy component.
	*/
	double energyError() const { return theMomentum.energyError(); }

	/**
	* Return the relative inconsistency in the spatial components.
	*/
	double rhoError() const { return theMomentum.rhoError(); }

	/**
	* Rescale energy, so that the invariant length/mass of the
	* LorentzVector agrees with the current one.
	*/
	void rescaleEnergy() { theMomentum.rescaleEnergy(); }

	/**
	* Rescale spatial component, so that the invariant length/mass of
	* the LorentzVector agrees with the current one.
	*/
	void rescaleRho() { theMomentum.rescaleRho(); }

	/**
	* Set the invariant length/mass member, so that it agrees with the
	* invariant length/mass of the LorentzVector.
	*/
	void rescaleMass() { theMomentum.rescaleMass(); }
	//@}

	/** @name Acces incormation about colour connections */
	//@{
	/**
	* True if this particle has colour information. To determine if
	* this particle is actually coloured, the coloured(), hasColour() or
	* hasAntiColour() methods should be used instead.
	*/
	bool hasColourInfo() const {
	return hasRep() && rep().theColourInfo;
	}

	/**
	* Return the colour lines to which this particles anti-colour is
	* connected.
	*/
	tColinePtr antiColourLine() const {
	return hasColourInfo() ? colourInfo()->antiColourLine() : tColinePtr();
	}

	/**
	* Return the colour lines to which this particles (\a anti-)colour
	* is connected.
	*/
	tColinePtr colourLine(bool anti = false) const {
	if ( anti ) return antiColourLine();
	return hasColourInfo() ? colourInfo()->colourLine() : tColinePtr();
	}

	/**
	* Return true if the particle is connected to the given (\a anti-)
	* colour \a line.
	*/
	bool hasColourLine(tcColinePtr line, bool anti = false) const {
	return hasColourInfo() ? colourInfo()->hasColourLine(line, anti) : false;
	}

	/**
	* Return true if the particle is connected to the given anti-colour
	* \a line.
	*/
	bool hasAntiColourLine(tcColinePtr line) const {
	return hasColourLine(line, true);
	}

	/**
	* True if this particle type is not a colour singlet.
	*/
	bool coloured() const { return data().coloured(); }

	/**
	* True if this particle type carries (\a anti-)colour.
	*/
	bool hasColour(bool anti = false) const { return data().hasColour(anti); }

	/**
	* True if this particle type carries anti-colour.
	*/
	bool hasAntiColour() const { return data().hasAntiColour(); }

	/**
	* Get the ColourBase object.
	*/
	tcCBPtr colourInfo() const {
	return hasRep() ? rep().theColourInfo : CBPtr();
	}

	/**
	* Get the ColourBase object.
	*/
	tCBPtr colourInfo() {
	if ( !rep().theColourInfo ) {
	switch(theData->iColour()) {
	case PDT::Colour6:
	case PDT::Colour6bar:
	rep().theColourInfo = new_ptr(MultiColour());
	break;
	default:
	rep().theColourInfo = new_ptr(ColourBase());
	}
	}
	return rep().theColourInfo;
	}

	/**
	* Set the ColourBase object.
	*/
	void colourInfo(tCBPtr c) {
	rep().theColourInfo = c;
	}

	/**
	* Get a pointer to the colour neighbor. Returns a particle in the
	* range \a first to \a last which colour is connected to the same
	* line as this particles anti-colour. If \a anti is true return
	* antiColourNeighbour().
	*/
	template <typename Iterator>
	typename std::iterator_traits<Iterator>::value_type
	colourNeighbour(Iterator first, Iterator last, bool anti = false) const;

	/**
	* Get a pointer to the anti-colour neighbor. Returns a particle in
	* the range \a first to \a last which anti-colour is
	* connected to the same line as this particles colour.
	*/
	template <typename Iterator>
	typename std::iterator_traits<Iterator>::value_type
	antiColourNeighbour(Iterator first, Iterator last) const {
	return colourNeighbour(first, last, true);
	}

	/**
	* Set the colour neighbor. Connects the given particles colour to
	* the same colour line as this particles anti-colour. If \a anti is
	* true call antiColourNeighbour(tPPtr).
	*/
	void colourNeighbour(tPPtr, bool anti = false);

	/**
	* Set the anti-colour neighbor. Connects the given particles
	* anti-colour to the same colour line as this particles colour.
	*/
	void antiColourNeighbour(tPPtr p) { colourNeighbour(p, true); }

	/**
	* Connect colour. Create a colour line connecting to it this
	* particles colour and the given particles anti-colour.
	*/
	void antiColourConnect(tPPtr neighbour) {
	colourConnect(neighbour, true);
	}

	/**
	* Connect colour. Create a colour line connecting to it this
	* particles anti-colour and the given particles colour. If \a anti
	* is true call antiColourConnect(tPPtr).
	*/
	void colourConnect(tPPtr neighbour, bool anti = false) {
	colourNeighbour(neighbour, anti);
	}

	/**
	* Incoming colour. Return the parent particle which colour is
	* connected to the same colour line as this particle. If \a anti is
	* true return incomingAntiColour().
	*/
	tPPtr incomingColour(bool anti = false) const;

	/**
	* Incoming anti-colour. Return the parent particle which
	* anti-colour is connected to the same colour line as this
	* particle.
	*/
	tPPtr incomingAntiColour() const { return incomingColour(true); }

	/**
	* Set incoming colour. Connect this particles colour to the same
	* colour line as the given particle. If \a anti
	* is true call incomingAntiColour(tPPtr).
	*/
	void incomingColour(tPPtr p, bool anti = false) { p->outgoingColour(this, anti); }

	/**
	* Set incoming anti-colour. Connect this particles anti colour to
	* the same colour line as the given particle.
	*/
	void incomingAntiColour(tPPtr p) { p->outgoingColour(this, true); }

	/**
	* Outgoing colour. Return the daughter particle which colour is
	* connected to the same colour line as this particle. If \a anti is
	* true return outgoingAntiColour().
	*/
	tPPtr outgoingColour(bool anti = false) const;
	/**
	* Outgoing anti-colour. Return the daughter particle which
	* anti-colour is connected to the same colour line as this
	* particle.
	*/
	tPPtr outgoingAntiColour() const { return outgoingColour(true); }

	/**
	* Set outgoing colour. Connect this particles colour to the same
	* colour line as the given particle. If \a anti
	* is true call outgoingAntiColour(tPPtr).
	*/
	void outgoingColour(tPPtr, bool anti = false);

	/**
	* Set outgoing anti-colour. Connect this particles anti-colour to
	* the same colour line as the given particle.
	*/
	void outgoingAntiColour(tPPtr p) { outgoingColour(p, true); }

	/**
	* Specify colour flow. Calls outgoingColour(tPPtr,bool).
	*/
	void colourFlow(tPPtr child, bool anti = false) {
	outgoingColour(child, anti);
	}

	/**
	* Specify anticolour flow. Calls outgoingAntiColour(tPPtr,bool).
	*/
	void antiColourFlow(tPPtr child) { colourFlow(child, true); }

	/**
	* Remove all colour information;
	*/
	void resetColour() {
	if ( hasColourInfo() ) rep().theColourInfo = CBPtr();
	}

	//@}

	/** @name Functions to access spin. */
	//@{
	/**
	* Return the Spin object.
	*/
	tcSpinPtr spinInfo() const {
	return hasRep() ? rep().theSpinInfo : SpinPtr();
	}

	/**
	* Return the Spin object.
	*/
	tSpinPtr spinInfo() {
	return hasRep() ? rep().theSpinInfo : SpinPtr();
	}

	/**
	* Set the Spin object.
	*/
	void spinInfo(tSpinPtr s) { rep().theSpinInfo = s; }
	//@}

	/** @name Accessing user-defined information. */
	//@{
	/**
	* Access user-defined information as a vector of EventInfoBase pointers.
	*/
	const EIVector & getInfo() const {
	static const EIVector null;
	return hasRep() ? rep().theExtraInfo : null;
	}

	/**
	* Access user-defined information as a vector of EventInfoBase pointers.
	*/
	EIVector & getInfo() { return rep().theExtraInfo; }
	//@}

	public:

	/** @name Accessing user-defined information. */
	//@{
	/**
	* True if this particle has instantiated the object with
	* information other than type and momentum.
	*/
	bool hasRep() const { return theRep; }

	/**
	* If this particle has only a type and momentum, instantiate the
	* rest of the information.
	*/
	void initFull();

	//@}

	public:

	/** @name Input and output functions. */
	//@{
	/**
	* Standard function for writing to a persistent stream.
	*/
	void persistentOutput(PersistentOStream &) const;

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

	//@}

	/**
	* Print particle info to a stream \a os. The \a step is used to
	* access information about colour neighbors and other struff.
	*/
	ostream & print(ostream & os, tcStepPtr step = tcStepPtr()) const;

	/**
	* Print a range of particles.
	*/
	template <typename Iterator>
	static void PrintParticles(ostream & os, Iterator first, Iterator last,
	tcStepPtr step = tcStepPtr());

	/**
	* Print a container of particles.
	*/
	template <typename Cont>
	static inline void PrintParticles(ostream & os, const Cont & c,
	tcStepPtr step = tcStepPtr()) {
	PrintParticles(os, c.begin(), c.end(), step);
	}

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

	/**
	* Specify how to print particles. The format string is analogous to
	* the one used by eg. the unix 'date' command as described above.
	*/
	static string outputFormat;

	private:

	/**
	* Standard clone function.
	*/
	virtual PPtr clone() const;

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

	/**
	* Set the order-number for this particle in the current event.
	*/
	void number(int n) { rep().theNumber = n; }

	/**
	* Remove the given particle from the list of children.
	*/
	void removeChild(tPPtr c) {
	if ( hasRep() )
	rep().theChildren.erase(remove(rep().theChildren.begin(),
	rep().theChildren.end(), c),
	rep().theChildren.end());
	}

	/**
	* Remove the given particle from the list of parents.
	*/
	void removeParent(tPPtr p) {
	if ( hasRep() )
	rep().theParents.erase(remove(rep().theParents.begin(),
	rep().theParents.end(), p),
	rep().theParents.end());
	}

	/**
	* Set the mass of this particle.
	*/
	void mass(Energy m) { theMomentum.setMass(m); }

	/**
	* Set the invaiant life time of this particle.
	*/
	void lifeTime(Length t) { rep().theLifeLength.setTau(t); }

	/**
	* Return a reference to the bulk information of this particle. if
	* no ParticleRep object exists, one is created.
	*/
	ParticleRep & rep() {
	if ( !hasRep() ) initFull();
	return *theRep;
	}

	/**
	* Return a reference to the bulk information of this particle. if
	* no ParticleRep object exists, we return the default values.
	*/
	const ParticleRep & rep() const {
	static const ParticleRep null;
	return hasRep() ? *theRep : null;
	}

	/**
	* The pointer to the ParticleData object
	*/
	cEventPDPtr theData;

	/**
	* The momentum.
	*/
	Lorentz5Momentum theMomentum;

	/**
	* The rest of the information in this particle is only instantiated
	* if needed.
	*/
	ParticleRep * theRep;

	+ /**
	+ * The status code of the particle
	+ */
	+ int theStatus;
	+
	public:

	/**
	* This class is used internally in the Particle class to represent
	* information besides momentum and type. A corresponding object
	* will only be instantiated if needed to save memory and time when
	* temporarily creating particles.
	*/
	struct ParticleRep {

	/**
	* Default constructor.
	*/
	ParticleRep() : theScale(-1.0GeV2), theVetoScale(-1.0GeV2), theNumber(0) {}

	/**
	* Copy constructor.
	*/
	ParticleRep(const ParticleRep &);

	/**
	* The pointers to the parents.
	*/
	tParticleVector theParents;

	/**
	* The pointers to the children.
	*/
	ParticleVector theChildren;

	/**
	* The pointer to the previous instance.
	*/
	tPPtr thePrevious;

	/**
	* The pointer to the next instance.
	*/
	PPtr theNext;

	/**
	* If this particle has decayed this is the pointer to the
	* corresponding decay mode.
	*/
	tDMPtr theDecayMode;

	/**
	* The pointer to the first step where this particle occurred.
	*/
	tStepPtr theBirthStep;

	/**
	* The creation point.
	*/
	LorentzPoint theVertex;

	/**
	* The life time/length.
	*/
	Lorentz5Distance theLifeLength;

	/**
	* the resolution scale.
	*/
	Energy2 theScale;

	/**
	* the veto scale.
	*/
	Energy2 theVetoScale;

	/**
	* The order-number for this particle in the current event.
	*/
	int theNumber;

	/**
	* A pointer to the colour information object.
	*/
	CBPtr theColourInfo;

	/**
	* Spin information
	*/
	SpinPtr theSpinInfo;

	/**
	* Additional used-defined information.
	*/
	EIVector theExtraInfo;

	};

	public:

	/**
	* Print out debugging information for this object on std::cerr. To
	* be called from within a debugger via the debug() function.
	*/
	virtual void debugme() const;

	protected:

	/**
	* Private default constructor must only be used by the
	* PersistentIStream class via the ClassTraits<Particle> class.
	*/
	- Particle() : theRep(0) {}
	+ Particle() : theRep(0), theStatus(0) {}

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

	private:

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

	/**
	* Describe concrete class with persistent data.
	*/
	static ClassDescription<Particle> initParticle;

	};

	/**
	* Write a Particle object to a stream.
	*/
	ostream & operator<<(ostream &, const Particle &);


	/** @cond TRAITSPECIALIZATIONS */

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

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

	/** @endcond */

	}

	#ifndef ThePEG_TEMPLATES_IN_CC_FILE
	#include "Particle.tcc"
	#endif

	#endif /* ThePEG_Particle_H */
	diff --git a/Vectors/HepMCConverter.tcc b/Vectors/HepMCConverter.tcc
	--- a/Vectors/HepMCConverter.tcc
	+++ b/Vectors/HepMCConverter.tcc
	@@ -1,337 +1,337 @@
	// -- C++ --
	//
	// HepMCConverter.tcc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2017 Leif Lonnblad
	//
	// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HepMCConverter class.
	//

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


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

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

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

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

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

	init(ev, nocopies);

	}

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

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

	init(ev, nocopies);

	}

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

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

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

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

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

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

	// Create GenParticle's and map them to the ThePEG particles.
	for ( int i = 0, N = all.size(); i < N; ++i ) {
	tcPPtr p = all[i];
	if ( nocopies && p->next() ) continue;
	if ( pmap.find(p) != pmap.end() ) continue;
	pmap[p] = createParticle(p);
	if ( !p->children().empty() \|\| p->next() ) {
	// If the particle has children it should have a decay vertex:
	vertices.push_back(Vertex());
	decv[p] = &vertices.back();
	vertices.back().in.insert(p);
	}

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

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

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

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

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

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

	// and the PDF info
	setPdfInfo(ev);

	// and the cross section info
	Traits::setCrossSection(*geneve,
	eh->integratedXSec()/picobarn,
	eh->integratedXSecErr()/picobarn);
	for ( int i = 0, N = all.size(); i < N; ++i ) {
	tcPPtr p = all[i];
	if ( pmap.find(p) == pmap.end() ) continue;
	GenParticlePtrT gp = pmap[p];
	if ( p->hasColourInfo() ) {
	// Check if the particle is connected to colour lines, in which
	// case the lines are mapped to an integer and set in the
	// GenParticle's Flow info.
	tcColinePtr l;
	if ( (l = p->colourLine()) ) {
	if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500;
	Traits::setColourLine(*gp, 1, flowmap[l]);
	}
	if ( (l = p->antiColourLine()) ) {
	if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500;
	Traits::setColourLine(*gp, 2, flowmap[l]);
	}
	}
	if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) {
	DPair pol = p->spinInfo()->polarization();
	Traits::setPolarization(*gp, pol.first, pol.second);
	}
	}
	}

	template <typename HepMCEventT, typename Traits>
	typename HepMCConverter<HepMCEventT,Traits>::GenParticlePtrT
	HepMCConverter<HepMCEventT,Traits>::createParticle(tcPPtr p) const {
	int status = 1;
	size_t nChildren = p->children().size();
	if ( nChildren > 0 \|\| p->next() ) status = 11;
	if ( nChildren > 1 ) {
	long id = p->data().id();
	if ( BaryonMatcher::Check(id) \|\| MesonMatcher::Check(id) \|\|
	id == ParticleID::muminus \|\| id == ParticleID::muplus \|\|
	id == ParticleID::tauminus \|\| id == ParticleID::tauplus ) {
	bool child = false;
	for(unsigned int ix=0;ix<nChildren;++ix) {
	if(p->children()[ix]->id()==id) {
	child = true;
	break;
	}
	}
	if ( !child ) {
	if(p->data().widthCut()!=ZERO) {
	if(p->mass() <= p->data().massMax() &&
	p->mass() >= p->data().massMin() )
	status = 2;
	}
	else {
	status = 2;
	}
	}
	}
	}
	GenParticlePtrT gp =
	- Traits::newParticle(p->momentum(), p->id(), status, energyUnit);
	+ Traits::newParticle(p->momentum(), p->id(), p->status() ? p->status() : status, energyUnit);

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

	return gp;

	}

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

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

	GenVertexPtrT gv = Traits::newVertex();

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

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

	return gv;
	}

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

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

	}



	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:43 AM (1 d, 8 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4983123
Default Alt Text: (59 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions