Page MenuHomeHEPForge
Files F8310305

No OneTemporary
Actions

View Options

Index: contrib/contribs/CMPPlugin/trunk/example-CMP.cc
===================================================================
--- contrib/contribs/CMPPlugin/trunk/example-CMP.cc (revision 1469)
+++ contrib/contribs/CMPPlugin/trunk/example-CMP.cc (revision 1470)
@@ -1,174 +1,176 @@
// To run this example, use the following command:
//
// ./example-CMP < ../data/pythia8_Zq_vshort.dat
//
// NB: the example file reads in a file with 6 light flavours
// (though for this algorithm we need to strip off all but one), and
// an extra high density of quarks, to help with "make check"
// test things more thoroughly
//----------------------------------------------------------------------
// $Id$
//
-// Copyright (c) 2023, Fabrizio Caola, Radoslaw Grabarczyk,
+// Copyright (c) 2024, Michal Czakon, Alexander Mitov, and Rene Poncelet
+//
+// based on initial version by Fabrizio Caola, Radoslaw Grabarczyk,
// Maxwell Hutt, Gavin P. Salam, Ludovic Scyboz, and Jesse Thaler
//
//----------------------------------------------------------------------
// 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 <iomanip>
#include "fastjet/PseudoJet.hh"
#include "fastjet/contrib/CMPPlugin.hh"
using namespace std;
using namespace fastjet;
using namespace fastjet::contrib;
// forward declaration to make things clearer
void read_event(vector<PseudoJet> &event);
//----------------------------------------------------------------------
int main(int iargc, char **argv){
// give user control over printout (mainly relevant for make check)
// usage: "./example [nevmax [njetmax]] < data/pythia8_Zq_vshort.dat"
unsigned int nevmax = 2;
unsigned int njetmax = 1;
if (iargc > 1) nevmax = stoi(argv[1]);
if (iargc > 2) njetmax = stoi(argv[2]);
// print banner for FastJet at the start, so it doesn't mix
// into the other output
ClusterSequence::print_banner();
// we start with a base jet definition (here antikt)
double R = 0.4;
JetDefinition base_jet_def(antikt_algorithm, R);
// enable it to track flavours (default is net flavour)
FlavRecombiner flav_recombiner;
base_jet_def.set_recombiner(&flav_recombiner);
// CMP parameters:
// CMP 'a' parameter in
// kappa_ij = 1/a * (kT_i^2 + kT_j^2) / (2*kT_max^2)
double CMP_a = 0.1;
// correction to original CMP algo: do not change this if you want IRC safety!
CMPPlugin::CorrectionType CMP_corr =
CMPPlugin::CorrectionType::SqrtCoshyCosPhiArgument_a2;
// Dynamic definition of ktmax
CMPPlugin::ClusteringType CMP_clust = CMPPlugin::ClusteringType::DynamicKtMax;
// CMP plugin
JetDefinition CMP_jet_def(new CMPPlugin(R, CMP_a, CMP_corr, CMP_clust));
// enable it to track flavours (default is net flavour)
CMP_jet_def.set_recombiner(&flav_recombiner);
CMP_jet_def.delete_plugin_when_unused();
cout << "! this analysis considers only b-flavour from input events !" << endl;
cout << "base jet definition: " << base_jet_def.description() << endl;
cout << "CMP jet definition: " << CMP_jet_def.description() << endl;
// loop over some number of events
unsigned int n_events = 10;
for (unsigned int iev = 0; iev < n_events && iev < nevmax; iev++) {
// read in input particles: see that routine for info
// on how to set up the PseudoJets with flavour information
vector<PseudoJet> event;
read_event(event);
cout << "\n#---------------------------------------------------------------\n";
cout << "# read event " << iev << " with " << event.size() << " particles" << endl;
// run the jet clustering with the base jet definition and the
// CMPPlugin-based jet definition
vector<PseudoJet> base_jets = base_jet_def(event);
vector<PseudoJet> CMP_jets = CMP_jet_def(event);
// ----------------------------------------------------
// loop over the two leading jets and print out their properties
for (unsigned int ijet = 0; ijet < base_jets.size() && ijet < njetmax; ijet++) {
// first print out the original anti-kt jets
const auto & base_jet = base_jets[ijet];
cout << endl;
cout << "base jet " << ijet << ": ";
cout << "pt=" << base_jet.pt() << " rap=" << base_jet.rap() << " phi=" << base_jet.phi();
cout << ", flav = " << FlavHistory::current_flavour_of(base_jet).description() << endl;
// for the first event, print out the jet constituents' pt and initial and final flavours
cout << "constituents:" << endl;
for (const auto & c: sorted_by_pt(base_jet.constituents())) {
cout << " pt = " << setw(10) << c.pt();
cout << ", orig. flav = " << setw(8) << FlavHistory::initial_flavour_of(c).description();
cout << ", final flav = " << setw(8) << FlavHistory::current_flavour_of(c).description();
cout << endl;
}
}
// ----------------------------------------------------
// loop over the two leading jets and print out their properties
for (unsigned int ijet = 0; ijet < CMP_jets.size() && ijet < njetmax; ijet++) {
// and the CMP jets
const auto & CMP_jet = CMP_jets [ijet];
cout << endl;
cout << "CMP jet " << ijet << ": ";
cout << "pt=" << CMP_jet.pt() << " rap=" << CMP_jet.rap() << " phi=" << CMP_jet.phi();
cout << ", flav = " << FlavHistory::current_flavour_of(CMP_jet).description() << endl;
// for the first event, print out the jet constituents' pt and initial and final flavours
cout << "constituents:" << endl;
for (const auto & c: sorted_by_pt(CMP_jet.constituents())) {
cout << " pt = " << setw(10) << c.pt();
cout << ", orig. flav = " << setw(8) << FlavHistory::initial_flavour_of(c).description();
cout << ", final flav = " << setw(8) << FlavHistory::current_flavour_of(c).description();
cout << endl;
}
}
}
return 0;
}
// read in input particles and set up PseudoJets with flavour information
void read_event(vector<PseudoJet> &event){
// read in the input particles and their PDG IDs
string line;
double px, py, pz, E;
int pdg_id;
event.resize(0);
while(getline(cin,line)) {
if(line[0] == '#') continue;
istringstream iss(line);
iss >> px >> py >> pz >> E >> pdg_id;
// create a fastjet::PseudoJet with these components and put it onto
// back of the input_particles vector
PseudoJet p(px,py,pz,E);
// assign information about flavour (will be deleted automatically)
FlavInfo flav_info_init(pdg_id);
// we take only b-quarks
flav_info_init.reset_all_but_flav(5);
// need to cast to const for constructor
p.set_user_info(new FlavHistory(const_cast<const FlavInfo&>(flav_info_init)));
event.push_back(p);
if (cin.peek() == '\n' || cin.peek() == EOF) {
getline(cin,line);
break;
}
}
}
Index: contrib/contribs/CMPPlugin/trunk/README
===================================================================
--- contrib/contribs/CMPPlugin/trunk/README (revision 1469)
+++ contrib/contribs/CMPPlugin/trunk/README (revision 1470)
@@ -1,66 +1,52 @@
-CMPPlugin Contrib
-===============
-
-This contrib provides a range of tools for studying Jet Flavour,
-provided jointly by several groups who are active on the subject (see
-[AUTHORS](AUTHORS)).
-
-It includes
-
-- Flavoured anti-kT algorithm by Czakon-Mitov-Poncelet ([CMPPlugin](#cmpplugin))
-
-**The CMPPlugin code needs FastJet ≥ 3.4.1 to compile**
-
CMPPlugin
=========
-/// Implementation of the
-
+/// Implementation of the CMP flavoured anti-kt algorithm
This code provides an implementation of the flavoured anti-kt algorithm
from
> Infrared-safe flavoured anti-kT jets,
> by Michal Czakon, Alexander Mitov and Rene Poncelet
> https://arxiv.org/pdf/2205.11879 (v2)
including modifications to ensure IR-safety suggested in
> Flavoured jets with exact anti-kt kinematics and tests of infrared and collinear safety,
> by Fabrizio Caola, Radoslaw Grabarczyk, Maxwell Hutt, Gavin P. Salam, Ludovic Scyboz, and Jesse Thaler
> https://arxiv.org/abs/2306.07314
To learn how to use the library, the [`example-CMP.cc`](example-CMP.cc) code is
a good place to get started.
Main principles of CMP
----------------------
The CMP algorithm provides a anti-kT-like flavour safe jet clustering algorithm.
IR safety is ensured through a modified pseudojet distance for "soft"
flavour-anti-flavour pairs. Through this modification the clustering sequence is
not exactly anti-kT. The modification can be steered by a parameter a (the limit
a->0 recovers the IR flavour unsafe anti-kT algorithm).
Code Structure
--------------
The main interface to the CMP algorithm is in [`CMPPlugin.hh`](CMPPlugin.hh):
```cpp
// create an CMPPlugin for a given a and jetradius R
double a = 0.1;
double R = 0.4;
JetDefinition jet_def(new CMPPlugin(R,a));
jet_def.delete_plugin_when_unused();
```
The plugin can be used as standard with the FastJet package, e.g.:
```cpp
vector<PseudoJet> particles;
// ... fill the particles ...
auto jets = jet_def(particles);
for (const auto & jet : jets) {
cout << "jet pt = " << jet.pt()
<< ", flav = " << FlavHistory::current_flavour_of(jet) << endl;
}
```
Index: contrib/contribs/CMPPlugin/trunk/include/fastjet/contrib/CMPPlugin.hh
===================================================================
--- contrib/contribs/CMPPlugin/trunk/include/fastjet/contrib/CMPPlugin.hh (revision 1469)
+++ contrib/contribs/CMPPlugin/trunk/include/fastjet/contrib/CMPPlugin.hh (revision 1470)
@@ -1,168 +1,171 @@
-
/// Implementation of the flavoured anti-kt algorithm
/// by Michal Czakon, Alexander Mitov and Rene Poncelet,
/// as described in https://arxiv.org/pdf/2205.11879 (v2)
///
+/// Authors: Michal Czakon, Alexander Mitov, and Rene Poncelet
+///
+/// based on initial version by
+///
/// Authors: Fabrizio Caola, Radoslaw Grabarczyk, Maxwell Hutt,
-// Rene Poncelet, Gavin P. Salam, Ludovic Scyboz,
-// and Jesse Thaler
+/// Gavin P. Salam, Ludovic Scyboz, and Jesse Thaler
+///
#ifndef __CMPPLUGIN_HH__
#define __CMPPLUGIN_HH__
#include "fastjet/contrib/FlavInfo.hh"
// to facilitate use with fjcore
#ifndef __FJC_FLAVINFO_USEFJCORE__
#include "fastjet/NNH.hh"
#endif
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
using namespace std;
using namespace contrib;
/// Plugin for algorithm proposed by Czakon, Mitov and Poncelet
/// as described in arXiv:2205.11879 with the modification
class CMPPlugin : public JetDefinition::Plugin {
public:
/// correction to original CMP
enum CorrectionType {
/// this uses the CMP distance measure factor as written in 2205.11879v1
/// Eq. 2.9:
///
/// Sij = 1 - Theta(1-κ)cos(πκ/2)
NoCorrection,
/// SqrtCoshyCosPhiArgument causes kappa in Eq.(2.9) to be
/// mutiplied by
///
/// sqrt(M), with M = 2*(cosh(drap) - cos(dphi))/deltaR2
///
/// which tends to 1 in the small deltaR^2 limit, and fixes
/// a joint infrared/collinear safety issue when deltaR is large
///
SqrtCoshyCosPhiArgument,
/// Default: same as above, but with damping a = 2
///
/// sqrt(M), with M = 2*(1/a^2*[cosh(a*drap)-1] - [cos(dphi)-1])/deltaR2
///
SqrtCoshyCosPhiArgument_a2,
// NOT RECOMMENDED: with M from above, uses
//
// Sij = 1 - Theta(1-κ)cos(πκ/2) * M
CoshyCosPhi,
// with M from above, uses
//
// Sij = (1 - Theta(1-κ)cos(πκ/2)) * M
OverAllCoshyCosPhi,
// same as above, but with damping a = 2
OverAllCoshyCosPhi_a2
};
/// definition of ktmax
/// 0 : dynamic ktmax (pT of the hardest pseudojet lying around, default)
/// 1 : fixed ktmax (pT of the hardest pseudojet clustered with anti-kt)
enum ClusteringType {DynamicKtMax, FixedKtMax};
/// Main constructor for the class
///
/// @param R: the jet radius parameter
///
/// @param a: CMP parameter a
///
/// @param correction_type: the type of correction we use in the CMP distance
///
/// @param clustering_type: the definition of ktmax in the CMP distance
///
/// @param spherical: true if using the e+e- version of the algorithm
CMPPlugin(double R, double a,
CorrectionType correction_type = SqrtCoshyCosPhiArgument_a2,
ClusteringType clustering_type = DynamicKtMax, bool spherical=false)
: _R(R), _a(a), _correction_type(correction_type),
_clustering_type(clustering_type), _spherical(spherical) {}
/// copy constructor
CMPPlugin(const CMPPlugin &plugin) { *this = plugin; }
/// CMP parameter
double a() const { return _a; }
/// Jet radius of antikt distance
double R() const { return _R; }
/// whether to use the e+e- version of the algo
bool is_spherical() const { return _spherical; }
// Required by base class:
virtual std::string description() const;
virtual double precise_squared_distance(const PseudoJet & j1, const PseudoJet & j2) const;
virtual double distance_opposite_flavour(const PseudoJet & j1, const PseudoJet & j2,
const double ktmax) const;
virtual void run_clustering(ClusterSequence &) const;
void cross_product(const PseudoJet & p1, const PseudoJet & p2,
PseudoJet & retp, bool lightlike=false) const {
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;
}
// is this exactly what we want? The masses seem to change, need to check!
retp.reset(px, py, pz, E);
return;
}
double dot_product_3d(const PseudoJet & a, const PseudoJet & b) const {
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
double one_minus_costheta(const PseudoJet & p1, const PseudoJet & p2) const {
if (p1.m2() == 0 && p2.m2() == 0) {
// use the 4-vector dot product.
// For massless particles it gives us E1*E2*(1-cos theta)
double res = dot_product(p1,p2) / (p1.E() * p2.E());
return res;
} else {
double p1mod = sqrt(p1.modp2());
double p2mod = sqrt(p2.modp2());
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, cross_result, 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
double res = -cross_result.m2()/(p1p2mod * (p1p2mod+dot));
return res;
}
return 1.0 - dot/p1p2mod;
}
}
private:
static const double _deltaR2_handover;
double _R, _a;
CorrectionType _correction_type;
ClusteringType _clustering_type;
bool _spherical;
template <typename NN>
void _NN_clustering(ClusterSequence &cs, NN &nn) const;
//template <typename NN, typename CMPNNInfo>
//void k_tmax_clustering(ClusterSequence &cs, NN &nn, CMPNNInfo info) const;
};
FASTJET_END_NAMESPACE
#endif // __CMPPLUGIN_HH__
Index: contrib/contribs/CMPPlugin/trunk/AUTHORS
===================================================================
--- contrib/contribs/CMPPlugin/trunk/AUTHORS (revision 1469)
+++ contrib/contribs/CMPPlugin/trunk/AUTHORS (revision 1470)
@@ -1,4 +1,4 @@
-The flavoured anti-kT algorithm by Czakon-Mitov-Poncelet
-(CMPPlugin) code has been provided by Fabrizio Caola, Radoslaw
-Grabarczyk, Maxwell Hutt, Rene Poncelet, Gavin P. Salam,
-Ludovic Scyboz, and Jesse Thaler
+The CMP flavoured anti-kT algorithm by Michal Czakon, Alexander Mitov,
+and Rene Poncelet (CMPPlugin) based on initial version by Fabrizio
+Caola, Radoslaw Grabarczyk, Maxwell Hutt, Gavin P. Salam, Ludovic Scyboz,
+and Jesse Thaler
Index: contrib/contribs/CMPPlugin/trunk/ChangeLog
===================================================================
--- contrib/contribs/CMPPlugin/trunk/ChangeLog (revision 1469)
+++ contrib/contribs/CMPPlugin/trunk/ChangeLog (revision 1470)
@@ -1,12 +1,26 @@
+2024-12-12 Rene Poncelet
+
+ * AUTHORS:
+ corrected attributions
+
+ * README:
+ Reflect corrected structure of a single plugin
+
+ * include/fastjet/contrib/CMPPlugin.hh:
+ Fixed authorship in headers
+
+ * example-CMP.cc:
+ Fixed authorship in headers
+
2024-12-11 Gavin Salam
* example-CMP.cc:
changed example line to reflect actual usage
* include/fastjet/contrib/CMPPlugin.hh:
added fastjet/contrib to #include "FlavInfo.hh"
2024-12-11 Rene Poncelet:
* [many files]:
Initial creation and import of CMPPlugin version 1.0.0 from
git@github.com:jetflav/CMPPlugin.git (7673166d9dc)

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:02 PM (7 h, 11 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4015537
Default Alt Text
(17 KB)

Event Timeline