Belle II Software light-2406-ragdoll
Thrust.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
10// //
11// Thrust.cc //
12// //
13// calculation of the thrust quantity //
14// from a vector of 3D momenta. //
15// Returns the thrust axes (vector). //
16// {cc,h} //
17// Brandt, S., Dahmen, H.D. Axes and //
18// scalar measures of two-jet and //
19// three-jet events. Z. Phys. C - //
20// Particles and Fields 1, 61–70 //
21// (1979). //
22// https://doi.org/10.1007/BF01450381 //
23// http://home.thep.lu.se/~torbjorn/ppp2011/lec8.pdf
24// https://home.fnal.gov/~mrenna/lu_update/lutp0613man2/node235.html
25// //
26// Last modified 2021 March 15th //
28
29
30/* ******************** NOTES ******************** *
31 * INPUT : std::vector<XYZVector> : momenta *
32 * OUTPUT: XYZVector : thrust axis *
33 * *
34 * Thrust Axis defined as: *
35 * thrust_axis = ( vec{n} ÷ ||n|| ) * thrust *
36 * *
37 * Thrust (Scalar) defined as: *
38 * Σ(p_i · vec{n}) *
39 * thrust = max( ÷ ) *
40 * Σ(||p_i||) * ||n|| *
41 * *
42 * *********************************************** */
43
44/* *************************** STEPS *************************** *
45 * Initialization: *
46 * 1. Initialize all variables *
47 * 2. Parse momenta vector to compute magnitude Σ(||p_i||) *
48 * *
49 * Pick a thrust axis and optimize (max Σ(p_i · vec{n})): *
50 * 3. Define p_k as thrust axis (initial conditions) *
51 * 4. Define p_j = Sum (± p_i) as thrust axis (trial) *
52 * - note: ± = to sign of (p_i · p_k) *
53 * 5. Check if ( p_i · p_j ) * ( p_i · p_k ) < 0 for all p_i *
54 * 6. While True, p_k = p_j and repeat 4 -> 5, else see go to 7 *
55 * 7. Compute thrust associated to selected trial axis *
56 * *
57 * Find best local Thrust (scalar) maxima possible: *
58 * 8. Keep trial axis as thrust axis if trial > base thrust *
59 * 9. Repeat from 3 -> 8 for each momentum in momenta vector *
60 * *
61 * Return thrust axis associated to best Thrust: *
62 * 10. Multiply normalized thrust axis by Thrust (<= 1) *
63 * *
64 * ************************************************************* */
65
66#include <analysis/ContinuumSuppression/Thrust.h>
67
68using namespace Belle2;
69using namespace ROOT::Math;
70
71XYZVector Thrust::calculateThrust(const std::vector<XYZVector>& momenta)
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}
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::XYZVector > &momenta)
calculates the thrust axis
Definition: Thrust.cc:71
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24