Fork me on GitHub

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

ImprovedOutputFile Timing dual_readout llp
Last change on this file since fc4e460 was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 8.1 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 677 2014-06-12 18:56:46Z 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
38Njettiness::Njettiness(const AxesDefinition & axes_def, const MeasureDefinition & measure_def)
39: _axes_def(axes_def.create()), _measure_def(measure_def.create()) {
40 setMeasureFunctionAndAxesFinder(); // call helper function to do the hard work
41}
42
43Njettiness::Njettiness(AxesMode axes_mode, const MeasureDefinition & measure_def)
44: _axes_def(createAxesDef(axes_mode)), _measure_def(measure_def.create()) {
45 setMeasureFunctionAndAxesFinder(); // call helper function to do the hard work
46}
47
48// Convert from MeasureMode enum to MeasureDefinition
49// This returns a pointer that will be claimed by a SharedPtr
50MeasureDefinition* Njettiness::createMeasureDef(MeasureMode measure_mode, int num_para, double para1, double para2, double para3) const {
51
52 // definition of maximum Rcutoff for non-cutoff measures, changed later by other measures
53 double Rcutoff = std::numeric_limits<double>::max(); //large number
54 // Most (but not all) measures have some kind of beta value
55 double beta = std::numeric_limits<double>::quiet_NaN();
56 // The normalized measures have an R0 value.
57 double R0 = std::numeric_limits<double>::quiet_NaN();
58
59 // Find the MeasureFunction and set the parameters.
60 switch (measure_mode) {
61 case normalized_measure:
62 beta = para1;
63 R0 = para2;
64 if(num_para == 2) {
65 return new NormalizedMeasure(beta,R0);
66 } else {
67 throw Error("normalized_measure needs 2 parameters (beta and R0)");
68 }
69 break;
70 case unnormalized_measure:
71 beta = para1;
72 if(num_para == 1) {
73 return new UnnormalizedMeasure(beta);
74 } else {
75 throw Error("unnormalized_measure needs 1 parameter (beta)");
76 }
77 break;
78 case geometric_measure:
79 beta = para1;
80 if (num_para == 1) {
81 return new GeometricMeasure(beta);
82 } else {
83 throw Error("geometric_measure needs 1 parameter (beta)");
84 }
85 break;
86 case normalized_cutoff_measure:
87 beta = para1;
88 R0 = para2;
89 Rcutoff = para3; //Rcutoff parameter is 3rd parameter in normalized_cutoff_measure
90 if (num_para == 3) {
91 return new NormalizedCutoffMeasure(beta,R0,Rcutoff);
92 } else {
93 throw Error("normalized_cutoff_measure has 3 parameters (beta, R0, Rcutoff)");
94 }
95 break;
96 case unnormalized_cutoff_measure:
97 beta = para1;
98 Rcutoff = para2; //Rcutoff parameter is 2nd parameter in normalized_cutoff_measure
99 if (num_para == 2) {
100 return new UnnormalizedCutoffMeasure(beta,Rcutoff);
101 } else {
102 throw Error("unnormalized_cutoff_measure has 2 parameters (beta, Rcutoff)");
103 }
104 break;
105 case geometric_cutoff_measure:
106 beta = para1;
107 Rcutoff = para2; //Rcutoff parameter is 2nd parameter in geometric_cutoff_measure
108 if(num_para == 2) {
109 return new GeometricCutoffMeasure(beta,Rcutoff);
110 } else {
111 throw Error("geometric_cutoff_measure has 2 parameters (beta, Rcutoff)");
112 }
113 break;
114 default:
115 assert(false);
116 break;
117 }
118 return NULL;
119}
120
121// Convert from AxesMode enum to AxesDefinition
122// This returns a pointer that will be claimed by a SharedPtr
123AxesDefinition* Njettiness::createAxesDef(Njettiness::AxesMode axes_mode) const {
124
125 switch (axes_mode) {
126 case wta_kt_axes:
127 return new WTA_KT_Axes();
128 case wta_ca_axes:
129 return new WTA_CA_Axes();
130 case kt_axes:
131 return new KT_Axes();
132 case ca_axes:
133 return new CA_Axes();
134 case antikt_0p2_axes:
135 return new AntiKT_Axes(0.2);
136 case onepass_wta_kt_axes:
137 return new OnePass_WTA_KT_Axes();
138 case onepass_wta_ca_axes:
139 return new OnePass_WTA_CA_Axes();
140 case onepass_kt_axes:
141 return new OnePass_KT_Axes();
142 case onepass_ca_axes:
143 return new OnePass_CA_Axes();
144 case onepass_antikt_0p2_axes:
145 return new OnePass_AntiKT_Axes(0.2);
146 case onepass_manual_axes:
147 return new OnePass_Manual_Axes();
148 case min_axes:
149 return new MultiPass_Axes(100);
150 case manual_axes:
151 return new Manual_Axes();
152 default:
153 assert(false);
154 return NULL;
155 }
156}
157
158
159// Parsing needed for constructor to set AxesFinder and MeasureFunction
160// All of the parameter handling is here, and checking that number of parameters is correct.
161void Njettiness::setMeasureFunctionAndAxesFinder() {
162 // Get the correct MeasureFunction and AxesFinders
163 _measureFunction.reset(_measure_def->createMeasureFunction());
164 _startingAxesFinder.reset(_axes_def->createStartingAxesFinder(*_measure_def));
165 _finishingAxesFinder.reset(_axes_def->createFinishingAxesFinder(*_measure_def));
166}
167
168// setAxes for Manual mode
169void Njettiness::setAxes(const std::vector<fastjet::PseudoJet> & myAxes) {
170 if (_axes_def->supportsManualAxes()) {
171 _currentAxes = myAxes;
172 } else {
173 throw Error("You can only use setAxes for manual AxesDefinitions");
174 }
175}
176
177// Calculates and returns all TauComponents that user would want.
178// This information is stored in _current_tau_components for later access as well.
179TauComponents Njettiness::getTauComponents(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets) const {
180 if (inputJets.size() <= n_jets) { //if not enough particles, return zero
181 _currentAxes = inputJets;
182 _currentAxes.resize(n_jets,fastjet::PseudoJet(0.0,0.0,0.0,0.0));
183 _current_tau_components = TauComponents();
184 _seedAxes = _currentAxes;
185 _currentJets = _currentAxes;
186 _currentBeam = PseudoJet(0.0,0.0,0.0,0.0);
187 } else {
188
189 _seedAxes = _startingAxesFinder->getAxes(n_jets,inputJets,_currentAxes); //sets starting point for minimization
190 if (_finishingAxesFinder) {
191 _currentAxes = _finishingAxesFinder->getAxes(n_jets,inputJets,_seedAxes);
192 } else {
193 _currentAxes = _seedAxes;
194 }
195
196 // Find partition and store information
197 // (jet information in _currentJets, beam in _currentBeam)
198 _currentJets = _measureFunction->get_partition(inputJets,_currentAxes,&_currentBeam);
199
200 // Find tau value and store information
201 _current_tau_components = _measureFunction->result_from_partition(_currentJets, _currentAxes,&_currentBeam); // sets current Tau Values
202 }
203 return _current_tau_components;
204}
205
206
207// Partition a list of particles according to which N-jettiness axis they are closest to.
208// Return a vector of length _currentAxes.size() (which should be N).
209// Each vector element is a list of ints corresponding to the indices in
210// particles of the particles belonging to that jet.
211std::vector<std::list<int> > Njettiness::getPartitionList(const std::vector<fastjet::PseudoJet> & particles) const {
212 // core code is in MeasureFunction
213 return _measureFunction->get_partition_list(particles,_currentAxes);
214}
215
216
217} // namespace contrib
218
219FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.