// Nsubjettiness Package // Questions/Comments? jthaler@jthaler.net // // Copyright (c) 2011-14 // Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason // // $Id: MeasureDefinition.cc 946 2016-06-14 19:11:27Z jthaler $ //---------------------------------------------------------------------- // This file is part of FastJet contrib. // // It is free software; you can redistribute it and/or modify it under // the terms of the GNU General Public License as published by the // Free Software Foundation; either version 2 of the License, or (at // your option) any later version. // // It is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public // License for more details. // // You should have received a copy of the GNU General Public License // along with this code. If not, see . //---------------------------------------------------------------------- // #include "AxesRefiner.hh" #include "MeasureDefinition.hh" #include FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh namespace contrib { /////// // // Measure Function // /////// //descriptions updated to include measure type std::string DefaultMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Default Measure (should not be used directly)"; return stream.str(); }; std::string NormalizedMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Normalized Measure (beta = " << _beta << ", R0 = " << _R0 << ")"; return stream.str(); }; std::string UnnormalizedMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Unnormalized Measure (beta = " << _beta << ", in GeV)"; return stream.str(); }; std::string NormalizedCutoffMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Normalized Cutoff Measure (beta = " << _beta << ", R0 = " << _R0 << ", Rcut = " << _Rcutoff << ")"; return stream.str(); }; std::string UnnormalizedCutoffMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Unnormalized Cutoff Measure (beta = " << _beta << ", Rcut = " << _Rcutoff << ", in GeV)"; return stream.str(); }; //std::string DeprecatedGeometricMeasure::description() const { // std::stringstream stream; // stream << std::fixed << std::setprecision(2) // << "Deprecated Geometric Measure (beta = " << _jet_beta << ", in GeV)"; // return stream.str(); //}; //std::string DeprecatedGeometricCutoffMeasure::description() const { // std::stringstream stream; // stream << std::fixed << std::setprecision(2) // << "Deprecated Geometric Cutoff Measure (beta = " << _jet_beta << ", Rcut = " << _Rcutoff << ", in GeV)"; // return stream.str(); //}; std::string ConicalMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Conical Measure (beta = " << _beta << ", Rcut = " << _Rcutoff << ", in GeV)"; return stream.str(); }; std::string OriginalGeometricMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Original Geometric Measure (Rcut = " << _Rcutoff << ", in GeV)"; return stream.str(); }; std::string ModifiedGeometricMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Modified Geometric Measure (Rcut = " << _Rcutoff << ", in GeV)"; return stream.str(); }; std::string ConicalGeometricMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "Conical Geometric Measure (beta = " << _jet_beta << ", gamma = " << _beam_gamma << ", Rcut = " << _Rcutoff << ", in GeV)"; return stream.str(); }; std::string XConeMeasure::description() const { std::stringstream stream; stream << std::fixed << std::setprecision(2) << "XCone Measure (beta = " << _jet_beta << ", Rcut = " << _Rcutoff << ", in GeV)"; return stream.str(); }; // Return all of the necessary TauComponents for specific input particles and axes TauComponents MeasureDefinition::component_result(const std::vector& particles, const std::vector& axes) const { // first find partition TauPartition partition = get_partition(particles,axes); // then return result calculated from partition return component_result_from_partition(partition,axes); } TauPartition MeasureDefinition::get_partition(const std::vector& particles, const std::vector& axes) const { TauPartition myPartition(axes.size()); // Figures out the partiting of the input particles into the various jet pieces // Based on which axis the parition is closest to for (unsigned i = 0; i < particles.size(); i++) { // find minimum distance; start with beam (-1) for reference int j_min = -1; double minRsq; if (has_beam()) minRsq = beam_distance_squared(particles[i]); else minRsq = std::numeric_limits::max(); // make it large value // check to see which axis the particle is closest to for (unsigned j = 0; j < axes.size(); j++) { double tempRsq = jet_distance_squared(particles[i],axes[j]); // delta R distance if (tempRsq < minRsq) { minRsq = tempRsq; j_min = j; } } if (j_min == -1) { assert(has_beam()); // should have beam for this to make sense. myPartition.push_back_beam(particles[i],i); } else { myPartition.push_back_jet(j_min,particles[i],i); } } return myPartition; } // Uses existing partition and calculates result // TODO: Can we cache this for speed up when doing area subtraction? TauComponents MeasureDefinition::component_result_from_partition(const TauPartition& partition, const std::vector& axes) const { std::vector jetPieces(axes.size(), 0.0); double beamPiece = 0.0; double tauDen = 0.0; if (!has_denominator()) tauDen = 1.0; // if no denominator, then 1.0 for no normalization factor // first find jet pieces for (unsigned j = 0; j < axes.size(); j++) { std::vector thisPartition = partition.jet(j).constituents(); for (unsigned i = 0; i < thisPartition.size(); i++) { jetPieces[j] += jet_numerator(thisPartition[i],axes[j]); //numerator jet piece if (has_denominator()) tauDen += denominator(thisPartition[i]); // denominator } } // then find beam piece if (has_beam()) { std::vector beamPartition = partition.beam().constituents(); for (unsigned i = 0; i < beamPartition.size(); i++) { beamPiece += beam_numerator(beamPartition[i]); //numerator beam piece if (has_denominator()) tauDen += denominator(beamPartition[i]); // denominator } } // create jets for storage in TauComponents std::vector jets = partition.jets(); return TauComponents(_tau_mode, jetPieces, beamPiece, tauDen, jets, axes); } // new methods added to generalize energy and angle squared for different measure types double DefaultMeasure::energy(const PseudoJet& jet) const { double energy; switch (_measure_type) { case pt_R : case perp_lorentz_dot : energy = jet.perp(); break; case E_theta : case lorentz_dot : energy = jet.e(); break; default : { assert(_measure_type == pt_R || _measure_type == E_theta || _measure_type == lorentz_dot || _measure_type == perp_lorentz_dot); energy = std::numeric_limits::quiet_NaN(); break; } } return energy; } double DefaultMeasure::angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const { double pseudoRsquared; switch(_measure_type) { case pt_R : { pseudoRsquared = jet1.squared_distance(jet2); break; } case E_theta : { // doesn't seem to be a fastjet built in for this double dot = jet1.px()*jet2.px() + jet1.py()*jet2.py() + jet1.pz()*jet2.pz(); double norm1 = sqrt(jet1.px()*jet1.px() + jet1.py()*jet1.py() + jet1.pz()*jet1.pz()); double norm2 = sqrt(jet2.px()*jet2.px() + jet2.py()*jet2.py() + jet2.pz()*jet2.pz()); double costheta = dot/(norm1 * norm2); if (costheta > 1.0) costheta = 1.0; // Need to handle case of numerical overflow double theta = acos(costheta); pseudoRsquared = theta*theta; break; } case lorentz_dot : { double dotproduct = dot_product(jet1,jet2); pseudoRsquared = 2.0 * dotproduct / (jet1.e() * jet2.e()); break; } case perp_lorentz_dot : { PseudoJet lightJet = lightFrom(jet2); // assuming jet2 is the axis double dotproduct = dot_product(jet1,lightJet); pseudoRsquared = 2.0 * dotproduct / (lightJet.pt() * jet1.pt()); break; } default : { assert(_measure_type == pt_R || _measure_type == E_theta || _measure_type == lorentz_dot || _measure_type == perp_lorentz_dot); pseudoRsquared = std::numeric_limits::quiet_NaN(); break; } } return pseudoRsquared; } /////// // // Axes Refining // /////// // uses minimization of N-jettiness to continually update axes until convergence. // The function returns the axes found at the (local) minimum // This is the general axes refiner that can be used for a generic measure (but is // overwritten in the case of the conical measure and the deprecated geometric measure) std::vector MeasureDefinition::get_one_pass_axes(int n_jets, const std::vector & particles, const std::vector& currentAxes, int nAttempts, double accuracy) const { assert(n_jets == (int)currentAxes.size()); std::vector seedAxes = currentAxes; std::vector temp_axes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0)); for (unsigned int k = 0; k < seedAxes.size(); k++) { seedAxes[k] = lightFrom(seedAxes[k]) * seedAxes[k].E(); // making light-like, but keeping energy } double seedTau = result(particles, seedAxes); std::vector bestAxesSoFar = seedAxes; double bestTauSoFar = seedTau; for (int i_att = 0; i_att < nAttempts; i_att++) { std::vector newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0)); std::vector summed_jets(seedAxes.size(), fastjet::PseudoJet(0,0,0,0)); // find closest axis and assign to that for (unsigned int i = 0; i < particles.size(); i++) { // start from unclustered beam measure int minJ = -1; double minDist = beam_distance_squared(particles[i]); // which axis am I closest to? for (unsigned int j = 0; j < seedAxes.size(); j++) { double tempDist = jet_distance_squared(particles[i],seedAxes[j]); if (tempDist < minDist) { minDist = tempDist; minJ = j; } } // if not unclustered, then cluster if (minJ != -1) { summed_jets[minJ] += particles[i]; // keep track of energy to use later. if (_useAxisScaling) { double pseudoMomentum = dot_product(lightFrom(seedAxes[minJ]),particles[i]) + accuracy; // need small offset to avoid potential divide by zero issues double axis_scaling = (double)jet_numerator(particles[i], seedAxes[minJ])/pseudoMomentum; newAxes[minJ] += particles[i]*axis_scaling; } } } if (!_useAxisScaling) newAxes = summed_jets; // convert the axes to LightLike and then back to PseudoJet for (unsigned int k = 0; k < newAxes.size(); k++) { if (newAxes[k].perp() > 0) { newAxes[k] = lightFrom(newAxes[k]); newAxes[k] *= summed_jets[k].E(); // scale by energy to get sensible result } } // calculate tau on new axes double newTau = result(particles, newAxes); // find the smallest value of tau (and the corresponding axes) so far if (newTau < bestTauSoFar) { bestAxesSoFar = newAxes; bestTauSoFar = newTau; } if (fabs(newTau - seedTau) < accuracy) {// close enough for jazz seedAxes = newAxes; seedTau = newTau; break; } seedAxes = newAxes; seedTau = newTau; } // return the axes corresponding to the smallest tau found throughout all iterations // this is to prevent the minimization from returning a non-minimized of tau due to potential oscillations around the minimum return bestAxesSoFar; } // One pass minimization for the DefaultMeasure // Given starting axes, update to find better axes by using Kmeans clustering around the old axes template std::vector DefaultMeasure::UpdateAxesFast(const std::vector & old_axes, const std::vector & inputJets, double accuracy ) const { assert(old_axes.size() == N); // some storage, declared static to save allocation/re-allocation costs static LightLikeAxis new_axes[N]; static fastjet::PseudoJet new_jets[N]; for (int n = 0; n < N; ++n) { new_axes[n].reset(0.0,0.0,0.0,0.0); new_jets[n].reset_momentum(0.0,0.0,0.0,0.0); } double precision = accuracy; //TODO: actually cascade this in /////////////// Assignment Step ////////////////////////////////////////////////////////// std::vector assignment_index(inputJets.size()); int k_assign = -1; for (unsigned i = 0; i < inputJets.size(); i++){ double smallestDist = std::numeric_limits::max(); //large number for (int k = 0; k < N; k++) { double thisDist = old_axes[k].DistanceSq(inputJets[i]); if (thisDist < smallestDist) { smallestDist = thisDist; k_assign = k; } } if (smallestDist > sq(_Rcutoff)) {k_assign = -1;} assignment_index[i] = k_assign; } //////////////// Update Step ///////////////////////////////////////////////////////////// double distPhi, old_dist; for (unsigned i = 0; i < inputJets.size(); i++) { int old_jet_i = assignment_index[i]; if (old_jet_i == -1) {continue;} const fastjet::PseudoJet& inputJet_i = inputJets[i]; LightLikeAxis& new_axis_i = new_axes[old_jet_i]; double inputPhi_i = inputJet_i.phi(); double inputRap_i = inputJet_i.rap(); // optimize pow() call // add noise (the precision term) to make sure we don't divide by zero if (_beta == 1.0) { double DR = std::sqrt(sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i)); old_dist = 1.0/DR; } else if (_beta == 2.0) { old_dist = 1.0; } else if (_beta == 0.0) { double DRSq = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i); old_dist = 1.0/DRSq; } else { old_dist = sq(precision) + old_axes[old_jet_i].DistanceSq(inputJet_i); old_dist = std::pow(old_dist, (0.5*_beta-1.0)); } // TODO: Put some of these addition functions into light-like axes // rapidity sum new_axis_i.set_rap(new_axis_i.rap() + inputJet_i.perp() * inputRap_i * old_dist); // phi sum distPhi = inputPhi_i - old_axes[old_jet_i].phi(); if (fabs(distPhi) <= M_PI){ new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * inputPhi_i * old_dist ); } else if (distPhi > M_PI) { new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (-2*M_PI + inputPhi_i) * old_dist ); } else if (distPhi < -M_PI) { new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (+2*M_PI + inputPhi_i) * old_dist ); } // weights sum new_axis_i.set_weight( new_axis_i.weight() + inputJet_i.perp() * old_dist ); // momentum magnitude sum new_jets[old_jet_i] += inputJet_i; } // normalize sums for (int k = 0; k < N; k++) { if (new_axes[k].weight() == 0) { // no particles were closest to this axis! Return to old axis instead of (0,0,0,0) new_axes[k] = old_axes[k]; } else { new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() ); new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() ); new_axes[k].set_phi( std::fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) ); new_axes[k].set_mom( std::sqrt(new_jets[k].modp2()) ); } } std::vector new_axes_vec(N); for (unsigned k = 0; k < N; ++k) new_axes_vec[k] = new_axes[k]; return new_axes_vec; } // Given N starting axes, this function updates all axes to find N better axes. // (This is just a wrapper for the templated version above.) // TODO: Consider removing this in a future version std::vector DefaultMeasure::UpdateAxes(const std::vector & old_axes, const std::vector & inputJets, double accuracy) const { int N = old_axes.size(); switch (N) { case 1: return UpdateAxesFast<1>(old_axes, inputJets, accuracy); case 2: return UpdateAxesFast<2>(old_axes, inputJets, accuracy); case 3: return UpdateAxesFast<3>(old_axes, inputJets, accuracy); case 4: return UpdateAxesFast<4>(old_axes, inputJets, accuracy); case 5: return UpdateAxesFast<5>(old_axes, inputJets, accuracy); case 6: return UpdateAxesFast<6>(old_axes, inputJets, accuracy); case 7: return UpdateAxesFast<7>(old_axes, inputJets, accuracy); case 8: return UpdateAxesFast<8>(old_axes, inputJets, accuracy); case 9: return UpdateAxesFast<9>(old_axes, inputJets, accuracy); case 10: return UpdateAxesFast<10>(old_axes, inputJets, accuracy); case 11: return UpdateAxesFast<11>(old_axes, inputJets, accuracy); case 12: return UpdateAxesFast<12>(old_axes, inputJets, accuracy); case 13: return UpdateAxesFast<13>(old_axes, inputJets, accuracy); case 14: return UpdateAxesFast<14>(old_axes, inputJets, accuracy); case 15: return UpdateAxesFast<15>(old_axes, inputJets, accuracy); case 16: return UpdateAxesFast<16>(old_axes, inputJets, accuracy); case 17: return UpdateAxesFast<17>(old_axes, inputJets, accuracy); case 18: return UpdateAxesFast<18>(old_axes, inputJets, accuracy); case 19: return UpdateAxesFast<19>(old_axes, inputJets, accuracy); case 20: return UpdateAxesFast<20>(old_axes, inputJets, accuracy); default: std::cout << "N-jettiness is hard-coded to only allow up to 20 jets!" << std::endl; return std::vector(); } } // uses minimization of N-jettiness to continually update axes until convergence. // The function returns the axes found at the (local) minimum std::vector DefaultMeasure::get_one_pass_axes(int n_jets, const std::vector & inputJets, const std::vector& seedAxes, int nAttempts, double accuracy ) const { // if the measure type doesn't use the pt_R metric, then the standard minimization scheme should be used if (_measure_type != pt_R) { return MeasureDefinition::get_one_pass_axes(n_jets, inputJets, seedAxes, nAttempts, accuracy); } // convert from PseudoJets to LightLikeAxes std::vector< LightLikeAxis > old_axes(n_jets, LightLikeAxis(0,0,0,0)); for (int k = 0; k < n_jets; k++) { old_axes[k].set_rap( seedAxes[k].rap() ); old_axes[k].set_phi( seedAxes[k].phi() ); old_axes[k].set_mom( seedAxes[k].modp() ); } // Find new axes by iterating (only one pass here) std::vector< LightLikeAxis > new_axes(n_jets, LightLikeAxis(0,0,0,0)); double cmp = std::numeric_limits::max(); //large number int h = 0; while (cmp > accuracy && h < nAttempts) { // Keep updating axes until near-convergence or too many update steps cmp = 0.0; h++; new_axes = UpdateAxes(old_axes, inputJets,accuracy); // Update axes for (int k = 0; k < n_jets; k++) { cmp += old_axes[k].Distance(new_axes[k]); } cmp = cmp / ((double) n_jets); old_axes = new_axes; } // Convert from internal LightLikeAxes to PseudoJet std::vector outputAxes; for (int k = 0; k < n_jets; k++) { fastjet::PseudoJet temp = old_axes[k].ConvertToPseudoJet(); outputAxes.push_back(temp); } // this is used to debug the minimization routine to make sure that it works. bool do_debug = false; if (do_debug) { // get this information to make sure that minimization is working properly double seed_tau = result(inputJets, seedAxes); double outputTau = result(inputJets, outputAxes); assert(outputTau <= seed_tau); } return outputAxes; } //// One-pass minimization for the Deprecated Geometric Measure //// Uses minimization of the geometric distance in order to find the minimum axes. //// It continually updates until it reaches convergence or it reaches the maximum number of attempts. //// This is essentially the same as a stable cone finder. //std::vector DeprecatedGeometricCutoffMeasure::get_one_pass_axes(int n_jets, // const std::vector & particles, // const std::vector& currentAxes, // int nAttempts, // double accuracy) const { // // assert(n_jets == (int)currentAxes.size()); //added int casting to get rid of compiler warning // // std::vector seedAxes = currentAxes; // double seedTau = result(particles, seedAxes); // // for (int i = 0; i < nAttempts; i++) { // // std::vector newAxes(seedAxes.size(),fastjet::PseudoJet(0,0,0,0)); // // // find closest axis and assign to that // for (unsigned int i = 0; i < particles.size(); i++) { // // // start from unclustered beam measure // int minJ = -1; // double minDist = beam_distance_squared(particles[i]); // // // which axis am I closest to? // for (unsigned int j = 0; j < seedAxes.size(); j++) { // double tempDist = jet_distance_squared(particles[i],seedAxes[j]); // if (tempDist < minDist) { // minDist = tempDist; // minJ = j; // } // } // // // if not unclustered, then cluster // if (minJ != -1) newAxes[minJ] += particles[i]; // } // // // calculate tau on new axes // seedAxes = newAxes; // double tempTau = result(particles, newAxes); // // // close enough to stop? // if (fabs(tempTau - seedTau) < accuracy) break; // seedTau = tempTau; // } // // return seedAxes; //} // Go from internal LightLikeAxis to PseudoJet fastjet::PseudoJet LightLikeAxis::ConvertToPseudoJet() { double px, py, pz, E; E = _mom; pz = (std::exp(2.0*_rap) - 1.0) / (std::exp(2.0*_rap) + 1.0) * E; px = std::cos(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) ); py = std::sin(_phi) * std::sqrt( std::pow(E,2) - std::pow(pz,2) ); return fastjet::PseudoJet(px,py,pz,E); } } //namespace contrib FASTJET_END_NAMESPACE