EventShapes.h
No OneTemporary
Actions

Size

14 KB

Referenced Files

None

Subscribers

None

EventShapes.h
View Options

	// -- C++ --
	//
	// EventShapes.h is a part of Herwig++ - A multi-purpose Monte Carlo
	// event generator Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for
	// details. Please respect the MCnet academic guidelines, see
	// GUIDELINES for details.
	//
	#ifndef HERWIG_EventShapes_H
	#define HERWIG_EventShapes_H
	//
	// This is the declaration of the EventShapes class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Handlers/AnalysisHandler.h"
	#include "ThePEG/Vectors/Lorentz5Vector.h"
	#include "ThePEG/Vectors/ThreeVector.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "EventShapes.fh"

	namespace Herwig {

	using namespace ThePEG;

	/** \ingroup Analysis
	*
	* The EventShapes class is designed so that certain event shapes, such
	* as the thrust are only calculated once per event given the speed of
	* the calculation.
	*
	* @see \ref EventShapesInterfaces "The interfaces" defined for
	* EventShapes.
	*/
	class EventShapes: public Interfaced {

	public:

	/**
	* The default constructor.
	*/
	EventShapes() : _thrustDone(false), _spherDone(false), _linTenDone(false),
	_hemDone(false), _useCmBoost(false),
	_mPlus(), _mMinus(), _bPlus(), _bMinus()
	{}

	/**
	* Member to reset the particles to be considered
	*/
	void reset(const tPVector &part) {
	_pv.resize(part.size());
	for(unsigned int ix=0;ix<part.size();++ix) _pv[ix]=part[ix]->momentum();
	_thrustDone = false;
	_spherDone = false;
	_linTenDone = false;
	_hemDone = false;
	_useCmBoost = false;
	}


	public:

	/**
	* Member functions to return thrust related shapes
	*/
	//@{
	/**
	* The thrust
	*/
	double thrust() {
	checkThrust();
	return _thrust[0];
	}

	/**
	* The major
	*/
	double thrustMajor() {
	checkThrust();
	return _thrust[1];
	}

	/**
	* The minor
	*/
	double thrustMinor() {
	checkThrust();
	return _thrust[2];
	}

	/**
	* The oblateness
	*/
	double oblateness() {
	checkThrust();
	return _thrust[1]-_thrust[2];
	}

	/**
	* The thrust axis
	*/
	Axis thrustAxis() {
	checkThrust();
	return _thrustAxis[0];
	}

	/**
	* The major axis
	*/
	Axis majorAxis() {
	checkThrust();
	return _thrustAxis[1];
	}

	/**
	* The minor axis
	*/
	Axis minorAxis() {
	checkThrust();
	return _thrustAxis[2];
	}
	//@}

	/**
	* Linear momentum tensor related event shapes
	*/
	//@{
	/**
	* The C parameter
	*/
	double CParameter() {
	checkLinTen();
	return 3.(_linTen[0]_linTen[1]+_linTen[1]*_linTen[2]
	+_linTen[2]*_linTen[0]);
	}

	/**
	* The D parameter
	*/
	double DParameter() {
	checkLinTen();
	return 27.(_linTen[0]_linTen[1]*_linTen[2]);
	}

	/**
	* The eigenvalues in descending order
	*/
	vector<double> linTenEigenValues() {
	checkLinTen();
	return _linTen;
	}


	/**
	* The eigenvectors in order of descending eigenvalue
	*/
	vector<Axis> linTenEigenVectors() {
	checkLinTen();
	return _linTenAxis;
	}

	//@}

	/**
	* Quadratic momentum tensor related variables
	*/
	//@{
	/**
	* The sphericity
	*/
	double sphericity() {
	checkSphericity();
	return 3./2.*(_spher[1]+_spher[2]);
	}

	/**
	* The aplanarity
	*/
	double aplanarity() {
	checkSphericity();
	return 3./2.*_spher[2];
	}


	/**
	* The planarity
	*/
	double planarity() {
	checkSphericity();
	return _spher[1]-_spher[2];
	}

	/**
	* The sphericity axis
	*/
	Axis sphericityAxis() {
	checkSphericity();
	return _spherAxis[0];
	}


	/**
	* The sphericity eigenvalues
	*/
	vector<double> sphericityEigenValues() {
	checkSphericity();
	return _spher;
	}

	/**
	* The sphericity eigenvectors
	*/
	vector<Axis> sphericityEigenVectors() {
	checkSphericity();
	return _spherAxis;
	} //@}

	/**
	* Jet mass related event shapes
	*/
	//@{
	/**
	* The high hemishpere mass squared divided by the visible energy
	* squared
	*/
	double Mhigh2() {
	checkHemispheres();
	return _mPlus;
	}

	/**
	* The low hemishpere mass squared divided by the visible energy
	* squared
	*/
	double Mlow2() {
	checkHemispheres();
	return _mMinus;
	}

	/**
	* The difference between the
	* hemishpere masses squared divided by the visible energy squared
	*/
	double Mdiff2() {
	checkHemispheres();
	return _mPlus-_mMinus;
	}

	//@}

	/**
	* Jet broadening related event shapes
	*/
	//@{
	/**
	* The wide jet broadening
	*/
	double Bmax() {
	checkHemispheres();
	return _bPlus;
	}

	/**
	* The narrow jet broadening
	*/
	double Bmin() {
	checkHemispheres();
	return _bMinus;
	}

	/**
	* The sum of the jet broadenings
	*/
	double Bsum() {
	checkHemispheres();
	return _bPlus+_bMinus;
	}


	/**
	* The difference of the jet broadenings
	*/
	double Bdiff() {
	checkHemispheres();
	return _bPlus-_bMinus;
	}
	//@}

	/**
	* Single particle variables which do not depend on event shapes axes
	*/
	//@{

	/**
	* The scaled momentum \f$\xi=-\log\left( p/E_{\rm beam}\right)\f$.
	*/
	double getXi(const Lorentz5Momentum & p,
	const Energy & Ebeam) {
	return((Ebeam > 0MeV && p.vect().mag() > 0MeV) ?
	log(Ebeam/p.vect().mag()) : -1.);
	}

	/**
	* Transverse momentum with respect to the beam
	*/
	Energy getPt(const Lorentz5Momentum & p) {
	return p.perp();
	}

	/**
	* Rapidity with respect to the beam direction
	*/
	double getRapidity(const Lorentz5Momentum & p) {
	return (p.t() > p.z() ? p.rapidity() : 1e99);
	}
	//@}

	/**
	* Single particle variables related to one of the shape axis.
	*/
	//@{
	/**
	* Transverse momentum with respect to the thrust axis in the event plane
	*/
	Energy ptInT(const Lorentz5Momentum & p) {
	checkThrust();
	return p.vect()*_thrustAxis[1];
	}

	/**
	* Transverse momentum with respect to the thrust axis out of the
	* event plane
	*/
	Energy ptOutT(const Lorentz5Momentum & p) {
	checkThrust();
	return p.vect()*_thrustAxis[2];
	}

	/**
	* Rapidity with respect to the thrust axis
	*/
	double yT(const Lorentz5Momentum & p) {
	checkThrust();
	return (p.t() > p.vect()*_thrustAxis[0] ?
	p.rapidity(_thrustAxis[0]) : 1e99);
	}

	/**
	* Transverse momentum with respect to the sphericity axis in the
	* event plane
	*/
	Energy ptInS(const Lorentz5Momentum & p) {
	checkSphericity();
	return p.vect()*_spherAxis[1];
	}

	/**
	* Transverse momentum with respect to the sphericity axis out of the
	* event plane
	*/
	Energy ptOutS(const Lorentz5Momentum & p) {
	checkSphericity();
	return p.vect()*_spherAxis[2];
	}

	/**
	* Rapidity with respect to the sphericity axis
	*/
	double yS(const Lorentz5Momentum & p) {
	checkSphericity();
	return (p.t() > p.vect()*_spherAxis[0] ?
	p.rapidity(_spherAxis[0]) : 1e99);
	}
	//@}


	/**
	* Energy-energy correlation (EEC) @param hi is the histogram and has
	* to be provided externally It is understood that the range of the
	* histogam is -1 < cos(chi) < 1. hi.front() contains the bin [-1 <
	* cos(chi) < -1+delta] and hi.back() the bin [1-delta < cos(chi) <
	* 1]. delta = 2/hi.size(). We use classical indices to access the
	* vector.
	*/
	void bookEEC(vector<double> & hi);

	/**
	* Before writing the histogram it has to be normalized according to
	* the number of events.
	*/
	void normalizeEEC(vector<double> & hi, long evts) {
	for (unsigned int bin = 0; bin < hi.size(); bin++) bin /= (hi.size()*evts);
	}

	/**
	* The asymmetry of EEC is calculated from a given \f$\cos\chi\f$ and
	* EEC histogram, which is a vector<double> as described above.
	*/
	double AEEC(vector<double> & hi, double& coschi) {
	if (coschi > 0. && coschi <= 1.) {
	int i = static_cast<int>( floor((-coschi+1.)/2.*hi.size()) );
	int j = static_cast<int>( floor(( coschi+1.)/2.*hi.size()) );
	return hi[i]-hi[j];
	} else {
	return 1e99;
	}
	}

	public:

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

	protected:

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

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

	private:

	/**
	* Check whether the initialization of a certain class of event shapes
	* has been calculated and if not do so
	*/
	//@{
	/**
	* Check if thrust related variables have been calculated and if not
	* do so
	*/
	void checkThrust() {
	if (!_thrustDone) {
	_thrustDone = true;
	calculateThrust();
	}
	}

	/**
	* Check if the linear tensor related variables have been calculated
	* and if not do so
	*/
	void checkLinTen() {
	if (!_linTenDone) {
	_linTenDone = true;
	diagonalizeTensors(true, _useCmBoost);
	}
	}

	/**
	* Check if the quadratic tensor related variables have been
	* calculated and if not do so
	*/
	void checkSphericity() {
	if (!_spherDone) {
	_spherDone = true;
	diagonalizeTensors(false, _useCmBoost);
	}
	}

	/**
	* Check if the hemisphere mass variables and jet broadenings have
	* been calculated and if not do so
	*/
	void checkHemispheres() {
	if (!_hemDone) {
	_hemDone = true;
	calcHemisphereMasses();
	}
	}
	//@}

	/**
	* Methods that actually calculate the event shapes
	*/
	//@{
	/**
	* Calculate the hemisphere masses and jet broadenings
	*/
	void calcHemisphereMasses();

	/**
	* Calculate the thrust and related axes
	*/
	void calculateThrust();

	/**
	* Diagonalize the tensors @param linear switch between
	* diagonalization of linear/quadratic tensor. @param cmboost tells
	* whether to boost into cm frame of all momenta first, or not
	* (default off, and no interface to this).
	*/
	void diagonalizeTensors(bool linear, bool cmboost);

	/**
	* Quite general diagonalization of a symmetric Matrix T, given as an
	* array of doubles. The symmetry is not checked explicitly as this
	* is clear in the context. It uses an explicit generic solution of
	* the eigenvalue problem and no numerical approximation, based on
	* Cardano's formula. @param T Matrix to be diagonalised
	*/
	vector<double> eigenvalues(const double T[3][3]);

	/**
	* The eigenvector of @param T to a given eigenvalue @param lam
	*/
	Axis eigenvector(const double T[3][3], const double &lam);

	/**
	* The eigenvectors of @param T corresponding to the eigenvectors
	* @param lam . The ordering of the vectors corresponds to the
	* ordering of the eigenvalues.
	*/
	vector<Axis> eigenvectors(const double T[3][3], const vector<double> &lam);

	/**
	* Member to calculate the thrust
	* @param p The three vectors
	* @param t The thrust-squared (up to an Energy scale factor)
	* @param taxis The thrust axis
	*/
	void calcT(const vector<Momentum3> &p, Energy2 &t, Axis &taxis);

	/**
	* Member to calculate the major
	* @param p The three vectors
	* @param m The major-squared (up to an Energy scale factor)
	* @param maxis The major axis
	*/
	void calcM(const vector<Momentum3> &p, Energy2 &m, Axis &maxis);
	//@}

	private:

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

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

	private:

	/**
	* Vector of particle momenta to be analysed
	*/
	vector<Lorentz5Momentum> _pv;

	/**
	* Various event shape axes
	*/
	//@{
	/**
	* The thrust related axes
	*/
	vector<Axis> _thrustAxis;

	/**
	* The sphericity related axes
	*/
	vector<Axis> _spherAxis;

	/**
	* The linearised tensor axes
	*/
	vector<Axis> _linTenAxis;
	//@}

	/**
	* Values of axis related event shapes
	*/
	//@{
	/**
	* Values of thrust related variables
	*/
	vector<double> _thrust;

	/**
	* Values of sphericity related variables
	*/
	vector<double> _spher;

	/**
	* Values of linearized tensor related variables
	*/
	vector<double> _linTen;
	//@}

	/**
	* Whether or not certain event axes have been calculated
	*/
	//@{
	/**
	* Whether or not the thrust is calculated
	*/
	bool _thrustDone;

	/**
	* Whether or not the sphericity is calculated
	*/
	bool _spherDone;

	/**
	* Whether or not the linearizes tensor is calculated
	*/
	bool _linTenDone;

	/**
	* Whether or not the hemisphere masses have been calculated
	*/
	bool _hemDone;
	//@}

	/**
	* Whether ot not to boost to the CMS frame for the tensor diagonalizations
	*/
	bool _useCmBoost;

	/**
	* Hemisphere masses
	*/
	//@{
	/**
	* The high hemisphere mass
	*/
	double _mPlus;

	/**
	* The low hemisphere mass
	*/
	double _mMinus;
	//@}

	/**
	* The jet broadenings
	*/
	//@{
	/**
	* The wide jet broadening
	*/
	double _bPlus;

	/**
	* The narrow jet broadening
	*/
	double _bMinus;
	//@}
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

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

	/** This template specialization informs ThePEG about the name of
	* the EventShapes class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::EventShapes>
	: public ClassTraitsBase<Herwig::EventShapes> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::EventShapes"; }
	/** Return the name(s) of the shared library (or libraries) be loaded to get
	* access to the EventShapes class and any other class on which it depends
	* (except the base class). */
	static string library() { return "HwAnalysis.so"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_EventShapes_H */

File Metadata

Mime Type: text/x-c++
Expires: Tue, Sep 30, 4:46 AM (15 h, 24 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6449935
Default Alt Text: EventShapes.h (14 KB)

EventShapes.hNo OneTemporaryActions

EventShapes.hView Options

File Metadata

Event Timeline

EventShapes.h
No OneTemporary
Actions

EventShapes.h
View Options