1#ifndef __FASTJET_CLUSTERSEQUENCE_HH__
2#define __FASTJET_CLUSTERSEQUENCE_HH__
37#include "fastjet/PseudoJet.hh"
44#include "fastjet/Error.hh"
45#include "fastjet/JetDefinition.hh"
46#include "fastjet/SharedPtr.hh"
47#include "fastjet/LimitedWarning.hh"
48#include "fastjet/FunctionOfPseudoJet.hh"
49#include "fastjet/ClusterSequenceStructure.hh"
50#include "fastjet/internal/thread_safety_helpers.hh"
51#include "fastjet/internal/deprecated.hh"
53FASTJET_BEGIN_NAMESPACE
57class ClusterSequenceStructure;
58class DynamicNearestNeighbours;
75 const std::vector<L> & pseudojets,
77 const bool & writeout_combinations =
false);
81 transfer_from_sequence(cs);
98 std::vector<PseudoJet> inclusive_jets (
const double ptmin = 0.0)
const;
103 int n_exclusive_jets (
const double dcut)
const;
108 std::vector<PseudoJet> exclusive_jets (
const double dcut)
const;
115 std::vector<PseudoJet> exclusive_jets (
const int njets)
const;
122 std::vector<PseudoJet> exclusive_jets_up_to (
const int njets)
const;
127 double exclusive_dmerge (
const int njets)
const;
133 double exclusive_dmerge_max (
const int njets)
const;
147 int njets = n_exclusive_jets_ycut(ycut);
148 return exclusive_jets(njets);
162 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
163 const double dcut)
const;
168 int n_exclusive_subjets(
const PseudoJet & jet,
169 const double dcut)
const;
176 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
184 std::vector<PseudoJet> exclusive_subjets_up_to (
const PseudoJet & jet,
191 double exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const;
198 double exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const;
205 double Q()
const {
return _Qtot;}
207 double Q2()
const {
return _Qtot*_Qtot;}
244 std::vector<PseudoJet> constituents (
const PseudoJet & jet)
const;
256 std::ostream & ostr = std::cout)
const;
260 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
261 const std::string & filename,
262 const std::string & comment =
"")
const;
269 void add_constituents (
const PseudoJet & jet,
270 std::vector<PseudoJet> & subjet_vector)
const;
279 std::string strategy_string (
Strategy strategy_in)
const;
296 void delete_self_when_unused();
307 void signal_imminent_self_deletion()
const;
314 double jet_scale_for_algorithm(
const PseudoJet & jet)
const;
328 assert(plugin_activated());
329 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
335 void plugin_record_ij_recombination(
int jet_i,
int jet_j,
double dij,
344 assert(plugin_activated());
345 _do_iB_recombination_step(jet_i, diB);
361 assert(plugin_activated());
375 virtual std::string description()
const {
return "This is a dummy extras class that contains no extra information! Derive from it if you want to use it to provide extra information from a plugin jet finder";}
383 _extras.reset(extras_in);
392#ifdef FASTJET_HAVE_AUTO_PTR_INTERFACE
393 FASTJET_DEPRECATED_MSG(
"Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead",
394 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in)){
395 _extras.reset(extras_in.release());
414 assert(plugin_activated());
415 _simple_N2_cluster<GBJ>();
469 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
487 const std::vector<PseudoJet> & jets()
const;
500 const std::vector<history_element> & history()
const;
506 unsigned int n_particles()
const;
512 std::vector<int> particle_jet_indices(
const std::vector<PseudoJet> &)
const;
529 std::vector<int> unique_history_order()
const;
534 std::vector<PseudoJet> unclustered_particles()
const;
540 std::vector<PseudoJet> childless_pseudojets()
const;
547 bool contains(
const PseudoJet &
object)
const;
566 return _structure_shared_ptr;
579 static void print_banner();
594 static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
608#ifdef FASTJET_HAVE_LIMITED_THREAD_SAFETY
609 static std::atomic<std::ostream*> _fastjet_banner_ostr;
611 static std::ostream * _fastjet_banner_ostr;
621 template<
class L>
void _transfer_input_jets(
622 const std::vector<L> & pseudojets);
628 const bool & writeout_combinations);
634 void _initialise_and_run_no_decant();
646 const bool & writeout_combinations);
652 void _decant_options_partial();
657 void _fill_initial_history();
662 void _do_ij_recombination_step(
const int jet_i,
const int jet_j,
663 const double dij,
int & newjet_k);
667 void _do_iB_recombination_step(
const int jet_i,
const double diB);
674 void _set_structure_shared_ptr(
PseudoJet & j);
678 void _update_structure_use_count();
692 _Parabola(
double a,
double b,
double c) : _a(a), _b(b), _c(c) {}
693 inline double operator()(
const double R)
const {
return _c*(_a*R*R + _b*R + 1);}
704 _Line(
double a,
double b) : _a(a), _b(b) {}
705 inline double operator()(
const double R)
const {
return _a*R + _b;}
729 void get_subhist_set(std::set<const history_element*> & subhist,
730 const PseudoJet & jet,
double dcut,
int maxjet)
const;
732 bool _writeout_combinations;
734 double _Rparam, _R2, _invR2;
740 int _structure_use_count_after_construction;
745#ifdef FASTJET_HAVE_THREAD_SAFETY
746 mutable std::atomic<bool> _deletes_self_when_unused;
753 bool _plugin_activated;
756 void _really_dumb_cluster ();
757 void _delaunay_cluster ();
759 template<
class BJ>
void _simple_N2_cluster ();
760 void _tiled_N2_cluster ();
761 void _faster_tiled_N2_cluster ();
764 void _minheap_faster_tiled_N2_cluster();
768 void _CP2DChan_cluster();
769 void _CP2DChan_cluster_2pi2R ();
770 void _CP2DChan_cluster_2piMultD ();
771 void _CP2DChan_limited_cluster(
double D);
772 void _do_Cambridge_inclusive_jets();
775 void _fast_NsqrtN_cluster();
777 void _add_step_to_history(
779 const int parent2,
const int jetp_index,
784 void _extract_tree_children(
int pos, std::valarray<bool> &,
785 const std::valarray<int> &, std::vector<int> &)
const;
789 void _extract_tree_parents (
int pos, std::valarray<bool> &,
790 const std::valarray<int> &, std::vector<int> &)
const;
794 typedef std::pair<int,int> TwoVertices;
795 typedef std::pair<double,TwoVertices> DijEntry;
796 typedef std::multimap<double,TwoVertices> DistMap;
799 void _add_ktdistance_to_map(
const int ii,
821 double eta, phi, kt2, NN_dist;
830 double eta, phi, kt2, NN_dist;
832 int _jets_index, tile_index, diJ_posn;
837 inline void label_minheap_update_needed() {diJ_posn = 1;}
838 inline void label_minheap_update_done() {diJ_posn = 0;}
839 inline bool minheap_update_needed()
const {
return diJ_posn==1;}
847 template <
class J>
void _bj_set_jetinfo( J *
const jet,
848 const int _jets_index)
const;
852 void _bj_remove_from_tiles( TiledJet *
const jet)
const;
855 template <
class J>
double _bj_dist(
const J *
const jeta,
856 const J *
const jetb)
const;
860 template <
class J>
double _bj_diJ(
const J *
const jeta)
const;
865 template <
class J>
inline J * _bj_of_hindex(
866 const int hist_index,
867 J *
const head, J *
const tail)
870 for(res = head; res<tail; res++) {
871 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {
break;}
884 template <
class J>
void _bj_set_NN_nocross(J *
const jeta,
885 J *
const head,
const J *
const tail)
const;
891 template <
class J>
void _bj_set_NN_crosscheck(J *
const jeta,
892 J *
const head,
const J *
const tail)
const;
898 static const int n_tile_neighbours = 9;
904 Tile * begin_tiles[n_tile_neighbours];
906 Tile ** surrounding_tiles;
916 std::vector<Tile> _tiles;
917 double _tiles_eta_min, _tiles_eta_max;
918 double _tile_size_eta, _tile_size_phi;
919 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
923 inline int _tile_index (
int ieta,
int iphi)
const {
926 return (ieta-_tiles_ieta_min)*_n_tiles_phi
927 + (iphi+_n_tiles_phi) % _n_tiles_phi;
932 int _tile_index(
const double eta,
const double phi)
const;
933 void _tj_set_jetinfo ( TiledJet *
const jet,
const int _jets_index);
934 void _bj_remove_from_tiles(TiledJet *
const jet);
935 void _initialise_tiles();
936 void _print_tiles(TiledJet * briefjets )
const;
937 void _add_neighbours_to_tile_union(
const int tile_index,
938 std::vector<int> & tile_union,
int & n_near_tiles)
const;
939 void _add_untagged_neighbours_to_tile_union(
const int tile_index,
940 std::vector<int> & tile_union,
int & n_near_tiles);
956 class EEAccurateBriefJet :
public EEBriefJet { };
963 void _simple_N2_cluster_BriefJet();
965 void _simple_N2_cluster_EEBriefJet();
966 void _simple_N2_cluster_EEAccurateBriefJet();
977template<
class L>
void ClusterSequence::_transfer_input_jets(
978 const std::vector<L> & pseudojets) {
982 _jets.reserve(pseudojets.size()*2);
988 for (
unsigned int i = 0; i < pseudojets.size(); i++) {
989 _jets.push_back(pseudojets[i]);}
1013template<
class L> ClusterSequence::ClusterSequence (
1014 const std::vector<L> & pseudojets,
1016 const bool & writeout_combinations) :
1017 _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
1028 _initialise_and_run_no_decant();
1053 std::vector<PseudoJet> jets;
1062 if (jets.size() != 0) {
1074template <
class J>
inline void ClusterSequence::_bj_set_jetinfo(
1075 J *
const jetA,
const int _jets_index)
const {
1076 jetA->eta =
_jets[_jets_index].rap();
1077 jetA->phi =
_jets[_jets_index].phi_02pi();
1079 jetA->_jets_index = _jets_index;
1081 jetA->NN_dist = _R2;
1089template <
class J>
inline double ClusterSequence::_bj_dist(
1090 const J *
const jetA,
const J *
const jetB)
const {
1092#ifndef FASTJET_NEW_DELTA_PHI
1094 double dphi = std::abs(jetA->phi - jetB->phi);
1095 double deta = (jetA->eta - jetB->eta);
1096 if (dphi > pi) {dphi = twopi - dphi;}
1099 double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1100 double deta = (jetA->eta - jetB->eta);
1102 return dphi*dphi + deta*deta;
1106template <
class J>
inline double ClusterSequence::_bj_diJ(
const J *
const jet)
const {
1107 double kt2 = jet->kt2;
1108 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1109 return jet->NN_dist * kt2;
1116template <
class J>
inline void ClusterSequence::_bj_set_NN_nocross(
1117 J *
const jet, J *
const head,
const J *
const tail)
const {
1118 double NN_dist = _R2;
1121 for (J * jetB = head; jetB != jet; jetB++) {
1122 double dist = _bj_dist(jet,jetB);
1123 if (dist < NN_dist) {
1130 for (J * jetB = jet+1; jetB != tail; jetB++) {
1131 double dist = _bj_dist(jet,jetB);
1132 if (dist < NN_dist) {
1139 jet->NN_dist = NN_dist;
1144template <
class J>
inline void ClusterSequence::_bj_set_NN_crosscheck(J *
const jet,
1145 J *
const head,
const J *
const tail)
const {
1146 double NN_dist = _R2;
1148 for (J * jetB = head; jetB != tail; jetB++) {
1149 double dist = _bj_dist(jet,jetB);
1150 if (dist < NN_dist) {
1154 if (dist < jetB->NN_dist) {
1155 jetB->NN_dist = dist;
1160 jet->NN_dist = NN_dist;
1163FASTJET_END_NAMESPACE
Contains any information related to the clustering that should be directly accessible to PseudoJet.
bool _deletes_self_when_unused
if true then the CS will delete itself when the last external object referring to it disappears.
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
bool plugin_activated() const
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
ClusterSequenceStructure StructureType
the structure type associated with a jet belonging to a ClusterSequence
void plugin_simple_N2_cluster()
allows a plugin to run a templated clustering (nearest-neighbour heuristic)
const JetDefinition & jet_def() const
return a reference to the jet definition
int n_exclusive_jets_ycut(double ycut) const
the number of exclusive jets at the given ycut
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
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,...
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
std::vector< PseudoJet > inclusive_jets(const double ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
std::string strategy_string() const
return the name of the strategy used to cluster the event
void _decant_options_partial()
assuming that the jet definition, writeout_combinations and _structure_shared_ptr have been set (e....
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
bool will_delete_self_when_unused() const
return true if the object has been told to delete itself when unused
void print_jets_for_root(const std::vector< PseudoJet > &jets, std::ostream &ostr=std::cout) const
output the supplied vector of jets in a format that can be read by an appropriate root script; the fo...
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
ClusterSequence(const ClusterSequence &cs)
copy constructor for a ClusterSequence
double exclusive_ymerge(int njets) const
return the ymin corresponding to the recombination that went from n+1 to n jets (sometimes known as y...
double exclusive_ymerge_max(int njets) const
same as exclusive_dmerge_max, but normalised to squared total energy
void plugin_associate_extras(Extras *extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
double Q() const
returns the sum of all energies in the event (relevant mainly for e+e-)
ClusterSequence()
default constructor
std::vector< PseudoJet > exclusive_jets_ycut(double ycut) const
the exclusive jets obtained at the given ycut
double Q2() const
return Q()^2
PseudoJet & plugin_non_const_jet(unsigned i)
return a non-const reference to the jets()[i], to allow plugins to modify its contents.
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
void _transfer_input_jets(const std::vector< L > &pseudojets)
transfer the vector<L> of input jets into our own vector<PseudoJet> _jets (with some reserved space f...
double jet_scale_for_algorithm(const PseudoJet &jet) const
returns the scale associated with a jet as required for this clustering algorithm (kt^2 for the kt-al...
const Extras * extras() const
returns a pointer to the extras object (may be null)
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],...
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
base class providing interface for a generic function of a PseudoJet
class that is intended to hold a full definition of the jet clusterer
bool is_spherical() const
returns true if the jet definition involves an algorithm intended for use on a spherical geometry (e....
std::vector< PseudoJet > operator()(const std::vector< L > &particles) const
cluster the supplied particles and returns a vector of resulting jets, sorted by pt (or energy in the...
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
structure analogous to BriefJet, but with the extra information needed for dealing with tiles
provides an object wich will return "true" the first time () is called and false afterwards
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
JetAlgorithm
the various families of jet-clustering algorithm
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
a single element in the clustering history
int parent1
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
int jetp_index
index in the _jets vector where we will find the PseudoJet object corresponding to this jet (i....
int parent2
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
int child
index in _history where the current jet is
double max_dij_so_far
the largest recombination distance seen so far in the clustering history.
double dij
the distance corresponding to the recombination at this stage of the clustering.