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: AxesDefinition.cc 833 2015-07-23 14:35:23Z 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 "AxesDefinition.hh"
|
---|
26 |
|
---|
27 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
---|
28 |
|
---|
29 | namespace contrib {
|
---|
30 |
|
---|
31 | // Repeatedly calls the one pass finder to try to find global minimum
|
---|
32 | std::vector<fastjet::PseudoJet> AxesDefinition::get_multi_pass_axes(int n_jets,
|
---|
33 | const std::vector <fastjet::PseudoJet> & inputJets,
|
---|
34 | const std::vector <fastjet::PseudoJet> & seedAxes,
|
---|
35 | const MeasureDefinition* measure) const {
|
---|
36 |
|
---|
37 | assert(n_jets == (int)seedAxes.size()); //added int casting to get rid of compiler warning
|
---|
38 |
|
---|
39 | // first iteration
|
---|
40 | std::vector<fastjet::PseudoJet> bestAxes = measure->get_one_pass_axes(n_jets, inputJets, seedAxes,_nAttempts,_accuracy);
|
---|
41 |
|
---|
42 | double bestTau = measure->result(inputJets,bestAxes);
|
---|
43 |
|
---|
44 | for (int l = 1; l < _Npass; l++) { // Do minimization procedure multiple times (l = 1 to start since first iteration is done already)
|
---|
45 |
|
---|
46 | // Add noise to current best axes
|
---|
47 | std::vector< PseudoJet > noiseAxes(n_jets, PseudoJet(0,0,0,0));
|
---|
48 | for (int k = 0; k < n_jets; k++) {
|
---|
49 | noiseAxes[k] = jiggle(bestAxes[k]);
|
---|
50 | }
|
---|
51 |
|
---|
52 | std::vector<fastjet::PseudoJet> testAxes = measure->get_one_pass_axes(n_jets, inputJets, noiseAxes,_nAttempts,_accuracy);
|
---|
53 | double testTau = measure->result(inputJets,testAxes);
|
---|
54 |
|
---|
55 | if (testTau < bestTau) {
|
---|
56 | bestTau = testTau;
|
---|
57 | bestAxes = testAxes;
|
---|
58 | }
|
---|
59 | }
|
---|
60 |
|
---|
61 | return bestAxes;
|
---|
62 | }
|
---|
63 |
|
---|
64 |
|
---|
65 | // Used by multi-pass minimization to jiggle axes with noise range.
|
---|
66 | PseudoJet AxesDefinition::jiggle(const PseudoJet& axis) const {
|
---|
67 | double phi_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
|
---|
68 | double rap_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
|
---|
69 |
|
---|
70 | double new_phi = axis.phi() + phi_noise;
|
---|
71 | if (new_phi >= 2.0*M_PI) new_phi -= 2.0*M_PI;
|
---|
72 | if (new_phi <= -2.0*M_PI) new_phi += 2.0*M_PI;
|
---|
73 |
|
---|
74 | PseudoJet newAxis(0,0,0,0);
|
---|
75 | newAxis.reset_PtYPhiM(axis.perp(),axis.rap() + rap_noise,new_phi);
|
---|
76 | return newAxis;
|
---|
77 | }
|
---|
78 |
|
---|
79 |
|
---|
80 | LimitedWarning HardestJetAxes::_too_few_axes_warning;
|
---|
81 | LimitedWarning ExclusiveJetAxes::_too_few_axes_warning;
|
---|
82 | LimitedWarning ExclusiveCombinatorialJetAxes::_too_few_axes_warning;
|
---|
83 |
|
---|
84 | std::vector<fastjet::PseudoJet> Manual_Axes::get_starting_axes(int,
|
---|
85 | const std::vector<fastjet::PseudoJet>&,
|
---|
86 | const MeasureDefinition *) const {
|
---|
87 | // This is a function dummy and should never be called
|
---|
88 | assert(false);
|
---|
89 | std::vector<fastjet::PseudoJet> dummy;
|
---|
90 | return dummy;
|
---|
91 | }
|
---|
92 |
|
---|
93 | } // namespace contrib
|
---|
94 |
|
---|
95 | FASTJET_END_NAMESPACE
|
---|