deals with clustering More...
#include <ClusterSequence.hh>
Classes | |
| struct | BriefJet |
| the fundamental structure which contains the minimal info about a jet, as needed for our plain N^2 algorithm -- the idea is to put all info that will be accessed N^2 times into an array of BriefJets... More... | |
| struct | EEBriefJet |
| fundamental structure for e+e- clustering More... | |
| class | Extras |
| a class intended to serve as a base in case a plugin needs to associate extra information with a ClusterSequence (see SISConePlugin. More... | |
| struct | history_element |
| a single element in the clustering history (see vector _history below). More... | |
| struct | Tile |
| The fundamental structures to be used for the tiled N^2 algorithm (see CCN27-44 for some discussion of pattern of tiling) More... | |
| class | TiledJet |
| structure analogous to BriefJet, but with the extra information needed for dealing with tiles More... | |
Public Types | |
| enum | JetType { Invalid = -3, InexistentParent = -2, BeamJet = -1 } |
Public Member Functions | |
| ClusterSequence () | |
| default constructor | |
| template<class L > | |
| ClusterSequence (const std::vector< L > &pseudojets, const double &R=1.0, const Strategy &strategy=Best, const bool &writeout_combinations=false) | |
| create a clustersequence starting from the supplied set of pseudojets and clustering them with the long-invariant kt algorithm (E-scheme recombination) with the supplied value for R. | |
| template<class L > | |
| ClusterSequence (const std::vector< L > &pseudojets, const JetDefinition &jet_def, const bool &writeout_combinations=false) | |
| create a clustersequence starting from the supplied set of pseudojets and clustering them with jet definition specified by jet_def (which also specifies the clustering strategy) | |
| virtual | ~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. | |
| int | n_exclusive_jets (const double &dcut) const |
| return the number of jets (in the sense of the exclusive algorithm) that would be obtained when running the algorithm with the given dcut. | |
| std::vector< PseudoJet > | exclusive_jets (const double &dcut) const |
| return a vector of all jets (in the sense of the exclusive algorithm) that would be obtained when running the algorithm with the given dcut. | |
| std::vector< PseudoJet > | exclusive_jets (const int &njets) const |
| return a vector of all jets when the event is clustered (in the exclusive sense) to exactly njets. | |
| double | exclusive_dmerge (const int &njets) const |
| return the dmin corresponding to the recombination that went from n+1 to n jets (sometimes known as d_{n n+1}). | |
| double | exclusive_dmerge_max (const int &njets) const |
| return the maximum of the dmin encountered during all recombinations up to the one that led to an n-jet final state; identical to exclusive_dmerge, except in cases where the dmin do not increase monotonically. | |
| 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_{n n+1}). | |
| double | exclusive_ymerge_max (int njets) const |
| same as exclusive_dmerge_max, but normalised to squared total energy | |
| int | n_exclusive_jets_ycut (double ycut) const |
| the number of exclusive jets at the given ycut | |
| std::vector< PseudoJet > | exclusive_jets_ycut (double ycut) const |
| the exclusive jets obtained at the given ycut | |
| std::vector< PseudoJet > | exclusive_subjets (const PseudoJet &jet, const double &dcut) const |
| return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that would be obtained when running the algorithm with the given dcut. | |
| int | n_exclusive_subjets (const PseudoJet &jet, const double &dcut) const |
| return the size of exclusive_subjets(...); still n ln n with same coefficient, but marginally more efficient than manually taking exclusive_subjets.size() | |
| std::vector< PseudoJet > | exclusive_subjets (const PseudoJet &jet, int nsub) const |
| return the list of subjets obtained by unclustering the supplied jet down to n subjets (or all constituents if there are fewer than n). | |
| double | exclusive_subdmerge (const PseudoJet &jet, int nsub) const |
| return the dij that was present in the merging nsub+1 -> nsub subjets inside this jet. | |
| double | exclusive_subdmerge_max (const PseudoJet &jet, int nsub) const |
| return the maximum dij that occurred in the whole event at the stage that the nsub+1 -> nsub merge of subjets occurred inside this jet. | |
| double | Q () const |
| returns the sum of all energies in the event (relevant mainly for e+e-) | |
| double | Q2 () const |
| return Q()^2 | |
| bool | object_in_jet (const PseudoJet &object, const PseudoJet &jet) const |
| returns true iff the object is included in the jet. | |
| bool | has_parents (const PseudoJet &jet, PseudoJet &parent1, PseudoJet &parent2) const |
| if the jet has parents in the clustering, it returns true and sets parent1 and parent2 equal to them. | |
| bool | has_child (const PseudoJet &jet, PseudoJet &child) const |
| if the jet has a child then return true and give the child jet otherwise return false and set the child to zero | |
| bool | has_child (const PseudoJet &jet, const PseudoJet *&childp) const |
| Version of has_child that sets a pointer to the child if the child exists;. | |
| bool | has_partner (const PseudoJet &jet, PseudoJet &partner) const |
| if this jet has a child (and so a partner) return true and give the partner, otherwise return false and set the partner to zero | |
| std::vector< PseudoJet > | constituents (const PseudoJet &jet) const |
| return a vector of the particles that make up jet | |
| 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 format is: jet-n jet-px jet-py jet-pz jet-E particle-n particle-rap particle-phi particle-pt particle-n particle-rap particle-phi particle-pt ... | |
| void | print_jets_for_root (const std::vector< PseudoJet > &jets, const std::string &filename, const std::string &comment="") const |
| print jets for root to the file labelled filename, with an optional comment at the beginning | |
| void | add_constituents (const PseudoJet &jet, std::vector< PseudoJet > &subjet_vector) const |
| add on to subjet_vector the constituents of jet (for internal use mainly) | |
| Strategy | strategy_used () const |
| return the enum value of the strategy used to cluster the event | |
| std::string | strategy_string () const |
| const JetDefinition & | jet_def () const |
| return a reference to the jet definition | |
| 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-algorithm, 1 for the Cambridge algorithm). | |
| 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 | |
| void | plugin_record_ij_recombination (int jet_i, int jet_j, double dij, const PseudoJet &newjet, int &newjet_k) |
| as for the simpler variant of plugin_record_ij_recombination, except that the new jet is attributed the momentum and user_index of newjet | |
| 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, with the specified diB; when looking for inclusive jets, any iB recombination will returned to the user as a jet. | |
| void | plugin_associate_extras (std::auto_ptr< Extras > extras_in) |
| the plugin can associated some extra information with the ClusterSequence object by calling this function | |
| bool | plugin_activated () const |
| returns true when the plugin is allowed to run the show. | |
| const Extras * | extras () const |
| returns a pointer to the extras object (may be null) | |
| template<class GBJ > | |
| void | plugin_simple_N2_cluster () |
| allows a plugin to run a templated clustering (nearest-neighbour heuristic) | |
| const std::vector< PseudoJet > & | jets () const |
| allow the user to access the internally stored _jets() array, which contains both the initial particles and the various intermediate and final stages of recombination. | |
| const std::vector < history_element > & | history () const |
| allow the user to access the raw internal history. | |
| unsigned int | n_particles () const |
| returns the number of particles that were provided to the clustering algorithm (helps the user find their way around the history and jets objects if they weren't paying attention beforehand). | |
| std::vector< int > | particle_jet_indices (const std::vector< PseudoJet > &) const |
| returns a vector of size n_particles() which indicates, for each of the initial particles (in the order in which they were supplied), which of the supplied jets it belongs to; if it does not belong to any of the supplied jets, the index is set to -1; | |
| std::vector< int > | unique_history_order () const |
| routine that returns an order in which to read the history such that clusterings that lead to identical jet compositions but different histories (because of degeneracies in the clustering order) will have matching constituents for each matching entry in the unique_history_order. | |
| std::vector< PseudoJet > | unclustered_particles () const |
| return the set of particles that have not been clustered. | |
| void | transfer_from_sequence (ClusterSequence &from_seq) |
| transfer the sequence contained in other_seq into our own; any plugin "extras" contained in the from_seq will be lost from there. | |
| template<> | |
| void | _bj_set_jetinfo (EEBriefJet *const jetA, const int _jets_index) const |
| template<> | |
| double | _bj_dist (const EEBriefJet *const jeta, const EEBriefJet *const jetb) const |
Static Public Member Functions | |
| static void | set_jet_algorithm (JetAlgorithm jet_algorithm) |
| set the default (static) jet finder across all current and future ClusterSequence objects -- deprecated and obsolescent (i.e. | |
| static void | set_jet_finder (JetAlgorithm jet_algorithm) |
| same as above for backward compatibility | |
Protected Member Functions | |
| bool | _potentially_valid (const PseudoJet &jet) const |
| returns true if the jet has a history index contained within the range of this CS | |
| template<class L > | |
| 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 for future growth). | |
| void | _initialise_and_run (const JetDefinition &jet_def, const bool &writeout_combinations) |
| This is the routine that will do all the initialisation and then run the clustering (may be called by various constructors). | |
| void | _initialise_and_run (const double &R, const Strategy &strategy, const bool &writeout_combinations) |
| This is an alternative routine for initialising and running the clustering, provided for legacy purposes. | |
| void | _decant_options (const JetDefinition &jet_def, const bool &writeout_combinations) |
| fills in the various member variables with "decanted" options from the jet_definition and writeout_combinations variables | |
| void | _fill_initial_history () |
| fill out the history (and jet cross refs) related to the initial set of jets (assumed already to have been "transferred"), without any clustering | |
| void | _do_ij_recombination_step (const int &jet_i, const int &jet_j, const double &dij, int &newjet_k) |
| carry out the recombination between the jets numbered jet_i and jet_j, at distance scale dij; return the index newjet_k of the result of the recombination of i and j. | |
| void | _do_iB_recombination_step (const int &jet_i, const double &diB) |
| carry out an recombination step in which _jets[jet_i] merges with the beam, | |
| void | get_subhist_set (std::set< const history_element * > &subhist, const PseudoJet &jet, double dcut, int maxjet) const |
set subhist to be a set pointers to history entries corresponding to the subjets of this jet; one stops going working down through the subjets either when
| |
Protected Attributes | |
| JetDefinition | _jet_def |
| std::vector< PseudoJet > | _jets |
| This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in the _history by looking at _jets[i].cluster_hist_index(). | |
| std::vector< history_element > | _history |
| this vector will contain the branching history; for each stage, _history[i].jetp_index indicates where to look in the _jets vector to get the physical PseudoJet. | |
| bool | _writeout_combinations |
| int | _initial_n |
| double | _Rparam |
| double | _R2 |
| double | _invR2 |
| double | _Qtot |
| Strategy | _strategy |
| JetAlgorithm | _jet_algorithm |
Static Protected Attributes | |
| static JetAlgorithm | _default_jet_algorithm = kt_algorithm |
Private Types | |
| typedef std::pair< int, int > | TwoVertices |
| typedef std::pair< double, TwoVertices > | DijEntry |
| typedef std::multimap< double, TwoVertices > | DistMap |
Private Member Functions | |
| void | _really_dumb_cluster () |
| Run the clustering in a very slow variant of the N^3 algorithm. | |
| void | _delaunay_cluster () |
| Run the clustering using a Hierarchical Delaunay triangulation and STL maps to achieve O(N*ln N) behaviour. | |
| template<class BJ > | |
| void | _simple_N2_cluster () |
| void | _tiled_N2_cluster () |
| run a tiled clustering | |
| void | _faster_tiled_N2_cluster () |
| run a tiled clustering | |
| void | _minheap_faster_tiled_N2_cluster () |
| run a tiled clustering, with our minheap for keeping track of the smallest dij | |
| void | _CP2DChan_cluster () |
| a 4pi variant of the closest pair clustering | |
| void | _CP2DChan_cluster_2pi2R () |
| a variant of the closest pair clustering which uses a region of size 2pi+2R in phi. | |
| void | _CP2DChan_cluster_2piMultD () |
| a variant of the closest pair clustering which uses a region of size 2pi + 2*0.3 and then carries on with 2pi+2*R | |
| void | _CP2DChan_limited_cluster (double D) |
| clusters only up to a distance Dlim -- does not deal with "inclusive" jets -- these are left to some other part of the program | |
| void | _do_Cambridge_inclusive_jets () |
| void | _add_step_to_history (const int &step_number, const int &parent1, const int &parent2, const int &jetp_index, const double &dij) |
| void | _extract_tree_children (int pos, std::valarray< bool > &, const std::valarray< int > &, std::vector< int > &) const |
| internal routine associated with the construction of the unique history order (following children in the tree) | |
| void | _extract_tree_parents (int pos, std::valarray< bool > &, const std::valarray< int > &, std::vector< int > &) const |
| internal routine associated with the construction of the unique history order (following parents in the tree) | |
| void | _add_ktdistance_to_map (const int &ii, DistMap &DijMap, const DynamicNearestNeighbours *DNN) |
| currently used only in the Voronoi based code | |
| void | _print_banner () |
| for making sure the user knows what it is they're running... | |
| template<class J > | |
| void | _bj_set_jetinfo (J *const jet, const int _jets_index) const |
| set the kinematic and labelling info for jeta so that it corresponds to _jets[_jets_index] | |
| void | _bj_remove_from_tiles (TiledJet *const jet) const |
| "remove" this jet, which implies updating links of neighbours and perhaps modifying the tile structure | |
| template<class J > | |
| double | _bj_dist (const J *const jeta, const J *const jetb) const |
| return the distance between two BriefJet objects | |
| template<class J > | |
| double | _bj_diJ (const J *const jeta) const |
| template<class J > | |
| J * | _bj_of_hindex (const int hist_index, J *const head, J *const tail) const |
| for testing purposes only: if in the range head--tail-1 there is a a jet which corresponds to hist_index in the history, then return a pointer to that jet; otherwise return tail. | |
| template<class J > | |
| void | _bj_set_NN_nocross (J *const jeta, J *const head, const J *const tail) const |
| updates (only towards smaller distances) the NN for jeta without checking whether in the process jeta itself might be a new NN of one of the jets being scanned -- span the range head to tail-1 with assumption that jeta is not contained in that range | |
| template<class J > | |
| void | _bj_set_NN_crosscheck (J *const jeta, J *const head, const J *const tail) const |
| reset the NN for jeta and DO check whether in the process jeta itself might be a new NN of one of the jets being scanned -- span the range head to tail-1 with assumption that jeta is not contained in that range | |
| int | _tile_index (int ieta, int iphi) const |
| int | _tile_index (const double &eta, const double &phi) const |
| return the tile index corresponding to the given eta,phi point | |
| void | _tj_set_jetinfo (TiledJet *const jet, const int _jets_index) |
| void | _bj_remove_from_tiles (TiledJet *const jet) |
| void | _initialise_tiles () |
Set up the tiles:
| |
| void | _print_tiles (TiledJet *briefjets) const |
| output the contents of the tiles | |
| void | _add_neighbours_to_tile_union (const int tile_index, std::vector< int > &tile_union, int &n_near_tiles) const |
| Add to the vector tile_union the tiles that are in the neighbourhood of the specified tile_index, including itself -- start adding from position n_near_tiles-1, and increase n_near_tiles as you go along (could have done it more C++ like with vector with reserved space, but fear is that it would have been slower, e.g. | |
| void | _add_untagged_neighbours_to_tile_union (const int tile_index, std::vector< int > &tile_union, int &n_near_tiles) |
| Like _add_neighbours_to_tile_union, but only adds neighbours if their "tagged" status is false; when a neighbour is added its tagged status is set to true. | |
| void | _simple_N2_cluster_BriefJet () |
| to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?) | |
| void | _simple_N2_cluster_EEBriefJet () |
| to avoid issues with template instantiation (OS X 10.3, gcc 3.3) | |
Private Attributes | |
| bool | _plugin_activated |
| std::auto_ptr< Extras > | _extras |
| std::vector< Tile > | _tiles |
| double | _tiles_eta_min |
| double | _tiles_eta_max |
| double | _tile_size_eta |
| double | _tile_size_phi |
| int | _n_tiles_phi |
| int | _tiles_ieta_min |
| int | _tiles_ieta_max |
Static Private Attributes | |
| static bool | _first_time = true |
| will be set by default to be true for the first run | |
| static int | _n_exclusive_warnings = 0 |
| record the number of warnings provided about the exclusive algorithm -- so that we don't print it out more than a few times. | |
| static const int | n_tile_neighbours = 9 |
| number of neighbours that a tile will have (rectangular geometry gives 9 neighbours). | |
Detailed Description
deals with clustering
Definition at line 66 of file ClusterSequence.hh.
Member Typedef Documentation
typedef std::pair<double,TwoVertices> ClusterSequence::DijEntry [private] |
Definition at line 601 of file ClusterSequence.hh.
typedef std::multimap<double,TwoVertices> ClusterSequence::DistMap [private] |
Definition at line 602 of file ClusterSequence.hh.
typedef std::pair<int,int> ClusterSequence::TwoVertices [private] |
Definition at line 600 of file ClusterSequence.hh.
Member Enumeration Documentation
Definition at line 408 of file ClusterSequence.hh.
{Invalid=-3, InexistentParent = -2, BeamJet = -1};
Constructor & Destructor Documentation
| ClusterSequence::ClusterSequence | ( | ) | [inline] |
| ClusterSequence::ClusterSequence | ( | const std::vector< L > & | pseudojets, |
| const double & | R = 1.0, |
||
| const Strategy & | strategy = Best, |
||
| const bool & | writeout_combinations = false |
||
| ) |
create a clustersequence starting from the supplied set of pseudojets and clustering them with the long-invariant kt algorithm (E-scheme recombination) with the supplied value for R.
If strategy=DumbN3 a very stupid N^3 algorithm is used for the clustering; otherwise strategy = NlnN* uses cylinders algorithms with some number of pi coverage. If writeout_combinations=true a summary of the recombination sequence is written out
Definition at line 796 of file ClusterSequence.hh.
References _initialise_and_run(), and _transfer_input_jets().
{
// transfer the initial jets (type L) into our own array
_transfer_input_jets(pseudojets);
// run the clustering
_initialise_and_run(R,strategy,writeout_combinations);
}
| ClusterSequence::ClusterSequence | ( | const std::vector< L > & | pseudojets, |
| const JetDefinition & | jet_def, | ||
| const bool & | writeout_combinations = false |
||
| ) |
create a clustersequence starting from the supplied set of pseudojets and clustering them with jet definition specified by jet_def (which also specifies the clustering strategy)
constructor of a jet-clustering sequence from a vector of four-momenta, with the jet definition specified by jet_def
Definition at line 813 of file ClusterSequence.hh.
References _initialise_and_run(), and _transfer_input_jets().
{
// transfer the initial jets (type L) into our own array
_transfer_input_jets(pseudojets);
// run the clustering
_initialise_and_run(jet_def,writeout_combinations);
}
| ClusterSequence::~ClusterSequence | ( | ) | [virtual] |
Definition at line 54 of file ClusterSequence.cc.
{}
Member Function Documentation
| void ClusterSequence::_add_ktdistance_to_map | ( | const int & | ii, |
| DistMap & | DijMap, | ||
| const DynamicNearestNeighbours * | DNN | ||
| ) | [private] |
currently used only in the Voronoi based code
Add the current kt distance for particle ii to the map (DijMap) using information from the DNN object.
Work as follows:
. if the kt is zero then it's nearest neighbour is taken to be the the beam jet and the distance is zero.
. if cylinder distance to nearest neighbour > _Rparam then it is yiB that is smallest and this is added to map.
. otherwise if the nearest neighbour jj has a larger kt then add dij to the map.
. otherwise do nothing
Definition at line 220 of file ClusterSequence_Delaunay.cc.
References DynamicNearestNeighbours::NearestNeighbourDistance(), and DynamicNearestNeighbours::NearestNeighbourIndex().
{
double yiB = jet_scale_for_algorithm(_jets[ii]);
if (yiB == 0.0) {
// in this case convention is that we do not worry about distances
// but directly state that nearest neighbour is beam
DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
} else {
double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
// Logic of following bit is: only add point to map if it has
// smaller kt2 than nearest neighbour j (if it has larger kt,
// then: either it is j's nearest neighbour and then we will
// include dij when we come to j; or it is not j's nearest
// neighbour and j will recombine with someone else).
// If DeltaR2 > 1.0 then in any case it will recombine with beam rather
// than with any neighbours.
// (put general normalisation here at some point)
if (DeltaR2 > 1.0) {
DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
} else {
double kt2i = jet_scale_for_algorithm(_jets[ii]);
int jj = DNN->NearestNeighbourIndex(ii);
if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
double dij = DeltaR2 * kt2i;
DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
}
}
}
}
| void ClusterSequence::_add_neighbours_to_tile_union | ( | const int | tile_index, |
| std::vector< int > & | tile_union, | ||
| int & | n_near_tiles | ||
| ) | const [private] |
Add to the vector tile_union the tiles that are in the neighbourhood of the specified tile_index, including itself -- start adding from position n_near_tiles-1, and increase n_near_tiles as you go along (could have done it more C++ like with vector with reserved space, but fear is that it would have been slower, e.g.
checking for end of vector at each stage to decide whether to resize it)
Definition at line 235 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::Tile::end_tiles.
| void ClusterSequence::_add_step_to_history | ( | const int & | step_number, |
| const int & | parent1, | ||
| const int & | parent2, | ||
| const int & | jetp_index, | ||
| const double & | dij | ||
| ) | [private] |
Definition at line 939 of file ClusterSequence.cc.
References ClusterSequence::history_element::child, ClusterSequence::history_element::dij, ClusterSequence::history_element::jetp_index, ClusterSequence::history_element::max_dij_so_far, ClusterSequence::history_element::parent1, and ClusterSequence::history_element::parent2.
{
history_element element;
element.parent1 = parent1;
element.parent2 = parent2;
element.jetp_index = jetp_index;
element.child = Invalid;
element.dij = dij;
element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
_history.push_back(element);
int local_step = _history.size()-1;
assert(local_step == step_number);
assert(parent1 >= 0);
_history[parent1].child = local_step;
if (parent2 >= 0) {_history[parent2].child = local_step;}
// get cross-referencing right from PseudoJets
if (jetp_index != Invalid) {
assert(jetp_index >= 0);
//cout << _jets.size() <<" "<<jetp_index<<"\n";
_jets[jetp_index].set_cluster_hist_index(local_step);
}
if (_writeout_combinations) {
cout << local_step << ": "
<< parent1 << " with " << parent2
<< "; y = "<< dij<<endl;
}
}
| void ClusterSequence::_add_untagged_neighbours_to_tile_union | ( | const int | tile_index, |
| std::vector< int > & | tile_union, | ||
| int & | n_near_tiles | ||
| ) | [inline, private] |
Like _add_neighbours_to_tile_union, but only adds neighbours if their "tagged" status is false; when a neighbour is added its tagged status is set to true.
Definition at line 250 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::Tile::end_tiles.
| double ClusterSequence::_bj_diJ | ( | const J *const | jeta | ) | const [inline, private] |
Definition at line 863 of file ClusterSequence.hh.
{
double kt2 = jet->kt2;
if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
return jet->NN_dist * kt2;
}
| double ClusterSequence::_bj_dist | ( | const J *const | jeta, |
| const J *const | jetb | ||
| ) | const [inline, private] |
return the distance between two BriefJet objects
Definition at line 854 of file ClusterSequence.hh.
References fastjet::pi, and twopi.
Referenced by _bj_set_NN_crosscheck(), and _bj_set_NN_nocross().
| double ClusterSequence::_bj_dist | ( | const EEBriefJet *const | jeta, |
| const EEBriefJet *const | jetb | ||
| ) | const |
Definition at line 95 of file ClusterSequence_N2.cc.
References ClusterSequence::EEBriefJet::nx, ClusterSequence::EEBriefJet::ny, and ClusterSequence::EEBriefJet::nz.
{
double dist = 1.0
- jeta->nx*jetb->nx
- jeta->ny*jetb->ny
- jeta->nz*jetb->nz;
dist *= 2; // distance is _2_*min(Ei^2,Ej^2)*(1-cos theta)
//cout << "Dist = " << dist << ": "
// << jeta->nx << " "
// << jeta->ny << " "
// << jeta->nz << " "
// <<endl;
return dist;
}
| J* ClusterSequence::_bj_of_hindex | ( | const int | hist_index, |
| J *const | head, | ||
| J *const | tail | ||
| ) | const [inline, private] |
for testing purposes only: if in the range head--tail-1 there is a a jet which corresponds to hist_index in the history, then return a pointer to that jet; otherwise return tail.
Definition at line 670 of file ClusterSequence.hh.
References _jets.
{
J * res;
for(res = head; res<tail; res++) {
if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
}
return res;
}
| void ClusterSequence::_bj_remove_from_tiles | ( | TiledJet *const | jet | ) | const [private] |
"remove" this jet, which implies updating links of neighbours and perhaps modifying the tile structure
| void ClusterSequence::_bj_remove_from_tiles | ( | TiledJet *const | jet | ) | [private] |
Definition at line 50 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::Tile::head, ClusterSequence::TiledJet::next, ClusterSequence::TiledJet::previous, and ClusterSequence::TiledJet::tile_index.
{
Tile * tile = & _tiles[jet->tile_index];
if (jet->previous == NULL) {
// we are at head of the tile, so reset it.
// If this was the only jet on the tile then tile->head will now be NULL
tile->head = jet->next;
} else {
// adjust link from previous jet in this tile
jet->previous->next = jet->next;
}
if (jet->next != NULL) {
// adjust backwards-link from next jet in this tile
jet->next->previous = jet->previous;
}
}
| void ClusterSequence::_bj_set_jetinfo | ( | J *const | jet, |
| const int | _jets_index | ||
| ) | const [inline, private] |
set the kinematic and labelling info for jeta so that it corresponds to _jets[_jets_index]
Definition at line 839 of file ClusterSequence.hh.
References _jets, _R2, and jet_scale_for_algorithm().
{
jetA->eta = _jets[_jets_index].rap();
jetA->phi = _jets[_jets_index].phi_02pi();
jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
jetA->_jets_index = _jets_index;
// initialise NN info as well
jetA->NN_dist = _R2;
jetA->NN = NULL;
}
| void ClusterSequence::_bj_set_jetinfo | ( | EEBriefJet *const | jetA, |
| const int | _jets_index | ||
| ) | const [inline] |
Definition at line 54 of file ClusterSequence_N2.cc.
References ClusterSequence::EEBriefJet::_jets_index, ee_genkt_algorithm, ee_kt_algorithm, ClusterSequence::EEBriefJet::kt2, ClusterSequence::EEBriefJet::NN, ClusterSequence::EEBriefJet::NN_dist, norm(), ClusterSequence::EEBriefJet::nx, ClusterSequence::EEBriefJet::ny, and ClusterSequence::EEBriefJet::nz.
{
double E = _jets[_jets_index].E();
double scale = E*E; // the default energy scale for the kt alg
double p = jet_def().extra_param(); // in case we're ee_genkt
switch (_jet_algorithm) {
case ee_kt_algorithm:
assert(_Rparam > 2.0); // force this to be true! [not best place, but works]
// recall that _invR2 is artificially set to 1 for this alg
// so that we automatically have dij = scale * 2(1-cos theta_ij)
// Normally, _Rparam should be automatically set to 4 from JetDefinition
break;
case ee_genkt_algorithm:
if (p <= 0 && scale < 1e-300) scale = 1e-300; // same dodgy safety as genkt
scale = pow(scale,p);
break;
default:
throw Error("Unrecognised jet algorithm");
}
jetA->kt2 = scale; // "kt2" might one day be renamed as "scale" or some such
double norm = _jets[_jets_index].modp2();
if (norm > 0) {
norm = 1.0/sqrt(norm);
jetA->nx = norm * _jets[_jets_index].px();
jetA->ny = norm * _jets[_jets_index].py();
jetA->nz = norm * _jets[_jets_index].pz();
} else {
jetA->nx = 0.0;
jetA->ny = 0.0;
jetA->nz = 1.0;
}
jetA->_jets_index = _jets_index;
// initialise NN info as well
jetA->NN_dist = _R2;
jetA->NN = NULL;
}
| void ClusterSequence::_bj_set_NN_crosscheck | ( | J *const | jeta, |
| J *const | head, | ||
| const J *const | tail | ||
| ) | const [inline, private] |
reset the NN for jeta and DO check whether in the process jeta itself might be a new NN of one of the jets being scanned -- span the range head to tail-1 with assumption that jeta is not contained in that range
Definition at line 901 of file ClusterSequence.hh.
References _bj_dist(), and _R2.
| void ClusterSequence::_bj_set_NN_nocross | ( | J *const | jeta, |
| J *const | head, | ||
| const J *const | tail | ||
| ) | const [inline, private] |
updates (only towards smaller distances) the NN for jeta without checking whether in the process jeta itself might be a new NN of one of the jets being scanned -- span the range head to tail-1 with assumption that jeta is not contained in that range
Definition at line 873 of file ClusterSequence.hh.
References _bj_dist(), and _R2.
{
double NN_dist = _R2;
J * NN = NULL;
if (head < jet) {
for (J * jetB = head; jetB != jet; jetB++) {
double dist = _bj_dist(jet,jetB);
if (dist < NN_dist) {
NN_dist = dist;
NN = jetB;
}
}
}
if (tail > jet) {
for (J * jetB = jet+1; jetB != tail; jetB++) {
double dist = _bj_dist(jet,jetB);
if (dist < NN_dist) {
NN_dist = dist;
NN = jetB;
}
}
}
jet->NN = NN;
jet->NN_dist = NN_dist;
}
| void ClusterSequence::_CP2DChan_cluster | ( | ) | [private] |
a 4pi variant of the closest pair clustering
Definition at line 223 of file ClusterSequence_CP2DChan.cc.
References cambridge_algorithm, ClosestPair2D::closest_pair(), d0::inline_maths::min(), ClosestPair2D::replace(), and twopi.
{
if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
unsigned int n = _jets.size();
vector<MirrorInfo> coordIDs(2*n); // link from original to mirror indices
vector<int> jetIDs(2*n); // link from mirror to original indices
vector<Coord2D> coords(2*n); // our coordinates (and copies)
// start things off...
double minrap = numeric_limits<double>::max();
double maxrap = -minrap;
int coord_index = 0;
for (unsigned i = 0; i < n; i++) {
// separate out points with infinite rapidity
if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
} else {
coordIDs[i].orig = coord_index;
coordIDs[i].mirror = coord_index+1;
coords[coord_index] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
jetIDs[coord_index] = i;
jetIDs[coord_index+1] = i;
minrap = min(coords[coord_index].x,minrap);
maxrap = max(coords[coord_index].x,maxrap);
coord_index += 2;
}
}
// label remaining "mirror info" as saying that there are no jets
for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
// set them to their actual size...
coords.resize(coord_index);
// establish limits (with some leeway on rapidity)
Coord2D left_edge(minrap-1.0, 0.0);
Coord2D right_edge(maxrap+1.0, 2*twopi);
// now create the closest pair search object
ClosestPair2D cp(coords, left_edge, right_edge);
// and start iterating...
vector<Coord2D> new_points(2);
vector<unsigned int> cIDs_to_remove(4);
vector<unsigned int> new_cIDs(2);
do {
// find out which pair we recombine
unsigned int cID1, cID2;
double distance2;
cp.closest_pair(cID1,cID2,distance2);
distance2 *= _invR2;
// if the closest distance > 1, we exit and go to "inclusive" stage
if (distance2 > 1.0) {break;}
// do the recombination...
int jet_i = jetIDs[cID1];
int jet_j = jetIDs[cID2];
int newjet_k;
_do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
// now prepare operations on CP structure
cIDs_to_remove[0] = coordIDs[jet_i].orig;
cIDs_to_remove[1] = coordIDs[jet_i].mirror;
cIDs_to_remove[2] = coordIDs[jet_j].orig;
cIDs_to_remove[3] = coordIDs[jet_j].mirror;
new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
// carry out the CP operation
//cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
// remarkable the following is faster...
new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
// signal that the following jets no longer exist
coordIDs[jet_i].orig = Invalid;
coordIDs[jet_j].orig = Invalid;
// and do bookkeeping for new jet
coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
jetIDs[new_cIDs[0]] = newjet_k;
jetIDs[new_cIDs[1]] = newjet_k;
// if we've reached one jet we should exit...
n--;
if (n == 1) {break;}
} while(true);
// now do "beam" recombinations
//for (unsigned int jet_i = 0; jet_i < coordIDs.size(); jet_i++) {
// // if we have a predeclared beam jet or a valid set of mirror
// // coordinates, recombine then recombine this jet with the beam
// if (coordIDs[jet_i].orig == BeamJet || coordIDs[jet_i].orig > 0) {
// // diB = 1 _always_ (for the cambridge alg.)
// _do_iB_recombination_step(jet_i, 1.0);
// }
//}
_do_Cambridge_inclusive_jets();
//n = _history.size();
//for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
// if (_history[hist_i].child == Invalid) {
// _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
// }
//}
}
| void ClusterSequence::_CP2DChan_cluster_2pi2R | ( | ) | [private] |
a variant of the closest pair clustering which uses a region of size 2pi+2R in phi.
Definition at line 193 of file ClusterSequence_CP2DChan.cc.
References cambridge_algorithm.
{
if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
// run the clustering with mirror copies kept such that only things
// within _Rparam of a border are mirrored
_CP2DChan_limited_cluster(_Rparam);
//
_do_Cambridge_inclusive_jets();
}
| void ClusterSequence::_CP2DChan_cluster_2piMultD | ( | ) | [private] |
a variant of the closest pair clustering which uses a region of size 2pi + 2*0.3 and then carries on with 2pi+2*R
Definition at line 209 of file ClusterSequence_CP2DChan.cc.
References d0::inline_maths::min().
{
// do a first run of clustering up to a small distance parameter,
if (_Rparam >= 0.39) {
_CP2DChan_limited_cluster(min(_Rparam/2,0.3));
}
// and then the final round of clustering
_CP2DChan_cluster_2pi2R ();
}
| void ClusterSequence::_CP2DChan_limited_cluster | ( | double | D | ) | [private] |
clusters only up to a distance Dlim -- does not deal with "inclusive" jets -- these are left to some other part of the program
Definition at line 69 of file ClusterSequence_CP2DChan.cc.
References ClosestPair2D::closest_pair(), Private::make_mirror(), d0::inline_maths::min(), fastjet::pi, and ClosestPair2D::replace_many().
{
unsigned int n = _initial_n;
vector<MirrorInfo> coordIDs(2*n); // coord IDs of a given jetID
vector<int> jetIDs(2*n); // jet ID for a given coord ID
vector<Coord2D> coords(2*n); // our coordinates (and copies)
// start things off...
double minrap = numeric_limits<double>::max();
double maxrap = -minrap;
int coord_index = -1;
int n_active = 0;
for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
// skip jets that already have children or that have infinite
// rapidity
if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
(_jets[jet_i].E() == abs(_jets[jet_i].pz()) &&
_jets[jet_i].perp2() == 0.0)
) {continue;}
n_active++;
coordIDs[jet_i].orig = ++coord_index;
coords[coord_index] = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
jetIDs[coord_index] = jet_i;
minrap = min(coords[coord_index].x,minrap);
maxrap = max(coords[coord_index].x,maxrap);
Coord2D mirror_point(coords[coord_index]);
if (make_mirror(mirror_point, Dlim)) {
coordIDs[jet_i].mirror = ++coord_index;
coords[coord_index] = mirror_point;
jetIDs[coord_index] = jet_i;
} else {
coordIDs[jet_i].mirror = Invalid;
}
}
// set them to their actual size...
coords.resize(coord_index+1);
// establish limits (with some leeway on rapidity)
Coord2D left_edge(minrap-1.0, -pi);
Coord2D right_edge(maxrap+1.0, 3*pi);
//cerr << "minrap, maxrap = " << minrap << " " << maxrap << endl;
// now create the closest pair search object
ClosestPair2D cp(coords, left_edge, right_edge);
// cross check to see what's going on...
//cerr << "Tree depths before: ";
//cp.print_tree_depths(cerr);
// and start iterating...
vector<Coord2D> new_points(2);
vector<unsigned int> cIDs_to_remove(4);
vector<unsigned int> new_cIDs(2);
do {
// find out which pair we recombine
unsigned int cID1, cID2;
double distance2;
cp.closest_pair(cID1,cID2,distance2);
// if the closest distance > Dlim, we exit and go to "inclusive" stage
if (distance2 > Dlim*Dlim) {break;}
// normalise distance as we like it
distance2 *= _invR2;
// do the recombination...
int jet_i = jetIDs[cID1];
int jet_j = jetIDs[cID2];
int newjet_k;
_do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
// don't bother with any further action if only one active particle
// is left (also avoid closest-pair error [cannot remove last particle]).
if (--n_active == 1) {break;}
// now prepare operations on CP structure
cIDs_to_remove.resize(0);
cIDs_to_remove.push_back(coordIDs[jet_i].orig);
cIDs_to_remove.push_back(coordIDs[jet_j].orig);
if (coordIDs[jet_i].mirror != Invalid)
cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
if (coordIDs[jet_j].mirror != Invalid)
cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
new_points.resize(0);
new_points.push_back(new_point);
if (make_mirror(new_point, Dlim)) new_points.push_back(new_point);
// carry out actions on search tree
cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
// now fill in info for new points...
coordIDs[newjet_k].orig = new_cIDs[0];
jetIDs[new_cIDs[0]] = newjet_k;
if (new_cIDs.size() == 2) {
coordIDs[newjet_k].mirror = new_cIDs[1];
jetIDs[new_cIDs[1]] = newjet_k;
} else {coordIDs[newjet_k].mirror = Invalid;}
//n_active--;
//if (n_active == 1) {break;}
} while(true);
// cross check to see what's going on...
//cerr << "Max tree depths after: ";
//cp.print_tree_depths(cerr);
}
| void ClusterSequence::_decant_options | ( | const JetDefinition & | jet_def, |
| const bool & | writeout_combinations | ||
| ) | [protected] |
fills in the various member variables with "decanted" options from the jet_definition and writeout_combinations variables
Definition at line 229 of file ClusterSequence.cc.
References JetDefinition::jet_algorithm(), JetDefinition::R(), and JetDefinition::strategy().
{
// let the user know what's going on
_print_banner();
// make a local copy of the jet definition (for future use?)
_jet_def = jet_def;
_writeout_combinations = writeout_combinations;
_jet_algorithm = jet_def.jet_algorithm();
_Rparam = jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
_strategy = jet_def.strategy();
// disallow interference from the plugin
_plugin_activated = false;
}
| void ClusterSequence::_delaunay_cluster | ( | ) | [private] |
Run the clustering using a Hierarchical Delaunay triangulation and STL maps to achieve O(N*ln N) behaviour.
There may be internally asserted assumptions about absence of points with coincident eta-phi coordinates.
Definition at line 58 of file ClusterSequence_Delaunay.cc.
References NlnN, NlnN3pi, NlnN4pi, DynamicNearestNeighbours::RemoveCombinedAddCombination(), DynamicNearestNeighbours::RemovePoint(), EtaPhi::sanitize(), twopi, and DynamicNearestNeighbours::Valid().
{
int n = _jets.size();
vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double>
for (int i = 0; i < n; i++) {
points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
points[i].sanitize(); // make sure things are in the right range
}
// initialise our DNN structure with the set of points
DynamicNearestNeighbours * DNN;
#ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
bool verbose = false;
bool ignore_nearest_is_mirror = (_Rparam < twopi);
if (_strategy == NlnN4pi) {
DNN = new Dnn4piCylinder(points,verbose);
} else if (_strategy == NlnN3pi) {
DNN = new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose);
} else if (_strategy == NlnN) {
DNN = new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose);
} else
#else
if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
ostringstream err;
err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
err << " supported because FastJet was compiled without CGAL"<<endl;
throw Error(err.str());
//assert(false);
}
#endif // DROP_CGAL
{
ostringstream err;
err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
assert(false);
throw Error(err.str());
}
// We will find nearest neighbour for each vertex, and include
// distance in map (NB DistMap is a typedef given in the .h file)
DistMap DijMap;
// fill the map with the minimal (as far as we know) subset of Dij
// distances (i.e. nearest neighbour ones).
for (int ii = 0; ii < n; ii++) {
_add_ktdistance_to_map(ii, DijMap, DNN);
}
// run the clustering (go up to i=n-1, but then will stop half-way down,
// when we reach that point -- it will be the final beam jet and there
// will be no nearest neighbours to find).
for (int i=0;i<n;i++) {
// find nearest vertices
// NB: skip cases where the point is not there anymore!
TwoVertices SmallestDijPair;
int jet_i, jet_j;
double SmallestDij;
bool Valid2;
bool recombine_with_beam;
do {
SmallestDij = DijMap.begin()->first;
SmallestDijPair = DijMap.begin()->second;
jet_i = SmallestDijPair.first;
jet_j = SmallestDijPair.second;
// distance is immediately removed regardless of whether or not
// it is used.
// Some temporary testing code relating to problems with the gcc-3.2 compiler
//cout << "got here and size is "<< DijMap.size()<< " and it is "<<SmallestDij <<"\n";
//cout << jet_i << " "<< jet_j<<"\n";
DijMap.erase(DijMap.begin());
//cout << "got beyond here\n";
// need to "prime" the validity of jet_j in such a way that
// if it corresponds to the beam then it is automatically valid.
recombine_with_beam = (jet_j == BeamJet);
if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
else {Valid2 = true;}
} while ( !DNN->Valid(jet_i) || !Valid2);
// The following part acts just on jet momenta and on the history.
// The action on the nearest-neighbour structures takes place
// later (only if at least 2 jets are around).
if (! recombine_with_beam) {
int nn; // will be index of new jet
_do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
//OBS // merge the two jets, add new jet, remove old ones
//OBS _jets.push_back(_jets[jet_i] + _jets[jet_j]);
//OBS
//OBS int nn = _jets.size()-1;
//OBS _jets[nn].set_cluster_hist_index(n+i);
//OBS
//OBS // get corresponding indices in history structure
//OBS int hist_i = _jets[jet_i].cluster_hist_index();
//OBS int hist_j = _jets[jet_j].cluster_hist_index();
//OBS
//OBS
//OBS _add_step_to_history(n+i,min(hist_i,hist_j), max(hist_i,hist_j),
//OBS _jets.size()-1, SmallestDij);
// add new point to points vector
EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
newpoint.sanitize(); // make sure it is in correct range
points.push_back(newpoint);
} else {
// recombine the jet with the beam
_do_iB_recombination_step(jet_i, SmallestDij);
//OBS _add_step_to_history(n+i,_jets[jet_i].cluster_hist_index(),BeamJet,
//OBS Invalid, SmallestDij);
}
// exit the loop because we do not want to look for nearest neighbours
// etc. of zero partons
if (i == n-1) {break;}
vector<int> updated_neighbours;
if (! recombine_with_beam) {
// update DNN
int point3;
DNN->RemoveCombinedAddCombination(jet_i, jet_j,
points[points.size()-1], point3,
updated_neighbours);
// C++ beginners' comment: static_cast to unsigned int is necessary
// to do away with warnings about type mismatch between point3 (int)
// and points.size (unsigned int)
if (static_cast<unsigned int> (point3) != points.size()-1) {
throw Error("INTERNAL ERROR: point3 != points.size()-1");}
} else {
// update DNN
DNN->RemovePoint(jet_i, updated_neighbours);
}
// update map
vector<int>::iterator it = updated_neighbours.begin();
for (; it != updated_neighbours.end(); ++it) {
int ii = *it;
_add_ktdistance_to_map(ii, DijMap, DNN);
}
} // end clustering loop
// remember to clean up!
delete DNN;
}
| void ClusterSequence::_do_Cambridge_inclusive_jets | ( | ) | [private] |
Definition at line 336 of file ClusterSequence_CP2DChan.cc.
{
unsigned int n = _history.size();
for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
if (_history[hist_i].child == Invalid) {
_do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
}
}
}
| void ClusterSequence::_do_iB_recombination_step | ( | const int & | jet_i, |
| const double & | diB | ||
| ) | [protected] |
carry out an recombination step in which _jets[jet_i] merges with the beam,
carries out the bookkeeping associated with the step of recombining jet_i with the beam
Definition at line 1116 of file ClusterSequence.cc.
Referenced by plugin_record_iB_recombination().
{
// get history index
int newstep_k = _history.size();
// recombine the jet with the beam
_add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
Invalid, diB);
}
| void ClusterSequence::_do_ij_recombination_step | ( | const int & | jet_i, |
| const int & | jet_j, | ||
| const double & | dij, | ||
| int & | newjet_k | ||
| ) | [protected] |
carry out the recombination between the jets numbered jet_i and jet_j, at distance scale dij; return the index newjet_k of the result of the recombination of i and j.
carries out the bookkeeping associated with the step of recombining jet_i and jet_j (assuming a distance dij) and returns the index of the recombined jet, newjet_k.
Definition at line 1083 of file ClusterSequence.cc.
References d0::inline_maths::min().
Referenced by plugin_record_ij_recombination().
{
// create the new jet by recombining the first two
PseudoJet newjet;
_jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
_jets.push_back(newjet);
// original version...
//_jets.push_back(_jets[jet_i] + _jets[jet_j]);
// get its index
newjet_k = _jets.size()-1;
// get history index
int newstep_k = _history.size();
// and provide jet with the info
_jets[newjet_k].set_cluster_hist_index(newstep_k);
// finally sort out the history
int hist_i = _jets[jet_i].cluster_hist_index();
int hist_j = _jets[jet_j].cluster_hist_index();
_add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
newjet_k, dij);
}
| void ClusterSequence::_extract_tree_children | ( | int | pos, |
| std::valarray< bool > & | , | ||
| const std::valarray< int > & | , | ||
| std::vector< int > & | |||
| ) | const [private] |
internal routine associated with the construction of the unique history order (following children in the tree)
Reimplemented in ClusterSequenceActiveArea.
| void ClusterSequence::_extract_tree_parents | ( | int | pos, |
| std::valarray< bool > & | , | ||
| const std::valarray< int > & | , | ||
| std::vector< int > & | |||
| ) | const [private] |
internal routine associated with the construction of the unique history order (following parents in the tree)
Reimplemented in ClusterSequenceActiveArea.
| void ClusterSequence::_faster_tiled_N2_cluster | ( | ) | [private] |
run a tiled clustering
Definition at line 510 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::TiledJet::_jets_index, ClusterSequence::Tile::begin_tiles, ClusterSequence::TiledJet::diJ_posn, ClusterSequence::Tile::end_tiles, ClusterSequence::Tile::head, ClusterSequence::TiledJet::next, ClusterSequence::TiledJet::NN, ClusterSequence::TiledJet::NN_dist, ClusterSequence::Tile::RH_tiles, ClusterSequence::Tile::tagged, and ClusterSequence::TiledJet::tile_index.
{
_initialise_tiles();
int n = _jets.size();
TiledJet * briefjets = new TiledJet[n];
TiledJet * jetA = briefjets, * jetB;
TiledJet oldB;
// will be used quite deep inside loops, but declare it here so that
// memory (de)allocation gets done only once
vector<int> tile_union(3*n_tile_neighbours);
// initialise the basic jet info
for (int i = 0; i< n; i++) {
_tj_set_jetinfo(jetA, i);
//cout << i<<": "<<jetA->tile_index<<"\n";
jetA++; // move on to next entry of briefjets
}
TiledJet * head = briefjets; // a nicer way of naming start
// set up the initial nearest neighbour information
vector<Tile>::const_iterator tile;
for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
// first do it on this tile
for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
double dist = _bj_dist(jetA,jetB);
if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
}
}
// then do it for RH tiles
for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
double dist = _bj_dist(jetA,jetB);
if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
}
}
}
// no need to do it for LH tiles, since they are implicitly done
// when we set NN for both jetA and jetB on the RH tiles.
}
// now create the diJ (where J is i's NN) table -- remember that
// we differ from standard normalisation here by a factor of R2
// (corrected for at the end).
struct diJ_plus_link {
double diJ; // the distance
TiledJet * jet; // the jet (i) for which we've found this distance
// (whose NN will the J).
};
diJ_plus_link * diJ = new diJ_plus_link[n];
jetA = head;
for (int i = 0; i < n; i++) {
diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
diJ[i].jet = jetA; // our compact diJ table will not be in
jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
// so set up bi-directional correspondence here.
jetA++; // have jetA follow i
}
// now run the recombination loop
int history_location = n-1;
while (n > 0) {
// find the minimum of the diJ on this round
diJ_plus_link * best, *stop; // pointers a bit faster than indices
// could use best to keep track of diJ min, but it turns out to be
// marginally faster to have a separate variable (avoids n
// dereferences at the expense of n/2 assignments).
double diJ_min = diJ[0].diJ; // initialise the best one here.
best = diJ; // and here
stop = diJ+n;
for (diJ_plus_link * here = diJ+1; here != stop; here++) {
if (here->diJ < diJ_min) {best = here; diJ_min = here->diJ;}
}
// do the recombination between A and B
history_location++;
jetA = best->jet;
jetB = jetA->NN;
// put the normalisation back in
diJ_min *= _invR2;
if (jetB != NULL) {
// jet-jet recombination
// If necessary relabel A & B to ensure jetB < jetA, that way if
// the larger of them == newtail then that ends up being jetA and
// the new jet that is added as jetB is inserted in a position that
// has a future!
if (jetA < jetB) {swap(jetA,jetB);}
int nn; // new jet index
_do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
//OBS// get the two history indices
//OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
//OBSint ihstry_b = _jets[jetB->_jets_index].cluster_hist_index();
//OBS// create the recombined jet
//OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
//OBSint nn = _jets.size() - 1;
//OBS_jets[nn].set_cluster_hist_index(history_location);
//OBS// update history
//OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
//OBS_add_step_to_history(history_location,
//OBS min(ihstry_a,ihstry_b),max(ihstry_a,ihstry_b),
//OBS nn, diJ_min);
// what was jetB will now become the new jet
_bj_remove_from_tiles(jetA);
oldB = * jetB; // take a copy because we will need it...
_bj_remove_from_tiles(jetB);
_tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
// (also registers the jet in the tiling)
} else {
// jet-beam recombination
// get the hist_index
_do_iB_recombination_step(jetA->_jets_index, diJ_min);
//OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
//OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
//OBS_add_step_to_history(history_location,ihstry_a,BeamJet,Invalid,diJ_min);
_bj_remove_from_tiles(jetA);
}
// first establish the set of tiles over which we are going to
// have to run searches for updated and new nearest-neighbours --
// basically a combination of vicinity of the tiles of the two old
// and one new jet.
int n_near_tiles = 0;
_add_untagged_neighbours_to_tile_union(jetA->tile_index,
tile_union, n_near_tiles);
if (jetB != NULL) {
if (jetB->tile_index != jetA->tile_index) {
_add_untagged_neighbours_to_tile_union(jetB->tile_index,
tile_union,n_near_tiles);
}
if (oldB.tile_index != jetA->tile_index &&
oldB.tile_index != jetB->tile_index) {
_add_untagged_neighbours_to_tile_union(oldB.tile_index,
tile_union,n_near_tiles);
}
}
// now update our nearest neighbour info and diJ table
// first reduce size of table
n--;
// then compactify the diJ by taking the last of the diJ and copying
// it to the position occupied by the diJ for jetA
diJ[n].jet->diJ_posn = jetA->diJ_posn;
diJ[jetA->diJ_posn] = diJ[n];
// Initialise jetB's NN distance as well as updating it for
// other particles.
// Run over all tiles in our union
for (int itile = 0; itile < n_near_tiles; itile++) {
Tile * tile = &_tiles[tile_union[itile]];
tile->tagged = false; // reset tag, since we're done with unions
// run over all jets in the current tile
for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
// see if jetI had jetA or jetB as a NN -- if so recalculate the NN
if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
jetI->NN_dist = _R2;
jetI->NN = NULL;
// now go over tiles that are neighbours of I (include own tile)
for (Tile ** near_tile = tile->begin_tiles;
near_tile != tile->end_tiles; near_tile++) {
// and then over the contents of that tile
for (TiledJet * jetJ = (*near_tile)->head;
jetJ != NULL; jetJ = jetJ->next) {
double dist = _bj_dist(jetI,jetJ);
if (dist < jetI->NN_dist && jetJ != jetI) {
jetI->NN_dist = dist; jetI->NN = jetJ;
}
}
}
diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
}
// check whether new jetB is closer than jetI's current NN and
// if jetI is closer than jetB's current (evolving) nearest
// neighbour. Where relevant update things
if (jetB != NULL) {
double dist = _bj_dist(jetI,jetB);
if (dist < jetI->NN_dist) {
if (jetI != jetB) {
jetI->NN_dist = dist;
jetI->NN = jetB;
diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
}
}
if (dist < jetB->NN_dist) {
if (jetI != jetB) {
jetB->NN_dist = dist;
jetB->NN = jetI;}
}
}
}
}
// finally, register the updated kt distance for B
if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
}
// final cleaning up;
delete[] diJ;
delete[] briefjets;
}
| void ClusterSequence::_fill_initial_history | ( | ) | [protected] |
fill out the history (and jet cross refs) related to the initial set of jets (assumed already to have been "transferred"), without any clustering
Definition at line 251 of file ClusterSequence.cc.
References ClusterSequence::history_element::child, ClusterSequence::history_element::dij, ClusterSequence::history_element::jetp_index, ClusterSequence::history_element::max_dij_so_far, ClusterSequence::history_element::parent1, and ClusterSequence::history_element::parent2.
{
//if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}
// reserve sufficient space for everything
_jets.reserve(_jets.size()*2);
_history.reserve(_jets.size()*2);
_Qtot = 0;
for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
history_element element;
element.parent1 = InexistentParent;
element.parent2 = InexistentParent;
element.child = Invalid;
element.jetp_index = i;
element.dij = 0.0;
element.max_dij_so_far = 0.0;
_history.push_back(element);
// do any momentum preprocessing needed by the recombination scheme
_jet_def.recombiner()->preprocess(_jets[i]);
// get cross-referencing right from PseudoJets
_jets[i].set_cluster_hist_index(i);
// determine the total energy in the event
_Qtot += _jets[i].E();
}
_initial_n = _jets.size();
}
| void ClusterSequence::_initialise_and_run | ( | const double & | R, |
| const Strategy & | strategy, | ||
| const bool & | writeout_combinations | ||
| ) | [protected] |
This is an alternative routine for initialising and running the clustering, provided for legacy purposes.
The jet finder is that specified in the static member _default_jet_algorithm.
Definition at line 57 of file ClusterSequence.cc.
{
JetDefinition jet_def(_default_jet_algorithm, R, strategy);
_initialise_and_run(jet_def, writeout_combinations);
}
| void ClusterSequence::_initialise_and_run | ( | const JetDefinition & | jet_def, |
| const bool & | writeout_combinations | ||
| ) | [protected] |
This is the routine that will do all the initialisation and then run the clustering (may be called by various constructors).
It assumes _jets contains the momenta to be clustered.
Definition at line 68 of file ClusterSequence.cc.
References antikt_algorithm, Best, cambridge_algorithm, ee_genkt_algorithm, ee_kt_algorithm, JetDefinition::jet_algorithm(), d0::inline_maths::min(), N2MinHeapTiled, N2Plain, N2PoorTiled, N2Tiled, N3Dumb, NlnN, NlnN3pi, NlnN4pi, NlnNCam, NlnNCam2pi2R, NlnNCam4pi, fastjet::pi, and plugin_algorithm.
Referenced by ClusterSequence().
{
// transfer all relevant info into internal variables
_decant_options(jet_def, writeout_combinations);
// set up the history entries for the initial particles (those
// currently in _jets)
_fill_initial_history();
// don't run anything if the event is empty
if (n_particles() == 0) return;
// ----- deal with special cases: plugins & e+e- ------
if (_jet_algorithm == plugin_algorithm) {
// allows plugin_xyz() functions to modify cluster sequence
_plugin_activated = true;
// let the plugin do its work here
_jet_def.plugin()->run_clustering( (*this) );
_plugin_activated = false;
return;
} else if (_jet_algorithm == ee_kt_algorithm ||
_jet_algorithm == ee_genkt_algorithm) {
// ignore requested strategy
_strategy = N2Plain;
if (_jet_algorithm == ee_kt_algorithm) {
// make sure that R is large enough so that "beam" recomb only
// occurs when a single particle is left
// Normally, this should be automatically set to 4 from JetDefinition
assert(_Rparam > 2.0);
// this is used to renormalise the dij to get a "standard" form
// and our convention in e+e- will be different from that
// in long.inv case; NB: _invR2 name should be changed -> _renorm_dij?
_invR2 = 1.0;
} else {
// as of 2009-01-09, choose R to be an angular distance, in
// radians. Since the algorithm uses 2(1-cos(theta)) as its
// squared angular measure, make sure that the _R2 is defined
// in a similar way.
if (_Rparam > pi) {
// choose a value that ensures that back-to-back particles will
// always recombine
//_R2 = 4.0000000000001;
_R2 = 2 * ( 3.0 + cos(_Rparam) );
} else {
_R2 = 2 * ( 1.0 - cos(_Rparam) );
}
_invR2 = 1.0/_R2;
}
//_simple_N2_cluster<EEBriefJet>();
_simple_N2_cluster_EEBriefJet();
return;
}
// automatically redefine the strategy according to N if that is
// what the user requested -- transition points (and especially
// their R-dependence) are based on empirical observations for a
// R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual
// core] with 2MB of cache).
if (_strategy == Best) {
int N = _jets.size();
if (N > 6200/pow(_Rparam,2.0)
&& jet_def.jet_algorithm() == cambridge_algorithm) {
_strategy = NlnNCam;}
else
#ifndef DROP_CGAL
if ((N > 16000/pow(_Rparam,1.15) && jet_def.jet_algorithm() != antikt_algorithm)
|| N > 35000/pow(_Rparam,1.15)) {
_strategy = NlnN; }
else
#endif // DROP_CGAL
if (N > 450) {
_strategy = N2MinHeapTiled;
}
else if (N > 55*max(0.5,min(1.0,_Rparam))) {// empirical scaling with R
_strategy = N2Tiled;
} else {
_strategy = N2Plain;
}
}
// run the code containing the selected strategy
if (_strategy == NlnN || _strategy == NlnN3pi
|| _strategy == NlnN4pi ) {
this->_delaunay_cluster();
} else if (_strategy == N3Dumb ) {
this->_really_dumb_cluster();
} else if (_strategy == N2Tiled) {
this->_faster_tiled_N2_cluster();
} else if (_strategy == N2PoorTiled) {
this->_tiled_N2_cluster();
} else if (_strategy == N2Plain) {
// BriefJet provides standard long.invariant kt alg.
//this->_simple_N2_cluster<BriefJet>();
this->_simple_N2_cluster_BriefJet();
} else if (_strategy == N2MinHeapTiled) {
this->_minheap_faster_tiled_N2_cluster();
} else if (_strategy == NlnNCam4pi) {
this->_CP2DChan_cluster();
} else if (_strategy == NlnNCam2pi2R) {
this->_CP2DChan_cluster_2pi2R();
} else if (_strategy == NlnNCam) {
this->_CP2DChan_cluster_2piMultD();
} else {
ostringstream err;
err << "Unrecognised value for strategy: "<<_strategy;
throw Error(err.str());
//assert(false);
}
}
| void ClusterSequence::_initialise_tiles | ( | ) | [private] |
Set up the tiles:
- decide the range in eta
- allocate the tiles
- set up the cross-referencing info between tiles.
The neighbourhood of a tile is set up as follows
LRR LXR LLR
such that tiles is an array containing XLLLLRRRR with pointers | \ RH_tiles \ surrounding_tiles
with appropriate precautions when close to the edge of the tiled region.
Definition at line 86 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::Tile::begin_tiles, ClusterSequence::Tile::end_tiles, ClusterSequence::Tile::head, ClusterSequence::Tile::RH_tiles, ClusterSequence::Tile::surrounding_tiles, ClusterSequence::Tile::tagged, and twopi.
{
// first decide tile sizes (with a lower bound to avoid huge memory use with
// very small R)
double default_size = max(0.1,_Rparam);
_tile_size_eta = default_size;
_n_tiles_phi = int(floor(twopi/default_size));
_tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
// always include zero rapidity in the tiling region
_tiles_eta_min = 0.0;
_tiles_eta_max = 0.0;
// but go no further than following
const double maxrap = 7.0;
// and find out how much further one should go
for(unsigned int i = 0; i < _jets.size(); i++) {
double eta = _jets[i].rap();
// first check if eta is in range -- to avoid taking into account
// very spurious rapidities due to particles with near-zero kt.
if (abs(eta) < maxrap) {
if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
}
}
// now adjust the values
_tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
_tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
_tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
_tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
// allocate the tiles
_tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
// now set up the cross-referencing between tiles
for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
Tile * tile = & _tiles[_tile_index(ieta,iphi)];
// no jets in this tile yet
tile->head = NULL; // first element of tiles points to itself
tile->begin_tiles[0] = tile;
Tile ** pptile = & (tile->begin_tiles[0]);
pptile++;
//
// set up L's in column to the left of X
tile->surrounding_tiles = pptile;
if (ieta > _tiles_ieta_min) {
// with the itile subroutine, we can safely run tiles from
// idphi=-1 to idphi=+1, because it takes care of
// negative and positive boundaries
for (int idphi = -1; idphi <=+1; idphi++) {
*pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
pptile++;
}
}
// now set up last L (below X)
*pptile = & _tiles[_tile_index(ieta,iphi-1)];
pptile++;
// set up first R (above X)
tile->RH_tiles = pptile;
*pptile = & _tiles[_tile_index(ieta,iphi+1)];
pptile++;
// set up remaining R's, to the right of X
if (ieta < _tiles_ieta_max) {
for (int idphi = -1; idphi <= +1; idphi++) {
*pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
pptile++;
}
}
// now put semaphore for end tile
tile->end_tiles = pptile;
// finally make sure tiles are untagged
tile->tagged = false;
}
}
}
| void ClusterSequence::_minheap_faster_tiled_N2_cluster | ( | ) | [private] |
run a tiled clustering, with our minheap for keeping track of the smallest dij
Definition at line 727 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::TiledJet::_jets_index, ClusterSequence::Tile::begin_tiles, ClusterSequence::Tile::end_tiles, ClusterSequence::Tile::head, ClusterSequence::TiledJet::label_minheap_update_done(), MinHeap::minloc(), MinHeap::minval(), ClusterSequence::TiledJet::next, ClusterSequence::TiledJet::NN, ClusterSequence::TiledJet::NN_dist, MinHeap::remove(), ClusterSequence::Tile::RH_tiles, ClusterSequence::Tile::tagged, ClusterSequence::TiledJet::tile_index, and MinHeap::update().
{
_initialise_tiles();
int n = _jets.size();
TiledJet * briefjets = new TiledJet[n];
TiledJet * jetA = briefjets, * jetB;
TiledJet oldB;
// will be used quite deep inside loops, but declare it here so that
// memory (de)allocation gets done only once
vector<int> tile_union(3*n_tile_neighbours);
// initialise the basic jet info
for (int i = 0; i< n; i++) {
_tj_set_jetinfo(jetA, i);
//cout << i<<": "<<jetA->tile_index<<"\n";
jetA++; // move on to next entry of briefjets
}
TiledJet * head = briefjets; // a nicer way of naming start
// set up the initial nearest neighbour information
vector<Tile>::const_iterator tile;
for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
// first do it on this tile
for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
double dist = _bj_dist(jetA,jetB);
if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
}
}
// then do it for RH tiles
for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
double dist = _bj_dist(jetA,jetB);
if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
}
}
}
// no need to do it for LH tiles, since they are implicitly done
// when we set NN for both jetA and jetB on the RH tiles.
}
//struct diJ_plus_link {
// double diJ; // the distance
// TiledJet * jet; // the jet (i) for which we've found this distance
// // (whose NN will the J).
//};
//diJ_plus_link * diJ = new diJ_plus_link[n];
//jetA = head;
//for (int i = 0; i < n; i++) {
// diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
// diJ[i].jet = jetA; // our compact diJ table will not be in
// jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
// // so set up bi-directional correspondence here.
// jetA++; // have jetA follow i
//}
vector<double> diJs(n);
for (int i = 0; i < n; i++) {
diJs[i] = _bj_diJ(&briefjets[i]);
briefjets[i].label_minheap_update_done();
}
MinHeap minheap(diJs);
// have a stack telling us which jets we'll have to update on the heap
vector<TiledJet *> jets_for_minheap;
jets_for_minheap.reserve(n);
// now run the recombination loop
int history_location = n-1;
while (n > 0) {
double diJ_min = minheap.minval() *_invR2;
jetA = head + minheap.minloc();
// do the recombination between A and B
history_location++;
jetB = jetA->NN;
if (jetB != NULL) {
// jet-jet recombination
// If necessary relabel A & B to ensure jetB < jetA, that way if
// the larger of them == newtail then that ends up being jetA and
// the new jet that is added as jetB is inserted in a position that
// has a future!
if (jetA < jetB) {swap(jetA,jetB);}
int nn; // new jet index
_do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
// what was jetB will now become the new jet
_bj_remove_from_tiles(jetA);
oldB = * jetB; // take a copy because we will need it...
_bj_remove_from_tiles(jetB);
_tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
// (also registers the jet in the tiling)
} else {
// jet-beam recombination
// get the hist_index
_do_iB_recombination_step(jetA->_jets_index, diJ_min);
_bj_remove_from_tiles(jetA);
}
// remove the minheap entry for jetA
minheap.remove(jetA-head);
// first establish the set of tiles over which we are going to
// have to run searches for updated and new nearest-neighbours --
// basically a combination of vicinity of the tiles of the two old
// and one new jet.
int n_near_tiles = 0;
_add_untagged_neighbours_to_tile_union(jetA->tile_index,
tile_union, n_near_tiles);
if (jetB != NULL) {
if (jetB->tile_index != jetA->tile_index) {
_add_untagged_neighbours_to_tile_union(jetB->tile_index,
tile_union,n_near_tiles);
}
if (oldB.tile_index != jetA->tile_index &&
oldB.tile_index != jetB->tile_index) {
_add_untagged_neighbours_to_tile_union(oldB.tile_index,
tile_union,n_near_tiles);
}
// indicate that we'll have to update jetB in the minheap
jetB->label_minheap_update_needed();
jets_for_minheap.push_back(jetB);
}
// Initialise jetB's NN distance as well as updating it for
// other particles.
// Run over all tiles in our union
for (int itile = 0; itile < n_near_tiles; itile++) {
Tile * tile = &_tiles[tile_union[itile]];
tile->tagged = false; // reset tag, since we're done with unions
// run over all jets in the current tile
for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
// see if jetI had jetA or jetB as a NN -- if so recalculate the NN
if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
jetI->NN_dist = _R2;
jetI->NN = NULL;
// label jetI as needing heap action...
if (!jetI->minheap_update_needed()) {
jetI->label_minheap_update_needed();
jets_for_minheap.push_back(jetI);}
// now go over tiles that are neighbours of I (include own tile)
for (Tile ** near_tile = tile->begin_tiles;
near_tile != tile->end_tiles; near_tile++) {
// and then over the contents of that tile
for (TiledJet * jetJ = (*near_tile)->head;
jetJ != NULL; jetJ = jetJ->next) {
double dist = _bj_dist(jetI,jetJ);
if (dist < jetI->NN_dist && jetJ != jetI) {
jetI->NN_dist = dist; jetI->NN = jetJ;
}
}
}
}
// check whether new jetB is closer than jetI's current NN and
// if jetI is closer than jetB's current (evolving) nearest
// neighbour. Where relevant update things
if (jetB != NULL) {
double dist = _bj_dist(jetI,jetB);
if (dist < jetI->NN_dist) {
if (jetI != jetB) {
jetI->NN_dist = dist;
jetI->NN = jetB;
// label jetI as needing heap action...
if (!jetI->minheap_update_needed()) {
jetI->label_minheap_update_needed();
jets_for_minheap.push_back(jetI);}
}
}
if (dist < jetB->NN_dist) {
if (jetI != jetB) {
jetB->NN_dist = dist;
jetB->NN = jetI;}
}
}
}
}
// deal with jets whose minheap entry needs updating
while (jets_for_minheap.size() > 0) {
TiledJet * jetI = jets_for_minheap.back();
jets_for_minheap.pop_back();
minheap.update(jetI-head, _bj_diJ(jetI));
jetI->label_minheap_update_done();
}
n--;
}
// final cleaning up;
delete[] briefjets;
}
| bool ClusterSequence::_potentially_valid | ( | const PseudoJet & | jet | ) | const [inline, protected] |
returns true if the jet has a history index contained within the range of this CS
Definition at line 487 of file ClusterSequence.hh.
References _history.
{
return jet.cluster_hist_index() >= 0
&& jet.cluster_hist_index() < int(_history.size());
}
| void ClusterSequence::_print_banner | ( | ) | [private] |
for making sure the user knows what it is they're running...
Definition at line 197 of file ClusterSequence.cc.
References fastjet_version.
{
if (!_first_time) {return;}
_first_time = false;
//Symp. Discr. Alg, p.472 (2002) and CGAL (http://www.cgal.org);
cout << "#--------------------------------------------------------------------------\n";
cout << "# FastJet release " << fastjet_version << endl;
cout << "# Written by M. Cacciari, G.P. Salam and G. Soyez \n";
cout << "# http://www.fastjet.fr \n";
cout << "# \n";
cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen \n";
cout << "# clustering using fast geometric algorithms, with jet areas and optional\n";
cout << "# external jet-finder plugins. If you use this code towards a scientific \n";
cout << "# publication please cite Phys. Lett. B641 (2006) [hep-ph/0512210] and \n";
cout << "# M. Cacciari, G.P. Salam and G. Soyez, http://fastjet.fr/ \n";
cout << "# \n";
cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code" ;
#ifndef DROP_CGAL
cout << endl << "# and CGAL: http://www.cgal.org/";
#endif // DROP_CGAL
cout << ".\n";
cout << "#-------------------------------------------------------------------------\n";
// make sure we really have the output done.
cout.flush();
}
| void ClusterSequence::_print_tiles | ( | TiledJet * | briefjets | ) | const [private] |
output the contents of the tiles
Definition at line 212 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::TiledJet::next.
{
for (vector<Tile>::const_iterator tile = _tiles.begin();
tile < _tiles.end(); tile++) {
cout << "Tile " << tile - _tiles.begin()<<" = ";
vector<int> list;
for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
list.push_back(jetI-briefjets);
//cout <<" "<<jetI-briefjets;
}
sort(list.begin(),list.end());
for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
cout <<"\n";
}
}
| void ClusterSequence::_really_dumb_cluster | ( | ) | [private] |
Run the clustering in a very slow variant of the N^3 algorithm.
The only thing this routine has going for it is that memory usage is O(N)!
Definition at line 50 of file ClusterSequence_DumbN3.cc.
References d0::inline_maths::min(), and d0::inline_maths::y().
{
// the array that will be overwritten here will be one
// of pointers to jets.
vector<PseudoJet *> jetsp(_jets.size());
vector<int> indices(_jets.size());
for (size_t i = 0; i<_jets.size(); i++) {
jetsp[i] = & _jets[i];
indices[i] = i;
}
for (int n = jetsp.size(); n > 0; n--) {
int ii, jj;
// find smallest beam distance [remember jet_scale_for_algorithm
// will return kt^2 for the kt-algorithm and 1 for the Cambridge/Aachen]
double ymin = jet_scale_for_algorithm(*(jetsp[0]));
ii = 0; jj = -2;
for (int i = 0; i < n; i++) {
double yiB = jet_scale_for_algorithm(*(jetsp[i]));
if (yiB < ymin) {
ymin = yiB; ii = i; jj = -2;}
}
// find smallest distance between pair of jetsp
for (int i = 0; i < n-1; i++) {
for (int j = i+1; j < n; j++) {
//double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2;
double y = min(jet_scale_for_algorithm(*(jetsp[i])),
jet_scale_for_algorithm(*(jetsp[j])))
* jetsp[i]->plain_distance(*jetsp[j])*_invR2;
if (y < ymin) {ymin = y; ii = i; jj = j;}
}
}
// output recombination sequence
// old "ktclus" way of labelling
//cout <<n<< " "<< ii+1 << " with " << jj+1 << "; y = "<< ymin<<endl;
// new delaunay way of labelling
int jjindex_or_beam, iiindex;
if (jj < 0) {jjindex_or_beam = BeamJet; iiindex = indices[ii];}
else {
jjindex_or_beam = max(indices[ii],indices[jj]);
iiindex = min(indices[ii],indices[jj]);
}
// now recombine
int newn = 2*jetsp.size() - n;
if (jj >= 0) {
// combine pair
int nn; // new jet index
_do_ij_recombination_step(jetsp[ii]-&_jets[0],
jetsp[jj]-&_jets[0], ymin, nn);
// sort out internal bookkeeping
jetsp[ii] = &_jets[nn];
// have jj point to jet that was pointed at by n-1
// (since original jj is no longer current, so put n-1 into jj)
jetsp[jj] = jetsp[n-1];
indices[ii] = newn;
indices[jj] = indices[n-1];
//OBS_jets.push_back(*jetsp[ii] + *jetsp[jj]);
//OBSjetsp[ii] = &_jets[_jets.size()-1];
//OBS// have jj point to jet that was pointed at by n-1
//OBS// (since original jj is no longer current, so put n-1 into jj)
//OBSjetsp[jj] = jetsp[n-1];
//OBS
//OBSindices[ii] = newn;
//OBSindices[jj] = indices[n-1];
//OBS_add_step_to_history(newn,iiindex,
//OBS jjindex_or_beam,_jets.size()-1,ymin);
} else {
// combine ii with beam
_do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin);
// put last jet (pointer) in place of ii (which has disappeared)
jetsp[ii] = jetsp[n-1];
indices[ii] = indices[n-1];
//OBS_add_step_to_history(newn,iiindex,jjindex_or_beam,Invalid, ymin);
}
}
}
| void ClusterSequence::_simple_N2_cluster | ( | ) | [private] |
| void ClusterSequence::_simple_N2_cluster_BriefJet | ( | ) | [private] |
to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?)
to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
Definition at line 115 of file ClusterSequence_N2.cc.
{
_simple_N2_cluster<BriefJet>();
}
| void ClusterSequence::_simple_N2_cluster_EEBriefJet | ( | ) | [private] |
to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
Definition at line 121 of file ClusterSequence_N2.cc.
{
_simple_N2_cluster<EEBriefJet>();
}
| int ClusterSequence::_tile_index | ( | int | ieta, |
| int | iphi | ||
| ) | const [inline, private] |
Definition at line 728 of file ClusterSequence.hh.
References _n_tiles_phi, and _tiles_ieta_min.
{
// note that (-1)%n = -1 so that we have to add _n_tiles_phi
// before performing modulo operation
return (ieta-_tiles_ieta_min)*_n_tiles_phi
+ (iphi+_n_tiles_phi) % _n_tiles_phi;
}
| int ClusterSequence::_tile_index | ( | const double & | eta, |
| const double & | phi | ||
| ) | const [private] |
return the tile index corresponding to the given eta,phi point
Definition at line 168 of file ClusterSequence_TiledN2.cc.
References twopi.
{
int ieta, iphi;
if (eta <= _tiles_eta_min) {ieta = 0;}
else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
else {
//ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
// following needed in case of rare but nasty rounding errors
if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
ieta = _tiles_ieta_max-_tiles_ieta_min;}
}
// allow for some extent of being beyond range in calculation of phi
// as well
//iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
// with just int and no floor, things run faster but beware
iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
return (iphi + ieta * _n_tiles_phi);
}
| void ClusterSequence::_tiled_N2_cluster | ( | ) | [private] |
run a tiled clustering
Definition at line 267 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::TiledJet::_jets_index, ClusterSequence::Tile::begin_tiles, ClusterSequence::Tile::end_tiles, ClusterSequence::Tile::head, ClusterSequence::TiledJet::next, ClusterSequence::TiledJet::NN, ClusterSequence::TiledJet::NN_dist, ClusterSequence::TiledJet::previous, ClusterSequence::Tile::RH_tiles, and ClusterSequence::TiledJet::tile_index.
{
_initialise_tiles();
int n = _jets.size();
TiledJet * briefjets = new TiledJet[n];
TiledJet * jetA = briefjets, * jetB;
TiledJet oldB;
// will be used quite deep inside loops, but declare it here so that
// memory (de)allocation gets done only once
vector<int> tile_union(3*n_tile_neighbours);
// initialise the basic jet info
for (int i = 0; i< n; i++) {
_tj_set_jetinfo(jetA, i);
//cout << i<<": "<<jetA->tile_index<<"\n";
jetA++; // move on to next entry of briefjets
}
TiledJet * tail = jetA; // a semaphore for the end of briefjets
TiledJet * head = briefjets; // a nicer way of naming start
// set up the initial nearest neighbour information
vector<Tile>::const_iterator tile;
for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
// first do it on this tile
for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
double dist = _bj_dist(jetA,jetB);
if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
}
}
// then do it for RH tiles
for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
double dist = _bj_dist(jetA,jetB);
if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
}
}
}
}
// now create the diJ (where J is i's NN) table -- remember that
// we differ from standard normalisation here by a factor of R2
double * diJ = new double[n];
jetA = head;
for (int i = 0; i < n; i++) {
diJ[i] = _bj_diJ(jetA);
jetA++; // have jetA follow i
}
// now run the recombination loop
int history_location = n-1;
while (tail != head) {
// find the minimum of the diJ on this round
double diJ_min = diJ[0];
int diJ_min_jet = 0;
for (int i = 1; i < n; i++) {
if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min = diJ[i];}
}
// do the recombination between A and B
history_location++;
jetA = & briefjets[diJ_min_jet];
jetB = jetA->NN;
// put the normalisation back in
diJ_min *= _invR2;
//if (n == 19) {cout << "Hello "<<jetA-head<<" "<<jetB-head<<"\n";}
//cout <<" WILL RECOMBINE "<< jetA-briefjets<<" "<<jetB-briefjets<<"\n";
if (jetB != NULL) {
// jet-jet recombination
// If necessary relabel A & B to ensure jetB < jetA, that way if
// the larger of them == newtail then that ends up being jetA and
// the new jet that is added as jetB is inserted in a position that
// has a future!
if (jetA < jetB) {swap(jetA,jetB);}
int nn; // new jet index
_do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
//OBS// get the two history indices
//OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
//OBSint hist_b = _jets[jetB->_jets_index].cluster_hist_index();
//OBS// create the recombined jet
//OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
//OBSint nn = _jets.size() - 1;
//OBS_jets[nn].set_cluster_hist_index(history_location);
//OBS// update history
//OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
//OBS_add_step_to_history(history_location,
//OBS min(hist_a,hist_b),max(hist_a,hist_b),
//OBS nn, diJ_min);
// what was jetB will now become the new jet
_bj_remove_from_tiles(jetA);
oldB = * jetB; // take a copy because we will need it...
_bj_remove_from_tiles(jetB);
_tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
} else {
// jet-beam recombination
_do_iB_recombination_step(jetA->_jets_index, diJ_min);
//OBS// get the hist_index
//OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
//OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
//OBS_add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min);
_bj_remove_from_tiles(jetA);
}
// first establish the set of tiles over which we are going to
// have to run searches for updated and new nearest-neighbours
int n_near_tiles = 0;
_add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
if (jetB != NULL) {
bool sort_it = false;
if (jetB->tile_index != jetA->tile_index) {
sort_it = true;
_add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
}
if (oldB.tile_index != jetA->tile_index &&
oldB.tile_index != jetB->tile_index) {
sort_it = true;
_add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
}
if (sort_it) {
// sort the tiles before then compressing the list
sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
// and now condense the list
int nnn = 1;
for (int i = 1; i < n_near_tiles; i++) {
if (tile_union[i] != tile_union[nnn-1]) {
tile_union[nnn] = tile_union[i];
nnn++;
}
}
n_near_tiles = nnn;
}
}
// now update our nearest neighbour info and diJ table
// first reduce size of table
tail--; n--;
if (jetA == tail) {
// there is nothing to be done
} else {
// Copy last jet contents and diJ info into position of jetA
*jetA = *tail;
diJ[jetA - head] = diJ[tail-head];
// IN the tiling fix pointers to tail and turn them into
// pointers to jetA (from predecessors, successors and the tile
// head if need be)
if (jetA->previous == NULL) {
_tiles[jetA->tile_index].head = jetA;
} else {
jetA->previous->next = jetA;
}
if (jetA->next != NULL) {jetA->next->previous = jetA;}
}
// Initialise jetB's NN distance as well as updating it for
// other particles.
for (int itile = 0; itile < n_near_tiles; itile++) {
Tile * tile = &_tiles[tile_union[itile]];
for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
// see if jetI had jetA or jetB as a NN -- if so recalculate the NN
if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
jetI->NN_dist = _R2;
jetI->NN = NULL;
// now go over tiles that are neighbours of I (include own tile)
for (Tile ** near_tile = tile->begin_tiles;
near_tile != tile->end_tiles; near_tile++) {
// and then over the contents of that tile
for (TiledJet * jetJ = (*near_tile)->head;
jetJ != NULL; jetJ = jetJ->next) {
double dist = _bj_dist(jetI,jetJ);
if (dist < jetI->NN_dist && jetJ != jetI) {
jetI->NN_dist = dist; jetI->NN = jetJ;
}
}
}
diJ[jetI-head] = _bj_diJ(jetI); // update diJ
}
// check whether new jetB is closer than jetI's current NN and
// if need to update things
if (jetB != NULL) {
double dist = _bj_dist(jetI,jetB);
if (dist < jetI->NN_dist) {
if (jetI != jetB) {
jetI->NN_dist = dist;
jetI->NN = jetB;
diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
}
}
if (dist < jetB->NN_dist) {
if (jetI != jetB) {
jetB->NN_dist = dist;
jetB->NN = jetI;}
}
}
}
}
if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
//cout << n<<" "<<briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
// remember to update pointers to tail
for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles;
near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
// and then the contents of that tile
for (TiledJet * jetJ = (*near_tile)->head;
jetJ != NULL; jetJ = jetJ->next) {
if (jetJ->NN == tail) {jetJ->NN = jetA;}
}
}
//for (int i = 0; i < n; i++) {
// if (briefjets[i].NN-briefjets >= n && briefjets[i].NN != NULL) {cout <<"YOU MUST BE CRAZY for n ="<<n<<", i = "<<i<<", NN = "<<briefjets[i].NN-briefjets<<"\n";}
//}
if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
//cout << briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
}
// final cleaning up;
delete[] diJ;
delete[] briefjets;
}
| void ClusterSequence::_tj_set_jetinfo | ( | TiledJet *const | jet, |
| const int | _jets_index | ||
| ) | [inline, private] |
Definition at line 191 of file ClusterSequence_TiledN2.cc.
References ClusterSequence::Tile::head, ClusterSequence::TiledJet::next, and ClusterSequence::TiledJet::previous.
{
// first call the generic setup
_bj_set_jetinfo<>(jet, _jets_index);
// Then do the setup specific to the tiled case.
// Find out which tile it belonds to
jet->tile_index = _tile_index(jet->eta, jet->phi);
// Insert it into the tile's linked list of jets
Tile * tile = &_tiles[jet->tile_index];
jet->previous = NULL;
jet->next = tile->head;
if (jet->next != NULL) {jet->next->previous = jet;}
tile->head = jet;
}
| void ClusterSequence::_transfer_input_jets | ( | const std::vector< L > & | pseudojets | ) | [protected] |
transfer the vector<L> of input jets into our own vector<PseudoJet> _jets (with some reserved space for future growth).
Definition at line 777 of file ClusterSequence.hh.
References _jets.
Referenced by ClusterSequence().
{
// this will ensure that we can point to jets without difficulties
// arising.
_jets.reserve(pseudojets.size()*2);
// insert initial jets this way so that any type L that can be
// converted to a pseudojet will work fine (basically PseudoJet
// and any type that has [] subscript access to the momentum
// components, such as CLHEP HepLorentzVector).
for (unsigned int i = 0; i < pseudojets.size(); i++) {
_jets.push_back(pseudojets[i]);}
}
| void ClusterSequence::add_constituents | ( | const PseudoJet & | jet, |
| std::vector< PseudoJet > & | subjet_vector | ||
| ) | const |
add on to subjet_vector the constituents of jet (for internal use mainly)
Definition at line 909 of file ClusterSequence.cc.
{
// find out position in cluster history
int i = jet.cluster_hist_index();
int parent1 = _history[i].parent1;
int parent2 = _history[i].parent2;
if (parent1 == InexistentParent) {
// It is an original particle (labelled by its parent having value
// InexistentParent), therefore add it on to the subjet vector
// Note: we add the initial particle and not simply 'jet' so that
// calling add_constituents with a subtracted jet containing
// only one particle will work.
subjet_vector.push_back(_jets[i]);
return;
}
// add parent 1
add_constituents(_jets[_history[parent1].jetp_index], subjet_vector);
// see if parent2 is a real jet; if it is then add its constituents
if (parent2 != BeamJet) {
add_constituents(_jets[_history[parent2].jetp_index], subjet_vector);
}
}
| vector< PseudoJet > ClusterSequence::constituents | ( | const PseudoJet & | jet | ) | const |
return a vector of the particles that make up jet
Definition at line 818 of file ClusterSequence.cc.
{
vector<PseudoJet> subjets;
add_constituents(jet, subjets);
return subjets;
}
| double ClusterSequence::exclusive_dmerge | ( | const int & | njets | ) | const |
return the dmin corresponding to the recombination that went from n+1 to n jets (sometimes known as d_{n n+1}).
return the dmin corresponding to the recombination that went from n+1 to n jets
Definition at line 538 of file ClusterSequence.cc.
Referenced by exclusive_ymerge().
{
assert(njets >= 0);
if (njets >= _initial_n) {return 0.0;}
return _history[2*_initial_n-njets-1].dij;
}
| double ClusterSequence::exclusive_dmerge_max | ( | const int & | njets | ) | const |
return the maximum of the dmin encountered during all recombinations up to the one that led to an n-jet final state; identical to exclusive_dmerge, except in cases where the dmin do not increase monotonically.
Definition at line 550 of file ClusterSequence.cc.
Referenced by exclusive_ymerge_max().
{
assert(njets >= 0);
if (njets >= _initial_n) {return 0.0;}
return _history[2*_initial_n-njets-1].max_dij_so_far;
}
| vector< PseudoJet > ClusterSequence::exclusive_jets | ( | const double & | dcut | ) | const |
return a vector of all jets (in the sense of the exclusive algorithm) that would be obtained when running the algorithm with the given dcut.
Definition at line 458 of file ClusterSequence.cc.
Referenced by exclusive_jets_ycut().
{
int njets = n_exclusive_jets(dcut);
return exclusive_jets(njets);
}
| vector< PseudoJet > ClusterSequence::exclusive_jets | ( | const int & | njets | ) | const |
return a vector of all jets when the event is clustered (in the exclusive sense) to exactly njets.
Definition at line 466 of file ClusterSequence.cc.
References cambridge_algorithm, ee_genkt_algorithm, ee_kt_algorithm, genkt_algorithm, kt_algorithm, and plugin_algorithm.
{
// make sure the user does not ask for more than jets than there
// were particles in the first place.
assert (njets <= _initial_n);
// provide a warning when extracting exclusive jets for algorithms
// that does not support it explicitly.
// Native algorithm that support it are: kt, ee_kt, cambridge,
// genkt and ee_genkt (both with p>=0)
// For plugins, we check Plugin::exclusive_sequence_meaningful()
if (( _jet_def.jet_algorithm() != kt_algorithm) &&
( _jet_def.jet_algorithm() != cambridge_algorithm) &&
( _jet_def.jet_algorithm() != ee_kt_algorithm) &&
(((_jet_def.jet_algorithm() != genkt_algorithm) &&
(_jet_def.jet_algorithm() != ee_genkt_algorithm)) ||
(_jet_def.extra_param() <0)) &&
((_jet_def.jet_algorithm() != plugin_algorithm) ||
(!_jet_def.plugin()->exclusive_sequence_meaningful())) &&
(_n_exclusive_warnings < 5)) {
_n_exclusive_warnings++;
cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
}
// calculate the point where we have to stop the clustering.
// relation between stop_point, njets assumes one extra jet disappears
// at each clustering.
int stop_point = 2*_initial_n - njets;
// some sanity checking to make sure that e+e- does not give us
// surprises (should we ever implement e+e-)...
if (2*_initial_n != static_cast<int>(_history.size())) {
ostringstream err;
err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
throw Error(err.str());
//assert(false);
}
// now go forwards and reconstitute the jets that we have --
// basically for any history element, see if the parent jets to
// which it refers were created before the stopping point -- if they
// were then add them to the list, otherwise they are subsequent
// recombinations of the jets that we are looking for.
vector<PseudoJet> jets;
for (unsigned int i = stop_point; i < _history.size(); i++) {
int parent1 = _history[i].parent1;
if (parent1 < stop_point) {
jets.push_back(_jets[_history[parent1].jetp_index]);
}
int parent2 = _history[i].parent2;
if (parent2 < stop_point && parent2 > 0) {
jets.push_back(_jets[_history[parent2].jetp_index]);
}
}
// sanity check...
if (static_cast<int>(jets.size()) != njets) {
ostringstream err;
err << "ClusterSequence::exclusive_jets: size of returned vector ("
<<jets.size()<<") does not coincide with requested number of jets ("
<<njets<<")";
throw Error(err.str());
}
return jets;
}
| std::vector<PseudoJet> ClusterSequence::exclusive_jets_ycut | ( | double | ycut | ) | const [inline] |
the exclusive jets obtained at the given ycut
Definition at line 145 of file ClusterSequence.hh.
References exclusive_jets(), and n_exclusive_jets_ycut().
{
int njets = n_exclusive_jets_ycut(ycut);
return exclusive_jets(njets);
}
| double ClusterSequence::exclusive_subdmerge | ( | const PseudoJet & | jet, |
| int | nsub | ||
| ) | const |
return the dij that was present in the merging nsub+1 -> nsub subjets inside this jet.
If the jet has nsub or fewer constituents, it will return 0.
will be zero if nconst <= nsub, since highest will be an original particle have zero dij
Definition at line 622 of file ClusterSequence.cc.
{
set<const history_element*> subhist;
// get the set of history elements that correspond to subjets at
// scale dcut
get_subhist_set(subhist, jet, -1.0, nsub);
set<const history_element*>::iterator highest = subhist.end();
highest--;
return (*highest)->dij;
}
| double ClusterSequence::exclusive_subdmerge_max | ( | const PseudoJet & | jet, |
| int | nsub | ||
| ) | const |
return the maximum dij that occurred in the whole event at the stage that the nsub+1 -> nsub merge of subjets occurred inside this jet.
If the jet has nsub or fewer constituents, it will return 0.
will be zero if nconst <= nsub, since highest will be an original particle have zero dij
Definition at line 643 of file ClusterSequence.cc.
{
set<const history_element*> subhist;
// get the set of history elements that correspond to subjets at
// scale dcut
get_subhist_set(subhist, jet, -1.0, nsub);
set<const history_element*>::iterator highest = subhist.end();
highest--;
return (*highest)->max_dij_so_far;
}
| std::vector< PseudoJet > ClusterSequence::exclusive_subjets | ( | const PseudoJet & | jet, |
| const double & | dcut | ||
| ) | const |
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that would be obtained when running the algorithm with the given dcut.
Time taken is O(m ln m), where m is the number of subjets that are found. If m gets to be of order of the total number of constituents in the jet, this could be substantially slower than just getting that list of constituents.
Definition at line 562 of file ClusterSequence.cc.
{
set<const history_element*> subhist;
// get the set of history elements that correspond to subjets at
// scale dcut
get_subhist_set(subhist, jet, dcut, 0);
// now transfer this into a sequence of jets
vector<PseudoJet> subjets;
subjets.reserve(subhist.size());
for (set<const history_element*>::iterator elem = subhist.begin();
elem != subhist.end(); elem++) {
subjets.push_back(_jets[(*elem)->jetp_index]);
}
return subjets;
}
| std::vector< PseudoJet > ClusterSequence::exclusive_subjets | ( | const PseudoJet & | jet, |
| int | n | ||
| ) | const |
return the list of subjets obtained by unclustering the supplied jet down to n subjets (or all constituents if there are fewer than n).
requires n ln n time
Definition at line 598 of file ClusterSequence.cc.
{
set<const history_element*> subhist;
// get the set of history elements that correspond to subjets at
// scale dcut
get_subhist_set(subhist, jet, -1.0, n);
// now transfer this into a sequence of jets
vector<PseudoJet> subjets;
subjets.reserve(subhist.size());
for (set<const history_element*>::iterator elem = subhist.begin();
elem != subhist.end(); elem++) {
subjets.push_back(_jets[(*elem)->jetp_index]);
}
return subjets;
}
| double ClusterSequence::exclusive_ymerge | ( | int | njets | ) | const [inline] |
return the ymin corresponding to the recombination that went from n+1 to n jets (sometimes known as y_{n n+1}).
Definition at line 136 of file ClusterSequence.hh.
References exclusive_dmerge(), and Q2().
{return exclusive_dmerge(njets) / Q2();}
| double ClusterSequence::exclusive_ymerge_max | ( | int | njets | ) | const [inline] |
same as exclusive_dmerge_max, but normalised to squared total energy
Definition at line 139 of file ClusterSequence.hh.
References exclusive_dmerge_max(), and Q2().
{return exclusive_dmerge_max(njets)/Q2();}
| const Extras* ClusterSequence::extras | ( | ) | const [inline] |
returns a pointer to the extras object (may be null)
Definition at line 320 of file ClusterSequence.hh.
References _extras.
{return _extras.get();}
| void ClusterSequence::get_subhist_set | ( | std::set< const history_element * > & | subhist, |
| const PseudoJet & | jet, | ||
| double | dcut, | ||
| int | maxjet | ||
| ) | const [protected] |
set subhist to be a set pointers to history entries corresponding to the subjets of this jet; one stops going working down through the subjets either when
- there is no further to go
- one has found maxjet entries
- max_dij_so_far <= dcut By setting maxjet=0 one can use just dcut; by setting dcut<0 one can use jet maxjet
return a set of pointers to history entries corresponding to the subjets of this jet; one stops going working down through the subjets either when
- there is no further to go
- one has found maxjet entries
- max_dij_so_far <= dcut
Definition at line 667 of file ClusterSequence.cc.
References ClusterSequence::history_element::max_dij_so_far, ClusterSequence::history_element::parent1, and ClusterSequence::history_element::parent2.
{
subhist.clear();
subhist.insert(&(_history[jet.cluster_hist_index()]));
// establish the set of jets that are relevant
int njet = 1;
while (true) {
// first find out if we need to probe deeper into jet.
// Get history element closest to end of sequence
set<const history_element*>::iterator highest = subhist.end();
assert (highest != subhist.begin());
highest--;
const history_element* elem = *highest;
// make sure we haven't got too many jets
if (njet == maxjet) break;
// make sure it has parents
if (elem->parent1 < 0) break;
// make sure that we still resolve it at scale dcut
if (elem->max_dij_so_far <= dcut) break;
// then do so: replace "highest" with its two parents
subhist.erase(highest);
subhist.insert(&(_history[elem->parent1]));
subhist.insert(&(_history[elem->parent2]));
njet++;
}
}
| bool ClusterSequence::has_child | ( | const PseudoJet & | jet, |
| const PseudoJet *& | childp | ||
| ) | const |
Version of has_child that sets a pointer to the child if the child exists;.
Definition at line 770 of file ClusterSequence.cc.
References ClusterSequence::history_element::child.
{
const history_element & hist = _history[jet.cluster_hist_index()];
// check that this jet has a child and that the child corresponds to
// a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
// behaviour if the child is the same jet but made inclusive...?]
if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
childp = &(_jets[_history[hist.child].jetp_index]);
return true;
} else {
childp = NULL;
return false;
}
}
| bool ClusterSequence::has_child | ( | const PseudoJet & | jet, |
| PseudoJet & | child | ||
| ) | const |
if the jet has a child then return true and give the child jet otherwise return false and set the child to zero
Definition at line 748 of file ClusterSequence.cc.
{
//const history_element & hist = _history[jet.cluster_hist_index()];
//
//if (hist.child >= 0) {
// child = _jets[_history[hist.child].jetp_index];
// return true;
//} else {
// child = PseudoJet(0.0,0.0,0.0,0.0);
// return false;
//}
const PseudoJet * childp;
bool res = has_child(jet, childp);
if (res) {
child = *childp;
return true;
} else {
child = PseudoJet(0.0,0.0,0.0,0.0);
return false;
}
}
| bool ClusterSequence::has_parents | ( | const PseudoJet & | jet, |
| PseudoJet & | parent1, | ||
| PseudoJet & | parent2 | ||
| ) | const |
if the jet has parents in the clustering, it returns true and sets parent1 and parent2 equal to them.
if it has no parents it returns false and sets parent1 and parent2 to zero
Definition at line 722 of file ClusterSequence.cc.
References ClusterSequence::history_element::parent1, and ClusterSequence::history_element::parent2.
{
const history_element & hist = _history[jet.cluster_hist_index()];
// make sure we do not run into any unexpected situations --
// i.e. both parents valid, or neither
assert ((hist.parent1 >= 0 && hist.parent2 >= 0) ||
(hist.parent1 < 0 && hist.parent2 < 0));
if (hist.parent1 < 0) {
parent1 = PseudoJet(0.0,0.0,0.0,0.0);
parent2 = parent1;
return false;
} else {
parent1 = _jets[_history[hist.parent1].jetp_index];
parent2 = _jets[_history[hist.parent2].jetp_index];
// order the parents in decreasing pt
if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
return true;
}
}
| bool ClusterSequence::has_partner | ( | const PseudoJet & | jet, |
| PseudoJet & | partner | ||
| ) | const |
if this jet has a child (and so a partner) return true and give the partner, otherwise return false and set the partner to zero
Definition at line 791 of file ClusterSequence.cc.
References ClusterSequence::history_element::child, ClusterSequence::history_element::parent1, and ClusterSequence::history_element::parent2.
{
const history_element & hist = _history[jet.cluster_hist_index()];
// make sure we have a child and that the child does not correspond
// to a clustering with the beam (or some other invalid quantity)
if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
const history_element & child_hist = _history[hist.child];
if (child_hist.parent1 == jet.cluster_hist_index()) {
// partner will be child's parent2 -- for iB clustering
// parent2 will not be valid
partner = _jets[_history[child_hist.parent2].jetp_index];
} else {
// partner will be child's parent1
partner = _jets[_history[child_hist.parent1].jetp_index];
}
return true;
} else {
partner = PseudoJet(0.0,0.0,0.0,0.0);
return false;
}
}
| const std::vector< ClusterSequence::history_element > & ClusterSequence::history | ( | ) | const [inline] |
allow the user to access the raw internal history.
This is present (as for jets()) in part so that protected derived classes can access this information about other ClusterSequences.
A user who wishes to follow the details of the ClusterSequence can also make use of this information (and should consult the history_element documentation for more information), but should be aware that these internal structures may evolve in future FastJet versions.
Definition at line 830 of file ClusterSequence.hh.
References _history.
Referenced by ClusterSequenceActiveArea::_transfer_areas(), ClusterSequenceActiveArea::_transfer_ghost_free_history(), and NestedDefsPlugin::run_clustering().
{
return _history;
}
| vector< PseudoJet > ClusterSequence::inclusive_jets | ( | const double & | ptmin = 0.0 | ) | const |
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
Time taken should be of the order of the number of jets returned.
Definition at line 387 of file ClusterSequence.cc.
References antikt_algorithm, cambridge_algorithm, cambridge_for_passive_algorithm, ee_genkt_algorithm, ee_kt_algorithm, genkt_algorithm, kt_algorithm, and plugin_algorithm.
{
double dcut = ptmin*ptmin;
int i = _history.size() - 1; // last jet
vector<PseudoJet> jets;
if (_jet_algorithm == kt_algorithm) {
while (i >= 0) {
// with our specific definition of dij and diB (i.e. R appears only in
// dij), then dij==diB is the same as the jet.perp2() and we can exploit
// this in selecting the jets...
if (_history[i].max_dij_so_far < dcut) {break;}
if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
// for beam jets
int parent1 = _history[i].parent1;
jets.push_back(_jets[_history[parent1].jetp_index]);}
i--;
}
} else if (_jet_algorithm == cambridge_algorithm) {
while (i >= 0) {
// inclusive jets are all at end of clustering sequence in the
// Cambridge algorithm -- so if we find a non-exclusive jet, then
// we can exit
if (_history[i].parent2 != BeamJet) {break;}
int parent1 = _history[i].parent1;
const PseudoJet & jet = _jets[_history[parent1].jetp_index];
if (jet.perp2() >= dcut) {jets.push_back(jet);}
i--;
}
} else if (_jet_algorithm == plugin_algorithm
|| _jet_algorithm == ee_kt_algorithm
|| _jet_algorithm == antikt_algorithm
|| _jet_algorithm == genkt_algorithm
|| _jet_algorithm == ee_genkt_algorithm
|| _jet_algorithm == cambridge_for_passive_algorithm) {
// for inclusive jets with a plugin algorithm, we make no
// assumptions about anything (relation of dij to momenta,
// ordering of the dij, etc.)
while (i >= 0) {
if (_history[i].parent2 == BeamJet) {
int parent1 = _history[i].parent1;
const PseudoJet & jet = _jets[_history[parent1].jetp_index];
if (jet.perp2() >= dcut) {jets.push_back(jet);}
}
i--;
}
} else {throw Error("cs::inclusive_jets(...): Unrecognized jet algorithm");}
return jets;
}
| const JetDefinition& ClusterSequence::jet_def | ( | ) | const [inline] |
return a reference to the jet definition
Definition at line 264 of file ClusterSequence.hh.
References _jet_def.
Referenced by ClusterSequenceArea::_warn_if_range_unsuitable(), and ClusterSequenceAreaBase::n_empty_jets().
{return _jet_def;}
| double ClusterSequence::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-algorithm, 1 for the Cambridge algorithm).
[May become virtual at some point]
Definition at line 322 of file ClusterSequence.cc.
References antikt_algorithm, cambridge_algorithm, cambridge_for_passive_algorithm, genkt_algorithm, and kt_algorithm.
Referenced by _bj_set_jetinfo().
{
if (_jet_algorithm == kt_algorithm) {return jet.kt2();}
else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
else if (_jet_algorithm == antikt_algorithm) {
double kt2=jet.kt2();
return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
} else if (_jet_algorithm == genkt_algorithm) {
double kt2 = jet.kt2();
double p = jet_def().extra_param();
if (p <= 0 && kt2 < 1e-300) kt2 = 1e-300; // dodgy safety check
return pow(kt2, p);
} else if (_jet_algorithm == cambridge_for_passive_algorithm) {
double kt2 = jet.kt2();
double lim = _jet_def.extra_param();
if (kt2 < lim*lim && kt2 != 0.0) {
return 1.0/kt2;
} else {return 1.0;}
} else {throw Error("Unrecognised jet algorithm");}
}
| const std::vector< PseudoJet > & ClusterSequence::jets | ( | ) | const [inline] |
allow the user to access the internally stored _jets() array, which contains both the initial particles and the various intermediate and final stages of recombination.
The first n_particles() entries are the original particles, in the order in which they were supplied to the ClusterSequence constructor. It can be useful to access them for example when examining whether a given input object is part of a specific jet, via the objects_in_jet(...) member function (which only takes PseudoJets that are registered in the ClusterSequence).
One of the other (internal uses) is related to the fact because we don't seem to be able to access protected elements of the class for an object that is not "this" (at least in case where "this" is of a slightly different kind from the object, both derived from ClusterSequence).
Definition at line 826 of file ClusterSequence.hh.
References _jets.
Referenced by ClusterSequenceActiveArea::_transfer_areas(), TrackJetPlugin::run_clustering(), SISConeSphericalPlugin::run_clustering(), SISConePlugin::run_clustering(), PxConePlugin::run_clustering(), NestedDefsPlugin::run_clustering(), JadePlugin::run_clustering(), EECambridgePlugin::run_clustering(), D0RunIIConePlugin::run_clustering(), CMSIterativeConePlugin::run_clustering(), CDFMidPointPlugin::run_clustering(), CDFJetCluPlugin::run_clustering(), and ATLASConePlugin::run_clustering().
{
return _jets;
}
| int ClusterSequence::n_exclusive_jets | ( | const double & | dcut | ) | const |
return the number of jets (in the sense of the exclusive algorithm) that would be obtained when running the algorithm with the given dcut.
Definition at line 439 of file ClusterSequence.cc.
Referenced by n_exclusive_jets_ycut().
{
// first locate the point where clustering would have stopped (i.e. the
// first time max_dij_so_far > dcut)
int i = _history.size() - 1; // last jet
while (i >= 0) {
if (_history[i].max_dij_so_far <= dcut) {break;}
i--;
}
int stop_point = i + 1;
// relation between stop_point, njets assumes one extra jet disappears
// at each clustering.
int njets = 2*_initial_n - stop_point;
return njets;
}
| int ClusterSequence::n_exclusive_jets_ycut | ( | double | ycut | ) | const [inline] |
the number of exclusive jets at the given ycut
Definition at line 142 of file ClusterSequence.hh.
References n_exclusive_jets(), and Q2().
Referenced by exclusive_jets_ycut().
{return n_exclusive_jets(ycut*Q2());}
| int ClusterSequence::n_exclusive_subjets | ( | const PseudoJet & | jet, |
| const double & | dcut | ||
| ) | const |
return the size of exclusive_subjets(...); still n ln n with same coefficient, but marginally more efficient than manually taking exclusive_subjets.size()
Definition at line 584 of file ClusterSequence.cc.
{
set<const history_element*> subhist;
// get the set of history elements that correspond to subjets at
// scale dcut
get_subhist_set(subhist, jet, dcut, 0);
return subhist.size();
}
| unsigned int ClusterSequence::n_particles | ( | ) | const [inline] |
returns the number of particles that were provided to the clustering algorithm (helps the user find their way around the history and jets objects if they weren't paying attention beforehand).
Definition at line 834 of file ClusterSequence.hh.
References _initial_n.
Referenced by ClusterSequenceActiveArea::_transfer_areas().
{return _initial_n;}
| bool ClusterSequence::object_in_jet | ( | const PseudoJet & | object, |
| const PseudoJet & | jet | ||
| ) | const |
returns true iff the object is included in the jet.
NB: this is only sensible if the object is already registered within the cluster sequence, so you cannot use it with an input particle to the CS (since the particle won't have the history index set properly).
For nice clustering structures it should run in O(ln(N)) time but in worst cases (certain cone plugins) it can take O(n) time, where n is the number of particles in the jet.
Definition at line 699 of file ClusterSequence.cc.
{
// make sure the object conceivably belongs to this clustering
// sequence
assert(_potentially_valid(object) && _potentially_valid(jet));
const PseudoJet * this_object = &object;
const PseudoJet * childp;
while(true) {
if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
return true;
} else if (has_child(*this_object, childp)) {this_object = childp;}
else {return false;}
}
}
| vector< int > ClusterSequence::particle_jet_indices | ( | const std::vector< PseudoJet > & | ) | const |
returns a vector of size n_particles() which indicates, for each of the initial particles (in the order in which they were supplied), which of the supplied jets it belongs to; if it does not belong to any of the supplied jets, the index is set to -1;
Definition at line 878 of file ClusterSequence.cc.
{
vector<int> indices(n_particles());
// first label all particles as not belonging to any jets
for (unsigned ipart = 0; ipart < n_particles(); ipart++)
indices[ipart] = -1;
// then for each of the jets relabel its consituents as belonging to
// that jet
for (unsigned ijet = 0; ijet < jets.size(); ijet++) {
vector<PseudoJet> jet_constituents(constituents(jets[ijet]));
for (unsigned ip = 0; ip < jet_constituents.size(); ip++) {
// a safe (if slightly redundant) way of getting the particle
// index (for initial particles it is actually safe to assume
// ipart=iclust).
unsigned iclust = jet_constituents[ip].cluster_hist_index();
unsigned ipart = history()[iclust].jetp_index;
indices[ipart] = ijet;
}
}
return indices;
}
| bool ClusterSequence::plugin_activated | ( | ) | const [inline] |
returns true when the plugin is allowed to run the show.
Definition at line 317 of file ClusterSequence.hh.
References _plugin_activated.
Referenced by plugin_record_iB_recombination(), plugin_record_ij_recombination(), and plugin_simple_N2_cluster().
{return _plugin_activated;}
| void ClusterSequence::plugin_associate_extras | ( | std::auto_ptr< Extras > | extras_in | ) | [inline] |
the plugin can associated some extra information with the ClusterSequence object by calling this function
Definition at line 312 of file ClusterSequence.hh.
References _extras.
Referenced by SISConeSphericalPlugin::run_clustering(), and SISConePlugin::run_clustering().
{
_extras = extras_in;
}
| void ClusterSequence::plugin_record_iB_recombination | ( | int | jet_i, |
| double | diB | ||
| ) | [inline] |
record the fact that there has been a recombination between jets()[jet_i] and the beam, with the specified diB; when looking for inclusive jets, any iB recombination will returned to the user as a jet.
Definition at line 296 of file ClusterSequence.hh.
References _do_iB_recombination_step(), and plugin_activated().
Referenced by TrackJetPlugin::run_clustering(), SISConeSphericalPlugin::run_clustering(), SISConePlugin::run_clustering(), PxConePlugin::run_clustering(), NestedDefsPlugin::run_clustering(), JadePlugin::run_clustering(), EECambridgePlugin::run_clustering(), D0RunIIConePlugin::run_clustering(), CMSIterativeConePlugin::run_clustering(), CDFMidPointPlugin::run_clustering(), CDFJetCluPlugin::run_clustering(), and ATLASConePlugin::run_clustering().
{
assert(plugin_activated());
_do_iB_recombination_step(jet_i, diB);
}
| void ClusterSequence::plugin_record_ij_recombination | ( | int | jet_i, |
| int | jet_j, | ||
| double | dij, | ||
| int & | newjet_k | ||
| ) | [inline] |
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
Definition at line 279 of file ClusterSequence.hh.
References _do_ij_recombination_step(), and plugin_activated().
Referenced by TrackJetPlugin::run_clustering(), SISConeSphericalPlugin::run_clustering(), SISConePlugin::run_clustering(), PxConePlugin::run_clustering(), NestedDefsPlugin::run_clustering(), JadePlugin::run_clustering(), EECambridgePlugin::run_clustering(), D0RunIIConePlugin::run_clustering(), CMSIterativeConePlugin::run_clustering(), CDFMidPointPlugin::run_clustering(), CDFJetCluPlugin::run_clustering(), and ATLASConePlugin::run_clustering().
{
assert(plugin_activated());
_do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
}
| void ClusterSequence::plugin_record_ij_recombination | ( | int | jet_i, |
| int | jet_j, | ||
| double | dij, | ||
| const PseudoJet & | newjet, | ||
| int & | newjet_k | ||
| ) |
as for the simpler variant of plugin_record_ij_recombination, except that the new jet is attributed the momentum and user_index of newjet
Definition at line 372 of file ClusterSequence.cc.
{
plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
// now transfer newjet into place
int tmp_index = _jets[newjet_k].cluster_hist_index();
_jets[newjet_k] = newjet;
_jets[newjet_k].set_cluster_hist_index(tmp_index);
}
| void ClusterSequence::plugin_simple_N2_cluster | ( | ) | [inline] |
allows a plugin to run a templated clustering (nearest-neighbour heuristic)
This has N^2 behaviour on "good" distance, but a worst case behaviour of N^3 (and many algs trigger the worst case behaviour)
For more details on how this works, see GenBriefJet below
Definition at line 329 of file ClusterSequence.hh.
References plugin_activated().
{
assert(plugin_activated());
_simple_N2_cluster<GBJ>();
}
| void ClusterSequence::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 format is: jet-n jet-px jet-py jet-pz jet-E particle-n particle-rap particle-phi particle-pt particle-n particle-rap particle-phi particle-pt ...
#END ... [i.e. above repeated]
| void ClusterSequence::print_jets_for_root | ( | const std::vector< PseudoJet > & | jets, |
| const std::string & | filename, | ||
| const std::string & | comment = "" |
||
| ) | const |
print jets for root to the file labelled filename, with an optional comment at the beginning
Definition at line 852 of file ClusterSequence.cc.
{
std::ofstream ostr(filename.c_str());
if (comment != "") ostr << "# " << comment << endl;
print_jets_for_root(jets, ostr);
}
| double ClusterSequence::Q | ( | ) | const [inline] |
returns the sum of all energies in the event (relevant mainly for e+e-)
Definition at line 192 of file ClusterSequence.hh.
References _Qtot.
{return _Qtot;}
| double ClusterSequence::Q2 | ( | ) | const [inline] |
return Q()^2
Definition at line 194 of file ClusterSequence.hh.
References _Qtot.
Referenced by exclusive_ymerge(), exclusive_ymerge_max(), n_exclusive_jets_ycut(), and EECambridgePlugin::run_clustering().
| static void ClusterSequence::set_jet_algorithm | ( | JetAlgorithm | jet_algorithm | ) | [inline, static] |
set the default (static) jet finder across all current and future ClusterSequence objects -- deprecated and obsolescent (i.e.
may be suppressed in a future release).
Definition at line 371 of file ClusterSequence.hh.
References _default_jet_algorithm.
{_default_jet_algorithm = jet_algorithm;}
| static void ClusterSequence::set_jet_finder | ( | JetAlgorithm | jet_algorithm | ) | [inline, static] |
same as above for backward compatibility
Definition at line 373 of file ClusterSequence.hh.
References _default_jet_algorithm.
{_default_jet_algorithm = jet_algorithm;}
| string ClusterSequence::strategy_string | ( | ) | const |
Definition at line 288 of file ClusterSequence.cc.
References N2MinHeapTiled, N2Plain, N2PoorTiled, N2Tiled, N3Dumb, NlnN, NlnN3pi, NlnN4pi, NlnNCam, NlnNCam2pi2R, NlnNCam4pi, and plugin_strategy.
{
string strategy;
switch(_strategy) {
case NlnN:
strategy = "NlnN"; break;
case NlnN3pi:
strategy = "NlnN3pi"; break;
case NlnN4pi:
strategy = "NlnN4pi"; break;
case N2Plain:
strategy = "N2Plain"; break;
case N2Tiled:
strategy = "N2Tiled"; break;
case N2MinHeapTiled:
strategy = "N2MinHeapTiled"; break;
case N2PoorTiled:
strategy = "N2PoorTiled"; break;
case N3Dumb:
strategy = "N3Dumb"; break;
case NlnNCam4pi:
strategy = "NlnNCam4pi"; break;
case NlnNCam2pi2R:
strategy = "NlnNCam2pi2R"; break;
case NlnNCam:
strategy = "NlnNCam"; break; // 2piMultD
case plugin_strategy:
strategy = "plugin strategy"; break;
default:
strategy = "Unrecognized";
}
return strategy;
}
| Strategy ClusterSequence::strategy_used | ( | ) | const [inline] |
return the enum value of the strategy used to cluster the event
Definition at line 260 of file ClusterSequence.hh.
References _strategy.
Referenced by ClusterSequenceActiveArea::_transfer_ghost_free_history().
{return _strategy;}
| void ClusterSequence::transfer_from_sequence | ( | ClusterSequence & | from_seq | ) |
transfer the sequence contained in other_seq into our own; any plugin "extras" contained in the from_seq will be lost from there.
Definition at line 348 of file ClusterSequence.cc.
References _extras, _history, _initial_n, _invR2, _jet_algorithm, _jet_def, _jets, _plugin_activated, _R2, _Rparam, _strategy, and _writeout_combinations.
Referenced by ClusterSequenceArea::initialize_and_run_cswa().
{
// the metadata
_jet_def = from_seq._jet_def ;
_writeout_combinations = from_seq._writeout_combinations ;
_initial_n = from_seq._initial_n ;
_Rparam = from_seq._Rparam ;
_R2 = from_seq._R2 ;
_invR2 = from_seq._invR2 ;
_strategy = from_seq._strategy ;
_jet_algorithm = from_seq._jet_algorithm ;
_plugin_activated = from_seq._plugin_activated ;
// the data
_jets = from_seq._jets;
_history = from_seq._history;
// the following transferse ownership of the extras from the from_seq
_extras = from_seq._extras;
}
| vector< PseudoJet > ClusterSequence::unclustered_particles | ( | ) | const |
return the set of particles that have not been clustered.
For kt and cam/aachen algorithms this should always be null, but for cone type algorithms it can be non-null;
Definition at line 1038 of file ClusterSequence.cc.
Referenced by ClusterSequenceActiveArea::_transfer_areas().
{
vector<PseudoJet> unclustered;
for (unsigned i = 0; i < n_particles() ; i++) {
if (_history[i].child == Invalid)
unclustered.push_back(_jets[_history[i].jetp_index]);
}
return unclustered;
}
| vector< int > ClusterSequence::unique_history_order | ( | ) | const |
routine that returns an order in which to read the history such that clusterings that lead to identical jet compositions but different histories (because of degeneracies in the clustering order) will have matching constituents for each matching entry in the unique_history_order.
The order has the property that an entry's parents will always appear prior to that entry itself.
Roughly speaking the order is such that we first provide all steps that lead to the final jet containing particle 1; then we have the steps that lead to reconstruction of the jet containing the next-lowest-numbered unclustered particle, etc... [see GPS CCN28-12 for more info -- of course a full explanation here would be better...]
Definition at line 983 of file ClusterSequence.cc.
References d0::inline_maths::min().
Referenced by ClusterSequenceActiveArea::_transfer_areas().
{
// first construct an array that will tell us the lowest constituent
// of a given jet -- this will always be one of the original
// particles, whose order is well defined and so will help us to
// follow the tree in a unique manner.
valarray<int> lowest_constituent(_history.size());
int hist_n = _history.size();
lowest_constituent = hist_n; // give it a large number
for (int i = 0; i < hist_n; i++) {
// sets things up for the initial partons
lowest_constituent[i] = min(lowest_constituent[i],i);
// propagates them through to the children of this parton
if (_history[i].child > 0) lowest_constituent[_history[i].child]
= min(lowest_constituent[_history[i].child],lowest_constituent[i]);
}
// establish an array for what we have and have not extracted so far
valarray<bool> extracted(_history.size()); extracted = false;
vector<int> unique_tree;
unique_tree.reserve(_history.size());
// now work our way through the tree
for (unsigned i = 0; i < n_particles(); i++) {
if (!extracted[i]) {
unique_tree.push_back(i);
extracted[i] = true;
_extract_tree_children(i, extracted, lowest_constituent, unique_tree);
}
}
return unique_tree;
}
Member Data Documentation
JetAlgorithm ClusterSequence::_default_jet_algorithm = kt_algorithm [static, protected] |
Definition at line 482 of file ClusterSequence.hh.
Referenced by set_jet_algorithm(), and set_jet_finder().
std::auto_ptr<Extras> ClusterSequence::_extras [private] |
Definition at line 564 of file ClusterSequence.hh.
Referenced by extras(), plugin_associate_extras(), and transfer_from_sequence().
bool ClusterSequence::_first_time = true [static, private] |
will be set by default to be true for the first run
Definition at line 612 of file ClusterSequence.hh.
std::vector<history_element> ClusterSequence::_history [protected] |
this vector will contain the branching history; for each stage, _history[i].jetp_index indicates where to look in the _jets vector to get the physical PseudoJet.
Definition at line 540 of file ClusterSequence.hh.
Referenced by _potentially_valid(), history(), and transfer_from_sequence().
int ClusterSequence::_initial_n [protected] |
Definition at line 554 of file ClusterSequence.hh.
Referenced by n_particles(), and transfer_from_sequence().
double ClusterSequence::_invR2 [protected] |
Definition at line 555 of file ClusterSequence.hh.
Referenced by transfer_from_sequence().
JetAlgorithm ClusterSequence::_jet_algorithm [protected] |
Definition at line 558 of file ClusterSequence.hh.
Referenced by transfer_from_sequence().
JetDefinition ClusterSequence::_jet_def [protected] |
Definition at line 483 of file ClusterSequence.hh.
Referenced by jet_def(), and transfer_from_sequence().
std::vector<PseudoJet> ClusterSequence::_jets [protected] |
This contains the physical PseudoJets; for each PseudoJet one can find the corresponding position in the _history by looking at _jets[i].cluster_hist_index().
Definition at line 534 of file ClusterSequence.hh.
Referenced by ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts(), _bj_of_hindex(), _bj_set_jetinfo(), _transfer_input_jets(), jets(), and transfer_from_sequence().
int ClusterSequence::_n_exclusive_warnings = 0 [static, private] |
record the number of warnings provided about the exclusive algorithm -- so that we don't print it out more than a few times.
Definition at line 617 of file ClusterSequence.hh.
int ClusterSequence::_n_tiles_phi [private] |
Definition at line 724 of file ClusterSequence.hh.
Referenced by _tile_index().
bool ClusterSequence::_plugin_activated [private] |
Definition at line 563 of file ClusterSequence.hh.
Referenced by plugin_activated(), and transfer_from_sequence().
double ClusterSequence::_Qtot [protected] |
Definition at line 556 of file ClusterSequence.hh.
double ClusterSequence::_R2 [protected] |
Definition at line 555 of file ClusterSequence.hh.
Referenced by _bj_set_jetinfo(), _bj_set_NN_crosscheck(), _bj_set_NN_nocross(), and transfer_from_sequence().
double ClusterSequence::_Rparam [protected] |
Definition at line 555 of file ClusterSequence.hh.
Referenced by transfer_from_sequence().
Strategy ClusterSequence::_strategy [protected] |
Definition at line 557 of file ClusterSequence.hh.
Referenced by strategy_used(), and transfer_from_sequence().
double ClusterSequence::_tile_size_eta [private] |
Definition at line 723 of file ClusterSequence.hh.
double ClusterSequence::_tile_size_phi [private] |
Definition at line 723 of file ClusterSequence.hh.
std::vector<Tile> ClusterSequence::_tiles [private] |
Definition at line 721 of file ClusterSequence.hh.
double ClusterSequence::_tiles_eta_max [private] |
Definition at line 722 of file ClusterSequence.hh.
double ClusterSequence::_tiles_eta_min [private] |
Definition at line 722 of file ClusterSequence.hh.
int ClusterSequence::_tiles_ieta_max [private] |
Definition at line 724 of file ClusterSequence.hh.
int ClusterSequence::_tiles_ieta_min [private] |
Definition at line 724 of file ClusterSequence.hh.
Referenced by _tile_index().
bool ClusterSequence::_writeout_combinations [protected] |
Definition at line 553 of file ClusterSequence.hh.
Referenced by transfer_from_sequence().
const int ClusterSequence::n_tile_neighbours = 9 [static, private] |
number of neighbours that a tile will have (rectangular geometry gives 9 neighbours).
Definition at line 703 of file ClusterSequence.hh.
The documentation for this class was generated from the following files:
- include/fastjet/ClusterSequence.hh
- src/ClusterSequence.cc
- src/ClusterSequence_CP2DChan.cc
- src/ClusterSequence_Delaunay.cc
- src/ClusterSequence_DumbN3.cc
- src/ClusterSequence_N2.cc
- src/ClusterSequence_TiledN2.cc
