Belle II Software light-2406-ragdoll
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.

72{
73
74 /* STEP 1: Initialization of Variables */
75 XYZVector thrust_axis, trial_axis, base_axis;
76 auto end = momenta.end();
77 auto begin = momenta.begin();
78 decltype(begin) itr;
79
80 double sum_magnitude_mom = 0.;
81 double thrust = 0.;
82
83 /*
84 STEP 2: Parse momenta vector to compute magnitude Σ(||p_i||)
85 */
86
87 for (auto const& momentum : momenta)
88 sum_magnitude_mom += momentum.R();
89
90 /*
91 STEP 3: For each momentum in momenta vector,
92 use momentum as initial axis to optimize
93 */
94
95 for (auto const& mom : momenta) {
96 // By convention, thrust axis in same direction as Z axis
97 trial_axis = (mom.Z() >= 0.) ? XYZVector(mom) : XYZVector(-mom);
98
99 // Normalize if magnitude != 0
100 trial_axis = trial_axis.Unit();
101
102 /*
103 STEP 4: Store the previous trial axis as a base axis and initialize
104 a new trial axis as the Z-aligned sum of the momentum vectors
105 */
106
107 itr = begin;
108 while (itr != end) {
109 base_axis = trial_axis;
110 trial_axis = XYZVector();
111
112 // Z-alignment of momenta and addition to trial axis
113 for (auto const& momentum : momenta)
114 trial_axis += (momentum.Dot(base_axis) >= 0.) ? \
115 XYZVector(momentum) : XYZVector(-momentum);
116
117 /*
118 STEP 5: Check ( p_i · trial_axis ) * ( p_i · base_axis ) < 0 ∀ p_i
119 */
120
121 for (itr = begin; itr != end; itr++)
122 if ((*itr).Dot(trial_axis) * (*itr).Dot(base_axis) < 0.) break;
123
124 /*
125 STEP 6: While condition True:
126 base_axis = trial_axis
127 repeat STEP 4 and STEP 5;
128 else:
129 go to STEP 7 (while from STEP 4 is done)
130 */
131 }
132
133 double trial_mag(trial_axis.R()); // pre-compute
134 // trial_axis = trial_axis.Unit();
135
136 /*
137 STEP 7: Compute the thrust associated to the selected trial axis
138 Σ(p_i · vec{n})
139 thrust = max( ÷ )
140 Σ(||p_i||) * ||n||
141 */
142 double trial_thrust(0.);
143 for (auto const& momentum : momenta)
144 trial_thrust += std::fabs(momentum.Dot(trial_axis));
145
146 trial_thrust /= (sum_magnitude_mom * trial_mag);
147 // trial_thrust /= sum_magnitude_mom;
148
149 /*
150 STEP 8: Keep trial axis as thrust axis better
151 */
152 if (trial_thrust > thrust) {
153 thrust = trial_thrust;
154 thrust_axis = trial_axis / trial_mag;
155 // thrust_axis = trial_axis;
156 }
157
158 /*
159 STEP 9: Repeat from STEP 3 to STEP 8 (select base axis)
160 to avoid falling in local maxima
161 */
162 }
163
164 /*
165 STEP 10: Multiply normalized thrust axis by Thrust (<= 1)
166 */
167 thrust_axis *= thrust;
168 return thrust_axis;
169}

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