| 1 |
ntran |
1.1.2.1 |
// Nsubjettiness Plugin
|
| 2 |
|
|
// Jesse Thaler and Ken Van Tilburg
|
| 3 |
|
|
// Version 0.3.2 (November 4, 2011)
|
| 4 |
|
|
// Questions/Comments? jthaler@jthaler.net
|
| 5 |
|
|
|
| 6 |
|
|
#ifndef __NJETTINESS_HH__
|
| 7 |
|
|
#define __NJETTINESS_HH__
|
| 8 |
|
|
|
| 9 |
|
|
|
| 10 |
|
|
#include "fastjet/PseudoJet.hh"
|
| 11 |
|
|
#include "fastjet/ClusterSequence.hh"
|
| 12 |
|
|
#include <cmath>
|
| 13 |
|
|
#include <vector>
|
| 14 |
|
|
#include <list>
|
| 15 |
|
|
|
| 16 |
|
|
///////
|
| 17 |
|
|
//
|
| 18 |
|
|
// Helper classes and enums
|
| 19 |
|
|
//
|
| 20 |
|
|
///////
|
| 21 |
|
|
|
| 22 |
|
|
// Choose Axes Mode
|
| 23 |
|
|
enum NsubAxesMode {
|
| 24 |
|
|
nsub_kt_axes, // exclusive kt axes
|
| 25 |
|
|
nsub_ca_axes, // exclusive ca axes
|
| 26 |
|
|
nsub_antikt_0p2_axes, // inclusive hardest axes with antikt-0.2
|
| 27 |
|
|
nsub_min_axes, // axes that minimize N-subjettiness (100 passes by default)
|
| 28 |
|
|
nsub_manual_axes, // set your own axes with setAxes()
|
| 29 |
|
|
nsub_1pass_from_kt_axes, // one-pass minimization from kt starting point
|
| 30 |
|
|
nsub_1pass_from_ca_axes, // one-pass minimization from ca starting point
|
| 31 |
|
|
nsub_1pass_from_antikt_0p2_axes, // one-pass minimization from antikt-0.2 starting point
|
| 32 |
|
|
nsub_1pass_from_manual_axes // one-pass minimization from manual starting point
|
| 33 |
|
|
};
|
| 34 |
|
|
|
| 35 |
|
|
// Parameters that define Nsubjettiness
|
| 36 |
|
|
class NsubParameters {
|
| 37 |
|
|
private:
|
| 38 |
|
|
double _beta; // angular weighting exponent
|
| 39 |
|
|
double _R0; // characteristic jet radius (for normalization)
|
| 40 |
|
|
double _Rcutoff; // Cutoff scale for cone jet finding (default is large number such that no boundaries are used)
|
| 41 |
|
|
|
| 42 |
|
|
public:
|
| 43 |
|
|
NsubParameters(double myBeta, double myR0, double myRcutoff=10000.0) :
|
| 44 |
|
|
_beta(myBeta), _R0(myR0), _Rcutoff(myRcutoff) {}
|
| 45 |
|
|
double beta() const {return _beta;}
|
| 46 |
|
|
double R0() const {return _R0;}
|
| 47 |
|
|
double Rcutoff() const {return _Rcutoff;}
|
| 48 |
|
|
};
|
| 49 |
|
|
|
| 50 |
|
|
// Parameters that change minimization procedure.
|
| 51 |
|
|
// Set automatically when you choose NsubAxesMode, but can be adjusted manually as well
|
| 52 |
|
|
class KmeansParameters {
|
| 53 |
|
|
private:
|
| 54 |
|
|
int _n_iterations; // Number of iterations to run (0 for no minimization, 1 for one-pass, >>1 for global minimum)
|
| 55 |
|
|
double _precision; // Desired precision in axes alignment
|
| 56 |
|
|
int _halt; // maximum number of steps per iteration
|
| 57 |
|
|
double _noise_range;// noise range for random initialization
|
| 58 |
|
|
|
| 59 |
|
|
public:
|
| 60 |
|
|
KmeansParameters() : _n_iterations(0), _precision(0.0), _halt(0.0), _noise_range(0.0) {}
|
| 61 |
|
|
KmeansParameters(int my_n_iterations, double my_precision, int my_halt, double my_noise_range) :
|
| 62 |
|
|
_n_iterations(my_n_iterations), _precision(my_precision), _halt(my_halt), _noise_range(my_noise_range) {}
|
| 63 |
|
|
int n_iterations() const { return _n_iterations;}
|
| 64 |
|
|
double precision() const {return _precision;}
|
| 65 |
|
|
int halt() const {return _halt;}
|
| 66 |
|
|
double noise_range() const {return _noise_range;}
|
| 67 |
|
|
};
|
| 68 |
|
|
|
| 69 |
|
|
// helper class for minimization
|
| 70 |
|
|
class LightLikeAxis {
|
| 71 |
|
|
private:
|
| 72 |
|
|
double _rap, _phi, _weight, _mom;
|
| 73 |
|
|
|
| 74 |
|
|
public:
|
| 75 |
|
|
LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
|
| 76 |
|
|
_rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}
|
| 77 |
|
|
double rap() {return _rap;}
|
| 78 |
|
|
double phi() {return _phi;}
|
| 79 |
|
|
double weight() {return _weight;}
|
| 80 |
|
|
double mom() {return _mom;}
|
| 81 |
|
|
void set_rap(double newRap) {_rap = newRap;}
|
| 82 |
|
|
void set_phi(double newPhi) {_phi = newPhi;}
|
| 83 |
|
|
void set_weight(double newWeight) {_weight = newWeight;}
|
| 84 |
|
|
void set_mom(double newMom) {_mom = newMom;}
|
| 85 |
|
|
};
|
| 86 |
|
|
|
| 87 |
|
|
///////
|
| 88 |
|
|
//
|
| 89 |
|
|
// Functions for minimization.
|
| 90 |
|
|
// TODO: Wrap these in N-subjettiness class
|
| 91 |
|
|
//
|
| 92 |
|
|
///////
|
| 93 |
|
|
|
| 94 |
|
|
// Calculates distance between two points in rapidity-azimuth plane
|
| 95 |
|
|
// TODO: Convert to built in fastjet DeltaR distance
|
| 96 |
|
|
double Distance(double rap1, double phi1, double rap2, double phi2){
|
| 97 |
|
|
double distRap = fabs(rap1-rap2); double distPhi = fabs(phi1-phi2);
|
| 98 |
|
|
if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
|
| 99 |
|
|
double distance = sqrt(distRap*distRap + distPhi*distPhi);
|
| 100 |
|
|
return distance;
|
| 101 |
|
|
}
|
| 102 |
|
|
|
| 103 |
|
|
// Given starting axes, update to find better axes
|
| 104 |
|
|
std::vector<LightLikeAxis> UpdateAxes(std::vector <LightLikeAxis> old_axes,
|
| 105 |
|
|
const std::vector <fastjet::PseudoJet> & inputJets, NsubParameters paraNsub) {
|
| 106 |
|
|
int n_jets = old_axes.size();
|
| 107 |
|
|
|
| 108 |
|
|
double beta = paraNsub.beta();
|
| 109 |
|
|
double Rcutoff = paraNsub.Rcutoff();
|
| 110 |
|
|
|
| 111 |
|
|
/////////////// Assignment Step //////////////////////////////////////////////////////////
|
| 112 |
|
|
std::vector<double> assignment_index(inputJets.size());
|
| 113 |
|
|
int k_assign = -1;
|
| 114 |
|
|
|
| 115 |
|
|
for (unsigned int i = 0; i < inputJets.size(); i++){
|
| 116 |
|
|
double smallestDist = 1000.0;
|
| 117 |
|
|
for (int k = 0; k < n_jets; k++) {
|
| 118 |
|
|
double thisDist = Distance(inputJets[i].rap(),inputJets[i].phi(),old_axes[k].rap(),old_axes[k].phi());
|
| 119 |
|
|
if (thisDist < smallestDist) {
|
| 120 |
|
|
smallestDist = thisDist;
|
| 121 |
|
|
k_assign = k;
|
| 122 |
|
|
}
|
| 123 |
|
|
}
|
| 124 |
|
|
if (smallestDist > Rcutoff) {k_assign = -1;}
|
| 125 |
|
|
assignment_index[i] = k_assign;
|
| 126 |
|
|
}
|
| 127 |
|
|
|
| 128 |
|
|
//////////////// Update Step /////////////////////////////////////////////////////////////
|
| 129 |
|
|
std::vector <LightLikeAxis> new_axes(n_jets,LightLikeAxis(0,0,0,0));
|
| 130 |
|
|
std::vector< std::vector<double> > fourVecJets(4, std::vector<double>(n_jets,0));
|
| 131 |
|
|
double distPhi, old_dist;
|
| 132 |
|
|
for (unsigned int i = 0; i < inputJets.size(); i++){
|
| 133 |
|
|
if (assignment_index[i] == -1) {continue;}
|
| 134 |
|
|
old_dist = Distance(inputJets[i].rap(),inputJets[i].phi(),old_axes[assignment_index[i]].rap(), old_axes[assignment_index[i]].phi());
|
| 135 |
|
|
old_dist = pow(old_dist, (beta-2));
|
| 136 |
|
|
// rapidity sum
|
| 137 |
|
|
new_axes[assignment_index[i]].set_rap(new_axes[assignment_index[i]].rap() + inputJets[i].perp() * inputJets[i].rap() * old_dist);
|
| 138 |
|
|
// phi sum
|
| 139 |
|
|
distPhi = inputJets[i].phi() - old_axes[assignment_index[i]].phi();
|
| 140 |
|
|
if (fabs(distPhi) <= M_PI){
|
| 141 |
|
|
new_axes[assignment_index[i]].set_phi( new_axes[assignment_index[i]].phi() + inputJets[i].perp() * inputJets[i].phi() * old_dist );
|
| 142 |
|
|
} else if (distPhi > M_PI) {
|
| 143 |
|
|
new_axes[assignment_index[i]].set_phi( new_axes[assignment_index[i]].phi() + inputJets[i].perp() * (-2*M_PI + inputJets[i].phi()) * old_dist );
|
| 144 |
|
|
} else if (distPhi < -M_PI) {
|
| 145 |
|
|
new_axes[assignment_index[i]].set_phi( new_axes[assignment_index[i]].phi() + inputJets[i].perp() * (+2*M_PI + inputJets[i].phi()) * old_dist );
|
| 146 |
|
|
}
|
| 147 |
|
|
// weights sum
|
| 148 |
|
|
new_axes[assignment_index[i]].set_weight( new_axes[assignment_index[i]].weight() + inputJets[i].perp() * old_dist );
|
| 149 |
|
|
// momentum magnitude sum
|
| 150 |
|
|
fourVecJets[0][assignment_index[i]] += inputJets[i].px();
|
| 151 |
|
|
fourVecJets[1][assignment_index[i]] += inputJets[i].py();
|
| 152 |
|
|
fourVecJets[2][assignment_index[i]] += inputJets[i].pz();
|
| 153 |
|
|
fourVecJets[3][assignment_index[i]] += inputJets[i].e();
|
| 154 |
|
|
}
|
| 155 |
|
|
// normalize sums
|
| 156 |
|
|
for (int k = 0; k < n_jets; k++){
|
| 157 |
|
|
if (new_axes[k].weight() == 0) {continue;}
|
| 158 |
|
|
new_axes[k].set_rap( new_axes[k].rap() / new_axes[k].weight() );
|
| 159 |
|
|
new_axes[k].set_phi( new_axes[k].phi() / new_axes[k].weight() );
|
| 160 |
|
|
new_axes[k].set_phi( fmod(new_axes[k].phi() + 2*M_PI, 2*M_PI) );
|
| 161 |
|
|
new_axes[k].set_mom( sqrt( pow(fourVecJets[0][k],2) + pow(fourVecJets[1][k],2) + pow(fourVecJets[2][k],2) ) );
|
| 162 |
|
|
}
|
| 163 |
|
|
return new_axes;
|
| 164 |
|
|
}
|
| 165 |
|
|
|
| 166 |
|
|
// Go from internal LightLikeAxis to PseudoJet
|
| 167 |
|
|
// TODO: Make part of LightLikeAxis class.
|
| 168 |
|
|
std::vector<fastjet::PseudoJet> ConvertToPseudoJet(std::vector <LightLikeAxis> axes) {
|
| 169 |
|
|
|
| 170 |
|
|
int n_jets = axes.size();
|
| 171 |
|
|
|
| 172 |
|
|
double px, py, pz, E;
|
| 173 |
|
|
std::vector<fastjet::PseudoJet> FourVecJets;
|
| 174 |
|
|
for (int k = 0; k < n_jets; k++) {
|
| 175 |
|
|
E = axes[k].mom();
|
| 176 |
|
|
pz = (exp(2.0*axes[k].rap()) - 1.0) / (exp(2.0*axes[k].rap()) + 1.0) * E;
|
| 177 |
|
|
px = cos(axes[k].phi()) * sqrt( pow(E,2) - pow(pz,2) );
|
| 178 |
|
|
py = sin(axes[k].phi()) * sqrt( pow(E,2) - pow(pz,2) );
|
| 179 |
|
|
fastjet::PseudoJet temp = fastjet::PseudoJet(px,py,pz,E);
|
| 180 |
|
|
FourVecJets.push_back(temp);
|
| 181 |
|
|
}
|
| 182 |
|
|
return FourVecJets;
|
| 183 |
|
|
}
|
| 184 |
|
|
|
| 185 |
|
|
// N-subjettiness pieces
|
| 186 |
|
|
std::vector<double> ConstituentTauValue(const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& axes,const NsubParameters paraNsub) {// Returns the sub-tau values, i.e. a std::vector of the contributions to tau_N of each Voronoi region (or region within R_0)
|
| 187 |
|
|
double beta = paraNsub.beta();
|
| 188 |
|
|
double R0 = paraNsub.R0();
|
| 189 |
|
|
double Rcutoff = paraNsub.Rcutoff();
|
| 190 |
|
|
|
| 191 |
|
|
std::vector<double> tauNum(axes.size(),0.0),tau(axes.size());
|
| 192 |
|
|
double tauDen = 0.0;
|
| 193 |
|
|
for (unsigned int i = 0; i < particles.size(); i++) {
|
| 194 |
|
|
// find minimum distance
|
| 195 |
|
|
int j_min = -1;
|
| 196 |
|
|
double minR = 10000.0; // large number
|
| 197 |
|
|
for (unsigned int j = 0; j < axes.size(); j++) {
|
| 198 |
|
|
double tempR = sqrt(particles[i].squared_distance(axes[j])); // delta R distance
|
| 199 |
|
|
if (tempR < minR) {minR = tempR; j_min = j;}
|
| 200 |
|
|
}
|
| 201 |
|
|
if (minR > Rcutoff) {minR = Rcutoff;}
|
| 202 |
|
|
tauNum[j_min] += particles[i].perp() * pow(minR,beta);
|
| 203 |
|
|
tauDen += particles[i].perp() * pow(R0,beta);
|
| 204 |
|
|
}
|
| 205 |
|
|
for (unsigned int j = 0; j < axes.size(); j++) {
|
| 206 |
|
|
tau[j] = tauNum[j]/tauDen;
|
| 207 |
|
|
}
|
| 208 |
|
|
return tau;
|
| 209 |
|
|
}
|
| 210 |
|
|
|
| 211 |
|
|
// N-subjettiness values
|
| 212 |
|
|
double TauValue(const std::vector <fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes,const NsubParameters paraNsub) {// Calculates tau_N
|
| 213 |
|
|
std::vector<double> tau_vec = ConstituentTauValue(particles, axes,paraNsub);
|
| 214 |
|
|
double tau = 0.0;
|
| 215 |
|
|
for (unsigned int j = 0; j < tau_vec.size(); j++) {tau += tau_vec[j];}
|
| 216 |
|
|
return tau;
|
| 217 |
|
|
}
|
| 218 |
|
|
|
| 219 |
|
|
// Get exclusive kT subjets
|
| 220 |
|
|
std::vector<fastjet::PseudoJet> GetKTAxes(const int n_jets, const std::vector <fastjet::PseudoJet> & inputJets) {
|
| 221 |
|
|
fastjet::JetDefinition jet_def = fastjet::JetDefinition(fastjet::kt_algorithm,M_PI/2.0,fastjet::E_scheme,fastjet::Best);
|
| 222 |
|
|
fastjet::ClusterSequence jet_clust_seq(inputJets, jet_def);
|
| 223 |
|
|
return jet_clust_seq.exclusive_jets(n_jets);
|
| 224 |
|
|
}
|
| 225 |
|
|
|
| 226 |
|
|
// Get exclusive CA subjets
|
| 227 |
|
|
std::vector<fastjet::PseudoJet> GetCAAxes(const int n_jets, const std::vector <fastjet::PseudoJet> & inputJets) {
|
| 228 |
|
|
fastjet::JetDefinition jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,M_PI/2.0,fastjet::E_scheme,fastjet::Best);
|
| 229 |
|
|
fastjet::ClusterSequence jet_clust_seq(inputJets, jet_def);
|
| 230 |
|
|
return jet_clust_seq.exclusive_jets(n_jets);
|
| 231 |
|
|
}
|
| 232 |
|
|
|
| 233 |
|
|
// Get inclusive anti kT hardest subjets
|
| 234 |
|
|
std::vector<fastjet::PseudoJet> GetAntiKTAxes(const int n_jets, const double R0, const std::vector <fastjet::PseudoJet> & inputJets) {
|
| 235 |
|
|
fastjet::JetDefinition jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,R0,fastjet::E_scheme,fastjet::Best);
|
| 236 |
|
|
fastjet::ClusterSequence jet_clust_seq(inputJets, jet_def);
|
| 237 |
|
|
std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
|
| 238 |
|
|
myJets.resize(n_jets); // only keep n hardest
|
| 239 |
|
|
return myJets;
|
| 240 |
|
|
}
|
| 241 |
|
|
|
| 242 |
|
|
|
| 243 |
|
|
// Get minimization axes
|
| 244 |
|
|
std::vector<fastjet::PseudoJet> GetMinimumAxes(std::vector <fastjet::PseudoJet> seedAxes, const std::vector <fastjet::PseudoJet> & inputJets, KmeansParameters para,
|
| 245 |
|
|
NsubParameters paraNsub) {
|
| 246 |
|
|
int n_jets = seedAxes.size();
|
| 247 |
|
|
double noise = 0, tau = 10000.0, tau_tmp, cmp;
|
| 248 |
|
|
std::vector< LightLikeAxis > new_axes(n_jets, LightLikeAxis(0,0,0,0)), old_axes(n_jets, LightLikeAxis(0,0,0,0));
|
| 249 |
|
|
std::vector<fastjet::PseudoJet> tmp_min_axes, min_axes;
|
| 250 |
|
|
for (int l = 0; l < para.n_iterations(); l++) { // Do minimization procedure multiple times
|
| 251 |
|
|
// Add noise to guess for the axes
|
| 252 |
|
|
for (int k = 0; k < n_jets; k++) {
|
| 253 |
|
|
if ( 0 == l ) { // Don't add noise on first pass
|
| 254 |
|
|
old_axes[k].set_rap( seedAxes[k].rap() + noise );
|
| 255 |
|
|
old_axes[k].set_phi( seedAxes[k].phi() + noise );
|
| 256 |
|
|
} else {
|
| 257 |
|
|
noise = ((double)rand()/(double)RAND_MAX) * para.noise_range() * 2 - para.noise_range();
|
| 258 |
|
|
old_axes[k].set_rap( seedAxes[k].rap() + noise );
|
| 259 |
|
|
noise = ((double)rand()/(double)RAND_MAX) * para.noise_range() * 2 - para.noise_range();
|
| 260 |
|
|
old_axes[k].set_phi( seedAxes[k].phi() + noise );
|
| 261 |
|
|
}
|
| 262 |
|
|
}
|
| 263 |
|
|
cmp = 100.0; int h = 0;
|
| 264 |
|
|
while (cmp > para.precision() && h < para.halt()) { // Keep updating axes until near-convergence or too many update steps
|
| 265 |
|
|
cmp = 0.0; h++;
|
| 266 |
|
|
new_axes = UpdateAxes(old_axes, inputJets, paraNsub); // Update axes
|
| 267 |
|
|
for (int k = 0; k < n_jets; k++) {
|
| 268 |
|
|
cmp += Distance(new_axes[k].rap(),new_axes[k].phi(),old_axes[k].rap(),old_axes[k].phi());
|
| 269 |
|
|
}
|
| 270 |
|
|
cmp = cmp / ((double) n_jets);
|
| 271 |
|
|
old_axes = new_axes;
|
| 272 |
|
|
}
|
| 273 |
|
|
tmp_min_axes = ConvertToPseudoJet(old_axes); // Convert axes directions into four-std::vectors
|
| 274 |
|
|
tau_tmp = TauValue(inputJets, tmp_min_axes,paraNsub);
|
| 275 |
|
|
if (tau_tmp < tau) {tau = tau_tmp; min_axes = tmp_min_axes;} // Keep axes and tau only if they are best so far
|
| 276 |
|
|
}
|
| 277 |
|
|
return min_axes;
|
| 278 |
|
|
}
|
| 279 |
|
|
|
| 280 |
|
|
///////
|
| 281 |
|
|
//
|
| 282 |
|
|
// Main Nsubjettiness Class
|
| 283 |
|
|
//
|
| 284 |
|
|
///////
|
| 285 |
|
|
|
| 286 |
|
|
class Nsubjettiness {
|
| 287 |
|
|
|
| 288 |
|
|
private:
|
| 289 |
|
|
NsubAxesMode _axes; //Axes mode choice
|
| 290 |
|
|
NsubParameters _paraNsub; //Parameters for Nsubjettiness
|
| 291 |
|
|
KmeansParameters _paraKmeans; //Parameters for Minimization Procedure (set by NsubAxesMode automatically, but can change manually if desired)
|
| 292 |
|
|
|
| 293 |
|
|
std::vector<fastjet::PseudoJet> _currentAxes;
|
| 294 |
|
|
|
| 295 |
|
|
void establishAxes(const int n_jets, const std::vector <fastjet::PseudoJet> & inputs);
|
| 296 |
|
|
|
| 297 |
|
|
public:
|
| 298 |
|
|
Nsubjettiness(NsubAxesMode axes, NsubParameters paraNsub);
|
| 299 |
|
|
|
| 300 |
|
|
void setParaKmeans(KmeansParameters newPara) {_paraKmeans = newPara;}
|
| 301 |
|
|
void setParaNsub(NsubParameters newPara) {_paraNsub = newPara;}
|
| 302 |
|
|
|
| 303 |
|
|
// setAxes for Manual mode
|
| 304 |
|
|
void setAxes(std::vector<fastjet::PseudoJet> myAxes) {
|
| 305 |
|
|
assert((_axes == nsub_manual_axes) || (_axes == nsub_1pass_from_manual_axes));
|
| 306 |
|
|
_currentAxes = myAxes;
|
| 307 |
|
|
}
|
| 308 |
|
|
|
| 309 |
|
|
// The value of N-subjettiness
|
| 310 |
|
|
double getTau(const int n_jets, const std::vector<fastjet::PseudoJet> & inputJets) {
|
| 311 |
|
|
if (n_jets >= (int) inputJets.size()) {
|
| 312 |
|
|
_currentAxes = inputJets;
|
| 313 |
|
|
_currentAxes.resize(n_jets,fastjet::PseudoJet(0.0,0.0,0.0,0.0));
|
| 314 |
|
|
return 0.0;
|
| 315 |
|
|
}
|
| 316 |
|
|
establishAxes(n_jets, inputJets); // sets current Axes
|
| 317 |
|
|
return TauValue(inputJets,_currentAxes,_paraNsub);
|
| 318 |
|
|
}
|
| 319 |
|
|
|
| 320 |
|
|
// get axes used by getTau.
|
| 321 |
|
|
std::vector<fastjet::PseudoJet> currentAxes() {
|
| 322 |
|
|
return _currentAxes;
|
| 323 |
|
|
}
|
| 324 |
|
|
|
| 325 |
|
|
// partition inputs by Voronoi (each vector stores indices corresponding to inputJets)
|
| 326 |
|
|
std::vector<std::list<int> > getPartition(const std::vector<fastjet::PseudoJet> & inputJets);
|
| 327 |
|
|
|
| 328 |
|
|
// partition inputs by Voronoi
|
| 329 |
|
|
std::vector<fastjet::PseudoJet> getJets(const std::vector<fastjet::PseudoJet> & inputJets);
|
| 330 |
|
|
|
| 331 |
|
|
};
|
| 332 |
|
|
|
| 333 |
|
|
|
| 334 |
|
|
//Use NsubAxesMode to pick which type of axes to use
|
| 335 |
|
|
void Nsubjettiness::establishAxes(const int n_jets, const std::vector <fastjet::PseudoJet> & inputs) {
|
| 336 |
|
|
switch (_axes) {
|
| 337 |
|
|
case nsub_kt_axes:
|
| 338 |
|
|
_currentAxes = GetKTAxes(n_jets,inputs);
|
| 339 |
|
|
break;
|
| 340 |
|
|
case nsub_ca_axes:
|
| 341 |
|
|
_currentAxes = GetCAAxes(n_jets,inputs);
|
| 342 |
|
|
break;
|
| 343 |
|
|
case nsub_antikt_0p2_axes:
|
| 344 |
|
|
_currentAxes = GetAntiKTAxes(n_jets,0.2,inputs);
|
| 345 |
|
|
break;
|
| 346 |
|
|
case nsub_1pass_from_kt_axes:
|
| 347 |
|
|
case nsub_min_axes:
|
| 348 |
|
|
_currentAxes = GetKTAxes(n_jets,inputs);
|
| 349 |
|
|
_currentAxes = GetMinimumAxes(_currentAxes, inputs, _paraKmeans, _paraNsub);
|
| 350 |
|
|
break;
|
| 351 |
|
|
case nsub_1pass_from_ca_axes:
|
| 352 |
|
|
_currentAxes = GetCAAxes(n_jets,inputs);
|
| 353 |
|
|
_currentAxes = GetMinimumAxes(_currentAxes, inputs, _paraKmeans, _paraNsub);
|
| 354 |
|
|
break;
|
| 355 |
|
|
case nsub_1pass_from_antikt_0p2_axes:
|
| 356 |
|
|
_currentAxes = GetAntiKTAxes(n_jets,0.2,inputs);
|
| 357 |
|
|
_currentAxes = GetMinimumAxes(_currentAxes, inputs, _paraKmeans, _paraNsub);
|
| 358 |
|
|
break;
|
| 359 |
|
|
case nsub_1pass_from_manual_axes:
|
| 360 |
|
|
assert((int) _currentAxes.size() == n_jets);
|
| 361 |
|
|
_currentAxes = GetMinimumAxes(_currentAxes, inputs, _paraKmeans, _paraNsub);
|
| 362 |
|
|
break;
|
| 363 |
|
|
case nsub_manual_axes:
|
| 364 |
|
|
assert((int) _currentAxes.size() == n_jets);
|
| 365 |
|
|
break;
|
| 366 |
|
|
default:
|
| 367 |
|
|
assert(false);
|
| 368 |
|
|
break;
|
| 369 |
|
|
}
|
| 370 |
|
|
}
|
| 371 |
|
|
|
| 372 |
|
|
//Constructor sets KmeansParameters from NsubAxesMode input
|
| 373 |
|
|
Nsubjettiness::Nsubjettiness(NsubAxesMode axes, NsubParameters paraNsub) : _axes(axes), _paraNsub(paraNsub), _paraKmeans() {
|
| 374 |
|
|
switch (_axes) {
|
| 375 |
|
|
case nsub_kt_axes:
|
| 376 |
|
|
case nsub_ca_axes:
|
| 377 |
|
|
case nsub_antikt_0p2_axes:
|
| 378 |
|
|
_paraKmeans = KmeansParameters(0,0.0,0.0,0.0);
|
| 379 |
|
|
break;
|
| 380 |
|
|
case nsub_1pass_from_kt_axes:
|
| 381 |
|
|
case nsub_1pass_from_ca_axes:
|
| 382 |
|
|
case nsub_1pass_from_antikt_0p2_axes:
|
| 383 |
|
|
case nsub_1pass_from_manual_axes:
|
| 384 |
|
|
_paraKmeans = KmeansParameters(1,0.0001,1000.0,0.8);
|
| 385 |
|
|
break;
|
| 386 |
|
|
case nsub_min_axes:
|
| 387 |
|
|
_paraKmeans = KmeansParameters(100,0.0001,1000.0,0.8);
|
| 388 |
|
|
break;
|
| 389 |
|
|
case nsub_manual_axes:
|
| 390 |
|
|
_paraKmeans = KmeansParameters(0,0.0,0.0,0.0);
|
| 391 |
|
|
break;
|
| 392 |
|
|
default:
|
| 393 |
|
|
assert(false);
|
| 394 |
|
|
break;
|
| 395 |
|
|
}
|
| 396 |
|
|
}
|
| 397 |
|
|
|
| 398 |
|
|
|
| 399 |
|
|
std::vector<std::list<int> > Nsubjettiness::getPartition(const std::vector<fastjet::PseudoJet> & particles) {
|
| 400 |
|
|
double Rcutoff = _paraNsub.Rcutoff();
|
| 401 |
|
|
|
| 402 |
|
|
std::vector<std::list<int> > partitions(_currentAxes.size());
|
| 403 |
|
|
|
| 404 |
|
|
for (unsigned int i = 0; i < particles.size(); i++) {
|
| 405 |
|
|
// find minimum distance
|
| 406 |
|
|
int j_min = -1;
|
| 407 |
|
|
double minR = 10000.0; // large number
|
| 408 |
|
|
for (unsigned int j = 0; j < _currentAxes.size(); j++) {
|
| 409 |
|
|
double tempR = sqrt(particles[i].squared_distance(_currentAxes[j])); // delta R distance
|
| 410 |
|
|
if (tempR < minR) {
|
| 411 |
|
|
minR = tempR;
|
| 412 |
|
|
j_min = j;
|
| 413 |
|
|
}
|
| 414 |
|
|
}
|
| 415 |
|
|
if (minR > Rcutoff) {
|
| 416 |
|
|
// do nothing
|
| 417 |
|
|
} else {
|
| 418 |
|
|
partitions[j_min].push_back(i);
|
| 419 |
|
|
}
|
| 420 |
|
|
}
|
| 421 |
|
|
return partitions;
|
| 422 |
|
|
}
|
| 423 |
|
|
|
| 424 |
|
|
// TODO: Make this call getPartition
|
| 425 |
|
|
std::vector<fastjet::PseudoJet> Nsubjettiness::getJets(const std::vector<fastjet::PseudoJet> & particles) {
|
| 426 |
|
|
double Rcutoff = _paraNsub.Rcutoff();
|
| 427 |
|
|
|
| 428 |
|
|
std::vector<fastjet::PseudoJet> jets(_currentAxes.size());
|
| 429 |
|
|
|
| 430 |
|
|
for (unsigned int i = 0; i < particles.size(); i++) {
|
| 431 |
|
|
// find minimum distance
|
| 432 |
|
|
int j_min = -1;
|
| 433 |
|
|
double minR = 10000.0; // large number
|
| 434 |
|
|
for (unsigned int j = 0; j < _currentAxes.size(); j++) {
|
| 435 |
|
|
double tempR = sqrt(particles[i].squared_distance(_currentAxes[j])); // delta R distance
|
| 436 |
|
|
if (tempR < minR) {
|
| 437 |
|
|
minR = tempR;
|
| 438 |
|
|
j_min = j;
|
| 439 |
|
|
}
|
| 440 |
|
|
}
|
| 441 |
|
|
if (minR > Rcutoff) {
|
| 442 |
|
|
// do nothing
|
| 443 |
|
|
} else {
|
| 444 |
|
|
jets[j_min] += particles[i];
|
| 445 |
|
|
}
|
| 446 |
|
|
}
|
| 447 |
|
|
return jets;
|
| 448 |
|
|
}
|
| 449 |
|
|
|
| 450 |
|
|
#endif
|
| 451 |
|
|
|
