FastJet 3.0.0
|
00001 //STARTHEADER 00002 // $Id: ClusterSequence.hh 2596 2011-09-23 10:09:17Z salam $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 00030 #ifndef __FASTJET_CLUSTERSEQUENCE_HH__ 00031 #define __FASTJET_CLUSTERSEQUENCE_HH__ 00032 00033 #include<vector> 00034 #include<map> 00035 #include "fastjet/internal/DynamicNearestNeighbours.hh" 00036 #include "fastjet/PseudoJet.hh" 00037 #include<memory> 00038 #include<cassert> 00039 #include<iostream> 00040 #include<string> 00041 #include<set> 00042 #include<cmath> // needed to get double std::abs(double) 00043 #include "fastjet/Error.hh" 00044 #include "fastjet/JetDefinition.hh" 00045 #include "fastjet/SharedPtr.hh" 00046 #include "fastjet/LimitedWarning.hh" 00047 #include "fastjet/FunctionOfPseudoJet.hh" 00048 #include "fastjet/ClusterSequenceStructure.hh" 00049 00050 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00051 00052 //---------------------------------------------------------------------- 00053 // here's where we put the main page for fastjet (as explained in the 00054 // Doxygen faq) 00055 // We put in inside te fastjet namespace to have the links without 00056 // having to specify (fastjet::) 00057 //...................................................................... 00058 /** \mainpage FastJet code documentation 00059 * 00060 * These pages provide automatically generated documentation for the 00061 * FastJet package. 00062 * 00063 * \section useful_classes The most useful classes 00064 * 00065 * Many of the facilities of FastJet can be accessed through the three 00066 * following classes: 00067 * 00068 * - PseudoJet: the basic class for holding the 4-momentum of a 00069 * particle or a jet. 00070 * 00071 * - JetDefinition: the combination of a #JetAlgorithm and its 00072 * associated parameters. 00073 * 00074 * - ClusterSequence: constructed with a vector of input (PseudoJet) 00075 * particles and a JetDefinition, it computes and stores the 00076 * information on how the input particles are clustered into jets. 00077 * 00078 * \section advanced_classes Selected more advanced classes 00079 * 00080 * - ClusterSequenceArea: with the help of an AreaDefinition, provides 00081 * jets that also contain information about their area. 00082 * 00083 * \section Tools Selected additional tools 00084 * 00085 * - JetMedianBackgroundEstimator: with the help of a Selector, a JetDefinition and 00086 * an AreaDefinition, allows one to estimate the background noise density in an event; for a simpler, quicker, effective alternative, use GridMedianBackgroundEstimator 00087 * 00088 * - Transformer: class from which are derived various tools for 00089 * manipulating jets and accessing their substructure. Examples are 00090 * Subtractor, Filter, Pruner and various taggers (e.g. JHTopTagger 00091 * and MassDropTagger). 00092 * 00093 * \section further_info Further information 00094 * 00095 * - Selected classes ordered by topics can be found under the <a 00096 * href="modules.html">modules</a> tab. 00097 * 00098 * - The complete list of classes is available under the <a 00099 * href="annotated.html">classes</a> tab. 00100 * 00101 * - For non-class material (<a href="namespacefastjet.html#enum-members">enums</a>, 00102 * <a href="namespacefastjet.html#typedef-members">typedefs</a>, 00103 * <a href="namespacefastjet.html#func-members">functions</a>), see the 00104 * #fastjet documentation 00105 * 00106 * - For further information and normal documentation, see the main <a 00107 * href="http://www.lpthe.jussieu.fr/~salam/fastjet">FastJet</a> page. 00108 * 00109 * \section examples Examples 00110 * See our \subpage Examples page 00111 */ 00112 //---------------------------------------------------------------------- 00113 00114 // forward declaration 00115 class ClusterSequenceStructure; 00116 00117 /// @ingroup basic_classes 00118 /// \class ClusterSequence 00119 /// deals with clustering 00120 class ClusterSequence { 00121 00122 00123 public: 00124 00125 /// default constructor 00126 ClusterSequence () : _deletes_self_when_unused(false) {} 00127 00128 // /// create a clustersequence starting from the supplied set 00129 // /// of pseudojets and clustering them with the long-invariant 00130 // /// kt algorithm (E-scheme recombination) with the supplied 00131 // /// value for R. 00132 // /// 00133 // /// If strategy=DumbN3 a very stupid N^3 algorithm is used for the 00134 // /// clustering; otherwise strategy = NlnN* uses cylinders algorithms 00135 // /// with some number of pi coverage. If writeout_combinations=true a 00136 // /// summary of the recombination sequence is written out 00137 // template<class L> ClusterSequence (const std::vector<L> & pseudojets, 00138 // const double & R = 1.0, 00139 // const Strategy & strategy = Best, 00140 // const bool & writeout_combinations = false); 00141 00142 00143 /// create a clustersequence starting from the supplied set 00144 /// of pseudojets and clustering them with jet definition specified 00145 /// by jet_def (which also specifies the clustering strategy) 00146 template<class L> ClusterSequence ( 00147 const std::vector<L> & pseudojets, 00148 const JetDefinition & jet_def, 00149 const bool & writeout_combinations = false); 00150 00151 /// copy constructor for a ClusterSequence 00152 ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) { 00153 transfer_from_sequence(cs); 00154 } 00155 00156 // virtual ClusterSequence destructor, in case any derived class 00157 // thinks of needing a destructor at some point 00158 virtual ~ClusterSequence (); //{} 00159 00160 // NB: in the routines that follow, for extracting lists of jets, a 00161 // list structure might be more efficient, if sometimes a little 00162 // more awkward to use (at least for old fortran hands). 00163 00164 /// return a vector of all jets (in the sense of the inclusive 00165 /// algorithm) with pt >= ptmin. Time taken should be of the order 00166 /// of the number of jets returned. 00167 std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const; 00168 00169 /// return the number of jets (in the sense of the exclusive 00170 /// algorithm) that would be obtained when running the algorithm 00171 /// with the given dcut. 00172 int n_exclusive_jets (const double & dcut) const; 00173 00174 /// return a vector of all jets (in the sense of the exclusive 00175 /// algorithm) that would be obtained when running the algorithm 00176 /// with the given dcut. 00177 std::vector<PseudoJet> exclusive_jets (const double & dcut) const; 00178 00179 /// return a vector of all jets when the event is clustered (in the 00180 /// exclusive sense) to exactly njets. 00181 /// 00182 /// If there are fewer than njets particles in the ClusterSequence 00183 /// an error is thrown 00184 std::vector<PseudoJet> exclusive_jets (const int & njets) const; 00185 00186 /// return a vector of all jets when the event is clustered (in the 00187 /// exclusive sense) to exactly njets. 00188 /// 00189 /// If there are fewer than njets particles in the ClusterSequence 00190 /// the function just returns however many particles there were. 00191 std::vector<PseudoJet> exclusive_jets_up_to (const int & njets) const; 00192 00193 /// return the dmin corresponding to the recombination that went 00194 /// from n+1 to n jets (sometimes known as d_{n n+1}). If the number 00195 /// of particles in the event is <= njets, the function returns 0. 00196 double exclusive_dmerge (const int & njets) const; 00197 00198 /// return the maximum of the dmin encountered during all recombinations 00199 /// up to the one that led to an n-jet final state; identical to 00200 /// exclusive_dmerge, except in cases where the dmin do not increase 00201 /// monotonically. 00202 double exclusive_dmerge_max (const int & njets) const; 00203 00204 /// return the ymin corresponding to the recombination that went from 00205 /// n+1 to n jets (sometimes known as y_{n n+1}). 00206 double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();} 00207 00208 /// same as exclusive_dmerge_max, but normalised to squared total energy 00209 double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();} 00210 00211 /// the number of exclusive jets at the given ycut 00212 int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());} 00213 00214 /// the exclusive jets obtained at the given ycut 00215 std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const { 00216 int njets = n_exclusive_jets_ycut(ycut); 00217 return exclusive_jets(njets); 00218 } 00219 00220 00221 //int n_exclusive_jets (const PseudoJet & jet, const double & dcut) const; 00222 00223 /// return a vector of all subjets of the current jet (in the sense 00224 /// of the exclusive algorithm) that would be obtained when running 00225 /// the algorithm with the given dcut. 00226 /// 00227 /// Time taken is O(m ln m), where m is the number of subjets that 00228 /// are found. If m gets to be of order of the total number of 00229 /// constituents in the jet, this could be substantially slower than 00230 /// just getting that list of constituents. 00231 std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 00232 const double & dcut) const; 00233 00234 /// return the size of exclusive_subjets(...); still n ln n with same 00235 /// coefficient, but marginally more efficient than manually taking 00236 /// exclusive_subjets.size() 00237 int n_exclusive_subjets(const PseudoJet & jet, 00238 const double & dcut) const; 00239 00240 /// return the list of subjets obtained by unclustering the supplied 00241 /// jet down to nsub subjets. Throws an error if there are fewer than 00242 /// nsub particles in the jet. 00243 /// 00244 /// This requires nsub ln nsub time 00245 std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet, 00246 int nsub) const; 00247 00248 /// return the list of subjets obtained by unclustering the supplied 00249 /// jet down to nsub subjets (or all constituents if there are fewer 00250 /// than nsub). 00251 /// 00252 /// This requires nsub ln nsub time 00253 std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet, 00254 int nsub) const; 00255 00256 /// return the dij that was present in the merging nsub+1 -> nsub 00257 /// subjets inside this jet. 00258 /// 00259 /// Returns 0 if there were nsub or fewer constituents in the jet. 00260 double exclusive_subdmerge(const PseudoJet & jet, int nsub) const; 00261 00262 /// return the maximum dij that occurred in the whole event at the 00263 /// stage that the nsub+1 -> nsub merge of subjets occurred inside 00264 /// this jet. 00265 /// 00266 /// Returns 0 if there were nsub or fewer constituents in the jet. 00267 double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const; 00268 00269 //std::vector<PseudoJet> exclusive_jets (const PseudoJet & jet, 00270 // const int & njets) const; 00271 //double exclusive_dmerge (const PseudoJet & jet, const int & njets) const; 00272 00273 /// returns the sum of all energies in the event (relevant mainly for e+e-) 00274 double Q() const {return _Qtot;} 00275 /// return Q()^2 00276 double Q2() const {return _Qtot*_Qtot;} 00277 00278 /// returns true iff the object is included in the jet. 00279 /// 00280 /// NB: this is only sensible if the object is already registered 00281 /// within the cluster sequence, so you cannot use it with an input 00282 /// particle to the CS (since the particle won't have the history 00283 /// index set properly). 00284 /// 00285 /// For nice clustering structures it should run in O(ln(N)) time 00286 /// but in worst cases (certain cone plugins) it can take O(n) time, 00287 /// where n is the number of particles in the jet. 00288 bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const; 00289 00290 /// if the jet has parents in the clustering, it returns true 00291 /// and sets parent1 and parent2 equal to them. 00292 /// 00293 /// if it has no parents it returns false and sets parent1 and 00294 /// parent2 to zero 00295 bool has_parents(const PseudoJet & jet, PseudoJet & parent1, 00296 PseudoJet & parent2) const; 00297 00298 /// if the jet has a child then return true and give the child jet 00299 /// otherwise return false and set the child to zero 00300 bool has_child(const PseudoJet & jet, PseudoJet & child) const; 00301 00302 /// Version of has_child that sets a pointer to the child if the child 00303 /// exists; 00304 bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const; 00305 00306 /// if this jet has a child (and so a partner) return true 00307 /// and give the partner, otherwise return false and set the 00308 /// partner to zero 00309 bool has_partner(const PseudoJet & jet, PseudoJet & partner) const; 00310 00311 00312 /// return a vector of the particles that make up jet 00313 std::vector<PseudoJet> constituents (const PseudoJet & jet) const; 00314 00315 00316 /// output the supplied vector of jets in a format that can be read 00317 /// by an appropriate root script; the format is: 00318 /// jet-n jet-px jet-py jet-pz jet-E 00319 /// particle-n particle-rap particle-phi particle-pt 00320 /// particle-n particle-rap particle-phi particle-pt 00321 /// ... 00322 /// #END 00323 /// ... [i.e. above repeated] 00324 void print_jets_for_root(const std::vector<PseudoJet> & jets, 00325 std::ostream & ostr = std::cout) const; 00326 00327 /// print jets for root to the file labelled filename, with an 00328 /// optional comment at the beginning 00329 void print_jets_for_root(const std::vector<PseudoJet> & jets, 00330 const std::string & filename, 00331 const std::string & comment = "") const; 00332 00333 // Not yet. Perhaps in a future release. 00334 // /// print out all inclusive jets with pt > ptmin 00335 // virtual void print_jets (const double & ptmin=0.0) const; 00336 00337 /// add on to subjet_vector the constituents of jet (for internal use mainly) 00338 void add_constituents (const PseudoJet & jet, 00339 std::vector<PseudoJet> & subjet_vector) const; 00340 00341 /// return the enum value of the strategy used to cluster the event 00342 inline Strategy strategy_used () const {return _strategy;} 00343 00344 /// return the name of the strategy used to cluster the event 00345 std::string strategy_string () const {return strategy_string(_strategy);} 00346 00347 /// return the name of the strategy associated with the enum strategy_in 00348 std::string strategy_string (Strategy strategy_in) const; 00349 00350 00351 /// return a reference to the jet definition 00352 const JetDefinition & jet_def() const {return _jet_def;} 00353 00354 /// by calling this routine you tell the ClusterSequence to delete 00355 /// itself when all the Pseudojets associated with it have gone out 00356 /// of scope. 00357 /// 00358 /// At the time you call this, there must be at least one jet or 00359 /// other object outside the CS that is associated with the CS 00360 /// (e.g. the result of inclusive_jets()). 00361 /// 00362 /// NB: after having made this call, the user is still allowed to 00363 /// delete the CS or let it go out of scope. Jets associated with it 00364 /// will then simply not be able to access their substructure after 00365 /// that point. 00366 void delete_self_when_unused(); 00367 00368 /// return true if the object has been told to delete itself 00369 /// when unused 00370 bool will_delete_self_when_unused() const {return _deletes_self_when_unused;} 00371 00372 /// tell the ClusterSequence it's about to be self deleted (internal use only) 00373 void signal_imminent_self_deletion() const; 00374 00375 /// returns the scale associated with a jet as required for this 00376 /// clustering algorithm (kt^2 for the kt-algorithm, 1 for the 00377 /// Cambridge algorithm). [May become virtual at some point] 00378 double jet_scale_for_algorithm(const PseudoJet & jet) const; 00379 00380 /// 00381 00382 //----- next follow functions designed specifically for plugins, which 00383 // may only be called when plugin_activated() returns true 00384 00385 /// record the fact that there has been a recombination between 00386 /// jets()[jet_i] and jets()[jet_k], with the specified dij, and 00387 /// return the index (newjet_k) allocated to the new jet, whose 00388 /// momentum is assumed to be the 4-vector sum of that of jet_i and 00389 /// jet_j 00390 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 00391 int & newjet_k) { 00392 assert(plugin_activated()); 00393 _do_ij_recombination_step(jet_i, jet_j, dij, newjet_k); 00394 } 00395 00396 /// as for the simpler variant of plugin_record_ij_recombination, 00397 /// except that the new jet is attributed the momentum and 00398 /// user_index of newjet 00399 void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, 00400 const PseudoJet & newjet, 00401 int & newjet_k); 00402 00403 /// record the fact that there has been a recombination between 00404 /// jets()[jet_i] and the beam, with the specified diB; when looking 00405 /// for inclusive jets, any iB recombination will returned to the user 00406 /// as a jet. 00407 void plugin_record_iB_recombination(int jet_i, double diB) { 00408 assert(plugin_activated()); 00409 _do_iB_recombination_step(jet_i, diB); 00410 } 00411 00412 /// @ingroup extra_info 00413 /// \class Extras 00414 /// base class to store extra information that plugins may provide 00415 /// 00416 /// a class intended to serve as a base in case a plugin needs to 00417 /// associate extra information with a ClusterSequence (see 00418 /// SISConePlugin.* for an example). 00419 class Extras { 00420 public: 00421 virtual ~Extras() {} 00422 virtual std::string description() const {return "This is a dummy extras class that contains no extra information! Derive from it if you want to use it to provide extra information from a plugin jet finder";} 00423 }; 00424 00425 /// the plugin can associate some extra information with the 00426 /// ClusterSequence object by calling this function 00427 inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) { 00428 //_extras = extras_in; 00429 _extras.reset(extras_in.release()); 00430 } 00431 00432 /// returns true when the plugin is allowed to run the show. 00433 inline bool plugin_activated() const {return _plugin_activated;} 00434 00435 /// returns a pointer to the extras object (may be null) 00436 const Extras * extras() const {return _extras.get();} 00437 00438 /// allows a plugin to run a templated clustering (nearest-neighbour heuristic) 00439 /// 00440 /// This has N^2 behaviour on "good" distance, but a worst case behaviour 00441 /// of N^3 (and many algs trigger the worst case behaviour) 00442 /// 00443 /// 00444 /// For more details on how this works, see GenBriefJet below 00445 template<class GBJ> void plugin_simple_N2_cluster () { 00446 assert(plugin_activated()); 00447 _simple_N2_cluster<GBJ>(); 00448 } 00449 00450 00451 public: 00452 //DEP /// set the default (static) jet finder across all current and future 00453 //DEP /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be 00454 //DEP /// suppressed in a future release). 00455 //DEP static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;} 00456 //DEP /// same as above for backward compatibility 00457 //DEP static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;} 00458 00459 00460 /// \ingroup extra_info 00461 /// \struct history_element 00462 /// a single element in the clustering history 00463 /// 00464 /// (see vector _history below). 00465 struct history_element{ 00466 int parent1; /// index in _history where first parent of this jet 00467 /// was created (InexistentParent if this jet is an 00468 /// original particle) 00469 00470 int parent2; /// index in _history where second parent of this jet 00471 /// was created (InexistentParent if this jet is an 00472 /// original particle); BeamJet if this history entry 00473 /// just labels the fact that the jet has recombined 00474 /// with the beam) 00475 00476 int child; /// index in _history where the current jet is 00477 /// recombined with another jet to form its child. It 00478 /// is Invalid if this jet does not further 00479 /// recombine. 00480 00481 int jetp_index; /// index in the _jets vector where we will find the 00482 /// PseudoJet object corresponding to this jet 00483 /// (i.e. the jet created at this entry of the 00484 /// history). NB: if this element of the history 00485 /// corresponds to a beam recombination, then 00486 /// jetp_index=Invalid. 00487 00488 double dij; /// the distance corresponding to the recombination 00489 /// at this stage of the clustering. 00490 00491 double max_dij_so_far; /// the largest recombination distance seen 00492 /// so far in the clustering history. 00493 }; 00494 00495 enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1}; 00496 00497 /// allow the user to access the internally stored _jets() array, 00498 /// which contains both the initial particles and the various 00499 /// intermediate and final stages of recombination. 00500 /// 00501 /// The first n_particles() entries are the original particles, 00502 /// in the order in which they were supplied to the ClusterSequence 00503 /// constructor. It can be useful to access them for example when 00504 /// examining whether a given input object is part of a specific 00505 /// jet, via the objects_in_jet(...) member function (which only takes 00506 /// PseudoJets that are registered in the ClusterSequence). 00507 /// 00508 /// One of the other (internal uses) is related to the fact 00509 /// because we don't seem to be able to access protected elements of 00510 /// the class for an object that is not "this" (at least in case where 00511 /// "this" is of a slightly different kind from the object, both 00512 /// derived from ClusterSequence). 00513 const std::vector<PseudoJet> & jets() const; 00514 00515 /// allow the user to access the raw internal history. 00516 /// 00517 /// This is present (as for jets()) in part so that protected 00518 /// derived classes can access this information about other 00519 /// ClusterSequences. 00520 /// 00521 /// A user who wishes to follow the details of the ClusterSequence 00522 /// can also make use of this information (and should consult the 00523 /// history_element documentation for more information), but should 00524 /// be aware that these internal structures may evolve in future 00525 /// FastJet versions. 00526 const std::vector<history_element> & history() const; 00527 00528 /// returns the number of particles that were provided to the 00529 /// clustering algorithm (helps the user find their way around the 00530 /// history and jets objects if they weren't paying attention 00531 /// beforehand). 00532 unsigned int n_particles() const; 00533 00534 /// returns a vector of size n_particles() which indicates, for 00535 /// each of the initial particles (in the order in which they were 00536 /// supplied), which of the supplied jets it belongs to; if it does 00537 /// not belong to any of the supplied jets, the index is set to -1; 00538 std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const; 00539 00540 /// routine that returns an order in which to read the history 00541 /// such that clusterings that lead to identical jet compositions 00542 /// but different histories (because of degeneracies in the 00543 /// clustering order) will have matching constituents for each 00544 /// matching entry in the unique_history_order. 00545 /// 00546 /// The order has the property that an entry's parents will always 00547 /// appear prior to that entry itself. 00548 /// 00549 /// Roughly speaking the order is such that we first provide all 00550 /// steps that lead to the final jet containing particle 1; then we 00551 /// have the steps that lead to reconstruction of the jet containing 00552 /// the next-lowest-numbered unclustered particle, etc... 00553 /// [see GPS CCN28-12 for more info -- of course a full explanation 00554 /// here would be better...] 00555 std::vector<int> unique_history_order() const; 00556 00557 /// return the set of particles that have not been clustered. For 00558 /// kt and cam/aachen algorithms this should always be null, but for 00559 /// cone type algorithms it can be non-null; 00560 std::vector<PseudoJet> unclustered_particles() const; 00561 00562 /// Return the list of pseudojets in the ClusterSequence that do not 00563 /// have children (and are not among the inclusive jets). They may 00564 /// result from a clustering step or may be one of the pseudojets 00565 /// returned by unclustered_particles(). 00566 std::vector<PseudoJet> childless_pseudojets() const; 00567 00568 /// returns true if the object (jet or particle) is contained by (ie 00569 /// belongs to) this cluster sequence. 00570 /// 00571 /// Tests performed: if thejet's interface is this cluster sequence 00572 /// and its cluster history index is in a consistent range. 00573 bool contains(const PseudoJet & object) const; 00574 00575 /// transfer the sequence contained in other_seq into our own; 00576 /// any plugin "extras" contained in the from_seq will be lost 00577 /// from there. 00578 /// 00579 /// It also sets the ClusterSequence pointers of the PseudoJets in 00580 /// the history to point to this ClusterSequence 00581 /// 00582 /// When specified, the second argument is an action that will be 00583 /// applied on every jets in the resulting ClusterSequence 00584 void transfer_from_sequence(const ClusterSequence & from_seq, 00585 const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0); 00586 00587 /// retrieve a shared pointer to the wrapper to this ClusterSequence 00588 /// 00589 /// this may turn useful if you want to track when this 00590 /// ClusterSequence goes out of scope 00591 const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const{ 00592 return _structure_shared_ptr; 00593 } 00594 00595 /// the structure type associated with a jet belonging to a ClusterSequence 00596 typedef ClusterSequenceStructure StructureType; 00597 00598 00599 protected: 00600 //DEP static JetAlgorithm _default_jet_algorithm; 00601 JetDefinition _jet_def; 00602 00603 /// transfer the vector<L> of input jets into our own vector<PseudoJet> 00604 /// _jets (with some reserved space for future growth). 00605 template<class L> void _transfer_input_jets( 00606 const std::vector<L> & pseudojets); 00607 00608 /// This is the routine that will do all the initialisation and 00609 /// then run the clustering (may be called by various constructors). 00610 /// It assumes _jets contains the momenta to be clustered. 00611 void _initialise_and_run (const JetDefinition & jet_def, 00612 const bool & writeout_combinations); 00613 00614 //DEP /// This is an alternative routine for initialising and running the 00615 //DEP /// clustering, provided for legacy purposes. The jet finder is that 00616 //DEP /// specified in the static member _default_jet_algorithm. 00617 //DEP void _initialise_and_run (const double & R, 00618 //DEP const Strategy & strategy, 00619 //DEP const bool & writeout_combinations); 00620 00621 /// fills in the various member variables with "decanted" options from 00622 /// the jet_definition and writeout_combinations variables 00623 void _decant_options(const JetDefinition & jet_def, 00624 const bool & writeout_combinations); 00625 00626 /// fill out the history (and jet cross refs) related to the initial 00627 /// set of jets (assumed already to have been "transferred"), 00628 /// without any clustering 00629 void _fill_initial_history(); 00630 00631 /// carry out the recombination between the jets numbered jet_i and 00632 /// jet_j, at distance scale dij; return the index newjet_k of the 00633 /// result of the recombination of i and j. 00634 void _do_ij_recombination_step(const int & jet_i, const int & jet_j, 00635 const double & dij, int & newjet_k); 00636 00637 /// carry out an recombination step in which _jets[jet_i] merges with 00638 /// the beam, 00639 void _do_iB_recombination_step(const int & jet_i, const double & diB); 00640 00641 /// every time a jet is added internally during clustering, this 00642 /// should be called to set the jet's structure shared ptr to point 00643 /// to the CS (and the count of internally associated objects is 00644 /// also updated). This should not be called outside construction of 00645 /// a CS object. 00646 void _set_structure_shared_ptr(PseudoJet & j); 00647 00648 /// make sure that the CS's internal tally of the use count matches 00649 /// that of the _structure_shared_ptr 00650 void _update_structure_use_count(); 00651 00652 00653 /// This contains the physical PseudoJets; for each PseudoJet one 00654 /// can find the corresponding position in the _history by looking 00655 /// at _jets[i].cluster_hist_index(). 00656 std::vector<PseudoJet> _jets; 00657 00658 00659 /// this vector will contain the branching history; for each stage, 00660 /// _history[i].jetp_index indicates where to look in the _jets 00661 /// vector to get the physical PseudoJet. 00662 std::vector<history_element> _history; 00663 00664 /// set subhist to be a set pointers to history entries corresponding to the 00665 /// subjets of this jet; one stops going working down through the 00666 /// subjets either when 00667 /// - there is no further to go 00668 /// - one has found maxjet entries 00669 /// - max_dij_so_far <= dcut 00670 /// By setting maxjet=0 one can use just dcut; by setting dcut<0 00671 /// one can use jet maxjet 00672 void get_subhist_set(std::set<const history_element*> & subhist, 00673 const PseudoJet & jet, double dcut, int maxjet) const; 00674 00675 bool _writeout_combinations; 00676 int _initial_n; 00677 double _Rparam, _R2, _invR2; 00678 double _Qtot; 00679 Strategy _strategy; 00680 JetAlgorithm _jet_algorithm; 00681 00682 SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure 00683 int _structure_use_count_after_construction; //< info of use when CS handles its own memory 00684 /// if true then the CS will delete itself when the last external 00685 /// object referring to it disappears. It is mutable so as to ensure 00686 /// that signal_imminent_self_deletion() [const] can make relevant 00687 /// changes. 00688 mutable bool _deletes_self_when_unused; 00689 00690 private: 00691 00692 bool _plugin_activated; 00693 //std::auto_ptr<Extras> _extras; // things the plugin might want to add 00694 SharedPtr<Extras> _extras; // things the plugin might want to add 00695 00696 void _really_dumb_cluster (); 00697 void _delaunay_cluster (); 00698 //void _simple_N2_cluster (); 00699 template<class BJ> void _simple_N2_cluster (); 00700 void _tiled_N2_cluster (); 00701 void _faster_tiled_N2_cluster (); 00702 00703 // 00704 void _minheap_faster_tiled_N2_cluster(); 00705 00706 // things needed specifically for Cambridge with Chan's 2D closest 00707 // pairs method 00708 void _CP2DChan_cluster(); 00709 void _CP2DChan_cluster_2pi2R (); 00710 void _CP2DChan_cluster_2piMultD (); 00711 void _CP2DChan_limited_cluster(double D); 00712 void _do_Cambridge_inclusive_jets(); 00713 00714 // NSqrtN method for C/A 00715 void _fast_NsqrtN_cluster(); 00716 00717 void _add_step_to_history(const int & step_number, const int & parent1, 00718 const int & parent2, const int & jetp_index, 00719 const double & dij); 00720 00721 /// internal routine associated with the construction of the unique 00722 /// history order (following children in the tree) 00723 void _extract_tree_children(int pos, std::valarray<bool> &, 00724 const std::valarray<int> &, std::vector<int> &) const; 00725 00726 /// internal routine associated with the construction of the unique 00727 /// history order (following parents in the tree) 00728 void _extract_tree_parents (int pos, std::valarray<bool> &, 00729 const std::valarray<int> &, std::vector<int> &) const; 00730 00731 00732 // these will be useful shorthands in the Voronoi-based code 00733 typedef std::pair<int,int> TwoVertices; 00734 typedef std::pair<double,TwoVertices> DijEntry; 00735 typedef std::multimap<double,TwoVertices> DistMap; 00736 00737 /// currently used only in the Voronoi based code 00738 void _add_ktdistance_to_map(const int & ii, 00739 DistMap & DijMap, 00740 const DynamicNearestNeighbours * DNN); 00741 00742 /// for making sure the user knows what it is they're running... 00743 void _print_banner(); 00744 00745 00746 /// will be set by default to be true for the first run 00747 static bool _first_time; 00748 00749 /// record the number of warnings provided about the exclusive 00750 /// algorithm -- so that we don't print it out more than a few 00751 /// times. 00752 static int _n_exclusive_warnings; 00753 00754 /// the limited warning member for notification of user that 00755 /// their requested strategy has been overridden (usually because 00756 /// they have R>2pi and not all strategies work then) 00757 static LimitedWarning _changed_strategy_warning; 00758 00759 //---------------------------------------------------------------------- 00760 /// the fundamental structure which contains the minimal info about 00761 /// a jet, as needed for our plain N^2 algorithm -- the idea is to 00762 /// put all info that will be accessed N^2 times into an array of 00763 /// BriefJets... 00764 struct BriefJet { 00765 double eta, phi, kt2, NN_dist; 00766 BriefJet * NN; 00767 int _jets_index; 00768 }; 00769 00770 00771 /// structure analogous to BriefJet, but with the extra information 00772 /// needed for dealing with tiles 00773 class TiledJet { 00774 public: 00775 double eta, phi, kt2, NN_dist; 00776 TiledJet * NN, *previous, * next; 00777 int _jets_index, tile_index, diJ_posn; 00778 // routines that are useful in the minheap version of tiled 00779 // clustering ("misuse" the otherwise unused diJ_posn, so as 00780 // to indicate whether jets need to have their minheap entries 00781 // updated). 00782 inline void label_minheap_update_needed() {diJ_posn = 1;} 00783 inline void label_minheap_update_done() {diJ_posn = 0;} 00784 inline bool minheap_update_needed() const {return diJ_posn==1;} 00785 }; 00786 00787 //-- some of the functions that follow are templates and will work 00788 //as well for briefjet and tiled jets 00789 00790 /// set the kinematic and labelling info for jeta so that it corresponds 00791 /// to _jets[_jets_index] 00792 template <class J> void _bj_set_jetinfo( J * const jet, 00793 const int _jets_index) const; 00794 00795 /// "remove" this jet, which implies updating links of neighbours and 00796 /// perhaps modifying the tile structure 00797 void _bj_remove_from_tiles( TiledJet * const jet) const; 00798 00799 /// return the distance between two BriefJet objects 00800 template <class J> double _bj_dist(const J * const jeta, 00801 const J * const jetb) const; 00802 00803 // return the diJ (multiplied by _R2) for this jet assuming its NN 00804 // info is correct 00805 template <class J> double _bj_diJ(const J * const jeta) const; 00806 00807 /// for testing purposes only: if in the range head--tail-1 there is a 00808 /// a jet which corresponds to hist_index in the history, then 00809 /// return a pointer to that jet; otherwise return tail. 00810 template <class J> inline J * _bj_of_hindex( 00811 const int hist_index, 00812 J * const head, J * const tail) 00813 const { 00814 J * res; 00815 for(res = head; res<tail; res++) { 00816 if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;} 00817 } 00818 return res; 00819 } 00820 00821 00822 //-- remaining functions are different in various cases, so we 00823 // will use templates but are not sure if they're useful... 00824 00825 /// updates (only towards smaller distances) the NN for jeta without checking 00826 /// whether in the process jeta itself might be a new NN of one of 00827 /// the jets being scanned -- span the range head to tail-1 with 00828 /// assumption that jeta is not contained in that range 00829 template <class J> void _bj_set_NN_nocross(J * const jeta, 00830 J * const head, const J * const tail) const; 00831 00832 /// reset the NN for jeta and DO check whether in the process jeta 00833 /// itself might be a new NN of one of the jets being scanned -- 00834 /// span the range head to tail-1 with assumption that jeta is not 00835 /// contained in that range 00836 template <class J> void _bj_set_NN_crosscheck(J * const jeta, 00837 J * const head, const J * const tail) const; 00838 00839 00840 00841 /// number of neighbours that a tile will have (rectangular geometry 00842 /// gives 9 neighbours). 00843 static const int n_tile_neighbours = 9; 00844 //---------------------------------------------------------------------- 00845 /// The fundamental structures to be used for the tiled N^2 algorithm 00846 /// (see CCN27-44 for some discussion of pattern of tiling) 00847 struct Tile { 00848 /// pointers to neighbouring tiles, including self 00849 Tile * begin_tiles[n_tile_neighbours]; 00850 /// neighbouring tiles, excluding self 00851 Tile ** surrounding_tiles; 00852 /// half of neighbouring tiles, no self 00853 Tile ** RH_tiles; 00854 /// just beyond end of tiles 00855 Tile ** end_tiles; 00856 /// start of list of BriefJets contained in this tile 00857 TiledJet * head; 00858 /// sometimes useful to be able to tag a tile 00859 bool tagged; 00860 }; 00861 std::vector<Tile> _tiles; 00862 double _tiles_eta_min, _tiles_eta_max; 00863 double _tile_size_eta, _tile_size_phi; 00864 int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max; 00865 00866 // reasonably robust return of tile index given ieta and iphi, in particular 00867 // it works even if iphi is negative 00868 inline int _tile_index (int ieta, int iphi) const { 00869 // note that (-1)%n = -1 so that we have to add _n_tiles_phi 00870 // before performing modulo operation 00871 return (ieta-_tiles_ieta_min)*_n_tiles_phi 00872 + (iphi+_n_tiles_phi) % _n_tiles_phi; 00873 } 00874 00875 // routines for tiled case, including some overloads of the plain 00876 // BriefJet cases 00877 int _tile_index(const double & eta, const double & phi) const; 00878 void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index); 00879 void _bj_remove_from_tiles(TiledJet * const jet); 00880 void _initialise_tiles(); 00881 void _print_tiles(TiledJet * briefjets ) const; 00882 void _add_neighbours_to_tile_union(const int tile_index, 00883 std::vector<int> & tile_union, int & n_near_tiles) const; 00884 void _add_untagged_neighbours_to_tile_union(const int tile_index, 00885 std::vector<int> & tile_union, int & n_near_tiles); 00886 00887 00888 //---------------------------------------------------------------------- 00889 /// fundamental structure for e+e- clustering 00890 struct EEBriefJet { 00891 double NN_dist; // obligatorily present 00892 double kt2; // obligatorily present == E^2 in general 00893 EEBriefJet * NN; // must be present too 00894 int _jets_index; // must also be present! 00895 //........................................................... 00896 double nx, ny, nz; // our internal storage for fast distance calcs 00897 }; 00898 00899 /// to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?) 00900 //void _dummy_N2_cluster_instantiation(); 00901 00902 00903 /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3) 00904 void _simple_N2_cluster_BriefJet(); 00905 /// to avoid issues with template instantiation (OS X 10.3, gcc 3.3) 00906 void _simple_N2_cluster_EEBriefJet(); 00907 }; 00908 00909 00910 //********************************************************************** 00911 //************** START OF INLINE MATERIAL ****************** 00912 //********************************************************************** 00913 00914 00915 //---------------------------------------------------------------------- 00916 // Transfer the initial jets into our internal structure 00917 template<class L> void ClusterSequence::_transfer_input_jets( 00918 const std::vector<L> & pseudojets) { 00919 00920 // this will ensure that we can point to jets without difficulties 00921 // arising. 00922 _jets.reserve(pseudojets.size()*2); 00923 00924 // insert initial jets this way so that any type L that can be 00925 // converted to a pseudojet will work fine (basically PseudoJet 00926 // and any type that has [] subscript access to the momentum 00927 // components, such as CLHEP HepLorentzVector). 00928 for (unsigned int i = 0; i < pseudojets.size(); i++) { 00929 _jets.push_back(pseudojets[i]);} 00930 00931 } 00932 00933 // //---------------------------------------------------------------------- 00934 // // initialise from some generic type... Has to be made available 00935 // // here in order for it the template aspect of it to work... 00936 // template<class L> ClusterSequence::ClusterSequence ( 00937 // const std::vector<L> & pseudojets, 00938 // const double & R, 00939 // const Strategy & strategy, 00940 // const bool & writeout_combinations) { 00941 // 00942 // // transfer the initial jets (type L) into our own array 00943 // _transfer_input_jets(pseudojets); 00944 // 00945 // // run the clustering 00946 // _initialise_and_run(R,strategy,writeout_combinations); 00947 // } 00948 00949 00950 //---------------------------------------------------------------------- 00951 /// constructor of a jet-clustering sequence from a vector of 00952 /// four-momenta, with the jet definition specified by jet_def 00953 template<class L> ClusterSequence::ClusterSequence ( 00954 const std::vector<L> & pseudojets, 00955 const JetDefinition & jet_def, 00956 const bool & writeout_combinations) { 00957 00958 // transfer the initial jets (type L) into our own array 00959 _transfer_input_jets(pseudojets); 00960 00961 // run the clustering 00962 _initialise_and_run(jet_def,writeout_combinations); 00963 } 00964 00965 00966 inline const std::vector<PseudoJet> & ClusterSequence::jets () const { 00967 return _jets; 00968 } 00969 00970 inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const { 00971 return _history; 00972 } 00973 00974 inline unsigned int ClusterSequence::n_particles() const {return _initial_n;} 00975 00976 00977 00978 //---------------------------------------------------------------------- 00979 template <class J> inline void ClusterSequence::_bj_set_jetinfo( 00980 J * const jetA, const int _jets_index) const { 00981 jetA->eta = _jets[_jets_index].rap(); 00982 jetA->phi = _jets[_jets_index].phi_02pi(); 00983 jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]); 00984 jetA->_jets_index = _jets_index; 00985 // initialise NN info as well 00986 jetA->NN_dist = _R2; 00987 jetA->NN = NULL; 00988 } 00989 00990 00991 00992 00993 //---------------------------------------------------------------------- 00994 template <class J> inline double ClusterSequence::_bj_dist( 00995 const J * const jetA, const J * const jetB) const { 00996 double dphi = std::abs(jetA->phi - jetB->phi); 00997 double deta = (jetA->eta - jetB->eta); 00998 if (dphi > pi) {dphi = twopi - dphi;} 00999 return dphi*dphi + deta*deta; 01000 } 01001 01002 //---------------------------------------------------------------------- 01003 template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const { 01004 double kt2 = jet->kt2; 01005 if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}} 01006 return jet->NN_dist * kt2; 01007 } 01008 01009 01010 //---------------------------------------------------------------------- 01011 // set the NN for jet without checking whether in the process you might 01012 // have discovered a new nearest neighbour for another jet 01013 template <class J> inline void ClusterSequence::_bj_set_NN_nocross( 01014 J * const jet, J * const head, const J * const tail) const { 01015 double NN_dist = _R2; 01016 J * NN = NULL; 01017 if (head < jet) { 01018 for (J * jetB = head; jetB != jet; jetB++) { 01019 double dist = _bj_dist(jet,jetB); 01020 if (dist < NN_dist) { 01021 NN_dist = dist; 01022 NN = jetB; 01023 } 01024 } 01025 } 01026 if (tail > jet) { 01027 for (J * jetB = jet+1; jetB != tail; jetB++) { 01028 double dist = _bj_dist(jet,jetB); 01029 if (dist < NN_dist) { 01030 NN_dist = dist; 01031 NN = jetB; 01032 } 01033 } 01034 } 01035 jet->NN = NN; 01036 jet->NN_dist = NN_dist; 01037 } 01038 01039 01040 //---------------------------------------------------------------------- 01041 template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet, 01042 J * const head, const J * const tail) const { 01043 double NN_dist = _R2; 01044 J * NN = NULL; 01045 for (J * jetB = head; jetB != tail; jetB++) { 01046 double dist = _bj_dist(jet,jetB); 01047 if (dist < NN_dist) { 01048 NN_dist = dist; 01049 NN = jetB; 01050 } 01051 if (dist < jetB->NN_dist) { 01052 jetB->NN_dist = dist; 01053 jetB->NN = jet; 01054 } 01055 } 01056 jet->NN = NN; 01057 jet->NN_dist = NN_dist; 01058 } 01059 01060 01061 01062 FASTJET_END_NAMESPACE 01063 01064 #endif // __FASTJET_CLUSTERSEQUENCE_HH__