Index: contrib/contribs/LundPlane/trunk/EEHelpers.hh
===================================================================
--- contrib/contribs/LundPlane/trunk/EEHelpers.hh (revision 0)
+++ contrib/contribs/LundPlane/trunk/EEHelpers.hh (revision 1286)
@@ -0,0 +1,262 @@
+// $Id: EEHelpers.hh 1153 2018-08-22 11:56:35Z akarlberg, lscyboz $
+//
+// Copyright (c) 2018-, Frederic A. Dreyer, A. Karlberg, Gavin P. Salam,
+// Ludovic Scyboz, Gregory Soyez
+//
+// 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 .
+//----------------------------------------------------------------------
+
+#ifndef __FASTJET_CONTRIB_EEHELPERS_HH__
+#define __FASTJET_CONTRIB_EEHELPERS_HH__
+
+#include "fastjet/PseudoJet.hh"
+#include
+
+FASTJET_BEGIN_NAMESPACE
+
+namespace contrib{
+
+//----------------------------------------------------------------------
+/// Returns the 3-vector cross-product of p1 and p2. If lightlike is false
+/// then the energy component is zero; if it's true the the energy component
+/// is arranged so that the vector is lighlike
+inline PseudoJet cross_product(const PseudoJet & p1, const PseudoJet & p2, bool lightlike=false) {
+ double px = p1.py() * p2.pz() - p2.py() * p1.pz();
+ double py = p1.pz() * p2.px() - p2.pz() * p1.px();
+ double pz = p1.px() * p2.py() - p2.px() * p1.py();
+
+ double E;
+ if (lightlike) {
+ E = sqrt(px*px + py*py + pz*pz);
+ } else {
+ E = 0.0;
+ }
+ return PseudoJet(px, py, pz, E);
+}
+
+/// Map angles to [-pi, pi]
+inline const double map_to_pi(const double &phi) {
+ if (phi < -M_PI) return phi + 2 * M_PI;
+ else if (phi > M_PI) return phi - 2 * M_PI;
+ else return phi;
+}
+
+inline double dot_product_3d(const PseudoJet & a, const PseudoJet & b) {
+ return a.px()*b.px() + a.py()*b.py() + a.pz()*b.pz();
+}
+
+/// Returns (1-cos theta) where theta is the angle between p1 and p2
+inline double one_minus_costheta(const PseudoJet & p1, const PseudoJet & p2) {
+
+ if (p1.m2() == 0 && p2.m2() == 0) {
+ // use the 4-vector dot product.
+ // For massless particles it gives us E1*E2*(1-cos theta)
+ return dot_product(p1,p2) / (p1.E() * p2.E());
+ } else {
+ double p1mod = p1.modp();
+ double p2mod = p2.modp();
+ double p1p2mod = p1mod*p2mod;
+ double dot = dot_product_3d(p1,p2);
+
+ if (dot > (1-std::numeric_limits::epsilon()) * p1p2mod) {
+ PseudoJet cross_result = cross_product(p1, p2, false);
+ // the mass^2 of cross_result is equal to
+ // -(px^2 + py^2 + pz^2) = (p1mod*p2mod*sintheta_ab)^2
+ // so we can get
+ return -cross_result.m2()/(p1p2mod * (p1p2mod+dot));
+ }
+
+ return 1.0 - dot/p1p2mod;
+
+ }
+}
+
+// Get the angle between two planes defined by normalized vectors
+// n1, n2. The sign is decided by the direction of a vector n.
+inline double signed_angle_between_planes(const PseudoJet& n1,
+ const PseudoJet& n2, const PseudoJet& n) {
+
+ // Two vectors passed as arguments should be normalised to 1.
+ assert(fabs(n1.modp()-1) < sqrt(std::numeric_limits::epsilon()) && fabs(n2.modp()-1) < sqrt(std::numeric_limits::epsilon()));
+
+ double omcost = one_minus_costheta(n1,n2);
+ double theta;
+
+ // If theta ~ pi, we return pi.
+ if(fabs(omcost-2) < sqrt(std::numeric_limits::epsilon())) {
+ theta = M_PI;
+ } else if (omcost > sqrt(std::numeric_limits::epsilon())) {
+ double cos_theta = 1.0 - omcost;
+ theta = acos(cos_theta);
+ } else {
+ // we are at small angles, so use small-angle formulas
+ theta = sqrt(2. * omcost);
+ }
+
+ PseudoJet cp = cross_product(n1,n2);
+ double sign = dot_product_3d(cp,n);
+
+ if (sign > 0) return theta;
+ else return -theta;
+}
+
+class Matrix3 {
+public:
+ /// constructs an empty matrix
+ Matrix3() : matrix_({{{{0,0,0}}, {{0,0,0}}, {{0,0,0}}}}) {}
+
+ /// constructs a diagonal matrix with "unit" along each diagonal entry
+ Matrix3(double unit) : matrix_({{{{unit,0,0}}, {{0,unit,0}}, {{0,0,unit}}}}) {}
+
+ /// constructs a matrix from the array,3> object
+ Matrix3(const std::array,3> & mat) : matrix_(mat) {}
+
+ /// returns the entry at row i, column j
+ inline double operator()(int i, int j) const {
+ return matrix_[i][j];
+ }
+
+ /// returns a matrix for carrying out azimuthal rotation
+ /// around the z direction by an angle phi
+ static Matrix3 azimuthal_rotation(double phi) {
+ double cos_phi = cos(phi);
+ double sin_phi = sin(phi);
+ Matrix3 phi_rot( {{ {{cos_phi, sin_phi, 0}},
+ {{-sin_phi, cos_phi, 0}},
+ {{0,0,1}}}});
+ return phi_rot;
+ }
+
+ /// returns a matrix for carrying out a polar-angle
+ /// rotation, in the z-x plane, by an angle theta
+ static Matrix3 polar_rotation(double theta) {
+ double cos_theta = cos(theta);
+ double sin_theta = sin(theta);
+ Matrix3 theta_rot( {{ {{cos_theta, 0, sin_theta}},
+ {{0,1,0}},
+ {{-sin_theta, 0, cos_theta}}}});
+ return theta_rot;
+ }
+
+ /// This provides a rotation matrix that takes the z axis to the
+ /// direction of p. With skip_pre_rotation = false (the default), it
+ /// has the characteristic that if p is close to the z axis then the
+ /// azimuthal angle of anything at much larger angle is conserved.
+ ///
+ /// If skip_pre_rotation is true, then the azimuthal angles are not
+ /// when p is close to the z axis.
+ template
+ static Matrix3 from_direction(const T & p, bool skip_pre_rotation = false) {
+ double pt = p.pt();
+ double modp = p.modp();
+ double cos_theta = p.pz() / modp;
+ double sin_theta = pt / modp;
+ double cos_phi, sin_phi;
+ if (pt > 0.0) {
+ cos_phi = p.px()/pt;
+ sin_phi = p.py()/pt;
+ } else {
+ cos_phi = 1.0;
+ sin_phi = 0.0;
+ }
+
+ Matrix3 phi_rot({{ {{ cos_phi,-sin_phi, 0 }},
+ {{ sin_phi, cos_phi, 0 }},
+ {{ 0, 0, 1 }} }});
+ Matrix3 theta_rot( {{ {{ cos_theta, 0, sin_theta }},
+ {{ 0, 1, 0 }},
+ {{-sin_theta, 0, cos_theta }} }});
+
+ // since we have orthogonal matrices, the inverse and transpose
+ // are identical; we use the transpose for the frontmost rotation
+ // because
+ if (skip_pre_rotation) {
+ return phi_rot * theta_rot;
+ } else {
+ return phi_rot * (theta_rot * phi_rot.transpose());
+ }
+ }
+
+ template
+ static Matrix3 from_direction_no_pre_rotn(const T & p) {
+ return from_direction(p,true);
+ }
+
+
+
+ /// returns the transposed matrix
+ Matrix3 transpose() const {
+ // 00 01 02
+ // 10 11 12
+ // 20 21 22
+ Matrix3 result = *this;
+ std::swap(result.matrix_[0][1],result.matrix_[1][0]);
+ std::swap(result.matrix_[0][2],result.matrix_[2][0]);
+ std::swap(result.matrix_[1][2],result.matrix_[2][1]);
+ return result;
+ }
+
+ // returns the product with another matrix
+ Matrix3 operator*(const Matrix3 & other) const {
+ Matrix3 result;
+ // r_{ij} = sum_k this_{ik} & other_{kj}
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ for (int k = 0; k < 3; k++) {
+ result.matrix_[i][j] += this->matrix_[i][k] * other.matrix_[k][j];
+ }
+ }
+ }
+ return result;
+ }
+
+ friend std::ostream & operator<<(std::ostream & ostr, const Matrix3 & mat);
+private:
+ std::array,3> matrix_;
+};
+
+inline std::ostream & operator<<(std::ostream & ostr, const Matrix3 & mat) {
+ ostr << mat.matrix_[0][0] << " " << mat.matrix_[0][1] << " " << mat.matrix_[0][2] << std::endl;
+ ostr << mat.matrix_[1][0] << " " << mat.matrix_[1][1] << " " << mat.matrix_[1][2] << std::endl;
+ ostr << mat.matrix_[2][0] << " " << mat.matrix_[2][1] << " " << mat.matrix_[2][2] << std::endl;
+ return ostr;
+}
+
+/// returns the project of this matrix with the PseudoJet,
+/// maintaining the 4th component of the PseudoJet unchanged
+inline PseudoJet operator*(const Matrix3 & mat, const PseudoJet & p) {
+ // r_{i} = m_{ij} p_j
+ std::array res3{{0,0,0}};
+ for (unsigned i = 0; i < 3; i++) {
+ for (unsigned j = 0; j < 3; j++) {
+ res3[i] += mat(i,j) * p[j];
+ }
+ }
+ // return a jet that maintains all internal pointers by
+ // initialising the result from the input jet and
+ // then resetting the momentum.
+ PseudoJet result(p);
+ // maintain the energy component as it was
+ result.reset_momentum(res3[0], res3[1], res3[2], p[3]);
+ return result;
+}
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+#endif // __FASTJET_CONTRIB_EEHELPERS_HH__
+
Index: contrib/contribs/LundPlane/trunk/NEWS
===================================================================
--- contrib/contribs/LundPlane/trunk/NEWS (revision 1285)
+++ contrib/contribs/LundPlane/trunk/NEWS (revision 1286)
@@ -1,4 +1,7 @@
+2021/11/08: release of version 1.0.4
+- Recursive e+e- generator with
+ spin-correlation observables
2020/02/23: release of version 1.0.3
- fixes issue with C++98 compatibility
2018/10/29: release of version 1.0.2
2018/08/30: release of version 1.0.1
Index: contrib/contribs/LundPlane/trunk/example_dpsi_slice.cc
===================================================================
--- contrib/contribs/LundPlane/trunk/example_dpsi_slice.cc (revision 0)
+++ contrib/contribs/LundPlane/trunk/example_dpsi_slice.cc (revision 1286)
@@ -0,0 +1,173 @@
+//----------------------------------------------------------------------
+/// \file example_dpsi_slice.cc
+///
+/// This example program is meant to illustrate how the
+/// fastjet::contrib::LundEEGenerator class is used.
+///
+/// Run this example with
+///
+/// \verbatim
+/// ./example_dpsi_slice < ../data/single-ee-event.dat
+/// \endverbatim
+
+//----------------------------------------------------------------------
+// $Id: example_dpsi_slice.cc 1164 2018-08-23 10:06:45Z gsalam $
+//
+// Copyright (c) 2018-, Frederic A. Dreyer, A. Karlberg, Gavin P. Salam,
+// Ludovic Scyboz, Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see .
+//----------------------------------------------------------------------
+
+#include
+#include
+#include
+
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/EECambridgePlugin.hh"
+#include
+#include "LundEEGenerator.hh" // In external code, this should be fastjet/contrib/LundEEGenerator.hh
+using namespace std;
+using namespace fastjet;
+
+// Definitions for the slice observable (ymax = 1, zcut = 0.1)
+double abs_rap_slice = 1;
+double z2_cut = 0.1;
+
+// forward declaration to make things clearer
+void read_event(vector &event);
+
+// returns true if PseudoJet p is in the central slice of half-width abs_rap_slice
+// w.r.t. the event axis pref
+bool in_slice(const PseudoJet & p, const PseudoJet & pref) {
+
+ // Rotate p, and use the rapidity to determine if p is in the slice
+ contrib::Matrix3 rotmat = contrib::Matrix3::from_direction(pref).transpose();
+ const PseudoJet p_rot = rotmat*p;
+
+ return fabs(p_rot.rap()) < abs_rap_slice;
+}
+
+//----------------------------------------------------------------------
+int main(){
+
+ //----------------------------------------------------------
+ // read in input particles
+ vector event;
+ read_event(event);
+
+ cout << "# read an event with " << event.size() << " particles" << endl;
+ //----------------------------------------------------------
+ // create an instance of RecursiveLundEEGenerator, with default options
+ int depth = -1;
+ bool dynamic_psi_reference = true;
+ fastjet::contrib::RecursiveLundEEGenerator lund(depth, dynamic_psi_reference);
+
+ // first get some C/A jets
+ double y3_cut = 1.0;
+ JetDefinition::Plugin* ee_plugin = new EECambridgePlugin(y3_cut);
+ JetDefinition jet_def(ee_plugin);
+ ClusterSequence cs(event, jet_def);
+
+ // Get the event axis (to be used for the slice)
+ std::vector excl = cs.exclusive_jets(2);
+ assert(excl.size() == 2);
+ // order the two jets according to momentum along z axis
+ if (excl[0].pz() < excl[1].pz()) {
+ std::swap(excl[0],excl[1]);
+ }
+ const PseudoJet ev_axis = excl[0]-excl[1];
+
+ // Get the list of primary declusterings
+ const vector declusts = lund.result(cs);
+
+ // find the declustering that throws something into a fixed slice
+ // (from a parent that was not in the slice) and that throws the
+ // largest pt (relative to the z axis) into that slice
+
+ int index_of_max_pt_in_slice = -1;
+ double max_pt_in_slice = 0.0;
+ double psi_1;
+ for (unsigned i = 0; i < declusts.size(); i++) {
+ const auto & decl = declusts[i];
+ if (!in_slice(decl.harder(), ev_axis) && in_slice(decl.softer(), ev_axis)) {
+ if (decl.softer().pt() > max_pt_in_slice) {
+ index_of_max_pt_in_slice = i;
+ max_pt_in_slice = decl.softer().pt();
+ psi_1 = decl.psibar();
+ }
+ }
+ }
+ if (index_of_max_pt_in_slice < 0) return 0;
+
+ // establish what we need to follow and the reference psi_1
+ int iplane_to_follow = declusts[index_of_max_pt_in_slice].leaf_iplane();
+
+ vector secondaries;
+ for (const auto & declust: declusts){
+ if (declust.iplane() == iplane_to_follow) secondaries.push_back(&declust);
+ }
+
+ int index_of_max_kt_secondary = -1;
+ double dpsi;
+ for (uint i_secondary=0; i_secondaryz() > z2_cut) {
+
+ index_of_max_kt_secondary = i_secondary;
+ double psi_2 = secondaries[i_secondary]->psibar();
+ dpsi = contrib::map_to_pi(psi_2 - psi_1);
+
+ break;
+ }
+ }
+ if (index_of_max_kt_secondary < 0) return 0;
+
+ cout << "Primary in the central slice, with Lund coordinates ( ln 1/Delta, ln kt, psibar ):" << endl;
+ pair coords = declusts[index_of_max_pt_in_slice].lund_coordinates();
+ cout << "index [" << index_of_max_pt_in_slice << "](" << coords.first << ", " << coords.second << ", "
+ << declusts[index_of_max_pt_in_slice].psibar() << ")" << endl;
+ cout << endl << "with Lund coordinates for the (highest-kT) secondary plane that passes the zcut of "
+ << z2_cut << endl;
+ coords = secondaries[index_of_max_kt_secondary]->lund_coordinates();
+ cout << "index [" << index_of_max_kt_secondary << "](" << coords.first << ", " << coords.second << ", "
+ << secondaries[index_of_max_kt_secondary]->psibar() << ")";
+
+ cout << " --> delta_psi,slice = " << dpsi << endl;
+
+ // Delete the EECambridge plugin
+ delete ee_plugin;
+
+ return 0;
+}
+
+// read in input particles
+void read_event(vector &event){
+ string line;
+ while (getline(cin, line)) {
+ istringstream linestream(line);
+ // take substrings to avoid problems when there is 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);
+ }
+}
Index: contrib/contribs/LundPlane/trunk/example_dpsi_collinear.ref
===================================================================
--- contrib/contribs/LundPlane/trunk/example_dpsi_collinear.ref (revision 0)
+++ contrib/contribs/LundPlane/trunk/example_dpsi_collinear.ref (revision 1286)
@@ -0,0 +1,19 @@
+# read an event with 70 particles
+#--------------------------------------------------------------------------
+# FastJet release 3.3.2
+# M. Cacciari, G.P. Salam and G. Soyez
+# A software package for jet finding and analysis at colliders
+# http://fastjet.fr
+#
+# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package
+# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210].
+#
+# FastJet is provided without warranty under the terms of the GNU GPLv2.
+# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code
+# and 3rd party plugin jet algorithms. See COPYING file for details.
+#--------------------------------------------------------------------------
+Primary that passes the zcut of 0.1, with Lund coordinates ( ln 1/Delta, ln kt, psibar ):
+index [0](0.461518, 2.74826, 2.32912)
+
+with Lund coordinates for the secondary plane that passes the zcut of 0.1
+index [0](0.860316, 1.51421, -2.65024) --> delta_psi = 1.30383
Index: contrib/contribs/LundPlane/trunk/Makefile
===================================================================
--- contrib/contribs/LundPlane/trunk/Makefile (revision 1285)
+++ contrib/contribs/LundPlane/trunk/Makefile (revision 1286)
@@ -1,84 +1,88 @@
# If you are using this Makefile standalone and fastjet-config is not
# in your path, edit this line to specify the full path
FASTJETCONFIG=fastjet-config
PREFIX=`$(FASTJETCONFIG) --prefix`
CXX=g++
CXXFLAGS= -O3 -Wall -g
install_script = $(SHELL) ../utils/install-sh
check_script = ../utils/check.sh
# global contrib-wide Makefile include may override some of the above
# variables (leading "-" means don't give an error if you can't find
# the file)
-include ../.Makefile.inc
#------------------------------------------------------------------------
# things that are specific to this contrib
NAME=LundPlane
-SRCS=LundGenerator.cc LundWithSecondary.cc SecondaryLund.cc
-EXAMPLES=example example_secondary
-INSTALLED_HEADERS=LundGenerator.hh LundWithSecondary.hh SecondaryLund.hh LundJSON.hh
+SRCS=LundGenerator.cc LundWithSecondary.cc SecondaryLund.cc LundEEGenerator.cc
+EXAMPLES=example example_secondary example_dpsi_collinear example_dpsi_slice
+INSTALLED_HEADERS=LundGenerator.hh LundWithSecondary.hh SecondaryLund.hh LundEEGenerator.hh LundJSON.hh
#------------------------------------------------------------------------
CXXFLAGS+= $(shell $(FASTJETCONFIG) --cxxflags)
-LDFLAGS += -lm $(shell $(FASTJETCONFIG) --libs)
+LDFLAGS += -lm $(shell $(FASTJETCONFIG) --libs) -lfastjetplugins
OBJS = $(SRCS:.cc=.o)
EXAMPLES_SRCS = $(EXAMPLES:=.cc)
install_HEADER = $(install_script) -c -m 644
install_LIB = $(install_script) -c -m 644
install_DIR = $(install_script) -d
install_DATA = $(install_script) -c -m 644
install_PROGRAM = $(install_script) -c -s
install_SCRIPT = $(install_script) -c
.PHONY: clean distclean examples check install
# compilation of the code (default target)
all: lib$(NAME).a
lib$(NAME).a: $(OBJS)
ar cru lib$(NAME).a $(OBJS)
ranlib lib$(NAME).a
# building the examples
examples: $(EXAMPLES)
# the following construct makes it possible to automatically build
# each of the examples listed in $EXAMPLES
$(EXAMPLES): % : %.o all
$(CXX) $(CXXFLAGS) -o $@ $< -L. -l$(NAME) $(LDFLAGS)
# check that everything went fine
check: examples
@for prog in $(EXAMPLES); do\
- $(check_script) $${prog} ../data/single-event.dat || exit 1; \
+ if [ "$${prog}" = "example_dpsi_collinear" ] || [ "$${prog}" = "example_dpsi_slice" ]; then \
+ $(check_script) $${prog} ../data/single-ee-event.dat || exit 1; \
+ else \
+ $(check_script) $${prog} ../data/single-event.dat || exit 1; \
+ fi; \
done
@echo "All tests successful"
# cleaning the directory
clean:
rm -f *~ *.o
distclean: clean
rm -f lib$(NAME).a $(EXAMPLES)
# install things in PREFIX/...
install: all
$(install_DIR) $(PREFIX)/include/fastjet/contrib
for header in $(INSTALLED_HEADERS); do\
$(install_HEADER) $$header $(PREFIX)/include/fastjet/contrib/;\
done
$(install_DIR) $(PREFIX)/lib
$(install_LIB) lib$(NAME).a $(PREFIX)/lib
depend:
makedepend -Y -- -- $(SRCS) $(EXAMPLES_SRCS)
# DO NOT DELETE
LundGenerator.o: LundGenerator.hh
LundWithSecondary.o: LundWithSecondary.hh LundGenerator.hh SecondaryLund.hh
SecondaryLund.o: SecondaryLund.hh LundGenerator.hh
example.o: LundGenerator.hh SecondaryLund.hh LundJSON.hh
example_secondary.o: LundWithSecondary.hh LundGenerator.hh SecondaryLund.hh
Index: contrib/contribs/LundPlane/trunk/LundEEGenerator.hh
===================================================================
--- contrib/contribs/LundPlane/trunk/LundEEGenerator.hh (revision 0)
+++ contrib/contribs/LundPlane/trunk/LundEEGenerator.hh (revision 1286)
@@ -0,0 +1,373 @@
+// $Id: LundEEGenerator.hh 1153 2018-08-22 11:56:35Z akarlberg, lscyboz $
+//
+// Copyright (c) 2018-, Frederic A. Dreyer, A. Karlberg, Gavin P. Salam,
+// Ludovic Scyboz, Gregory Soyez
+//
+//----------------------------------------------------------------------
+// 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 .
+//----------------------------------------------------------------------
+
+#ifndef __FASTJET_CONTRIB_LUNDEEGENERATOR_HH__
+#define __FASTJET_CONTRIB_LUNDEEGENERATOR_HH__
+
+#include "EEHelpers.hh"
+
+#include
+#include "fastjet/tools/Recluster.hh"
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/PseudoJet.hh"
+#include
+#include
+#include
+#include
+
+using namespace std;
+
+FASTJET_BEGIN_NAMESPACE
+
+namespace contrib{
+
+class LundEEGenerator;
+
+//----------------------------------------------------------------------
+/// \class LundEEDeclustering
+/// Contains the declustering variables associated with a single qnode
+/// on the LundEE plane
+class LundEEDeclustering {
+public:
+
+ /// return the pair PseudoJet, i.e. sum of the two subjets
+ const PseudoJet & pair() const {return pair_;}
+ /// returns the subjet with larger transverse momentum
+ const PseudoJet & harder() const {return harder_;}
+ /// returns the subjet with smaller transverse momentum
+ const PseudoJet & softer() const {return softer_;}
+
+
+ /// returns pair().m() [cached]
+ double m() const {return m_;}
+
+ /// returns the effective pseudorapidity of the emission [cached]
+ double eta() const {return eta_;}
+
+ /// returns sin(theta) of the branching [cached]
+ double sin_theta() const {return sin_theta_;}
+
+ /// returns softer().modp() / (softer().modp() + harder().modp()) [cached]
+ double z() const {return z_;}
+
+ /// returns softer().modp() * sin(theta()) [cached]
+ double kt() const {return kt_;}
+
+ /// returns ln(softer().modp() * sin(theta())) [cached]
+ double lnkt() const {return lnkt_;}
+
+ /// returns z() * Delta() [cached]
+ double kappa() const {return kappa_;}
+
+ /// returns the index of the plane to which this branching belongs
+ int iplane() const {return iplane_;}
+
+ /// returns the depth of the plane on which this declustering
+ /// occurred. 0 is the primary plane, 1 is the first set of leaves, etc.
+ int depth() const {return depth_;}
+
+ /// returns iplane (plane index) of the leaf associated with the
+ /// potential further declustering of the softer of the objects in
+ /// this splitting
+ int leaf_iplane() const {return leaf_iplane_;}
+
+ /// Returns sign_s, indicating the initial parent jet index of this splitting
+ int sign_s() const {return sign_s_;}
+
+ /// returns an azimuthal type angle between this declustering plane and the previous one
+ double psi() const {return psi_;}
+
+ /// update the azimuthal angle (deprecated)
+ void set_psi(double psi) {psi_ = psi;}
+
+ /// returns the azimuthal angle psibar between this declustering plane and the previous one
+ double psibar() const {return psibar_;}
+
+
+ /// returns the coordinates in the Lund plane
+ std::pair const lund_coordinates() const {
+ return std::pair(eta_,lnkt_);
+ }
+
+ virtual ~LundEEDeclustering() {}
+
+private:
+ int iplane_;
+ double psi_, psibar_, lnkt_, eta_;
+ double m_, z_, kt_, kappa_, sin_theta_;
+ PseudoJet pair_, harder_, softer_;
+ int depth_ = -1, leaf_iplane_ = -1;
+ int sign_s_;
+
+protected:
+ /// the constructor is private, because users will not generally be
+ /// constructing a LundEEDeclustering element themselves.
+ LundEEDeclustering(const PseudoJet& pair,
+ const PseudoJet& j1, const PseudoJet& j2,
+ int iplane = -1, double psi = 0.0, double psibar = 0.0, int depth = -1, int leaf_iplane = -1, int sign_s = 1);
+
+ friend class LundEEGenerator;
+ friend class RecursiveLundEEGenerator;
+
+};
+
+
+/// Default comparison operator for LundEEDeclustering, using kt as the ordering.
+/// Useful when including declusterings in structures like priority queues
+inline bool operator<(const LundEEDeclustering& d1, const LundEEDeclustering& d2) {
+ return d1.kt() < d2.kt();
+}
+
+//----------------------------------------------------------------------
+/// Class to carry out Lund declustering to get anything from the
+/// primary Lund plane declusterings to the full Lund diagram with all
+/// its leaves, etc.
+class RecursiveLundEEGenerator {
+ public:
+ /// constructs a RecursiveLundEEGenerator with the specified depth.
+ /// - depth = 0 means only primary declusterings are registered
+ /// - depth = 1 means the first set of leaves are declustered
+ /// - ...
+ /// - depth < 0 means no limit, i.e. recurse through all leaves
+ RecursiveLundEEGenerator(int max_depth = 0, bool dynamical_psi_ref = false) :
+ max_depth_(max_depth), nx_(1,0,0,0), ny_(0,1,0,0), dynamical_psi_ref_(dynamical_psi_ref)
+ {}
+
+ /// destructor
+ virtual ~RecursiveLundEEGenerator() {}
+
+ /// This takes a cluster sequence with an e+e- C/A style algorithm, e.g.
+ /// EECambridgePlugin(ycut=1.0).
+ ///
+ /// The output is a vector of LundEEDeclustering objects, ordered
+ /// according to kt
+ virtual std::vector result(const ClusterSequence & cs) const {
+ std::vector exclusive_jets = cs.exclusive_jets(2);
+ assert(exclusive_jets.size() == 2);
+
+ // order the two jets according to momentum along z axis
+ if (exclusive_jets[0].pz() < exclusive_jets[1].pz()) {
+ std::swap(exclusive_jets[0],exclusive_jets[1]);
+ }
+
+ PseudoJet d_ev = exclusive_jets[0] - exclusive_jets[1];
+ Matrix3 rotmat = Matrix3::from_direction(d_ev);
+
+ std::vector declusterings;
+ int depth = 0;
+ int max_iplane_sofar = 1;
+ for (unsigned ijet = 0; ijet < exclusive_jets.size(); ijet++) {
+
+ // reference direction for psibar calculation
+ PseudoJet axis = d_ev/sqrt(d_ev.modp2());
+ PseudoJet ref_plane = axis;
+
+ int sign_s = ijet == 0? +1 : -1;
+ // We can pass a vector normal to a plane of reference for phi definitions
+ append_to_vector(declusterings, exclusive_jets[ijet], depth, ijet, max_iplane_sofar,
+ rotmat, sign_s, exclusive_jets[0], exclusive_jets[1], ref_plane, 0., true);
+ }
+
+ // a typedef to save typing below
+ typedef LundEEDeclustering LD;
+ // sort so that declusterings are returned in order of decreasing
+ // kt (if result of the lambda is true, then first object appears
+ // before the second one in the final sorted list)
+ sort(declusterings.begin(), declusterings.end(),
+ [](const LD & d1, LD & d2){return d1.kt() > d2.kt();});
+
+ return declusterings;
+ }
+
+ private:
+
+ /// internal routine to recursively carry out the declusterings,
+ /// adding each one to the declusterings vector; the primary
+ /// ones are dealt with first (from large to small angle),
+ /// and then secondary ones take place.
+ void append_to_vector(std::vector & declusterings,
+ const PseudoJet & jet, int depth,
+ int iplane, int & max_iplane_sofar,
+ const Matrix3 & rotmat, int sign_s,
+ const PseudoJet & harder,
+ const PseudoJet & softer,
+ const PseudoJet & psibar_ref_plane,
+ const double & last_psibar, bool first_time) const {
+ PseudoJet j1, j2;
+ if (!jet.has_parents(j1, j2)) return;
+ if (j1.modp2() < j2.modp2()) std::swap(j1,j2);
+
+ // calculation of azimuth psi
+ Matrix3 new_rotmat;
+ if (dynamical_psi_ref_) {
+ new_rotmat = Matrix3::from_direction(rotmat.transpose()*(sign_s*jet)) * rotmat;
+ } else {
+ new_rotmat = rotmat;
+ }
+ PseudoJet rx = new_rotmat * nx_;
+ PseudoJet ry = new_rotmat * ny_;
+ PseudoJet u1 = j1/j1.modp(), u2 = j2/j2.modp();
+ PseudoJet du = u2 - u1;
+ double x = du.px() * rx.px() + du.py() * rx.py() + du.pz() * rx.pz();
+ double y = du.px() * ry.px() + du.py() * ry.py() + du.pz() * ry.pz();
+ double psi = atan2(y,x);
+
+ // calculation of psibar
+ double psibar = 0.;
+ PseudoJet n1, n2;
+
+ // First psibar for this jet
+ if (first_time) {
+
+ // Compute the angle between the planes spanned by (some axis,j1) and by (j1,j2)
+ n1 = cross_product(psibar_ref_plane,j1);
+ n2 = cross_product(j1,j2);
+
+ double signed_angle = 0.;
+ n2 /= n2.modp();
+ if (n1.modp() != 0) {
+ n1 /= n1.modp();
+ signed_angle = signed_angle_between_planes(n1,n2,j1);
+ }
+
+ psibar = map_to_pi(j1.phi() + signed_angle);
+ }
+ // Else take the value of psibar_i and the plane from the last splitting to define psibar_{i+1}
+ else {
+ n2 = cross_product(j1,j2);
+ n2 /= n2.modp();
+ psibar = map_to_pi(last_psibar + signed_angle_between_planes(psibar_ref_plane, n2, j1));
+ }
+
+ int leaf_iplane = -1;
+ // we will recurse into the softer "parent" only if the depth is
+ // not yet at its limit or if there is no limit on the depth (max_depth<0)
+ bool recurse_into_softer = (depth < max_depth_ || max_depth_ < 0);
+ if (recurse_into_softer) {
+ max_iplane_sofar += 1;
+ leaf_iplane = max_iplane_sofar;
+ }
+
+ LundEEDeclustering declust(jet, j1, j2, iplane, psi, psibar, depth, leaf_iplane, sign_s);
+ declusterings.push_back(declust);
+
+ // now recurse
+ // for the definition of psibar, we recursively pass the last splitting plane (normal to n2) and the last value
+ // of psibar
+ append_to_vector(declusterings, j1, depth, iplane, max_iplane_sofar, new_rotmat, sign_s, u1, u2, n2, psibar, false);
+ if (recurse_into_softer) {
+ append_to_vector(declusterings, j2, depth+1, leaf_iplane, max_iplane_sofar, new_rotmat, sign_s, u1, u2, n2, psibar, false);
+ }
+ }
+
+ int max_depth_ = 0;
+ /// vectors used to help define psi
+ PseudoJet nx_;
+ PseudoJet ny_;
+ bool dynamical_psi_ref_;
+};
+
+//----------------------------------------------------------------------
+/// \class LundEEGenerator
+/// Generates vector of LundEEDeclustering for a given jet
+/// corresponding to its LundEE primary plane.
+class LundEEGenerator : public FunctionOfPseudoJet< std::vector > {
+public:
+ LundEEGenerator() {}
+
+ /// destructor
+ virtual ~LundEEGenerator() {}
+
+ // definitions needed for comparison of declusterings
+ struct CompareDeclustWithKt {
+ bool operator ()(const LundEEDeclustering& d1, const LundEEDeclustering& d2) const {
+ return d1.kt() < d2.kt();
+ }
+ };
+
+ /// return a vector of the declusterings of the primary plane of the
+ /// jet, registering the supplied information (iplane) about how to label
+ /// the plane
+ virtual std::vector result_iplane(const PseudoJet& jet,
+ int iplane) const;
+
+ /// create the primary declustering sequence of a jet (with a default plane
+ /// label of -1)
+ virtual std::vector result(const PseudoJet& jet) const override{
+ return this->result_iplane(jet, -1);
+ }
+
+protected:
+ /// create a declustering from a pseudojet and its parents
+ LundEEDeclustering single_declustering(const PseudoJet& pair,
+ const PseudoJet& j1, const PseudoJet& j2,
+ PseudoJet& cross, int iplane) const;
+};
+
+//----------------------------------------------------------------------
+class LundEEGeneratorOrdered : public LundEEGenerator {
+public:
+ LundEEGeneratorOrdered() {}
+ virtual ~LundEEGeneratorOrdered() {}
+
+ /// return vector of declusterings ordered according to the defined measure
+ virtual std::vector result_iplane(const PseudoJet& jet,
+ int iplane) const override;
+
+ /// create an ordered declustering sequence of two jets
+ virtual std::vector result_twojets(const PseudoJet& jet1,
+ const PseudoJet& jet2) const;
+
+ /// create an ordered declustering sequence of multiple jets and n secondary planes
+ virtual std::vector result_jets(const std::vector& jets,
+ int n = 0) const;
+
+protected:
+ /// measure to use for ordering
+ virtual double ordering_measure(const LundEEDeclustering& d) const = 0;
+
+ /// fill a priority queue with ordered declusterings
+ virtual void fill_pq(const PseudoJet& jet, int iplane,
+ std::priority_queue< std::pair >& pq) const;
+
+ virtual std::vector pq_to_vector(std::priority_queue< std::pair >& pq) const;
+
+};
+
+class LundEEGeneratorKtOrdered : public LundEEGeneratorOrdered {
+protected:
+ /// measure to use for ordering
+ virtual double ordering_measure(const LundEEDeclustering& d) const override {return d.lnkt();};
+};
+
+} // namespace contrib
+
+FASTJET_END_NAMESPACE
+
+/// for output of declustering information
+std::ostream & operator<<(std::ostream & ostr, const fastjet::contrib::LundEEDeclustering & d);
+
+
+
+#endif // __FASTJET_CONTRIB_LUNDEEGENERATOR_HH__
+
Index: contrib/contribs/LundPlane/trunk/AUTHORS
===================================================================
--- contrib/contribs/LundPlane/trunk/AUTHORS (revision 1285)
+++ contrib/contribs/LundPlane/trunk/AUTHORS (revision 1286)
@@ -1,12 +1,30 @@
The LundPlane FastJet contrib is developed and maintained by:
Frederic Dreyer
+ Alexander Karlberg
Gavin P. Salam
+ Ludovic Scyboz
Gregory Soyez
The physics is based on:
The Lund Jet Plane.
Frederic A. Dreyer, Gavin P. Salam, and Gregory Soyez
arXiv:1807.04758
+The definition of azimuthal angles, and examples of
+all-order observables sensitive to spin correlations
+(example_dpsi_collinear and example_dpsi_slice), are from:
+
+ Spin correlations in final-state parton showers and
+ jet observables
+ A. Karlberg, G. P. Salam, L. Scyboz, and R. Verheyen
+ Eur. Phys. J. C 81, 681 (2021)
+ https://doi.org/10.1140/epjc/s10052-021-09378-0
+
+and
+
+ Soft spin correlations in final-state parton showers
+ K. Hamilton, A. Karlberg, G. P. Salam, L. Scyboz,
+ and R. Verheyen
+ arXiv:2111.01161
\ No newline at end of file
Index: contrib/contribs/LundPlane/trunk/VERSION
===================================================================
--- contrib/contribs/LundPlane/trunk/VERSION (revision 1285)
+++ contrib/contribs/LundPlane/trunk/VERSION (revision 1286)
@@ -1 +1 @@
-1.0.3
+1.0.4
Index: contrib/contribs/LundPlane/trunk/example_dpsi_collinear.cc
===================================================================
--- contrib/contribs/LundPlane/trunk/example_dpsi_collinear.cc (revision 0)
+++ contrib/contribs/LundPlane/trunk/example_dpsi_collinear.cc (revision 1286)
@@ -0,0 +1,145 @@
+//----------------------------------------------------------------------
+/// \file example_dpsi_collinear.cc
+///
+/// This example program is meant to illustrate how the
+/// fastjet::contrib::LundEEGenerator class is used.
+///
+/// Run this example with
+///
+/// \verbatim
+/// ./example_dpsi_collinear < ../data/single-ee-event.dat
+/// \endverbatim
+
+//----------------------------------------------------------------------
+// $Id: example_dpsi_collinear.cc 1164 2018-08-23 10:06:45Z gsalam $
+//
+// Copyright (c) 2018-, Frederic A. Dreyer, A. Karlberg, Gavin P. Salam,
+// Ludovic Scyboz, Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see .
+//----------------------------------------------------------------------
+
+#include
+#include
+#include
+
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/EECambridgePlugin.hh"
+#include
+#include "LundEEGenerator.hh" // In external code, this should be fastjet/contrib/LundEEGenerator.hh
+using namespace std;
+using namespace fastjet;
+
+// Definitions for the collinear observable
+double z1_cut = 0.1;
+double z2_cut = 0.1;
+
+// forward declaration to make things clearer
+void read_event(vector &event);
+
+//----------------------------------------------------------------------
+int main(){
+
+ //----------------------------------------------------------
+ // read in input particles
+ vector event;
+ read_event(event);
+
+ cout << "# read an event with " << event.size() << " particles" << endl;
+ //----------------------------------------------------------
+ // create an instance of RecursiveLundEEGenerator, with default options
+ int depth = -1;
+ bool dynamic_psi_reference = true;
+ fastjet::contrib::RecursiveLundEEGenerator lund(depth, dynamic_psi_reference);
+
+ // first get some C/A jets
+ double y3_cut = 1.0;
+ JetDefinition::Plugin* ee_plugin = new EECambridgePlugin(y3_cut);
+ JetDefinition jet_def(ee_plugin);
+ ClusterSequence cs(event, jet_def);
+
+ // Get the list of primary declusterings
+ const vector declusts = lund.result(cs);
+
+ // Find the highest-kt primary declustering (ordered in kt by default)
+ int i_primary = -1;
+ double psi_1;
+ for (int i=0; i z1_cut) {
+ i_primary = i;
+ psi_1 = declusts[i].psibar();
+ break;
+ }
+ }
+ if (i_primary < 0) return 0;
+
+ // Find the highest-kt secondary associated to that Lund leaf
+ int iplane_to_follow = declusts[i_primary].leaf_iplane();
+ vector secondaries;
+ for (const auto & declust: declusts){
+ if (declust.iplane() == iplane_to_follow) secondaries.push_back(&declust);
+ }
+ if(secondaries.size() < 1) return 0;
+
+ int i_secondary = -1;
+ double dpsi;
+ for (int i=0; iz() > z2_cut) {
+ i_secondary = i;
+ double psi_2 = secondaries[i]->psibar();
+ dpsi = contrib::map_to_pi(psi_2-psi_1);
+ break;
+ }
+ }
+ if (i_secondary < 0) return 0;
+
+ cout << "Primary that passes the zcut of "
+ << z1_cut << ", with Lund coordinates ( ln 1/Delta, ln kt, psibar ):" << endl;
+ pair coords = declusts[i_primary].lund_coordinates();
+ cout << "index [" << i_primary << "](" << coords.first << ", " << coords.second << ", "
+ << declusts[i_primary].psibar() << ")" << endl;
+ cout << endl << "with Lund coordinates for the secondary plane that passes the zcut of "
+ << z2_cut << endl;
+ coords = secondaries[i_secondary]->lund_coordinates();
+ cout << "index [" << i_secondary << "](" << coords.first << ", " << coords.second << ", "
+ << secondaries[i_secondary]->psibar() << ")";
+
+ cout << " --> delta_psi = " << dpsi << endl;
+
+ // Delete the EECambridge plugin
+ delete ee_plugin;
+
+ return 0;
+}
+
+// read in input particles
+void read_event(vector &event){
+ string line;
+ while (getline(cin, line)) {
+ istringstream linestream(line);
+ // take substrings to avoid problems when there is 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);
+ }
+}
Index: contrib/contribs/LundPlane/trunk/README
===================================================================
--- contrib/contribs/LundPlane/trunk/README (revision 1285)
+++ contrib/contribs/LundPlane/trunk/README (revision 1286)
@@ -1,46 +1,62 @@
------------------------------------------------------------------------
LundPlane FastJet contrib
------------------------------------------------------------------------
The LundPlane FastJet contrib provides tools necessary to construct
the Lund Plane of a jet.
It is based on arXiv:1807.04758 by Frederic Dreyer, Gavin Salam, and
Gregory Soyez.
+Azimuthal-angle definitions are taken from arXiv:2103.16526, by
+Alexander Karlberg, Gavin Salam, Ludovic Scyboz and Rob Verheyen, and
+arXiv:2111.01161 by the authors above and Keith Hamilton.
+
It provides the following C++ classes:
- LundDeclustering
Contains the declustering variables associated with a node on
the Lund plane.
- LundGenerator
A tool to generate a vector of LundDeclustering for a given jet
corresponding to its primary Lund plane.
+- LundEEGenerator
+ A recursive tool to extract declusterings for a given jet in
+ e+e- collisions, i.e. its full Lund structure down to a user-given
+ depth
+
- LundWithSecondary
A tool to generate primary, and secondary, lund planes of a jet
according to the definition given as input.
- SecondaryLund, SecondaryLund_{mMDT|dotmMDT|Mass}
Provides a definition for the leading emission.
-There are two C++ examples which can be compiled with
+There are four C++ examples which can be compiled with
> make example example_secondary
The example program can be run using
> ./example < ../data/single-event.dat
and will output the primary lund plane declusterings as well as
generating a jets.json file, which can be read by a python script:
> python3 ./example.py
Some python classes to assist with reading and processing the json file
are provided in the read_lund_json.py module.
The example_secondary program outputs both the primary and secondary
Lund plane declusterings to screen, and can be run with:
-> ./example_secondary < ../data/single-event.dat
\ No newline at end of file
+> ./example_secondary < ../data/single-event.dat
+
+The example_dpsi_collinear and example_dpsi_slice programs compute the
+contribution to the collinear azimuthal delta_psi observable from
+2103.16526, and to the slice observable delta_psi_slice from 2111.01161,
+respectively. It can be run with
+
+> ./example_dpsi_collinear < ../data/single-ee-event.dat
Index: contrib/contribs/LundPlane/trunk/example_dpsi_slice.ref
===================================================================
--- contrib/contribs/LundPlane/trunk/example_dpsi_slice.ref (revision 0)
+++ contrib/contribs/LundPlane/trunk/example_dpsi_slice.ref (revision 1286)
@@ -0,0 +1,19 @@
+# read an event with 70 particles
+#--------------------------------------------------------------------------
+# FastJet release 3.3.2
+# M. Cacciari, G.P. Salam and G. Soyez
+# A software package for jet finding and analysis at colliders
+# http://fastjet.fr
+#
+# Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package
+# for scientific work and optionally PLB641(2006)57 [hep-ph/0512210].
+#
+# FastJet is provided without warranty under the terms of the GNU GPLv2.
+# It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code
+# and 3rd party plugin jet algorithms. See COPYING file for details.
+#--------------------------------------------------------------------------
+Primary in the central slice, with Lund coordinates ( ln 1/Delta, ln kt, psibar ):
+index [0](0.461518, 2.74826, 2.32912)
+
+with Lund coordinates for the (highest-kT) secondary plane that passes the zcut of 0.1
+index [0](0.860316, 1.51421, -2.65024) --> delta_psi,slice = 1.30383
Index: contrib/contribs/LundPlane/trunk/LundEEGenerator.cc
===================================================================
--- contrib/contribs/LundPlane/trunk/LundEEGenerator.cc (revision 0)
+++ contrib/contribs/LundPlane/trunk/LundEEGenerator.cc (revision 1286)
@@ -0,0 +1,222 @@
+// $Id: LundEEGenerator.cc 1153 2018-08-22 11:56:35Z frdreyer $
+//
+// Copyright (c) 2018-, Frederic A. Dreyer, A. Karlberg, Gavin P. Salam,
+// Ludovic Scyboz, Gregory Soyez
+//
+//----------------------------------------------------------------------
+// This file is part of FastJet contrib.
+//
+// It is free software; you can redistribute it and/or modify it under
+// the terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2 of the License, or (at
+// your option) any later version.
+//
+// It is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+// License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this code. If not, see .
+//----------------------------------------------------------------------
+
+#include
+#include "LundEEGenerator.hh"
+
+using namespace std;
+
+FASTJET_BEGIN_NAMESPACE
+
+namespace contrib{
+
+LundEEDeclustering::LundEEDeclustering(const PseudoJet& pair,
+ const PseudoJet& j1, const PseudoJet& j2,
+ int iplane, double psi, double psibar,
+ int depth, int leaf_iplane, int sign_s)
+ : iplane_(iplane), psi_(psi), psibar_(psibar), m_(pair.m()), pair_(pair), depth_(depth), leaf_iplane_(leaf_iplane), sign_s_(sign_s) {
+
+ double omc = one_minus_costheta(j1,j2);
+
+ if (omc > sqrt(numeric_limits::epsilon())) {
+ double cos_theta = 1.0 - omc;
+ double theta = acos(cos_theta);
+ sin_theta_ = sin(theta);
+ eta_ = -log(tan(theta/2.0));
+ } else {
+ // we are at small angles, so use small-angle formulas
+ double theta = sqrt(2. * omc);
+ sin_theta_ = theta;
+ eta_ = -log(theta/2);
+ }
+
+ // establish which of j1 and j2 is softer
+ if (j1.modp2() > j2.modp2()) {
+ harder_ = j1;
+ softer_ = j2;
+ } else {
+ harder_ = j2;
+ softer_ = j1;
+ }
+
+ // now work out the various LundEE declustering variables
+ double softer_modp = softer_.modp();
+ z_ = softer_modp / (softer_modp + harder_.modp());
+ kt_ = softer_modp * sin_theta_;
+ lnkt_ = log(kt_);
+ kappa_ = z_ * sin_theta_;
+}
+
+//----------------------------------------------------------------------
+/// return a vector of the declusterings of the primary plane of the
+/// jet, registering information about the plane from which this
+/// comes
+std::vector LundEEGenerator::result_iplane(const PseudoJet& jet, int iplane) const {
+ std::vector result;
+ PseudoJet j = jet;
+ PseudoJet pair, j1, j2, cross;
+ pair = j;
+ while (pair.has_parents(j1, j2)) {
+ // order the jets in their 3-mom norm
+ if (j1.modp2() < j2.modp2()) std::swap(j1,j2);
+ // on the first invocation cross has its default value of 0;
+ result.push_back(single_declustering(pair, j1, j2, cross, iplane));
+ pair = j1;
+ }
+ return result;
+}
+
+//----------------------------------------------------------------------
+/// get declustering corresponding to a given pseudojet and its parents.
+///
+/// If cross is a null vector, then the psi value for the declustering
+/// will be set to zero and the cross vector will be defined as
+/// j1 x j2.
+///
+/// If cross is non-zero, then the psi value is determined as the
+/// angle between this j1 x j2 and the existing cross.
+LundEEDeclustering LundEEGenerator::single_declustering(const PseudoJet& pair,
+ const PseudoJet& j1,
+ const PseudoJet& j2,
+ PseudoJet& cross,
+ int iplane) const {
+ // get the angle between current declustering plane and previous one
+ PseudoJet j1xj2 = cross_product(j1,j2,false);
+ // the first time one encounters the next line cross == 0
+ // (default PJ initialisation), so one simply defines the direction
+ // meant by psi=0
+ double psi;
+ if (cross == double(0)) {
+ psi = 0;
+ cross = j1xj2;
+ // normalise it
+ cross /= cross.modp();
+ } else {
+ psi = theta(cross, j1xj2);
+ }
+
+ // prepare for next iteration
+ // return declustering element
+ return LundEEDeclustering(pair, j1, j2, iplane, psi);
+}
+
+//----------------------------------------------------------------------
+/// create an ordered declustering sequence for one jet
+std::vector LundEEGeneratorOrdered::result_iplane(const PseudoJet& jet,
+ int iplane) const {
+ std::priority_queue< std::pair > pq_declusts;
+ this->fill_pq(jet, iplane, pq_declusts);
+ return this->pq_to_vector(pq_declusts);
+}
+
+//----------------------------------------------------------------------
+/// create an ordered declustering sequence for two jets
+std::vector LundEEGeneratorOrdered::result_twojets(const PseudoJet& jet1,
+ const PseudoJet& jet2) const {
+ std::priority_queue< std::pair > pq_declusts;
+ this->fill_pq(jet1, 0, pq_declusts);
+ this->fill_pq(jet2, 1, pq_declusts);
+ return this->pq_to_vector(pq_declusts);
+}
+
+//----------------------------------------------------------------------
+/// create an ordered declustering sequence for multiple jets, along with n secondary planes
+std::vector LundEEGeneratorOrdered::result_jets(const std::vector& jets,
+ int n) const{
+ std::priority_queue< std::pair > pq_declusts;
+ unsigned int i=0;
+ for(; i < jets.size(); ++i)
+ this->fill_pq(jets[i], i, pq_declusts);
+ std::vector< std::pair > tmp;
+ // get the n elements from which we want secondary planes
+ for(int isec=0; isec < n; ++isec) {
+ tmp.push_back(pq_declusts.top());
+ pq_declusts.pop();
+ }
+ // get the n secondary declustering sequences and insert them into priority queue
+ for(int isec=0; isec < n; ++isec) {
+ this->fill_pq(tmp[isec].second.softer(), isec+i, pq_declusts);
+ pq_declusts.push(tmp[isec]);
+ }
+ // return sorted vector of all the declusterings
+ return this->pq_to_vector(pq_declusts);
+}
+
+//----------------------------------------------------------------------
+/// fill a priority queue with ordered declusterings
+void LundEEGeneratorOrdered::fill_pq(const PseudoJet& jet, int iplane,
+ std::priority_queue< std::pair >& pq) const {
+
+ PseudoJet j = jet;
+ PseudoJet pair, j1, j2, cross;
+ pair = j;
+ while (pair.has_parents(j1, j2)) {
+ // order the jets in their 3-mom norm
+ if (j1.modp2() < j2.modp2()) std::swap(j1,j2);
+ LundEEDeclustering d = single_declustering(pair, j1, j2, cross, iplane);
+ pq.push(std::make_pair(this->ordering_measure(d), d));
+ pair = j1;
+ }
+}
+
+//----------------------------------------------------------------------
+/// turn a priority queue to an ordered vector
+std::vector LundEEGeneratorOrdered::pq_to_vector(std::priority_queue< std::pair >& pq) const {
+ std::vector result;
+ PseudoJet cross;
+ while (!pq.empty()) {
+ LundEEDeclustering d = pq.top().second;
+ // update the psi angle according to the new ordering
+ PseudoJet newcross = PseudoJet(d.harder().py()*d.softer().pz()
+ - d.harder().pz()*d.softer().py(),
+ d.harder().pz()*d.softer().px()
+ - d.harder().px()*d.softer().pz(),
+ d.harder().px()*d.softer().py()
+ - d.harder().py()*d.softer().px(), 0.0);
+ if (cross != double(0)) d.set_psi(theta(cross, newcross));
+ cross = newcross;
+ // add declustering to result vector
+ result.push_back(d);
+ pq.pop();
+ }
+ return result;
+}
+
+} // namespace contrib
+
+
+FASTJET_END_NAMESPACE
+
+std::ostream & operator<<(std::ostream & ostr, const fastjet::contrib::LundEEDeclustering & d) {
+ ostr << "kt = " << d.kt()
+ << " z = " << d.z()
+ << " eta = " << d.eta()
+ << " psi = " << d.psi()
+ << " psibar = " << d.psibar()
+ << " m = " << d.m()
+ << " iplane = " << d.iplane()
+ << " depth = " << d.depth()
+ << " leaf_iplane = " << d.leaf_iplane();
+ return ostr;
+}
+