| 1 |
// Dear emacs, this is -*- c++ -*-
|
| 2 |
|
| 3 |
#ifndef Njettiness_ICC
|
| 4 |
#define Njettiness_ICC
|
| 5 |
|
| 6 |
#include "Njettiness.h"
|
| 7 |
|
| 8 |
// Given starting axes, update to find better axes
|
| 9 |
template <int N>
|
| 10 |
std::vector<LightLikeAxis> Njettiness::UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
|
| 11 |
const std::vector <fastjet::PseudoJet> & inputJets,
|
| 12 |
NsubParameters paraNsub, double precision) {
|
| 13 |
assert(old_axes.size() == N);
|
| 14 |
|
| 15 |
// some storage, declared static to save allocation/re-allocation costs
|
| 16 |
static LightLikeAxis new_axes[N];
|
| 17 |
static fastjet::PseudoJet new_jets[N];
|
| 18 |
for (int n = 0; n < N; ++n) {
|
| 19 |
new_axes[n].reset(0.0,0.0,0.0,0.0);
|
| 20 |
#ifdef FASTJET2
|
| 21 |
new_jets[n].reset(0.0,0.0,0.0,0.0);
|
| 22 |
#else
|
| 23 |
// use cheaper reset if available
|
| 24 |
new_jets[n].reset_momentum(0.0,0.0,0.0,0.0);
|
| 25 |
#endif
|
| 26 |
}
|
| 27 |
|
| 28 |
|
| 29 |
double beta = paraNsub.beta();
|
| 30 |
double Rcutoff = paraNsub.Rcutoff();
|
| 31 |
|
| 32 |
/////////////// Assignment Step //////////////////////////////////////////////////////////
|
| 33 |
std::vector<int> assignment_index(inputJets.size());
|
| 34 |
int k_assign = -1;
|
| 35 |
|
| 36 |
for (unsigned i = 0; i < inputJets.size(); i++){
|
| 37 |
double smallestDist = 1000000.0;
|
| 38 |
for (int k = 0; k < N; k++) {
|
| 39 |
double thisDist = DistanceSq(inputJets[i].rap(),inputJets[i].phi(),old_axes[k].rap(),old_axes[k].phi());
|
| 40 |
if (thisDist < smallestDist) {
|
| 41 |
smallestDist = thisDist;
|
| 42 |
k_assign = k;
|
| 43 |
}
|
| 44 |
}
|
| 45 |
if (smallestDist > sq(Rcutoff)) {k_assign = -1;}
|
| 46 |
assignment_index[i] = k_assign;
|
| 47 |
}
|
| 48 |
|
| 49 |
//////////////// Update Step /////////////////////////////////////////////////////////////
|
| 50 |
double distPhi, old_dist;
|
| 51 |
for (unsigned i = 0; i < inputJets.size(); i++) {
|
| 52 |
int old_jet_i = assignment_index[i];
|
| 53 |
if (old_jet_i == -1) {continue;}
|
| 54 |
|
| 55 |
const fastjet::PseudoJet& inputJet_i = inputJets[i];
|
| 56 |
LightLikeAxis& new_axis_i = new_axes[old_jet_i];
|
| 57 |
double inputPhi_i = inputJet_i.phi();
|
| 58 |
double inputRap_i = inputJet_i.rap();
|
| 59 |
|
| 60 |
// optimize pow() call
|
| 61 |
// add noise (the precision term) to make sure we don't divide by zero
|
| 62 |
if (beta == 1.0) {
|
| 63 |
double DR = std::sqrt(sq(precision) + DistanceSq(inputRap_i, inputPhi_i, old_axes[old_jet_i].rap(), old_axes[old_jet_i].phi()));
|
| 64 |
old_dist = 1.0/DR;
|
| 65 |
} else if (beta == 2.0) {
|
| 66 |
old_dist = 1.0;
|
| 67 |
} else if (beta == 0.0) {
|
| 68 |
double DRSq = sq(precision) + DistanceSq(inputRap_i, inputPhi_i, old_axes[old_jet_i].rap(), old_axes[old_jet_i].phi());
|
| 69 |
old_dist = 1.0/DRSq;
|
| 70 |
} else {
|
| 71 |
old_dist = sq(precision) + DistanceSq(inputRap_i, inputPhi_i, old_axes[old_jet_i].rap(), old_axes[old_jet_i].phi());
|
| 72 |
old_dist = std::pow(old_dist, (0.5*beta-1.0));
|
| 73 |
}
|
| 74 |
|
| 75 |
// rapidity sum
|
| 76 |
new_axis_i.set_rap(new_axis_i.rap() + inputJet_i.perp() * inputRap_i * old_dist);
|
| 77 |
// phi sum
|
| 78 |
distPhi = inputPhi_i - old_axes[old_jet_i].phi();
|
| 79 |
if (fabs(distPhi) <= M_PI){
|
| 80 |
new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * inputPhi_i * old_dist );
|
| 81 |
} else if (distPhi > M_PI) {
|
| 82 |
new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (-2*M_PI + inputPhi_i) * old_dist );
|
| 83 |
} else if (distPhi < -M_PI) {
|
| 84 |
new_axis_i.set_phi( new_axis_i.phi() + inputJet_i.perp() * (+2*M_PI + inputPhi_i) * old_dist );
|
| 85 |
}
|
| 86 |
// weights sum
|
| 87 |
new_axis_i.set_weight( new_axis_i.weight() + inputJet_i.perp() * old_dist );
|
| 88 |
// momentum magnitude sum
|
| 89 |
new_jets[old_jet_i] += inputJet_i;
|
| 90 |
}
|
| 91 |
// normalize sums
|
| 92 |
for (int k = 0; k < N; k++) {
|
| 93 |
if (new_axes[k].weight() == 0) {
|
| 94 |
// no particles were closest to this axis! Return to old axis instead of (0,0,0,0)
|
| 95 |
new_axes[k] = old_axes[k];
|
| 96 |
} else {
|
| 97 |
new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() );
|
| 98 |
new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() );
|
| 99 |
new_axes[k].set_phi( std::fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) );
|
| 100 |
new_axes[k].set_mom( std::sqrt(new_jets[k].modp2()) );
|
| 101 |
}
|
| 102 |
}
|
| 103 |
std::vector<LightLikeAxis> new_axes_vec(N);
|
| 104 |
for (unsigned k = 0; k < N; ++k) new_axes_vec[k] = new_axes[k];
|
| 105 |
return new_axes_vec;
|
| 106 |
}
|
| 107 |
|
| 108 |
#endif
|
