CandUtils/src/Thrust.cc

0001 #include "PhysicsTools/CandUtils/interface/Thrust.h"
0002 #include <cmath>
0003 using namespace reco;
0004 const double pi = M_PI, pi2 = 2 * pi, pi_2 = pi / 2, pi_4 = pi / 4;
0005
0006 void Thrust::init(const std::vector<const Candidate*>& cands) {
0007   int i = 0;
0008   for (std::vector<const Candidate*>::const_iterator t = cands.begin(); t != cands.end(); ++t, ++i)
0009     pSum_ += (p_[i] = (*t)->momentum()).r();
0010   axis_ = axis(finalAxis(initialAxis()));
0011   if (axis_.z() < 0)
0012     axis_ *= -1;
0013   thrust_ = thrust(axis_);
0014 }
0015
0016 Thrust::ThetaPhi Thrust::initialAxis() const {
0017   static const int nSegsTheta = 10, nSegsPhi = 10, nSegs = nSegsTheta * nSegsPhi;
0018   int i, j;
0019   double thr[nSegs], max = 0;
0020   int indI = 0, indJ = 0, index = -1;
0021   for (i = 0; i < nSegsTheta; ++i) {
0022     double z = cos(pi * i / (nSegsTheta - 1));
0023     double r = sqrt(1 - z * z);
0024     for (j = 0; j < nSegsPhi; ++j) {
0025       double phi = pi2 * j / nSegsPhi;
0026       thr[i * nSegsPhi + j] = thrust(Vector(r * cos(phi), r * sin(phi), z));
0027       if (thr[i * nSegsPhi + j] > max) {
0028         index = i * nSegsPhi + j;
0029         indI = i;
0030         indJ = j;
0031         max = thr[index];
0032       }
0033     }
0034   }
0035
0036   // take max and one point on either size, fitting to a parabola and

0037   // extrapolating to the real max.  Do this separately for each dimension.

0038   // y = a x^2 + b x + c.  At the max, x = 0, on either side, x = +/-1.

0039   // do phi first

0040   double a, b, c = max;
0041   int ind1 = indJ + 1;
0042   if (ind1 >= nSegsPhi)
0043     ind1 -= nSegsPhi;
0044   int ind2 = indJ - 1;
0045   if (ind2 < 0)
0046     ind2 += nSegsPhi;
0047   a = (thr[ind1] + thr[ind2] - 2 * c) / 2;
0048   b = thr[ind1] - a - c;
0049   double maxPhiInd = 0;
0050   if (a != 0)
0051     maxPhiInd = -b / (2 * a);
0052   double maxThetaInd;
0053   if (indI == 0 || indI == (nSegsTheta - 1))
0054     maxThetaInd = indI;
0055   else {
0056     ind1 = indI + 1;
0057     ind2 = indI - 1;
0058     a = (thr[ind1] + thr[ind2] - 2 * c) / 2;
0059     b = thr[ind1] - a - c;
0060     maxThetaInd = 0;
0061     if (a != 0)
0062       maxThetaInd = -b / (2 * a);
0063   }
0064   return ThetaPhi(pi * (maxThetaInd + indI) / (nSegsTheta - 1), pi2 * (maxPhiInd + indJ) / nSegsPhi);
0065 }
0066
0067 Thrust::ThetaPhi Thrust::finalAxis(ThetaPhi best) const {
0068   static const double epsilon = 0.0001;
0069   double maxChange1 = 0.0, maxChange2 = 0.0, a = 0.0, b = 0.0, c = 0.0, thr = 0.0;
0070   int mandCt = 3, maxCt = 1000;
0071   bool done;
0072   do {
0073     parabola(a, b, c, axis(best), axis(best.theta + epsilon, best.phi), axis(best.theta - epsilon, best.phi));
0074     maxChange1 = 10 * (b < 0 ? -1 : 1);
0075     if (a != 0)
0076       maxChange1 = -b / (2 * a);
0077     while (fabs(maxChange1 * epsilon) > pi_4)
0078       maxChange1 /= 2;
0079     if (maxChange1 == 0 && (best.theta == 0 || best.theta == pi)) {
0080       best.phi += pi_2;
0081       if (best.phi > pi2)
0082         best.phi -= pi2;
0083       parabola(a, b, c, axis(best), axis(best.theta + epsilon, best.phi), axis(best.theta - epsilon, best.phi));
0084       maxChange1 = 10 * (b < 0 ? -1 : 1);
0085       if (a != 0)
0086         maxChange1 = -b / (2 * a);
0087     }
0088     do {
0089       // L.L.: fixed odd behavoir adding epsilon (???)

0090       thr = thrust(axis(best.theta + maxChange1 * epsilon, best.phi)) + epsilon;
0091       if (thr < c)
0092         maxChange1 /= 2;
0093     } while (thr < c);
0094
0095     best.theta += maxChange1 * epsilon;
0096     if (best.theta > pi) {
0097       best.theta = pi - (best.theta - pi);
0098       best.phi += pi;
0099       if (best.phi > pi2)
0100         best.phi -= pi2;
0101     }
0102     if (best.theta < 0) {
0103       best.theta *= -1;
0104       best.phi += pi;
0105       if (best.phi > pi2)
0106         best.phi -= pi2;
0107     }
0108     parabola(a, b, c, axis(best), axis(best.theta, best.phi + epsilon), axis(best.theta, best.phi - epsilon));
0109     maxChange2 = 10 * (b < 0 ? -1 : 1);
0110     if (a != 0)
0111       maxChange2 = -b / (2 * a);
0112     while (fabs(maxChange2 * epsilon) > pi_4) {
0113       maxChange2 /= 2;
0114     }
0115     do {
0116       // L.L.: fixed odd behavoir adding epsilon

0117       thr = thrust(axis(best.theta, best.phi + maxChange2 * epsilon)) + epsilon;
0118       if (thr < c)
0119         maxChange2 /= 2;
0120     } while (thr < c);
0121     best.phi += maxChange2 * epsilon;
0122     if (best.phi > pi2)
0123       best.phi -= pi2;
0124     if (best.phi < 0)
0125       best.phi += pi2;
0126     if (mandCt > 0)
0127       mandCt--;
0128     maxCt--;
0129     done = (fabs(maxChange1) > 1 || fabs(maxChange2) > 1 || mandCt) && (maxCt > 0);
0130   } while (done);
0131
0132   return best;
0133 }
0134
0135 void Thrust::parabola(double& a, double& b, double& c, const Vector& a1, const Vector& a2, const Vector& a3) const {
0136   double t1 = thrust(a1), t2 = thrust(a2), t3 = thrust(a3);
0137   a = (t2 - 2 * c + t3) / 2;
0138   b = t2 - a - c;
0139   c = t1;
0140 }
0141
0142 Thrust::Vector Thrust::axis(double theta, double phi) const {
0143   double theSin = sin(theta);
0144   return Vector(theSin * cos(phi), theSin * sin(phi), cos(theta));
0145 }
0146
0147 double Thrust::thrust(const Vector& axis) const {
0148   double result = 0;
0149   double sum = 0;
0150   for (unsigned int i = 0; i < n_; ++i)
0151     sum += fabs(axis.Dot(p_[i]));
0152   if (pSum_ > 0)
0153     result = sum / pSum_;
0154   return result;
0155 }