fastjet::ClusterSequence Class Reference

deals with clustering More...

#include <ClusterSequence.hh>

Inheritance diagram for fastjet::ClusterSequence:

[legend]Collaboration diagram for fastjet::ClusterSequence:

[legend]List of all members.

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)
	constructor of a jet-clustering sequence from a vector of four-momenta, with the jet definition specified by jet_def
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
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.
std::vector< PseudoJet >	constituents (const PseudoJet &jet) const
	return a vector of the particles that make up jet
void	add_constituents (const PseudoJet &jet, std::vector< PseudoJet > &subjet_vector) const
	add on to subjet_vector the subjets of jet.
Strategy	strategy_used () const
	return the enum value of the strategy used to cluster the event
std::string	strategy_string () const
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)
const std::vector< PseudoJet > &	jets () const
	allow the user to access the jets in this raw manner (needed 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).
const std::vector< history_element > &	history () const
	allow the user to access the history in this raw manner (see above for motivation).
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.
Static Public Member Functions
static void	set_jet_finder (JetFinder jet_finder)
	set the default (static) jet finder across all current and future ClusterSequence objects -- deprecated and obsolescent (i.e.
Protected Member Functions
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)
	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.
void	_do_iB_recombination_step (const int &jet_i, const double &diB)
	carries out the bookkeeping associated with the step of recombining jet_i with the beam
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
Strategy	_strategy
JetFinder	_jet_finder
Static Protected Attributes
static JetFinder	_default_jet_finder = 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.
void	_simple_N2_cluster ()
	Order(N^2) clustering.
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)
	Add the current kt distance for particle ii to the map (DijMap) using information from the DNN object.
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: decide the range in eta allocate the tiles set up the cross-referencing info between 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.
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).
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...
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...

---

Detailed Description

deals with clustering

Definition at line 65 of file ClusterSequence.hh.

---

Member Typedef Documentation

typedef std::pair<double,TwoVertices> fastjet::ClusterSequence::DijEntry [private]

Definition at line 383 of file ClusterSequence.hh.

typedef std::multimap<double,TwoVertices> fastjet::ClusterSequence::DistMap [private]

Definition at line 384 of file ClusterSequence.hh.

typedef std::pair<int,int> fastjet::ClusterSequence::TwoVertices [private]

Definition at line 382 of file ClusterSequence.hh.

---

Member Enumeration Documentation

enum fastjet::ClusterSequence::JetType

Enumeration values: Invalid  InexistentParent  BeamJet  Definition at line 236 of file ClusterSequence.hh. {Invalid=-3, InexistentParent = -2, BeamJet = -1};

---

Constructor & Destructor Documentation

fastjet::ClusterSequence::ClusterSequence (   )  [inline]

default constructor Definition at line 71 of file ClusterSequence.hh. {}

template<class L> fastjet::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 558 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); }

template<class L> fastjet::ClusterSequence::ClusterSequence (  const std::vector< L > &  pseudojets, const JetDefinition &  jet_def, const bool &  writeout_combinations = false )

constructor of a jet-clustering sequence from a vector of four-momenta, with the jet definition specified by jet_def Definition at line 575 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); }

---

Member Function Documentation

void fastjet::ClusterSequence::_add_ktdistance_to_map (  const int &  ii, DistMap &  DijMap, const DynamicNearestNeighbours *  DNN )  [private]

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 _invR2, _jets, jet_scale_for_algorithm(), fastjet::DynamicNearestNeighbours::NearestNeighbourDistance(), and fastjet::DynamicNearestNeighbours::NearestNeighbourIndex(). Referenced by _delaunay_cluster(). { 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 fastjet::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 232 of file ClusterSequence_TiledN2.cc. References _tiles. Referenced by _tiled_N2_cluster(). { for (Tile * const * near_tile = _tiles[tile_index].begin_tiles; near_tile != _tiles[tile_index].end_tiles; near_tile++){ // get the tile number tile_union[n_near_tiles] = *near_tile - & _tiles[0]; n_near_tiles++; } }

void fastjet::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 504 of file ClusterSequence.cc. References _history, _jets, _writeout_combinations, fastjet::ClusterSequence::history_element::child, fastjet::ClusterSequence::history_element::dij, Invalid, fastjet::ClusterSequence::history_element::jetp_index, fastjet::ClusterSequence::history_element::max_dij_so_far, fastjet::ClusterSequence::history_element::parent1, and fastjet::ClusterSequence::history_element::parent2. Referenced by _do_iB_recombination_step(), and _do_ij_recombination_step(). { 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 fastjet::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 247 of file ClusterSequence_TiledN2.cc. References _tiles. Referenced by _faster_tiled_N2_cluster(), and _minheap_faster_tiled_N2_cluster(). { for (Tile ** near_tile = _tiles[tile_index].begin_tiles; near_tile != _tiles[tile_index].end_tiles; near_tile++){ if (! (*near_tile)->tagged) { (*near_tile)->tagged = true; // get the tile number tile_union[n_near_tiles] = *near_tile - & _tiles[0]; n_near_tiles++; } } }

template<class J> double fastjet::ClusterSequence::_bj_diJ (  const J *const   jeta  )  const [inline, private]

Definition at line 633 of file ClusterSequence.hh. Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), _simple_N2_cluster(), and _tiled_N2_cluster(). { double kt2 = jet->kt2; if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}} return jet->NN_dist * kt2; }

template<class J> double fastjet::ClusterSequence::_bj_dist (  const J *const   jeta, const J *const   jetb )  const [inline, private]

return the distance between two BriefJet objects Definition at line 624 of file ClusterSequence.hh. References fastjet::pi, and fastjet::twopi. Referenced by _bj_set_NN_crosscheck(), _bj_set_NN_nocross(), _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), _simple_N2_cluster(), and _tiled_N2_cluster(). { double dphi = std::abs(jetA->phi - jetB->phi); double deta = (jetA->eta - jetB->eta); if (dphi > pi) {dphi = twopi - dphi;} return dphi*dphi + deta*deta; }

template<class J> J* fastjet::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 450 of file ClusterSequence.hh. { J * res; for(res = head; res<tail; res++) { if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;} } return res; }

void fastjet::ClusterSequence::_bj_remove_from_tiles (  TiledJet *const   jet  )  [private]

Definition at line 49 of file ClusterSequence_TiledN2.cc. References _tiles, fastjet::ClusterSequence::Tile::head, fastjet::ClusterSequence::TiledJet::next, fastjet::ClusterSequence::TiledJet::previous, and fastjet::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 fastjet::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 Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster().

template<class J> void fastjet::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 609 of file ClusterSequence.hh. References _jets, _R2, and jet_scale_for_algorithm(). Referenced by _simple_N2_cluster(). { 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; }

template<class J> void fastjet::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 671 of file ClusterSequence.hh. References _bj_dist(), and _R2. Referenced by _simple_N2_cluster(). { double NN_dist = _R2; J * NN = NULL; for (J * jetB = head; jetB != tail; jetB++) { double dist = _bj_dist(jet,jetB); if (dist < NN_dist) { NN_dist = dist; NN = jetB; } if (dist < jetB->NN_dist) { jetB->NN_dist = dist; jetB->NN = jet; } } jet->NN = NN; jet->NN_dist = NN_dist; }

template<class J> void fastjet::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 643 of file ClusterSequence.hh. References _bj_dist(), and _R2. Referenced by _simple_N2_cluster(). { 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 fastjet::ClusterSequence::_CP2DChan_cluster (   )  [private]

a 4pi variant of the closest pair clustering Definition at line 223 of file ClusterSequence_CP2DChan.cc. References _do_Cambridge_inclusive_jets(), _do_ij_recombination_step(), _invR2, _jet_finder, _jets, BeamJet, fastjet::cambridge_algorithm, fastjet::ClosestPair2D::closest_pair(), Invalid, fastjet::ClosestPair2D::replace(), and fastjet::twopi. Referenced by _initialise_and_run(). { if (_jet_finder != 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 fastjet::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 _CP2DChan_limited_cluster(), _do_Cambridge_inclusive_jets(), _jet_finder, _Rparam, and fastjet::cambridge_algorithm. Referenced by _CP2DChan_cluster_2piMultD(), and _initialise_and_run(). { if (_jet_finder != 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 fastjet::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 _CP2DChan_cluster_2pi2R(), _CP2DChan_limited_cluster(), and _Rparam. Referenced by _initialise_and_run(). { // 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 fastjet::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 _do_ij_recombination_step(), _history, _initial_n, _invR2, _jets, fastjet::ClosestPair2D::closest_pair(), Invalid, fastjet::Private::make_mirror(), fastjet::pi, and fastjet::ClosestPair2D::replace_many(). Referenced by _CP2DChan_cluster_2pi2R(), and _CP2DChan_cluster_2piMultD(). { 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 fastjet::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 173 of file ClusterSequence.cc. References _invR2, _jet_def, _jet_finder, _plugin_activated, _print_banner(), _R2, _Rparam, _strategy, and _writeout_combinations. Referenced by _initialise_and_run(), and fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA(). { // 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_finder = jet_def.jet_finder(); _Rparam = jet_def.R(); _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2; _strategy = jet_def.strategy(); // disallow interference from the plugin _plugin_activated = false; }

void fastjet::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 _add_ktdistance_to_map(), _do_iB_recombination_step(), _do_ij_recombination_step(), _jets, _Rparam, _strategy, BeamJet, fastjet::NlnN, fastjet::NlnN3pi, fastjet::NlnN4pi, fastjet::DynamicNearestNeighbours::RemoveCombinedAddCombination(), fastjet::DynamicNearestNeighbours::RemovePoint(), fastjet::EtaPhi::sanitize(), strategy_string(), fastjet::twopi, and fastjet::DynamicNearestNeighbours::Valid(). Referenced by _initialise_and_run(). { 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 fastjet::ClusterSequence::_do_Cambridge_inclusive_jets (   )  [private]

Definition at line 336 of file ClusterSequence_CP2DChan.cc. References _do_iB_recombination_step(), _history, and Invalid. Referenced by _CP2DChan_cluster(), and _CP2DChan_cluster_2pi2R(). { 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 fastjet::ClusterSequence::_do_iB_recombination_step (  const int &  jet_i, const double &  diB )  [protected]

carries out the bookkeeping associated with the step of recombining jet_i with the beam Definition at line 681 of file ClusterSequence.cc. References _add_step_to_history(), _history, _jets, BeamJet, and Invalid. Referenced by _delaunay_cluster(), _do_Cambridge_inclusive_jets(), _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), _really_dumb_cluster(), _simple_N2_cluster(), _tiled_N2_cluster(), and fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(). { // 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 fastjet::ClusterSequence::_do_ij_recombination_step (  const int &  jet_i, const int &  jet_j, const double &  dij, int &  newjet_k )  [protected]

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 648 of file ClusterSequence.cc. References _add_step_to_history(), _history, _jet_def, _jets, and fastjet::JetDefinition::recombiner(). Referenced by _CP2DChan_cluster(), _CP2DChan_limited_cluster(), _delaunay_cluster(), _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), _really_dumb_cluster(), _simple_N2_cluster(), _tiled_N2_cluster(), and fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(). { // 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 fastjet::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) Definition at line 584 of file ClusterSequence.cc. References _extract_tree_parents(), and _history. Referenced by unique_history_order(). { if (!extracted[position]) { // that means we may have unidentified parents around, so go and // collect them (extracted[position]) will then be made true) _extract_tree_parents(position,extracted,lowest_constituent,unique_tree); } // now look after the children... int child = _history[position].child; if (child >= 0) _extract_tree_children(child,extracted,lowest_constituent,unique_tree); }

void fastjet::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) Definition at line 616 of file ClusterSequence.cc. References _history. Referenced by _extract_tree_children(). { if (!extracted[position]) { int parent1 = _history[position].parent1; int parent2 = _history[position].parent2; // where relevant order parents so that we will first treat the // one containing the smaller "lowest_constituent" if (parent1 >= 0 && parent2 >= 0) { if (lowest_constituent[parent1] > lowest_constituent[parent2]) swap(parent1, parent2); } // then actually run through the parents to extract the constituents... if (parent1 >= 0 && !extracted[parent1]) _extract_tree_parents(parent1,extracted,lowest_constituent,unique_tree); if (parent2 >= 0 && !extracted[parent2]) _extract_tree_parents(parent2,extracted,lowest_constituent,unique_tree); // finally declare this position to be accounted for and push it // onto our list. unique_tree.push_back(position); extracted[position] = true; } }

void fastjet::ClusterSequence::_faster_tiled_N2_cluster (   )  [private]

run a tiled clustering Definition at line 507 of file ClusterSequence_TiledN2.cc. References _add_untagged_neighbours_to_tile_union(), _bj_diJ(), _bj_dist(), _bj_remove_from_tiles(), _do_iB_recombination_step(), _do_ij_recombination_step(), _initialise_tiles(), _invR2, _jets, fastjet::ClusterSequence::TiledJet::_jets_index, _R2, _tiles, _tj_set_jetinfo(), fastjet::ClusterSequence::Tile::begin_tiles, fastjet::ClusterSequence::TiledJet::diJ_posn, fastjet::ClusterSequence::Tile::end_tiles, fastjet::ClusterSequence::Tile::head, n_tile_neighbours, fastjet::ClusterSequence::TiledJet::next, fastjet::ClusterSequence::TiledJet::NN, fastjet::ClusterSequence::TiledJet::NN_dist, fastjet::ClusterSequence::Tile::tagged, and fastjet::ClusterSequence::TiledJet::tile_index. Referenced by _initialise_and_run(). { _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 fastjet::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 195 of file ClusterSequence.cc. References _history, _initial_n, _jet_def, _jets, fastjet::ClusterSequence::history_element::child, fastjet::ClusterSequence::history_element::dij, InexistentParent, Invalid, fastjet::ClusterSequence::history_element::jetp_index, fastjet::ClusterSequence::history_element::max_dij_so_far, fastjet::ClusterSequence::history_element::parent1, fastjet::ClusterSequence::history_element::parent2, and fastjet::JetDefinition::recombiner(). Referenced by _initialise_and_run(), and fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA(). { 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); 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); } _initial_n = _jets.size(); }

void fastjet::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_finder. Definition at line 52 of file ClusterSequence.cc. References _default_jet_finder, and _initialise_and_run(). { JetDefinition jet_def(_default_jet_finder, R, strategy); _initialise_and_run(jet_def, writeout_combinations); }

void fastjet::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 63 of file ClusterSequence.cc. References _CP2DChan_cluster(), _CP2DChan_cluster_2pi2R(), _CP2DChan_cluster_2piMultD(), _decant_options(), _delaunay_cluster(), _faster_tiled_N2_cluster(), _fill_initial_history(), _jet_def, _jet_finder, _jets, _minheap_faster_tiled_N2_cluster(), _plugin_activated, _really_dumb_cluster(), _Rparam, _simple_N2_cluster(), _strategy, _tiled_N2_cluster(), fastjet::Best, fastjet::cambridge_algorithm, fastjet::JetDefinition::jet_finder(), fastjet::N2MinHeapTiled, fastjet::N2Plain, fastjet::N2PoorTiled, fastjet::N2Tiled, fastjet::N3Dumb, fastjet::NlnN, fastjet::NlnN3pi, fastjet::NlnN4pi, fastjet::NlnNCam, fastjet::NlnNCam2pi2R, fastjet::NlnNCam4pi, fastjet::JetDefinition::plugin(), and fastjet::plugin_algorithm. Referenced by _initialise_and_run(), fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA(), and 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(); // run the plugin if that's what's decreed if (_jet_finder == plugin_algorithm) { _plugin_activated = true; _jet_def.plugin()->run_clustering( (*this) ); _plugin_activated = false; 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_finder() == cambridge_algorithm) { _strategy = NlnNCam;} else #ifndef DROP_CGAL if (N > 16000/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) { this->_simple_N2_cluster(); } 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 fastjet::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 85 of file ClusterSequence_TiledN2.cc. References _jets, _n_tiles_phi, _Rparam, _tile_index(), _tile_size_eta, _tile_size_phi, _tiles, _tiles_eta_max, _tiles_eta_min, _tiles_ieta_max, _tiles_ieta_min, fastjet::ClusterSequence::Tile::begin_tiles, fastjet::ClusterSequence::Tile::end_tiles, fastjet::ClusterSequence::Tile::head, fastjet::ClusterSequence::Tile::RH_tiles, fastjet::ClusterSequence::Tile::surrounding_tiles, fastjet::ClusterSequence::Tile::tagged, and fastjet::twopi. Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster(). { // first decide tile sizes _tile_size_eta = _Rparam; _n_tiles_phi = int(floor(twopi/_Rparam)); _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 fastjet::ClusterSequence::_minheap_faster_tiled_N2_cluster (   )  [private]

run a tiled clustering, with our minheap for keeping track of the smallest dij Definition at line 724 of file ClusterSequence_TiledN2.cc. References _add_untagged_neighbours_to_tile_union(), _bj_diJ(), _bj_dist(), _bj_remove_from_tiles(), _do_iB_recombination_step(), _do_ij_recombination_step(), _initialise_tiles(), _invR2, _jets, fastjet::ClusterSequence::TiledJet::_jets_index, _R2, _tiles, _tj_set_jetinfo(), fastjet::ClusterSequence::Tile::begin_tiles, fastjet::ClusterSequence::Tile::end_tiles, fastjet::ClusterSequence::Tile::head, fastjet::ClusterSequence::TiledJet::label_minheap_update_done(), fastjet::MinHeap::minloc(), fastjet::MinHeap::minval(), n_tile_neighbours, fastjet::ClusterSequence::TiledJet::next, fastjet::ClusterSequence::TiledJet::NN, fastjet::ClusterSequence::TiledJet::NN_dist, fastjet::MinHeap::remove(), fastjet::ClusterSequence::Tile::tagged, fastjet::ClusterSequence::TiledJet::tile_index, and fastjet::MinHeap::update(). Referenced by _initialise_and_run(). { _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; }

void fastjet::ClusterSequence::_print_banner (   )  [private]

for making sure the user knows what it is they're running... Definition at line 145 of file ClusterSequence.cc. References _first_time, and fastjet_version. Referenced by _decant_options(). { 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 Matteo Cacciari and Gavin Salam \n"; cout << "# http://www.lpthe.jussieu.fr/~salam/fastjet \n"; cout << "# \n"; cout << "# Longitudinally invariant Kt, and inclusive Cambridge/Aachen clustering\n"; cout << "# using fast geometric algorithms, with optional external jet-finder \n"; cout << "# plugins. Please cite hep-ph/0512210 if you use this code. \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)"; #ifndef DROP_CGAL cout << " as well as CGAL: http://www.cgal.org/"; #endif // DROP_CGAL cout << ".\n"; cout << "#-----------------------------------------------------------------------\n"; }

void fastjet::ClusterSequence::_print_tiles (  TiledJet *  briefjets  )  const [private]

output the contents of the tiles Definition at line 209 of file ClusterSequence_TiledN2.cc. References _tiles. { 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 fastjet::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 _do_iB_recombination_step(), _do_ij_recombination_step(), _invR2, _jets, BeamJet, and jet_scale_for_algorithm(). Referenced by _initialise_and_run(). { // 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 fastjet::ClusterSequence::_simple_N2_cluster (   )  [private]

Order(N^2) clustering. Definition at line 48 of file ClusterSequence_N2.cc. References _bj_diJ(), _bj_dist(), _bj_set_jetinfo(), _bj_set_NN_crosscheck(), _bj_set_NN_nocross(), _do_iB_recombination_step(), _do_ij_recombination_step(), _invR2, _jets, fastjet::ClusterSequence::BriefJet::_jets_index, and fastjet::ClusterSequence::BriefJet::NN. Referenced by _initialise_and_run(). { int n = _jets.size(); BriefJet * briefjets = new BriefJet[n]; BriefJet * jetA = briefjets, * jetB; // initialise the basic jet info for (int i = 0; i< n; i++) { _bj_set_jetinfo(jetA, i); jetA++; // move on to next entry of briefjets } BriefJet * tail = jetA; // a semaphore for the end of briefjets BriefJet * head = briefjets; // a nicer way of naming start // now initialise the NN distances: jetA will run from 1..n-1; and // jetB from 0..jetA-1 for (jetA = head + 1; jetA != tail; jetA++) { // set NN info for jetA based on jets running from head..jetA-1, // checking in the process whether jetA itself is an undiscovered // NN of one of those jets. _bj_set_NN_crosscheck(jetA, head, 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 (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 //OBS int hist_a = _jets[jetA->_jets_index].cluster_hist_index(); //OBS int 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]); //OBS int nn = _jets.size() - 1; //OBS _jets[nn].set_cluster_hist_index(history_location); //OBS // update history //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_set_jetinfo(jetB, nn); } else { // jet-beam recombination _do_iB_recombination_step(jetA->_jets_index, diJ_min); //OBS // get the hist_index //OBS int hist_a = _jets[jetA->_jets_index].cluster_hist_index(); //OBS _add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); } // now update our nearest neighbour info and diJ table // first reduce size of table tail--; n--; // Copy last jet contents and diJ info into position of jetA *jetA = *tail; diJ[jetA - head] = diJ[tail-head]; // Initialise jetB's NN distance as well as updating it for // other particles. // NB: by having different loops for jetB == or != NULL we could // perhaps save a few percent (usually avoid one if inside loop), // but will not do it for now because on laptop fluctuations are // too large to reliably measure a few percent difference... for (BriefJet * jetI = head; jetI != tail; jetI++) { // see if jetI had jetA or jetB as a NN -- if so recalculate the NN if (jetI->NN == jetA || jetI->NN == jetB) { _bj_set_NN_nocross(jetI, head, tail); diJ[jetI-head] = _bj_diJ(jetI); // update diJ } // check whether new jetB is closer than jetI's current NN and // if need be 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 jetI's NN is the new tail then relabel it so that it becomes jetA if (jetI->NN == tail) {jetI->NN = jetA;} } if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);} } // final cleaning up; delete[] diJ; delete[] briefjets; }

int fastjet::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 165 of file ClusterSequence_TiledN2.cc. References _n_tiles_phi, _tile_size_eta, _tile_size_phi, _tiles_eta_max, _tiles_eta_min, _tiles_ieta_max, _tiles_ieta_min, and fastjet::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); }

int fastjet::ClusterSequence::_tile_index (  int  ieta, int  iphi )  const [inline, private]

Definition at line 508 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tj_set_jetinfo(). { // 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; }

void fastjet::ClusterSequence::_tiled_N2_cluster (   )  [private]

run a tiled clustering Definition at line 264 of file ClusterSequence_TiledN2.cc. References _add_neighbours_to_tile_union(), _bj_diJ(), _bj_dist(), _bj_remove_from_tiles(), _do_iB_recombination_step(), _do_ij_recombination_step(), _initialise_tiles(), _invR2, _jets, fastjet::ClusterSequence::TiledJet::_jets_index, _R2, _tiles, _tj_set_jetinfo(), fastjet::ClusterSequence::Tile::begin_tiles, fastjet::ClusterSequence::Tile::end_tiles, fastjet::ClusterSequence::Tile::head, n_tile_neighbours, fastjet::ClusterSequence::TiledJet::next, fastjet::ClusterSequence::TiledJet::NN, fastjet::ClusterSequence::TiledJet::NN_dist, fastjet::ClusterSequence::TiledJet::previous, and fastjet::ClusterSequence::TiledJet::tile_index. Referenced by _initialise_and_run(). { _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 fastjet::ClusterSequence::_tj_set_jetinfo (  TiledJet *const   jet, const int  _jets_index )  [inline, private]

Definition at line 188 of file ClusterSequence_TiledN2.cc. References _tile_index(), _tiles, and fastjet::ClusterSequence::Tile::head. Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster(). { // 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; }

template<class L> void fastjet::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 539 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 fastjet::ClusterSequence::add_constituents (  const PseudoJet &  jet, std::vector< PseudoJet > &  subjet_vector )  const

add on to subjet_vector the subjets of jet. Definition at line 477 of file ClusterSequence.cc. References _history, _jets, BeamJet, fastjet::PseudoJet::cluster_hist_index(), and InexistentParent. Referenced by constituents(). { // 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 subjet_vector.push_back(jet); 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 > fastjet::ClusterSequence::constituents (  const PseudoJet &  jet  )  const

return a vector of the particles that make up jet Definition at line 434 of file ClusterSequence.cc. References add_constituents(). Referenced by main(), particle_jet_indices(), and print_jets(). { vector<PseudoJet> subjets; add_constituents(jet, subjets); return subjets; }

double fastjet::ClusterSequence::exclusive_dmerge (  const int &  njets  )  const

return the dmin corresponding to the recombination that went from n+1 to n jets Definition at line 413 of file ClusterSequence.cc. References _history, and _initial_n. { assert(njets > 0); if (njets >= _initial_n) {return 0.0;} return _history[2*_initial_n-njets-1].dij; }

double fastjet::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 425 of file ClusterSequence.cc. References _history, and _initial_n. { assert(njets > 0); if (njets >= _initial_n) {return 0.0;} return _history[2*_initial_n-njets-1].max_dij_so_far; }

vector< PseudoJet > fastjet::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 353 of file ClusterSequence.cc. References _history, _initial_n, _jet_def, _jets, _n_exclusive_warnings, fastjet::JetDefinition::jet_finder(), jets(), and fastjet::kt_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 if (_jet_def.jet_finder() != kt_algorithm && _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; }

vector< PseudoJet > fastjet::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 345 of file ClusterSequence.cc. References n_exclusive_jets(). Referenced by main(). { int njets = n_exclusive_jets(dcut); return exclusive_jets(njets); }

const Extras* fastjet::ClusterSequence::extras (   )  const [inline]

returns a pointer to the extras object (may be null) Definition at line 195 of file ClusterSequence.hh. Referenced by main(). {return _extras.get();}

const std::vector< ClusterSequence::history_element > & fastjet::ClusterSequence::history (   )  const [inline]

allow the user to access the history in this raw manner (see above for motivation). Definition at line 592 of file ClusterSequence.hh. References _history. Referenced by fastjet::ClusterSequenceActiveArea::_transfer_areas(), fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(), main(), and particle_jet_indices(). { return _history; }

vector< PseudoJet > fastjet::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 279 of file ClusterSequence.cc. References _history, _jet_finder, _jets, BeamJet, fastjet::cambridge_algorithm, jets(), fastjet::kt_algorithm, fastjet::PseudoJet::perp2(), and fastjet::plugin_algorithm. Referenced by main(), fastjet::ClusterSequenceActiveArea::parabolic_pt_per_unit_area(), fastjet::ClusterSequenceActiveArea::pt_per_unit_area(), and run_jet_finder(). { double dcut = ptmin*ptmin; int i = _history.size() - 1; // last jet vector<PseudoJet> jets; if (_jet_finder == 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_finder == 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_finder == plugin_algorithm) { // for inclusive jets with a plugin algorith, 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("Unrecognized jet algorithm");} return jets; }

double fastjet::ClusterSequence::jet_scale_for_algorithm (  const PseudoJet &  jet  )  const [inline]

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 600 of file ClusterSequence.hh. References _jet_finder, fastjet::cambridge_algorithm, fastjet::PseudoJet::kt2(), and fastjet::kt_algorithm. Referenced by _add_ktdistance_to_map(), _bj_set_jetinfo(), and _really_dumb_cluster(). { if (_jet_finder == kt_algorithm) {return jet.kt2();} else if (_jet_finder == cambridge_algorithm) {return 1.0;} else {throw Error("Unrecognised jet finder");} }

const std::vector< PseudoJet > & fastjet::ClusterSequence::jets (   )  const [inline]

allow the user to access the jets in this raw manner (needed 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 588 of file ClusterSequence.hh. References _jets. Referenced by fastjet::ClusterSequenceActiveArea::_transfer_areas(), exclusive_jets(), inclusive_jets(), fastjet::SISConePlugin::run_clustering(), fastjet::PxConePlugin::run_clustering(), fastjet::CDFMidPointPlugin::run_clustering(), and fastjet::CDFJetCluPlugin::run_clustering(). { return _jets; }

int fastjet::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 326 of file ClusterSequence.cc. References _history, and _initial_n. Referenced by exclusive_jets(). { // 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; }

unsigned int fastjet::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 596 of file ClusterSequence.hh. References _initial_n. Referenced by fastjet::ClusterSequenceActiveArea::_transfer_areas(), particle_jet_indices(), unclustered_particles(), and unique_history_order(). {return _initial_n;}

vector< int > fastjet::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 446 of file ClusterSequence.cc. References constituents(), history(), and n_particles(). { 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 fastjet::ClusterSequence::plugin_activated (   )  const [inline]

returns true when the plugin is allowed to run the show. Definition at line 192 of file ClusterSequence.hh. {return _plugin_activated;}

void fastjet::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 187 of file ClusterSequence.hh. Referenced by fastjet::SISConePlugin::run_clustering(). { _extras = extras_in; }

void fastjet::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 171 of file ClusterSequence.hh. Referenced by fastjet::SISConePlugin::run_clustering(), fastjet::PxConePlugin::run_clustering(), fastjet::CDFMidPointPlugin::run_clustering(), and fastjet::CDFJetCluPlugin::run_clustering(). { assert(plugin_activated()); _do_iB_recombination_step(jet_i, diB); }

void fastjet::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 264 of file ClusterSequence.cc. References _jets, and plugin_record_ij_recombination(). { 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 fastjet::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 154 of file ClusterSequence.hh. Referenced by plugin_record_ij_recombination(), fastjet::SISConePlugin::run_clustering(), fastjet::PxConePlugin::run_clustering(), fastjet::CDFMidPointPlugin::run_clustering(), and fastjet::CDFJetCluPlugin::run_clustering(). { assert(plugin_activated()); _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k); }

static void fastjet::ClusterSequence::set_jet_finder (  JetFinder  jet_finder  )  [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 201 of file ClusterSequence.hh. {_default_jet_finder = jet_finder;}

string fastjet::ClusterSequence::strategy_string (   )  const

Definition at line 227 of file ClusterSequence.cc. References _strategy, fastjet::N2MinHeapTiled, fastjet::N2Plain, fastjet::N2PoorTiled, fastjet::N2Tiled, fastjet::N3Dumb, fastjet::NlnN, fastjet::NlnN3pi, fastjet::NlnN4pi, fastjet::NlnNCam, fastjet::NlnNCam2pi2R, fastjet::NlnNCam4pi, and fastjet::plugin_strategy. Referenced by _delaunay_cluster(), and main(). { 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 fastjet::ClusterSequence::strategy_used (   )  const [inline]

return the enum value of the strategy used to cluster the event Definition at line 137 of file ClusterSequence.hh. Referenced by fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(). {return _strategy;}

vector< PseudoJet > fastjet::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 603 of file ClusterSequence.cc. References _history, _jets, Invalid, and n_particles(). { 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 > fastjet::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 548 of file ClusterSequence.cc. References _extract_tree_children(), _history, and n_particles(). Referenced by fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), and main(). { // 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

JetFinder fastjet::ClusterSequence::_default_jet_finder = kt_algorithm [static, protected]

Definition at line 47 of file ClusterSequence.cc. Referenced by _initialise_and_run().

std::auto_ptr<Extras> fastjet::ClusterSequence::_extras [private]

Definition at line 347 of file ClusterSequence.hh.

bool fastjet::ClusterSequence::_first_time = true [static, private]

will be set by default to be true for the first run Definition at line 140 of file ClusterSequence.cc. Referenced by _print_banner().

std::vector<history_element> fastjet::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 336 of file ClusterSequence.hh. Referenced by _add_step_to_history(), _CP2DChan_limited_cluster(), _do_Cambridge_inclusive_jets(), _do_iB_recombination_step(), _do_ij_recombination_step(), _extract_tree_children(), _extract_tree_parents(), _fill_initial_history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(), add_constituents(), exclusive_dmerge(), exclusive_dmerge_max(), exclusive_jets(), history(), inclusive_jets(), n_exclusive_jets(), unclustered_particles(), and unique_history_order().

int fastjet::ClusterSequence::_initial_n [protected]

Definition at line 339 of file ClusterSequence.hh. Referenced by _CP2DChan_limited_cluster(), _fill_initial_history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), exclusive_dmerge(), exclusive_dmerge_max(), exclusive_jets(), n_exclusive_jets(), and n_particles().

double fastjet::ClusterSequence::_invR2 [protected]

Definition at line 340 of file ClusterSequence.hh. Referenced by _add_ktdistance_to_map(), _CP2DChan_cluster(), _CP2DChan_limited_cluster(), _decant_options(), _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), _really_dumb_cluster(), _simple_N2_cluster(), and _tiled_N2_cluster().

JetDefinition fastjet::ClusterSequence::_jet_def [protected]

Definition at line 286 of file ClusterSequence.hh. Referenced by _decant_options(), _do_ij_recombination_step(), _fill_initial_history(), _initialise_and_run(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), and exclusive_jets().

JetFinder fastjet::ClusterSequence::_jet_finder [protected]

Definition at line 342 of file ClusterSequence.hh. Referenced by _CP2DChan_cluster(), _CP2DChan_cluster_2pi2R(), _decant_options(), _initialise_and_run(), inclusive_jets(), and jet_scale_for_algorithm().

std::vector<PseudoJet> fastjet::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 330 of file ClusterSequence.hh. Referenced by fastjet::ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts(), _add_ktdistance_to_map(), _add_step_to_history(), _bj_set_jetinfo(), _CP2DChan_cluster(), _CP2DChan_limited_cluster(), _delaunay_cluster(), _do_iB_recombination_step(), _do_ij_recombination_step(), _faster_tiled_N2_cluster(), _fill_initial_history(), _initialise_and_run(), fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA(), _initialise_tiles(), _minheap_faster_tiled_N2_cluster(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), _really_dumb_cluster(), _simple_N2_cluster(), _tiled_N2_cluster(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), _transfer_input_jets(), add_constituents(), exclusive_jets(), inclusive_jets(), jets(), plugin_record_ij_recombination(), and unclustered_particles().

int fastjet::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 141 of file ClusterSequence.cc. Referenced by exclusive_jets().

int fastjet::ClusterSequence::_n_tiles_phi [private]

Definition at line 504 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

bool fastjet::ClusterSequence::_plugin_activated [private]

Definition at line 346 of file ClusterSequence.hh. Referenced by _decant_options(), and _initialise_and_run().

double fastjet::ClusterSequence::_R2 [protected]

Definition at line 340 of file ClusterSequence.hh. Referenced by _bj_set_jetinfo(), _bj_set_NN_crosscheck(), _bj_set_NN_nocross(), _decant_options(), _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster().

double fastjet::ClusterSequence::_Rparam [protected]

Definition at line 340 of file ClusterSequence.hh. Referenced by _CP2DChan_cluster_2pi2R(), _CP2DChan_cluster_2piMultD(), _decant_options(), _delaunay_cluster(), _initialise_and_run(), and _initialise_tiles().

Strategy fastjet::ClusterSequence::_strategy [protected]

Definition at line 341 of file ClusterSequence.hh. Referenced by _decant_options(), _delaunay_cluster(), _initialise_and_run(), fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(), and strategy_string().

double fastjet::ClusterSequence::_tile_size_eta [private]

Definition at line 503 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

double fastjet::ClusterSequence::_tile_size_phi [private]

Definition at line 503 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

std::vector<Tile> fastjet::ClusterSequence::_tiles [private]

Definition at line 501 of file ClusterSequence.hh. Referenced by _add_neighbours_to_tile_union(), _add_untagged_neighbours_to_tile_union(), _bj_remove_from_tiles(), _faster_tiled_N2_cluster(), _initialise_tiles(), _minheap_faster_tiled_N2_cluster(), _print_tiles(), _tiled_N2_cluster(), and _tj_set_jetinfo().

double fastjet::ClusterSequence::_tiles_eta_max [private]

Definition at line 502 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

double fastjet::ClusterSequence::_tiles_eta_min [private]

Definition at line 502 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

int fastjet::ClusterSequence::_tiles_ieta_max [private]

Definition at line 504 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

int fastjet::ClusterSequence::_tiles_ieta_min [private]

Definition at line 504 of file ClusterSequence.hh. Referenced by _initialise_tiles(), and _tile_index().

bool fastjet::ClusterSequence::_writeout_combinations [protected]

Definition at line 338 of file ClusterSequence.hh. Referenced by _add_step_to_history(), and _decant_options().

const int fastjet::ClusterSequence::n_tile_neighbours = 9 [static, private]

number of neighbours that a tile will have (rectangular geometry gives 9 neighbours). Definition at line 483 of file ClusterSequence.hh. Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster().

---

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

---