CandUtils/interface/Thrust.h

0001 #ifndef CandUtils_Thrust_h
0002 #define CandUtils_Thrust_h
0003 /** \class Thrust

0004  *

0005  * Utility to compute thrust value and axis from

0006  * a collection of candidates.

0007  *

0008  * Ported from BaBar implementation.

0009  *

0010  * The thrust axis is the vector T which maximises

0011  * the following expression:

0012  *

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

0014  *  t = ---------------------------------

0015  *      sum_i=1...N ( (p_i . _p_i)^(1/2) )

0016  *

0017  * where p_i, i=1...N are the particle momentum vectors.

0018  * The thrust value is the maximum value of t.

0019  *

0020  * The thrust axis has a two-fold ambiguity due to taking the

0021  * absolute value of the dot product. This computation returns

0022  * by convention a thurs axis with a positive component along

0023  * the z-direction is positive.

0024  *

0025  * The thrust measure the alignment of a collection of particles

0026  * along a common axis.  The lower the thrust, the more spherical

0027  * the event is.  The higher the thrust, the more jet-like the

0028  *

0029  * \author Luca Lista, INFN

0030  *

0031  *

0032  *

0033  */
0034 #include "DataFormats/Math/interface/Vector3D.h"
0035 #include "DataFormats/Candidate/interface/Candidate.h"
0036 #include <vector>
0037
0038 class Thrust {
0039 public:
0040   /// spatial vector

0041   typedef math::XYZVector Vector;
0042   /// constructor from first and last iterators

0043   template <typename const_iterator>
0044   Thrust(const_iterator begin, const_iterator end) : thrust_(0), axis_(0, 0, 0), pSum_(0), n_(end - begin), p_(n_) {
0045     if (n_ == 0)
0046       return;
0047     std::vector<const reco::Candidate *> cands;
0048     for (const_iterator i = begin; i != end; ++i) {
0049       cands.push_back(&*i);
0050     }
0051     init(cands);
0052   }
0053   /// thrust value (in the range [0.5, 1.0])

0054   double thrust() const { return thrust_; }
0055   /// thrust axis (with magnitude = 1)

0056   const Vector &axis() const { return axis_; }
0057
0058 private:
0059   double thrust_;
0060   Vector axis_;
0061   double pSum_;
0062   const unsigned int n_;
0063   std::vector<Vector> p_;
0064
0065   struct ThetaPhi {
0066     ThetaPhi(double t, double p) : theta(t), phi(p) {}
0067     double theta, phi;
0068   };
0069   double thrust(const Vector &theAxis) const;
0070   ThetaPhi initialAxis() const;
0071   ThetaPhi finalAxis(ThetaPhi) const;
0072   Vector axis(double theta, double phi) const;
0073   Vector axis(const ThetaPhi &tp) const { return axis(tp.theta, tp.phi); }
0074   void parabola(double &a, double &b, double &c, const Vector &, const Vector &, const Vector &) const;
0075   void init(const std::vector<const reco::Candidate *> &);
0076 };
0077
0078 #endif