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