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 <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#ifndef __FASTJET_CONTRIB_RECURSIVELUNDEEGENERATOR_HH__
#define __FASTJET_CONTRIB_RECURSIVELUNDEEGENERATOR_HH__
#include "LundEEHelpers.hh"
#include <fastjet/internal/base.hh>
#include "fastjet/tools/Recluster.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/PseudoJet.hh"
#include <string>
#include <vector>
#include <utility>
#include <queue>
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<double,double> const lund_coordinates() const {
return std::pair<double,double>(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<LundEEDeclustering> result(const ClusterSequence & cs) const {
std::vector<PseudoJet> 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<LundEEDeclustering> 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<LundEEDeclustering> & 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__
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 <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#include <iostream>
#include <fstream>
#include <sstream>
#include "fastjet/PseudoJet.hh"
#include "fastjet/EECambridgePlugin.hh"
#include <string>
#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<PseudoJet> &event);
//----------------------------------------------------------------------
int main(){
//----------------------------------------------------------
// read in input particles
vector<PseudoJet> 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<contrib::LundEEDeclustering> 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<declusts.size(); ++i) {
if (declusts[i].depth() == 0 && declusts[i].z() > 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<const contrib::LundEEDeclustering *> 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; i<secondaries.size(); ++i) {
- if (secondaries[i]->z() > 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; i<secondaries.size(); ++i) {
+ if (secondaries[i]->z() > 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<double,double> 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<double,double> 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<PseudoJet> &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 <gavin.salam@physics.ox.ac.uk>
* example_dpsi_slice.cc:
replaced a uint -> unsigned int, to resolve a compilation issue
with gcc-13 on MacOS
2023-09-22 Fri <gavin.salam@physics.ox.ac.uk>
* 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 <soyez@fastjet.fr>
* 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 <gavin.salam@physics.ox.ac.uk> + 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 <gavin.salam@physics.ox.ac.uk>
* VERSION:
* NEWS:
prepared for release 2.0.2
* Makefile:
updated some dependencies
* EEHelpers.hh:
added #include <limits>, as per request from Andy Buckley
for compilation with g++-12 on some systems
2021-12-06 Gavin Salam <gavin.salam@physics.ox.ac.uk>
* NEWS:
* VERSION:
prepared for release 2.0.1
* AUTHORS:
fixed missing names and publication info
2021-12-06 Gregory Soyez <soyez@fastjet.fr>
* SecondaryLund.cc:
fixed int v. unsigned int in loop over vector indices
2021-12-06 Gavin Salam <gavin.salam@physics.ox.ac.uk>
* example.py:
fixed name of executable in comments about how to execute
this (thanks to Matteo Cacciari)
2021-11-09 Ludovic Scyboz <ludovic.scyboz@physics.ox.ac.uk>
* 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 <gavin.salam@cern.ch>
* 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 <gavin.salam@cern.ch>
* read_lund_json.py:
removed extraneous normalisation of zeroth bin in
the LundImage class.
Added documentation.
2018-08-30 Gavin Salam <gavin.salam@cern.ch>
* VERSION:
* NEWS:
Release of version 1.0.1
2018-08-23 Gavin Salam <gavin.salam@cern.ch>
* 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 <fdreyer@mit.edu>
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 <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#ifndef __FASTJET_CONTRIB_EEHELPERS_HH__
#define __FASTJET_CONTRIB_EEHELPERS_HH__
#include "fastjet/PseudoJet.hh"
#include <array>
#include <limits>
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<double>::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<double>::epsilon()) && fabs(n2.modp()-1) < sqrt(std::numeric_limits<double>::epsilon()));
double omcost = one_minus_costheta(n1,n2);
double theta;
// If theta ~ pi, we return pi.
if(fabs(omcost-2) < sqrt(std::numeric_limits<double>::epsilon())) {
theta = M_PI;
} else if (omcost > sqrt(std::numeric_limits<double>::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<array<...,3>,3> object
Matrix3(const std::array<std::array<double,3>,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<class T>
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<class T>
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<std::array<double,3>,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<double,3> 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