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.