Index: contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.hh =================================================================== --- contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.hh (revision 1365) +++ contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.hh (revision 1366) @@ -1,349 +1,362 @@ // $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 + /// The actual value of psibar is arbitrary (but see below), while differences in psibar /// values between different clusterings are meaningful. + /// + /// The absolute value of psibar should in general not be used, but + /// it has the following properties. (1) the highest kt splitting ("ref") in the + /// positive-z hemisphere defined psibar=0 (if not then the highest kt + /// splitting in the negative-z hemisphere does). If the +z jet aligns along + /// the z axis, then psibar of collinear splitting i in that jet is + /// + /// psibar_i = phi_i - phi_ref + /// + /// For the negative-z jet, the psibar of collinear splitting i is + /// + /// psibar_i = phi_i - pi - phi_ref + /// 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 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, 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, 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} // The signed angle is multiplied by sign_s (+1 for the +z hemisphere, -1 for the -z hemisphere) // to take the correct orientation into account. else { n2 = lund_plane::cross_product(j1,j2); n2 /= n2.modp(); psibar = lund_plane::map_to_pi(last_psibar + sign_s*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 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, 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__ Index: contrib/contribs/LundPlane/trunk/python/gen-test-event.py =================================================================== --- contrib/contribs/LundPlane/trunk/python/gen-test-event.py (revision 1365) +++ contrib/contribs/LundPlane/trunk/python/gen-test-event.py (revision 1366) @@ -1,56 +1,87 @@ #!/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 * +import argparse # produce an event with a collinear branching on each side, # some angle phi between the two branches, and an additional # particle in the z>0 hemisphere at phi!=0 (to test whether # signs are set properly across the two hemispheres). # # (That particle generates the highest-kt primary on that side, # but that primary does not pass the z-cut. Hence the two # opposite-side primaries will have psibar != 0) # # - smaller value of scale make the event more collinear # - phi = 0 causes the soft particles to be on opposite sides # - mirror_event rotates the event by pi (to test whether # seeding the declustering with the other jet makes a difference) -scale = 0.001 -phi = 0.1 -mirror_event = False + +argparser = argparse.ArgumentParser() +argparser.add_argument("--scale", type=float, default=0.001) +argparser.add_argument("--phi", type=float, default=0.1) +argparser.add_argument("--mirror", action="store_true") +argparser.add_argument("--phi-ref", type=float, default=1.15) +args = argparser.parse_args() + +scale = args.scale +phi = args.phi +mirror_event = args.mirror + + +rap = lambda pt,z: -log(pt/z/2) + +z1 = 0.12 +pt1 = 0.007*scale +y1 = rap(pt1, z1) +phi1 = args.phi + +# if args.phi_ref == 0, then we expect that +# psibar2 = phi2 - pi +z2 = 0.17 +pt2 = 0.004*scale +y2 = -rap(pt2, z2) +phi2 = pi + +# highest kt emission, which will set psibar = 0 +z3 = 0.05 +pt3 = scale +y3 = rap(pt3, z3) +phi3 = args.phi_ref particles_in = [ fj.PseudoJet(0,0,+1,1), fj.PseudoJet(0,0,-1,1), - fj.PtYPhiM(0.12*scale, 1.0*(1 - log(scale)), -phi), - fj.PtYPhiM(0.10*scale, -1.0*(1 - log(scale)), pi), - fj.PtYPhiM(0.15*scale, 0.92*(1 - log(scale)), 1.15), + fj.PtYPhiM(pt1, y1, phi1), + # the psibar that will come out here should be + fj.PtYPhiM(pt2, y2, phi2), + fj.PtYPhiM(pt3, y3, phi3), ] # 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_in)[0] #print(psum) # then boost to balance the event for p in particles_in: p.unboost(psum) # and mirror the event if requested particles = [] for p in particles_in: if mirror_event: p = fj.PseudoJet(-p.px(), -p.py(), -p.pz(), p.E()) particles.append(p) for p in particles: - print(p.px(), p.py(), p.pz(), p.E()) + print(p.px(), p.py(), p.pz(), p.E(), p.rap()) Index: contrib/contribs/LundPlane/trunk/ChangeLog =================================================================== --- contrib/contribs/LundPlane/trunk/ChangeLog (revision 1365) +++ contrib/contribs/LundPlane/trunk/ChangeLog (revision 1366) @@ -1,178 +1,193 @@ +2024-01-18 Ludo Scyboz (+ minor contribs from Gregory Soyez + Gavin Salam) + + * RecursiveLundEEGenerator.hh: + * example_dpsi_collinear.ref: + * example_dpsi_slice.ref: + * python/gen-test-event.py: + * NEWS: + another bug-fix for opposite hemispheres differences of psibar, where + a missing sign needed to be added in the second hemisphere. This + changes the reference files for opposite hemisphere dpsibar. For + same-hemisphere dpsibar, the sign can change (depending on the hemisphere). + Also updated python/gen-test-event.py to have command-line options + for more flexibility in examining output and for testing explanation + of psibar absolute value in RecursiveLundEEGenerator comments. + 2024-01-10 Gavin Salam + Alexander Karlberg * NEWS: * VERSION: preparing for release 2.1.0 (NEWS still needs correct date to be added) * AUTHORS: added publication info for soft spin paper 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.