| 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: AxesFinder.cc 670 2014-06-06 01:24:42Z 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 "AxesFinder.hh"
|
|---|
| 26 |
|
|---|
| 27 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 28 |
|
|---|
| 29 | namespace contrib{
|
|---|
| 30 |
|
|---|
| 31 | ///////
|
|---|
| 32 | //
|
|---|
| 33 | // Functions for minimization.
|
|---|
| 34 | //
|
|---|
| 35 | ///////
|
|---|
| 36 |
|
|---|
| 37 | // Given starting axes, update to find better axes by using Kmeans clustering around the old axes
|
|---|
| 38 | template <int N>
|
|---|
| 39 | std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
|
|---|
| 40 | const std::vector <fastjet::PseudoJet> & inputJets) const {
|
|---|
| 41 | assert(old_axes.size() == N);
|
|---|
| 42 |
|
|---|
| 43 | // some storage, declared static to save allocation/re-allocation costs
|
|---|
| 44 | static LightLikeAxis new_axes[N];
|
|---|
| 45 | static fastjet::PseudoJet new_jets[N];
|
|---|
| 46 | for (int n = 0; n < N; ++n) {
|
|---|
| 47 | new_axes[n].reset(0.0,0.0,0.0,0.0);
|
|---|
| 48 | new_jets[n].reset_momentum(0.0,0.0,0.0,0.0);
|
|---|
| 49 | }
|
|---|
| 50 |
|
|---|
| 51 | double precision = _precision;
|
|---|
| 52 |
|
|---|
| 53 | /////////////// Assignment Step //////////////////////////////////////////////////////////
|
|---|
| 54 | std::vector<int> assignment_index(inputJets.size());
|
|---|
| 55 | int k_assign = -1;
|
|---|
| 56 |
|
|---|
| 57 | for (unsigned i = 0; i < inputJets.size(); i++){
|
|---|
| 58 | double smallestDist = std::numeric_limits<double>::max(); //large number
|
|---|
| 59 | for (int k = 0; k < N; k++) {
|
|---|
| 60 | double thisDist = old_axes[k].DistanceSq(inputJets[i]);
|
|---|
| 61 | if (thisDist < smallestDist) {
|
|---|
| 62 | smallestDist = thisDist;
|
|---|
| 63 | k_assign = k;
|
|---|
| 64 | }
|
|---|
| 65 | }
|
|---|
| 66 | if (smallestDist > sq(_Rcutoff)) {k_assign = -1;}
|
|---|
| 67 | assignment_index[i] = k_assign;
|
|---|
| 68 | }
|
|---|
| 69 |
|
|---|
| 70 | //////////////// Update Step /////////////////////////////////////////////////////////////
|
|---|
| 71 | double distPhi, old_dist;
|
|---|
| 72 | for (unsigned i = 0; i < inputJets.size(); i++) {
|
|---|
| 73 | int old_jet_i = assignment_index[i];
|
|---|
| 74 | if (old_jet_i == -1) {continue;}
|
|---|
| 75 |
|
|---|
| 76 | const fastjet::PseudoJet& inputJet_i = inputJets[i];
|
|---|
| 77 | LightLikeAxis& new_axis_i = new_axes[old_jet_i];
|
|---|
| 78 | double inputPhi_i = inputJet_i.phi();
|
|---|
| 79 | double inputRap_i = inputJet_i.rap();
|
|---|
| 80 |
|
|---|
| 81 | // optimize pow() call
|
|---|
| 82 | // add noise (the precision term) to make sure we don't divide by zero
|
|---|
| 83 | if (_beta == 1.0) {
|
|---|
| 84 | double DR = std::sqrt(sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i));
|
|---|
| 85 | old_dist = 1.0/DR;
|
|---|
| 86 | } else if (_beta == 2.0) {
|
|---|
| 87 | old_dist = 1.0;
|
|---|
| 88 | } else if (_beta == 0.0) {
|
|---|
| 89 | double DRSq = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
|
|---|
| 90 | old_dist = 1.0/DRSq;
|
|---|
| 91 | } else {
|
|---|
| 92 | old_dist = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i);
|
|---|
| 93 | old_dist = std::pow(old_dist, (0.5*_beta-1.0));
|
|---|
| 94 | }
|
|---|
| 95 |
|
|---|
| 96 | // TODO: Put some of these addition functions into light-like axes
|
|---|
| 97 | // rapidity sum
|
|---|
| 98 | new_axis_i.set_rap(new_axis_i.rap() + inputJet_i.perp() * inputRap_i * old_dist);
|
|---|
| 99 | // phi sum
|
|---|
| 100 | distPhi = inputPhi_i - old_axes[old_jet_i].phi();
|
|---|
| 101 | if (fabs(distPhi) <= M_PI){
|
|---|
| 102 | new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * inputPhi_i * old_dist );
|
|---|
| 103 | } else if (distPhi > M_PI) {
|
|---|
| 104 | new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (-2*M_PI + inputPhi_i) * old_dist );
|
|---|
| 105 | } else if (distPhi < -M_PI) {
|
|---|
| 106 | new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (+2*M_PI + inputPhi_i) * old_dist );
|
|---|
| 107 | }
|
|---|
| 108 | // weights sum
|
|---|
| 109 | new_axis_i.set_weight( new_axis_i.weight() + inputJet_i.perp() * old_dist );
|
|---|
| 110 | // momentum magnitude sum
|
|---|
| 111 | new_jets[old_jet_i] += inputJet_i;
|
|---|
| 112 | }
|
|---|
| 113 | // normalize sums
|
|---|
| 114 | for (int k = 0; k < N; k++) {
|
|---|
| 115 | if (new_axes[k].weight() == 0) {
|
|---|
| 116 | // no particles were closest to this axis! Return to old axis instead of (0,0,0,0)
|
|---|
| 117 | new_axes[k] = old_axes[k];
|
|---|
| 118 | } else {
|
|---|
| 119 | new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() );
|
|---|
| 120 | new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() );
|
|---|
| 121 | new_axes[k].set_phi( std::fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) );
|
|---|
| 122 | new_axes[k].set_mom( std::sqrt(new_jets[k].modp2()) );
|
|---|
| 123 | }
|
|---|
| 124 | }
|
|---|
| 125 | std::vector<LightLikeAxis> new_axes_vec(N);
|
|---|
| 126 | for (unsigned k = 0; k < N; ++k) new_axes_vec[k] = new_axes[k];
|
|---|
| 127 | return new_axes_vec;
|
|---|
| 128 | }
|
|---|
| 129 |
|
|---|
| 130 | // Given N starting axes, this function updates all axes to find N better axes.
|
|---|
| 131 | // (This is just a wrapper for the templated version above.)
|
|---|
| 132 | std::vector<LightLikeAxis> AxesFinderFromOnePassMinimization::UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
|
|---|
| 133 | const std::vector <fastjet::PseudoJet> & inputJets) const {
|
|---|
| 134 | int N = old_axes.size();
|
|---|
| 135 | switch (N) {
|
|---|
| 136 | case 1: return UpdateAxesFast<1>(old_axes, inputJets);
|
|---|
| 137 | case 2: return UpdateAxesFast<2>(old_axes, inputJets);
|
|---|
| 138 | case 3: return UpdateAxesFast<3>(old_axes, inputJets);
|
|---|
| 139 | case 4: return UpdateAxesFast<4>(old_axes, inputJets);
|
|---|
| 140 | case 5: return UpdateAxesFast<5>(old_axes, inputJets);
|
|---|
| 141 | case 6: return UpdateAxesFast<6>(old_axes, inputJets);
|
|---|
| 142 | case 7: return UpdateAxesFast<7>(old_axes, inputJets);
|
|---|
| 143 | case 8: return UpdateAxesFast<8>(old_axes, inputJets);
|
|---|
| 144 | case 9: return UpdateAxesFast<9>(old_axes, inputJets);
|
|---|
| 145 | case 10: return UpdateAxesFast<10>(old_axes, inputJets);
|
|---|
| 146 | case 11: return UpdateAxesFast<11>(old_axes, inputJets);
|
|---|
| 147 | case 12: return UpdateAxesFast<12>(old_axes, inputJets);
|
|---|
| 148 | case 13: return UpdateAxesFast<13>(old_axes, inputJets);
|
|---|
| 149 | case 14: return UpdateAxesFast<14>(old_axes, inputJets);
|
|---|
| 150 | case 15: return UpdateAxesFast<15>(old_axes, inputJets);
|
|---|
| 151 | case 16: return UpdateAxesFast<16>(old_axes, inputJets);
|
|---|
| 152 | case 17: return UpdateAxesFast<17>(old_axes, inputJets);
|
|---|
| 153 | case 18: return UpdateAxesFast<18>(old_axes, inputJets);
|
|---|
| 154 | case 19: return UpdateAxesFast<19>(old_axes, inputJets);
|
|---|
| 155 | case 20: return UpdateAxesFast<20>(old_axes, inputJets);
|
|---|
| 156 | default: std::cout << "N-jettiness is hard-coded to only allow up to 20 jets!" << std::endl;
|
|---|
| 157 | return std::vector<LightLikeAxis>();
|
|---|
| 158 | }
|
|---|
| 159 |
|
|---|
| 160 | }
|
|---|
| 161 |
|
|---|
| 162 | // uses minimization of N-jettiness to continually update axes until convergence.
|
|---|
| 163 | // The function returns the axes found at the (local) minimum
|
|---|
| 164 | std::vector<fastjet::PseudoJet> AxesFinderFromOnePassMinimization::getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) const {
|
|---|
| 165 |
|
|---|
| 166 | // convert from PseudoJets to LightLikeAxes
|
|---|
| 167 | std::vector< LightLikeAxis > old_axes(n_jets, LightLikeAxis(0,0,0,0));
|
|---|
| 168 | for (int k = 0; k < n_jets; k++) {
|
|---|
| 169 | old_axes[k].set_rap( seedAxes[k].rap() );
|
|---|
| 170 | old_axes[k].set_phi( seedAxes[k].phi() );
|
|---|
| 171 | }
|
|---|
| 172 |
|
|---|
| 173 | // Find new axes by iterating (only one pass here)
|
|---|
| 174 | std::vector< LightLikeAxis > new_axes(n_jets, LightLikeAxis(0,0,0,0));
|
|---|
| 175 | double cmp = std::numeric_limits<double>::max(); //large number
|
|---|
| 176 | int h = 0;
|
|---|
| 177 | while (cmp > _precision && h < _halt) { // Keep updating axes until near-convergence or too many update steps
|
|---|
| 178 | cmp = 0.0;
|
|---|
| 179 | h++;
|
|---|
| 180 | new_axes = UpdateAxes(old_axes, inputJets); // Update axes
|
|---|
| 181 | for (int k = 0; k < n_jets; k++) {
|
|---|
| 182 | cmp += old_axes[k].Distance(new_axes[k]);
|
|---|
| 183 | }
|
|---|
| 184 | cmp = cmp / ((double) n_jets);
|
|---|
| 185 | old_axes = new_axes;
|
|---|
| 186 | }
|
|---|
| 187 |
|
|---|
| 188 | // Convert from internal LightLikeAxes to PseudoJet
|
|---|
| 189 | std::vector<fastjet::PseudoJet> outputAxes;
|
|---|
| 190 | for (int k = 0; k < n_jets; k++) {
|
|---|
| 191 | fastjet::PseudoJet temp = old_axes[k].ConvertToPseudoJet();
|
|---|
| 192 | outputAxes.push_back(temp);
|
|---|
| 193 | }
|
|---|
| 194 |
|
|---|
| 195 | // this is used to debug the minimization routine to make sure that it works.
|
|---|
| 196 | bool do_debug = false;
|
|---|
| 197 | if (do_debug) {
|
|---|
| 198 | // get this information to make sure that minimization is working properly
|
|---|
| 199 | TauComponents seed_tau_components = _measureFunction.result(inputJets, seedAxes);
|
|---|
| 200 | double seed_tau = seed_tau_components.tau();
|
|---|
| 201 | TauComponents tau_components = _measureFunction.result(inputJets, outputAxes);
|
|---|
| 202 | double outputTau = tau_components.tau();
|
|---|
| 203 | assert(outputTau <= seed_tau);
|
|---|
| 204 | }
|
|---|
| 205 |
|
|---|
| 206 | return outputAxes;
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | PseudoJet AxesFinderFromKmeansMinimization::jiggle(const PseudoJet& axis) const {
|
|---|
| 210 | double phi_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
|
|---|
| 211 | double rap_noise = ((double)rand()/(double)RAND_MAX) * _noise_range * 2.0 - _noise_range;
|
|---|
| 212 |
|
|---|
| 213 | double new_phi = axis.phi() + phi_noise;
|
|---|
| 214 | if (new_phi >= 2.0*M_PI) new_phi -= 2.0*M_PI;
|
|---|
| 215 | if (new_phi <= -2.0*M_PI) new_phi += 2.0*M_PI;
|
|---|
| 216 |
|
|---|
| 217 | PseudoJet newAxis(0,0,0,0);
|
|---|
| 218 | newAxis.reset_PtYPhiM(axis.perp(),axis.rap() + rap_noise,new_phi);
|
|---|
| 219 | return newAxis;
|
|---|
| 220 | }
|
|---|
| 221 |
|
|---|
| 222 |
|
|---|
| 223 | // Repeatedly calls the one pass finder to try to find global minimum
|
|---|
| 224 | std::vector<fastjet::PseudoJet> AxesFinderFromKmeansMinimization::getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& seedAxes) const {
|
|---|
| 225 |
|
|---|
| 226 | // first iteration
|
|---|
| 227 | std::vector<fastjet::PseudoJet> bestAxes = _onePassFinder.getAxes(n_jets, inputJets, seedAxes);
|
|---|
| 228 |
|
|---|
| 229 | double bestTau = (_measureFunction.result(inputJets,bestAxes)).tau();
|
|---|
| 230 |
|
|---|
| 231 | for (int l = 1; l < _n_iterations; l++) { // Do minimization procedure multiple times (l = 1 to start since first iteration is done already)
|
|---|
| 232 |
|
|---|
| 233 | // Add noise to current best axes
|
|---|
| 234 | std::vector< PseudoJet > noiseAxes(n_jets, PseudoJet(0,0,0,0));
|
|---|
| 235 | for (int k = 0; k < n_jets; k++) {
|
|---|
| 236 | noiseAxes[k] = jiggle(bestAxes[k]);
|
|---|
| 237 | }
|
|---|
| 238 |
|
|---|
| 239 | std::vector<fastjet::PseudoJet> testAxes = _onePassFinder.getAxes(n_jets, inputJets, noiseAxes);
|
|---|
| 240 | double testTau = (_measureFunction.result(inputJets,testAxes)).tau();
|
|---|
| 241 |
|
|---|
| 242 | if (testTau < bestTau) {
|
|---|
| 243 | bestTau = testTau;
|
|---|
| 244 | bestAxes = testAxes;
|
|---|
| 245 | }
|
|---|
| 246 | }
|
|---|
| 247 |
|
|---|
| 248 | return bestAxes;
|
|---|
| 249 | }
|
|---|
| 250 |
|
|---|
| 251 | // Uses minimization of the geometric distance in order to find the minimum axes.
|
|---|
| 252 | // It continually updates until it reaches convergence or it reaches the maximum number of attempts.
|
|---|
| 253 | // This is essentially the same as a stable cone finder.
|
|---|
| 254 | std::vector<fastjet::PseudoJet> AxesFinderFromGeometricMinimization::getAxes(int /*n_jets*/, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes) const {
|
|---|
| 255 |
|
|---|
| 256 | std::vector<fastjet::PseudoJet> seedAxes = currentAxes;
|
|---|
| 257 | double seedTau = _function.tau(particles, seedAxes);
|
|---|
| 258 |
|
|---|
| 259 | for (int i = 0; i < _nAttempts; i++) {
|
|---|
| 260 |
|
|---|
| 261 | std::vector<fastjet::PseudoJet> newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0));
|
|---|
| 262 |
|
|---|
| 263 | // find closest axis and assign to that
|
|---|
| 264 | for (unsigned int i = 0; i < particles.size(); i++) {
|
|---|
| 265 |
|
|---|
| 266 | // start from unclustered beam measure
|
|---|
| 267 | int minJ = -1;
|
|---|
| 268 | double minDist = _function.beam_distance_squared(particles[i]);
|
|---|
| 269 |
|
|---|
| 270 | // which axis am I closest to?
|
|---|
| 271 | for (unsigned int j = 0; j < seedAxes.size(); j++) {
|
|---|
| 272 | double tempDist = _function.jet_distance_squared(particles[i],seedAxes[j]);
|
|---|
| 273 | if (tempDist < minDist) {
|
|---|
| 274 | minDist = tempDist;
|
|---|
| 275 | minJ = j;
|
|---|
| 276 | }
|
|---|
| 277 | }
|
|---|
| 278 |
|
|---|
| 279 | // if not unclustered, then cluster
|
|---|
| 280 | if (minJ != -1) newAxes[minJ] += particles[i];
|
|---|
| 281 | }
|
|---|
| 282 |
|
|---|
| 283 | // calculate tau on new axes
|
|---|
| 284 | seedAxes = newAxes;
|
|---|
| 285 | double tempTau = _function.tau(particles, newAxes);
|
|---|
| 286 |
|
|---|
| 287 | // close enough to stop?
|
|---|
| 288 | if (fabs(tempTau - seedTau) < _accuracy) break;
|
|---|
| 289 | seedTau = tempTau;
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 | return seedAxes;
|
|---|
| 293 | }
|
|---|
| 294 |
|
|---|
| 295 | // Go from internal LightLikeAxis to PseudoJet
|
|---|
| 296 | fastjet::PseudoJet LightLikeAxis::ConvertToPseudoJet() {
|
|---|
| 297 | double px, py, pz, E;
|
|---|
| 298 | E = _mom;
|
|---|
| 299 | pz = (std::exp(2.0*_rap) - 1.0) / (std::exp(2.0*_rap) + 1.0) * E;
|
|---|
| 300 | px = std::cos(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
|
|---|
| 301 | py = std::sin(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) );
|
|---|
| 302 | return fastjet::PseudoJet(px,py,pz,E);
|
|---|
| 303 | }
|
|---|
| 304 |
|
|---|
| 305 | } //namespace contrib
|
|---|
| 306 |
|
|---|
| 307 | FASTJET_END_NAMESPACE
|
|---|
