Thrust Class Reference

Class to calculate the thrust axis. More...

#include <Thrust.h>

Static Public Member Functions
static ROOT::Math::XYZVector	calculateThrust (const std::vector< ROOT::Math::XYZVector > &momenta)
	calculates the thrust axis

Detailed Description

Class to calculate the thrust axis.

Definition at line 22 of file Thrust.h.

Member Function Documentation

calculateThrust()

XYZVector calculateThrust ( const std::vector< ROOT::Math::XYZVector > &  momenta )

static

calculates the thrust axis

Definition at line 71 of file Thrust.cc.

{
 
  /* STEP 1: Initialization of Variables */
  XYZVector thrust_axis, trial_axis, base_axis;
  auto end = momenta.end();
  auto begin = momenta.begin();
  decltype(begin) itr;
 
  double sum_magnitude_mom = 0.;
  double thrust = 0.;
 
  /*
    STEP 2: Parse momenta vector to compute magnitude Σ(||p_i||)
  */
 
  for (auto const& momentum : momenta)
    sum_magnitude_mom += momentum.R();
 
  /*
    STEP 3: For each momentum in momenta vector,
            use momentum as initial axis to optimize
  */
 
  for (auto const& mom : momenta) {
    // By convention, thrust axis in same direction as Z axis
    trial_axis = (mom.Z() >= 0.) ? XYZVector(mom) : XYZVector(-mom);
 
    // Normalize if magnitude != 0
    trial_axis = trial_axis.Unit();
 
    /*
      STEP 4: Store the previous trial axis as a base axis and initialize
              a new trial axis as the Z-aligned sum of the momentum vectors
    */
 
    itr = begin;
    while (itr != end) {
      base_axis = trial_axis;
      trial_axis = XYZVector();
 
      // Z-alignment of momenta and addition to trial axis
      for (auto const& momentum : momenta)
        trial_axis += (momentum.Dot(base_axis) >= 0.) ? \
                      XYZVector(momentum) : XYZVector(-momentum);
 
      /*
        STEP 5: Check ( p_i · trial_axis ) * ( p_i · base_axis ) < 0 ∀ p_i
      */
 
      for (itr = begin; itr != end; itr++)
        if ((*itr).Dot(trial_axis) * (*itr).Dot(base_axis) < 0.) break;
 
      /*
        STEP 6: While condition True:
                    base_axis = trial_axis
                    repeat STEP 4 and STEP 5;
                else:
                    go to STEP 7 (while from STEP 4 is done)
      */
    }
 
    double trial_mag(trial_axis.R()); // pre-compute
    // trial_axis = trial_axis.Unit();
 
    /*
      STEP 7: Compute the thrust associated to the selected trial axis
                                Σ(p_i · vec{n})
                  thrust = max(       ÷           )
                                Σ(||p_i||) * ||n||
    */
    double trial_thrust(0.);
    for (auto const& momentum : momenta)
      trial_thrust += std::fabs(momentum.Dot(trial_axis));
 
    trial_thrust /= (sum_magnitude_mom * trial_mag);
    // trial_thrust /= sum_magnitude_mom;
 
    /*
      STEP 8: Keep trial axis as thrust axis better
    */
    if (trial_thrust > thrust) {
      thrust = trial_thrust;
      thrust_axis = trial_axis / trial_mag;
      // thrust_axis = trial_axis;
    }
 
    /*
      STEP 9: Repeat from STEP 3 to STEP 8 (select base axis)
              to avoid falling in local maxima
    */
  }
 
  /*
    STEP 10: Multiply normalized thrust axis by Thrust (<= 1)
  */
  thrust_axis *= thrust;
  return thrust_axis;
}

---

The documentation for this class was generated from the following files:

• analysis/ContinuumSuppression/include/Thrust.h
• analysis/ContinuumSuppression/src/Thrust.cc