Page MenuHomeHEPForge
Files F7878320

No OneTemporary
Actions

View Options

Index: contrib/contribs/EnergyCorrelator/trunk/EnergyCorrelator.cc
===================================================================
--- contrib/contribs/EnergyCorrelator/trunk/EnergyCorrelator.cc (revision 1096)
+++ contrib/contribs/EnergyCorrelator/trunk/EnergyCorrelator.cc (revision 1097)
@@ -1,1081 +1,1116 @@
// EnergyCorrelator Package
// Questions/Comments? Email the authors:
// larkoski@mit.edu, lnecib@mit.edu,
// gavin.salam@cern.ch jthaler@jthaler.net
//
// Copyright (c) 2013-2016
// Andrew Larkoski, Lina Necib, Gavin Salam, and Jesse Thaler
//
// $Id$
//----------------------------------------------------------------------
// 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 <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#include "EnergyCorrelator.hh"
#include <sstream>
#include <limits>
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib {
double EnergyCorrelator::result(const PseudoJet& jet) const {
// if jet does not have constituents, throw error
if (!jet.has_constituents()) throw Error("EnergyCorrelator called on jet with no constituents.");
// get N = 0 case out of the way
if (_N == 0) return 1.0;
// find constituents
std::vector<fastjet::PseudoJet> particles = jet.constituents();
// return zero if the number of constituents is less than _N
if (particles.size() < _N) return 0.0 ;
double answer = 0.0;
// take care of N = 1 case.
if (_N == 1) {
for (unsigned int i = 0; i < particles.size(); i++) {
answer += energy(particles[i]);
}
return answer;
}
double half_beta = _beta/2.0;
// take care of N = 2 case.
if (_N == 2) {
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) { //note offset by one so that angle is never called on identical pairs
answer += energy(particles[i])
* energy(particles[j])
* pow(angleSquared(particles[i],particles[j]), half_beta);
}
}
return answer;
}
// if N > 5, then throw error
if (_N > 5) {
throw Error("EnergyCorrelator is only hard coded for N = 0,1,2,3,4,5");
}
// Now deal with N = 3,4,5. Different options if storage array is used or not.
if (_strategy == storage_array) {
// For N > 2, fill static storage array to save computation time.
// Make energy storage
std::vector<double> energyStore;
energyStore.resize(particles.size());
// Make angular storage
std::vector< std::vector<double> > angleStore;
angleStore.resize(particles.size());
for (unsigned int i = 0; i < angleStore.size(); i++) {
angleStore[i].resize(i);
}
// Fill storage with energy/angle information
for (unsigned int i = 0; i < particles.size(); i++) {
energyStore[i] = energy(particles[i]);
for (unsigned int j = 0; j < i; j++) {
angleStore[i][j] = pow(angleSquared(particles[i],particles[j]), half_beta);
}
}
// now do recursion
if (_N == 3) {
for (unsigned int i = 2; i < particles.size(); i++) {
for (unsigned int j = 1; j < i; j++) {
double ans_ij = energyStore[i]
* energyStore[j]
* angleStore[i][j];
for (unsigned int k = 0; k < j; k++) {
answer += ans_ij
* energyStore[k]
* angleStore[i][k]
* angleStore[j][k];
}
}
}
} else if (_N == 4) {
for (unsigned int i = 3; i < particles.size(); i++) {
for (unsigned int j = 2; j < i; j++) {
double ans_ij = energyStore[i]
* energyStore[j]
* angleStore[i][j];
for (unsigned int k = 1; k < j; k++) {
double ans_ijk = ans_ij
* energyStore[k]
* angleStore[i][k]
* angleStore[j][k];
for (unsigned int l = 0; l < k; l++) {
answer += ans_ijk
* energyStore[l]
* angleStore[i][l]
* angleStore[j][l]
* angleStore[k][l];
}
}
}
}
} else if (_N == 5) {
for (unsigned int i = 4; i < particles.size(); i++) {
for (unsigned int j = 3; j < i; j++) {
double ans_ij = energyStore[i]
* energyStore[j]
* angleStore[i][j];
for (unsigned int k = 2; k < j; k++) {
double ans_ijk = ans_ij
* energyStore[k]
* angleStore[i][k]
* angleStore[j][k];
for (unsigned int l = 1; l < k; l++) {
double ans_ijkl = ans_ijk
* energyStore[l]
* angleStore[i][l]
* angleStore[j][l]
* angleStore[k][l];
for (unsigned int m = 0; m < l; m++) {
answer += ans_ijkl
* energyStore[m]
* angleStore[i][m]
* angleStore[j][m]
* angleStore[k][m]
* angleStore[l][m];
}
}
}
}
}
} else {
assert(_N <= 5);
}
} else if (_strategy == slow) {
if (_N == 3) {
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
double ans_ij = energy(particles[i])
* energy(particles[j])
* pow(angleSquared(particles[i],particles[j]), half_beta);
for (unsigned int k = j + 1; k < particles.size(); k++) {
answer += ans_ij
* energy(particles[k])
* pow(angleSquared(particles[i],particles[k]), half_beta)
* pow(angleSquared(particles[j],particles[k]), half_beta);
}
}
}
} else if (_N == 4) {
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
double ans_ij = energy(particles[i])
* energy(particles[j])
* pow(angleSquared(particles[i],particles[j]), half_beta);
for (unsigned int k = j + 1; k < particles.size(); k++) {
double ans_ijk = ans_ij
* energy(particles[k])
* pow(angleSquared(particles[i],particles[k]), half_beta)
* pow(angleSquared(particles[j],particles[k]), half_beta);
for (unsigned int l = k + 1; l < particles.size(); l++) {
answer += ans_ijk
* energy(particles[l])
* pow(angleSquared(particles[i],particles[l]), half_beta)
* pow(angleSquared(particles[j],particles[l]), half_beta)
* pow(angleSquared(particles[k],particles[l]), half_beta);
}
}
}
}
} else if (_N == 5) {
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
double ans_ij = energy(particles[i])
* energy(particles[j])
* pow(angleSquared(particles[i],particles[j]), half_beta);
for (unsigned int k = j + 1; k < particles.size(); k++) {
double ans_ijk = ans_ij
* energy(particles[k])
* pow(angleSquared(particles[i],particles[k]), half_beta)
* pow(angleSquared(particles[j],particles[k]), half_beta);
for (unsigned int l = k + 1; l < particles.size(); l++) {
double ans_ijkl = ans_ijk
* energy(particles[l])
* pow(angleSquared(particles[i],particles[l]), half_beta)
* pow(angleSquared(particles[j],particles[l]), half_beta)
* pow(angleSquared(particles[k],particles[l]), half_beta);
for (unsigned int m = l + 1; m < particles.size(); m++) {
answer += ans_ijkl
* energy(particles[m])
* pow(angleSquared(particles[i],particles[m]), half_beta)
* pow(angleSquared(particles[j],particles[m]), half_beta)
* pow(angleSquared(particles[k],particles[m]), half_beta)
* pow(angleSquared(particles[l],particles[m]), half_beta);
}
}
}
}
}
} else {
assert(_N <= 5);
}
} else {
assert(_strategy == slow || _strategy == storage_array);
}
return answer;
}
double EnergyCorrelator::energy(const PseudoJet& jet) const {
if (_measure == pt_R) {
return jet.perp();
} else if (_measure == E_theta || _measure == E_inv) {
return jet.e();
} else {
assert(_measure==pt_R || _measure==E_theta || _measure==E_inv);
return std::numeric_limits<double>::quiet_NaN();
}
}
double EnergyCorrelator::angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const {
if (_measure == pt_R) {
return jet1.squared_distance(jet2);
} else if (_measure == 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 = jet1.px()*jet1.px() + jet1.py()*jet1.py() + jet1.pz()*jet1.pz();
double norm2 = jet2.px()*jet2.px() + jet2.py()*jet2.py() + jet2.pz()*jet2.pz();
double costheta = dot/(sqrt(norm1 * norm2));
if (costheta > 1.0) costheta = 1.0; // Need to handle case of numerical overflow
double theta = acos(costheta);
return theta*theta;
} else if (_measure == E_inv) {
if (jet1.E() < 0.0000001 || jet2.E() < 0.0000001) return 0.0;
else {
double dot4 = max(jet1.E()*jet2.E() - jet1.px()*jet2.px() - jet1.py()*jet2.py() - jet1.pz()*jet2.pz(),0.0);
return 2.0 * dot4 / jet1.E() / jet2.E();
}
} else {
assert(_measure==pt_R || _measure==E_theta || _measure==E_inv);
return std::numeric_limits<double>::quiet_NaN();
}
}
+ double EnergyCorrelatorGeneralized::multiply_angles(double angle_list[], int n_angles, unsigned int N_total) const {
+ // Compute the product of the n_angles smallest angles.
+ // std::partial_sort could also work, but since angle_list contains
+ // less than 10 elements, this way is usually faster.
+ double product = 1;
+
+ for (int a = 0; a < n_angles; a++) {
+ double cur_min = angle_list[0];
+ int cur_min_pos = 0;
+ for (unsigned int b = 1; b < N_total; b++) {
+ if (angle_list[b] < cur_min) {
+ cur_min = angle_list[b];
+ cur_min_pos = b;
+ }
+ }
+
+ // multiply it by the next smallest
+ product *= cur_min;
+ angle_list[cur_min_pos] = INT_MAX;
+ }
+ return product;
+ }
+
+ void EnergyCorrelatorGeneralized::precompute_energies_and_angles(std::vector<fastjet::PseudoJet> const &particles, double* energyStore, double** angleStore) const {
+ // Fill storage with energy/angle information
+ unsigned int nC = particles.size();
+ for (unsigned int i = 0; i < nC; i++) {
+ angleStore[i] = new double[i];
+ }
+
+ double half_beta = _beta/2.0;
+ for (unsigned int i = 0; i < particles.size(); i++) {
+ energyStore[i] = energy(particles[i]);
+ for (unsigned int j = 0; j < i; j++) {
+ if (half_beta == 1.0){
+ angleStore[i][j] = angleSquared(particles[i], particles[j]);
+ } else {
+ angleStore[i][j] = pow(angleSquared(particles[i], particles[j]), half_beta);
+ }
+ }
+ }
+ }
+
+ double EnergyCorrelatorGeneralized::evaluate_n3(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const {
+ unsigned int N_total = 3;
+ double angle1, angle2, angle3;
+ double angle;
+ double answer = 0;
+
+ for (unsigned int i = 2; i < nC; i++) {
+ for (unsigned int j = 1; j < i; j++) {
+ double mult_energy_i_j = energyStore[i] * energyStore[j];
+
+ for (unsigned int k = 0; k < j; k++) {
+ angle1 = angleStore[i][j];
+ angle2 = angleStore[i][k];
+ angle3 = angleStore[j][k];
+
+ double angle_list[] = {angle1, angle2, angle3};
+
+ if (n_angles == N_total) {
+ angle = angle1 * angle2 * angle3;
+ } else {
+ angle = multiply_angles(angle_list, n_angles, N_total);
+ }
+
+ answer += mult_energy_i_j
+ * energyStore[k]
+ * angle;
+ }
+ }
+ }
+ return answer;
+ }
+
+ double EnergyCorrelatorGeneralized::evaluate_n4(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const {
+ double answer = 0;
+ double angle1, angle2, angle3, angle4, angle5, angle6;
+ unsigned int N_total = 6;
+ double angle;
+
+ for (unsigned int i = 3; i < nC; i++) {
+ for (unsigned int j = 2; j < i; j++) {
+ for (unsigned int k = 1; k < j; k++) {
+ for (unsigned int l = 0; l < k; l++) {
+
+ angle1 = angleStore[i][j];
+ angle2 = angleStore[i][k];
+ angle3 = angleStore[i][l];
+ angle4 = angleStore[j][k];
+ angle5 = angleStore[j][l];
+ angle6 = angleStore[k][l];
+
+ double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6};
+
+ if (n_angles == N_total) {
+ angle = angle1 * angle2 * angle3 * angle4 * angle5 * angle6;
+ } else {
+ angle = multiply_angles(angle_list, n_angles, N_total);
+ }
+
+ answer += energyStore[i]
+ * energyStore[j]
+ * energyStore[k]
+ * energyStore[l]
+ * angle;
+ }
+ }
+ }
+ }
+ return answer;
+ }
+
+ double EnergyCorrelatorGeneralized::evaluate_n5(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const {
+
+ double answer = 0;
+ double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10;
+ unsigned int N_total = 10;
+ double angle;
+
+ for (unsigned int i = 4; i < nC; i++) {
+ for (unsigned int j = 3; j < i; j++) {
+ for (unsigned int k = 2; k < j; k++) {
+ for (unsigned int l = 1; l < k; l++) {
+ for (unsigned int m = 0; m < l; m++) {
+
+ angle1 = angleStore[i][j];
+ angle2 = angleStore[i][k];
+ angle3 = angleStore[i][l];
+ angle4 = angleStore[i][m];
+ angle5 = angleStore[j][k];
+ angle6 = angleStore[j][l];
+ angle7 = angleStore[j][m];
+ angle8 = angleStore[k][l];
+ angle9 = angleStore[k][m];
+ angle10 = angleStore[l][m];
+
+ double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8,
+ angle9, angle10};
+
+ angle = multiply_angles(angle_list, n_angles, N_total);
+
+ answer += energyStore[i]
+ * energyStore[j]
+ * energyStore[k]
+ * energyStore[l]
+ * energyStore[m]
+ * angle;
+ }
+ }
+ }
+ }
+ }
+ return answer;
+ }
+
double EnergyCorrelatorGeneralized::result(const PseudoJet& jet) const {
// if jet does not have constituents, throw error
if (!jet.has_constituents()) throw Error("EnergyCorrelator called on jet with no constituents.");
// Throw an error if N < 0
// Not needed if N is unsigned integer
//if (_N < 0 ) throw Error("N cannot be negative");
-
// get N = 0 case out of the way
if (_N == 0) return 1.0;
// take care of N = 1 case.
if (_N == 1) return 1.0;
// find constituents
std::vector<fastjet::PseudoJet> particles = jet.constituents();
double answer = 0.0;
// return zero if the number of constituents is less than _N for the ECFG
if (particles.size() < _N) return 0.0 ;
// The normalization is the energy or pt of the jet, which is also ECF(1, beta)
double EJ = _helper_correlator.result(jet);
// The overall normalization
double norm = pow(EJ, _N);
// Find the max number of angles and throw an error if unsuitable
int N_total = int(_N*(_N-1)/2);
if (_angles > N_total) throw Error("Requested number of angles for EnergyCorrelatorGeneralized is larger than number of angles available");
if (_angles < -1) throw Error("Negative number of angles called for EnergyCorrelatorGeneralized");
double half_beta = _beta/2.0;
// take care of N = 2 case.
if (_N == 2) {
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) { //note offset by one so that angle is never called on identical pairs
answer += energy(particles[i])
* energy(particles[j])
* pow(angleSquared(particles[i],particles[j]), half_beta)/norm;
}
}
return answer;
}
// if N > 4, then throw error
if (_N > 5) {
throw Error("EnergyCorrelatorGeneralized is only hard coded for N = 0,1,2,3,4,5");
}
// Now deal with N = 3,4,5. Different options if storage array is used or not.
if (_strategy == EnergyCorrelator::storage_array) {
// For N > 2, fill static storage array to save computation time.
+ unsigned int nC = particles.size();
// Make energy storage
- std::vector<double> energyStore;
- energyStore.resize(particles.size());
+ double *energyStore = new double[nC];
// Make angular storage
- std::vector < std::vector<double> > angleStore;
- angleStore.resize(particles.size());
- for (unsigned int i = 0; i < angleStore.size(); i++) {
- angleStore[i].resize(i);
- }
+ double **angleStore = new double*[nC];
- // Fill storage with energy/angle information
- for (unsigned int i = 0; i < particles.size(); i++) {
- energyStore[i] = energy(particles[i]);
- for (unsigned int j = 0; j < i; j++) {
- if (half_beta == 1.0){
- angleStore[i][j] = angleSquared(particles[i], particles[j]);
- } else {
- angleStore[i][j] = pow(angleSquared(particles[i], particles[j]), half_beta);
- }
- }
+ precompute_energies_and_angles(particles, energyStore, angleStore);
+
+ unsigned int n_angles = _angles;
+ if (_angles < 0) {
+ n_angles = N_total;
}
// now do recursion
if (_N == 3) {
- unsigned int N_total = 3;
- double angle1, angle2, angle3;
- double angle;
-
- for (unsigned int i = 2; i < particles.size(); i++) {
- for (unsigned int j = 1; j < i; j++) {
- for (unsigned int k = 0; k < j; k++) {
-
- angle1 = angleStore[i][j];
- angle2 = angleStore[i][k];
- angle3 = angleStore[j][k];
-
- if (_angles == -1){
- // When _angles = -1, it defaults to considering all angles.
- angle = angle1*angle2*angle3;
- } else {
- double angle_list[] = {angle1, angle2, angle3};
- std::vector<double> angle_vector(angle_list, angle_list + N_total);
- std::partial_sort(angle_vector.begin(), angle_vector.begin() + _angles, angle_vector.begin() + N_total);
-
- angle = angle_vector[0];
- for ( int l = 1; l < _angles; l++) { angle *= angle_vector[l]; }
- }
-
- answer += energyStore[i]
- * energyStore[j]
- * energyStore[k]
- * angle /norm;
- }
- }
- }
+ answer = evaluate_n3(nC, n_angles, energyStore, angleStore) / norm;
} else if (_N == 4) {
- double angle1, angle2, angle3, angle4, angle5, angle6;
- unsigned int N_total = 6;
- double angle;
-
- for (unsigned int i = 3; i < particles.size(); i++) {
- for (unsigned int j = 2; j < i; j++) {
- for (unsigned int k = 1; k < j; k++) {
- for (unsigned int l = 0; l < k; l++) {
-
- angle1 = angleStore[i][j];
- angle2 = angleStore[i][k];
- angle3 = angleStore[i][l];
- angle4 = angleStore[j][k];
- angle5 = angleStore[j][l];
- angle6 = angleStore[k][l];
-
- if(_angles == -1) {
- angle = angle1*angle2*angle3*angle4*angle5*angle6;
- } else {
-
- double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6};
- std::vector<double> angle_vector(angle_list, angle_list + N_total);
- std::partial_sort(angle_vector.begin(), angle_vector.begin() + _angles, angle_vector.begin() + N_total);
-
- angle = angle_vector[0];
- for ( int s = 1; s < _angles; s++) { angle *= angle_vector[s]; }
-
- }
- answer += energyStore[i]
- * energyStore[j]
- * energyStore[k]
- * energyStore[l]
- * angle /norm;
-
- }
- }
- }
- }
+ answer = evaluate_n4(nC, n_angles, energyStore, angleStore) / norm;
} else if (_N == 5) {
- double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10;
- unsigned int N_total = 10;
- double angle;
-
- for (unsigned int i = 4; i < particles.size(); i++) {
- for (unsigned int j = 3; j < i; j++) {
- for (unsigned int k = 2; k < j; k++) {
- for (unsigned int l = 1; l < k; l++) {
- for (unsigned int m = 0; m < l; m++) {
-
- angle1 = angleStore[i][j];
- angle2 = angleStore[i][k];
- angle3 = angleStore[i][l];
- angle4 = angleStore[i][m];
- angle5 = angleStore[j][k];
- angle6 = angleStore[j][l];
- angle7 = angleStore[j][m];
- angle8 = angleStore[k][l];
- angle9 = angleStore[k][m];
- angle10 = angleStore[l][m];
-
- if (_angles == -1){
- angle = angle1*angle2*angle3*angle4*angle5*angle6*angle7*angle8*angle9*angle10;
- } else {
- double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6,
- angle7, angle8, angle9, angle10};
- std::vector<double> angle_vector(angle_list, angle_list + N_total);
- std::partial_sort(angle_vector.begin(), angle_vector.begin() + _angles, angle_vector.begin() + N_total);
-
- angle = angle_vector[0];
- for ( int s = 1; s < _angles; s++) { angle *= angle_vector[s]; }
- }
-
- answer += energyStore[i]
- * energyStore[j]
- * energyStore[k]
- * energyStore[l]
- * energyStore[m]
- * angle /norm;
- }
- }
- }
- }
- }
-
+ answer = evaluate_n5(nC, n_angles, energyStore, angleStore) / norm;
} else {
assert(_N <= 5);
}
} else if (_strategy == EnergyCorrelator::slow) {
if (_N == 3) {
unsigned int N_total = 3;
double angle1, angle2, angle3;
double angle;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
for (unsigned int k = j + 1; k < particles.size(); k++) {
angle1 = angleSquared(particles[i], particles[j]);
angle2 = angleSquared(particles[i], particles[k]);
angle3 = angleSquared(particles[j], particles[k]);
if (_angles == -1){
angle = angle1*angle2*angle3;
} else {
double angle_list[] = {angle1, angle2, angle3};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
angle = angle_vector[0];
for ( int l = 1; l < _angles; l++) { angle = angle * angle_vector[l]; }
}
answer += energy(particles[i])
* energy(particles[j])
* energy(particles[k])
* pow(angle, half_beta) /norm;
}
}
}
} else if (_N == 4) {
double angle1, angle2, angle3, angle4, angle5, angle6;
unsigned int N_total = 6;
double angle;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
for (unsigned int k = j + 1; k < particles.size(); k++) {
for (unsigned int l = k + 1; l < particles.size(); l++) {
angle1 = angleSquared(particles[i], particles[j]);
angle2 = angleSquared(particles[i], particles[k]);
angle3 = angleSquared(particles[i], particles[l]);
angle4 = angleSquared(particles[j], particles[k]);
angle5 = angleSquared(particles[j], particles[l]);
angle6 = angleSquared(particles[k], particles[l]);
if(_angles == -1) {
angle = angle1*angle2*angle3*angle4*angle5*angle6;
} else {
double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
angle = angle_vector[0];
for ( int s = 1; s < _angles; s++) { angle = angle * angle_vector[s]; }
}
answer += energy(particles[i])
* energy(particles[j])
* energy(particles[k])
* energy(particles[l])
* pow(angle, half_beta)/norm;
}
}
}
}
} else if (_N == 5) {
double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10;
unsigned int N_total = 10;
double angle;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
for (unsigned int k = j + 1; k < particles.size(); k++) {
for (unsigned int l = k + 1; l < particles.size(); l++) {
for (unsigned int m = l + 1; m < particles.size(); m++) {
angle1 = angleSquared(particles[i], particles[j]);
angle2 = angleSquared(particles[i], particles[k]);
angle3 = angleSquared(particles[i], particles[l]);
angle4 = angleSquared(particles[j], particles[k]);
angle5 = angleSquared(particles[j], particles[l]);
angle6 = angleSquared(particles[k], particles[l]);
angle7 = angleSquared(particles[m], particles[i]);
angle8 = angleSquared(particles[m], particles[j]);
angle9 = angleSquared(particles[m], particles[k]);
angle10 = angleSquared(particles[m], particles[l]);
if (_angles == -1){
angle = angle1*angle2*angle3*angle4*angle5*angle6*angle7*angle8*angle9*angle10;
} else {
double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6,
angle7, angle8, angle9, angle10};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
angle = angle_vector[0];
for ( int s = 1; s < _angles; s++) { angle = angle * angle_vector[s]; }
}
answer += energy(particles[i])
* energy(particles[j])
* energy(particles[k])
* energy(particles[l])
* energy(particles[m])
* pow(angle, half_beta) /norm;
}
}
}
}
}
} else {
assert(_N <= 5);
}
} else {
assert(_strategy == EnergyCorrelator::slow || _strategy == EnergyCorrelator::storage_array);
}
return answer;
}
std::vector<double> EnergyCorrelatorGeneralized::result_all_angles(const PseudoJet& jet) const {
// if jet does not have constituents, throw error
if (!jet.has_constituents()) throw Error("EnergyCorrelator called on jet with no constituents.");
// Throw an error if N < 1
if (_N < 1 ) throw Error("N cannot be negative or zero");
// get the N = 1 case out of the way
if (_N == 1) {
std::vector<double> ans (1, 1.0);
return ans;
}
// find constituents
std::vector<fastjet::PseudoJet> particles = jet.constituents();
// return zero if the number of constituents is less than _N for the ECFG
if (particles.size() < _N) {
std::vector<double> ans (_N, 0.0);
return ans;
}
// The normalization is the energy or pt of the jet, which is also ECF(1, beta)
double EJ = _helper_correlator.result(jet);
// The overall normalization
double norm = pow(EJ, _N);
// Find the max number of angles and throw an error if it unsuitable
int N_total = _N * (_N - 1)/2;
double half_beta = _beta/2.0;
// take care of N = 2 case.
if (_N == 2) {
double answer = 0.0;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) { //note offset by one so that angle is never called on identical pairs
answer += energy(particles[i])
* energy(particles[j])
* pow(angleSquared(particles[i],particles[j]), half_beta)/norm;
}
}
std::vector<double> ans(N_total, answer);
return ans;
}
// Prepare the answer vector
std::vector<double> ans (N_total, 0.0);
// if N > 4, then throw error
if (_N > 5) {
throw Error("EnergyCorrelatorGeneralized is only hard coded for N = 0,1,2,3,4,5");
}
// Now deal with N = 3,4,5. Different options if storage array is used or not.
if (_strategy == EnergyCorrelator::storage_array) {
// For N > 2, fill static storage array to save computation time.
// Make energy storage
std::vector<double> energyStore;
energyStore.resize(particles.size());
// Make angular storage
std::vector < std::vector<double> > angleStore;
angleStore.resize(particles.size());
for (unsigned int i = 0; i < angleStore.size(); i++) {
angleStore[i].resize(i);
}
// Fill storage with energy/angle information
for (unsigned int i = 0; i < particles.size(); i++) {
energyStore[i] = energy(particles[i]);
for (unsigned int j = 0; j < i; j++) {
if (half_beta == 1){
angleStore[i][j] = angleSquared(particles[i], particles[j]);
} else {
angleStore[i][j] = pow(angleSquared(particles[i], particles[j]), half_beta);
}
}
}
// now do recursion
if (_N == 3) {
double angle1, angle2, angle3;
for (unsigned int i = 2; i < particles.size(); i++) {
for (unsigned int j = 1; j < i; j++) {
for (unsigned int k = 0; k < j; k++) {
angle1 = angleStore[i][j];
angle2 = angleStore[i][k];
angle3 = angleStore[j][k];
double angle_list[] = {angle1, angle2, angle3};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
std::vector<double> final_angles (N_total, angle_vector[0]);
double z_product = energyStore[i] * energyStore[j] * energyStore[k]/norm;
ans[0] += z_product * final_angles[0];
for ( int s=1 ; s<N_total; s++){
final_angles[s] = final_angles[s-1]*angle_vector[s];
ans[s] += z_product * final_angles[s];
}
}
}
}
} else if (_N == 4) {
double angle1, angle2, angle3, angle4, angle5, angle6;
for (unsigned int i = 3; i < particles.size(); i++) {
for (unsigned int j = 2; j < i; j++) {
for (unsigned int k = 1; k < j; k++) {
for (unsigned int l = 0; l < k; l++) {
angle1 = angleStore[i][j];
angle2 = angleStore[i][k];
angle3 = angleStore[i][l];
angle4 = angleStore[j][k];
angle5 = angleStore[j][l];
angle6 = angleStore[k][l];
double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
std::vector<double> final_angles (N_total, angle_vector[0]);
double z_product = energyStore[i] * energyStore[j] * energyStore[k] * energyStore [l]/norm;
ans[0] += z_product * final_angles[0];
for ( int s=1 ; s<N_total; s++){
final_angles[s] = final_angles[s-1]*angle_vector[s];
ans[s] += z_product * final_angles[s];
}
}
}
}
}
} else if (_N == 5) {
double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10;
for (unsigned int i = 4; i < particles.size(); i++) {
for (unsigned int j = 3; j < i; j++) {
for (unsigned int k = 2; k < j; k++) {
for (unsigned int l = 1; l < k; l++) {
for (unsigned int m = 0; m < l; m++) {
angle1 = angleStore[i][j];
angle2 = angleStore[i][k];
angle3 = angleStore[i][l];
angle4 = angleStore[i][m];
angle5 = angleStore[j][k];
angle6 = angleStore[j][l];
angle7 = angleStore[j][m];
angle8 = angleStore[k][l];
angle9 = angleStore[k][m];
angle10 = angleStore[l][m];
double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6,
- angle7, angle8, angle9, angle10};
+ angle7, angle8, angle9, angle10};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
std::vector<double> final_angles (N_total, angle_vector[0]);
double z_product = energyStore[i] * energyStore[j] * energyStore[k]
* energyStore[l] * energyStore[m]/norm;
ans[0] += z_product * final_angles[0];
for ( int s=1 ; s<N_total; s++){
final_angles[s] = final_angles[s-1]*angle_vector[s];
ans[s] += z_product * final_angles[s];
}
}
}
}
}
}
} else {
assert(_N <= 5);
}
} else if (_strategy == EnergyCorrelator::slow) {
if (_N == 3) {
double angle1, angle2, angle3;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
for (unsigned int k = j + 1; k < particles.size(); k++) {
angle1 = pow(angleSquared(particles[i], particles[j]), half_beta);
angle2 = pow(angleSquared(particles[i], particles[k]), half_beta);
angle3 = pow(angleSquared(particles[j], particles[k]), half_beta);
double angle_list[] = {angle1, angle2, angle3};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
std::vector<double> final_angles (N_total, angle_vector[0]);
double z_product = energy(particles[i])
* energy(particles[j])
* energy(particles[k])/norm;
ans[0] += z_product * final_angles[0];
for ( int s=1 ; s<N_total; s++){
final_angles[s] = final_angles[s-1]*angle_vector[s];
ans[s] += z_product * final_angles[s];
}
}
}
}
} else if (_N == 4) {
double angle1, angle2, angle3, angle4, angle5, angle6;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
for (unsigned int k = j + 1; k < particles.size(); k++) {
for (unsigned int l = k + 1; l < particles.size(); l++) {
angle1 = pow(angleSquared(particles[i], particles[j]), half_beta);
angle2 = pow(angleSquared(particles[i], particles[k]), half_beta);
angle3 = pow(angleSquared(particles[i], particles[l]), half_beta);
angle4 = pow(angleSquared(particles[j], particles[k]), half_beta);
angle5 = pow(angleSquared(particles[j], particles[l]), half_beta);
angle6 = pow(angleSquared(particles[k], particles[l]), half_beta);
double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
std::vector<double> final_angles (N_total, angle_vector[0]);
double z_product = energy(particles[i])
* energy(particles[j])
* energy(particles[k])
* energy(particles[l])/norm;
ans[0] += z_product * final_angles[0];
for ( int s=1 ; s<N_total; s++){
final_angles[s] = final_angles[s-1]*angle_vector[s];
ans[s] += z_product * final_angles[s];
}
}
}
}
}
} else if (_N == 5) {
double angle1, angle2, angle3, angle4, angle5, angle6, angle7, angle8, angle9, angle10;
for (unsigned int i = 0; i < particles.size(); i++) {
for (unsigned int j = i + 1; j < particles.size(); j++) {
for (unsigned int k = j + 1; k < particles.size(); k++) {
for (unsigned int l = k + 1; l < particles.size(); l++) {
for (unsigned int m = l + 1; m < particles.size(); m++) {
angle1 = pow(angleSquared(particles[i], particles[j]), half_beta);
angle2 = pow(angleSquared(particles[i], particles[k]), half_beta);
angle3 = pow(angleSquared(particles[i], particles[l]), half_beta);
angle4 = pow(angleSquared(particles[j], particles[k]), half_beta);
angle5 = pow(angleSquared(particles[j], particles[l]), half_beta);
angle6 = pow(angleSquared(particles[k], particles[l]), half_beta);
angle7 = pow(angleSquared(particles[m], particles[i]), half_beta);
angle8 = pow(angleSquared(particles[m], particles[j]), half_beta);
angle9 = pow(angleSquared(particles[m], particles[k]), half_beta);
angle10 = pow(angleSquared(particles[m], particles[l]), half_beta);
double angle_list[] = {angle1, angle2, angle3, angle4, angle5, angle6,
angle7, angle8, angle9, angle10};
std::vector<double> angle_vector(angle_list, angle_list + N_total);
std::sort(angle_vector.begin(), angle_vector.begin() + N_total);
std::vector<double> final_angles (N_total, angle_vector[0]);
double z_product = energy(particles[i])
* energy(particles[j])
* energy(particles[k])
* energy(particles[l])
* energy(particles[m])/norm;
ans[0] += z_product * final_angles[0];
for ( int s=1 ; s<N_total; s++){
final_angles[s] = final_angles[s-1]*angle_vector[s];
ans[s] += z_product * final_angles[s];
}
}
}
}
}
}
} else {
assert(_N <= 5);
}
} else {
assert(_strategy == EnergyCorrelator::slow || _strategy == EnergyCorrelator::storage_array);
}
return ans;
}
// call _helper_correlator to get energy information
double EnergyCorrelatorGeneralized::energy(const PseudoJet& jet) const {
return _helper_correlator.energy(jet);
}
-
+
// call _helper_correlator to get angle information
double EnergyCorrelatorGeneralized::angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const {
return _helper_correlator.angleSquared(jet1, jet2);
}
string EnergyCorrelator::description_parameters() const {
ostringstream oss;
oss << "N=" << _N << ", beta=" << _beta;
if (_measure == pt_R) oss << ", pt_R measure";
else if (_measure == E_theta) oss << ", E_theta measure";
else if (_measure == E_inv) oss << ", E_inv measure";
else throw Error("unrecognized measure");
if (_strategy == slow) oss << " and 'slow' strategy";
else if (_strategy == storage_array) oss << " and 'storage_array' strategy";
else throw Error("unrecognized strategy");
return oss.str();
}
string EnergyCorrelator::description_no_N() const {
ostringstream oss;
oss << "beta=" << _beta;
if (_measure == pt_R) oss << ", pt_R measure";
else if (_measure == E_theta) oss << ", E_theta measure";
else if (_measure == E_inv) oss << ", E_inv measure";
else throw Error("unrecognized measure");
if (_strategy == slow) oss << " and 'slow' strategy";
else if (_strategy == storage_array) oss << " and 'storage_array' strategy";
else throw Error("unrecognized strategy");
return oss.str();
}
string EnergyCorrelator::description() const {
ostringstream oss;
oss << "Energy Correlator ECF(N,beta) for ";
oss << description_parameters();
return oss.str();
}
string EnergyCorrelatorRatio::description() const {
ostringstream oss;
oss << "Energy Correlator ratio ECF(N+1,beta)/ECF(N,beta) for ";
oss << EnergyCorrelator(_N,_beta,_measure,_strategy).description_parameters();
return oss.str();
}
string EnergyCorrelatorDoubleRatio::description() const {
ostringstream oss;
oss << "Energy Correlator double ratio ECF(N-1,beta)ECF(N+1,beta)/ECF(N,beta)^2 for ";
oss << EnergyCorrelator(_N,_beta,_measure,_strategy).description_parameters();
return oss.str();
}
string EnergyCorrelatorC1::description() const {
ostringstream oss;
oss << "Energy Correlator observable C1 ECF(2,beta)/ECF(1,beta)^2 for ";
oss << EnergyCorrelator(2,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorC2::description() const {
ostringstream oss;
oss << "Energy Correlator observable C2 ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorD2::description() const {
ostringstream oss;
oss << "Energy Correlator observable D2 ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorGeneralizedD2::description() const {
ostringstream oss;
oss << "Energy Correlator observable D2 ECFN(3,alpha)/ECFN(2,beta)^(3 alpha/beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorNseries::description() const {
ostringstream oss;
oss << "Energy Correlator observable N_n ECFG(2,n+1,beta)/ECFG(1,n,beta)^2 for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorMseries::description() const {
ostringstream oss;
oss << "Energy Correlator observable M_n ECFG(1,n+1,beta)/ECFG(1,n,beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorN2::description() const {
ostringstream oss;
oss << "Energy Correlator observable N2 ECFG(2,3,beta)/ECFG(1,2,beta)^2 for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorN3::description() const {
ostringstream oss;
oss << "Energy Correlator observable N3 ECFG(2,4,beta)/ECFG(1,3,beta)^2 for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorM2::description() const {
ostringstream oss;
oss << "Energy Correlator observable M2 ECFG(1,3,beta)/ECFG(1,2,beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorCseries::description() const {
ostringstream oss;
oss << "Energy Correlator double ratio ECFN(N-1,beta)ECFN(N+1,beta)/ECFN(N,beta)^2 for ";
oss << EnergyCorrelator(_N,_beta,_measure,_strategy).description_parameters();
return oss.str();
}
string EnergyCorrelatorUseries::description() const {
ostringstream oss;
oss << "Energy Correlator observable U_n ECFG(1,n+1,beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorU1::description() const {
ostringstream oss;
oss << "Energy Correlator observable U_1 ECFG(1,2,beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorU2::description() const {
ostringstream oss;
oss << "Energy Correlator observable U_2 ECFG(1,3,beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
string EnergyCorrelatorU3::description() const {
ostringstream oss;
oss << "Energy Correlator observable U_3 ECFG(1,4,beta) for ";
oss << EnergyCorrelator(3,_beta,_measure,_strategy).description_no_N();
return oss.str();
}
} // namespace contrib
FASTJET_END_NAMESPACE
Index: contrib/contribs/EnergyCorrelator/trunk/example.cc
===================================================================
--- contrib/contribs/EnergyCorrelator/trunk/example.cc (revision 1096)
+++ contrib/contribs/EnergyCorrelator/trunk/example.cc (revision 1097)
@@ -1,655 +1,686 @@
// Example showing usage of energy correlator classes.
//
// Compile it with "make example" and run it with
//
// ./example < ../data/single-event.dat
//
// Copyright (c) 2013-2016
// Andrew Larkoski, Lina Necib, Gavin Salam, and Jesse Thaler
//
// $Id$
//----------------------------------------------------------------------
// 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 <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#include <iomanip>
#include <stdlib.h>
#include <stdio.h>
#include <ctime>
#include <iostream>
#include <istream>
#include <fstream>
#include <sstream>
#include <string>
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include <sstream>
#include "EnergyCorrelator.hh" // In external code, this should be fastjet/contrib/EnergyCorrelator.hh
using namespace std;
using namespace fastjet;
using namespace fastjet::contrib;
// forward declaration to make things clearer
void read_event(vector<PseudoJet> &event);
void analyze(const vector<PseudoJet> & input_particles);
//----------------------------------------------------------------------
int main(){
//----------------------------------------------------------
// read in input particles
vector<PseudoJet> event;
read_event(event);
cout << "# read an event with " << event.size() << " particles" << endl;
//----------------------------------------------------------
// illustrate how this EnergyCorrelator contrib works
analyze(event);
return 0;
}
// read in input particles
void read_event(vector<PseudoJet> &event){
string line;
while (getline(cin, line)) {
istringstream linestream(line);
// take substrings to avoid problems when there are extra "pollution"
// characters (e.g. line-feed).
if (line.substr(0,4) == "#END") {return;}
if (line.substr(0,1) == "#") {continue;}
double px,py,pz,E;
linestream >> px >> py >> pz >> E;
PseudoJet particle(px,py,pz,E);
// push event onto back of full_event vector
event.push_back(particle);
}
}
////////
//
// Main Routine for Analysis
//
///////
void analyze(const vector<PseudoJet> & input_particles) {
/////// EnergyCorrelator /////////////////////////////
// Initial clustering with anti-kt algorithm
JetAlgorithm algorithm = antikt_algorithm;
double jet_rad = 1.00; // jet radius for anti-kt algorithm
JetDefinition jetDef = JetDefinition(algorithm,jet_rad,E_scheme,Best);
ClusterSequence clust_seq(input_particles,jetDef);
vector<PseudoJet> antikt_jets = sorted_by_pt(clust_seq.inclusive_jets());
for (int j = 0; j < 2; j++) { // Two hardest jets per event
if (antikt_jets[j].perp() > 200) {
PseudoJet myJet = antikt_jets[j];
// various values of beta
vector<double> betalist;
betalist.push_back(0.1);
betalist.push_back(0.2);
betalist.push_back(0.5);
betalist.push_back(1.0);
betalist.push_back(1.5);
betalist.push_back(2.0);
// various values of alpha
vector<double> alphalist;
alphalist.push_back(0.1);
alphalist.push_back(0.2);
alphalist.push_back(0.5);
alphalist.push_back(1.0);
// checking the two energy/angle modes
vector<EnergyCorrelator::Measure> measurelist;
measurelist.push_back(EnergyCorrelator::pt_R);
measurelist.push_back(EnergyCorrelator::E_theta);
//measurelist.push_back(EnergyCorrelator::E_inv);
vector<string> modename;
modename.push_back("pt_R");
modename.push_back("E_theta");
//modename.push_back("E_inv");
for (unsigned int M = 0; M < measurelist.size(); M++) {
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelator: ECF(N,beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s %14s %14s %14s %15s\n","beta", "N=1 (GeV)", "N=2 (GeV^2)", "N=3 (GeV^3)", "N=4 (GeV^4)", "N=5 (GeV^5)");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelator ECF0(0,beta,measurelist[M]);
EnergyCorrelator ECF1(1,beta,measurelist[M]);
EnergyCorrelator ECF2(2,beta,measurelist[M]);
EnergyCorrelator ECF3(3,beta,measurelist[M]);
EnergyCorrelator ECF4(4,beta,measurelist[M]);
EnergyCorrelator ECF5(5,beta,measurelist[M]);
printf("%7.3f %14.2f %14.2f %14.2f %14.2f %15.2f \n",beta,ECF1(myJet),ECF2(myJet),ECF3(myJet),ECF4(myJet),ECF5(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorRatio: r_N^(beta) = ECF(N+1,beta)/ECF(N,beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s %14s %14s %14s %15s \n","beta", "N=0 (GeV)", "N=1 (GeV)", "N=2 (GeV)", "N=3 (GeV)","N=4 (GeV)");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorRatio r0(0,beta,measurelist[M]);
EnergyCorrelatorRatio r1(1,beta,measurelist[M]);
EnergyCorrelatorRatio r2(2,beta,measurelist[M]);
EnergyCorrelatorRatio r3(3,beta,measurelist[M]);
EnergyCorrelatorRatio r4(4,beta,measurelist[M]);
printf("%7.3f %14.4f %14.4f %14.4f %14.4f %15.4f \n",beta,r0(myJet),r1(myJet),r2(myJet),r3(myJet),r4(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorDoubleRatio: C_N^(beta) = r_N^(beta)/r_{N-1}^(beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s %14s %14s %14s \n","beta", "N=1", "N=2", "N=3", "N=4");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorDoubleRatio C1(1,beta,measurelist[M]);
EnergyCorrelatorDoubleRatio C2(2,beta,measurelist[M]);
EnergyCorrelatorDoubleRatio C3(3,beta,measurelist[M]);
EnergyCorrelatorDoubleRatio C4(4,beta,measurelist[M]);
printf("%7.3f %14.6f %14.6f %14.6f %14.6f \n",beta,C1(myJet),C2(myJet),C3(myJet),C4(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorC1: C_1^(beta) = ECF(2,beta)/ECF(1,beta)^2 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta","C1 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorC1 c1(beta,measurelist[M]);
printf("%7.3f %14.6f \n",beta,c1(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorC2: C_2^(beta) = ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta","C2 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorC2 c2(beta,measurelist[M]);
printf("%7.3f %14.6f \n",beta,c2(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorD2: D_2^(beta) = ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta","D2 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorD2 d2(beta,measurelist[M]);
printf("%7.3f %14.6f \n",beta,d2(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorGeneralizedD2: D_2^(alpha, beta) = ECFN(3,alpha)/ECFN(2,beta)^(3*alpha/beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %18s %18s %18s %18s\n","beta","alpha = 0.100","alpha = 0.200","alpha = 0.500","alpha = 1.000");
for (unsigned int B = 1; B < betalist.size(); B++) {
double beta = betalist[B];
printf("%7.3f ", beta);
for (unsigned int A = 0; A < alphalist.size(); A++) {
double alpha = alphalist[A];
EnergyCorrelatorGeneralizedD2 d2(alpha, beta, measurelist[M]);
printf("%18.4f ", d2(myJet));
}
printf("\n");
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorGeneralized (angles = N Choose 2): ECFN(N, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %7s %14s %14s %14s\n","beta", "N=1", "N=2", "N=3", "N=4");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorGeneralized ECF1(-1,1, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF2(-1,2, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF3(-1,3, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF4(-1,4, beta, measurelist[M]);
//EnergyCorrelatorGeneralized ECF5(-1, 5, beta, measurelist[M]);
printf("%7.3f %7.2f %14.10f %14.10f %14.10f \n", beta, ECF1(myJet), ECF2(myJet), ECF3(myJet),
ECF4(myJet));
}
cout << "-------------------------------------------------------------------------------------" <<
endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorGeneralized: ECFG(angles, N, beta=1) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %7s %14s %14s %14s\n","angles", "N=1", "N=2", "N=3", "N=4");
double beta = 1.0;
for (unsigned int A = 1; A < 2; A++) {
double angle = A;
EnergyCorrelatorGeneralized ECF1(angle, 1, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF2(angle, 2, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF3(angle, 3, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF4(angle, 4, beta, measurelist[M], EnergyCorrelator::slow);
printf("%7.0f %7.2f %14.10f %14.10f %14.10f \n", angle, ECF1(myJet), ECF2(myJet), ECF3(myJet),
ECF4(myJet));
}
for (unsigned int A = 2; A < 4; A++) {
double angle = A;
EnergyCorrelatorGeneralized ECF3(angle, 3, beta, measurelist[M]);
EnergyCorrelatorGeneralized ECF4(angle, 4, beta, measurelist[M]);
printf("%7.0f %7s %14s %14.10f %14.10f \n", angle, " " , " " ,ECF3(myJet), ECF4(myJet));
}
for (unsigned int A = 4; A < 7; A++) {
double angle = A;
EnergyCorrelatorGeneralized ECF4(angle, 4, beta, measurelist[M]);
printf("%7.0f %7s %14s %14s %14.10f \n", angle, " ", " ", " ", ECF4(myJet) );
}
cout << "-------------------------------------------------------------------------------------" <<
endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorNseries: N_i(beta) = ECFG(i+1, 2, beta)/ECFG(i, 1, beta)^2 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s %14s %14s \n","beta", "N=1", "N=2", "N=3");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorNseries N1(1,beta,measurelist[M]);
EnergyCorrelatorNseries N2(2,beta,measurelist[M]);
EnergyCorrelatorNseries N3(3,beta,measurelist[M]);
-
printf("%7.3f %14.6f %14.6f %14.6f \n",beta,N1(myJet),N2(myJet),N3(myJet));
+
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorN2: N2(beta) = ECFG(3, 2, beta)/ECFG(2, 1, beta)^2 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta", "N2 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorN2 N2(beta,measurelist[M]);
printf("%7.3f %14.6f \n",beta,N2(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorN3: N3(beta) = ECFG(4, 2, beta)/ECFG(3, 1, beta)^2 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta", "N3 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorN3 N3(beta,measurelist[M]);
printf("%7.3f %14.6f \n",beta,N3(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorMseries: M_i(beta) = ECFG(i+1, 1, beta)/ECFN(i, 1, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s %14s %14s \n","beta", "N=1", "N=2", "N=3");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorMseries M1(1,beta,measurelist[M]);
EnergyCorrelatorMseries M2(2,beta,measurelist[M]);
EnergyCorrelatorMseries M3(3,beta,measurelist[M]);
printf("%7.3f %14.6f %14.6f %14.6f \n",beta,M1(myJet),M2(myJet),M3(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorM2: M2(beta) = ECFG(3, 1, beta)/ECFG(3, 1, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta", "M2 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorM2 M2(beta,measurelist[M]);
printf("%7.3f %14.6f \n",beta,M2(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorCseries: C_i(beta) = ECFN(i-1, beta)*ECFN(i+1, beta)/ECFN(i, beta)^2 with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %20s %20s %20s \n","beta", "N=1", "N=2", "N=3");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorCseries C1(1,beta,measurelist[M]);
EnergyCorrelatorCseries C2(2,beta,measurelist[M]);
EnergyCorrelatorCseries C3(3,beta,measurelist[M]);
printf("%7.3f %20.10f %20.10f %20.10f \n",beta,C1(myJet),C2(myJet),C3(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorUseries: U_i(beta) = ECFG(i+1, 1, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %20s %20s %20s \n","beta", "N=1", "N=2", "N=3");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorUseries U1(1,beta,measurelist[M]);
EnergyCorrelatorUseries U2(2,beta,measurelist[M]);
EnergyCorrelatorUseries U3(3,beta,measurelist[M]);
printf("%7.3f %20.10f %20.10f %20.10f \n",beta,U1(myJet),U2(myJet),U3(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorU1: U1(beta) = ECFG(2, 1, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta", "U1 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorU1 U1(beta,measurelist[M]);
printf("%7.3f %14.10f \n",beta,U1(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorU2: U2(beta) = ECFG(3, 1, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta", "U2 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorU2 U2(beta,measurelist[M]);
printf("%7.3f %14.10f \n",beta,U2(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
cout << "EnergyCorrelatorU3: U3(beta) = ECFG(4, 1, beta) with " << modename[M] << endl;
cout << "-------------------------------------------------------------------------------------" << endl;
printf("%7s %14s \n","beta", "U3 obs");
for (unsigned int B = 0; B < betalist.size(); B++) {
double beta = betalist[B];
EnergyCorrelatorU3 U3(beta,measurelist[M]);
printf("%7.3f %14.10f \n",beta,U3(myJet));
}
cout << "-------------------------------------------------------------------------------------" << endl << endl;
// timing tests for the developers
double do_timing_test = false;
if (do_timing_test) {
cout << "jet with pt = " << myJet.pt() << " and " << myJet.constituents().size() << " constituents" << endl;
clock_t clock_begin, clock_end;
double num_iter;
double beta = 0.5;
cout << setprecision(6);
// test C1
num_iter = 20000;
clock_begin = clock();
EnergyCorrelatorDoubleRatio C1s(1,beta,measurelist[M],EnergyCorrelator::slow);
EnergyCorrelatorDoubleRatio C1f(1,beta,measurelist[M],EnergyCorrelator::storage_array);
cout << "timing " << C1s.description() << endl;
cout << "timing " << C1f.description() << endl;
for (int t = 0; t < num_iter; t++) {
C1s(myJet);
}
clock_end = clock();
cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C1"<< endl;
num_iter = 20000;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
C1f(myJet);
}
clock_end = clock();
cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C1"<< endl;
// test C2
num_iter = 1000;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorDoubleRatio C2(2,beta,measurelist[M],EnergyCorrelator::slow);
C2(myJet);
}
clock_end = clock();
cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C2"<< endl;
num_iter = 10000;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorDoubleRatio C2(2,beta,measurelist[M],EnergyCorrelator::storage_array);
C2(myJet);
}
clock_end = clock();
cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C2"<< endl;
// test C3
num_iter = 100;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorDoubleRatio C3(3,beta,measurelist[M],EnergyCorrelator::slow);
C3(myJet);
}
clock_end = clock();
cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C3"<< endl;
num_iter = 3000;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorDoubleRatio C3(3,beta,measurelist[M],EnergyCorrelator::storage_array);
C3(myJet);
}
clock_end = clock();
cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C3"<< endl;
// test C4
num_iter = 10;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorDoubleRatio C4(4,beta,measurelist[M],EnergyCorrelator::slow);
C4(myJet);
}
clock_end = clock();
cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C4"<< endl;
num_iter = 300;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorDoubleRatio C4(4,beta,measurelist[M],EnergyCorrelator::storage_array);
C4(myJet);
}
clock_end = clock();
cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per C4"<< endl;
// test N2
num_iter = 10;
clock_begin = clock();
+ num_iter = 300;
+ clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
- EnergyCorrelatorN2 N2(beta,measurelist[M],EnergyCorrelator::slow);
+ EnergyCorrelatorN2 N2(beta,measurelist[M],EnergyCorrelator::storage_array);
N2(myJet);
}
clock_end = clock();
- cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N2"<< endl;
+ EnergyCorrelatorN2 N2test(beta,measurelist[M],EnergyCorrelator::storage_array);
+ cout << "Beta is: "<< beta << endl;
+ cout << "Result of N2: "<< N2test(myJet) << endl;
+ cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N2"<< endl;
+
num_iter = 300;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
- EnergyCorrelatorN2 N2(beta,measurelist[M],EnergyCorrelator::storage_array);
- N2(myJet);
+ EnergyCorrelatorN3 N3(beta,measurelist[M],EnergyCorrelator::storage_array);
+ N3(myJet);
}
clock_end = clock();
- cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N2"<< endl;
+ EnergyCorrelatorN3 N3test(beta,measurelist[M],EnergyCorrelator::storage_array);
+ cout << "Beta is: "<< beta << endl;
+ cout << "Result of N3: "<< N3test(myJet) << endl;
+ cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N3"<< endl;
+
- // test N3
- num_iter = 10;
- clock_begin = clock();
+
+ num_iter = 300;
+ clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
- EnergyCorrelatorN3 N3(beta,measurelist[M],EnergyCorrelator::slow);
- N3(myJet);
+ EnergyCorrelatorGeneralized ECF1(2,3, beta, measurelist[M]);
+ ECF1(myJet);
}
clock_end = clock();
- cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N3"<< endl;
+ EnergyCorrelatorGeneralized ECF1test(2,3, beta, measurelist[M]);
+ cout << "Beta is: "<< beta << endl;
+ cout << "Result of 2e3: "<< ECF1test(myJet) << endl;
+ cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per 2e3"<< endl;
+
num_iter = 300;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
- EnergyCorrelatorN3 N3(beta,measurelist[M],EnergyCorrelator::storage_array);
- N3(myJet);
+ EnergyCorrelatorGeneralized ECF3(2,4, beta, measurelist[M]);
+ ECF3(myJet);
}
clock_end = clock();
- cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per N3"<< endl;
+ EnergyCorrelatorGeneralized ECF2test(2,4, beta, measurelist[M]);
+ cout << "Beta is: "<< beta << endl;
+ cout << "Result of 2e4: "<< ECF2test(myJet) << endl;
+ cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per 2e4"<< endl;
+
+
+// num_iter = 300;
+// clock_begin = clock();
+// for (int t = 0; t < num_iter; t++) {
+// EnergyCorrelatorGeneralized ECF5(2,5, beta, measurelist[M]);
+// ECF5(myJet);
+// }
+// clock_end = clock();
+// EnergyCorrelatorGeneralized ECF5test(2,5, beta, measurelist[M]);
+// cout << "Beta is: "<< beta << endl;
+// cout << "Result of 2e5: "<< ECF5test(myJet) << endl;
+// cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per 2e5"<< endl;
+//
// test M2
num_iter = 10;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorM2 M2(beta,measurelist[M],EnergyCorrelator::slow);
M2(myJet);
}
clock_end = clock();
cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M2"<< endl;
num_iter = 300;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorM2 M2(beta,measurelist[M],EnergyCorrelator::storage_array);
M2(myJet);
}
clock_end = clock();
cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M2"<< endl;
// test M3
num_iter = 10;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorMseries M3(3,beta,measurelist[M],EnergyCorrelator::slow);
M3(myJet);
}
clock_end = clock();
cout << "Slow method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M3"<< endl;
num_iter = 300;
clock_begin = clock();
for (int t = 0; t < num_iter; t++) {
EnergyCorrelatorMseries M3(3,beta,measurelist[M],EnergyCorrelator::storage_array);
M3(myJet);
}
clock_end = clock();
cout << "Storage array method: " << (clock_end-clock_begin)/double(CLOCKS_PER_SEC*num_iter)*1000 << " ms per M3"<< endl;
}
}
}
}
}
Index: contrib/contribs/EnergyCorrelator/trunk/EnergyCorrelator.hh
===================================================================
--- contrib/contribs/EnergyCorrelator/trunk/EnergyCorrelator.hh (revision 1096)
+++ contrib/contribs/EnergyCorrelator/trunk/EnergyCorrelator.hh (revision 1097)
@@ -1,1070 +1,1074 @@
#ifndef __FASTJET_CONTRIB_ENERGYCORRELATOR_HH__
#define __FASTJET_CONTRIB_ENERGYCORRELATOR_HH__
// EnergyCorrelator Package
// Questions/Comments? Email the authors.
// larkoski@mit.edu, lnecib@mit.edu,
// gavin.salam@cern.ch jthaler@jthaler.net
//
// Copyright (c) 2013-2016
// Andrew Larkoski, Lina Necib Gavin Salam, and Jesse Thaler
//
// $Id$
//----------------------------------------------------------------------
// 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 <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#include <fastjet/internal/base.hh>
#include "fastjet/FunctionOfPseudoJet.hh"
#include <string>
#include <climits>
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib{
/// \mainpage EnergyCorrelator contrib
///
/// The EnergyCorrelator contrib provides an implementation of energy
/// correlators and their ratios as described in arXiv:1305.0007 by
/// Larkoski, Salam and Thaler. Additionally, the ratio observable
/// D2 described in arXiv:1409.6298 by Larkoski, Moult and Neill
/// is also included in this contrib. Finally, a generalized version of
/// the energy correlation functions is added, defined in
/// arXiv:1609.07483 by Moult, Necib and Thaler, which allow the
/// definition of the M series, N series, and U series observables.
/// There is also a generalized version of D2.
///
///
/// <p>There are 4 main classes:
///
/// - EnergyCorrelator
/// - EnergyCorrelatorRatio
/// - EnergyCorrelatorDoubleRatio
/// - EnergyCorrelatorGeneralized
///
/// <p>There are five classes that define useful combinations of the ECFs.
///
/// - EnergyCorrelatorNseries
/// - EnergyCorrelatorMseries
/// - EnergyCorrelatorUseries
/// - EnergyCorrelatorD2
/// - EnergyCorrelatorGeneralizedD2
///
/// <p> There are also aliases for easier access:
/// - EnergyCorrelatorCseries (same as EnergyCorrelatorDoubleRatio)
/// - EnergyCorrelatorC1 (EnergyCorrelatorCseries with i=1)
/// - EnergyCorrelatorC2 (EnergyCorrelatorCseries with i=2)
/// - EnergyCorrelatorN2 (EnergyCorrelatorNseries with i=2)
/// - EnergyCorrelatorN3 (EnergyCorrelatorNseries with i=3)
/// - EnergyCorrelatorM2 (EnergyCorrelatorMseries with i=2)
/// - EnergyCorrelatorU1 (EnergyCorrelatorUseries with i=1)
/// - EnergyCorrelatorU2 (EnergyCorrelatorUseries with i=2)
/// - EnergyCorrelatorU3 (EnergyCorrelatorUseries with i=3)
///
/// Each of these classes is a FastJet FunctionOfPseudoJet.
/// EnergyCorrelatorDoubleRatio (which is equivalent to EnergyCorrelatorCseries)
/// is in particular is useful for quark/gluon discrimination and boosted
/// object tagging.
///
/// Using the original 2- and 3-point correlators, EnergyCorrelationD2 has
/// been shown to be the optimal combination for boosted 2-prong tagging.
///
/// The EnergyCorrelatorNseries and EnergyCorrelatorMseries use
/// generalized correlation functions with different angular scaling,
/// and are intended for use on 2-prong and 3-prong jets.
/// The EnergyCorrelatorUseries is useful for quark/gluon discrimimation.
///
/// See the file example.cc for an illustration of usage and
/// example_basic_usage.cc for the most commonly used functions.
//------------------------------------------------------------------------
/// \class EnergyCorrelator
/// ECF(N,beta) is the N-point energy correlation function, with an angular exponent beta.
///
/// It is defined as follows
///
/// - \f$ \mathrm{ECF}(1,\beta) = \sum_i E_i \f$
/// - \f$ \mathrm{ECF}(2,\beta) = \sum_{i<j} E_i E_j \theta_{ij}^\beta \f$
/// - \f$ \mathrm{ECF}(3,\beta) = \sum_{i<j<k} E_i E_j E_k (\theta_{ij} \theta_{ik} \theta_{jk})^\beta \f$
/// - \f$ \mathrm{ECF}(4,\beta) = \sum_{i<j<k<l} E_i E_j E_k E_l (\theta_{ij} \theta_{ik} \theta_{il} \theta_{jk} \theta_{jl} \theta_{kl})^\beta \f$
/// - ...
///
/// The correlation can be determined with energies and angles (as
/// given above) or with transverse momenta and boost invariant angles
/// (the code's default). The choice is controlled by
/// EnergyCorrelator::Measure provided in the constructor.
///
/// The current implementation handles values of N up to and including 5.
/// Run times scale as n^N/N!, where n is the number of particles in a jet.
class EnergyCorrelator : public FunctionOfPseudoJet<double> {
friend class EnergyCorrelatorGeneralized; ///< This allow ECFG to access the energy and angle definitions
///< of this class, which are otherwise private.
public:
enum Measure {
pt_R, ///< use transverse momenta and boost-invariant angles,
///< eg \f$\mathrm{ECF}(2,\beta) = \sum_{i<j} p_{ti} p_{tj} \Delta R_{ij}^{\beta} \f$
E_theta, /// use energies and angles,
/// eg \f$\mathrm{ECF}(2,\beta) = \sum_{i<j} E_{i} E_{j} \theta_{ij}^{\beta} \f$
E_inv /// use energies and invariant mass,
/// eg \f$\mathrm{ECF}(2,\beta) = \sum_{i<j} E_{i} E_{j} (\frac{2 p_{i} \cdot p_{j}}{E_{i} E_{j}})^{\beta/2} \f$
};
enum Strategy {
slow, ///< interparticle angles are not cached.
///< For N>=3 this leads to many expensive recomputations,
///< but has only O(n) memory usage for n particles
storage_array /// the interparticle angles are cached. This gives a significant speed
/// improvement for N>=3, but has a memory requirement of (4n^2) bytes.
};
public:
/// constructs an N-point correlator with angular exponent beta,
/// using the specified choice of energy and angular measure as well
/// one of two possible underlying computational Strategy
EnergyCorrelator(unsigned int N,
double beta,
Measure measure = pt_R,
Strategy strategy = storage_array) :
_N(N), _beta(beta), _measure(measure), _strategy(strategy) {};
/// destructor
virtual ~EnergyCorrelator(){}
/// returns the value of the energy correlator for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
/// returns the the part of the description related to the parameters
std::string description_parameters() const;
std::string description_no_N() const;
private:
unsigned int _N;
double _beta;
Measure _measure;
Strategy _strategy;
double energy(const PseudoJet& jet) const;
double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const;
};
// core EnergyCorrelator::result code in .cc file.
//------------------------------------------------------------------------
/// \class EnergyCorrelatorRatio
/// A class to calculate the ratio of (N+1)-point to N-point energy correlators,
/// ECF(N+1,beta)/ECF(N,beta),
/// called \f$ r_N^{(\beta)} \f$ in the publication.
class EnergyCorrelatorRatio : public FunctionOfPseudoJet<double> {
public:
/// constructs an (N+1)-point to N-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorRatio(unsigned int N,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _N(N), _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorRatio() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
unsigned int _N;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorRatio::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelator(_N + 1, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelator(_N, _beta, _measure, _strategy).result(jet);
return numerator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorDoubleRatio
/// Calculates the double ratio of energy correlators, ECF(N-1,beta)*ECF(N+1)/ECF(N,beta)^2.
///
/// A class to calculate a double ratio of energy correlators,
/// ECF(N-1,beta)*ECF(N+1,beta)/ECF(N,beta)^2,
/// called \f$C_N^{(\beta)}\f$ in the publication, and equal to
/// \f$ r_N^{(\beta)}/r_{N-1}^{(\beta)} \f$.
///
class EnergyCorrelatorDoubleRatio : public FunctionOfPseudoJet<double> {
public:
EnergyCorrelatorDoubleRatio(unsigned int N,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _N(N), _beta(beta), _measure(measure), _strategy(strategy) {
if (_N < 1) throw Error("EnergyCorrelatorDoubleRatio: N must be 1 or greater.");
};
virtual ~EnergyCorrelatorDoubleRatio() {}
/// returns the value of the energy correlator double-ratio for a
/// jet's constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
unsigned int _N;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorDoubleRatio::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelator(_N - 1, _beta, _measure, _strategy).result(jet) * EnergyCorrelator(_N + 1, _beta, _measure, _strategy).result(jet);
double denominator = pow(EnergyCorrelator(_N, _beta, _measure, _strategy).result(jet), 2.0);
return numerator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorC1
/// A class to calculate the normalized 2-point energy correlators,
/// ECF(2,beta)/ECF(1,beta)^2,
/// called \f$ C_1^{(\beta)} \f$ in the publication.
class EnergyCorrelatorC1 : public FunctionOfPseudoJet<double> {
public:
/// constructs a 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorC1(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorC1() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorC1::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelator(2, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelator(1, _beta, _measure, _strategy).result(jet);
return numerator/denominator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorC2
/// A class to calculate the double ratio of 3-point to 2-point
/// energy correlators,
/// ECF(3,beta)*ECF(1,beta)/ECF(2,beta)^2,
/// called \f$ C_2^{(\beta)} \f$ in the publication.
class EnergyCorrelatorC2 : public FunctionOfPseudoJet<double> {
public:
/// constructs a 3-point to 2-point correlator double ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorC2(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorC2() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorC2::result(const PseudoJet& jet) const {
double numerator3 = EnergyCorrelator(3, _beta, _measure, _strategy).result(jet);
double numerator1 = EnergyCorrelator(1, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelator(2, _beta, _measure, _strategy).result(jet);
return numerator3*numerator1/denominator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorD2
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// ECF(3,beta)*ECF(1,beta)^3/ECF(2,beta)^3,
/// called \f$ D_2^{(\beta)} \f$ in the publication.
class EnergyCorrelatorD2 : public FunctionOfPseudoJet<double> {
public:
/// constructs an 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorD2(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorD2() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorD2::result(const PseudoJet& jet) const {
double numerator3 = EnergyCorrelator(3, _beta, _measure, _strategy).result(jet);
double numerator1 = EnergyCorrelator(1, _beta, _measure, _strategy).result(jet);
double denominator2 = EnergyCorrelator(2, _beta, _measure, _strategy).result(jet);
return numerator3*numerator1*numerator1*numerator1/denominator2/denominator2/denominator2;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorGeneralized
/// A generalized and normalized version of the N-point energy correlators, with
/// angular exponent beta and v number of pairwise angles. When \f$v = {N \choose 2}\f$
/// (or, for convenience, \f$v = -1\f$), EnergyCorrelatorGeneralized just gives normalized
/// versions of EnergyCorrelator:
/// - \f$ \mathrm{ECFG}(-1,1,\beta) = \mathrm{ECFN}(N,\beta) = \mathrm{ECF}(N,\beta)/\mathrm{ECF}(1,\beta)\f$
///
/// Note that there is no separate class that implements ECFN, though it is a
/// notation that we will use in this documentation. Examples of the low-point normalized
/// correlators are:
/// - \f$\mathrm{ECFN}(1,\beta) = 1\f$
/// - \f$\mathrm{ECFN}(2,\beta) = \sum_{i<j} z_i z_j \theta_{ij}^\beta \f$
/// - \f$\mathrm{ECFN}(3,\beta) = \sum_{i<j<k} z_i z_j z_k (\theta_{ij} \theta_{ik} \theta_{jk})^\beta \f$
/// - \f$\mathrm{ECFN}(4,\beta) = \sum_{i<j<k<l} z_i z_j z_k z_l (\theta_{ij} \theta_{ik} \theta_{il} \theta_{jk} \theta_{jl} \theta_{kl})^\beta \f$
/// - ...
/// where the \f$z_i\f$'s are the energy fractions.
///
/// When a new value of v is given, the generalized energy correlators are defined as
/// - \f$\mathrm{ECFG}(0,1,\beta) = 1\f$
/// - \f$\mathrm{ECFG}(1,2,\beta) = \sum_{i<j} z_i z_j \theta_{ij}^\beta \f$
/// - \f$\mathrm{ECFG}(v,3,\beta) = \sum_{i<j<k} z_i z_j z_k \prod_{m = 1}^{v} \min^{(m)} \{ \theta_{ij}, \theta_{jk}, \theta_{ki} \}^\beta \f$
/// - \f$\mathrm{ECFG}(v,4,\beta) = \sum_{i<j<k<l} z_i z_j z_k z_l \prod_{m = 1}^{v} \min^{(m)} \{ \theta_{ij}, \theta_{ik}, \theta_{il}, \theta_{jk}, \theta_{jl}, \theta_{kl} \}^\beta \f$
/// - \f$\mathrm{ECFG}(v,n,\beta) = \sum_{i_1 < i_2 < \dots < i_n} z_{i_1} z_{i_2} \dots z_{i_n} \prod_{m = 1}^{v} \min^{(m)}_{s < t \in \{i_1, i_2 , \dots, i_n \}} \{ \theta_{st}^{\beta} \}\f$,
///
/// where \f$\min^{(m)}\f$ means the m-th smallest element of the list.
///
/// The correlation can be determined with energies and angles (as
/// given above) or with transverse momenta and boost invariant angles
/// (the code's default). The choice is controlled by
/// EnergyCorrelator::Measure provided in the constructor.
///
/// The current implementation handles values of N up to and including 5.
///
class EnergyCorrelatorGeneralized : public FunctionOfPseudoJet<double> {
public:
/// constructs an N-point correlator with v_angles pairwise angles
/// and angular exponent beta,
/// using the specified choice of energy and angular measure as well
/// one of two possible underlying computational Strategy
EnergyCorrelatorGeneralized(int v_angles,
unsigned int N,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _angles(v_angles), _N(N), _beta(beta), _measure(measure), _strategy(strategy), _helper_correlator(1,_beta, _measure, _strategy) {};
/// destructor
virtual ~EnergyCorrelatorGeneralized(){}
/// returns the value of the normalized energy correlator for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::vector<double> result_all_angles(const PseudoJet& jet) const;
private:
int _angles;
unsigned int _N;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
EnergyCorrelator _helper_correlator;
double energy(const PseudoJet& jet) const;
double angleSquared(const PseudoJet& jet1, const PseudoJet& jet2) const;
-
+ double multiply_angles(double angles[], int n_angles, unsigned int N_total) const;
+ void precompute_energies_and_angles(std::vector<fastjet::PseudoJet> const &particles, double* energyStore, double** angleStore) const;
+ double evaluate_n3(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const;
+ double evaluate_n4(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const;
+ double evaluate_n5(unsigned int nC, unsigned int n_angles, double* energyStore, double** angleStore) const;
};
//------------------------------------------------------------------------
/// \class EnergyCorrelatorGeneralizedD2
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// ECFN(3,alpha)/ECFN(2,beta)^3 alpha/beta,
/// called \f$ D_2^{(\alpha, \beta)} \f$ in the publication.
class EnergyCorrelatorGeneralizedD2 : public FunctionOfPseudoJet<double> {
public:
/// constructs an 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorGeneralizedD2(
double alpha,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _alpha(alpha), _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorGeneralizedD2() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _alpha;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorGeneralizedD2::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelatorGeneralized(-1, 3, _alpha, _measure, _strategy).result(jet);
double denominator = EnergyCorrelatorGeneralized(-1, 2, _beta, _measure, _strategy).result(jet);
return numerator/pow(denominator, 3.0*_alpha/_beta);
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorNseries
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// N_n = ECFG(2,n+1,beta)/ECFG(1,n,beta)^2,
/// called \f$ N_i^{(\alpha, \beta)} \f$ in the publication.
/// By definition, N_1^{beta} = ECFG(1, 2, 2*beta), where the angular exponent
/// is twice as big since the N series should involve two pairwise angles.
class EnergyCorrelatorNseries : public FunctionOfPseudoJet<double> {
public:
/// constructs a n 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorNseries(
unsigned int n,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _n(n), _beta(beta), _measure(measure), _strategy(strategy) {
if (_n < 1) throw Error("EnergyCorrelatorNseries: n must be 1 or greater.");
};
virtual ~EnergyCorrelatorNseries() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
unsigned int _n;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorNseries::result(const PseudoJet& jet) const {
if (_n == 1) return EnergyCorrelatorGeneralized(1, 2, 2*_beta, _measure, _strategy).result(jet);
// By definition, N1 = ECFN(2, 2 beta)
double numerator = EnergyCorrelatorGeneralized(2, _n + 1, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelatorGeneralized(1, _n, _beta, _measure, _strategy).result(jet);
return numerator/denominator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorN2
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// ECFG(2,3,beta)/ECFG(1,2,beta)^2,
/// called \f$ N_2^{(\beta)} \f$ in the publication.
class EnergyCorrelatorN2 : public FunctionOfPseudoJet<double> {
public:
/// constructs an 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorN2(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorN2() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorN2::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelatorGeneralized(2, 3, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet);
return numerator/denominator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorN3
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// ECFG(2,4,beta)/ECFG(1,3,beta)^2,
/// called \f$ N_3^{(\beta)} \f$ in the publication.
class EnergyCorrelatorN3 : public FunctionOfPseudoJet<double> {
public:
/// constructs an 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorN3(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorN3() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorN3::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelatorGeneralized(2, 4, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelatorGeneralized(1, 3, _beta, _measure, _strategy).result(jet);
return numerator/denominator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorMseries
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// M_n = ECFG(1,n+1,beta)/ECFG(1,n,beta),
/// called \f$ M_i^{(\alpha, \beta)} \f$ in the publication.
/// By definition, M_1^{beta} = ECFG(1,2,beta)
class EnergyCorrelatorMseries : public FunctionOfPseudoJet<double> {
public:
/// constructs a n 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorMseries(
unsigned int n,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _n(n), _beta(beta), _measure(measure), _strategy(strategy) {
if (_n < 1) throw Error("EnergyCorrelatorMseries: n must be 1 or greater.");
};
virtual ~EnergyCorrelatorMseries() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
unsigned int _n;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorMseries::result(const PseudoJet& jet) const {
if (_n == 1) return EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet);
double numerator = EnergyCorrelatorGeneralized(1, _n + 1, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelatorGeneralized(1, _n, _beta, _measure, _strategy).result(jet);
return numerator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorM2
/// A class to calculate the observable formed from the ratio of the
/// 3-point and 2-point energy correlators,
/// ECFG(1,3,beta)/ECFG(1,2,beta),
/// called \f$ M_2^{(\beta)} \f$ in the publication.
class EnergyCorrelatorM2 : public FunctionOfPseudoJet<double> {
public:
/// constructs an 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorM2(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorM2() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorM2::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelatorGeneralized(1, 3, _beta, _measure, _strategy).result(jet);
double denominator = EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet);
return numerator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorCseries
/// Calculates the C series energy correlators, ECFN(N-1,beta)*ECFN(N+1,beta)/ECFN(N,beta)^2.
/// This is equivalent to EnergyCorrelatorDoubleRatio
///
/// A class to calculate a double ratio of energy correlators,
/// ECFN(N-1,beta)*ECFN(N+1,beta)/ECFN(N,beta)^2,
/// called \f$C_N^{(\beta)}\f$ in the publication, and equal to
/// \f$ r_N^{(\beta)}/r_{N-1}^{(\beta)} \f$.
///
class EnergyCorrelatorCseries : public FunctionOfPseudoJet<double> {
public:
EnergyCorrelatorCseries(unsigned int N,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _N(N), _beta(beta), _measure(measure), _strategy(strategy) {
if (_N < 1) throw Error("EnergyCorrelatorCseries: N must be 1 or greater.");
};
virtual ~EnergyCorrelatorCseries() {}
/// returns the value of the energy correlator double-ratio for a
/// jet's constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
unsigned int _N;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorCseries::result(const PseudoJet& jet) const {
double numerator = EnergyCorrelatorGeneralized(-1, _N - 1, _beta, _measure, _strategy).result(jet) * EnergyCorrelatorGeneralized(-1, _N + 1, _beta, _measure, _strategy).result(jet);
double denominator = pow(EnergyCorrelatorGeneralized(-1, _N, _beta, _measure, _strategy).result(jet), 2.0);
return numerator/denominator;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorUseries
/// A class to calculate the observable used for quark versus gluon discrimination
/// U_n = ECFG(1,n+1,beta),
/// called \f$ U_i^{(\beta)} \f$ in the publication.
class EnergyCorrelatorUseries : public FunctionOfPseudoJet<double> {
public:
/// constructs a n 3-point to 2-point correlator ratio with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorUseries(
unsigned int n,
double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _n(n), _beta(beta), _measure(measure), _strategy(strategy) {
if (_n < 1) throw Error("EnergyCorrelatorUseries: n must be 1 or greater.");
};
virtual ~EnergyCorrelatorUseries() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
unsigned int _n;
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorUseries::result(const PseudoJet& jet) const {
double answer = EnergyCorrelatorGeneralized(1, _n + 1, _beta, _measure, _strategy).result(jet);
return answer;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorU1
/// A class to calculate the observable formed from
/// ECFG(1,2,beta),
/// called \f$ U_1^{(\beta)} \f$ in the publication.
class EnergyCorrelatorU1 : public FunctionOfPseudoJet<double> {
public:
/// constructs a 2-point correlator with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorU1(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorU1() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorU1::result(const PseudoJet& jet) const {
double answer = EnergyCorrelatorGeneralized(1, 2, _beta, _measure, _strategy).result(jet);
return answer;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorU2
/// A class to calculate the observable formed from
/// ECFG(1,3,beta),
/// called \f$ U_2^{(\beta)} \f$ in the publication.
class EnergyCorrelatorU2 : public FunctionOfPseudoJet<double> {
public:
/// constructs a 3-point correlator with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorU2(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorU2() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorU2::result(const PseudoJet& jet) const {
double answer = EnergyCorrelatorGeneralized(1, 3, _beta, _measure, _strategy).result(jet);
return answer;
}
//------------------------------------------------------------------------
/// \class EnergyCorrelatorU3
/// A class to calculate the observable formed from
/// ECFG(1,4,beta),
/// called \f$ U_3^{(\beta)} \f$ in the publication.
class EnergyCorrelatorU3 : public FunctionOfPseudoJet<double> {
public:
/// constructs a 4-point correlator with
/// angular exponent beta, using the specified choice of energy and
/// angular measure as well one of two possible underlying
/// computational strategies
EnergyCorrelatorU3(double beta,
EnergyCorrelator::Measure measure = EnergyCorrelator::pt_R,
EnergyCorrelator::Strategy strategy = EnergyCorrelator::storage_array)
: _beta(beta), _measure(measure), _strategy(strategy) {};
virtual ~EnergyCorrelatorU3() {}
/// returns the value of the energy correlator ratio for a jet's
/// constituents. (Normally accessed by the parent class's
/// operator()).
double result(const PseudoJet& jet) const;
std::string description() const;
private:
double _beta;
EnergyCorrelator::Measure _measure;
EnergyCorrelator::Strategy _strategy;
};
inline double EnergyCorrelatorU3::result(const PseudoJet& jet) const {
double answer = EnergyCorrelatorGeneralized(1, 4, _beta, _measure, _strategy).result(jet);
return answer;
}
} // namespace contrib
FASTJET_END_NAMESPACE
#endif // __FASTJET_CONTRIB_ENERGYCORRELATOR_HH__
Index: contrib/contribs/EnergyCorrelator/trunk/VERSION
===================================================================
--- contrib/contribs/EnergyCorrelator/trunk/VERSION (revision 1096)
+++ contrib/contribs/EnergyCorrelator/trunk/VERSION (revision 1097)
@@ -1 +1 @@
-1.2.1-devel
\ No newline at end of file
+1.2.2-devel
\ No newline at end of file
Index: contrib/contribs/EnergyCorrelator/trunk/ChangeLog
===================================================================
--- contrib/contribs/EnergyCorrelator/trunk/ChangeLog (revision 1096)
+++ contrib/contribs/EnergyCorrelator/trunk/ChangeLog (revision 1097)
@@ -1,249 +1,262 @@
+2018-01-04 Lina Necib <lnecib@caltech.edu>
+
+ * EnergyCorrelator.hh/cc
+ Sped up evaluations of ECFG by a factor of 4.
+ Refactored some of the code in functions
+
+ * VERSION
+ Making 1.2.2-devel
+
+ * example.cc
+ Added some timings tests for ECFGs, N2, and N3.
+
+
2017-01-25 Jesse Thaler <jthaler@jthaler.net>
* EnergyCorrelator.hh/cc
Converting all _N to unsigned integers, removing _N < 0 warning
Added warning to constructors for N < 1 in some cases.
* VERSION
Making 1.2.1-devel
* Makefile
Added -Wextra to compiler flags
2016-10-07 Jesse Thaler <jthaler@jthaler.net>
* AUTHORS/COPYING:
Updated the publication arXiv number to 1609.07483.
* EnergyCorrelator.hh
Fixed typo in EnergyCorrelatorGeneralized description
* NEWS
Added "USeries" to news.
* VERSION
Changed version to 1.2.0
2016-09-27 Lina Necib <lnecib@mit.edu>
* EnergyCorrelator.hh/README
Updated the publication arXiv number to 1609.07483.
2016-09-23 Lina Necib <lnecib@mit.edu>
* EnergyCorrelator.cc
Made the evaluation of ECFG faster by replacing sort -> partial_sort and simplified multiplication syntax
* example.cc/ref
Fixed minor typo
* VERSION
Changed version to 1.2.0-rc3
2016-09-23 Lina Necib <lnecib@mit.edu>
* EnergyCorrelator.hh
Fixed minor doxygen typo
* example.cc/ref
Changed EnergyCorrelatorNormalized -> EnergyCorrelatorGeneralized in function explanations
* VERSION
Changed version to 1.2.0-rc2
2016-09-22 Jesse Thaler <jthaler@jthaler.net>
* EnergyCorrelator.cc/hh
Renamed EnergyCorrelatorNormalized -> EnergyCorrelatorGeneralized
Changed order of arguments for EnergyCorrelatorGeneralized
Updated doxygen documentation
* example_basic_usage.cc
Changed EnergyCorrelatorNormalized -> EnergyCorrelatorGeneralized
* README
Updated for EnergyCorrelatorGeneralized
2016-09-15 Lina Necib <lnecib@mit.edu>
*VERSION:
1.2.0-rc1
2016-08-24 Jesse Thaler <jthaler@jthaler.net>
Minor comment fixes throughout.
* EnergyCorrelator.cc/hh
Put in _helper_correlator when defining overall normalization
Removed angle error for result_all_angles()
*NEWS:
Made this more readable.
*README:
Clarified some of the documentation.
2016-08-23 Lina Necib <lnecib@mit.edu>
*VERSION:
*NEWS:
*AUTHORS:
*COPYING:
*README:
*EnergyCorrelator.cc
Added if statement that the ECF and ECFNs return exactly zero
if the number of particles is less than N of the ECF.
*EnergyCorrelator.hh
*example.cc
*example_basic_usage.cc
Updated this example.
2016-08-21 Lina Necib <lnecib@mit.edu>
*VERSION:
*NEWS:
*AUTHORS:
*COPYING:
*README:
*EnergyCorrelator.cc
Added Cseries.
*EnergyCorrelator.hh
Added Cseries.
*example.cc
Added example of Cseries
*example_basic_usage.cc
Simplified examples.
2016-08-17 Lina Necib <lnecib@mit.edu>
*VERSION:
*NEWS:
*AUTHORS:
Added placeholder for new reference
*COPYING:
Added placeholder for new reference
*README:
added information on different measures E_inv
*EnergyCorrelator.cc
Minor text edits + added comments
*EnergyCorrelator.hh
Minor text edits + added comments
*example_basic_usage.cc
added a simplified example file that shows the use of the
different observables.
2016-06-23 Lina Necib <lnecib@mit.edu>
*VERSION:
1.2.0-alpha0
*NEWS:
*AUTHORS:
Lina Necib
*COPYING:
*README:
added descriptions of functions that will appear shortly
arXiv:160X.XXXXX
*EnergyCorrelator.cc
added code to calculate normalized ECFS.
*EnergyCorrelator.hh
added code to calculate normalized ECFS, Nseries, generalized
D2, N2, N3, and M2.
*example.cc
added calculation of normalized ECFS, Nseries, generalized
D2, N2, N3, and M2 to example file.
2014-11-20 Jesse Thaler <jthaler@mit.edu>
* README:
Typos fixed
2014-11-13 Andrew Larkoski <larkoski@mit.edu>
*VERSION:
*NEWS:
release of version 1.1.0
*AUTHORS:
*COPYING:
*README:
added reference to arXiv:1409.6298.
*EnergyCorrelator.cc
*EnergyCorrelator.hh
added code to calculate C1, C2, and D2 observables.
*example.cc
added calculation of C1, C2, and D2 to example file.
2013-05-01 Gavin Salam <gavin.salam@cern.ch>
* VERSION:
* NEWS:
release of version 1.0.1
* README:
updated to reflect v1.0 interface.
2013-05-01 Jesse Thaler
* EnergyCorrelator.cc
Switched from NAN to std::numeric_limits<double>::quiet_NaN()
2013-04-30 Jesse Thaler
* EnergyCorrelator.cc
Implemented Gregory Soyez's suggestions on errors/asserts
2013-04-30 Gavin Salam <gavin.salam@cern.ch>
* VERSION:
* NEWS:
released v. 1.0.0
* EnergyCorrelator.hh:
* example.cc
small changes to documentation to reflect the change below + an
gave an explicit command-line in the example.
2013-04-29 Jesse Thaler
* EnergyCorrelator.cc
- Added support for N = 5
* example.cc|.ref
- Added N = 5 to test suite.
2013-04-29 Gavin Salam <gavin.salam@cern.ch>
* EnergyCorrelator.hh|cc:
- reduced memory usage from roughly 8n^2 to 4n^2 bytes (note that
sums are now performed in a different order, so results may
change to within rounding error).
- Implemented description() for all three classes.
- Worked on doxygen-style comments and moved some bits of code into
the .cc file.
* Doxyfile: *** ADDED ***
* example.cc:
developers' timing component now uses clock(), to get
finer-grained timing, and also outputs a description for some of
the correlators being used.
2013-04-26 + 04-27 Jesse Thaler <jthaler@MIT.EDU>
* EnergyCorrelator.hh|cc:
added temporary storage array for interparticle angles, to speed
up the energy correlator calculation for N>2
* example.cc
has internal option for printing out timing info.
2013-04-26 Gavin Salam <gavin.salam@cern.ch> + Jesse & Andrew
* Makefile:
added explicit dependencies
* example.cc (analyze):
added a little bit of commented code for timing tests.
2013-03-10 <jthaler>
Initial creation

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:47 PM (1 d, 16 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805450
Default Alt Text
(142 KB)

Event Timeline