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 | #include "Njettiness.hh"
|
---|
25 |
|
---|
26 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
---|
27 |
|
---|
28 | namespace contrib {
|
---|
29 |
|
---|
30 |
|
---|
31 | ///////
|
---|
32 | //
|
---|
33 | // Main Njettiness Class
|
---|
34 | //
|
---|
35 | ///////
|
---|
36 |
|
---|
37 | // Helper function to correlate one pass minimization with appropriate measure
|
---|
38 | void Njettiness::setOnePassAxesFinder(MeasureMode measure_mode, AxesFinder* startingFinder, double beta, double Rcutoff) {
|
---|
39 | if (measure_mode == normalized_measure || measure_mode == unnormalized_measure || measure_mode == normalized_cutoff_measure || measure_mode == unnormalized_cutoff_measure) {
|
---|
40 | _axesFinder = new AxesFinderFromOnePassMinimization(startingFinder, beta, Rcutoff);
|
---|
41 | }
|
---|
42 | else if (measure_mode == geometric_measure || measure_mode == geometric_cutoff_measure) {
|
---|
43 | _axesFinder = new AxesFinderFromGeometricMinimization(startingFinder, beta, Rcutoff);
|
---|
44 | }
|
---|
45 | else {
|
---|
46 | std::cerr << "Minimization only set up for normalized_measure, unnormalized_measure, normalized_cutoff_measure, unnormalized_cutoff_measure, geometric_measure, geometric_cutoff_measure" << std::endl;
|
---|
47 | exit(1); }
|
---|
48 | }
|
---|
49 |
|
---|
50 | // Parsing needed for constructor to set AxesFinder and MeasureFunction
|
---|
51 | // All of the parameter handling is here, and checking that number of parameters is correct.
|
---|
52 | void Njettiness::setMeasureFunctionandAxesFinder(AxesMode axes_mode, MeasureMode measure_mode, double para1, double para2, double para3, double para4) {
|
---|
53 |
|
---|
54 | // definition of maximum Rcutoff for non-cutoff measures, changed later by other measures
|
---|
55 | double Rcutoff = std::numeric_limits<double>::max(); //large number
|
---|
56 | // Most (but all measures have some kind of beta value)
|
---|
57 | double beta = NAN;
|
---|
58 | // The normalized measures have an R0 value.
|
---|
59 | double R0 = NAN;
|
---|
60 |
|
---|
61 | // Find the MeasureFunction and set the parameters.
|
---|
62 | switch (measure_mode) {
|
---|
63 | case normalized_measure:
|
---|
64 | beta = para1;
|
---|
65 | R0 = para2;
|
---|
66 | if(correctParameterCount(2, para1, para2, para3, para4))
|
---|
67 | _measureFunction = new DefaultNormalizedMeasure(beta, R0, Rcutoff); //normalized_measure requires 2 parameters, beta and R0
|
---|
68 | else {
|
---|
69 | std::cerr << "normalized_measure needs 2 parameters (beta and R0)" << std::endl;
|
---|
70 | exit(1); }
|
---|
71 | break;
|
---|
72 | case unnormalized_measure:
|
---|
73 | beta = para1;
|
---|
74 | if(correctParameterCount(1, para1, para2, para3, para4))
|
---|
75 | _measureFunction = new DefaultUnnormalizedMeasure(beta, Rcutoff); //unnormalized_measure requires 1 parameter, beta
|
---|
76 | else {
|
---|
77 | std::cerr << "unnormalized_measure needs 1 parameter (beta)" << std::endl;
|
---|
78 | exit(1); }
|
---|
79 | break;
|
---|
80 | case geometric_measure:
|
---|
81 | beta = para1;
|
---|
82 | if(correctParameterCount(1, para1, para2, para3, para4))
|
---|
83 | _measureFunction = new GeometricMeasure(beta,Rcutoff); //geometric_measure requires 1 parameter, beta
|
---|
84 | else {
|
---|
85 | std::cerr << "geometric_measure needs 1 parameter (beta)" << std::endl;
|
---|
86 | exit(1); }
|
---|
87 | break;
|
---|
88 | case normalized_cutoff_measure:
|
---|
89 | beta = para1;
|
---|
90 | R0 = para2;
|
---|
91 | Rcutoff = para3; //Rcutoff parameter is 3rd parameter in normalized_cutoff_measure
|
---|
92 | if(correctParameterCount(3, para1, para2, para3, para4))
|
---|
93 | _measureFunction = new DefaultNormalizedMeasure(beta, R0, Rcutoff); //normalized_cutoff_measure requires 3 parameters, beta, R0, and Rcutoff
|
---|
94 | else {
|
---|
95 | std::cerr << "normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)" << std::endl;
|
---|
96 | exit(1); }
|
---|
97 | break;
|
---|
98 | case unnormalized_cutoff_measure:
|
---|
99 | beta = para1;
|
---|
100 | Rcutoff = para2; //Rcutoff parameter is 2nd parameter in normalized_cutoff_measure
|
---|
101 | if (correctParameterCount(2, para1, para2, para3, para4))
|
---|
102 | _measureFunction = new DefaultUnnormalizedMeasure(beta, Rcutoff); //unnormalized_cutoff_measure requires 2 parameters, beta and Rcutoff
|
---|
103 | else {
|
---|
104 | std::cerr << "unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)" << std::endl;
|
---|
105 | exit(1); }
|
---|
106 | break;
|
---|
107 | case geometric_cutoff_measure:
|
---|
108 | beta = para1;
|
---|
109 | Rcutoff = para2; //Rcutoff parameter is 2nd parameter in geometric_cutoff_measure
|
---|
110 | if(correctParameterCount(2, para1, para2, para3, para4))
|
---|
111 | _measureFunction = new GeometricMeasure(beta,Rcutoff); //geometric_cutoff_measure requires 2 parameters, beta and Rcutoff
|
---|
112 | else {
|
---|
113 | std::cerr << "geometric_cutoff_measure has 2 parameters (beta,Rcutoff)" << std::endl;
|
---|
114 | exit(1); }
|
---|
115 | break;
|
---|
116 | default:
|
---|
117 | assert(false);
|
---|
118 | break;
|
---|
119 | }
|
---|
120 |
|
---|
121 | // Choose which AxesFinder from user input.
|
---|
122 | // Uses setOnePassAxesFinder helpful function to use beta and Rcutoff values about (if needed)
|
---|
123 | switch (axes_mode) {
|
---|
124 | case wta_kt_axes:
|
---|
125 | _axesFinder = new AxesFinderFromWTA_KT();
|
---|
126 | break;
|
---|
127 | case wta_ca_axes:
|
---|
128 | _axesFinder = new AxesFinderFromWTA_CA();
|
---|
129 | break;
|
---|
130 | case kt_axes:
|
---|
131 | _axesFinder = new AxesFinderFromKT();
|
---|
132 | break;
|
---|
133 | case ca_axes:
|
---|
134 | _axesFinder = new AxesFinderFromCA();
|
---|
135 | break;
|
---|
136 | case antikt_0p2_axes:
|
---|
137 | _axesFinder = new AxesFinderFromAntiKT(0.2);
|
---|
138 | break;
|
---|
139 | case onepass_wta_kt_axes:
|
---|
140 | setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_KT(), beta, Rcutoff);
|
---|
141 | break;
|
---|
142 | case onepass_wta_ca_axes:
|
---|
143 | setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA_CA(), beta, Rcutoff);
|
---|
144 | break;
|
---|
145 | case onepass_kt_axes:
|
---|
146 | setOnePassAxesFinder(measure_mode, new AxesFinderFromKT(), beta, Rcutoff);
|
---|
147 | break;
|
---|
148 | case onepass_ca_axes:
|
---|
149 | setOnePassAxesFinder(measure_mode, new AxesFinderFromCA(), beta, Rcutoff);
|
---|
150 | break;
|
---|
151 | case onepass_antikt_0p2_axes:
|
---|
152 | setOnePassAxesFinder(measure_mode, new AxesFinderFromAntiKT(0.2), beta, Rcutoff);
|
---|
153 | break;
|
---|
154 | case onepass_manual_axes:
|
---|
155 | setOnePassAxesFinder(measure_mode, new AxesFinderFromUserInput(), beta, Rcutoff);
|
---|
156 | break;
|
---|
157 | case min_axes: //full minimization is not defined for geometric_measure.
|
---|
158 | if (measure_mode == normalized_measure || measure_mode == unnormalized_measure || measure_mode == normalized_cutoff_measure || measure_mode == unnormalized_cutoff_measure)
|
---|
159 | //Defaults to 100 iteration to find minimum
|
---|
160 | _axesFinder = new AxesFinderFromKmeansMinimization(new AxesFinderFromKT(), beta, Rcutoff, 100);
|
---|
161 | else {
|
---|
162 | std::cerr << "Multi-pass minimization only set up for normalized_measure, unnormalized_measure, normalized_cutoff_measure, unnormalized_cutoff_measure." << std::endl;
|
---|
163 | exit(1);
|
---|
164 | }
|
---|
165 | break;
|
---|
166 | case manual_axes:
|
---|
167 | _axesFinder = new AxesFinderFromUserInput();
|
---|
168 | break;
|
---|
169 | // These options have been commented out because they have not been fully tested
|
---|
170 | // case wta2_kt_axes: // option for alpha = 2 added
|
---|
171 | // _axesFinder = new AxesFinderFromWTA2_KT();
|
---|
172 | // break;
|
---|
173 | // case wta2_ca_axes: // option for alpha = 2 added
|
---|
174 | // _axesFinder = new AxesFinderFromWTA2_CA();
|
---|
175 | // break;
|
---|
176 | // case onepass_wta2_kt_axes: // option for alpha = 2 added
|
---|
177 | // setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA2_KT(), beta, Rcutoff);
|
---|
178 | // break;
|
---|
179 | // case onepass_wta2_ca_axes: // option for alpha = 2 added
|
---|
180 | // setOnePassAxesFinder(measure_mode, new AxesFinderFromWTA2_CA(), beta, Rcutoff);
|
---|
181 | // break;
|
---|
182 | default:
|
---|
183 | assert(false);
|
---|
184 | break;
|
---|
185 | }
|
---|
186 |
|
---|
187 | }
|
---|
188 |
|
---|
189 | // setAxes for Manual mode
|
---|
190 | void Njettiness::setAxes(std::vector<fastjet::PseudoJet> myAxes) {
|
---|
191 | if (_current_axes_mode == manual_axes || _current_axes_mode == onepass_manual_axes) {
|
---|
192 | _currentAxes = myAxes;
|
---|
193 | }
|
---|
194 | else {
|
---|
195 | std::cerr << "You can only use setAxes if using manual_axes or onepass_manual_axes measure mode" << std::endl;
|
---|
196 | exit(1);
|
---|
197 | }
|
---|
198 | }
|
---|
199 |
|
---|
200 | // Calculates and returns all TauComponents that user would want.
|
---|
201 | // This information is stored in _current_tau_components for later access as well.
|
---|
202 | TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) {
|
---|
203 | if (inputJets.size() <= n_jets) { //if not enough particles, return zero
|
---|
204 | _currentAxes = inputJets;
|
---|
205 | _currentAxes.resize(n_jets,fastjet::PseudoJet(0.0,0.0,0.0,0.0));
|
---|
206 | _current_tau_components = TauComponents();
|
---|
207 | _seedAxes = _currentAxes;
|
---|
208 | } else {
|
---|
209 | _currentAxes = _axesFinder->getAxes(n_jets,inputJets,_currentAxes); // sets current Axes
|
---|
210 | _seedAxes = _axesFinder->seedAxes(); // sets seed Axes (if one pass minimization was used)
|
---|
211 | _current_tau_components = _measureFunction->result(inputJets, _currentAxes); // sets current Tau Values
|
---|
212 | }
|
---|
213 | return _current_tau_components;
|
---|
214 | }
|
---|
215 |
|
---|
216 |
|
---|
217 | // Partition a list of particles according to which N-jettiness axis they are closest to.
|
---|
218 | // Return a vector of length _currentAxes.size() (which should be N).
|
---|
219 | // Each vector element is a list of ints corresponding to the indices in
|
---|
220 | // particles of the particles belonging to that jet.
|
---|
221 | // TODO: Consider moving to MeasureFunction
|
---|
222 | std::vector<std::list<int> > Njettiness::getPartition(const std::vector<fastjet::PseudoJet> & particles) {
|
---|
223 | std::vector<std::list<int> > partitions(_currentAxes.size());
|
---|
224 |
|
---|
225 | for (unsigned i = 0; i < particles.size(); i++) {
|
---|
226 |
|
---|
227 | int j_min = -1;
|
---|
228 | // find minimum distance
|
---|
229 | double minR = std::numeric_limits<double>::max(); //large number
|
---|
230 | for (unsigned j = 0; j < _currentAxes.size(); j++) {
|
---|
231 | double tempR = _measureFunction->jet_distance_squared(particles[i],_currentAxes[j]); // delta R distance
|
---|
232 | if (tempR < minR) {
|
---|
233 | minR = tempR;
|
---|
234 | j_min = j;
|
---|
235 | }
|
---|
236 | }
|
---|
237 | if (_measureFunction->do_cluster(particles[i],_currentAxes[j_min])) partitions[j_min].push_back(i);
|
---|
238 | }
|
---|
239 | return partitions;
|
---|
240 | }
|
---|
241 |
|
---|
242 | // Having found axes, assign each particle in particles to an axis, and return a set of jets.
|
---|
243 | // Each jet is the sum of particles closest to an axis (Njet = Naxes).
|
---|
244 | // TODO: Consider moving to MeasureFunction
|
---|
245 | std::vector<fastjet::PseudoJet> Njettiness::getJets(const std::vector<fastjet::PseudoJet> & particles) {
|
---|
246 |
|
---|
247 | std::vector<fastjet::PseudoJet> jets(_currentAxes.size());
|
---|
248 |
|
---|
249 | std::vector<std::list<int> > partition = getPartition(particles);
|
---|
250 | for (unsigned j = 0; j < partition.size(); ++j) {
|
---|
251 | std::list<int>::const_iterator it, itE;
|
---|
252 | for (it = partition[j].begin(), itE = partition[j].end(); it != itE; ++it) {
|
---|
253 | jets[j] += particles[*it];
|
---|
254 | }
|
---|
255 | }
|
---|
256 | return jets;
|
---|
257 | }
|
---|
258 |
|
---|
259 | } // namespace contrib
|
---|
260 |
|
---|
261 | FASTJET_END_NAMESPACE
|
---|