fastjet: fastjet::ClusterSequence Class Reference

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

Definition at line 447 of file ClusterSequence.hh.

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

Definition at line 448 of file ClusterSequence.hh.

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

Definition at line 449 of file ClusterSequence.hh.

Enumerator:

Invalid
InexistentParent
BeamJet

Definition at line 290 of file ClusterSequence.hh.

00290 {Invalid=-3, InexistentParent = -2, BeamJet = -1};

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

default constructor

Definition at line 71 of file ClusterSequence.hh.

00071 {}

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`
	)			`[inline]`

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 623 of file ClusterSequence.hh.

References _initialise_and_run(), and _transfer_input_jets().

00627                                                                       {
00628 
00629   // transfer the initial jets (type L) into our own array
00630   _transfer_input_jets(pseudojets);
00631 
00632   // run the clustering
00633   _initialise_and_run(R,strategy,writeout_combinations);
00634 }

template<class L>

fastjet::ClusterSequence::ClusterSequence	(	const std::vector< L > &	pseudojets,
		const JetDefinition &	jet_def,
		const bool &	writeout_combinations = `false`
	)			`[inline]`

constructor of a jet-clustering sequence from a vector of four-momenta, with the jet definition specified by jet_def

Definition at line 640 of file ClusterSequence.hh.

References _initialise_and_run(), _transfer_input_jets(), and jet_def().

00643                                                                       {
00644 
00645   // transfer the initial jets (type L) into our own array
00646   _transfer_input_jets(pseudojets);
00647 
00648   // run the clustering
00649   _initialise_and_run(jet_def,writeout_combinations);
00650 }

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 330 of file ClusterSequence.cc.

References _history, _jet_algorithm, _jets, fastjet::antikt_algorithm, BeamJet, fastjet::cambridge_algorithm, fastjet::cambridge_for_passive_algorithm, jets(), fastjet::kt_algorithm, fastjet::PseudoJet::perp2(), and fastjet::plugin_algorithm.

Referenced by fastjet::ClusterSequenceAreaBase::empty_area(), fastjet::ClusterSequenceAreaBase::get_median_rho_and_sigma(), main(), fastjet::ClusterSequenceAreaBase::parabolic_pt_per_unit_area(), fastjet::ClusterSequenceActiveArea::parabolic_pt_per_unit_area(), fastjet::ClusterSequenceActiveArea::pt_per_unit_area(), run_jet_finder(), and fastjet::ClusterSequenceAreaBase::subtracted_jets().

00330                                                                             {
00331   double dcut = ptmin*ptmin;
00332   int i = _history.size() - 1; // last jet
00333   vector<PseudoJet> jets;
00334   if (_jet_algorithm == kt_algorithm) {
00335     while (i >= 0) {
00336       // with our specific definition of dij and diB (i.e. R appears only in 
00337       // dij), then dij==diB is the same as the jet.perp2() and we can exploit
00338       // this in selecting the jets...
00339       if (_history[i].max_dij_so_far < dcut) {break;}
00340       if (_history[i].parent2 == BeamJet && _history[i].dij >= dcut) {
00341         // for beam jets
00342         int parent1 = _history[i].parent1;
00343         jets.push_back(_jets[_history[parent1].jetp_index]);}
00344       i--;
00345     }
00346   } else if (_jet_algorithm == cambridge_algorithm) {
00347     while (i >= 0) {
00348       // inclusive jets are all at end of clustering sequence in the
00349       // Cambridge algorithm -- so if we find a non-exclusive jet, then
00350       // we can exit
00351       if (_history[i].parent2 != BeamJet) {break;}
00352       int parent1 = _history[i].parent1;
00353       const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00354       if (jet.perp2() >= dcut) {jets.push_back(jet);}
00355       i--;
00356     }
00357   } else if (_jet_algorithm == plugin_algorithm 
00358              || _jet_algorithm == antikt_algorithm
00359              || _jet_algorithm == cambridge_for_passive_algorithm) {
00360     // for inclusive jets with a plugin algorithm, we make no
00361     // assumptions about anything (relation of dij to momenta,
00362     // ordering of the dij, etc.)
00363     while (i >= 0) {
00364       if (_history[i].parent2 == BeamJet) {
00365         int parent1 = _history[i].parent1;
00366         const PseudoJet & jet = _jets[_history[parent1].jetp_index];
00367         if (jet.perp2() >= dcut) {jets.push_back(jet);}
00368       }
00369       i--;
00370     }
00371   } else {throw Error("Unrecognized jet algorithm");}
00372   return jets;
00373 }

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 379 of file ClusterSequence.cc.

References _history, and _initial_n.

Referenced by exclusive_jets().

00379                                                                 {
00380 
00381   // first locate the point where clustering would have stopped (i.e. the
00382   // first time max_dij_so_far > dcut)
00383   int i = _history.size() - 1; // last jet
00384   while (i >= 0) {
00385     if (_history[i].max_dij_so_far <= dcut) {break;}
00386     i--;
00387   }
00388   int stop_point = i + 1;
00389   // relation between stop_point, njets assumes one extra jet disappears
00390   // at each clustering.
00391   int njets = 2*_initial_n - stop_point;
00392   return njets;
00393 }

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 398 of file ClusterSequence.cc.

References n_exclusive_jets().

Referenced by main().

00398                                                                             {
00399   int njets = n_exclusive_jets(dcut);
00400   return exclusive_jets(njets);
00401 }

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 406 of file ClusterSequence.cc.

References _history, _initial_n, _jet_def, _jets, _n_exclusive_warnings, fastjet::JetDefinition::jet_algorithm(), jets(), and fastjet::kt_algorithm.

00406                                                                           {
00407 
00408   // make sure the user does not ask for more than jets than there
00409   // were particles in the first place.
00410   assert (njets <= _initial_n);
00411 
00412   // provide a warning when extracting exclusive jets 
00413   if (_jet_def.jet_algorithm() != kt_algorithm && _n_exclusive_warnings < 5) {
00414     _n_exclusive_warnings++;
00415     cerr << "FastJet WARNING: dcut and exclusive jets for jet-finders other than kt should be interpreted with care." << endl;
00416   }
00417 
00418 
00419   // calculate the point where we have to stop the clustering.
00420   // relation between stop_point, njets assumes one extra jet disappears
00421   // at each clustering.
00422   int stop_point = 2*_initial_n - njets;
00423 
00424   // some sanity checking to make sure that e+e- does not give us
00425   // surprises (should we ever implement e+e-)...
00426   if (2*_initial_n != static_cast<int>(_history.size())) {
00427     ostringstream err;
00428     err << "2*_initial_n != _history.size() -- this endangers internal assumptions!\n";
00429     throw Error(err.str());
00430     //assert(false);
00431   }
00432 
00433   // now go forwards and reconstitute the jets that we have --
00434   // basically for any history element, see if the parent jets to
00435   // which it refers were created before the stopping point -- if they
00436   // were then add them to the list, otherwise they are subsequent
00437   // recombinations of the jets that we are looking for.
00438   vector<PseudoJet> jets;
00439   for (unsigned int i = stop_point; i < _history.size(); i++) {
00440     int parent1 = _history[i].parent1;
00441     if (parent1 < stop_point) {
00442       jets.push_back(_jets[_history[parent1].jetp_index]);
00443     }
00444     int parent2 = _history[i].parent2;
00445     if (parent2 < stop_point && parent2 > 0) {
00446       jets.push_back(_jets[_history[parent2].jetp_index]);
00447     }
00448     
00449   }
00450 
00451   // sanity check...
00452   if (static_cast<int>(jets.size()) != njets) {
00453     ostringstream err;
00454     err << "ClusterSequence::exclusive_jets: size of returned vector ("
00455          <<jets.size()<<") does not coincide with requested number of jets ("
00456          <<njets<<")";
00457     throw Error(err.str());
00458   }
00459 
00460   return jets;
00461 }

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 466 of file ClusterSequence.cc.

References _history, and _initial_n.

00466                                                                  {
00467   assert(njets > 0);
00468   if (njets >= _initial_n) {return 0.0;}
00469   return _history[2*_initial_n-njets-1].dij;
00470 }

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 478 of file ClusterSequence.cc.

References _history, and _initial_n.

00478                                                                      {
00479   assert(njets > 0);
00480   if (njets >= _initial_n) {return 0.0;}
00481   return _history[2*_initial_n-njets-1].max_dij_so_far;
00482 }

bool fastjet::ClusterSequence::object_in_jet	(	const PseudoJet &	object,
		const PseudoJet &	jet
	)			const

returns true iff the object is included in the jet.

NB: this is only sensible if the object is already registered within the cluster sequence, so you cannot use it with an input particle to the CS (since the particle won't have the history index set properly).

For nice clustering structures it should run in O(ln(N)) time but in worst cases (certain cone plugins) it can take O(n) time, where n is the number of particles in the jet.

Definition at line 487 of file ClusterSequence.cc.

References _potentially_valid(), fastjet::PseudoJet::cluster_hist_index(), and has_child().

00488                                                                  {
00489 
00490   // make sure the object conceivably belongs to this clustering
00491   // sequence
00492   assert(_potentially_valid(object) && _potentially_valid(jet));
00493 
00494   const PseudoJet * this_object = &object;
00495   const PseudoJet * childp;
00496   while(true) {
00497     if (this_object->cluster_hist_index() == jet.cluster_hist_index()) {
00498       return true;
00499     } else if (has_child(*this_object, childp)) {this_object = childp;}
00500     else {return false;}
00501   }
00502 }

bool fastjet::ClusterSequence::has_parents	(	const PseudoJet &	jet,
		PseudoJet &	parent1,
		PseudoJet &	parent2
	)			const

if the jet has parents in the clustering, it returns true and sets parent1 and parent2 equal to them.

if it has no parents it returns false and sets parent1 and parent2 to zero

Definition at line 510 of file ClusterSequence.cc.

References _history, _jets, fastjet::PseudoJet::cluster_hist_index(), fastjet::ClusterSequence::history_element::parent1, fastjet::ClusterSequence::history_element::parent2, and fastjet::PseudoJet::perp2().

00511                                                          {
00512 
00513   const history_element & hist = _history[jet.cluster_hist_index()];
00514 
00515   // make sure we do not run into any unexpected situations --
00516   // i.e. both parents valid, or neither
00517   assert ((hist.parent1 >= 0 && hist.parent2 >= 0) || 
00518           (hist.parent1 < 0 && hist.parent2 < 0));
00519 
00520   if (hist.parent1 < 0) {
00521     parent1 = PseudoJet(0.0,0.0,0.0,0.0);
00522     parent2 = parent1;
00523     return false;
00524   } else {
00525     parent1 = _jets[_history[hist.parent1].jetp_index];
00526     parent2 = _jets[_history[hist.parent2].jetp_index];
00527     // order the parents in decreasing pt
00528     if (parent1.perp2() < parent2.perp2()) swap(parent1,parent2);
00529     return true;
00530   }
00531 }

bool fastjet::ClusterSequence::has_child	(	const PseudoJet &	jet,
		PseudoJet &	child
	)			const

if the jet has a child then return true and give the child jet otherwise return false and set the child to zero

Definition at line 536 of file ClusterSequence.cc.

Referenced by object_in_jet().

00536                                                                               {
00537 
00538   //const history_element & hist = _history[jet.cluster_hist_index()];
00539   //
00540   //if (hist.child >= 0) {
00541   //  child = _jets[_history[hist.child].jetp_index];
00542   //  return true;
00543   //} else {
00544   //  child = PseudoJet(0.0,0.0,0.0,0.0);
00545   //  return false;
00546   //}
00547   const PseudoJet * childp;
00548   bool res = has_child(jet, childp);
00549   if (res) {
00550     child = *childp;
00551     return true;
00552   } else {
00553     child = PseudoJet(0.0,0.0,0.0,0.0);
00554     return false;
00555   }
00556 }

bool fastjet::ClusterSequence::has_child	(	const PseudoJet &	jet,
		const PseudoJet *&	childp
	)			const

Version of has_child that sets a pointer to the child if the child exists;.

Definition at line 558 of file ClusterSequence.cc.

References _history, _jets, fastjet::ClusterSequence::history_element::child, and fastjet::PseudoJet::cluster_hist_index().

00558                                                                                        {
00559 
00560   const history_element & hist = _history[jet.cluster_hist_index()];
00561 
00562   // check that this jet has a child and that the child corresponds to
00563   // a true jet [RETHINK-IF-CHANGE-NUMBERING: what is the right
00564   // behaviour if the child is the same jet but made inclusive...?]
00565   if (hist.child >= 0 && _history[hist.child].jetp_index >= 0) {
00566     childp = &(_jets[_history[hist.child].jetp_index]);
00567     return true;
00568   } else {
00569     childp = NULL;
00570     return false;
00571   }
00572 }

bool fastjet::ClusterSequence::has_partner	(	const PseudoJet &	jet,
		PseudoJet &	partner
	)			const

if this jet has a child (and so a partner) return true and give the partner, otherwise return false and set the partner to zero

Definition at line 579 of file ClusterSequence.cc.

References _history, _jets, fastjet::ClusterSequence::history_element::child, fastjet::PseudoJet::cluster_hist_index(), fastjet::ClusterSequence::history_element::parent1, and fastjet::ClusterSequence::history_element::parent2.

00580                                                          {
00581 
00582   const history_element & hist = _history[jet.cluster_hist_index()];
00583 
00584   // make sure we have a child and that the child does not correspond
00585   // to a clustering with the beam (or some other invalid quantity)
00586   if (hist.child >= 0 && _history[hist.child].parent2 >= 0) {
00587     const history_element & child_hist = _history[hist.child];
00588     if (child_hist.parent1 == jet.cluster_hist_index()) {
00589       // partner will be child's parent2 -- for iB clustering
00590       // parent2 will not be valid
00591       partner = _jets[_history[child_hist.parent2].jetp_index];
00592     } else {
00593       // partner will be child's parent1
00594       partner = _jets[_history[child_hist.parent1].jetp_index];
00595     }
00596     return true;
00597   } else {
00598     partner = PseudoJet(0.0,0.0,0.0,0.0);
00599     return false;
00600   }
00601 }

vector< PseudoJet > fastjet::ClusterSequence::constituents ( const PseudoJet & jet ) const

return a vector of the particles that make up jet

Definition at line 606 of file ClusterSequence.cc.

References add_constituents().

Referenced by print_jets().

00606                                                                             {
00607   vector<PseudoJet> subjets;
00608   add_constituents(jet, subjets);
00609   return subjets;
00610 }

void fastjet::ClusterSequence::print_jets_for_root	(	const std::vector< PseudoJet > &	jets,
		std::ostream &	ostr = `std::cout`
	)			const

output the supplied vector of jets in a format that can be read by an appropriate root script; the format is: jet-n jet-px jet-py jet-pz jet-E particle-n particle-rap particle-phi particle-pt particle-n particle-rap particle-phi particle-pt

.. END ... [i.e. above repeated]

void fastjet::ClusterSequence::add_constituents	(	const PseudoJet &	jet,
		std::vector< PseudoJet > &	subjet_vector
	)			const

add on to subjet_vector the subjets of jet (for internal use mainly)

Referenced by constituents().

Strategy fastjet::ClusterSequence::strategy_used ( ) const [inline]

return the enum value of the strategy used to cluster the event

Definition at line 187 of file ClusterSequence.hh.

Referenced by fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history().

00187 {return _strategy;}

string fastjet::ClusterSequence::strategy_string ( ) const

Definition at line 236 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().

00236                                                 {
00237   string strategy;
00238   switch(_strategy) {
00239   case NlnN:
00240     strategy = "NlnN"; break;
00241   case NlnN3pi:
00242     strategy = "NlnN3pi"; break;
00243   case NlnN4pi:
00244     strategy = "NlnN4pi"; break;
00245   case N2Plain:
00246     strategy = "N2Plain"; break;
00247   case N2Tiled:
00248     strategy = "N2Tiled"; break;
00249   case N2MinHeapTiled:
00250     strategy = "N2MinHeapTiled"; break;
00251   case N2PoorTiled:
00252     strategy = "N2PoorTiled"; break;
00253   case N3Dumb:
00254     strategy = "N3Dumb"; break;
00255   case NlnNCam4pi:
00256     strategy = "NlnNCam4pi"; break;
00257   case NlnNCam2pi2R:
00258     strategy = "NlnNCam2pi2R"; break;
00259   case NlnNCam:
00260     strategy = "NlnNCam"; break; // 2piMultD
00261   case plugin_strategy:
00262     strategy = "plugin strategy"; break;
00263   default:
00264     strategy = "Unrecognized";
00265   }
00266   return strategy;
00267 }

const JetDefinition& fastjet::ClusterSequence::jet_def ( ) const [inline]

return a reference to the jet definition

Definition at line 191 of file ClusterSequence.hh.

00191 {return _jet_def;}

double fastjet::ClusterSequence::jet_scale_for_algorithm ( const PseudoJet & jet ) const

returns the scale associated with a jet as required for this clustering algorithm (kt^2 for the kt-algorithm, 1 for the Cambridge algorithm).

[May become virtual at some point]

Definition at line 270 of file ClusterSequence.cc.

References _jet_algorithm, _jet_def, fastjet::antikt_algorithm, fastjet::cambridge_algorithm, fastjet::cambridge_for_passive_algorithm, fastjet::JetDefinition::extra_param(), fastjet::PseudoJet::kt2(), and fastjet::kt_algorithm.

Referenced by _add_ktdistance_to_map(), _bj_set_jetinfo(), and _really_dumb_cluster().

00271                                                                {
00272   if (_jet_algorithm == kt_algorithm)             {return jet.kt2();}
00273   else if (_jet_algorithm == cambridge_algorithm) {return 1.0;}
00274   else if (_jet_algorithm == antikt_algorithm) {
00275     double kt2=jet.kt2();
00276     return kt2 > 1e-300 ? 1.0/kt2 : 1e300;
00277   } else if (_jet_algorithm == cambridge_for_passive_algorithm) {
00278     double kt2 = jet.kt2();
00279     double lim = _jet_def.extra_param();
00280     if (kt2 < lim*lim && kt2 != 0.0) {
00281       return 1.0/kt2;
00282     } else {return 1.0;}
00283   } else {throw Error("Unrecognised jet finder");}
00284 }

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 206 of file ClusterSequence.hh.

Referenced by plugin_record_ij_recombination().

00207                                                       {
00208     assert(plugin_activated());
00209     _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
00210   }

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 315 of file ClusterSequence.cc.

References _jets, and plugin_record_ij_recombination().

00317                                                      {
00318 
00319   plugin_record_ij_recombination(jet_i, jet_j, dij, newjet_k);
00320 
00321   // now transfer newjet into place
00322   int tmp_index = _jets[newjet_k].cluster_hist_index();
00323   _jets[newjet_k] = newjet;
00324   _jets[newjet_k].set_cluster_hist_index(tmp_index);
00325 }

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 223 of file ClusterSequence.hh.

00223                                                              {
00224     assert(plugin_activated());
00225     _do_iB_recombination_step(jet_i, diB);
00226   }

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 239 of file ClusterSequence.hh.

00239                                                                      {
00240     _extras = extras_in;
00241   }

bool fastjet::ClusterSequence::plugin_activated ( ) const [inline]

returns true when the plugin is allowed to run the show.

Definition at line 244 of file ClusterSequence.hh.

00244 {return _plugin_activated;}

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

returns a pointer to the extras object (may be null)

Definition at line 247 of file ClusterSequence.hh.

00247 {return _extras.get();}

static void fastjet::ClusterSequence::set_jet_algorithm ( JetAlgorithm jet_algorithm ) [inline, static]

set the default (static) jet finder across all current and future ClusterSequence objects -- deprecated and obsolescent (i.e.

may be suppressed in a future release).

Definition at line 253 of file ClusterSequence.hh.

00253 {_default_jet_algorithm = jet_algorithm;}

static void fastjet::ClusterSequence::set_jet_finder ( JetAlgorithm jet_algorithm ) [inline, static]

same as above for backward compatibility

Definition at line 255 of file ClusterSequence.hh.

00255 {_default_jet_algorithm = jet_algorithm;}

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 653 of file ClusterSequence.hh.

References _jets.

Referenced by fastjet::ClusterSequenceActiveArea::_transfer_areas(), exclusive_jets(), inclusive_jets(), and fastjet::ClusterSequenceAreaBase::subtracted_jets().

00653                                                                  {
00654   return _jets;
00655 }

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 657 of file ClusterSequence.hh.

References _history.

Referenced by fastjet::ClusterSequenceActiveArea::_transfer_areas(), and fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history().

00657                                                                                          {
00658   return _history;
00659 }

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 661 of file ClusterSequence.hh.

References _initial_n.

Referenced by fastjet::ClusterSequenceVoronoiArea::_initializeVA(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), unclustered_particles(), and unique_history_order().

00661 {return _initial_n;}

std::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;

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 759 of file ClusterSequence.cc.

References _extract_tree_children(), _history, and n_particles().

Referenced by fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), fastjet::ClusterSequenceActiveArea::_run_AA(), and fastjet::ClusterSequenceActiveArea::_transfer_areas().

00759                                                         {
00760 
00761   // first construct an array that will tell us the lowest constituent
00762   // of a given jet -- this will always be one of the original
00763   // particles, whose order is well defined and so will help us to
00764   // follow the tree in a unique manner.
00765   valarray<int> lowest_constituent(_history.size());
00766   int hist_n = _history.size();
00767   lowest_constituent = hist_n; // give it a large number
00768   for (int i = 0; i < hist_n; i++) {
00769     // sets things up for the initial partons
00770     lowest_constituent[i] = min(lowest_constituent[i],i); 
00771     // propagates them through to the children of this parton
00772     if (_history[i].child > 0) lowest_constituent[_history[i].child] 
00773       = min(lowest_constituent[_history[i].child],lowest_constituent[i]);
00774   }
00775 
00776   // establish an array for what we have and have not extracted so far
00777   valarray<bool> extracted(_history.size()); extracted = false;
00778   vector<int> unique_tree;
00779   unique_tree.reserve(_history.size());
00780 
00781   // now work our way through the tree
00782   for (unsigned i = 0; i < n_particles(); i++) {
00783     if (!extracted[i]) {
00784       unique_tree.push_back(i);
00785       extracted[i] = true;
00786       _extract_tree_children(i, extracted, lowest_constituent, unique_tree);
00787     }
00788   }
00789 
00790   return unique_tree;
00791 }

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 814 of file ClusterSequence.cc.

References _history, _jets, Invalid, and n_particles().

Referenced by fastjet::ClusterSequenceActiveArea::_transfer_areas(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::empty_area(), main(), and run_jet_finder().

00814                                                                {
00815   vector<PseudoJet> unclustered;
00816   for (unsigned i = 0; i < n_particles() ; i++) {
00817     if (_history[i].child == Invalid) 
00818       unclustered.push_back(_jets[_history[i].jetp_index]);
00819   }
00820   return unclustered;
00821 }

void fastjet::ClusterSequence::transfer_from_sequence ( ClusterSequence & from_seq )

transfer the sequence contained in other_seq into our own; any plugin "extras" contained in the from_seq will be lost from there.

Definition at line 291 of file ClusterSequence.cc.

References _extras, _history, _initial_n, _invR2, _jet_algorithm, _jet_def, _jets, _plugin_activated, _R2, _Rparam, _strategy, and _writeout_combinations.

Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), and fastjet::ClusterSequenceArea::initialize_and_run_cswa().

00291                                                                        {
00292 
00293   // the metadata
00294   _jet_def                 = from_seq._jet_def                ;
00295   _writeout_combinations   = from_seq._writeout_combinations  ;
00296   _initial_n               = from_seq._initial_n              ;
00297   _Rparam                  = from_seq._Rparam                 ;
00298   _R2                      = from_seq._R2                     ;
00299   _invR2                   = from_seq._invR2                  ;
00300   _strategy                = from_seq._strategy               ;
00301   _jet_algorithm           = from_seq._jet_algorithm          ;
00302   _plugin_activated        = from_seq._plugin_activated       ;
00303 
00304   // the data
00305   _jets     = from_seq._jets;
00306   _history  = from_seq._history;
00307   // the following transferse ownership of the extras from the from_seq
00308   _extras   = from_seq._extras;
00309 
00310 }

bool fastjet::ClusterSequence::_potentially_valid ( const PseudoJet & jet ) const [inline, protected]

returns true if the jet has a history index contained within the range of this CS

Definition at line 348 of file ClusterSequence.hh.

References fastjet::PseudoJet::cluster_hist_index().

Referenced by object_in_jet().

00348                                                        {
00349     return jet.cluster_hist_index() >= 0 
00350       && jet.cluster_hist_index() < int(_history.size());
00351   }

template<class L>

void fastjet::ClusterSequence::_transfer_input_jets ( const std::vector< L > & pseudojets ) [inline, protected]

transfer the vector<L> of input jets into our own vector<PseudoJet> _jets (with some reserved space for future growth).

Definition at line 604 of file ClusterSequence.hh.

References _jets.

Referenced by ClusterSequence().

00605                                                                         {
00606 
00607   // this will ensure that we can point to jets without difficulties
00608   // arising.
00609   _jets.reserve(pseudojets.size()*2);
00610 
00611   // insert initial jets this way so that any type L that can be
00612   // converted to a pseudojet will work fine (basically PseudoJet
00613   // and any type that has [] subscript access to the momentum
00614   // components, such as CLHEP HepLorentzVector).
00615   for (unsigned int i = 0; i < pseudojets.size(); i++) {
00616     _jets.push_back(pseudojets[i]);}
00617   
00618 }

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_algorithm, _jet_def, _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_algorithm(), jet_def(), 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 fastjet::ClusterSequenceActiveArea::_initialise_AA(), _initialise_and_run(), and ClusterSequence().

00065                                                                       {
00066 
00067   // transfer all relevant info into internal variables
00068   _decant_options(jet_def, writeout_combinations);
00069 
00070   // set up the history entries for the initial particles (those
00071   // currently in _jets)
00072   _fill_initial_history();
00073 
00074   // run the plugin if that's what's decreed
00075   if (_jet_algorithm == plugin_algorithm) {
00076     _plugin_activated = true;
00077     _jet_def.plugin()->run_clustering( (*this) );
00078     _plugin_activated = false;
00079     return;
00080   }
00081 
00082 
00083   // automatically redefine the strategy according to N if that is
00084   // what the user requested -- transition points (and especially
00085   // their R-dependence) are based on empirical observations for a
00086   // R=0.4, 0.7 and 1.0, running on toth (3.4GHz, Pentium IV D [dual
00087   // core] with 2MB of cache).
00088   if (_strategy == Best) {
00089     int N = _jets.size();
00090     if (N > 6200/pow(_Rparam,2.0) 
00091         && jet_def.jet_algorithm() == cambridge_algorithm) {
00092       _strategy = NlnNCam;}
00093     else
00094 #ifndef DROP_CGAL
00095     if (N > 16000/pow(_Rparam,1.15)) {
00096       _strategy = NlnN; }   
00097     else                    
00098 #endif  // DROP_CGAL
00099       if (N > 450) {
00100       _strategy = N2MinHeapTiled;
00101     }
00102     else if (N > 55*max(0.5,min(1.0,_Rparam))) {// empirical scaling with R
00103       _strategy = N2Tiled;
00104     } else {
00105       _strategy = N2Plain;
00106     }
00107   }
00108 
00109 
00110   // run the code containing the selected strategy
00111   if (_strategy == NlnN || _strategy == NlnN3pi 
00112       || _strategy == NlnN4pi ) {
00113     this->_delaunay_cluster();
00114   } else if (_strategy ==  N3Dumb ) {
00115     this->_really_dumb_cluster();
00116   } else if (_strategy == N2Tiled) {
00117     this->_faster_tiled_N2_cluster();
00118   } else if (_strategy == N2PoorTiled) {
00119     this->_tiled_N2_cluster();
00120   } else if (_strategy == N2Plain) {
00121     this->_simple_N2_cluster();
00122   } else if (_strategy == N2MinHeapTiled) {
00123     this->_minheap_faster_tiled_N2_cluster();
00124   } else if (_strategy == NlnNCam4pi) {
00125     this->_CP2DChan_cluster();
00126   } else if (_strategy == NlnNCam2pi2R) {
00127     this->_CP2DChan_cluster_2pi2R();
00128   } else if (_strategy == NlnNCam) {
00129     this->_CP2DChan_cluster_2piMultD();
00130   } else {
00131     ostringstream err;
00132     err << "Unrecognised value for strategy: "<<_strategy;
00133     throw Error(err.str());
00134     //assert(false);
00135   }
00136 }

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_algorithm.

Definition at line 52 of file ClusterSequence.cc.

References _default_jet_algorithm, _initialise_and_run(), and jet_def().

00055                                                                       {
00056 
00057   JetDefinition jet_def(_default_jet_algorithm, R, strategy);
00058   _initialise_and_run(jet_def, writeout_combinations);
00059 }

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 182 of file ClusterSequence.cc.

References _invR2, _jet_algorithm, _jet_def, _plugin_activated, _print_banner(), _R2, _Rparam, _strategy, _writeout_combinations, fastjet::JetDefinition::jet_algorithm(), jet_def(), fastjet::JetDefinition::R(), and fastjet::JetDefinition::strategy().

Referenced by fastjet::ClusterSequenceActiveArea::_initialise_AA(), and _initialise_and_run().

00183                                                                           {
00184 
00185   // let the user know what's going on
00186   _print_banner();
00187 
00188   // make a local copy of the jet definition (for future use?)
00189   _jet_def = jet_def;
00190   
00191   _writeout_combinations = writeout_combinations;
00192   _jet_algorithm = jet_def.jet_algorithm();
00193   _Rparam = jet_def.R();  _R2 = _Rparam*_Rparam; _invR2 = 1.0/_R2;
00194   _strategy = jet_def.strategy();
00195 
00196   // disallow interference from the plugin
00197   _plugin_activated = false;
00198   
00199 }

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 204 of file ClusterSequence.cc.

References _history, _initial_n, _jet_def, _jets, InexistentParent, Invalid, and fastjet::JetDefinition::recombiner().

Referenced by fastjet::ClusterSequenceActiveArea::_initialise_AA(), and _initialise_and_run().

00204                                              {
00205 
00206   if (_jets.size() == 0) {throw Error("Cannot run jet-finder on empty event");}
00207 
00208   // reserve sufficient space for everything
00209   _jets.reserve(_jets.size()*2);
00210   _history.reserve(_jets.size()*2);
00211 
00212   for (int i = 0; i < static_cast<int>(_jets.size()) ; i++) {
00213     history_element element;
00214     element.parent1 = InexistentParent;
00215     element.parent2 = InexistentParent;
00216     element.child   = Invalid;
00217     element.jetp_index = i;
00218     element.dij     = 0.0;
00219     element.max_dij_so_far = 0.0;
00220 
00221     _history.push_back(element);
00222     
00223     // do any momentum preprocessing needed by the recombination scheme
00224     _jet_def.recombiner()->preprocess(_jets[i]);
00225 
00226     // get cross-referencing right from PseudoJets
00227     _jets[i].set_cluster_hist_index(i);
00228   }
00229   _initial_n = _jets.size();
00230 }

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 859 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().

00862                                                {
00863 
00864   // create the new jet by recombining the first two
00865   PseudoJet newjet;
00866   _jet_def.recombiner()->recombine(_jets[jet_i], _jets[jet_j], newjet);
00867   _jets.push_back(newjet);
00868   // original version...
00869   //_jets.push_back(_jets[jet_i] + _jets[jet_j]);
00870 
00871   // get its index
00872   newjet_k = _jets.size()-1;
00873 
00874   // get history index
00875   int newstep_k = _history.size();
00876   // and provide jet with the info
00877   _jets[newjet_k].set_cluster_hist_index(newstep_k);
00878 
00879   // finally sort out the history 
00880   int hist_i = _jets[jet_i].cluster_hist_index();
00881   int hist_j = _jets[jet_j].cluster_hist_index();
00882 
00883   _add_step_to_history(newstep_k, min(hist_i, hist_j), max(hist_i,hist_j),
00884                        newjet_k, dij);
00885 
00886 }

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 892 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().

00893                                                                          {
00894   // get history index
00895   int newstep_k = _history.size();
00896 
00897   // recombine the jet with the beam
00898   _add_step_to_history(newstep_k,_jets[jet_i].cluster_hist_index(),BeamJet,
00899                        Invalid, diB);
00900 
00901 }

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().

00050                                             {
00051 
00052   // the array that will be overwritten here will be one
00053   // of pointers to jets.
00054   vector<PseudoJet *> jetsp(_jets.size());
00055   vector<int>         indices(_jets.size());
00056 
00057   for (size_t i = 0; i<_jets.size(); i++) {
00058     jetsp[i] = & _jets[i];
00059     indices[i] = i;
00060   }
00061 
00062   for (int n = jetsp.size(); n > 0; n--) {
00063     int ii, jj;
00064     // find smallest beam distance [remember jet_scale_for_algorithm 
00065     // will return kt^2 for the kt-algorithm and 1 for the Cambridge/Aachen]
00066     double ymin = jet_scale_for_algorithm(*(jetsp[0]));
00067     ii = 0; jj = -2;
00068     for (int i = 0; i < n; i++) {
00069       double yiB = jet_scale_for_algorithm(*(jetsp[i]));
00070       if (yiB < ymin) {
00071         ymin = yiB; ii = i; jj = -2;}
00072     }
00073 
00074     // find smallest distance between pair of jetsp
00075     for (int i = 0; i < n-1; i++) {
00076       for (int j = i+1; j < n; j++) {
00077         //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2;
00078         double y = min(jet_scale_for_algorithm(*(jetsp[i])), 
00079                        jet_scale_for_algorithm(*(jetsp[j])))
00080                     * jetsp[i]->plain_distance(*jetsp[j])*_invR2;
00081         if (y < ymin) {ymin = y; ii = i; jj = j;}
00082       }
00083     }
00084 
00085     // output recombination sequence
00086     // old "ktclus" way of labelling
00087     //cout <<n<< " "<< ii+1 << " with " << jj+1 << "; y = "<< ymin<<endl;
00088     // new delaunay way of labelling
00089     int jjindex_or_beam, iiindex;
00090     if (jj < 0) {jjindex_or_beam = BeamJet; iiindex = indices[ii];} 
00091     else {
00092       jjindex_or_beam = max(indices[ii],indices[jj]);
00093       iiindex =         min(indices[ii],indices[jj]);
00094     }
00095 
00096     // now recombine
00097     int newn = 2*jetsp.size() - n;
00098     if (jj >= 0) {
00099       // combine pair
00100       int nn; // new jet index
00101       _do_ij_recombination_step(jetsp[ii]-&_jets[0], 
00102                                 jetsp[jj]-&_jets[0], ymin, nn);
00103       
00104       // sort out internal bookkeeping
00105       jetsp[ii] = &_jets[nn];
00106       // have jj point to jet that was pointed at by n-1 
00107       // (since original jj is no longer current, so put n-1 into jj)
00108       jetsp[jj] = jetsp[n-1];
00109       indices[ii] = newn;
00110       indices[jj] = indices[n-1];
00111 
00112       //OBS_jets.push_back(*jetsp[ii] + *jetsp[jj]);
00113       //OBSjetsp[ii] = &_jets[_jets.size()-1];
00114       //OBS// have jj point to jet that was pointed at by n-1 
00115       //OBS// (since original jj is no longer current, so put n-1 into jj)
00116       //OBSjetsp[jj] = jetsp[n-1];
00117       //OBS
00118       //OBSindices[ii] = newn;
00119       //OBSindices[jj] = indices[n-1];
00120       //OBS_add_step_to_history(newn,iiindex,
00121       //OBS                           jjindex_or_beam,_jets.size()-1,ymin);
00122     } else {
00123       // combine ii with beam
00124       _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin);
00125       // put last jet (pointer) in place of ii (which has disappeared)
00126       jetsp[ii] = jetsp[n-1];
00127       indices[ii] = indices[n-1];
00128       //OBS_add_step_to_history(newn,iiindex,jjindex_or_beam,Invalid, ymin);
00129     }
00130   }
00131 
00132 }

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().

00058                                          {
00059 
00060   int n = _jets.size();
00061 
00062   vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double>
00063   for (int i = 0; i < n; i++) {
00064     points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
00065     points[i].sanitize(); // make sure things are in the right range
00066   }
00067 
00068   // initialise our DNN structure with the set of points
00069   DynamicNearestNeighbours * DNN;
00070 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
00071   bool verbose = false;
00072   bool ignore_nearest_is_mirror = (_Rparam < twopi);
00073   if (_strategy == NlnN4pi) {
00074     DNN = new Dnn4piCylinder(points,verbose);
00075   } else if (_strategy == NlnN3pi) {
00076     DNN = new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose);
00077   } else if (_strategy == NlnN) {
00078     DNN = new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose);
00079   } else 
00080 #else
00081   if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
00082     ostringstream err;
00083     err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
00084     err << "       supported because FastJet was compiled without CGAL"<<endl;
00085     throw Error(err.str());
00086     //assert(false);
00087   }
00088 #endif // DROP_CGAL
00089   {
00090     ostringstream err;
00091     err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
00092     assert(false);
00093     throw Error(err.str());
00094   }
00095 
00096   // We will find nearest neighbour for each vertex, and include
00097   // distance in map (NB DistMap is a typedef given in the .h file)
00098   DistMap DijMap;
00099 
00100   // fill the map with the minimal (as far as we know) subset of Dij
00101   // distances (i.e. nearest neighbour ones).
00102   for (int ii = 0; ii < n; ii++) {
00103     _add_ktdistance_to_map(ii, DijMap, DNN);
00104   }
00105 
00106   // run the clustering (go up to i=n-1, but then will stop half-way down,
00107   // when we reach that point -- it will be the final beam jet and there
00108   // will be no nearest neighbours to find).
00109   for (int i=0;i<n;i++) {
00110     // find nearest vertices
00111     // NB: skip cases where the point is not there anymore!
00112     TwoVertices SmallestDijPair;
00113     int jet_i, jet_j;
00114     double SmallestDij;
00115     bool Valid2;
00116     bool recombine_with_beam;
00117     do { 
00118       SmallestDij = DijMap.begin()->first;
00119       SmallestDijPair = DijMap.begin()->second;
00120       jet_i = SmallestDijPair.first;
00121       jet_j = SmallestDijPair.second;
00122       // distance is immediately removed regardless of whether or not
00123       // it is used.
00124       // Some temporary testing code relating to problems with the gcc-3.2 compiler
00125       //cout << "got here and size is "<< DijMap.size()<< " and it is "<<SmallestDij <<"\n";
00126       //cout <<  jet_i << " "<< jet_j<<"\n";
00127       DijMap.erase(DijMap.begin());
00128       //cout << "got beyond here\n";
00129 
00130       // need to "prime" the validity of jet_j in such a way that 
00131       // if it corresponds to the beam then it is automatically valid.
00132       recombine_with_beam = (jet_j == BeamJet);
00133       if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);} 
00134       else {Valid2 = true;}
00135 
00136     } while ( !DNN->Valid(jet_i) || !Valid2);
00137 
00138 
00139     // The following part acts just on jet momenta and on the history.
00140     // The action on the nearest-neighbour structures takes place
00141     // later (only if at least 2 jets are around).
00142     if (! recombine_with_beam) {
00143       int nn; // will be index of new jet
00144       _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
00145       //OBS // merge the two jets, add new jet, remove old ones
00146       //OBS _jets.push_back(_jets[jet_i] + _jets[jet_j]);
00147       //OBS 
00148       //OBS int nn = _jets.size()-1;
00149       //OBS _jets[nn].set_cluster_hist_index(n+i);
00150       //OBS 
00151       //OBS // get corresponding indices in history structure
00152       //OBS int hist_i = _jets[jet_i].cluster_hist_index();
00153       //OBS int hist_j = _jets[jet_j].cluster_hist_index();
00154       //OBS 
00155       //OBS 
00156       //OBS _add_step_to_history(n+i,min(hist_i,hist_j), max(hist_i,hist_j),
00157       //OBS                   _jets.size()-1, SmallestDij);
00158 
00159       // add new point to points vector
00160       EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
00161       newpoint.sanitize(); // make sure it is in correct range
00162       points.push_back(newpoint);
00163     } else {
00164       // recombine the jet with the beam
00165       _do_iB_recombination_step(jet_i, SmallestDij);
00166       //OBS _add_step_to_history(n+i,_jets[jet_i].cluster_hist_index(),BeamJet,
00167       //OBS                        Invalid, SmallestDij);
00168     }
00169 
00170     // exit the loop because we do not want to look for nearest neighbours
00171     // etc. of zero partons
00172     if (i == n-1) {break;}
00173 
00174     vector<int> updated_neighbours;
00175     if (! recombine_with_beam) {
00176       // update DNN
00177       int point3;
00178       DNN->RemoveCombinedAddCombination(jet_i, jet_j, 
00179                                        points[points.size()-1], point3,
00180                                        updated_neighbours);
00181       // C++ beginners' comment: static_cast to unsigned int is necessary
00182       // to do away with warnings about type mismatch between point3 (int) 
00183       // and points.size (unsigned int)
00184       if (static_cast<unsigned int> (point3) != points.size()-1) {
00185         throw Error("INTERNAL ERROR: point3 != points.size()-1");}
00186     } else {
00187       // update DNN
00188       DNN->RemovePoint(jet_i, updated_neighbours);
00189     }
00190 
00191     // update map
00192     vector<int>::iterator it = updated_neighbours.begin();
00193     for (; it != updated_neighbours.end(); ++it) {
00194       int ii = *it;
00195       _add_ktdistance_to_map(ii, DijMap, DNN);
00196     }
00197       
00198   } // end clustering loop 
00199   
00200   // remember to clean up!
00201   delete DNN;
00202 }

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().

00048                                          {
00049   int n = _jets.size();
00050   BriefJet * briefjets = new BriefJet[n];
00051   BriefJet * jetA = briefjets, * jetB;
00052   
00053   // initialise the basic jet info 
00054   for (int i = 0; i< n; i++) {
00055     _bj_set_jetinfo(jetA, i);
00056     jetA++; // move on to next entry of briefjets
00057   }
00058   BriefJet * tail = jetA; // a semaphore for the end of briefjets
00059   BriefJet * head = briefjets; // a nicer way of naming start
00060 
00061   // now initialise the NN distances: jetA will run from 1..n-1; and
00062   // jetB from 0..jetA-1
00063   for (jetA = head + 1; jetA != tail; jetA++) {
00064     // set NN info for jetA based on jets running from head..jetA-1,
00065     // checking in the process whether jetA itself is an undiscovered
00066     // NN of one of those jets.
00067     _bj_set_NN_crosscheck(jetA, head, jetA);
00068   }
00069 
00070 
00071   // now create the diJ (where J is i's NN) table -- remember that 
00072   // we differ from standard normalisation here by a factor of R2
00073   double * diJ = new double[n];
00074   jetA = head;
00075   for (int i = 0; i < n; i++) {
00076     diJ[i] = _bj_diJ(jetA);
00077     jetA++; // have jetA follow i
00078   }
00079 
00080   // now run the recombination loop
00081   int history_location = n-1;
00082   while (tail != head) {
00083 
00084     // find the minimum of the diJ on this round
00085     double diJ_min = diJ[0];
00086     int diJ_min_jet = 0;
00087     for (int i = 1; i < n; i++) {
00088       if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min  = diJ[i];}
00089     }
00090     
00091     // do the recombination between A and B
00092     history_location++;
00093     jetA = & briefjets[diJ_min_jet];
00094     jetB = jetA->NN;
00095     // put the normalisation back in
00096     diJ_min *= _invR2; 
00097     if (jetB != NULL) {
00098       // jet-jet recombination
00099       // If necessary relabel A & B to ensure jetB < jetA, that way if
00100       // the larger of them == newtail then that ends up being jetA and 
00101       // the new jet that is added as jetB is inserted in a position that
00102       // has a future!
00103       if (jetA < jetB) {swap(jetA,jetB);}
00104 
00105       int nn; // new jet index
00106       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00107 
00108       //OBS // get the two history indices
00109       //OBS int hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00110       //OBS int hist_b = _jets[jetB->_jets_index].cluster_hist_index();
00111       //OBS // create the recombined jet
00112       //OBS _jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
00113       //OBS int nn = _jets.size() - 1;
00114       //OBS _jets[nn].set_cluster_hist_index(history_location);
00115       //OBS // update history
00116       //OBS _add_step_to_history(history_location, 
00117       //OBS                min(hist_a,hist_b),max(hist_a,hist_b),
00118       //OBS                        nn, diJ_min);
00119       // what was jetB will now become the new jet
00120       _bj_set_jetinfo(jetB, nn);
00121     } else {
00122       // jet-beam recombination
00123       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00124       //OBS // get the hist_index
00125       //OBS int hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00126       //OBS _add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); 
00127     }
00128 
00129     // now update our nearest neighbour info and diJ table
00130     // first reduce size of table
00131     tail--; n--;
00132     // Copy last jet contents and diJ info into position of jetA
00133     *jetA = *tail;
00134     diJ[jetA - head] = diJ[tail-head];
00135 
00136     // Initialise jetB's NN distance as well as updating it for 
00137     // other particles.
00138     // NB: by having different loops for jetB == or != NULL we could
00139     //     perhaps save a few percent (usually avoid one if inside loop),
00140     //     but will not do it for now because on laptop fluctuations are
00141     //     too large to reliably measure a few percent difference...
00142     for (BriefJet * jetI = head; jetI != tail; jetI++) {
00143       // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00144       if (jetI->NN == jetA || jetI->NN == jetB) {
00145         _bj_set_NN_nocross(jetI, head, tail);
00146         diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
00147       } 
00148       // check whether new jetB is closer than jetI's current NN and
00149       // if need be update things
00150       if (jetB != NULL) {
00151         double dist = _bj_dist(jetI,jetB);
00152         if (dist < jetI->NN_dist) {
00153           if (jetI != jetB) {
00154             jetI->NN_dist = dist;
00155             jetI->NN = jetB;
00156             diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
00157           }
00158         }
00159         if (dist < jetB->NN_dist) {
00160           if (jetI != jetB) {
00161             jetB->NN_dist = dist;
00162             jetB->NN      = jetI;}
00163         }
00164       }
00165       // if jetI's NN is the new tail then relabel it so that it becomes jetA
00166       if (jetI->NN == tail) {jetI->NN = jetA;}
00167     }
00168 
00169     
00170     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00171 
00172 
00173   }
00174   
00175 
00176   // final cleaning up;
00177   delete[] diJ;
00178   delete[] briefjets;
00179 }

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(), 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().

00264                                         {
00265 
00266   _initialise_tiles();
00267 
00268   int n = _jets.size();
00269   TiledJet * briefjets = new TiledJet[n];
00270   TiledJet * jetA = briefjets, * jetB;
00271   TiledJet oldB;
00272   
00273 
00274   // will be used quite deep inside loops, but declare it here so that
00275   // memory (de)allocation gets done only once
00276   vector<int> tile_union(3*n_tile_neighbours);
00277   
00278   // initialise the basic jet info 
00279   for (int i = 0; i< n; i++) {
00280     _tj_set_jetinfo(jetA, i);
00281     //cout << i<<": "<<jetA->tile_index<<"\n";
00282     jetA++; // move on to next entry of briefjets
00283   }
00284   TiledJet * tail = jetA; // a semaphore for the end of briefjets
00285   TiledJet * head = briefjets; // a nicer way of naming start
00286 
00287   // set up the initial nearest neighbour information
00288   vector<Tile>::const_iterator tile;
00289   for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00290     // first do it on this tile
00291     for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00292       for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00293         double dist = _bj_dist(jetA,jetB);
00294         if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00295         if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00296       }
00297     }
00298     // then do it for RH tiles
00299     for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00300       for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00301         for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00302           double dist = _bj_dist(jetA,jetB);
00303           if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00304           if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00305         }
00306       }
00307     }
00308   }
00309   
00310   // now create the diJ (where J is i's NN) table -- remember that 
00311   // we differ from standard normalisation here by a factor of R2
00312   double * diJ = new double[n];
00313   jetA = head;
00314   for (int i = 0; i < n; i++) {
00315     diJ[i] = _bj_diJ(jetA);
00316     jetA++; // have jetA follow i
00317   }
00318 
00319   // now run the recombination loop
00320   int history_location = n-1;
00321   while (tail != head) {
00322 
00323     // find the minimum of the diJ on this round
00324     double diJ_min = diJ[0];
00325     int diJ_min_jet = 0;
00326     for (int i = 1; i < n; i++) {
00327       if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min  = diJ[i];}
00328     }
00329 
00330     // do the recombination between A and B
00331     history_location++;
00332     jetA = & briefjets[diJ_min_jet];
00333     jetB = jetA->NN;
00334     // put the normalisation back in
00335     diJ_min *= _invR2; 
00336 
00337     //if (n == 19) {cout << "Hello "<<jetA-head<<" "<<jetB-head<<"\n";}
00338 
00339     //cout <<" WILL RECOMBINE "<< jetA-briefjets<<" "<<jetB-briefjets<<"\n";
00340 
00341     if (jetB != NULL) {
00342       // jet-jet recombination
00343       // If necessary relabel A & B to ensure jetB < jetA, that way if
00344       // the larger of them == newtail then that ends up being jetA and 
00345       // the new jet that is added as jetB is inserted in a position that
00346       // has a future!
00347       if (jetA < jetB) {swap(jetA,jetB);}
00348 
00349       int nn; // new jet index
00350       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00351 
00352       //OBS// get the two history indices
00353       //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00354       //OBSint hist_b = _jets[jetB->_jets_index].cluster_hist_index();
00355       //OBS// create the recombined jet
00356       //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
00357       //OBSint nn = _jets.size() - 1;
00358       //OBS_jets[nn].set_cluster_hist_index(history_location);
00359       //OBS// update history
00360       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
00361       //OBS_add_step_to_history(history_location, 
00362       //OBS                min(hist_a,hist_b),max(hist_a,hist_b),
00363       //OBS                        nn, diJ_min);
00364 
00365       // what was jetB will now become the new jet
00366       _bj_remove_from_tiles(jetA);
00367       oldB = * jetB;  // take a copy because we will need it...
00368       _bj_remove_from_tiles(jetB);
00369       _tj_set_jetinfo(jetB, nn); // also registers the jet in the tiling
00370     } else {
00371       // jet-beam recombination
00372       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00373             
00374       //OBS// get the hist_index
00375       //OBSint hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00376       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
00377       //OBS_add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); 
00378       _bj_remove_from_tiles(jetA);
00379     }
00380 
00381     // first establish the set of tiles over which we are going to 
00382     // have to run searches for updated and new nearest-neighbours
00383     int n_near_tiles = 0;
00384     _add_neighbours_to_tile_union(jetA->tile_index, tile_union, n_near_tiles);
00385     if (jetB != NULL) {
00386       bool sort_it = false;
00387       if (jetB->tile_index != jetA->tile_index) {
00388         sort_it = true;
00389         _add_neighbours_to_tile_union(jetB->tile_index,tile_union,n_near_tiles);
00390       }
00391       if (oldB.tile_index != jetA->tile_index && 
00392           oldB.tile_index != jetB->tile_index) {
00393         sort_it = true;
00394         _add_neighbours_to_tile_union(oldB.tile_index,tile_union,n_near_tiles);
00395       }
00396 
00397       if (sort_it) {
00398         // sort the tiles before then compressing the list
00399         sort(tile_union.begin(), tile_union.begin()+n_near_tiles);
00400         // and now condense the list
00401         int nnn = 1;
00402         for (int i = 1; i < n_near_tiles; i++) {
00403           if (tile_union[i] != tile_union[nnn-1]) {
00404             tile_union[nnn] = tile_union[i]; 
00405             nnn++;
00406           }
00407         }
00408         n_near_tiles = nnn;
00409       }
00410     }
00411 
00412     // now update our nearest neighbour info and diJ table
00413     // first reduce size of table
00414     tail--; n--;
00415     if (jetA == tail) {
00416       // there is nothing to be done
00417     } else {
00418       // Copy last jet contents and diJ info into position of jetA
00419       *jetA = *tail;
00420       diJ[jetA - head] = diJ[tail-head];
00421       // IN the tiling fix pointers to tail and turn them into
00422       // pointers to jetA (from predecessors, successors and the tile
00423       // head if need be)
00424       if (jetA->previous == NULL) {
00425         _tiles[jetA->tile_index].head = jetA;
00426       } else {
00427         jetA->previous->next = jetA;
00428       }
00429       if (jetA->next != NULL) {jetA->next->previous = jetA;}
00430     }
00431 
00432     // Initialise jetB's NN distance as well as updating it for 
00433     // other particles.
00434     for (int itile = 0; itile < n_near_tiles; itile++) {
00435       Tile * tile = &_tiles[tile_union[itile]];
00436       for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00437         // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00438         if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00439           jetI->NN_dist = _R2;
00440           jetI->NN      = NULL;
00441           // now go over tiles that are neighbours of I (include own tile)
00442           for (Tile ** near_tile  = tile->begin_tiles; 
00443                        near_tile != tile->end_tiles; near_tile++) {
00444             // and then over the contents of that tile
00445             for (TiledJet * jetJ  = (*near_tile)->head; 
00446                             jetJ != NULL; jetJ = jetJ->next) {
00447               double dist = _bj_dist(jetI,jetJ);
00448               if (dist < jetI->NN_dist && jetJ != jetI) {
00449                 jetI->NN_dist = dist; jetI->NN = jetJ;
00450               }
00451             }
00452           }
00453           diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
00454         }
00455         // check whether new jetB is closer than jetI's current NN and
00456         // if need to update things
00457         if (jetB != NULL) {
00458           double dist = _bj_dist(jetI,jetB);
00459           if (dist < jetI->NN_dist) {
00460             if (jetI != jetB) {
00461               jetI->NN_dist = dist;
00462               jetI->NN = jetB;
00463               diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
00464             }
00465           }
00466           if (dist < jetB->NN_dist) {
00467             if (jetI != jetB) {
00468               jetB->NN_dist = dist;
00469               jetB->NN      = jetI;}
00470           }
00471         }
00472       }
00473     }
00474 
00475 
00476     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00477     //cout << n<<" "<<briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
00478 
00479     // remember to update pointers to tail
00480     for (Tile ** near_tile = _tiles[tail->tile_index].begin_tiles; 
00481                  near_tile!= _tiles[tail->tile_index].end_tiles; near_tile++){
00482       // and then the contents of that tile
00483       for (TiledJet * jetJ = (*near_tile)->head; 
00484                      jetJ != NULL; jetJ = jetJ->next) {
00485         if (jetJ->NN == tail) {jetJ->NN = jetA;}
00486       }
00487     }
00488 
00489     //for (int i = 0; i < n; i++) {
00490     //  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";}
00491     //}
00492 
00493 
00494     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00495     //cout << briefjets[95].NN-briefjets<<" "<<briefjets[95].NN_dist <<"\n";
00496 
00497   }
00498 
00499   // final cleaning up;
00500   delete[] diJ;
00501   delete[] briefjets;
00502 }

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::TiledJet::diJ_posn, n_tile_neighbours, fastjet::ClusterSequence::TiledJet::next, fastjet::ClusterSequence::TiledJet::NN, fastjet::ClusterSequence::TiledJet::NN_dist, and fastjet::ClusterSequence::TiledJet::tile_index.

Referenced by _initialise_and_run().

00507                                                {
00508 
00509   _initialise_tiles();
00510 
00511   int n = _jets.size();
00512   TiledJet * briefjets = new TiledJet[n];
00513   TiledJet * jetA = briefjets, * jetB;
00514   TiledJet oldB;
00515   
00516 
00517   // will be used quite deep inside loops, but declare it here so that
00518   // memory (de)allocation gets done only once
00519   vector<int> tile_union(3*n_tile_neighbours);
00520   
00521   // initialise the basic jet info 
00522   for (int i = 0; i< n; i++) {
00523     _tj_set_jetinfo(jetA, i);
00524     //cout << i<<": "<<jetA->tile_index<<"\n";
00525     jetA++; // move on to next entry of briefjets
00526   }
00527   TiledJet * head = briefjets; // a nicer way of naming start
00528 
00529   // set up the initial nearest neighbour information
00530   vector<Tile>::const_iterator tile;
00531   for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00532     // first do it on this tile
00533     for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00534       for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00535         double dist = _bj_dist(jetA,jetB);
00536         if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00537         if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00538       }
00539     }
00540     // then do it for RH tiles
00541     for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00542       for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00543         for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00544           double dist = _bj_dist(jetA,jetB);
00545           if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00546           if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00547         }
00548       }
00549     }
00550     // no need to do it for LH tiles, since they are implicitly done
00551     // when we set NN for both jetA and jetB on the RH tiles.
00552   }
00553 
00554   
00555   // now create the diJ (where J is i's NN) table -- remember that 
00556   // we differ from standard normalisation here by a factor of R2
00557   // (corrected for at the end). 
00558   struct diJ_plus_link {
00559     double     diJ; // the distance
00560     TiledJet * jet; // the jet (i) for which we've found this distance
00561                     // (whose NN will the J).
00562   };
00563   diJ_plus_link * diJ = new diJ_plus_link[n];
00564   jetA = head;
00565   for (int i = 0; i < n; i++) {
00566     diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
00567     diJ[i].jet = jetA;  // our compact diJ table will not be in      
00568     jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
00569                         // so set up bi-directional correspondence here.
00570     jetA++; // have jetA follow i 
00571   }
00572 
00573   // now run the recombination loop
00574   int history_location = n-1;
00575   while (n > 0) {
00576 
00577     // find the minimum of the diJ on this round
00578     diJ_plus_link * best, *stop; // pointers a bit faster than indices
00579     // could use best to keep track of diJ min, but it turns out to be
00580     // marginally faster to have a separate variable (avoids n
00581     // dereferences at the expense of n/2 assignments).
00582     double diJ_min = diJ[0].diJ; // initialise the best one here.
00583     best = diJ;                  // and here
00584     stop = diJ+n;
00585     for (diJ_plus_link * here = diJ+1; here != stop; here++) {
00586       if (here->diJ < diJ_min) {best = here; diJ_min  = here->diJ;}
00587     }
00588 
00589     // do the recombination between A and B
00590     history_location++;
00591     jetA = best->jet;
00592     jetB = jetA->NN;
00593     // put the normalisation back in
00594     diJ_min *= _invR2; 
00595 
00596     if (jetB != NULL) {
00597       // jet-jet recombination
00598       // If necessary relabel A & B to ensure jetB < jetA, that way if
00599       // the larger of them == newtail then that ends up being jetA and 
00600       // the new jet that is added as jetB is inserted in a position that
00601       // has a future!
00602       if (jetA < jetB) {swap(jetA,jetB);}
00603 
00604       int nn; // new jet index
00605       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00606       
00607       //OBS// get the two history indices
00608       //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
00609       //OBSint ihstry_b = _jets[jetB->_jets_index].cluster_hist_index();
00610       //OBS// create the recombined jet
00611       //OBS_jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
00612       //OBSint nn = _jets.size() - 1;
00613       //OBS_jets[nn].set_cluster_hist_index(history_location);
00614       //OBS// update history
00615       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<jetB-head<<"; ";
00616       //OBS_add_step_to_history(history_location, 
00617       //OBS                min(ihstry_a,ihstry_b),max(ihstry_a,ihstry_b),
00618       //OBS                        nn, diJ_min);
00619       // what was jetB will now become the new jet
00620       _bj_remove_from_tiles(jetA);
00621       oldB = * jetB;  // take a copy because we will need it...
00622       _bj_remove_from_tiles(jetB);
00623       _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
00624                                  // (also registers the jet in the tiling)
00625     } else {
00626       // jet-beam recombination
00627       // get the hist_index
00628       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00629       //OBSint ihstry_a = _jets[jetA->_jets_index].cluster_hist_index();
00630       //OBS//cout <<n-1<<" "<<jetA-head<<" "<<-1<<"; ";
00631       //OBS_add_step_to_history(history_location,ihstry_a,BeamJet,Invalid,diJ_min); 
00632       _bj_remove_from_tiles(jetA);
00633     }
00634 
00635     // first establish the set of tiles over which we are going to
00636     // have to run searches for updated and new nearest-neighbours --
00637     // basically a combination of vicinity of the tiles of the two old
00638     // and one new jet.
00639     int n_near_tiles = 0;
00640     _add_untagged_neighbours_to_tile_union(jetA->tile_index, 
00641                                            tile_union, n_near_tiles);
00642     if (jetB != NULL) {
00643       if (jetB->tile_index != jetA->tile_index) {
00644         _add_untagged_neighbours_to_tile_union(jetB->tile_index,
00645                                                tile_union,n_near_tiles);
00646       }
00647       if (oldB.tile_index != jetA->tile_index && 
00648           oldB.tile_index != jetB->tile_index) {
00649         _add_untagged_neighbours_to_tile_union(oldB.tile_index,
00650                                                tile_union,n_near_tiles);
00651       }
00652     }
00653 
00654     // now update our nearest neighbour info and diJ table
00655     // first reduce size of table
00656     n--;
00657     // then compactify the diJ by taking the last of the diJ and copying
00658     // it to the position occupied by the diJ for jetA
00659     diJ[n].jet->diJ_posn = jetA->diJ_posn;
00660     diJ[jetA->diJ_posn] = diJ[n];
00661 
00662     // Initialise jetB's NN distance as well as updating it for 
00663     // other particles.
00664     // Run over all tiles in our union 
00665     for (int itile = 0; itile < n_near_tiles; itile++) {
00666       Tile * tile = &_tiles[tile_union[itile]];
00667       tile->tagged = false; // reset tag, since we're done with unions
00668       // run over all jets in the current tile
00669       for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00670         // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00671         if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00672           jetI->NN_dist = _R2;
00673           jetI->NN      = NULL;
00674           // now go over tiles that are neighbours of I (include own tile)
00675           for (Tile ** near_tile  = tile->begin_tiles; 
00676                        near_tile != tile->end_tiles; near_tile++) {
00677             // and then over the contents of that tile
00678             for (TiledJet * jetJ  = (*near_tile)->head; 
00679                             jetJ != NULL; jetJ = jetJ->next) {
00680               double dist = _bj_dist(jetI,jetJ);
00681               if (dist < jetI->NN_dist && jetJ != jetI) {
00682                 jetI->NN_dist = dist; jetI->NN = jetJ;
00683               }
00684             }
00685           }
00686           diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ kt-dist
00687         }
00688         // check whether new jetB is closer than jetI's current NN and
00689         // if jetI is closer than jetB's current (evolving) nearest
00690         // neighbour. Where relevant update things
00691         if (jetB != NULL) {
00692           double dist = _bj_dist(jetI,jetB);
00693           if (dist < jetI->NN_dist) {
00694             if (jetI != jetB) {
00695               jetI->NN_dist = dist;
00696               jetI->NN = jetB;
00697               diJ[jetI->diJ_posn].diJ = _bj_diJ(jetI); // update diJ...
00698             }
00699           }
00700           if (dist < jetB->NN_dist) {
00701             if (jetI != jetB) {
00702               jetB->NN_dist = dist;
00703               jetB->NN      = jetI;}
00704           }
00705         }
00706       }
00707     }
00708 
00709     // finally, register the updated kt distance for B
00710     if (jetB != NULL) {diJ[jetB->diJ_posn].diJ = _bj_diJ(jetB);}
00711 
00712   }
00713 
00714   // final cleaning up;
00715   delete[] diJ;
00716   delete[] briefjets;
00717 }

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::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(), and fastjet::ClusterSequence::TiledJet::tile_index.

Referenced by _initialise_and_run().

00724                                                        {
00725 
00726   _initialise_tiles();
00727 
00728   int n = _jets.size();
00729   TiledJet * briefjets = new TiledJet[n];
00730   TiledJet * jetA = briefjets, * jetB;
00731   TiledJet oldB;
00732   
00733 
00734   // will be used quite deep inside loops, but declare it here so that
00735   // memory (de)allocation gets done only once
00736   vector<int> tile_union(3*n_tile_neighbours);
00737   
00738   // initialise the basic jet info 
00739   for (int i = 0; i< n; i++) {
00740     _tj_set_jetinfo(jetA, i);
00741     //cout << i<<": "<<jetA->tile_index<<"\n";
00742     jetA++; // move on to next entry of briefjets
00743   }
00744   TiledJet * head = briefjets; // a nicer way of naming start
00745 
00746   // set up the initial nearest neighbour information
00747   vector<Tile>::const_iterator tile;
00748   for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
00749     // first do it on this tile
00750     for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00751       for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
00752         double dist = _bj_dist(jetA,jetB);
00753         if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00754         if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00755       }
00756     }
00757     // then do it for RH tiles
00758     for (Tile ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
00759       for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
00760         for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
00761           double dist = _bj_dist(jetA,jetB);
00762           if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
00763           if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
00764         }
00765       }
00766     }
00767     // no need to do it for LH tiles, since they are implicitly done
00768     // when we set NN for both jetA and jetB on the RH tiles.
00769   }
00770 
00771   
00775   //struct diJ_plus_link {
00776   //  double     diJ; // the distance
00777   //  TiledJet * jet; // the jet (i) for which we've found this distance
00778   //                  // (whose NN will the J).
00779   //};
00780   //diJ_plus_link * diJ = new diJ_plus_link[n];
00781   //jetA = head;
00782   //for (int i = 0; i < n; i++) {
00783   //  diJ[i].diJ = _bj_diJ(jetA); // kt distance * R^2
00784   //  diJ[i].jet = jetA;  // our compact diJ table will not be in            
00785   //  jetA->diJ_posn = i; // one-to-one corresp. with non-compact jets,
00786   //                      // so set up bi-directional correspondence here.
00787   //  jetA++; // have jetA follow i 
00788   //}
00789 
00790   vector<double> diJs(n);
00791   for (int i = 0; i < n; i++) {
00792     diJs[i] = _bj_diJ(&briefjets[i]);
00793     briefjets[i].label_minheap_update_done();
00794   }
00795   MinHeap minheap(diJs);
00796   // have a stack telling us which jets we'll have to update on the heap
00797   vector<TiledJet *> jets_for_minheap;
00798   jets_for_minheap.reserve(n); 
00799 
00800   // now run the recombination loop
00801   int history_location = n-1;
00802   while (n > 0) {
00803 
00804     double diJ_min = minheap.minval() *_invR2;
00805     jetA = head + minheap.minloc();
00806 
00807     // do the recombination between A and B
00808     history_location++;
00809     jetB = jetA->NN;
00810 
00811     if (jetB != NULL) {
00812       // jet-jet recombination
00813       // If necessary relabel A & B to ensure jetB < jetA, that way if
00814       // the larger of them == newtail then that ends up being jetA and 
00815       // the new jet that is added as jetB is inserted in a position that
00816       // has a future!
00817       if (jetA < jetB) {swap(jetA,jetB);}
00818 
00819       int nn; // new jet index
00820       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00821       
00822       // what was jetB will now become the new jet
00823       _bj_remove_from_tiles(jetA);
00824       oldB = * jetB;  // take a copy because we will need it...
00825       _bj_remove_from_tiles(jetB);
00826       _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
00827                                  // (also registers the jet in the tiling)
00828     } else {
00829       // jet-beam recombination
00830       // get the hist_index
00831       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00832       _bj_remove_from_tiles(jetA);
00833     }
00834 
00835     // remove the minheap entry for jetA
00836     minheap.remove(jetA-head);
00837 
00838     // first establish the set of tiles over which we are going to
00839     // have to run searches for updated and new nearest-neighbours --
00840     // basically a combination of vicinity of the tiles of the two old
00841     // and one new jet.
00842     int n_near_tiles = 0;
00843     _add_untagged_neighbours_to_tile_union(jetA->tile_index, 
00844                                            tile_union, n_near_tiles);
00845     if (jetB != NULL) {
00846       if (jetB->tile_index != jetA->tile_index) {
00847         _add_untagged_neighbours_to_tile_union(jetB->tile_index,
00848                                                tile_union,n_near_tiles);
00849       }
00850       if (oldB.tile_index != jetA->tile_index && 
00851           oldB.tile_index != jetB->tile_index) {
00852         _add_untagged_neighbours_to_tile_union(oldB.tile_index,
00853                                                tile_union,n_near_tiles);
00854       }
00855       // indicate that we'll have to update jetB in the minheap
00856       jetB->label_minheap_update_needed();
00857       jets_for_minheap.push_back(jetB);
00858     }
00859 
00860 
00861     // Initialise jetB's NN distance as well as updating it for 
00862     // other particles.
00863     // Run over all tiles in our union 
00864     for (int itile = 0; itile < n_near_tiles; itile++) {
00865       Tile * tile = &_tiles[tile_union[itile]];
00866       tile->tagged = false; // reset tag, since we're done with unions
00867       // run over all jets in the current tile
00868       for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00869         // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00870         if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
00871           jetI->NN_dist = _R2;
00872           jetI->NN      = NULL;
00873           // label jetI as needing heap action...
00874           if (!jetI->minheap_update_needed()) {
00875             jetI->label_minheap_update_needed();
00876             jets_for_minheap.push_back(jetI);}
00877           // now go over tiles that are neighbours of I (include own tile)
00878           for (Tile ** near_tile  = tile->begin_tiles; 
00879                        near_tile != tile->end_tiles; near_tile++) {
00880             // and then over the contents of that tile
00881             for (TiledJet * jetJ  = (*near_tile)->head; 
00882                             jetJ != NULL; jetJ = jetJ->next) {
00883               double dist = _bj_dist(jetI,jetJ);
00884               if (dist < jetI->NN_dist && jetJ != jetI) {
00885                 jetI->NN_dist = dist; jetI->NN = jetJ;
00886               }
00887             }
00888           }
00889         }
00890         // check whether new jetB is closer than jetI's current NN and
00891         // if jetI is closer than jetB's current (evolving) nearest
00892         // neighbour. Where relevant update things
00893         if (jetB != NULL) {
00894           double dist = _bj_dist(jetI,jetB);
00895           if (dist < jetI->NN_dist) {
00896             if (jetI != jetB) {
00897               jetI->NN_dist = dist;
00898               jetI->NN = jetB;
00899               // label jetI as needing heap action...
00900               if (!jetI->minheap_update_needed()) {
00901                 jetI->label_minheap_update_needed();
00902                 jets_for_minheap.push_back(jetI);}
00903             }
00904           }
00905           if (dist < jetB->NN_dist) {
00906             if (jetI != jetB) {
00907               jetB->NN_dist = dist;
00908               jetB->NN      = jetI;}
00909           }
00910         }
00911       }
00912     }
00913 
00914     // deal with jets whose minheap entry needs updating
00915     while (jets_for_minheap.size() > 0) {
00916       TiledJet * jetI = jets_for_minheap.back(); 
00917       jets_for_minheap.pop_back();
00918       minheap.update(jetI-head, _bj_diJ(jetI));
00919       jetI->label_minheap_update_done();
00920     }
00921     n--;
00922   }
00923 
00924   // final cleaning up;
00925   delete[] briefjets;
00926 }

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_algorithm, _jets, BeamJet, fastjet::cambridge_algorithm, fastjet::ClosestPair2D::closest_pair(), Invalid, fastjet::ClosestPair2D::replace(), and fastjet::twopi.

Referenced by _initialise_and_run().

00223                                          {
00224 
00225   if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
00226 
00227   unsigned int n = _jets.size();
00228 
00229   vector<MirrorInfo>   coordIDs(2*n);  // link from original to mirror indices
00230   vector<int>          jetIDs(2*n);     // link from mirror to original indices
00231   vector<Coord2D>      coords(2*n);   // our coordinates (and copies)
00232 
00233   // start things off...
00234   double minrap = numeric_limits<double>::max();
00235   double maxrap = -minrap;
00236   int coord_index = 0;
00237   for (unsigned i = 0; i < n; i++) {
00238     // separate out points with infinite rapidity
00239     if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
00240       coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
00241     } else {
00242       coordIDs[i].orig   = coord_index;
00243       coordIDs[i].mirror = coord_index+1;
00244       coords[coord_index]   = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
00245       coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
00246       jetIDs[coord_index]   = i;
00247       jetIDs[coord_index+1] = i;
00248       minrap = min(coords[coord_index].x,minrap);
00249       maxrap = max(coords[coord_index].x,maxrap);
00250       coord_index += 2;
00251     }
00252   }
00253   // label remaining "mirror info" as saying that there are no jets
00254   for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
00255 
00256   // set them to their actual size...
00257   coords.resize(coord_index);
00258 
00259   // establish limits (with some leeway on rapidity)
00260   Coord2D left_edge(minrap-1.0, 0.0);
00261   Coord2D right_edge(maxrap+1.0, 2*twopi);
00262 
00263   // now create the closest pair search object
00264   ClosestPair2D cp(coords, left_edge, right_edge);
00265 
00266   // and start iterating...
00267   vector<Coord2D> new_points(2);
00268   vector<unsigned int> cIDs_to_remove(4);
00269   vector<unsigned int> new_cIDs(2);
00270   do {
00271     // find out which pair we recombine
00272     unsigned int cID1, cID2;
00273     double distance2;
00274     cp.closest_pair(cID1,cID2,distance2);
00275     distance2 *= _invR2;
00276 
00277     // if the closest distance > 1, we exit and go to "inclusive" stage
00278     if (distance2 > 1.0) {break;}
00279 
00280     // do the recombination...
00281     int jet_i = jetIDs[cID1];
00282     int jet_j = jetIDs[cID2];
00283     int newjet_k;
00284     _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
00285 
00286     // now prepare operations on CP structure
00287     cIDs_to_remove[0] = coordIDs[jet_i].orig;
00288     cIDs_to_remove[1] = coordIDs[jet_i].mirror;
00289     cIDs_to_remove[2] = coordIDs[jet_j].orig;
00290     cIDs_to_remove[3] = coordIDs[jet_j].mirror;
00291     new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
00292     new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
00293     // carry out the CP operation
00294     //cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
00295     // remarkable the following is faster...
00296     new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
00297     new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
00298     // signal that the following jets no longer exist
00299     coordIDs[jet_i].orig = Invalid;
00300     coordIDs[jet_j].orig = Invalid;
00301     // and do bookkeeping for new jet
00302     coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
00303     jetIDs[new_cIDs[0]] = newjet_k;
00304     jetIDs[new_cIDs[1]] = newjet_k;
00305 
00306     // if we've reached one jet we should exit...
00307     n--;
00308     if (n == 1) {break;}
00309 
00310   } while(true);
00311   
00312 
00313   // now do "beam" recombinations 
00314   //for (unsigned int jet_i = 0; jet_i < coordIDs.size(); jet_i++) {
00315   //  // if we have a predeclared beam jet or a valid set of mirror
00316   //  // coordinates, recombine then recombine this jet with the beam
00317   //  if (coordIDs[jet_i].orig == BeamJet || coordIDs[jet_i].orig > 0) {
00318   //    // diB = 1 _always_ (for the cambridge alg.)
00319   //    _do_iB_recombination_step(jet_i, 1.0);
00320   //  }
00321   //}
00322 
00323   _do_Cambridge_inclusive_jets();
00324 
00325   //n = _history.size();
00326   //for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
00327   //  if (_history[hist_i].child == Invalid) {
00328   //    _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
00329   //  }
00330   //}
00331 
00332 }

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_algorithm, _Rparam, and fastjet::cambridge_algorithm.

Referenced by _CP2DChan_cluster_2piMultD(), and _initialise_and_run().

00193                                                {
00194 
00195   if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
00196 
00197   // run the clustering with mirror copies kept such that only things
00198   // within _Rparam of a border are mirrored
00199   _CP2DChan_limited_cluster(_Rparam);
00200 
00201   //
00202   _do_Cambridge_inclusive_jets();
00203 }

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().

00209                                                   {
00210 
00211   // do a first run of clustering up to a small distance parameter,
00212   if (_Rparam >= 0.39) {
00213     _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
00214   }
00215 
00216   // and then the final round of clustering
00217   _CP2DChan_cluster_2pi2R ();
00218 }

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().

00069                                                             {
00070   
00071   unsigned int n = _initial_n;
00072 
00073   vector<MirrorInfo>   coordIDs(2*n); // coord IDs of a given jetID
00074   vector<int>          jetIDs(2*n);   // jet ID for a given coord ID
00075   vector<Coord2D>      coords(2*n);   // our coordinates (and copies)
00076 
00077   // start things off...
00078   double minrap = numeric_limits<double>::max();
00079   double maxrap = -minrap;
00080   int coord_index = -1;
00081   int n_active = 0;
00082   for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
00083 
00084     // skip jets that already have children or that have infinite
00085     // rapidity
00086     if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
00087         (_jets[jet_i].E() == abs(_jets[jet_i].pz()) && 
00088          _jets[jet_i].perp2() == 0.0)
00089         ) {continue;}
00090 
00091     n_active++;
00092 
00093     coordIDs[jet_i].orig = ++coord_index;
00094     coords[coord_index]  = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
00095     jetIDs[coord_index]  = jet_i;
00096     minrap = min(coords[coord_index].x,minrap);
00097     maxrap = max(coords[coord_index].x,maxrap);
00098 
00099     Coord2D mirror_point(coords[coord_index]);
00100     if (make_mirror(mirror_point, Dlim)) {
00101       coordIDs[jet_i].mirror = ++coord_index;
00102       coords[coord_index] = mirror_point;
00103       jetIDs[coord_index] = jet_i;
00104     } else {
00105       coordIDs[jet_i].mirror = Invalid;
00106     }
00107   }
00108 
00109   // set them to their actual size...
00110   coords.resize(coord_index+1);
00111 
00112   // establish limits (with some leeway on rapidity)
00113   Coord2D left_edge(minrap-1.0, -pi);
00114   Coord2D right_edge(maxrap+1.0, 3*pi);
00115 
00116   //cerr << "minrap, maxrap = " << minrap << " " << maxrap << endl;
00117 
00118   // now create the closest pair search object
00119   ClosestPair2D cp(coords, left_edge, right_edge);
00120 
00121   // cross check to see what's going on...
00122   //cerr << "Tree depths before: ";
00123   //cp.print_tree_depths(cerr);
00124 
00125   // and start iterating...
00126   vector<Coord2D> new_points(2);
00127   vector<unsigned int> cIDs_to_remove(4);
00128   vector<unsigned int> new_cIDs(2);
00129 
00130   do {
00131     // find out which pair we recombine
00132     unsigned int cID1, cID2;
00133     double distance2;
00134     cp.closest_pair(cID1,cID2,distance2);
00135 
00136     // if the closest distance > Dlim, we exit and go to "inclusive" stage
00137     if (distance2 > Dlim*Dlim) {break;}
00138 
00139     // normalise distance as we like it
00140     distance2 *= _invR2;
00141 
00142     // do the recombination...
00143     int jet_i = jetIDs[cID1];
00144     int jet_j = jetIDs[cID2];
00145     int newjet_k;
00146     _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
00147 
00148     // don't bother with any further action if only one active particle
00149     // is left (also avoid closest-pair error [cannot remove last particle]).
00150     if (--n_active == 1) {break;}
00151 
00152     // now prepare operations on CP structure
00153     cIDs_to_remove.resize(0);
00154     cIDs_to_remove.push_back(coordIDs[jet_i].orig);
00155     cIDs_to_remove.push_back(coordIDs[jet_j].orig);
00156     if (coordIDs[jet_i].mirror != Invalid) 
00157       cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
00158     if (coordIDs[jet_j].mirror != Invalid) 
00159       cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
00160 
00161     Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
00162     new_points.resize(0);
00163     new_points.push_back(new_point);
00164     if (make_mirror(new_point, Dlim)) new_points.push_back(new_point);
00165     
00166     // carry out actions on search tree
00167     cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
00168 
00169     // now fill in info for new points...
00170     coordIDs[newjet_k].orig = new_cIDs[0];
00171     jetIDs[new_cIDs[0]]       = newjet_k;
00172     if (new_cIDs.size() == 2) {
00173       coordIDs[newjet_k].mirror = new_cIDs[1];
00174       jetIDs[new_cIDs[1]]         = newjet_k;
00175     } else {coordIDs[newjet_k].mirror = Invalid;}
00176     
00178     //n_active--;
00179     //if (n_active == 1) {break;}
00180 
00181   } while(true);
00182   
00183   // cross check to see what's going on...
00184   //cerr << "Max tree depths after: ";
00185   //cp.print_tree_depths(cerr);
00186 
00187 }

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().

00336                                                     {
00337   unsigned int n = _history.size();
00338   for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
00339     if (_history[hist_i].child == Invalid) {
00340       _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
00341     }
00342   }
00343 }

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 715 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().

00718                                    {
00719 
00720   history_element element;
00721   element.parent1 = parent1;
00722   element.parent2 = parent2;
00723   element.jetp_index = jetp_index;
00724   element.child = Invalid;
00725   element.dij   = dij;
00726   element.max_dij_so_far = max(dij,_history[_history.size()-1].max_dij_so_far);
00727   _history.push_back(element);
00728 
00729   int local_step = _history.size()-1;
00730   assert(local_step == step_number);
00731 
00732   assert(parent1 >= 0);
00733   _history[parent1].child = local_step;
00734   if (parent2 >= 0) {_history[parent2].child = local_step;}
00735 
00736   // get cross-referencing right from PseudoJets
00737   if (jetp_index != Invalid) {
00738     assert(jetp_index >= 0);
00739     //cout << _jets.size() <<" "<<jetp_index<<"\n";
00740     _jets[jetp_index].set_cluster_hist_index(local_step);
00741   }
00742 
00743   if (_writeout_combinations) {
00744     cout << local_step << ": " 
00745          << parent1 << " with " << parent2
00746          << "; y = "<< dij<<endl;
00747   }
00748 
00749 }

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)

Referenced by unique_history_order().

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)

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().

00223                                                                 {
00224   
00225   double yiB = jet_scale_for_algorithm(_jets[ii]);
00226   if (yiB == 0.0) {
00227     // in this case convention is that we do not worry about distances
00228     // but directly state that nearest neighbour is beam
00229     DijMap.insert(DijEntry(yiB,  TwoVertices(ii,-1)));
00230   } else {
00231     double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
00232     // Logic of following bit is: only add point to map if it has
00233     // smaller kt2 than nearest neighbour j (if it has larger kt,
00234     // then: either it is j's nearest neighbour and then we will
00235     // include dij when we come to j; or it is not j's nearest
00236     // neighbour and j will recombine with someone else).
00237     
00238     // If DeltaR2 > 1.0 then in any case it will recombine with beam rather
00239     // than with any neighbours.
00240     // (put general normalisation here at some point)
00241     if (DeltaR2 > 1.0) {
00242       DijMap.insert(DijEntry(yiB,  TwoVertices(ii,-1)));
00243     } else {
00244       double kt2i = jet_scale_for_algorithm(_jets[ii]);
00245       int jj = DNN->NearestNeighbourIndex(ii);
00246       if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
00247         double dij = DeltaR2 * kt2i;
00248         DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
00249       }
00250     }
00251   }
00252 }

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

for making sure the user knows what it is they're running...

Definition at line 153 of file ClusterSequence.cc.

References _first_time, and fastjet_version.

Referenced by _decant_options().

00153                                     {
00154 
00155   if (!_first_time) {return;}
00156   _first_time = false;
00157   
00158   
00159   //Symp. Discr. Alg, p.472 (2002) and  CGAL (http://www.cgal.org);
00160 
00161   cout << "#--------------------------------------------------------------------------\n";
00162   cout << "#                      FastJet release " << fastjet_version << endl;
00163   cout << "#            Written by M. Cacciari, G.P. Salam and G. Soyez            \n"; 
00164   cout << "#              http://www.lpthe.jussieu.fr/~salam/fastjet               \n"; 
00165   cout << "#                                                                       \n";
00166   cout << "# Longitudinally invariant Kt, anti-Kt, and inclusive Cambridge/Aachen  \n";
00167   cout << "# clustering using fast geometric algorithms, with area measures and optional\n";
00168   cout << "# external jet-finder plugins.                                          \n";
00169   cout << "# Please cite Phys. Lett. B641 (2006) [hep-ph/0512210] if you use this code.\n";
00170   cout << "#                                                                       \n";
00171   cout << "# This package uses T.Chan's closest pair algorithm, Proc.13th ACM-SIAM \n";
00172   cout << "# Symp. Discr. Alg, p.472 (2002), S.Fortune's Voronoi algorithm and code " ;
00173 #ifndef DROP_CGAL
00174   cout << endl << "# and CGAL: http://www.cgal.org/";
00175 #endif  // DROP_CGAL
00176   cout << ".\n";
00177   cout << "#-------------------------------------------------------------------------\n";
00178 }

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 666 of file ClusterSequence.hh.

References _jets, _R2, and jet_scale_for_algorithm().

Referenced by _simple_N2_cluster().

00667                                                                          {
00668     jetA->eta  = _jets[_jets_index].rap();
00669     jetA->phi  = _jets[_jets_index].phi_02pi();
00670     jetA->kt2  = jet_scale_for_algorithm(_jets[_jets_index]);
00671     jetA->_jets_index = _jets_index;
00672     // initialise NN info as well
00673     jetA->NN_dist = _R2;
00674     jetA->NN      = NULL;
00675 }

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>

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 681 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().

00682                                                                   {
00683   double dphi = std::abs(jetA->phi - jetB->phi);
00684   double deta = (jetA->eta - jetB->eta);
00685   if (dphi > pi) {dphi = twopi - dphi;}
00686   return dphi*dphi + deta*deta;
00687 }

template<class J>

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

Definition at line 690 of file ClusterSequence.hh.

Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), _simple_N2_cluster(), and _tiled_N2_cluster().

00690                                                                                    {
00691   double kt2 = jet->kt2;
00692   if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
00693   return jet->NN_dist * kt2;
00694 }

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 515 of file ClusterSequence.hh.

00518           {
00519     J * res;
00520     for(res = head; res<tail; res++) {
00521       if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
00522     }
00523     return res;
00524   }

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 700 of file ClusterSequence.hh.

References _bj_dist(), and _R2.

Referenced by _simple_N2_cluster().

00701                                                                             {
00702   double NN_dist = _R2;
00703   J * NN  = NULL;
00704   if (head < jet) {
00705     for (J * jetB = head; jetB != jet; jetB++) {
00706       double dist = _bj_dist(jet,jetB);
00707       if (dist < NN_dist) {
00708         NN_dist = dist;
00709         NN = jetB;
00710       }
00711     }
00712   }
00713   if (tail > jet) {
00714     for (J * jetB = jet+1; jetB != tail; jetB++) {
00715       double dist = _bj_dist(jet,jetB);
00716       if (dist < NN_dist) {
00717         NN_dist = dist;
00718         NN = jetB;
00719       }
00720     }
00721   }
00722   jet->NN = NN;
00723   jet->NN_dist = NN_dist;
00724 }

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 728 of file ClusterSequence.hh.

References _bj_dist(), and _R2.

Referenced by _simple_N2_cluster().

00729                                                                 {
00730   double NN_dist = _R2;
00731   J * NN  = NULL;
00732   for (J * jetB = head; jetB != tail; jetB++) {
00733     double dist = _bj_dist(jet,jetB);
00734     if (dist < NN_dist) {
00735       NN_dist = dist;
00736       NN = jetB;
00737     }
00738     if (dist < jetB->NN_dist) {
00739       jetB->NN_dist = dist;
00740       jetB->NN = jet;
00741     }
00742   }
00743   jet->NN = NN;
00744   jet->NN_dist = NN_dist;
00745 }

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

Definition at line 573 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tj_set_jetinfo().

00573                                                     {
00574     // note that (-1)%n = -1 so that we have to add _n_tiles_phi
00575     // before performing modulo operation
00576     return (ieta-_tiles_ieta_min)*_n_tiles_phi
00577                   + (iphi+_n_tiles_phi) % _n_tiles_phi;
00578   }

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.

00165                                                                              {
00166   int ieta, iphi;
00167   if      (eta <= _tiles_eta_min) {ieta = 0;}
00168   else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
00169   else {
00170     //ieta = int(floor((eta - _tiles_eta_min) / _tile_size_eta));
00171     ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
00172     // following needed in case of rare but nasty rounding errors
00173     if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
00174       ieta = _tiles_ieta_max-_tiles_ieta_min;} 
00175   }
00176   // allow for some extent of being beyond range in calculation of phi
00177   // as well
00178   //iphi = (int(floor(phi/_tile_size_phi)) + _n_tiles_phi) % _n_tiles_phi;
00179   // with just int and no floor, things run faster but beware
00180   iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
00181   return (iphi + ieta * _n_tiles_phi);
00182 }

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, fastjet::ClusterSequence::Tile::head, fastjet::ClusterSequence::TiledJet::next, and fastjet::ClusterSequence::TiledJet::previous.

Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster().

00189                                                                      {
00190   // first call the generic setup
00191   _bj_set_jetinfo<>(jet, _jets_index);
00192 
00193   // Then do the setup specific to the tiled case.
00194 
00195   // Find out which tile it belonds to
00196   jet->tile_index = _tile_index(jet->eta, jet->phi);
00197 
00198   // Insert it into the tile's linked list of jets
00199   Tile * tile = &_tiles[jet->tile_index];
00200   jet->previous   = NULL;
00201   jet->next       = tile->head;
00202   if (jet->next != NULL) {jet->next->previous = jet;}
00203   tile->head      = jet;
00204 }

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.

00049                                                                 {
00050   Tile * tile = & _tiles[jet->tile_index];
00051 
00052   if (jet->previous == NULL) {
00053     // we are at head of the tile, so reset it.
00054     // If this was the only jet on the tile then tile->head will now be NULL
00055     tile->head = jet->next;
00056   } else {
00057     // adjust link from previous jet in this tile
00058     jet->previous->next = jet->next;
00059   }
00060   if (jet->next != NULL) {
00061     // adjust backwards-link from next jet in this tile
00062     jet->next->previous = jet->previous;
00063   }
00064 }

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::end_tiles, 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().

00085                                         {
00086 
00087   // first decide tile sizes
00088   _tile_size_eta = _Rparam;
00089   _n_tiles_phi   = int(floor(twopi/_Rparam));
00090   _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
00091 
00092   // always include zero rapidity in the tiling region
00093   _tiles_eta_min = 0.0;
00094   _tiles_eta_max = 0.0;
00095   // but go no further than following
00096   const double maxrap = 7.0;
00097 
00098   // and find out how much further one should go
00099   for(unsigned int i = 0; i < _jets.size(); i++) {
00100     double eta = _jets[i].rap();
00101     // first check if eta is in range -- to avoid taking into account
00102     // very spurious rapidities due to particles with near-zero kt.
00103     if (abs(eta) < maxrap) {
00104       if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
00105       if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
00106     }
00107   }
00108 
00109   // now adjust the values
00110   _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
00111   _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
00112   _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
00113   _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
00114 
00115   // allocate the tiles
00116   _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
00117 
00118   // now set up the cross-referencing between tiles
00119   for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
00120     for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
00121       Tile * tile = & _tiles[_tile_index(ieta,iphi)];
00122       // no jets in this tile yet
00123       tile->head = NULL; // first element of tiles points to itself
00124       tile->begin_tiles[0] =  tile;
00125       Tile ** pptile = & (tile->begin_tiles[0]);
00126       pptile++;
00127       //
00128       // set up L's in column to the left of X
00129       tile->surrounding_tiles = pptile;
00130       if (ieta > _tiles_ieta_min) {
00131         // with the itile subroutine, we can safely run tiles from
00132         // idphi=-1 to idphi=+1, because it takes care of
00133         // negative and positive boundaries
00134         for (int idphi = -1; idphi <=+1; idphi++) {
00135           *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
00136           pptile++;
00137         }       
00138       }
00139       // now set up last L (below X)
00140       *pptile = & _tiles[_tile_index(ieta,iphi-1)];
00141       pptile++;
00142       // set up first R (above X)
00143       tile->RH_tiles = pptile;
00144       *pptile = & _tiles[_tile_index(ieta,iphi+1)];
00145       pptile++;
00146       // set up remaining R's, to the right of X
00147       if (ieta < _tiles_ieta_max) {
00148         for (int idphi = -1; idphi <= +1; idphi++) {
00149           *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
00150           pptile++;
00151         }       
00152       }
00153       // now put semaphore for end tile
00154       tile->end_tiles = pptile;
00155       // finally make sure tiles are untagged
00156       tile->tagged = false;
00157     }
00158   }
00159 
00160 }

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.

00209                                                               {
00210   for (vector<Tile>::const_iterator tile = _tiles.begin(); 
00211        tile < _tiles.end(); tile++) {
00212     cout << "Tile " << tile - _tiles.begin()<<" = ";
00213     vector<int> list;
00214     for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
00215       list.push_back(jetI-briefjets);
00216       //cout <<" "<<jetI-briefjets;
00217     }
00218     sort(list.begin(),list.end());
00219     for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
00220     cout <<"\n";
00221   }
00222 }

void fastjet::ClusterSequence::_add_neighbours_to_tile_union	(	const int	tile_index,
		std::vector< int > &	tile_union,
		int &	n_near_tiles
	)			const `[private]`

Referenced by _tiled_N2_cluster().

void fastjet::ClusterSequence::_add_untagged_neighbours_to_tile_union	(	const int	tile_index,
		std::vector< int > &	tile_union,
		int &	n_near_tiles
	)			`[private]`

Referenced by _faster_tiled_N2_cluster(), and _minheap_faster_tiled_N2_cluster().

JetAlgorithm fastjet::ClusterSequence::_default_jet_algorithm = kt_algorithm [static, protected]

Definition at line 343 of file ClusterSequence.hh.

Referenced by _initialise_and_run().

JetDefinition fastjet::ClusterSequence::_jet_def [protected]

Definition at line 344 of file ClusterSequence.hh.

Referenced by _decant_options(), _do_ij_recombination_step(), _fill_initial_history(), _initialise_and_run(), fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), fastjet::ClusterSequenceVoronoiArea::_initializeVA(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), exclusive_jets(), jet_scale_for_algorithm(), and transfer_from_sequence().

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 395 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::ClusterSequencePassiveArea::_initialise_and_run_PA(), _initialise_tiles(), fastjet::ClusterSequenceVoronoiArea::_initializeVA(), _minheap_faster_tiled_N2_cluster(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), _really_dumb_cluster(), fastjet::ClusterSequenceActiveArea::_resize_and_zero_AA(), fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), fastjet::ClusterSequenceActiveArea::_run_AA(), _simple_N2_cluster(), _tiled_N2_cluster(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), _transfer_input_jets(), exclusive_jets(), has_child(), has_parents(), has_partner(), inclusive_jets(), jets(), plugin_record_ij_recombination(), transfer_from_sequence(), and unclustered_particles().

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 401 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(), _fill_initial_history(), fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), fastjet::ClusterSequenceVoronoiArea::_initializeVA(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::_post_process(), fastjet::ClusterSequenceActiveArea::_transfer_areas(), fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(), exclusive_dmerge(), exclusive_dmerge_max(), exclusive_jets(), has_child(), has_parents(), has_partner(), history(), inclusive_jets(), n_exclusive_jets(), transfer_from_sequence(), unclustered_particles(), and unique_history_order().

bool fastjet::ClusterSequence::_writeout_combinations [protected]

Definition at line 403 of file ClusterSequence.hh.

Referenced by _add_step_to_history(), _decant_options(), and transfer_from_sequence().

int fastjet::ClusterSequence::_initial_n [protected]

Definition at line 404 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(), n_particles(), and transfer_from_sequence().

double fastjet::ClusterSequence::_Rparam [protected]

Definition at line 405 of file ClusterSequence.hh.

Referenced by _CP2DChan_cluster_2pi2R(), _CP2DChan_cluster_2piMultD(), _decant_options(), _delaunay_cluster(), _initialise_and_run(), _initialise_tiles(), and transfer_from_sequence().

double fastjet::ClusterSequence::_R2 [protected]

Definition at line 405 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(), _tiled_N2_cluster(), and transfer_from_sequence().

double fastjet::ClusterSequence::_invR2 [protected]

Definition at line 405 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(), _tiled_N2_cluster(), and transfer_from_sequence().

Strategy fastjet::ClusterSequence::_strategy [protected]

Definition at line 406 of file ClusterSequence.hh.

Referenced by _decant_options(), _delaunay_cluster(), _initialise_and_run(), fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history(), strategy_string(), and transfer_from_sequence().

JetAlgorithm fastjet::ClusterSequence::_jet_algorithm [protected]

Definition at line 407 of file ClusterSequence.hh.

Referenced by _CP2DChan_cluster(), _CP2DChan_cluster_2pi2R(), _decant_options(), _initialise_and_run(), inclusive_jets(), jet_scale_for_algorithm(), and transfer_from_sequence().

bool fastjet::ClusterSequence::_plugin_activated [private]

Definition at line 411 of file ClusterSequence.hh.

Referenced by _decant_options(), _initialise_and_run(), and transfer_from_sequence().

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

Definition at line 412 of file ClusterSequence.hh.

Referenced by transfer_from_sequence().

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

will be set by default to be true for the first run

Definition at line 459 of file ClusterSequence.hh.

Referenced by _print_banner().

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 464 of file ClusterSequence.hh.

Referenced by exclusive_jets().

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 548 of file ClusterSequence.hh.

Referenced by _faster_tiled_N2_cluster(), _minheap_faster_tiled_N2_cluster(), and _tiled_N2_cluster().

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

Definition at line 566 of file ClusterSequence.hh.

Referenced by _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_min [private]

Definition at line 567 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

double fastjet::ClusterSequence::_tiles_eta_max [private]

Definition at line 567 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

double fastjet::ClusterSequence::_tile_size_eta [private]

Definition at line 568 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

double fastjet::ClusterSequence::_tile_size_phi [private]

Definition at line 568 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

int fastjet::ClusterSequence::_n_tiles_phi [private]

Definition at line 569 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

int fastjet::ClusterSequence::_tiles_ieta_min [private]

Definition at line 569 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

int fastjet::ClusterSequence::_tiles_ieta_max [private]

Definition at line 569 of file ClusterSequence.hh.

Referenced by _initialise_tiles(), and _tile_index().

fastjet::ClusterSequence Class Reference

Public Types

Public Member Functions

Static Public Member Functions

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Types

Private Member Functions

Private Attributes

Static Private Attributes

Classes

Detailed Description

Member Typedef Documentation

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


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.
bool	object_in_jet (const PseudoJet &object, const PseudoJet &jet) const
	returns true iff the object is included in the jet.
bool	has_parents (const PseudoJet &jet, PseudoJet &parent1, PseudoJet &parent2) const
	if the jet has parents in the clustering, it returns true and sets parent1 and parent2 equal to them.
bool	has_child (const PseudoJet &jet, PseudoJet &child) const
	if the jet has a child then return true and give the child jet otherwise return false and set the child to zero
bool	has_child (const PseudoJet &jet, const PseudoJet *&childp) const
	Version of has_child that sets a pointer to the child if the child exists;.
bool	has_partner (const PseudoJet &jet, PseudoJet &partner) const
	if this jet has a child (and so a partner) return true and give the partner, otherwise return false and set the partner to zero
std::vector< PseudoJet >	constituents (const PseudoJet &jet) const
	return a vector of the particles that make up jet
void	print_jets_for_root (const std::vector< PseudoJet > &jets, std::ostream &ostr=std::cout) const
	output the supplied vector of jets in a format that can be read by an appropriate root script; the format is: jet-n jet-px jet-py jet-pz jet-E particle-n particle-rap particle-phi particle-pt particle-n particle-rap particle-phi particle-pt
void	add_constituents (const PseudoJet &jet, std::vector< PseudoJet > &subjet_vector) const
	add on to subjet_vector the subjets of jet (for internal use mainly)
Strategy	strategy_used () const
	return the enum value of the strategy used to cluster the event
std::string	strategy_string () const
const JetDefinition &	jet_def () const
	return a reference to the jet definition
double	jet_scale_for_algorithm (const PseudoJet &jet) const
	returns the scale associated with a jet as required for this clustering algorithm (kt^2 for the kt-algorithm, 1 for the Cambridge algorithm).
void	plugin_record_ij_recombination (int jet_i, int jet_j, double dij, int &newjet_k)
	record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
void	plugin_record_ij_recombination (int jet_i, int jet_j, double dij, const PseudoJet &newjet, int &newjet_k)
	as for the simpler variant of plugin_record_ij_recombination, except that the new jet is attributed the momentum and user_index of newjet
void	plugin_record_iB_recombination (int jet_i, double diB)
	record the fact that there has been a recombination between jets()[jet_i] and the beam, with the specified diB; when looking for inclusive jets, any iB recombination will returned to the user as a jet.
void	plugin_associate_extras (std::auto_ptr< Extras > extras_in)
	the plugin can associated some extra information with the ClusterSequence object by calling this function
bool	plugin_activated () const
	returns true when the plugin is allowed to run the show.
const Extras *	extras () const
	returns a pointer to the extras object (may be null)
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.
void	transfer_from_sequence (ClusterSequence &from_seq)
	transfer the sequence contained in other_seq into our own; any plugin "extras" contained in the from_seq will be lost from there.
Static Public Member Functions
static void	set_jet_algorithm (JetAlgorithm jet_algorithm)
	set the default (static) jet finder across all current and future ClusterSequence objects -- deprecated and obsolescent (i.e.
static void	set_jet_finder (JetAlgorithm jet_algorithm)
	same as above for backward compatibility
Protected Member Functions
bool	_potentially_valid (const PseudoJet &jet) const
	returns true if the jet has a history index contained within the range of this CS
template<class L>
void	_transfer_input_jets (const std::vector< L > &pseudojets)
	transfer the vector<L> of input jets into our own vector<PseudoJet> _jets (with some reserved space for future growth).
void	_initialise_and_run (const JetDefinition &jet_def, const bool &writeout_combinations)
	This is the routine that will do all the initialisation and then run the clustering (may be called by various constructors).
void	_initialise_and_run (const double &R, const Strategy &strategy, const bool &writeout_combinations)
	This is an alternative routine for initialising and running the clustering, provided for legacy purposes.
void	_decant_options (const JetDefinition &jet_def, const bool &writeout_combinations)
	fills in the various member variables with "decanted" options from the jet_definition and writeout_combinations variables
void	_fill_initial_history ()
	fill out the history (and jet cross refs) related to the initial set of jets (assumed already to have been "transferred"), without any clustering
void	_do_ij_recombination_step (const int &jet_i, const int &jet_j, const double &dij, int &newjet_k)
	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
JetAlgorithm	_jet_algorithm
Static Protected Attributes
static JetAlgorithm	_default_jet_algorithm = kt_algorithm
Private Types
typedef std::pair< int, int >	TwoVertices
typedef std::pair< double, TwoVertices >	DijEntry
typedef std::multimap< double, TwoVertices >	DistMap
Private Member Functions
void	_really_dumb_cluster ()
	Run the clustering in a very slow variant of the N^3 algorithm.
void	_delaunay_cluster ()
	Run the clustering using a Hierarchical Delaunay triangulation and STL maps to achieve O(N*ln N) behaviour.
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 + 20.3 and then carries on with 2pi+2R
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
void	_add_untagged_neighbours_to_tile_union (const int tile_index, std::vector< int > &tile_union, int &n_near_tiles)
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...