| 1 |
rkogler |
1.1 |
// Njettiness Package
|
| 2 |
|
|
// Version 0.4.1 (January 22, 2012)
|
| 3 |
|
|
// Questions/Comments? jthaler@jthaler.net
|
| 4 |
|
|
|
| 5 |
|
|
// Copyright (c) 2011-12, Jesse Thaler, Ken Van Tilburg, and Christopher K.
|
| 6 |
|
|
// Vermilion
|
| 7 |
|
|
//
|
| 8 |
|
|
//----------------------------------------------------------------------
|
| 9 |
|
|
// This file is part of the N-jettiness package ("N-jettiness").
|
| 10 |
|
|
//
|
| 11 |
|
|
// This file has been modified to run with the SFrame analysis.
|
| 12 |
|
|
// It's original version comes from SpartyJet.
|
| 13 |
|
|
//----------------------------------------------------------------------
|
| 14 |
|
|
|
| 15 |
|
|
|
| 16 |
|
|
#ifndef __NJETTINESS_H__
|
| 17 |
|
|
#define __NJETTINESS_H__
|
| 18 |
|
|
|
| 19 |
|
|
#include "fastjet/PseudoJet.hh"
|
| 20 |
|
|
#include "fastjet/ClusterSequence.hh"
|
| 21 |
|
|
#include <cmath>
|
| 22 |
|
|
#include <vector>
|
| 23 |
|
|
#include <list>
|
| 24 |
|
|
|
| 25 |
|
|
///////
|
| 26 |
|
|
//
|
| 27 |
|
|
// Helper classes
|
| 28 |
|
|
//
|
| 29 |
|
|
///////
|
| 30 |
|
|
|
| 31 |
|
|
|
| 32 |
|
|
// Parameters that define Njettiness
|
| 33 |
|
|
class NsubParameters {
|
| 34 |
|
|
private:
|
| 35 |
|
|
double _beta; // angular weighting exponent
|
| 36 |
|
|
double _R0; // characteristic jet radius (for normalization)
|
| 37 |
|
|
double _Rcutoff; // Cutoff scale for cone jet finding (default is large number such that no boundaries are used)
|
| 38 |
|
|
|
| 39 |
|
|
public:
|
| 40 |
|
|
NsubParameters(const double mybeta, const double myR0, const double myRcutoff=10000.0) :
|
| 41 |
|
|
_beta(mybeta), _R0(myR0), _Rcutoff(myRcutoff) {}
|
| 42 |
|
|
double beta() const {return _beta;}
|
| 43 |
|
|
double R0() const {return _R0;}
|
| 44 |
|
|
double Rcutoff() const {return _Rcutoff;}
|
| 45 |
|
|
};
|
| 46 |
|
|
|
| 47 |
|
|
// Parameters that change minimization procedure.
|
| 48 |
|
|
// Set automatically when you choose NsubAxesMode, but can be adjusted manually as well
|
| 49 |
|
|
class KmeansParameters {
|
| 50 |
|
|
private:
|
| 51 |
|
|
int _n_iterations; // Number of iterations to run (0 for no minimization, 1 for one-pass, >>1 for global minimum)
|
| 52 |
|
|
double _precision; // Desired precision in axes alignment
|
| 53 |
|
|
int _halt; // maximum number of steps per iteration
|
| 54 |
|
|
double _noise_range;// noise range for random initialization
|
| 55 |
|
|
|
| 56 |
|
|
public:
|
| 57 |
|
|
KmeansParameters() : _n_iterations(0), _precision(0.0), _halt(0), _noise_range(0.0) {}
|
| 58 |
|
|
KmeansParameters(const int my_n_iterations, double my_precision, int my_halt, double my_noise_range) :
|
| 59 |
|
|
_n_iterations(my_n_iterations), _precision(my_precision), _halt(my_halt), _noise_range(my_noise_range) {}
|
| 60 |
|
|
int n_iterations() const { return _n_iterations;}
|
| 61 |
|
|
double precision() const {return _precision;}
|
| 62 |
|
|
int halt() const {return _halt;}
|
| 63 |
|
|
double noise_range() const {return _noise_range;}
|
| 64 |
|
|
};
|
| 65 |
|
|
|
| 66 |
|
|
// helper class for minimization
|
| 67 |
|
|
class LightLikeAxis {
|
| 68 |
|
|
private:
|
| 69 |
|
|
double _rap, _phi, _weight, _mom;
|
| 70 |
|
|
|
| 71 |
|
|
public:
|
| 72 |
|
|
LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
|
| 73 |
|
|
LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
|
| 74 |
|
|
_rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}
|
| 75 |
|
|
double rap() const {return _rap;}
|
| 76 |
|
|
double phi() const {return _phi;}
|
| 77 |
|
|
double weight() const {return _weight;}
|
| 78 |
|
|
double mom() const {return _mom;}
|
| 79 |
|
|
void set_rap(double my_set_rap) {_rap = my_set_rap;}
|
| 80 |
|
|
void set_phi(double my_set_phi) {_phi = my_set_phi;}
|
| 81 |
|
|
void set_weight(double my_set_weight) {_weight = my_set_weight;}
|
| 82 |
|
|
void set_mom(double my_set_mom) {_mom = my_set_mom;}
|
| 83 |
|
|
void reset(double my_rap, double my_phi, double my_weight, double my_mom) {_rap=my_rap; _phi=my_phi; _weight=my_weight; _mom=my_mom;}
|
| 84 |
|
|
};
|
| 85 |
|
|
|
| 86 |
|
|
|
| 87 |
|
|
///////
|
| 88 |
|
|
//
|
| 89 |
|
|
// Main Njettiness Class
|
| 90 |
|
|
//
|
| 91 |
|
|
///////
|
| 92 |
|
|
|
| 93 |
|
|
class Njettiness {
|
| 94 |
|
|
public:
|
| 95 |
|
|
enum AxesMode {
|
| 96 |
|
|
kt_axes, // exclusive kt axes
|
| 97 |
|
|
ca_axes, // exclusive ca axes
|
| 98 |
|
|
antikt_0p2_axes, // inclusive hardest axes with antikt-0.2
|
| 99 |
|
|
min_axes, // axes that minimize N-subjettiness (100 passes by default)
|
| 100 |
|
|
manual_axes, // set your own axes with setAxes()
|
| 101 |
|
|
onepass_kt_axes, // one-pass minimization from kt starting point
|
| 102 |
|
|
onepass_ca_axes, // one-pass minimization from ca starting point
|
| 103 |
|
|
onepass_antikt_0p2_axes, // one-pass minimization from antikt-0.2 starting point
|
| 104 |
|
|
onepass_manual_axes // one-pass minimization from manual starting point
|
| 105 |
|
|
};
|
| 106 |
|
|
|
| 107 |
|
|
private:
|
| 108 |
|
|
AxesMode _axes; //Axes mode choice
|
| 109 |
|
|
NsubParameters _paraNsub; //Parameters for Njettiness
|
| 110 |
|
|
KmeansParameters _paraKmeans; //Parameters for Minimization Procedure (set by NsubAxesMode automatically, but can change manually if desired)
|
| 111 |
|
|
|
| 112 |
|
|
std::vector<fastjet::PseudoJet> _currentAxes;
|
| 113 |
|
|
|
| 114 |
|
|
void establishAxes(unsigned n_jets, const std::vector <fastjet::PseudoJet> & inputs);
|
| 115 |
|
|
|
| 116 |
|
|
public:
|
| 117 |
|
|
Njettiness(AxesMode axes, NsubParameters paraNsub);
|
| 118 |
|
|
|
| 119 |
|
|
void setParaKmeans(KmeansParameters newPara) {_paraKmeans = newPara;}
|
| 120 |
|
|
void setParaNsub(NsubParameters newPara) {_paraNsub = newPara;}
|
| 121 |
|
|
|
| 122 |
|
|
// setAxes for Manual mode
|
| 123 |
|
|
void setAxes(std::vector<fastjet::PseudoJet> myAxes);
|
| 124 |
|
|
|
| 125 |
|
|
// The value of N-subjettiness
|
| 126 |
|
|
double getTau(unsigned n_jets, const std::vector<fastjet::PseudoJet> & inputJets);
|
| 127 |
|
|
|
| 128 |
|
|
// get axes used by getTau.
|
| 129 |
|
|
std::vector<fastjet::PseudoJet> currentAxes();
|
| 130 |
|
|
|
| 131 |
|
|
// partition inputs by Voronoi (each vector stores indices corresponding to inputJets)
|
| 132 |
|
|
std::vector<std::list<int> > getPartition(const std::vector<fastjet::PseudoJet> & inputJets);
|
| 133 |
|
|
|
| 134 |
|
|
// partition inputs by Voronoi
|
| 135 |
|
|
std::vector<fastjet::PseudoJet> getJets(const std::vector<fastjet::PseudoJet> & inputJets);
|
| 136 |
|
|
|
| 137 |
|
|
static double sq(double x);
|
| 138 |
|
|
|
| 139 |
|
|
// Calculates distance between two points in rapidity-azimuth plane
|
| 140 |
|
|
static double DistanceSq(double rap1, double phi1, double rap2, double phi2);
|
| 141 |
|
|
|
| 142 |
|
|
static double Distance(double rap1, double phi1, double rap2, double phi2);
|
| 143 |
|
|
|
| 144 |
|
|
template <int N> static std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
|
| 145 |
|
|
const std::vector <fastjet::PseudoJet> & inputJets,
|
| 146 |
|
|
NsubParameters paraNsub, double precision);
|
| 147 |
|
|
|
| 148 |
|
|
// Given starting axes, update to find better axes
|
| 149 |
|
|
// (This is just a wrapper for the templated version above.)
|
| 150 |
|
|
static std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
|
| 151 |
|
|
const std::vector <fastjet::PseudoJet> & inputJets, NsubParameters paraNsub, double precision);
|
| 152 |
|
|
|
| 153 |
|
|
|
| 154 |
|
|
// Go from internal LightLikeAxis to PseudoJet
|
| 155 |
|
|
static std::vector<fastjet::PseudoJet> ConvertToPseudoJet(const std::vector <LightLikeAxis>& axes);
|
| 156 |
|
|
|
| 157 |
|
|
// N-subjettiness pieces
|
| 158 |
|
|
// 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)
|
| 159 |
|
|
static std::vector<double> ConstituentTauValue(const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& axes, const NsubParameters& paraNsub);
|
| 160 |
|
|
|
| 161 |
|
|
// N-subjettiness values
|
| 162 |
|
|
// Calculates tau_N
|
| 163 |
|
|
static double TauValue(const std::vector <fastjet::PseudoJet>& particles, const std::vector<fastjet::PseudoJet>& axes,const NsubParameters& paraNsub);
|
| 164 |
|
|
|
| 165 |
|
|
// Get exclusive kT subjets
|
| 166 |
|
|
static std::vector<fastjet::PseudoJet> GetKTAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets);
|
| 167 |
|
|
|
| 168 |
|
|
// Get exclusive CA subjets
|
| 169 |
|
|
static std::vector<fastjet::PseudoJet> GetCAAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets);
|
| 170 |
|
|
|
| 171 |
|
|
// Get inclusive anti kT hardest subjets
|
| 172 |
|
|
static std::vector<fastjet::PseudoJet> GetAntiKTAxes(int n_jets, double R0, const std::vector <fastjet::PseudoJet> & inputJets);
|
| 173 |
|
|
|
| 174 |
|
|
// Get minimization axes
|
| 175 |
|
|
static std::vector<fastjet::PseudoJet> GetMinimumAxes(const std::vector <fastjet::PseudoJet> & seedAxes, const std::vector <fastjet::PseudoJet> & inputJets, KmeansParameters para,
|
| 176 |
|
|
NsubParameters paraNsub);
|
| 177 |
|
|
|
| 178 |
|
|
};
|
| 179 |
|
|
|
| 180 |
|
|
#endif
|
| 181 |
|
|
|
| 182 |
|
|
#ifndef __CINT__
|
| 183 |
|
|
#include "Njettiness.icc"
|
| 184 |
|
|
#endif
|
| 185 |
|
|
|
