Fork me on GitHub

source: git/external/fastjet/contribs/Nsubjettiness/AxesFinder.cc@ 49234af

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 49234af was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 12.5 KB
Line 
1// Nsubjettiness Package
2// Questions/Comments? jthaler@jthaler.net
3//
4// Copyright (c) 2011-14
5// Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
6//
7// $Id: AxesFinder.cc 670 2014-06-06 01:24:42Z jthaler $
8//----------------------------------------------------------------------
9// This file is part of FastJet contrib.
10//
11// It is free software; you can redistribute it and/or modify it under
12// the terms of the GNU General Public License as published by the
13// Free Software Foundation; either version 2 of the License, or (at
14// your option) any later version.
15//
16// It is distributed in the hope that it will be useful, but WITHOUT
17// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19// License for more details.
20//
21// You should have received a copy of the GNU General Public License
22// along with this code. If not, see <http://www.gnu.org/licenses/>.
23//----------------------------------------------------------------------
24
25#include "AxesFinder.hh"
26
27FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
28
29namespace contrib{
30
31///////
32//
33// Functions for minimization.
34//
35///////
36
37// Given starting axes, update to find better axes by using Kmeans clustering around the old axes
38template <int N>
39std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
40 const std::vector <fastjet::PseudoJet> & inputJets) const {
41 assert(old_axes.size() == N);
42
43 // some storage, declared static to save allocation/re-allocation costs
44 static LightLikeAxis new_axes[N];
45 static fastjet::PseudoJet new_jets[N];
46 for (int n = 0; n < N; ++n) {
47 new_axes[n].reset(0.0,0.0,0.0,0.0);
48 new_jets[n].reset_momentum(0.0,0.0,0.0,0.0);
49 }
50
51 double precision = _precision;
52
53 /////////////// Assignment Step //////////////////////////////////////////////////////////
54 std::vector<int> assignment_index(inputJets.size());
55 int k_assign = -1;
56
57 for (unsigned i = 0; i < inputJets.size(); i++){
58 double smallestDist = std::numeric_limits<double>::max(); //large number
59 for (int k = 0; k < N; k++) {
60 double thisDist = old_axes[k].DistanceSq(inputJets[i]);
61 if (thisDist < smallestDist) {
62 smallestDist = thisDist;
63 k_assign = k;
64 }
65 }
66 if (smallestDist > sq(_Rcutoff)) {k_assign = -1;}
67 assignment_index[i] = k_assign;
68 }
69
70 //////////////// Update Step /////////////////////////////////////////////////////////////
71 double distPhi, old_dist;
72 for (unsigned i = 0; i < inputJets.size(); i++) {
73 int old_jet_i = assignment_index[i];
74 if (old_jet_i == -1) {continue;}
75
76 const fastjet::PseudoJet& inputJet_i = inputJets[i];
77 LightLikeAxis& new_axis_i = new_axes[old_jet_i];
78 double inputPhi_i = inputJet_i.phi();
79 double inputRap_i = inputJet_i.rap();
80
81 // optimize pow() call
82 // add noise (the precision term) to make sure we don't divide by zero
83 if (_beta == 1.0) {
84 double DR = std::sqrt(sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i));
85 old_dist = 1.0/DR;
86 } else if (_beta == 2.0) {
87 old_dist = 1.0;
88 } else if (_beta == 0.0) {
89 double DRSq = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
90 old_dist = 1.0/DRSq;
91 } else {
92 old_dist = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
93 old_dist = std::pow(old_dist, (0.5*_beta-1.0));
94 }
95
96 // TODO: Put some of these addition functions into light-like axes
97 // rapidity sum
98 new_axis_i.set_rap(new_axis_i.rap() + inputJet_i.perp() * inputRap_i * old_dist);
99 // phi sum
100 distPhi = inputPhi_i - old_axes[old_jet_i].phi();
101 if (fabs(distPhi) <= M_PI){
102 new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * inputPhi_i * old_dist );
103 } else if (distPhi > M_PI) {
104 new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (-2*M_PI + inputPhi_i) * old_dist );
105 } else if (distPhi < -M_PI) {
106 new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (+2*M_PI + inputPhi_i) * old_dist );
107 }
108 // weights sum
109 new_axis_i.set_weight( new_axis_i.weight() + inputJet_i.perp() * old_dist );
110 // momentum magnitude sum
111 new_jets[old_jet_i] += inputJet_i;
112 }
113 // normalize sums
114 for (int k = 0; k < N; k++) {
115 if (new_axes[k].weight() == 0) {
116 // no particles were closest to this axis! Return to old axis instead of (0,0,0,0)
117 new_axes[k] = old_axes[k];
118 } else {
119 new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() );
120 new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() );
121 new_axes[k].set_phi( std::fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) );
122 new_axes[k].set_mom( std::sqrt(new_jets[k].modp2()) );
123 }
124 }
125 std::vector<LightLikeAxis> new_axes_vec(N);
126 for (unsigned k = 0; k < N; ++k) new_axes_vec[k] = new_axes[k];
127 return new_axes_vec;
128}
129
130// Given N starting axes, this function updates all axes to find N better axes.
131// (This is just a wrapper for the templated version above.)
132std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
133 const std::vector <fastjet::PseudoJet> & inputJets) const {
134 int N = old_axes.size();
135 switch (N) {
136 case 1: return UpdateAxesFast<1>(old_axes, inputJets);
137 case 2: return UpdateAxesFast<2>(old_axes, inputJets);
138 case 3: return UpdateAxesFast<3>(old_axes, inputJets);
139 case 4: return UpdateAxesFast<4>(old_axes, inputJets);
140 case 5: return UpdateAxesFast<5>(old_axes, inputJets);
141 case 6: return UpdateAxesFast<6>(old_axes, inputJets);
142 case 7: return UpdateAxesFast<7>(old_axes, inputJets);
143 case 8: return UpdateAxesFast<8>(old_axes, inputJets);
144 case 9: return UpdateAxesFast<9>(old_axes, inputJets);
145 case 10: return UpdateAxesFast<10>(old_axes, inputJets);
146 case 11: return UpdateAxesFast<11>(old_axes, inputJets);
147 case 12: return UpdateAxesFast<12>(old_axes, inputJets);
148 case 13: return UpdateAxesFast<13>(old_axes, inputJets);
149 case 14: return UpdateAxesFast<14>(old_axes, inputJets);
150 case 15: return UpdateAxesFast<15>(old_axes, inputJets);
151 case 16: return UpdateAxesFast<16>(old_axes, inputJets);
152 case 17: return UpdateAxesFast<17>(old_axes, inputJets);
153 case 18: return UpdateAxesFast<18>(old_axes, inputJets);
154 case 19: return UpdateAxesFast<19>(old_axes, inputJets);
155 case 20: return UpdateAxesFast<20>(old_axes, inputJets);
156 default: std::cout << "N-jettiness is hard-coded to only allow up to 20 jets!" << std::endl;
157 return std::vector<LightLikeAxis>();
158 }
159
160}
161
162// uses minimization of N-jettiness to continually update axes until convergence.
163// The function returns the axes found at the (local) minimum
164std::vector<fastjet::PseudoJet> AxesFinderFromOnePassMinimization::getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) const {
165
166 // convert from PseudoJets to LightLikeAxes
167 std::vector< LightLikeAxis > old_axes(n_jets, LightLikeAxis(0,0,0,0));
168 for (int k = 0; k < n_jets; k++) {
169 old_axes[k].set_rap( seedAxes[k].rap() );
170 old_axes[k].set_phi( seedAxes[k].phi() );
171 }
172
173 // Find new axes by iterating (only one pass here)
174 std::vector< LightLikeAxis > new_axes(n_jets, LightLikeAxis(0,0,0,0));
175 double cmp = std::numeric_limits<double>::max(); //large number
176 int h = 0;
177 while (cmp > _precision && h < _halt) { // Keep updating axes until near-convergence or too many update steps
178 cmp = 0.0;
179 h++;
180 new_axes = UpdateAxes(old_axes, inputJets); // Update axes
181 for (int k = 0; k < n_jets; k++) {
182 cmp += old_axes[k].Distance(new_axes[k]);
183 }
184 cmp = cmp / ((double) n_jets);
185 old_axes = new_axes;
186 }
187
188 // Convert from internal LightLikeAxes to PseudoJet
189 std::vector<fastjet::PseudoJet> outputAxes;
190 for (int k = 0; k < n_jets; k++) {
191 fastjet::PseudoJet temp = old_axes[k].ConvertToPseudoJet();
192 outputAxes.push_back(temp);
193 }
194
195 // this is used to debug the minimization routine to make sure that it works.
196 bool do_debug = false;
197 if (do_debug) {
198 // get this information to make sure that minimization is working properly
199 TauComponents seed_tau_components = _measureFunction.result(inputJets, seedAxes);
200 double seed_tau = seed_tau_components.tau();
201 TauComponents tau_components = _measureFunction.result(inputJets, outputAxes);
202 double outputTau = tau_components.tau();
203 assert(outputTau <= seed_tau);
204 }
205
206 return outputAxes;
207}
208
209PseudoJet AxesFinderFromKmeansMinimization::jiggle(const PseudoJet& axis) const {
210 double phi_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
211 double rap_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
212
213 double new_phi = axis.phi() + phi_noise;
214 if (new_phi >= 2.0*M_PI) new_phi -= 2.0*M_PI;
215 if (new_phi <= -2.0*M_PI) new_phi += 2.0*M_PI;
216
217 PseudoJet newAxis(0,0,0,0);
218 newAxis.reset_PtYPhiM(axis.perp(),axis.rap() + rap_noise,new_phi);
219 return newAxis;
220}
221
222
223// Repeatedly calls the one pass finder to try to find global minimum
224std::vector<fastjet::PseudoJet> AxesFinderFromKmeansMinimization::getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) const {
225
226 // first iteration
227 std::vector<fastjet::PseudoJet> bestAxes = _onePassFinder.getAxes(n_jets, inputJets, seedAxes);
228
229 double bestTau = (_measureFunction.result(inputJets,bestAxes)).tau();
230
231 for (int l = 1; l < _n_iterations; l++) { // Do minimization procedure multiple times (l = 1 to start since first iteration is done already)
232
233 // Add noise to current best axes
234 std::vector< PseudoJet > noiseAxes(n_jets, PseudoJet(0,0,0,0));
235 for (int k = 0; k < n_jets; k++) {
236 noiseAxes[k] = jiggle(bestAxes[k]);
237 }
238
239 std::vector<fastjet::PseudoJet> testAxes = _onePassFinder.getAxes(n_jets, inputJets, noiseAxes);
240 double testTau = (_measureFunction.result(inputJets,testAxes)).tau();
241
242 if (testTau < bestTau) {
243 bestTau = testTau;
244 bestAxes = testAxes;
245 }
246 }
247
248 return bestAxes;
249}
250
251// Uses minimization of the geometric distance in order to find the minimum axes.
252// It continually updates until it reaches convergence or it reaches the maximum number of attempts.
253// This is essentially the same as a stable cone finder.
254std::vector<fastjet::PseudoJet> AxesFinderFromGeometricMinimization::getAxes(int /*n_jets*/, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes) const {
255
256 std::vector<fastjet::PseudoJet> seedAxes = currentAxes;
257 double seedTau = _function.tau(particles, seedAxes);
258
259 for (int i = 0; i < _nAttempts; i++) {
260
261 std::vector<fastjet::PseudoJet> newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0));
262
263 // find closest axis and assign to that
264 for (unsigned int i = 0; i < particles.size(); i++) {
265
266 // start from unclustered beam measure
267 int minJ = -1;
268 double minDist = _function.beam_distance_squared(particles[i]);
269
270 // which axis am I closest to?
271 for (unsigned int j = 0; j < seedAxes.size(); j++) {
272 double tempDist = _function.jet_distance_squared(particles[i],seedAxes[j]);
273 if (tempDist < minDist) {
274 minDist = tempDist;
275 minJ = j;
276 }
277 }
278
279 // if not unclustered, then cluster
280 if (minJ != -1) newAxes[minJ] += particles[i];
281 }
282
283 // calculate tau on new axes
284 seedAxes = newAxes;
285 double tempTau = _function.tau(particles, newAxes);
286
287 // close enough to stop?
288 if (fabs(tempTau - seedTau) < _accuracy) break;
289 seedTau = tempTau;
290 }
291
292 return seedAxes;
293}
294
295// Go from internal LightLikeAxis to PseudoJet
296fastjet::PseudoJet LightLikeAxis::ConvertToPseudoJet() {
297 double px, py, pz, E;
298 E = _mom;
299 pz = (std::exp(2.0*_rap) - 1.0) / (std::exp(2.0*_rap) + 1.0) * E;
300 px = std::cos(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
301 py = std::sin(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
302 return fastjet::PseudoJet(px,py,pz,E);
303}
304
305} //namespace contrib
306
307FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.