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" 51 #include "fastjet/internal/deprecated.hh" 53 FASTJET_BEGIN_NAMESPACE
57 class ClusterSequenceStructure;
58 class 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;
255 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
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();
303 void signal_imminent_self_deletion()
const;
309 double jet_scale_for_algorithm(
const PseudoJet & jet)
const;
323 assert(plugin_activated());
324 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
330 void plugin_record_ij_recombination(
int jet_i,
int jet_j,
double dij,
339 assert(plugin_activated());
340 _do_iB_recombination_step(jet_i, diB);
353 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";}
361 _extras.reset(extras_in);
369 #ifdef FASTJET_HAVE_AUTO_PTR_INTERFACE 370 FASTJET_DEPRECATED_MSG(
"Please use ClusterSequence::plugin_associate_extras(Extras * extras_in)) instead")
371 inline
void plugin_associate_extras(
std::auto_ptr<
Extras> extras_in){
372 _extras.reset(extras_in.release());
390 assert(plugin_activated());
391 _simple_N2_cluster<GBJ>();
439 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
457 const std::vector<PseudoJet> & jets()
const;
470 const std::vector<history_element> & history()
const;
476 unsigned int n_particles()
const;
482 std::vector<int> particle_jet_indices(
const std::vector<PseudoJet> &)
const;
499 std::vector<int> unique_history_order()
const;
504 std::vector<PseudoJet> unclustered_particles()
const;
510 std::vector<PseudoJet> childless_pseudojets()
const;
517 bool contains(
const PseudoJet &
object)
const;
536 return _structure_shared_ptr;
549 static void print_banner();
564 static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
578 static std::ostream * _fastjet_banner_ostr;
588 template<
class L>
void _transfer_input_jets(
589 const std::vector<L> & pseudojets);
595 const bool & writeout_combinations);
601 void _initialise_and_run_no_decant();
613 const bool & writeout_combinations);
619 void _decant_options_partial();
624 void _fill_initial_history();
629 void _do_ij_recombination_step(
const int jet_i,
const int jet_j,
630 const double dij,
int & newjet_k);
634 void _do_iB_recombination_step(
const int jet_i,
const double diB);
641 void _set_structure_shared_ptr(
PseudoJet & j);
645 void _update_structure_use_count();
659 _Parabola(
double a,
double b,
double c) : _a(a), _b(b), _c(c) {}
660 inline double operator()(
const double R)
const {
return _c*(_a*R*R + _b*R + 1);}
671 _Line(
double a,
double b) : _a(a), _b(b) {}
672 inline double operator()(
const double R)
const {
return _a*R + _b;}
696 void get_subhist_set(std::set<const history_element*> & subhist,
697 const PseudoJet & jet,
double dcut,
int maxjet)
const;
699 bool _writeout_combinations;
701 double _Rparam, _R2, _invR2;
707 int _structure_use_count_after_construction;
716 bool _plugin_activated;
719 void _really_dumb_cluster ();
720 void _delaunay_cluster ();
722 template<
class BJ>
void _simple_N2_cluster ();
723 void _tiled_N2_cluster ();
724 void _faster_tiled_N2_cluster ();
727 void _minheap_faster_tiled_N2_cluster();
731 void _CP2DChan_cluster();
732 void _CP2DChan_cluster_2pi2R ();
733 void _CP2DChan_cluster_2piMultD ();
734 void _CP2DChan_limited_cluster(
double D);
735 void _do_Cambridge_inclusive_jets();
738 void _fast_NsqrtN_cluster();
740 void _add_step_to_history(
742 const int parent2,
const int jetp_index,
747 void _extract_tree_children(
int pos, std::valarray<bool> &,
748 const std::valarray<int> &, std::vector<int> &)
const;
752 void _extract_tree_parents (
int pos, std::valarray<bool> &,
753 const std::valarray<int> &, std::vector<int> &)
const;
757 typedef std::pair<int,int> TwoVertices;
758 typedef std::pair<double,TwoVertices> DijEntry;
759 typedef std::multimap<double,TwoVertices> DistMap;
762 void _add_ktdistance_to_map(
const int ii,
764 const DynamicNearestNeighbours * DNN);
768 static bool _first_time;
784 double eta, phi, kt2, NN_dist;
793 double eta, phi, kt2, NN_dist;
794 TiledJet * NN, *previous, * next;
795 int _jets_index, tile_index, diJ_posn;
800 inline void label_minheap_update_needed() {diJ_posn = 1;}
801 inline void label_minheap_update_done() {diJ_posn = 0;}
802 inline bool minheap_update_needed()
const {
return diJ_posn==1;}
810 template <
class J>
void _bj_set_jetinfo( J *
const jet,
811 const int _jets_index)
const;
815 void _bj_remove_from_tiles( TiledJet *
const jet)
const;
818 template <
class J>
double _bj_dist(
const J *
const jeta,
819 const J *
const jetb)
const;
823 template <
class J>
double _bj_diJ(
const J *
const jeta)
const;
828 template <
class J>
inline J * _bj_of_hindex(
829 const int hist_index,
830 J *
const head, J *
const tail)
833 for(res = head; res<tail; res++) {
834 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {
break;}
847 template <
class J>
void _bj_set_NN_nocross(J *
const jeta,
848 J *
const head,
const J *
const tail)
const;
854 template <
class J>
void _bj_set_NN_crosscheck(J *
const jeta,
855 J *
const head,
const J *
const tail)
const;
861 static const int n_tile_neighbours = 9;
867 Tile * begin_tiles[n_tile_neighbours];
869 Tile ** surrounding_tiles;
879 std::vector<Tile> _tiles;
880 double _tiles_eta_min, _tiles_eta_max;
881 double _tile_size_eta, _tile_size_phi;
882 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
886 inline int _tile_index (
int ieta,
int iphi)
const {
889 return (ieta-_tiles_ieta_min)*_n_tiles_phi
890 + (iphi+_n_tiles_phi) % _n_tiles_phi;
895 int _tile_index(
const double eta,
const double phi)
const;
896 void _tj_set_jetinfo ( TiledJet *
const jet,
const int _jets_index);
897 void _bj_remove_from_tiles(TiledJet *
const jet);
898 void _initialise_tiles();
899 void _print_tiles(TiledJet * briefjets )
const;
900 void _add_neighbours_to_tile_union(
const int tile_index,
901 std::vector<int> & tile_union,
int & n_near_tiles)
const;
902 void _add_untagged_neighbours_to_tile_union(
const int tile_index,
903 std::vector<int> & tile_union,
int & n_near_tiles);
921 void _simple_N2_cluster_BriefJet();
923 void _simple_N2_cluster_EEBriefJet();
934 template<
class L>
void ClusterSequence::_transfer_input_jets(
935 const std::vector<L> & pseudojets) {
939 _jets.reserve(pseudojets.size()*2);
945 for (
unsigned int i = 0; i < pseudojets.size(); i++) {
946 _jets.push_back(pseudojets[i]);}
970 template<
class L> ClusterSequence::ClusterSequence (
971 const std::vector<L> & pseudojets,
973 const bool & writeout_combinations) :
974 _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
985 _initialise_and_run_no_decant();
1004 std::vector<PseudoJet> JetDefinition::operator()(
const std::vector<L> & particles)
const {
1010 std::vector<PseudoJet>
jets;
1011 if (is_spherical()) {
1019 if (jets.size() != 0) {
1031 template <
class J>
inline void ClusterSequence::_bj_set_jetinfo(
1032 J *
const jetA,
const int _jets_index)
const {
1033 jetA->eta = _jets[_jets_index].rap();
1034 jetA->phi = _jets[_jets_index].phi_02pi();
1036 jetA->_jets_index = _jets_index;
1038 jetA->NN_dist = _R2;
1046 template <
class J>
inline double ClusterSequence::_bj_dist(
1047 const J *
const jetA,
const J *
const jetB)
const {
1049 #ifndef FASTJET_NEW_DELTA_PHI 1051 double dphi = std::abs(jetA->phi - jetB->phi);
1052 double deta = (jetA->eta - jetB->eta);
1053 if (dphi > pi) {dphi = twopi - dphi;}
1056 double dphi = pi-std::abs(pi-std::abs(jetA->phi - jetB->phi));
1057 double deta = (jetA->eta - jetB->eta);
1059 return dphi*dphi + deta*deta;
1063 template <
class J>
inline double ClusterSequence::_bj_diJ(
const J *
const jet)
const {
1064 double kt2 = jet->kt2;
1065 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
1066 return jet->NN_dist * kt2;
1073 template <
class J>
inline void ClusterSequence::_bj_set_NN_nocross(
1074 J *
const jet, J *
const head,
const J *
const tail)
const {
1075 double NN_dist = _R2;
1078 for (J * jetB = head; jetB != jet; jetB++) {
1079 double dist = _bj_dist(jet,jetB);
1080 if (dist < NN_dist) {
1087 for (J * jetB = jet+1; jetB != tail; jetB++) {
1088 double dist = _bj_dist(jet,jetB);
1089 if (dist < NN_dist) {
1096 jet->NN_dist = NN_dist;
1101 template <
class J>
inline void ClusterSequence::_bj_set_NN_crosscheck(J *
const jet,
1102 J *
const head,
const J *
const tail)
const {
1103 double NN_dist = _R2;
1105 for (J * jetB = head; jetB != tail; jetB++) {
1106 double dist = _bj_dist(jet,jetB);
1107 if (dist < NN_dist) {
1111 if (dist < jetB->NN_dist) {
1112 jetB->NN_dist = dist;
1117 jet->NN_dist = NN_dist;
1120 FASTJET_END_NAMESPACE
1122 #endif // __FASTJET_CLUSTERSEQUENCE_HH__ bool will_delete_self_when_unused() const
return true if the object has been told to delete itself when unused
ClusterSequenceStructure StructureType
the structure type associated with a jet belonging to a ClusterSequence
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
void plugin_associate_extras(Extras *extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
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...
const Extras * extras() const
returns a pointer to the extras object (may be null)
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
double Q() const
returns the sum of all energies in the event (relevant mainly for e+e-)
void plugin_simple_N2_cluster()
allows a plugin to run a templated clustering (nearest-neighbour heuristic)
ClusterSequence()
default constructor
ClusterSequence(const ClusterSequence &cs)
copy constructor for a ClusterSequence
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.
double Q2() const
return Q()^2
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
Contains any information related to the clustering that should be directly accessible to PseudoJet...
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], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
std::vector< PseudoJet > _jets
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in ...
const JetDefinition & jet_def() const
return a reference to the jet definition
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
int parent2
index in _history where first parent of this jet was created (InexistentParent if this jet is an orig...
double dij
index in the _jets vector where we will find the
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...
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
an implementation of C++0x shared pointers (or boost's)
a single element in the clustering history
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
bool plugin_activated() const
returns true when the plugin is allowed to run the show.
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...
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
void _decant_options_partial()
assuming that the jet definition, writeout_combinations and _structure_shared_ptr have been set (e...
JetAlgorithm
the various families of jet-clustering algorithm
std::vector< PseudoJet > exclusive_jets_ycut(double ycut) const
the exclusive jets obtained at the given ycut
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
bool _deletes_self_when_unused
if true then the CS will delete itself when the last external object referring to it disappears...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
double max_dij_so_far
the distance corresponding to the recombination at this stage of the clustering.
int n_exclusive_jets_ycut(double ycut) const
the number of exclusive jets at the given ycut
double exclusive_ymerge_max(int njets) const
same as exclusive_dmerge_max, but normalised to squared total energy
class that is intended to hold a full definition of the jet clusterer
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...
std::string strategy_string() const
return the name of the strategy used to cluster the event