Loading...
Searching...
No Matches
JadePlugin.cc
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2007-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet 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.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet 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 FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
// fastjet stuff
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JadePlugin.hh"
#include <iostream>
//#include "fastjet/internal/ClusterSequence_N2.icc"
#include "fastjet/NNH.hh"
#include "fastjet/NNFJN2Plain.hh"
// other stuff
#include <vector>
#include <sstream>
#include <limits>
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
//----------------------------------------------------------------------
/// class to help run a JADE algorithm
///
/// This class works both with NNH and NNFJN2Plain clustering
/// helpers. They both use the same init(...) call, but for the
/// clustering:
///
/// - NNH uses distance(...) and beam_distance()
/// - NNFJPlainN2 uses geometrical_distance(...), momentum_factor()
/// and geometrical_beam_distance()
///
/// For NNFJPlainN2 the 2 E_i E_j (1-cos theta_{ij}) factor
/// gets broken up into
///
/// sqrt(2)*min(E_i,E_j) * [sqrt(2)*max(E_i,E_j) (1 - cos \theta_{ij})]
///
/// The second factor is what we call the "geometrical_distance" even
/// though it isn't actually purely geometrical. But the fact that it
/// gets multiplied by min(E_i,E_j) to get the full distance is
/// sufficient for the validity of the FJ lemma, allowing for the use
/// of NNFJN2Plain.
class JadeBriefJet {
public:
void init(const PseudoJet & jet) {
double norm = 1.0/sqrt(jet.modp2());
nx = jet.px() * norm;
ny = jet.py() * norm;
nz = jet.pz() * norm;
rt2E = sqrt(2.0)*jet.E();
}
double distance(const JadeBriefJet * jet) const {
double dij = 1 - nx*jet->nx
- ny*jet->ny
- nz*jet->nz;
dij *= rt2E*jet->rt2E;
return dij;
}
double geometrical_distance(const JadeBriefJet * jet) const {
double dij = 1 - nx*jet->nx
- ny*jet->ny
- nz*jet->nz;
dij *= max(rt2E,jet->rt2E);
return dij;
}
double momentum_factor() const {
return rt2E;
}
double beam_distance() const {
return numeric_limits<double>::max();
}
double geometrical_beam_distance() const {
// get a number that is almost the same as max(), just a little
// smaller so as to ensure that when we divide it by rt2E and then
// multiply it again, we won't get an overflow.
// Watch out for cases where rt2E < 1.0 (cf. bug fix from
// andrii.verbytskyi@mpp.mpg.de on 2019-02-14)
const double almost_max = numeric_limits<double>::max() * (1 - 1e-13);
if (rt2E>1.0) return almost_max / rt2E;
else return almost_max;
}
private:
double rt2E, nx, ny, nz;
};
//----------------------------------------------------------------------
string JadePlugin::description () const {
ostringstream desc;
desc << "e+e- JADE algorithm plugin";
switch(_strategy) {
case strategy_NNH:
desc << ", using NNH strategy"; break;
case strategy_NNFJN2Plain:
desc << ", using NNFJN2Plain strategy"; break;
default:
throw Error("Unrecognized strategy in JadePlugin");
}
return desc.str();
}
// //----------------------------------------------------------------------
// void JadePlugin::run_clustering(ClusterSequence & cs) const {
// int njets = cs.jets().size();
//
// //SharedPtr<NNBase<> > nn;
// NNBase<> * nn;
// switch(_strategy) {
// case strategy_NNH:
// //nn.reset(new NNH<JadeBriefJet>(cs.jets()));
// nn = new NNH<JadeBriefJet>(cs.jets());
// break;
// case strategy_NNFJN2Plain:
// //nn.reset(new NNFJN2Plain<JadeBriefJet>(cs.jets()));
// nn = new NNFJN2Plain<JadeBriefJet>(cs.jets());
// break;
// default:
// throw Error("Unrecognized strategy in JadePlugin");
// }
// //NNH<JadeBriefJet> nnh(cs.jets());
// //NNFJN2Plain<JadeBriefJet> nnh(cs.jets());
//
// // if testing against Hoeth's implementation, need to rescale the
// // dij by Q^2.
// //double Q2 = cs.Q2();
//
// while (njets > 0) {
// int i, j, k;
// double dij = nn->dij_min(i, j);
//
// if (j >= 0) {
// cs.plugin_record_ij_recombination(i, j, dij, k);
// nn->merge_jets(i, j, cs.jets()[k], k);
// } else {
// double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
// cs.plugin_record_iB_recombination(i, diB);
// nn->remove_jet(i);
// }
// njets--;
// }
// delete nn;
// }
template<class N> void JadePlugin::_actual_run_clustering(ClusterSequence & cs) const {
int njets = cs.jets().size();
N nn(cs.jets());
// if testing against Hoeth's implementation, need to rescale the
// dij by Q^2.
//double Q2 = cs.Q2();
while (njets > 0) {
int i, j, k;
double dij = nn.dij_min(i, j);
if (j >= 0) {
cs.plugin_record_ij_recombination(i, j, dij, k);
nn.merge_jets(i, j, cs.jets()[k], k);
} else {
double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
cs.plugin_record_iB_recombination(i, diB);
nn.remove_jet(i);
}
njets--;
}
}
//----------------------------------------------------------------------
void JadePlugin::run_clustering(ClusterSequence & cs) const {
switch(_strategy) {
case strategy_NNH:
_actual_run_clustering<NNH<JadeBriefJet> >(cs);
break;
case strategy_NNFJN2Plain:
_actual_run_clustering<NNFJN2Plain<JadeBriefJet> >(cs);
break;
default:
throw Error("Unrecognized strategy in JadePlugin");
}
}
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
fastjet::ClusterSequence
deals with clustering
Definition ClusterSequence.hh:63
fastjet::ClusterSequence::plugin_record_iB_recombination
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam,...
Definition ClusterSequence.hh:343
fastjet::ClusterSequence::jets
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
Definition ClusterSequence.hh:1003
fastjet::ClusterSequence::plugin_record_ij_recombination
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k],...
Definition ClusterSequence.hh:326
fastjet::Error
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
fastjet::JadePlugin::description
virtual std::string description() const
return a textual description of the jet-definition implemented in this plugin
Definition JadePlugin.cc:124
fastjet::JadePlugin::run_clustering
virtual void run_clustering(ClusterSequence &) const
given a ClusterSequence that has been filled up with initial particles, the following function should...
Definition JadePlugin.cc:210
fastjet::norm
double norm(const VPoint p)
norm of a vector
Definition Voronoi.hh:155