#include <Thrust.h>

Classes
struct	ThetaPhi

Public Types
typedef math::XYZVector	Vector
	spatial vector More...

Public Member Functions
const Vector &	axis () const
	thrust axis (with magnitude = 1) More...

double	thrust () const
	thrust value (in the range [0.5, 1.0]) More...

template<typename const_iterator >
	Thrust (const_iterator begin, const_iterator end)
	constructor from first and last iterators More...

Private Member Functions
Vector	axis (const ThetaPhi &tp) const

Vector	axis (double theta, double phi) const

ThetaPhi	finalAxis (ThetaPhi) const

void	init (const std::vector< const reco::Candidate * > &)

ThetaPhi	initialAxis () const

void	parabola (double &a, double &b, double &c, const Vector &, const Vector &, const Vector &) const

double	thrust (const Vector &theAxis) const

Private Attributes
Vector	axis_

const unsigned int	n_

std::vector< Vector >	p_

double	pSum_

double	thrust_

Detailed Description

Utility to compute thrust value and axis from a collection of candidates.

Ported from BaBar implementation.

The thrust axis is the vector T which maximises the following expression:

  sum_i=1...N ( | T . p_i | )

t = ------------------------------— sum_i=1...N ( (p_i . _p_i)^(1/2) )

where p_i, i=1...N are the particle momentum vectors. The thrust value is the maximum value of t.

The thrust axis has a two-fold ambiguity due to taking the absolute value of the dot product. This computation returns by convention a thurs axis with a positive component along the z-direction is positive.

The thrust measure the alignment of a collection of particles along a common axis. The lower the thrust, the more spherical the event is. The higher the thrust, the more jet-like the

Author: Luca Lista, INFN

Definition at line 38 of file Thrust.h.

Member Typedef Documentation

◆ Vector

typedef math::XYZVector Thrust::Vector

spatial vector

Definition at line 41 of file Thrust.h.

Constructor & Destructor Documentation

◆ Thrust()

template<typename const_iterator >

Thrust::Thrust	(	const_iterator	begin,
		const_iterator	end
	)

inline

constructor from first and last iterators

Definition at line 44 of file Thrust.h.

                                                    : thrust_(0), axis_(0, 0, 0), pSum_(0), n_(end - begin), p_(n_) {
     if (n_ == 0)
       return;
     std::vector<const reco::Candidate *> cands;
     for (const_iterator i = begin; i != end; ++i) {
       cands.push_back(&*i);
     }
     init(cands);
   }

References HLT_FULL_cff::cands, mps_fire::end, mps_fire::i, init(), and n_.

Member Function Documentation

◆ axis() [1/3]

const Vector& Thrust::axis ( ) const

inline

thrust axis (with magnitude = 1)

Definition at line 56 of file Thrust.h.

56 { return axis_; }

References axis_.

◆ axis() [2/3]

Vector Thrust::axis ( const ThetaPhi & tp ) const

inlineprivate

Definition at line 73 of file Thrust.h.

73 { return axis(tp.theta, tp.phi); }

References axis(), and cmsswSequenceInfo::tp.

Referenced by axis().

◆ axis() [3/3]

Thrust::Vector Thrust::axis	(	double	theta,
		double	phi
	)		const

private

Definition at line 142 of file Thrust.cc.

                                                         {
   double theSin = sin(theta);
   return Vector(theSin * cos(phi), theSin * sin(phi), cos(theta));
 }

References funct::cos(), funct::sin(), and theta().

◆ finalAxis()

Thrust::ThetaPhi Thrust::finalAxis ( ThetaPhi best ) const

private

Definition at line 67 of file Thrust.cc.

                                                     {
   static const double epsilon = 0.0001;
   double maxChange1 = 0.0, maxChange2 = 0.0, a = 0.0, b = 0.0, c = 0.0, thr = 0.0;
   int mandCt = 3, maxCt = 1000;
   bool done;
   do {
     parabola(a, b, c, axis(best), axis(best.theta + epsilon, best.phi), axis(best.theta - epsilon, best.phi));
     maxChange1 = 10 * (b < 0 ? -1 : 1);
     if (a != 0)
       maxChange1 = -b / (2 * a);
     while (fabs(maxChange1 * epsilon) > pi_4)
       maxChange1 /= 2;
     if (maxChange1 == 0 && (best.theta == 0 || best.theta == pi)) {
       best.phi += pi_2;
       if (best.phi > pi2)
         best.phi -= pi2;
       parabola(a, b, c, axis(best), axis(best.theta + epsilon, best.phi), axis(best.theta - epsilon, best.phi));
       maxChange1 = 10 * (b < 0 ? -1 : 1);
       if (a != 0)
         maxChange1 = -b / (2 * a);
     }
     do {
       // L.L.: fixed odd behavoir adding epsilon (???)
       thr = thrust(axis(best.theta + maxChange1 * epsilon, best.phi)) + epsilon;
       if (thr < c)
         maxChange1 /= 2;
     } while (thr < c);
  
     best.theta += maxChange1 * epsilon;
     if (best.theta > pi) {
       best.theta = pi - (best.theta - pi);
       best.phi += pi;
       if (best.phi > pi2)
         best.phi -= pi2;
     }
     if (best.theta < 0) {
       best.theta *= -1;
       best.phi += pi;
       if (best.phi > pi2)
         best.phi -= pi2;
     }
     parabola(a, b, c, axis(best), axis(best.theta, best.phi + epsilon), axis(best.theta, best.phi - epsilon));
     maxChange2 = 10 * (b < 0 ? -1 : 1);
     if (a != 0)
       maxChange2 = -b / (2 * a);
     while (fabs(maxChange2 * epsilon) > pi_4) {
       maxChange2 /= 2;
     }
     do {
       // L.L.: fixed odd behavoir adding epsilon
       thr = thrust(axis(best.theta, best.phi + maxChange2 * epsilon)) + epsilon;
       if (thr < c)
         maxChange2 /= 2;
     } while (thr < c);
     best.phi += maxChange2 * epsilon;
     if (best.phi > pi2)
       best.phi -= pi2;
     if (best.phi < 0)
       best.phi += pi2;
     if (mandCt > 0)
       mandCt--;
     maxCt--;
     done = (fabs(maxChange1) > 1 || fabs(maxChange2) > 1 || mandCt) && (maxCt > 0);
   } while (done);
  
   return best;
 }

References a, b, c, fileCollector::done, geometryDiff::epsilon, Thrust::ThetaPhi::phi, pi, pi2, pi_2, pi_4, and Thrust::ThetaPhi::theta.

◆ init()

void Thrust::init ( const std::vector< const reco::Candidate * > & cands )

private

Definition at line 6 of file Thrust.cc.

                                                           {
   int i = 0;
   for (std::vector<const Candidate*>::const_iterator t = cands.begin(); t != cands.end(); ++t, ++i)
     pSum_ += (p_[i] = (*t)->momentum()).r();
   axis_ = axis(finalAxis(initialAxis()));
   if (axis_.z() < 0)
     axis_ *= -1;
   thrust_ = thrust(axis_);
 }

References HLT_FULL_cff::cands, mps_fire::i, alignCSCRings::r, and submitPVValidationJobs::t.

Referenced by Thrust().

◆ initialAxis()

Thrust::ThetaPhi Thrust::initialAxis ( ) const

private

Definition at line 16 of file Thrust.cc.

                                          {
   static const int nSegsTheta = 10, nSegsPhi = 10, nSegs = nSegsTheta * nSegsPhi;
   int i, j;
   double thr[nSegs], max = 0;
   int indI = 0, indJ = 0, index = -1;
   for (i = 0; i < nSegsTheta; ++i) {
     double z = cos(pi * i / (nSegsTheta - 1));
     double r = sqrt(1 - z * z);
     for (j = 0; j < nSegsPhi; ++j) {
       double phi = pi2 * j / nSegsPhi;
       thr[i * nSegsPhi + j] = thrust(Vector(r * cos(phi), r * sin(phi), z));
       if (thr[i * nSegsPhi + j] > max) {
         index = i * nSegsPhi + j;
         indI = i;
         indJ = j;
         max = thr[index];
       }
     }
   }
  
   // take max and one point on either size, fitting to a parabola and
   // extrapolating to the real max.  Do this separately for each dimension.
   // y = a x^2 + b x + c.  At the max, x = 0, on either side, x = +/-1.
   // do phi first
   double a, b, c = max;
   int ind1 = indJ + 1;
   if (ind1 >= nSegsPhi)
     ind1 -= nSegsPhi;
   int ind2 = indJ - 1;
   if (ind2 < 0)
     ind2 += nSegsPhi;
   a = (thr[ind1] + thr[ind2] - 2 * c) / 2;
   b = thr[ind1] - a - c;
   double maxPhiInd = 0;
   if (a != 0)
     maxPhiInd = -b / (2 * a);
   double maxThetaInd;
   if (indI == 0 || indI == (nSegsTheta - 1))
     maxThetaInd = indI;
   else {
     ind1 = indI + 1;
     ind2 = indI - 1;
     a = (thr[ind1] + thr[ind2] - 2 * c) / 2;
     b = thr[ind1] - a - c;
     maxThetaInd = 0;
     if (a != 0)
       maxThetaInd = -b / (2 * a);
   }
   return ThetaPhi(pi * (maxThetaInd + indI) / (nSegsTheta - 1), pi2 * (maxPhiInd + indJ) / nSegsPhi);
 }

References a, b, c, funct::cos(), mps_fire::i, dqmiolumiharvest::j, SiStripPI::max, pi, pi2, alignCSCRings::r, funct::sin(), and mathSSE::sqrt().

◆ parabola()

void Thrust::parabola	(	double &	a,
		double &	b,
		double &	c,
		const Vector &	a1,
		const Vector &	a2,
		const Vector &	a3
	)		const

private

Definition at line 135 of file Thrust.cc.

                                                                                                                  {
   double t1 = thrust(a1), t2 = thrust(a2), t3 = thrust(a3);
   a = (t2 - 2 * c + t3) / 2;
   b = t2 - a - c;
   c = t1;
 }

References a, testProducerWithPsetDescEmpty_cfi::a2, b, c, RandomServiceHelper::t1, RandomServiceHelper::t2, and RandomServiceHelper::t3.

◆ thrust() [1/2]

double Thrust::thrust ( ) const

inline

thrust value (in the range [0.5, 1.0])

Definition at line 54 of file Thrust.h.

54 { return thrust_; }

References thrust_.

◆ thrust() [2/2]

double Thrust::thrust ( const Vector & theAxis ) const

private

Definition at line 147 of file Thrust.cc.

                                               {
   double result = 0;
   double sum = 0;
   for (unsigned int i = 0; i < n_; ++i)
     sum += fabs(axis.Dot(p_[i]));
   if (pSum_ > 0)
     result = sum / pSum_;
   return result;
 }

References mps_fire::i, and mps_fire::result.

Member Data Documentation

◆ axis_

Vector Thrust::axis_

private

Definition at line 60 of file Thrust.h.

Referenced by axis().

◆ n_

const unsigned int Thrust::n_

private

Definition at line 62 of file Thrust.h.

Referenced by Thrust().

◆ p_

std::vector<Vector> Thrust::p_

private

Definition at line 63 of file Thrust.h.

◆ pSum_

double Thrust::pSum_

private

Definition at line 61 of file Thrust.h.

◆ thrust_

double Thrust::thrust_

private

Definition at line 59 of file Thrust.h.

Referenced by thrust().

Classes

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ Vector

Constructor & Destructor Documentation

◆ Thrust()

Member Function Documentation

◆ axis() [1/3]

◆ axis() [2/3]

◆ axis() [3/3]

◆ finalAxis()

◆ init()

◆ initialAxis()

◆ parabola()

◆ thrust() [1/2]

◆ thrust() [2/2]

Member Data Documentation

◆ axis_

◆ n_

◆ p_

◆ pSum_

◆ thrust_