Fork me on GitHub

source: git/external/fastjet/contribs/Nsubjettiness/AxesFinder.hh@ 5babc8b

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 5babc8b was 9687203, checked in by mselvaggi <mselvaggi@…>, 11 years ago

added nsubjettiness

  • Property mode set to 100644
File size: 16.2 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//----------------------------------------------------------------------
8// This file is part of FastJet contrib.
9//
10// It is free software; you can redistribute it and/or modify it under
11// the terms of the GNU General Public License as published by the
12// Free Software Foundation; either version 2 of the License, or (at
13// your option) any later version.
14//
15// It is distributed in the hope that it will be useful, but WITHOUT
16// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
18// License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with this code. If not, see <http://www.gnu.org/licenses/>.
22//----------------------------------------------------------------------
23
24
25#ifndef __FASTJET_CONTRIB_AXESFINDER_HH__
26#define __FASTJET_CONTRIB_AXESFINDER_HH__
27
28#include "WinnerTakeAllRecombiner.hh"
29#include "MeasureFunction.hh"
30
31#include "fastjet/PseudoJet.hh"
32#include "fastjet/ClusterSequence.hh"
33#include "fastjet/JetDefinition.hh"
34
35#include <cmath>
36#include <vector>
37#include <list>
38
39FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40
41namespace contrib{
42
43///////
44//
45// Axes Finder Options
46//
47///////
48
49//------------------------------------------------------------------------
50/// \class AxesFinder
51// This is the base class for all axes finders. These axes are used along with the MeasureFunctions to calculate
52// tau_N. There are different implementations of axes finding that are defined in derived classes below.
53class AxesFinder {
54
55protected:
56 AxesFinder* _startingFinder; // storing a possible starting finder if needed
57 std::vector<fastjet::PseudoJet> _seedAxes;
58
59 AxesFinder(AxesFinder* startingFinder = NULL) : _startingFinder(startingFinder) {}
60
61public:
62 virtual ~AxesFinder(){
63 if (_startingFinder) delete _startingFinder; //TODO: Convert to smart pointers to avoid this.
64 }
65
66 // Allow setting of seedAxes from a starting finder
67 std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
68 if (_startingFinder) {
69 _seedAxes = _startingFinder->getAxes(n_jets,inputs,currentAxes);
70 return getBetterAxes(n_jets,inputs,_seedAxes);
71 } else {
72 _seedAxes = getBetterAxes(n_jets,inputs,currentAxes);
73 return _seedAxes;
74 }
75 }
76
77 // say what the current seed axes are
78 std::vector<fastjet::PseudoJet> seedAxes() const {
79 return _seedAxes;
80 }
81
82 // This function should be overloaded, and updates the seedAxes
83 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& seedAxes) = 0;
84
85};
86
87//------------------------------------------------------------------------
88/// \class AxesFinderFromExclusiveJetDefinition
89// This class finds axes by clustering the particles and then finding the exclusive jets. This can be implemented
90// with different jet algorithms.
91class AxesFinderFromExclusiveJetDefinition : public AxesFinder {
92
93 private:
94 fastjet::JetDefinition _def;
95
96 public:
97 AxesFinderFromExclusiveJetDefinition(fastjet::JetDefinition def) : _def(def) {}
98
99 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
100 fastjet::ClusterSequence jet_clust_seq(inputs, _def);
101 return jet_clust_seq.exclusive_jets(n_jets);
102 }
103};
104
105//------------------------------------------------------------------------
106/// \class AxesFinderFromWTA_KT
107// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
108// winner take all recombination scheme.
109class AxesFinderFromWTA_KT : public AxesFinderFromExclusiveJetDefinition {
110 private:
111 const WinnerTakeAllRecombiner *recomb;
112 public:
113 AxesFinderFromWTA_KT() : AxesFinderFromExclusiveJetDefinition(
114 fastjet::JetDefinition(fastjet::kt_algorithm,
115 fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
116 recomb = new WinnerTakeAllRecombiner(),
117 fastjet::Best)) {}
118 ~AxesFinderFromWTA_KT() {delete recomb;}
119 };
120
121//------------------------------------------------------------------------
122/// \class AxesFinderFromWTA_CA
123// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
124// winner take all recombination scheme.
125class AxesFinderFromWTA_CA : public AxesFinderFromExclusiveJetDefinition {
126 private:
127 const WinnerTakeAllRecombiner *recomb;
128 public:
129 AxesFinderFromWTA_CA() : AxesFinderFromExclusiveJetDefinition(
130 fastjet::JetDefinition(fastjet::cambridge_algorithm,
131 fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
132 recomb = new WinnerTakeAllRecombiner(),
133 fastjet::Best)) {}
134 ~AxesFinderFromWTA_CA() {delete recomb;}
135};
136
137// The following classes are for testing, and are commented out for initial release
138//
139////------------------------------------------------------------------------
140///// \class AxesFinderFromWTA2_KT
141//// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
142//// winner take all recombination scheme with alpha = 2.
143//class AxesFinderFromWTA2_KT : public AxesFinderFromExclusiveJetDefinition {
144// private:
145// const WinnerTakeAllRecombiner *recomb;
146// public:
147// AxesFinderFromWTA2_KT() : AxesFinderFromExclusiveJetDefinition(
148// fastjet::JetDefinition(fastjet::kt_algorithm,
149// fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
150// recomb = new WinnerTakeAllRecombiner(2), // uses alpha = 2 here
151// fastjet::Best)) {}
152// ~AxesFinderFromWTA2_KT() {delete recomb;}
153// };
154//
155////------------------------------------------------------------------------
156///// \class AxesFinderFromWTA2_CA
157//// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
158//// winner take all recombination scheme with alpha = 2.
159//class AxesFinderFromWTA2_CA : public AxesFinderFromExclusiveJetDefinition {
160// private:
161// const WinnerTakeAllRecombiner *recomb;
162// public:
163// AxesFinderFromWTA2_CA() : AxesFinderFromExclusiveJetDefinition(
164// fastjet::JetDefinition(fastjet::cambridge_algorithm,
165// fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
166// recomb = new WinnerTakeAllRecombiner(2), //uses alpha = 2 here
167// fastjet::Best)) {}
168// ~AxesFinderFromWTA2_CA() {delete recomb;}
169//};
170
171//------------------------------------------------------------------------
172/// \class AxesFinderFromKT
173// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
174// E_scheme recombination.
175class AxesFinderFromKT : public AxesFinderFromExclusiveJetDefinition {
176 public:
177 AxesFinderFromKT() : AxesFinderFromExclusiveJetDefinition(
178 fastjet::JetDefinition(fastjet::kt_algorithm,
179 fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
180 fastjet::E_scheme,
181 fastjet::Best)) {}
182};
183
184//------------------------------------------------------------------------
185/// \class AxesFinderFromCA
186// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
187// E_scheme recombination.
188class AxesFinderFromCA : public AxesFinderFromExclusiveJetDefinition {
189 public:
190 AxesFinderFromCA() : AxesFinderFromExclusiveJetDefinition(
191 fastjet::JetDefinition(fastjet::cambridge_algorithm,
192 fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
193 fastjet::E_scheme,
194 fastjet::Best)) {}
195};
196
197
198//------------------------------------------------------------------------
199/// \class AxesFinderFromHardestJetDefinition
200// This class finds axes by clustering the particles and then finding the n hardest inclusive jets.
201// This can be implemented with different jet algorithms.
202class AxesFinderFromHardestJetDefinition : public AxesFinder {
203
204 private:
205 fastjet::JetDefinition _def;
206
207 public:
208 AxesFinderFromHardestJetDefinition(fastjet::JetDefinition def) : _def(def) {}
209
210 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
211 fastjet::ClusterSequence jet_clust_seq(inputs, _def);
212 std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
213 myJets.resize(n_jets); // only keep n hardest
214 return myJets;
215 }
216};
217
218//------------------------------------------------------------------------
219/// \class AxesFinderFromAntiKT
220// This class finds axes by finding the n hardest jets after clustering the particles according
221// to an anti kT algorithm and E_scheme.
222class AxesFinderFromAntiKT : public AxesFinderFromHardestJetDefinition {
223 public:
224 AxesFinderFromAntiKT(double R0) : AxesFinderFromHardestJetDefinition(fastjet::JetDefinition(fastjet::antikt_algorithm,R0,fastjet::E_scheme,fastjet::Best)) {}
225};
226
227
228//------------------------------------------------------------------------
229/// \class AxesFinderFromUserInput
230// This class allows the user to manually define the axes.
231class AxesFinderFromUserInput : public AxesFinder {
232
233 public:
234 AxesFinderFromUserInput() {}
235
236 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
237 assert(currentAxes.size() == (unsigned int) n_jets);
238 return currentAxes;
239 }
240};
241
242//This is a helper class for the Minimum Axes Finders. It is defined later.
243class LightLikeAxis;
244
245
246//------------------------------------------------------------------------
247/// \class AxesFinderFromOnePassMinimization
248// This class defines an AxesFinder that uses Kmeans minimization, but only on a single pass.
249class AxesFinderFromOnePassMinimization : public AxesFinder {
250
251 private:
252 double _precision; // Desired precision in axes alignment
253 int _halt; // maximum number of steps per iteration
254
255 double _beta;
256 double _Rcutoff;
257
258 DefaultUnnormalizedMeasure _measureFunction;
259
260 public:
261
262 // From a startingFinder, try to minimize the unnormalized_measure
263 AxesFinderFromOnePassMinimization(AxesFinder* startingFinder, double beta, double Rcutoff)
264 : AxesFinder(startingFinder),
265 _precision(0.0001), //hard coded for now
266 _halt(1000), //hard coded for now
267 _beta(beta),
268 _Rcutoff(Rcutoff),
269 _measureFunction(beta, Rcutoff)
270 {}
271
272 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);
273
274 template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
275 const std::vector <fastjet::PseudoJet> & inputJets);
276
277 std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
278 const std::vector <fastjet::PseudoJet> & inputJets);
279
280};
281
282
283//------------------------------------------------------------------------
284/// \class AxesFinderFromKmeansMinimization
285// This class finds finds axes by using Kmeans clustering to minimizaiton N-jettiness. Given a first set of
286// starting axes, it updates n times to get as close to the global minimum as possible. This class calls OnePass many times,
287// added noise to the axes.
288class AxesFinderFromKmeansMinimization : public AxesFinder{
289
290 private:
291 int _n_iterations; // Number of iterations to run (0 for no minimization, 1 for one-pass, >>1 for global minimum)
292 double _noise_range; // noise range for random initialization
293
294 DefaultUnnormalizedMeasure _measureFunction; //function to test whether minimum is reached
295
296 AxesFinderFromOnePassMinimization _onePassFinder; //one pass finder for minimization
297
298 PseudoJet jiggle(const PseudoJet& axis);
299
300 public:
301 AxesFinderFromKmeansMinimization(AxesFinder *startingFinder, double beta, double Rcutoff, int n_iterations) :
302 AxesFinder(startingFinder),
303 _n_iterations(n_iterations),
304 _noise_range(1.0), // hard coded for the time being
305 _measureFunction(beta, Rcutoff),
306 _onePassFinder(NULL, beta, Rcutoff)
307 {}
308
309 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);
310
311};
312
313//------------------------------------------------------------------------
314/// \class AxesFinderFromGeometricMinimization
315// This class finds axes by minimizing the Lorentz dot product distance between axes and particles. Given a first set of starting axes,
316// it essentially does stable cone finxing.
317class AxesFinderFromGeometricMinimization : public AxesFinder {
318
319 private:
320 MeasureFunction* _function;
321 double _Rcutoff;
322 double _nAttempts;
323 double _accuracy;
324
325
326 public:
327 AxesFinderFromGeometricMinimization(AxesFinder* startingFinder, double beta, double Rcutoff) : AxesFinder(startingFinder), _Rcutoff(Rcutoff) {
328 if (beta != 2.0) {
329 std::cerr << "Geometric minimization is currently only defined for beta = 2.0." << std::endl;
330 exit(1);
331 }
332
333 _nAttempts = 100;
334 _accuracy = 0.000000001;
335 _function = new GeometricMeasure(beta,_Rcutoff);
336 }
337
338 ~AxesFinderFromGeometricMinimization() {
339 delete _function;
340 }
341
342 virtual std::vector<fastjet::PseudoJet> getBetterAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes);
343};
344
345//------------------------------------------------------------------------
346/// \class LightLikeAxis
347// This is a helper class for the minimum Axes Finders classes above. It creates a convenient way of defining axes
348// in order to better facilitate calculations.
349class LightLikeAxis {
350private:
351 double _rap, _phi, _weight, _mom;
352
353 double DistanceSq(double rap2, double phi2) const {
354 double rap1 = _rap;
355 double phi1 = _phi;
356
357 double distRap = rap1-rap2;
358 double distPhi = std::fabs(phi1-phi2);
359 if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
360 return sq(distRap) + sq(distPhi);
361 }
362
363 double Distance(double rap2, double phi2) const {
364 return std::sqrt(DistanceSq(rap2,phi2));
365 }
366
367
368public:
369 LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
370 LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
371 _rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}
372 double rap() const {return _rap;}
373 double phi() const {return _phi;}
374 double weight() const {return _weight;}
375 double mom() const {return _mom;}
376 void set_rap(double my_set_rap) {_rap = my_set_rap;}
377 void set_phi(double my_set_phi) {_phi = my_set_phi;}
378 void set_weight(double my_set_weight) {_weight = my_set_weight;}
379 void set_mom(double my_set_mom) {_mom = my_set_mom;}
380 void reset(double my_rap, double my_phi, double my_weight, double my_mom) {_rap=my_rap; _phi=my_phi; _weight=my_weight; _mom=my_mom;}
381
382 // return PseudoJet with information
383 fastjet::PseudoJet ConvertToPseudoJet();
384
385 double DistanceSq(const fastjet::PseudoJet& input) const {
386 return DistanceSq(input.rap(),input.phi());
387 }
388
389 double Distance(const fastjet::PseudoJet& input) const {
390 return std::sqrt(DistanceSq(input));
391 }
392
393 double DistanceSq(const LightLikeAxis& input) const {
394 return DistanceSq(input.rap(),input.phi());
395 }
396
397 double Distance(const LightLikeAxis& input) const {
398 return std::sqrt(DistanceSq(input));
399 }
400
401};
402
403} //namespace contrib
404
405FASTJET_END_NAMESPACE
406
407#endif // __FASTJET_CONTRIB_AXESFINDER_HH__
Note: See TracBrowser for help on using the repository browser.