Thrust Class Reference

#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...
template<typename const_iterator >
	Thrust (const_iterator begin, const_iterator end)
	constructor from first and last iterators More...
double	thrust () const
	thrust value (in the range [0.5, 1.0]) More...

Private Member Functions
Vector	axis (double theta, double phi) const
Vector	axis (const ThetaPhi &tp) 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

typedef math::XYZVector Thrust::Vector

spatial vector

Definition at line 41 of file Thrust.h.

Constructor & Destructor Documentation

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.

References end, i, init(), and n_.

                                                   :
    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);
  } 

Member Function Documentation

const Vector& Thrust::axis ( ) const

inline

thrust axis (with magnitude = 1)

Definition at line 57 of file Thrust.h.

References axis_.

Referenced by axis().

{ return axis_; } 

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

private

Definition at line 134 of file Thrust.cc.

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

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

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

inlineprivate

Definition at line 74 of file Thrust.h.

References axis(), Thrust::ThetaPhi::phi, and Thrust::ThetaPhi::theta.

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

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

private

Definition at line 63 of file Thrust.cc.

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

                                                    {
  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;
}

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

private

Definition at line 6 of file Thrust.cc.

References axis_, i, alignCSCRings::r, and edmStreamStallGrapher::t.

Referenced by Thrust().

                                                          {
  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_);
}

Thrust::ThetaPhi Thrust::initialAxis ( ) const

private

Definition at line 16 of file Thrust.cc.

References a, b, trackerHits::c, funct::cos(), i, cmsHarvester::index, j, bookConverter::max, phi, pi, pi2, alignCSCRings::r, funct::sin(), mathSSE::sqrt(), and detailsBasic3DVector::z.

                                         {
  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);
}

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

private

Definition at line 126 of file Thrust.cc.

References trackerHits::c.

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

double Thrust::thrust ( ) const

inline

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

Definition at line 55 of file Thrust.h.

References thrust_.

Referenced by Rivet::CMS_EWK_11_021::analyze().

{ return thrust_; } 

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

private

Definition at line 139 of file Thrust.cc.

References i, and query::result.

                                               {
  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;
}

Member Data Documentation

Vector Thrust::axis_

private

Definition at line 61 of file Thrust.h.

Referenced by axis().

const unsigned int Thrust::n_

private

Definition at line 63 of file Thrust.h.

Referenced by Thrust().

std::vector<Vector> Thrust::p_

private

Definition at line 64 of file Thrust.h.

double Thrust::pSum_

private

Definition at line 62 of file Thrust.h.

double Thrust::thrust_

private

Definition at line 60 of file Thrust.h.

Referenced by thrust().