No OneTemporary
Actions

Size

26 KB

Subscribers

None

View Options

	diff --git a/EventRecord/SpinInfo.cc b/EventRecord/SpinInfo.cc
	--- a/EventRecord/SpinInfo.cc
	+++ b/EventRecord/SpinInfo.cc
	@@ -1,182 +1,187 @@
	// -- C++ --
	//
	// SpinInfo.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 2003-2017 Peter Richardson, 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 SpinInfo class.
	//
	// Author: Peter Richardson
	//
	#include "SpinInfo.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "Event.h"

	using namespace ThePEG;

	const double SpinInfo::_eps=1.0e-8;

	SpinInfo::SpinInfo(const SpinInfo & x)
	: EventInfoBase(x), _production(x._production), _decay(x._decay),
	_timelike(x._timelike),
	_prodloc(x._prodloc), _decayloc(x._decayloc),
	_decayed(x._decayed), _developed(x._developed),_rhomatrix(x._rhomatrix),
	_Dmatrix(x._Dmatrix),_spin(x._spin),
	_productionmomentum(x._productionmomentum),
	_decaymomentum(x._decaymomentum),
	_currentmomentum(x._currentmomentum)
	{
	x._production=VertexPtr();
	x._decay=VertexPtr();
	// set the vertex so it now points to the copy
	if(_production) {
	// for timelike
	if(_timelike) _production->resetOutgoing(this,_prodloc);
	// for spacelike
	else _production->resetIncoming(this,_prodloc);
	}
	}

	EIPtr SpinInfo::clone() const {
	tcSpinPtr temp=this;
	return const_ptr_cast<SpinPtr>(temp);
	}

	void SpinInfo::rebind(const EventTranslationMap & trans) {
	if(_production) _production = trans.translate(_production);
	if(_decay) _decay = trans.translate(_decay);
	EventInfoBase::rebind(trans);
	}


	NoPIOClassDescription<SpinInfo> SpinInfo::initSpinInfo;
	// Definition of the static class description member.

	void SpinInfo::Init() {}

	void SpinInfo::update() const {
	// number of instances fo this object
	int nref=referenceCount();
	if(nref<2\|\|nref>3) return;
	// work out the number of references there should be
	int nmin=0;
	// check the production pointers
	if(_production) {
	if(_timelike) { if(_production->outgoing()[_prodloc]==this) ++nmin;}
	else { if(_production->incoming()[_prodloc]==this) ++nmin;}
	}
	// check the decay pointers
	if(_decay) {
	if(_decay->incoming()[_decayloc]==this) ++nmin;
	}
	// delete the pointers
	SpinPtr temp;
	if(nmin+1==nref) {
	// delete the production pointers
	if(_production) {
	if(_timelike) {
	if(_production->outgoing()[_prodloc]==this)
	_production->resetOutgoing(SpinPtr(),_prodloc);
	}
	else {
	if(_production->incoming()[_prodloc]==this)
	_production->resetIncoming(SpinPtr(),_prodloc);
	}
	}
	// delete the decay pointers
	if(_decay) {
	if(_decay->incoming()[_decayloc]==this)
	_decay->resetIncoming(SpinPtr(),_decayloc);
	}
	}
	}

	void SpinInfo::decay(bool recursive) const {
	// if the particle has already been decayed do nothing
	if(_decayed) return;
	// otherwise we need to obtain the correct rho (timelike) or D (spacelike) matrix
	assert(_developed!=NeedsUpdate);
	if(_timelike) {
	if(_developed==Developed&&iSpin()!=PDT::Spin0) _developed=NeedsUpdate;
	if(productionVertex()) {
	if(recursive) redecay();
	else _rhomatrix = productionVertex()->getRhoMatrix(_prodloc,true);
	}
	}
	else {
	if(_developed==Developed&&iSpin()!=PDT::Spin0) _developed=NeedsUpdate;
	if(_production) _Dmatrix = _production->getDMatrix(_prodloc);
	}
	_decaymomentum = _currentmomentum;
	_decayed=true;
	}

	void SpinInfo::redevelop() const {
	assert(developed()==NeedsUpdate);
	// calculate rho/D matrix
	if(_timelike) {
	_Dmatrix = decayVertex() ?
	decayVertex()->getDMatrix(decayLocation()) : RhoDMatrix(iSpin());
	}
	else {
	_rhomatrix = decayVertex() ?
	decayVertex()->getRhoMatrix(decayLocation(),false) : RhoDMatrix(iSpin());
	}
	// update the D matrix of this spininfo
	_developed = Developed;
	// update the parent if needed
	if(productionVertex() &&
	productionVertex()->incoming().size()==1) {
	tcSpinPtr parent = _timelike ?
	productionVertex()->incoming()[0] : productionVertex()->outgoing()[0];
	- parent->needsUpdate();
	- parent->redevelop();
	+ if ( parent->developed() != StopUpdate ) {
	+ parent->needsUpdate();
	+ parent->redevelop();
	+ }
	}
	}

	void SpinInfo::develop() const {
	// if the particle has already been developed do nothing
	switch(_developed) {
	case Developed:
	return;
	+ case StopUpdate:
	+ return;
	case NeedsUpdate:
	redevelop();
	return;
	case Undeveloped:
	if(_timelike) {
	if(_decay) _Dmatrix = _decay->getDMatrix(_decayloc);
	else _Dmatrix = RhoDMatrix(iSpin());
	}
	else {
	if(_decay) _rhomatrix = _decay->getRhoMatrix(_decayloc,false);
	else _rhomatrix = RhoDMatrix(iSpin());
	}
	_developed=Developed;
	return;
	}
	}

	void SpinInfo::redecay() const {
	if(!productionVertex()) return;
	if(productionVertex()->incoming().size()==1) {
	tcSpinPtr parent;
	if(productionVertex()->incoming()[0]->timelike())
	parent = productionVertex()->incoming()[0];
	else {
	if(productionVertex()->outgoing()[0]!=this)
	parent = productionVertex()->outgoing()[0];
	else
	parent = productionVertex()->outgoing()[1];
	}
	- parent->redecay();
	+ if ( parent->developed() != StopUpdate )
	+ parent->redecay();
	}
	if(timelike())
	_rhomatrix = productionVertex()->getRhoMatrix(_prodloc,true);
	else
	_Dmatrix = productionVertex()->getDMatrix(_prodloc);
	}
	diff --git a/EventRecord/SpinInfo.h b/EventRecord/SpinInfo.h
	--- a/EventRecord/SpinInfo.h
	+++ b/EventRecord/SpinInfo.h
	@@ -1,450 +1,457 @@
	// -- C++ --
	//
	// SpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 2003-2017 Peter Richardson, 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_SpinInfo_H
	#define ThePEG_SpinInfo_H
	// This is the declaration of the SpinInfo class.

	#include "ThePEG/EventRecord/EventInfoBase.h"
	#include "ThePEG/PDT/PDT.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "HelicityVertex.h"

	namespace ThePEG {

	/**
	* The SpinInfo is the base class for the spin information for the
	* spin correlation algorithm. The implementations for different spin
	* states inherit from this.
	*
	* The class contains pointers to the vertex where the particle is
	* produced and where it decays, together with methods to set/get
	* these.
	*
	* There are two flags decayed which store information on the state
	* of the particle.
	*
	* The decayed() members provides access to the _decay data member
	* which is true if the spin density matrix required to perform the
	* decay of a timelike particle has been calculated (this would be a
	* decay matrix for a spacelike particle.) This is set by the
	* decay() method which calls a method from the production vertex to
	* calculate this matrix. The decay() method should be called by a
	* decayer which uses spin correlation method before it uses the
	* spin density matrix to calculate the matrix element for the
	* decay.
	*
	* The developed() member provides access to the _developed data
	* member which is true if the decay matrix required to perform the
	* decays of the siblings of a particle has been calculated (this
	* would a spin density matrix for a spacelike particle.) This is
	* set by the developed() method which calls a method from the decay
	* vertex to calculate the matrix. The developed() method is called
	* by a DecayHandler which is capable of performing spin
	* correlations after all the unstable particles produced by a
	* decaying particle are decayed.
	*
	* Methods are also provided to access the spin density and decay
	* matrices for a particle.
	*
	* @author Peter Richardson
	*
	*/
	class SpinInfo: public EventInfoBase {

	public:

	/**
	* Status for the implementation of spin correlations
	*/
	enum DevelopedStatus {
	Undeveloped=0, /*< Not developed. /
	Developed=1, /*< Developed. /
	- NeedsUpdate=2 /*< Developed but needs recalculating due to some change. /
	+ NeedsUpdate=2, /*< Developed but needs recalculating due to some change. /
	+ StopUpdate=3 /*< Stop recalculating at this spin info. /
	};

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* Default constructor.
	*/
	SpinInfo()
	: _timelike(false), _prodloc(-1), _decayloc(-1),
	_decayed(false), _developed(Undeveloped) {}

	/**
	* Standard Constructor.
	* @param s the spin.
	* @param p the production momentum.
	* @param time true if the particle is time-like.
	*/
	SpinInfo(PDT::Spin s,
	const Lorentz5Momentum & p = Lorentz5Momentum(),
	bool time = false)
	: _timelike(time), _prodloc(-1), _decayloc(-1),
	_decayed(false), _developed(Undeveloped),
	_rhomatrix(s), _Dmatrix(s), _spin(s),
	_productionmomentum(p), _currentmomentum(p) {}

	/**
	* Copy-constructor.
	*/
	SpinInfo(const SpinInfo &);
	//@}

	public:

	/**
	* Returns true if the polarization() has been implemented in a
	* subclass. This default version returns false.
	*/
	virtual bool hasPolarization() const { return false; }

	/**
	* Return the angles of the polarization vector as a pair of
	* doubles. first is the polar angle and second is the azimuth
	* wrt. the particles direction. This default version of the
	* function returns 0,0, and if a subclass implements a proper
	* function it should also implement 'hasPolarization()' to return
	* true.
	*/
	virtual DPair polarization() const { return DPair(); }

	public:

	/**
	* Standard Init function.
	*/
	static void Init();

	/**
	* Rebind to cloned objects. If a FermionSpinInfo is cloned together
	* with a whole Event and this has pointers to other event record
	* objects, these should be rebound to their clones in this
	* function.
	*/
	virtual void rebind(const EventTranslationMap & trans);

	/**
	* Standard clone method.
	*/
	virtual EIPtr clone() const;

	/**
	* Method to handle the delelation
	*/
	void update() const;

	/**
	* Perform a lorentz rotation of the spin information
	*/
	virtual void transform(const LorentzMomentum & m, const LorentzRotation & r) {
	_currentmomentum = m;
	_currentmomentum.transform(r);
	}

	public:


	/** @name Access the vertices. */
	//@{
	/**
	* Set the vertex at which the particle was produced.
	*/
	void productionVertex(VertexPtr in) const {
	_production=in;
	// add to the list of outgoing if timelike
	int temp(-1);
	if(_timelike) in->addOutgoing(this,temp);
	// or incoming if spacelike
	else in->addIncoming(this,temp);
	_prodloc=temp;
	}

	/**
	* Get the vertex at which the particle was produced.
	*/
	tcVertexPtr productionVertex() const { return _production; }

	/**
	* Set the vertex at which the particle decayed or branched.
	*/
	void decayVertex(VertexPtr in) const {
	if(in) {
	_decay=in;
	if(_timelike) {
	int temp(-1);
	in->addIncoming(this,temp);
	_decayloc=temp;
	assert(temp==0);
	}
	else {
	int temp(-1);
	in->addOutgoing(this,temp);
	_decayloc=temp;
	}
	}
	else {
	_decay=VertexPtr();
	_decayloc=-1;
	}
	}

	/**
	* Get the vertex at which the particle decayed or branched.
	*/
	tcVertexPtr decayVertex() const { return _decay; }
	//@}

	/** @name Access information about the associated particle. */
	//@{
	/**
	* Has the particle decayed?
	*/
	bool decayed() const { return _decayed; }

	/**
	* Set if the particle has decayed.
	*/
	void decayed(bool b) const { _decayed = b; }

	/**
	* Return true if the decay matrix required to perform the decays of
	* the siblings of a particle has been calculated.
	*/
	DevelopedStatus developed() const { return _developed; }

	/**
	* Calculate the rho matrix for the decay if not already done.
	*/
	void decay(bool recursive=false) const ;

	/**
	* Set the developed flag and calculate the D matrix for the decay.
	*/
	void develop() const ;

	/**
	* Needs update
	*/
	void needsUpdate() const {_developed=NeedsUpdate;}

	/**
	+ * Used for an unstable particle to temporarily stop
	+ * redevelop and redecay at that particle
	+ */
	+ void stopUpdate() const {_developed=StopUpdate;}
	+
	+ /**
	* Return 2s+1 for the particle
	*/
	PDT::Spin iSpin() const { return _spin; }

	/**
	* Return the momentum of the particle when it was produced.
	*/
	const Lorentz5Momentum & productionMomentum() const {
	return _productionmomentum;
	}

	/**
	* The current momentum of the particle
	*/
	const Lorentz5Momentum & currentMomentum() const {
	return _currentmomentum;
	}

	/**
	* Return true if particle is timelike (rather than spacelike).
	*/
	bool timelike() const { return _timelike; }
	//@}

	/**
	* Access to the locations
	*/
	//@{
	/**
	* Production Location
	*/
	int productionLocation() const {return _prodloc;}

	/**
	* Decay Location
	*/
	int decayLocation() const {return _decayloc;}
	//@}

	public:

	/** @name Access the rho and D matrices. */
	//@{
	/**
	* Access the rho matrix.
	*/
	RhoDMatrix rhoMatrix() const { return _rhomatrix; }

	/**
	* Access the rho matrix.
	*/
	RhoDMatrix & rhoMatrix() { return _rhomatrix; }

	/**
	* Access the D matrix.
	*/
	RhoDMatrix DMatrix() const { return _Dmatrix; }

	/**
	* Access the D matrix.
	*/
	RhoDMatrix & DMatrix() { return _Dmatrix; }
	//@}

	protected:

	/**
	* Check if momentum is near to the current momentum
	*/
	bool isNear(const Lorentz5Momentum & p) {
	return currentMomentum().isNear(p,_eps);
	}

	private:

	/**
	* Describe a concrete class without persistent data.
	*/
	static NoPIOClassDescription<SpinInfo> initSpinInfo;

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

	private:

	/**
	* Set the developed flag and calculate the D matrix for the decay,
	* and all decays further up the chain.
	*/
	void redevelop() const ;

	/**
	* Recursively recalulate all the rho matrices from the top of the chain
	*/
	void redecay() const ;

	private:

	/**
	* Pointer to the production vertex for the particle
	*/
	mutable VertexPtr _production;

	/**
	* Pointers to the decay vertex for the particle
	*/
	mutable VertexPtr _decay;

	/**
	* Is this is timelike (true) or spacelike (false ) particle? This
	* is used to decide if the particle is incoming or outgoing at the
	* production vertex
	*/
	bool _timelike;

	/**
	* Location in the hard vertex array at production.
	*/
	mutable int _prodloc;

	/**
	* Location in the hard vertex array at decay.
	*/
	mutable int _decayloc;

	/**
	* Has the particle been decayed? (I.e. has the rho matrix for the
	* decay been calculated.)
	*/
	mutable bool _decayed;

	/**
	* Has the particle been developed? (I.e. has the D matrix encoding
	* the info about the decay been calculated)
	*/
	mutable DevelopedStatus _developed;

	/**
	* Storage of the rho matrix.
	*/
	mutable RhoDMatrix _rhomatrix;

	/**
	* Storage of the decay matrix
	*/
	mutable RhoDMatrix _Dmatrix;

	/**
	* The spin of the particle
	*/
	PDT::Spin _spin;

	/**
	* Momentum of the particle when it was produced
	*/
	Lorentz5Momentum _productionmomentum;

	/**
	* Momentum of the particle when it decayed
	*/
	mutable Lorentz5Momentum _decaymomentum;

	/**
	* Current momentum of the particle
	*/
	Lorentz5Momentum _currentmomentum;

	/**
	* A small energy for comparing momenta to check if Lorentz Transformations
	* should be performed
	*/
	static const double _eps;
	};

	}


	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/**
	* This template specialization informs ThePEG about the name of the
	* SpinInfo class and the shared object where it is defined.
	*/
	template <>
	struct ClassTraits<ThePEG::SpinInfo>
	: public ClassTraitsBase<ThePEG::SpinInfo> {
	/**
	* Return the class name.
	*/
	static string className() { return "ThePEG::SpinInfo"; }
	};

	/** @endcond */

	}

	#endif /* ThePEG_SpinInfo_H */
	diff --git a/Utilities/SimplePhaseSpace.h b/Utilities/SimplePhaseSpace.h
	--- a/Utilities/SimplePhaseSpace.h
	+++ b/Utilities/SimplePhaseSpace.h
	@@ -1,232 +1,232 @@
	// -- C++ --
	//
	// SimplePhaseSpace.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_SimplePhaseSpace_H
	#define ThePEG_SimplePhaseSpace_H

	#include "ThePEG/Config/ThePEG.h"

	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/Vectors/LorentzRotation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/EventRecord/ParticleTraits.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "SimplePhaseSpace.xh"
	#include <numeric>

	namespace ThePEG {

	/**
	* SimplePhaseSpace defines a set of static functions to be used for
	* distributing momenta evenly in phase space. In most cases pointers
	* and references to both particle and momentum objects can be used as
	* arguments as long as the ParticleTraits class is specialized
	* properly. When needed, random numbers are generated with the
	* generator given by the static UseRandom class.
	*/
	namespace SimplePhaseSpace {

	/**
	* Set two momenta in their center of mass system. Their total
	* invariant mass squared is given by s, and their direction is
	* distributed isotropically.
	* @param s the total invariant mass squared.
	* @param p1 pointer or reference to the first momentum. Its
	* invariant mass will be preserved.
	* @param p2 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PType>
	void CMS(Energy2 s, PType & p1, PType & p2);

	/**
	* Set two momenta in their center of mass system. Their total
	* invariant mass squared is given by s, and their direction is
	* given in terms of the polar and azimuth angle of the first
	* momenta.
	* @param s the total invariant mass squared.
	* @param p1 pointer or reference to the first momentum. Its
	* invariant mass will be preserved.
	* @param p2 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param cosTheta cosine of the azimuth angle of the first momentum.
	* @param phi azimuth angle of the first momentum.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PType>
	void CMS(PType & p1, PType & p2, Energy2 s,
	double cosTheta, double phi);

	/**
	* Set two momenta in their center of mass system. Their total
	* invariant mass squared is given by s. The helper momentum p0 is
	* used so that afterwards \f$t=(p0-p1)^2\f$ and p1 has the azimuth
	* angle phi around p0.
	* @param p1 pointer or reference to the first momentum. Its
	* invariant mass will be preserved.
	* @param p2 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param s the total invariant mass squared.
	* @param t \f$=(p0-p1)^2\f$.
	* @param phi azimuth angle of the first momentum around p0.
	* @param p0 pointer or reference to an auxiliary momentum.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PType>
	void CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi,
	const PType & p0);

	/**
	* Set two momenta in their center of mass system. Their total
	* invariant mass squared is given by s. p1 will be along the z-axis.
	* @param p1 pointer or reference to the first momentum. Its
	* invariant mass will be preserved.
	* @param p2 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param s the total invariant mass squared.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PType>
	void CMS(PType & p1, PType & p2, Energy2 s);

	/**
	* Set two momenta in their center of mass system. Their total
	* invariant mass squared is given by s. The first will be along the
	* z-axis.
	* @param p a pair of pointers or references to the two momenta. Their
	* invariant masses will be preserved.
	* @param s the total invariant mass squared.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PPairType>
	void CMS(const PPairType & p, Energy2 s)
	{
	CMS(p.first, p.second, s);
	}

	/**
	* Set three momenta in their center of mass system. Their total
	* invariant mass squared is given by s. The energy fraction of
	* particle p1(3) is x1(3) of the total energy and the angles of the
	* system is distributed isotropically.
	* @param p1 pointer or reference to the first momentum. Its
	* invariant mass will be preserved.
	* @param p2 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param p3 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param s the total invariant mass squared.
	* @param x1 the energy fraction \f$2e_1/\sqrt{s}\f$.
	* @param x3 the energy fraction \f$2e_3/\sqrt{s}\f$.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PType>
	void CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
	double x1, double x3);

	/**
	* Set three momenta in their center of mass system. Their total
	* invariant mass squared is given by s. The energy fraction of
	* particle p1(3) is x1(3) of the total energy. Particle p1 is
	* initially placed along the z-axis and particle p2 is given
	* azimuth angle phii. Then the system is then rotated with
	* theta and phi respectively.
	* @param p1 pointer or reference to the first momentum. Its
	* invariant mass will be preserved.
	* @param p2 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param p3 pointer or reference to the second momentum. Its
	* invariant mass will be preserved.
	* @param s the total invariant mass squared.
	* @param x1 the energy fraction \f$2e_1/\sqrt{s}\f$.
	* @param x3 the energy fraction \f$2e_3/\sqrt{s}\f$.
	* @param phii the azimuth angle of p2 around p1.
	* @param theta the polar angle of p1.
	* @param phi the azimuth angle of p1.
	* @throw ImpossibleKinematics if the sum of the invariant masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename PType>
	void CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
	double x1, double x3, double phii = 0.0,
	double theta = 0.0, double phi = 0.0);

	/**
	* Calculate the absolute magnitude of the momenta of two particles
	* with masses m1 and m2 when put in their CMS of total invariant
	* mass squared s.
	* @param s the total invariant mass squared.
	* @param m1 the mass of particle 1.
	* @param m2 the mass of particle 2.
	* @throw ImpossibleKinematics if the sum of the masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	Energy getMagnitude(Energy2 s, Energy m1, Energy m2);

	/**
	* Return a three-vector given the absolute momentum, cos(theta) and
	* phi.
	* @param p the magnitude of the momentum.
	* @param costheta the cosine of the polar angle.
	* @param phi the azimuth angle.
	*/
	inline Momentum3 polar3Vector(Energy p, double costheta, double phi)
	{
	return Momentum3(psqrt(1.0 - sqr(costheta))sin(phi),
	psqrt(1.0 - sqr(costheta))cos(phi),
	p*costheta);
	}

	/**
	* Get a number of randomly distributed momenta.
	* Given a number specified invariant masses and a
	* total invariant mass m0, return corresponding four-momenta
	* randomly distributed according to phase space.
	* @param m0 the
	* total invariant mass of the resulting momenta.
	* @param m a vector
	* of invariant masses of the resulting momenta.
	* @return a vector
	* of momenta with the given masses randomly distributed.
	* @throw ImpossibleKinematics if the sum of the masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	vector<LorentzMomentum>
	CMSn(Energy m0, const vector<Energy> & m);

	/**
	* Set the momentum of a number of particles. Given a number of
	* particles and a total invariant mass m0, distribute their
	* four-momenta randomly according to phase space.
	* @param particles a container of particles or pointers to
	* particles. The invariant mass of these particles will not be
	* chaned.
	* @param m0 the
	* total invariant mass of the resulting momenta.
	* @throw ImpossibleKinematics if the sum of the masses was
	* larger than the given invariant mass (\f$\sqrt{s}\f$).
	*/
	template <typename Container>
	void CMSn(Container & particles, Energy m0);

	-};
	+}

	}

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

	#endif /* ThePEG_SimplePhaseSpace_H */

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:37 PM (1 d, 4 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4230692
Default Alt Text: (26 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions