Fork me on GitHub

source: git/external/fastjet/contribs/Nsubjettiness/Njettiness.cc@ 3b99510

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 3b99510 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: 7.9 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: Njettiness.cc 821 2015-06-15 18:50:53Z 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 "Njettiness.hh"
26
27FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
28
29namespace contrib {
30
31
32///////
33//
34// Main Njettiness Class
35//
36///////
37
38LimitedWarning Njettiness::_old_measure_warning;
39LimitedWarning Njettiness::_old_axes_warning;
40
41
42// Constructor
43Njettiness::Njettiness(const AxesDefinition & axes_def, const MeasureDefinition & measure_def)
44: _axes_def(axes_def.create()), _measure_def(measure_def.create()) {}
45
46// setAxes for Manual mode
47void Njettiness::setAxes(const std::vector<fastjet::PseudoJet> & myAxes) {
48 if (_axes_def()->needsManualAxes()) {
49 _currentAxes = myAxes;
50 } else {
51 throw Error("You can only use setAxes for manual AxesDefinitions");
52 }
53}
54
55// Calculates and returns all TauComponents that user would want.
56// This information is stored in _current_tau_components for later access as well.
57TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) const {
58 if (inputJets.size() <= n_jets) { //if not enough particles, return zero
59 _currentAxes = inputJets;
60 _currentAxes.resize(n_jets,fastjet::PseudoJet(0.0,0.0,0.0,0.0));
61 _current_tau_components = TauComponents();
62 _seedAxes = _currentAxes;
63 _currentPartition = TauPartition(n_jets); // empty partition
64 } else {
65 assert(_axes_def()); // this should never fail.
66
67 if (_axes_def()->needsManualAxes()) { // if manual mode
68 // take current axes as seeds
69 _seedAxes = _currentAxes;
70
71 // refine axes if requested
72 _currentAxes = _axes_def->get_refined_axes(n_jets,inputJets,_seedAxes, _measure_def());
73 } else { // non-manual axes
74
75 //set starting point for minimization
76 _seedAxes = _axes_def->get_starting_axes(n_jets,inputJets,_measure_def());
77
78 // refine axes as needed
79 _currentAxes = _axes_def->get_refined_axes(n_jets,inputJets,_seedAxes, _measure_def());
80
81 // NOTE: The above two function calls are combined in "AxesDefinition::get_axes"
82 // but are separated here to allow seed axes to be stored.
83 }
84
85 // Find and store partition
86 _currentPartition = _measure_def->get_partition(inputJets,_currentAxes);
87
88 // Find and store tau value
89 _current_tau_components = _measure_def->component_result_from_partition(_currentPartition, _currentAxes); // sets current Tau Values
90 }
91 return _current_tau_components;
92}
93
94
95///////
96//
97// Below is code for backward compatibility to use the old interface.
98// May be deleted in a future version
99//
100///////
101
102Njettiness::Njettiness(AxesMode axes_mode, const MeasureDefinition & measure_def)
103: _axes_def(createAxesDef(axes_mode)), _measure_def(measure_def.create()) {}
104
105// Convert from MeasureMode enum to MeasureDefinition
106// This returns a pointer that will be claimed by a SharedPtr
107MeasureDefinition* Njettiness::createMeasureDef(MeasureMode measure_mode, int num_para, double para1, double para2, double para3) const {
108
109 _old_measure_warning.warn("Njettiness::createMeasureDef: You are using the old MeasureMode way of specifying N-subjettiness measures. This is deprecated as of v2.1 and will be removed in v3.0. Please use MeasureDefinition instead.");
110
111 // definition of maximum Rcutoff for non-cutoff measures, changed later by other measures
112 double Rcutoff = std::numeric_limits<double>::max(); //large number
113 // Most (but not all) measures have some kind of beta value
114 double beta = std::numeric_limits<double>::quiet_NaN();
115 // The normalized measures have an R0 value.
116 double R0 = std::numeric_limits<double>::quiet_NaN();
117
118 // Find the MeasureFunction and set the parameters.
119 switch (measure_mode) {
120 case normalized_measure:
121 beta = para1;
122 R0 = para2;
123 if(num_para == 2) {
124 return new NormalizedMeasure(beta,R0);
125 } else {
126 throw Error("normalized_measure needs 2 parameters (beta and R0)");
127 }
128 break;
129 case unnormalized_measure:
130 beta = para1;
131 if(num_para == 1) {
132 return new UnnormalizedMeasure(beta);
133 } else {
134 throw Error("unnormalized_measure needs 1 parameter (beta)");
135 }
136 break;
137 case geometric_measure:
138 throw Error("This class has been removed. Please use OriginalGeometricMeasure, ModifiedGeometricMeasure, or ConicalGeometricMeasure with the new Njettiness constructor.");
139 break;
140 case normalized_cutoff_measure:
141 beta = para1;
142 R0 = para2;
143 Rcutoff = para3; //Rcutoff parameter is 3rd parameter in normalized_cutoff_measure
144 if (num_para == 3) {
145 return new NormalizedCutoffMeasure(beta,R0,Rcutoff);
146 } else {
147 throw Error("normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)");
148 }
149 break;
150 case unnormalized_cutoff_measure:
151 beta = para1;
152 Rcutoff = para2; //Rcutoff parameter is 2nd parameter in normalized_cutoff_measure
153 if (num_para == 2) {
154 return new UnnormalizedCutoffMeasure(beta,Rcutoff);
155 } else {
156 throw Error("unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)");
157 }
158 break;
159 case geometric_cutoff_measure:
160 throw Error("This class has been removed. Please use OriginalGeometricMeasure, ModifiedGeometricMeasure, or ConicalGeometricMeasure with the new Njettiness constructor.");
161 default:
162 assert(false);
163 break;
164 }
165 return NULL;
166}
167
168// Convert from AxesMode enum to AxesDefinition
169// This returns a pointer that will be claimed by a SharedPtr
170AxesDefinition* Njettiness::createAxesDef(Njettiness::AxesMode axes_mode) const {
171
172 _old_axes_warning.warn("Njettiness::createAxesDef: You are using the old AxesMode way of specifying N-subjettiness axes. This is deprecated as of v2.1 and will be removed in v3.0. Please use AxesDefinition instead.");
173
174
175 switch (axes_mode) {
176 case wta_kt_axes:
177 return new WTA_KT_Axes();
178 case wta_ca_axes:
179 return new WTA_CA_Axes();
180 case kt_axes:
181 return new KT_Axes();
182 case ca_axes:
183 return new CA_Axes();
184 case antikt_0p2_axes:
185 return new AntiKT_Axes(0.2);
186 case onepass_wta_kt_axes:
187 return new OnePass_WTA_KT_Axes();
188 case onepass_wta_ca_axes:
189 return new OnePass_WTA_CA_Axes();
190 case onepass_kt_axes:
191 return new OnePass_KT_Axes();
192 case onepass_ca_axes:
193 return new OnePass_CA_Axes();
194 case onepass_antikt_0p2_axes:
195 return new OnePass_AntiKT_Axes(0.2);
196 case onepass_manual_axes:
197 return new OnePass_Manual_Axes();
198 case min_axes:
199 return new MultiPass_Axes(100);
200 case manual_axes:
201 return new Manual_Axes();
202 default:
203 assert(false);
204 return NULL;
205 }
206}
207
208
209
210
211
212} // namespace contrib
213
214FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.