Fork me on GitHub

source: git/external/fastjet/contribs/Nsubjettiness/MeasureDefinition.cc@ 0bbf248

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0bbf248 was 973b92a, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

update FastJet library to 3.1.3 and Nsubjettiness library to 2.2.1
  • Property mode set to 100644
File size: 24.6 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: MeasureDefinition.cc 819 2015-06-12 21:23:24Z 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
26// #include "AxesRefiner.hh"
27#include "MeasureDefinition.hh"
28
29#include <iomanip>
30
31
32FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
33
34namespace contrib {
35
36///////
37//
38// Measure Function
39//
40///////
41
42
43//descriptions updated to include measure type
44std::string DefaultMeasure::description() const {
45 std::stringstream stream;
46 stream << std::fixed << std::setprecision(2)
47 << "Default Measure (should not be used directly)";
48 return stream.str();
49};
50
51std::string NormalizedMeasure::description() const {
52 std::stringstream stream;
53 stream << std::fixed << std::setprecision(2)
54 << "Normalized Measure (beta = " << _beta << ", R0 = " << _R0 << ")";
55 return stream.str();
56};
57
58std::string UnnormalizedMeasure::description() const {
59 std::stringstream stream;
60 stream << std::fixed << std::setprecision(2)
61 << "Unnormalized Measure (beta = " << _beta << ", in GeV)";
62 return stream.str();
63};
64
65
66std::string NormalizedCutoffMeasure::description() const {
67 std::stringstream stream;
68 stream << std::fixed << std::setprecision(2)
69 << "Normalized Cutoff Measure (beta = " << _beta << ", R0 = " << _R0 << ", Rcut = " << _Rcutoff << ")";
70 return stream.str();
71};
72
73std::string UnnormalizedCutoffMeasure::description() const {
74 std::stringstream stream;
75 stream << std::fixed << std::setprecision(2)
76 << "Unnormalized Cutoff Measure (beta = " << _beta << ", Rcut = " << _Rcutoff << ", in GeV)";
77 return stream.str();
78};
79
80//std::string DeprecatedGeometricMeasure::description() const {
81// std::stringstream stream;
82// stream << std::fixed << std::setprecision(2)
83// << "Deprecated Geometric Measure (beta = " << _jet_beta << ", in GeV)";
84// return stream.str();
85//};
86
87//std::string DeprecatedGeometricCutoffMeasure::description() const {
88// std::stringstream stream;
89// stream << std::fixed << std::setprecision(2)
90// << "Deprecated Geometric Cutoff Measure (beta = " << _jet_beta << ", Rcut = " << _Rcutoff << ", in GeV)";
91// return stream.str();
92//};
93
94std::string ConicalMeasure::description() const {
95 std::stringstream stream;
96 stream << std::fixed << std::setprecision(2)
97 << "Conical Measure (beta = " << _beta << ", Rcut = " << _Rcutoff << ", in GeV)";
98 return stream.str();
99};
100
101std::string OriginalGeometricMeasure::description() const {
102 std::stringstream stream;
103 stream << std::fixed << std::setprecision(2)
104 << "Original Geometric Measure (Rcut = " << _Rcutoff << ", in GeV)";
105 return stream.str();
106};
107
108std::string ModifiedGeometricMeasure::description() const {
109 std::stringstream stream;
110 stream << std::fixed << std::setprecision(2)
111 << "Modified Geometric Measure (Rcut = " << _Rcutoff << ", in GeV)";
112 return stream.str();
113};
114
115std::string ConicalGeometricMeasure::description() const {
116 std::stringstream stream;
117 stream << std::fixed << std::setprecision(2)
118 << "Conical Geometric Measure (beta = " << _jet_beta << ", gamma = " << _beam_gamma << ", Rcut = " << _Rcutoff << ", in GeV)";
119 return stream.str();
120};
121
122
123std::string XConeMeasure::description() const {
124 std::stringstream stream;
125 stream << std::fixed << std::setprecision(2)
126 << "XCone Measure (beta = " << _jet_beta << ", Rcut = " << _Rcutoff << ", in GeV)";
127 return stream.str();
128};
129
130// Return all of the necessary TauComponents for specific input particles and axes
131TauComponents MeasureDefinition::component_result(const std::vector<fastjet::PseudoJet>& particles,
132 const std::vector<fastjet::PseudoJet>& axes) const {
133
134 // first find partition
135 TauPartition partition = get_partition(particles,axes);
136
137 // then return result calculated from partition
138 return component_result_from_partition(partition,axes);
139}
140
141TauPartition MeasureDefinition::get_partition(const std::vector<fastjet::PseudoJet>& particles,
142 const std::vector<fastjet::PseudoJet>& axes) const {
143
144 TauPartition myPartition(axes.size());
145
146 // Figures out the partiting of the input particles into the various jet pieces
147 // Based on which axis the parition is closest to
148 for (unsigned i = 0; i < particles.size(); i++) {
149
150 // find minimum distance; start with beam (-1) for reference
151 int j_min = -1;
152 double minRsq;
153 if (has_beam()) minRsq = beam_distance_squared(particles[i]);
154 else minRsq = std::numeric_limits<double>::max(); // make it large value
155
156
157 // check to see which axis the particle is closest to
158 for (unsigned j = 0; j < axes.size(); j++) {
159 double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance
160
161 if (tempRsq < minRsq) {
162 minRsq = tempRsq;
163 j_min = j;
164 }
165 }
166
167 if (j_min == -1) {
168 assert(has_beam()); // should have beam for this to make sense.
169 myPartition.push_back_beam(particles[i],i);
170 } else {
171 myPartition.push_back_jet(j_min,particles[i],i);
172 }
173 }
174
175 return myPartition;
176}
177
178// Uses existing partition and calculates result
179// TODO: Can we cache this for speed up when doing area subtraction?
180TauComponents MeasureDefinition::component_result_from_partition(const TauPartition& partition,
181 const std::vector<fastjet::PseudoJet>& axes) const {
182
183 std::vector<double> jetPieces(axes.size(), 0.0);
184 double beamPiece = 0.0;
185
186 double tauDen = 0.0;
187 if (!has_denominator()) tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor
188
189 // first find jet pieces
190 for (unsigned j = 0; j < axes.size(); j++) {
191 std::vector<PseudoJet> thisPartition = partition.jet(j).constituents();
192 for (unsigned i = 0; i < thisPartition.size(); i++) {
193 jetPieces[j] += jet_numerator(thisPartition[i],axes[j]); //numerator jet piece
194 if (has_denominator()) tauDen += denominator(thisPartition[i]); // denominator
195 }
196 }
197
198 // then find beam piece
199 if (has_beam()) {
200 std::vector<PseudoJet> beamPartition = partition.beam().constituents();
201
202 for (unsigned i = 0; i < beamPartition.size(); i++) {
203 beamPiece += beam_numerator(beamPartition[i]); //numerator beam piece
204 if (has_denominator()) tauDen += denominator(beamPartition[i]); // denominator
205 }
206 }
207
208 // create jets for storage in TauComponents
209 std::vector<PseudoJet> jets = partition.jets();
210
211 return TauComponents(_tau_mode, jetPieces, beamPiece, tauDen, jets, axes);
212}
213
214// new methods added to generalize energy and angle squared for different measure types
215double DefaultMeasure::energy(const PseudoJet& jet) const {
216 double energy;
217 switch (_measure_type) {
218 case pt_R :
219 case perp_lorentz_dot :
220 energy = jet.perp();
221 break;
222 case E_theta :
223 case lorentz_dot :
224 energy = jet.e();
225 break;
226 default : {
227 assert(_measure_type == pt_R || _measure_type == E_theta || _measure_type == lorentz_dot || _measure_type == perp_lorentz_dot);
228 energy = std::numeric_limits<double>::quiet_NaN();
229 break;
230 }
231 }
232 return energy;
233}
234
235double DefaultMeasure::angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const {
236 double pseudoRsquared;
237 switch(_measure_type) {
238 case pt_R : {
239 pseudoRsquared = jet1.squared_distance(jet2);
240 break;
241 }
242 case E_theta : {
243 // doesn't seem to be a fastjet built in for this
244 double dot = jet1.px()*jet2.px() + jet1.py()*jet2.py() + jet1.pz()*jet2.pz();
245 double norm1 = sqrt(jet1.px()*jet1.px() + jet1.py()*jet1.py() + jet1.pz()*jet1.pz());
246 double norm2 = sqrt(jet2.px()*jet2.px() + jet2.py()*jet2.py() + jet2.pz()*jet2.pz());
247
248 double costheta = dot/(norm1 * norm2);
249 if (costheta > 1.0) costheta = 1.0; // Need to handle case of numerical overflow
250 double theta = acos(costheta);
251 pseudoRsquared = theta*theta;
252 break;
253 }
254 case lorentz_dot : {
255 double dotproduct = dot_product(jet1,jet2);
256 pseudoRsquared = 2.0 * dotproduct / (jet1.e() * jet2.e());
257 break;
258 }
259 case perp_lorentz_dot : {
260 PseudoJet lightJet = lightFrom(jet2); // assuming jet2 is the axis
261 double dotproduct = dot_product(jet1,lightJet);
262 pseudoRsquared = 2.0 * dotproduct / (lightJet.pt() * jet1.pt());
263 break;
264 }
265 default : {
266 assert(_measure_type == pt_R || _measure_type == E_theta || _measure_type == lorentz_dot || _measure_type == perp_lorentz_dot);
267 pseudoRsquared = std::numeric_limits<double>::quiet_NaN();
268 break;
269 }
270 }
271
272 return pseudoRsquared;
273
274}
275
276
277///////
278//
279// Axes Refining
280//
281///////
282
283// uses minimization of N-jettiness to continually update axes until convergence.
284// The function returns the axes found at the (local) minimum
285// This is the general axes refiner that can be used for a generic measure (but is
286// overwritten in the case of the conical measure and the deprecated geometric measure)
287std::vector<fastjet::PseudoJet> MeasureDefinition::get_one_pass_axes(int n_jets,
288 const std::vector <fastjet::PseudoJet> & particles,
289 const std::vector<fastjet::PseudoJet>& currentAxes,
290 int nAttempts,
291 double accuracy) const {
292
293 assert(n_jets == (int)currentAxes.size());
294
295 std::vector<fastjet::PseudoJet> seedAxes = currentAxes;
296
297 std::vector<fastjet::PseudoJet> temp_axes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0));
298 for (unsigned int k = 0; k < seedAxes.size(); k++) {
299 seedAxes[k] = lightFrom(seedAxes[k]) * seedAxes[k].E(); // making light-like, but keeping energy
300 }
301
302 double seedTau = result(particles, seedAxes);
303
304 std::vector<fastjet::PseudoJet> bestAxesSoFar = seedAxes;
305 double bestTauSoFar = seedTau;
306
307 for (int i_att = 0; i_att < nAttempts; i_att++) {
308
309 std::vector<fastjet::PseudoJet> newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0));
310 std::vector<fastjet::PseudoJet> summed_jets(seedAxes.size(), fastjet::PseudoJet(0,0,0,0));
311
312 // find closest axis and assign to that
313 for (unsigned int i = 0; i < particles.size(); i++) {
314
315 // start from unclustered beam measure
316 int minJ = -1;
317 double minDist = beam_distance_squared(particles[i]);
318
319 // which axis am I closest to?
320 for (unsigned int j = 0; j < seedAxes.size(); j++) {
321 double tempDist = jet_distance_squared(particles[i],seedAxes[j]);
322 if (tempDist < minDist) {
323 minDist = tempDist;
324 minJ = j;
325 }
326 }
327
328 // if not unclustered, then cluster
329 if (minJ != -1) {
330 summed_jets[minJ] += particles[i]; // keep track of energy to use later.
331 if (_useAxisScaling) {
332 double pseudoMomentum = dot_product(lightFrom(seedAxes[minJ]),particles[i]) + accuracy; // need small offset to avoid potential divide by zero issues
333 double axis_scaling = (double)jet_numerator(particles[i], seedAxes[minJ])/pseudoMomentum;
334
335 newAxes[minJ] += particles[i]*axis_scaling;
336 }
337 }
338 }
339 if (!_useAxisScaling) newAxes = summed_jets;
340
341 // convert the axes to LightLike and then back to PseudoJet
342 for (unsigned int k = 0; k < newAxes.size(); k++) {
343 if (newAxes[k].perp() > 0) {
344 newAxes[k] = lightFrom(newAxes[k]);
345 newAxes[k] *= summed_jets[k].E(); // scale by energy to get sensible result
346 }
347 }
348
349 // calculate tau on new axes
350 double newTau = result(particles, newAxes);
351
352 // find the smallest value of tau (and the corresponding axes) so far
353 if (newTau < bestTauSoFar) {
354 bestAxesSoFar = newAxes;
355 bestTauSoFar = newTau;
356 }
357
358 if (fabs(newTau - seedTau) < accuracy) {// close enough for jazz
359 seedAxes = newAxes;
360 seedTau = newTau;
361 break;
362 }
363
364 seedAxes = newAxes;
365 seedTau = newTau;
366
367}
368
369 // return the axes corresponding to the smallest tau found throughout all iterations
370 // this is to prevent the minimization from returning a non-minimized of tau due to potential oscillations around the minimum
371 return bestAxesSoFar;
372
373}
374
375
376// One pass minimization for the DefaultMeasure
377
378// Given starting axes, update to find better axes by using Kmeans clustering around the old axes
379template <int N>
380std::vector<LightLikeAxis> DefaultMeasure::UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
381 const std::vector <fastjet::PseudoJet> & inputJets,
382 double accuracy
383 ) const {
384 assert(old_axes.size() == N);
385
386 // some storage, declared static to save allocation/re-allocation costs
387 static LightLikeAxis new_axes[N];
388 static fastjet::PseudoJet new_jets[N];
389 for (int n = 0; n < N; ++n) {
390 new_axes[n].reset(0.0,0.0,0.0,0.0);
391 new_jets[n].reset_momentum(0.0,0.0,0.0,0.0);
392 }
393
394 double precision = accuracy; //TODO: actually cascade this in
395
396 /////////////// Assignment Step //////////////////////////////////////////////////////////
397 std::vector<int> assignment_index(inputJets.size());
398 int k_assign = -1;
399
400 for (unsigned i = 0; i < inputJets.size(); i++){
401 double smallestDist = std::numeric_limits<double>::max(); //large number
402 for (int k = 0; k < N; k++) {
403 double thisDist = old_axes[k].DistanceSq(inputJets[i]);
404 if (thisDist < smallestDist) {
405 smallestDist = thisDist;
406 k_assign = k;
407 }
408 }
409 if (smallestDist > sq(_Rcutoff)) {k_assign = -1;}
410 assignment_index[i] = k_assign;
411 }
412
413 //////////////// Update Step /////////////////////////////////////////////////////////////
414 double distPhi, old_dist;
415 for (unsigned i = 0; i < inputJets.size(); i++) {
416 int old_jet_i = assignment_index[i];
417 if (old_jet_i == -1) {continue;}
418
419 const fastjet::PseudoJet& inputJet_i = inputJets[i];
420 LightLikeAxis& new_axis_i = new_axes[old_jet_i];
421 double inputPhi_i = inputJet_i.phi();
422 double inputRap_i = inputJet_i.rap();
423
424 // optimize pow() call
425 // add noise (the precision term) to make sure we don't divide by zero
426 if (_beta == 1.0) {
427 double DR = std::sqrt(sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i));
428 old_dist = 1.0/DR;
429 } else if (_beta == 2.0) {
430 old_dist = 1.0;
431 } else if (_beta == 0.0) {
432 double DRSq = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
433 old_dist = 1.0/DRSq;
434 } else {
435 old_dist = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
436 old_dist = std::pow(old_dist, (0.5*_beta-1.0));
437 }
438
439 // TODO: Put some of these addition functions into light-like axes
440 // rapidity sum
441 new_axis_i.set_rap(new_axis_i.rap() + inputJet_i.perp() * inputRap_i * old_dist);
442 // phi sum
443 distPhi = inputPhi_i - old_axes[old_jet_i].phi();
444 if (fabs(distPhi) <= M_PI){
445 new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * inputPhi_i * old_dist );
446 } else if (distPhi > M_PI) {
447 new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (-2*M_PI + inputPhi_i) * old_dist );
448 } else if (distPhi < -M_PI) {
449 new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (+2*M_PI + inputPhi_i) * old_dist );
450 }
451 // weights sum
452 new_axis_i.set_weight( new_axis_i.weight() + inputJet_i.perp() * old_dist );
453 // momentum magnitude sum
454 new_jets[old_jet_i] += inputJet_i;
455 }
456 // normalize sums
457 for (int k = 0; k < N; k++) {
458 if (new_axes[k].weight() == 0) {
459 // no particles were closest to this axis! Return to old axis instead of (0,0,0,0)
460 new_axes[k] = old_axes[k];
461 } else {
462 new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() );
463 new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() );
464 new_axes[k].set_phi( std::fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) );
465 new_axes[k].set_mom( std::sqrt(new_jets[k].modp2()) );
466 }
467 }
468 std::vector<LightLikeAxis> new_axes_vec(N);
469 for (unsigned k = 0; k < N; ++k) new_axes_vec[k] = new_axes[k];
470 return new_axes_vec;
471}
472
473// Given N starting axes, this function updates all axes to find N better axes.
474// (This is just a wrapper for the templated version above.)
475// TODO: Consider removing this in a future version
476std::vector<LightLikeAxis> DefaultMeasure::UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
477 const std::vector <fastjet::PseudoJet> & inputJets,
478 double accuracy) const {
479 int N = old_axes.size();
480 switch (N) {
481 case 1: return UpdateAxesFast<1>(old_axes, inputJets, accuracy);
482 case 2: return UpdateAxesFast<2>(old_axes, inputJets, accuracy);
483 case 3: return UpdateAxesFast<3>(old_axes, inputJets, accuracy);
484 case 4: return UpdateAxesFast<4>(old_axes, inputJets, accuracy);
485 case 5: return UpdateAxesFast<5>(old_axes, inputJets, accuracy);
486 case 6: return UpdateAxesFast<6>(old_axes, inputJets, accuracy);
487 case 7: return UpdateAxesFast<7>(old_axes, inputJets, accuracy);
488 case 8: return UpdateAxesFast<8>(old_axes, inputJets, accuracy);
489 case 9: return UpdateAxesFast<9>(old_axes, inputJets, accuracy);
490 case 10: return UpdateAxesFast<10>(old_axes, inputJets, accuracy);
491 case 11: return UpdateAxesFast<11>(old_axes, inputJets, accuracy);
492 case 12: return UpdateAxesFast<12>(old_axes, inputJets, accuracy);
493 case 13: return UpdateAxesFast<13>(old_axes, inputJets, accuracy);
494 case 14: return UpdateAxesFast<14>(old_axes, inputJets, accuracy);
495 case 15: return UpdateAxesFast<15>(old_axes, inputJets, accuracy);
496 case 16: return UpdateAxesFast<16>(old_axes, inputJets, accuracy);
497 case 17: return UpdateAxesFast<17>(old_axes, inputJets, accuracy);
498 case 18: return UpdateAxesFast<18>(old_axes, inputJets, accuracy);
499 case 19: return UpdateAxesFast<19>(old_axes, inputJets, accuracy);
500 case 20: return UpdateAxesFast<20>(old_axes, inputJets, accuracy);
501 default: std::cout << "N-jettiness is hard-coded to only allow up to 20 jets!" << std::endl;
502 return std::vector<LightLikeAxis>();
503 }
504
505}
506
507// uses minimization of N-jettiness to continually update axes until convergence.
508// The function returns the axes found at the (local) minimum
509std::vector<fastjet::PseudoJet> DefaultMeasure::get_one_pass_axes(int n_jets,
510 const std::vector <fastjet::PseudoJet> & inputJets,
511 const std::vector<fastjet::PseudoJet>& seedAxes,
512 int nAttempts,
513 double accuracy
514 ) const {
515
516 // if the measure type doesn't use the pt_R metric, then the standard minimization scheme should be used
517 if (_measure_type != pt_R) {
518 return MeasureDefinition::get_one_pass_axes(n_jets, inputJets, seedAxes, nAttempts, accuracy);
519 }
520
521 // convert from PseudoJets to LightLikeAxes
522 std::vector< LightLikeAxis > old_axes(n_jets, LightLikeAxis(0,0,0,0));
523 for (int k = 0; k < n_jets; k++) {
524 old_axes[k].set_rap( seedAxes[k].rap() );
525 old_axes[k].set_phi( seedAxes[k].phi() );
526 }
527
528 // Find new axes by iterating (only one pass here)
529 std::vector< LightLikeAxis > new_axes(n_jets, LightLikeAxis(0,0,0,0));
530 double cmp = std::numeric_limits<double>::max(); //large number
531 int h = 0;
532
533 while (cmp > accuracy && h < nAttempts) { // Keep updating axes until near-convergence or too many update steps
534 cmp = 0.0;
535 h++;
536 new_axes = UpdateAxes(old_axes, inputJets,accuracy); // Update axes
537 for (int k = 0; k < n_jets; k++) {
538 cmp += old_axes[k].Distance(new_axes[k]);
539 }
540 cmp = cmp / ((double) n_jets);
541 old_axes = new_axes;
542 }
543
544 // Convert from internal LightLikeAxes to PseudoJet
545 std::vector<fastjet::PseudoJet> outputAxes;
546 for (int k = 0; k < n_jets; k++) {
547 fastjet::PseudoJet temp = old_axes[k].ConvertToPseudoJet();
548 outputAxes.push_back(temp);
549 }
550
551 // this is used to debug the minimization routine to make sure that it works.
552 bool do_debug = false;
553 if (do_debug) {
554 // get this information to make sure that minimization is working properly
555 double seed_tau = result(inputJets, seedAxes);
556 double outputTau = result(inputJets, outputAxes);
557 assert(outputTau <= seed_tau);
558 }
559
560 return outputAxes;
561}
562
563//// One-pass minimization for the Deprecated Geometric Measure
564//// Uses minimization of the geometric distance in order to find the minimum axes.
565//// It continually updates until it reaches convergence or it reaches the maximum number of attempts.
566//// This is essentially the same as a stable cone finder.
567//std::vector<fastjet::PseudoJet> DeprecatedGeometricCutoffMeasure::get_one_pass_axes(int n_jets,
568// const std::vector <fastjet::PseudoJet> & particles,
569// const std::vector<fastjet::PseudoJet>& currentAxes,
570// int nAttempts,
571// double accuracy) const {
572//
573// assert(n_jets == (int)currentAxes.size()); //added int casting to get rid of compiler warning
574//
575// std::vector<fastjet::PseudoJet> seedAxes = currentAxes;
576// double seedTau = result(particles, seedAxes);
577//
578// for (int i = 0; i < nAttempts; i++) {
579//
580// std::vector<fastjet::PseudoJet> newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0));
581//
582// // find closest axis and assign to that
583// for (unsigned int i = 0; i < particles.size(); i++) {
584//
585// // start from unclustered beam measure
586// int minJ = -1;
587// double minDist = beam_distance_squared(particles[i]);
588//
589// // which axis am I closest to?
590// for (unsigned int j = 0; j < seedAxes.size(); j++) {
591// double tempDist = jet_distance_squared(particles[i],seedAxes[j]);
592// if (tempDist < minDist) {
593// minDist = tempDist;
594// minJ = j;
595// }
596// }
597//
598// // if not unclustered, then cluster
599// if (minJ != -1) newAxes[minJ] += particles[i];
600// }
601//
602// // calculate tau on new axes
603// seedAxes = newAxes;
604// double tempTau = result(particles, newAxes);
605//
606// // close enough to stop?
607// if (fabs(tempTau - seedTau) < accuracy) break;
608// seedTau = tempTau;
609// }
610//
611// return seedAxes;
612//}
613
614
615// Go from internal LightLikeAxis to PseudoJet
616fastjet::PseudoJet LightLikeAxis::ConvertToPseudoJet() {
617 double px, py, pz, E;
618 E = _mom;
619 pz = (std::exp(2.0*_rap) - 1.0) / (std::exp(2.0*_rap) + 1.0) * E;
620 px = std::cos(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
621 py = std::sin(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
622 return fastjet::PseudoJet(px,py,pz,E);
623}
624
625} //namespace contrib
626
627FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.