30 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__ 31 #define __FASTJET_CLUSTERSEQUENCE_HH__ 35 #include "fastjet/PseudoJet.hh" 42 #include "fastjet/Error.hh" 43 #include "fastjet/JetDefinition.hh" 44 #include "fastjet/SharedPtr.hh" 45 #include "fastjet/LimitedWarning.hh" 46 #include "fastjet/FunctionOfPseudoJet.hh" 47 #include "fastjet/ClusterSequenceStructure.hh" 49 FASTJET_BEGIN_NAMESPACE
53 class ClusterSequenceStructure;
54 class DynamicNearestNeighbours;
86 const std::vector<L> & pseudojets,
88 const bool & writeout_combinations =
false);
92 transfer_from_sequence(cs);
106 std::vector<PseudoJet> inclusive_jets (
const double & ptmin = 0.0)
const;
111 int n_exclusive_jets (
const double & dcut)
const;
116 std::vector<PseudoJet> exclusive_jets (
const double & dcut)
const;
123 std::vector<PseudoJet> exclusive_jets (
const int & njets)
const;
130 std::vector<PseudoJet> exclusive_jets_up_to (
const int & njets)
const;
135 double exclusive_dmerge (
const int & njets)
const;
141 double exclusive_dmerge_max (
const int & njets)
const;
155 int njets = n_exclusive_jets_ycut(ycut);
156 return exclusive_jets(njets);
170 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
171 const double & dcut)
const;
176 int n_exclusive_subjets(
const PseudoJet & jet,
177 const double & dcut)
const;
184 std::vector<PseudoJet> exclusive_subjets (
const PseudoJet & jet,
192 std::vector<PseudoJet> exclusive_subjets_up_to (
const PseudoJet & jet,
199 double exclusive_subdmerge(
const PseudoJet & jet,
int nsub)
const;
206 double exclusive_subdmerge_max(
const PseudoJet & jet,
int nsub)
const;
213 double Q()
const {
return _Qtot;}
215 double Q2()
const {
return _Qtot*_Qtot;}
252 std::vector<PseudoJet> constituents (
const PseudoJet & jet)
const;
263 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
264 std::ostream & ostr = std::cout)
const;
268 void print_jets_for_root(
const std::vector<PseudoJet> & jets,
269 const std::string & filename,
270 const std::string & comment =
"")
const;
277 void add_constituents (
const PseudoJet & jet,
278 std::vector<PseudoJet> & subjet_vector)
const;
287 std::string strategy_string (
Strategy strategy_in)
const;
305 void delete_self_when_unused();
312 void signal_imminent_self_deletion()
const;
317 double jet_scale_for_algorithm(
const PseudoJet & jet)
const;
331 assert(plugin_activated());
332 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
338 void plugin_record_ij_recombination(
int jet_i,
int jet_j,
double dij,
347 assert(plugin_activated());
348 _do_iB_recombination_step(jet_i, diB);
361 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";}
368 _extras.reset(extras_in.release());
385 assert(plugin_activated());
386 _simple_N2_cluster<GBJ>();
434 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
452 const std::vector<PseudoJet> & jets()
const;
465 const std::vector<history_element> & history()
const;
471 unsigned int n_particles()
const;
477 std::vector<int> particle_jet_indices(
const std::vector<PseudoJet> &)
const;
494 std::vector<int> unique_history_order()
const;
499 std::vector<PseudoJet> unclustered_particles()
const;
505 std::vector<PseudoJet> childless_pseudojets()
const;
512 bool contains(
const PseudoJet &
object)
const;
531 return _structure_shared_ptr;
544 static void print_banner();
559 static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
573 static std::ostream * _fastjet_banner_ostr;
583 template<
class L>
void _transfer_input_jets(
584 const std::vector<L> & pseudojets);
590 const bool & writeout_combinations);
596 void _initialise_and_run_no_decant();
608 const bool & writeout_combinations);
614 void _decant_options_partial();
619 void _fill_initial_history();
624 void _do_ij_recombination_step(
const int & jet_i,
const int & jet_j,
625 const double & dij,
int & newjet_k);
629 void _do_iB_recombination_step(
const int & jet_i,
const double & diB);
636 void _set_structure_shared_ptr(
PseudoJet & j);
640 void _update_structure_use_count();
662 void get_subhist_set(std::set<const history_element*> & subhist,
663 const PseudoJet & jet,
double dcut,
int maxjet)
const;
665 bool _writeout_combinations;
667 double _Rparam, _R2, _invR2;
673 int _structure_use_count_after_construction;
682 bool _plugin_activated;
686 void _really_dumb_cluster ();
687 void _delaunay_cluster ();
689 template<
class BJ>
void _simple_N2_cluster ();
690 void _tiled_N2_cluster ();
691 void _faster_tiled_N2_cluster ();
694 void _minheap_faster_tiled_N2_cluster();
698 void _CP2DChan_cluster();
699 void _CP2DChan_cluster_2pi2R ();
700 void _CP2DChan_cluster_2piMultD ();
701 void _CP2DChan_limited_cluster(
double D);
702 void _do_Cambridge_inclusive_jets();
705 void _fast_NsqrtN_cluster();
707 void _add_step_to_history(
const int & step_number,
const int & parent1,
708 const int & parent2,
const int & jetp_index,
713 void _extract_tree_children(
int pos, std::valarray<bool> &,
714 const std::valarray<int> &, std::vector<int> &)
const;
718 void _extract_tree_parents (
int pos, std::valarray<bool> &,
719 const std::valarray<int> &, std::vector<int> &)
const;
723 typedef std::pair<int,int> TwoVertices;
724 typedef std::pair<double,TwoVertices> DijEntry;
725 typedef std::multimap<double,TwoVertices> DistMap;
728 void _add_ktdistance_to_map(
const int & ii,
730 const DynamicNearestNeighbours * DNN);
734 static bool _first_time;
739 static int _n_exclusive_warnings;
752 double eta, phi, kt2, NN_dist;
762 double eta, phi, kt2, NN_dist;
763 TiledJet * NN, *previous, * next;
764 int _jets_index, tile_index, diJ_posn;
769 inline void label_minheap_update_needed() {diJ_posn = 1;}
770 inline void label_minheap_update_done() {diJ_posn = 0;}
771 inline bool minheap_update_needed()
const {
return diJ_posn==1;}
779 template <
class J>
void _bj_set_jetinfo( J *
const jet,
780 const int _jets_index)
const;
784 void _bj_remove_from_tiles( TiledJet *
const jet)
const;
787 template <
class J>
double _bj_dist(
const J *
const jeta,
788 const J *
const jetb)
const;
792 template <
class J>
double _bj_diJ(
const J *
const jeta)
const;
797 template <
class J>
inline J * _bj_of_hindex(
798 const int hist_index,
799 J *
const head, J *
const tail)
802 for(res = head; res<tail; res++) {
803 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {
break;}
816 template <
class J>
void _bj_set_NN_nocross(J *
const jeta,
817 J *
const head,
const J *
const tail)
const;
823 template <
class J>
void _bj_set_NN_crosscheck(J *
const jeta,
824 J *
const head,
const J *
const tail)
const;
830 static const int n_tile_neighbours = 9;
836 Tile * begin_tiles[n_tile_neighbours];
838 Tile ** surrounding_tiles;
848 std::vector<Tile> _tiles;
849 double _tiles_eta_min, _tiles_eta_max;
850 double _tile_size_eta, _tile_size_phi;
851 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
855 inline int _tile_index (
int ieta,
int iphi)
const {
858 return (ieta-_tiles_ieta_min)*_n_tiles_phi
859 + (iphi+_n_tiles_phi) % _n_tiles_phi;
864 int _tile_index(
const double & eta,
const double & phi)
const;
865 void _tj_set_jetinfo ( TiledJet *
const jet,
const int _jets_index);
866 void _bj_remove_from_tiles(TiledJet *
const jet);
867 void _initialise_tiles();
868 void _print_tiles(TiledJet * briefjets )
const;
869 void _add_neighbours_to_tile_union(
const int tile_index,
870 std::vector<int> & tile_union,
int & n_near_tiles)
const;
871 void _add_untagged_neighbours_to_tile_union(
const int tile_index,
872 std::vector<int> & tile_union,
int & n_near_tiles);
891 void _simple_N2_cluster_BriefJet();
893 void _simple_N2_cluster_EEBriefJet();
904 template<
class L>
void ClusterSequence::_transfer_input_jets(
905 const std::vector<L> & pseudojets) {
909 _jets.reserve(pseudojets.size()*2);
915 for (
unsigned int i = 0; i < pseudojets.size(); i++) {
916 _jets.push_back(pseudojets[i]);}
940 template<
class L> ClusterSequence::ClusterSequence (
941 const std::vector<L> & pseudojets,
943 const bool & writeout_combinations) :
944 _jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
955 _initialise_and_run_no_decant();
972 template <
class J>
inline void ClusterSequence::_bj_set_jetinfo(
973 J *
const jetA,
const int _jets_index)
const {
974 jetA->eta = _jets[_jets_index].rap();
975 jetA->phi = _jets[_jets_index].phi_02pi();
977 jetA->_jets_index = _jets_index;
987 template <
class J>
inline double ClusterSequence::_bj_dist(
988 const J *
const jetA,
const J *
const jetB)
const {
989 double dphi = std::abs(jetA->phi - jetB->phi);
990 double deta = (jetA->eta - jetB->eta);
991 if (dphi > pi) {dphi = twopi - dphi;}
992 return dphi*dphi + deta*deta;
996 template <
class J>
inline double ClusterSequence::_bj_diJ(
const J *
const jet)
const {
997 double kt2 = jet->kt2;
998 if (jet->NN != NULL) {
if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
999 return jet->NN_dist * kt2;
1006 template <
class J>
inline void ClusterSequence::_bj_set_NN_nocross(
1007 J *
const jet, J *
const head,
const J *
const tail)
const {
1008 double NN_dist = _R2;
1011 for (J * jetB = head; jetB != jet; jetB++) {
1012 double dist = _bj_dist(jet,jetB);
1013 if (dist < NN_dist) {
1020 for (J * jetB = jet+1; jetB != tail; jetB++) {
1021 double dist = _bj_dist(jet,jetB);
1022 if (dist < NN_dist) {
1029 jet->NN_dist = NN_dist;
1034 template <
class J>
inline void ClusterSequence::_bj_set_NN_crosscheck(J *
const jet,
1035 J *
const head,
const J *
const tail)
const {
1036 double NN_dist = _R2;
1038 for (J * jetB = head; jetB != tail; jetB++) {
1039 double dist = _bj_dist(jet,jetB);
1040 if (dist < NN_dist) {
1044 if (dist < jetB->NN_dist) {
1045 jetB->NN_dist = dist;
1050 jet->NN_dist = NN_dist;
1053 FASTJET_END_NAMESPACE
1055 #endif // __FASTJET_CLUSTERSEQUENCE_HH__ const JetDefinition & jet_def() const
return a reference to the jet definition
ClusterSequenceStructure StructureType
the structure type associated with a jet belonging to 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...
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)
const std::vector< history_element > & history() const
allow the user to access the raw internal history.
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
int jetp_index
index in _history where the current jet is recombined with another jet to form its child...
double Q() const
returns the sum of all energies in the event (relevant mainly for e+e-)
double Q2() const
return Q()^2
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 ...
std::vector< history_element > _history
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates wher...
bool plugin_activated() const
returns true when the plugin is allowed to run the show.
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...
std::vector< PseudoJet > exclusive_jets_ycut(double ycut) const
the exclusive jets obtained at the given ycut
double exclusive_ymerge_max(int njets) const
same as exclusive_dmerge_max, but normalised to squared total energy
Strategy strategy_used() const
return the enum value of the strategy used to cluster the event
int child
index in _history where second parent of this jet was created (InexistentParent if this jet is an ori...
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
an implementation of C++0x shared pointers (or boost's)
std::string strategy_string() const
return the name of the strategy used to cluster the event
void plugin_associate_extras(std::auto_ptr< Extras > extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
a single element in the clustering history
int n_exclusive_jets_ycut(double ycut) const
the number of exclusive jets at the given ycut
bool will_delete_self_when_unused() const
return true if the object has been told to delete itself when unused
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
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.
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
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...
class that is intended to hold a full definition of the jet clusterer
unsigned int n_particles() const
returns the number of particles that were provided to the clustering algorithm (helps the user find t...