Index: contrib/contribs/LundPlane/trunk/example_dpsi_collinear.ref =================================================================== --- contrib/contribs/LundPlane/trunk/example_dpsi_collinear.ref (revision 1362) +++ contrib/contribs/LundPlane/trunk/example_dpsi_collinear.ref (revision 1363) @@ -1,19 +1,22 @@ # read an event with 70 particles #-------------------------------------------------------------------------- -# FastJet release 3.3.2 +# FastJet release 3.4.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. +# FastJet is provided without warranty under the GNU GPL v2 or higher. # 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) +Highest-kt primary that passes the zcut of 0.1, with Lund coordinates ( ln 1/Delta, ln kt, psibar ): +index [0](0.461518, 2.74826, -2.25293) 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 [0](0.860316, 1.51421, -0.949104) --> delta_psi(primary,secondary) = 1.30383 + +with Lund coordinates for the second primary plane that passes the zcut of 0.1 +index [2](1.99207, 1.49759, -1.91726) --> delta_psi(primary_1, primary_2) = 0.335672 Index: contrib/contribs/LundPlane/trunk/python/gen-test-event.py =================================================================== --- contrib/contribs/LundPlane/trunk/python/gen-test-event.py (revision 0) +++ contrib/contribs/LundPlane/trunk/python/gen-test-event.py (revision 1363) @@ -0,0 +1,38 @@ +#!/usr/bin/env python3 +# +# Script to help produce simple small events for testing the +# azimuthal angle aspects of the LundPlane contrib +# +# Usage (from the LundPlane directory): +# +# python3 python/gen-test-event.py | ./example_dpsi_collinear +# + +import fastjet as fj +from math import * + +# produce an event with a collinear branching on each side +# and some angle phi between the two branches. +# +# - smaller value of scale make the event more collinear +# - phi = 0 causes the soft particles to be on opposite sides +scale = 0.01 +phi = 0.1 +particles = [ + fj.PseudoJet(0,0,+1,1), + fj.PseudoJet(0,0,-1,1), + fj.PtYPhiM(0.12*scale, 1.0*(1 - log(scale)), 0.0), + fj.PtYPhiM(0.10*scale, -1.0*(1 - log(scale)), -pi - phi), +] + +# dumb way to get sum of all momenta -- clustering with e+e- Cambridge/Aachen +# with a large radius +psum = fj.JetDefinition(fj.ee_genkt_algorithm, 2*pi, 0.0)(particles)[0] +#print(psum) + +# then boost to balance the event +for p in particles: + p.unboost(psum) + +for p in particles: + print(p.px(), p.py(), p.pz(), p.E()) \ No newline at end of file Property changes on: contrib/contribs/LundPlane/trunk/python/gen-test-event.py ___________________________________________________________________ Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: contrib/contribs/LundPlane/trunk/example_dpsi_collinear.cc =================================================================== --- contrib/contribs/LundPlane/trunk/example_dpsi_collinear.cc (revision 1362) +++ contrib/contribs/LundPlane/trunk/example_dpsi_collinear.cc (revision 1363) @@ -1,145 +1,166 @@ //---------------------------------------------------------------------- /// \file example_dpsi_collinear.cc /// /// This example program is meant to illustrate how the /// fastjet::contrib::RecursiveLundEEGenerator class is used. /// /// Run this example with /// /// \verbatim /// ./example_dpsi_collinear < ../data/single-ee-event.dat /// \endverbatim //---------------------------------------------------------------------- // $Id$ // // Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, // Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // 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 "RecursiveLundEEGenerator.hh" // In external code, this should be fastjet/contrib/RecursiveLundEEGenerator.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); + fastjet::contrib::RecursiveLundEEGenerator lund(depth); // 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; + // Find the two highest-kt primary declustering (ordered in kt by default) + int i_primary_1 = -1, i_primary_2 = -1; + double psi_primary_1, psi_primary_2; + double dpsi_11; for (unsigned int i=0; i z1_cut) { - i_primary = i; - psi_1 = declusts[i].psibar(); - break; + if (i_primary_1 < 0) { + i_primary_1 = i; + psi_primary_1 = declusts[i].psibar(); + } else { + i_primary_2 = i; + psi_primary_2 = declusts[i].psibar(); + dpsi_11 = contrib::lund_plane::map_to_pi(psi_primary_2-psi_primary_1); + break; + } } } - if (i_primary < 0) return 0; + // If no primary at all, just return + if (i_primary_1 < 0) return 0; - // Find the highest-kt secondary associated to that Lund leaf - int iplane_to_follow = declusts[i_primary].leaf_iplane(); + // For the highest-kt primary: find the highest-kt secondary associated to that Lund leaf + // that passes the zcut of z2_cut + int iplane_to_follow = declusts[i_primary_1].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 (unsigned int i=0; iz() > z2_cut) { - i_secondary = i; - double psi_2 = secondaries[i]->psibar(); - dpsi = contrib::lund_plane::map_to_pi(psi_2-psi_1); - break; + double dpsi_12; + if (secondaries.size() > 0) { + for (unsigned int i=0; iz() > z2_cut) { + i_secondary = i; + double psi_2 = secondaries[i]->psibar(); + dpsi_12 = contrib::lund_plane::map_to_pi(psi_2-psi_primary_1); + break; + } } } - if (i_secondary < 0) return 0; - cout << "Primary that passes the zcut of " + cout << "Highest-kt 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() << ")"; + pair coords = declusts[i_primary_1].lund_coordinates(); + cout << "index [" << i_primary_1 << "](" << coords.first << ", " << coords.second << ", " + << declusts[i_primary_1].psibar() << ")" << endl; + + // dpsi_12 is the angle between the primary and the secondary + if (i_secondary >= 0) { + 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(primary,secondary) = " << dpsi_12 << endl; + } + // dpsi_11 is the angle between the two primaries + if (i_primary_2 >= 0) { + cout << endl << "with Lund coordinates for the second primary plane that passes the zcut of " + << z1_cut << endl; + coords = declusts[i_primary_2].lund_coordinates(); + cout << "index [" << i_primary_2 << "](" << coords.first << ", " << coords.second << ", " + << declusts[i_primary_2].psibar() << ")"; - cout << " --> delta_psi = " << dpsi << endl; + cout << " --> delta_psi(primary_1, primary_2) = " << dpsi_11 << 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/ChangeLog =================================================================== --- contrib/contribs/LundPlane/trunk/ChangeLog (revision 1362) +++ contrib/contribs/LundPlane/trunk/ChangeLog (revision 1363) @@ -1,143 +1,169 @@ +2024-01-10 Gavin Salam + Ludo Scyboz + Alexander Karlberg + + * LundEEHelpers.hh: + extended the comments in the signed_angle_between_planes function + + * RecursiveLundEEGenerator.hh: + changed definition of psibar, which fixes a bug in differences + of psibar values between opposite hemispheres. Differences + within a same hemisphere should stay the same. + + * example_dpsi_collinear.cc: + added extra output showing the primary declustering in the + second hemisphere. Also eliminated deprecated + dynamic_psi_reference argument in setting up the RecursiveLundEEGenerator. + + * example_dpsi_collinear.ref: + * example_dpsi_slice.ref: + updated these reference files so that are consistent with + the new definition of psibar. Note the delta psibar values + do not change. The collinear reference has acquired extra + output showing a primary declustering in the second hemisphere. + + * python/gen-test-event.py: *** ADDED *** + small script to generate test events for verifying output from + example_dpsi_collinear.ref + 2024-01-08 Gavin Salam * example_dpsi_slice.cc: replaced a uint -> unsigned int, to resolve a compilation issue with gcc-13 on MacOS 2023-09-22 Fri * NEWS: * VERSION: preparing 2.0.4 release * RecursiveLundEEGenerator.hh: fixed const issue in in ternal lambda function (reported by Seryubin Seraphim) * LundPlane.hh: *** ADDED *** added this to provide a generic include for the main classes (absence of such a standard-looking reported by Seryubin Seraphim) 2022-10-04 Gregory Soyez * example.cc: * example_secondary.cc: * example_dpsi_collinear.cc: * example_dpsi_slice.cc: removed (harmless) compilation warnings (+ minor uniformisation of indentation) 2022-10-04 Tue + Ludo Scyboz * VERSION: * NEWS: prepared for release 2.0.3 * EEHelpers.hh -> LundEEHelpers.hh: EEHelpers.hh was not being installed; given that it's a potentially common name, renamed it before including among the installation targets. Also placed whole contents in a new contrib::lund_plane namespace, because of certain potentially common names like cross_product * Makefile: replaced EEHelpers.hh -> LundEEHelpers.hh and made sure that LundEEHelpers.hh gets installed (without which RecursiveLundEEGenerator.hh cannot be used) * example_dpsi_collinear.cc: * example_dpsi_slice.cc: use of things from LundEEHelpers.hh updated to use of new namespace * RecursiveLundEEGenerator.hh: * RecursiveLundEEGenerator.cc: updated to use LundEEHelpers.hh (and associated namespace); also added a (protected) default constructor to LundEEDeclustering 2022-08-20 Gavin Salam * VERSION: * NEWS: prepared for release 2.0.2 * Makefile: updated some dependencies * EEHelpers.hh: added #include , as per request from Andy Buckley for compilation with g++-12 on some systems 2021-12-06 Gavin Salam * NEWS: * VERSION: prepared for release 2.0.1 * AUTHORS: fixed missing names and publication info 2021-12-06 Gregory Soyez * SecondaryLund.cc: fixed int v. unsigned int in loop over vector indices 2021-12-06 Gavin Salam * example.py: fixed name of executable in comments about how to execute this (thanks to Matteo Cacciari) 2021-11-09 Ludovic Scyboz * VERSION: preparing for release of 2.0.0 * RecursiveLundEEGenerator.hh: * RecursiveLundEEGenerator.cc: class for recursive Lund declustering in e+e- * example_dpsi_collinear.cc: spin-sensitive collinear observable from 2103.16526 * example_dpsi_slice.cc: spin-sensitive non-global observable from 2111.01161 2020-02-23 Gavin Salam * NEWS: * VERSION: preparing for release of 1.0.3 * example.cc: changed outfile open(filename) to outfile.open(filename.c_str()); to attempt to solve issue reported by Steven Schramm. 2018-10-26 Gavin Salam * read_lund_json.py: removed extraneous normalisation of zeroth bin in the LundImage class. Added documentation. 2018-08-30 Gavin Salam * VERSION: * NEWS: Release of version 1.0.1 2018-08-23 Gavin Salam * LundWithSecondary.hh: * LundWithSecondary.cc: added secondary_index(...), removed virtual qualifier from various functions * example_secondary.cc: * example_secondary.ref: example now prints out index of the primary declustering being used for the secondary. Referemce file updated accordingly. 2018-08-09 Frédéric Dreyer First version of LundPlane. Index: contrib/contribs/LundPlane/trunk/LundEEHelpers.hh =================================================================== --- contrib/contribs/LundPlane/trunk/LundEEHelpers.hh (revision 1362) +++ contrib/contribs/LundPlane/trunk/LundEEHelpers.hh (revision 1363) @@ -1,265 +1,270 @@ // $Id$ // // Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, // Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // // 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 #include FASTJET_BEGIN_NAMESPACE namespace contrib{ namespace lund_plane { //---------------------------------------------------------------------- /// 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. +/// Get the angle between two planes defined by normalized vectors +/// n1, n2. The sign is decided by the direction of a vector n, such that +/// if n1 x n2 points in the same direction as n, the angle is positive. +/// The result is in the range -pi < theta < pi. +/// +/// For example, labelling (x,y,z), and taking n1 = (1,0,0), +/// n2 = (0,1,0), n = (0,0,1), then the angle is +pi/2. 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 lund_plane } // namespace contrib FASTJET_END_NAMESPACE #endif // __FASTJET_CONTRIB_EEHELPERS_HH__ Index: contrib/contribs/LundPlane/trunk/example_dpsi_slice.ref =================================================================== --- contrib/contribs/LundPlane/trunk/example_dpsi_slice.ref (revision 1362) +++ contrib/contribs/LundPlane/trunk/example_dpsi_slice.ref (revision 1363) @@ -1,19 +1,19 @@ # read an event with 70 particles #-------------------------------------------------------------------------- -# FastJet release 3.3.2 +# FastJet release 3.4.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. +# FastJet is provided without warranty under the GNU GPL v2 or higher. # 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) +index [0](0.461518, 2.74826, -2.25293) 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 [0](0.860316, 1.51421, -0.949104) --> delta_psi,slice = 1.30383 Index: contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.hh =================================================================== --- contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.hh (revision 1362) +++ contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.hh (revision 1363) @@ -1,301 +1,347 @@ // $Id$ // // Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, // Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen // //---------------------------------------------------------------------- // 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_RECURSIVELUNDEEGENERATOR_HH__ #define __FASTJET_CONTRIB_RECURSIVELUNDEEGENERATOR_HH__ #include "LundEEHelpers.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 LundEEDeclustering /// Contains the declustering variables associated with a single node /// 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 angle psibar associated with this declustering. + /// The actual value of psibar is arbitrary, but differences in psibar + /// values between different clusterings are meaningful. + double psibar() const {return psibar_;} + /// (DEPRECATED) /// returns an azimuthal type angle between this declustering plane and the previous one /// Note: one should use psibar() instead, since we found that this definition of psi is /// not invariant under rotations of the event 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: /// default ctor (protected, should not normally be needed by users, /// but can be useful for derived classes) LundEEDeclustering() {} /// the constructor is protected, 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 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 + /// + /// The psibar values that are set in the result Lund tree have the + /// following property: + /// + /// - if the jet with the larger pz has splittings, then its + /// first splitting has psibar = 0 + /// - otherwise the first splitting of the other jet has psibar = 0 + /// + /// Note that this makes psibar IR unsafe (because an arbitrarily soft + /// splitting can be the one that gets the reference psibar=0 value), + /// but differences between psibar values are IR safe. + /// + /// NB: The dynamical_psi_ref option relates to the deprecated definition of psi + /// New code should use the psibar() function and dynamical_psi_ref + /// is irrelevant. 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]; lund_plane::Matrix3 rotmat = lund_plane::Matrix3::from_direction(d_ev); std::vector declusterings; int depth = 0; int max_iplane_sofar = 1; + +// 2024-01: new code, that fixes up issue of psibar differences +// between hemispheres. If RLEEG_NEWPSIBAR is false, answers +// will come out wrong. +#define RLEEG_NEWPSIBAR +#ifdef RLEEG_NEWPSIBAR + // 2024-01 -- attempt at new definition of psibar + PseudoJet ref_plane; + double last_psibar = 0.; + bool first_time = true; + + for (unsigned ijet = 0; ijet < exclusive_jets.size(); ijet++) { + int sign_s = ijet == 0? +1 : -1; + append_to_vector(declusterings, exclusive_jets[ijet], depth, ijet, max_iplane_sofar, + rotmat, sign_s, ref_plane, last_psibar, first_time); + } +#else 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; + bool // 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); + rotmat, sign_s, ref_plane, 0., true); } - +#endif + // 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, const 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 lund_plane::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 & 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 lund_plane::Matrix3 new_rotmat; if (dynamical_psi_ref_) { new_rotmat = lund_plane::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) { +#ifdef RLEEG_NEWPSIBAR +// 2024-01: new code, that fixes up issue of psibar differences +// between hemispheres + assert(last_psibar == 0.0); + psibar = 0.0; + n2 = lund_plane::cross_product(j1,j2); + n2 /= n2.modp(); + psibar_ref_plane = n2; + first_time = false; +#else // Compute the angle between the planes spanned by (some axis,j1) and by (j1,j2) n1 = lund_plane::cross_product(psibar_ref_plane,j1); n2 = lund_plane::cross_product(j1,j2); double signed_angle = 0.; n2 /= n2.modp(); if (n1.modp() != 0) { n1 /= n1.modp(); signed_angle = lund_plane::signed_angle_between_planes(n1,n2,j1); } psibar = lund_plane::map_to_pi(j1.phi() + signed_angle); +#endif } // Else take the value of psibar_i and the plane from the last splitting to define psibar_{i+1} else { n2 = lund_plane::cross_product(j1,j2); n2 /= n2.modp(); psibar = lund_plane::map_to_pi(last_psibar + lund_plane::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); + bool lcl_first_time = false; + append_to_vector(declusterings, j1, depth, iplane, max_iplane_sofar, new_rotmat, sign_s, n2, psibar, lcl_first_time); 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); + append_to_vector(declusterings, j2, depth+1, leaf_iplane, max_iplane_sofar, new_rotmat, sign_s, n2, psibar, lcl_first_time); } } int max_depth_ = 0; /// vectors used to help define psi PseudoJet nx_; PseudoJet ny_; bool dynamical_psi_ref_; }; } // namespace contrib FASTJET_END_NAMESPACE /// for output of declustering information std::ostream & operator<<(std::ostream & ostr, const fastjet::contrib::LundEEDeclustering & d); #endif // __FASTJET_CONTRIB_RECURSIVELUNDEEGENERATOR_HH__