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