/usr/include/fastjet/ClusterSequence.hh is in libfastjet-dev 3.0.6+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 | //STARTHEADER
// $Id: ClusterSequence.hh 3114 2013-05-04 08:46:00Z salam $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development and are described in hep-ph/0512210. If you use
// FastJet as part of work towards a scientific publication, please
// include a citation to the FastJet paper.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
#ifndef __FASTJET_CLUSTERSEQUENCE_HH__
#define __FASTJET_CLUSTERSEQUENCE_HH__
#include<vector>
#include<map>
#include "fastjet/PseudoJet.hh"
#include<memory>
#include<cassert>
#include<iostream>
#include<string>
#include<set>
#include<cmath> // needed to get double std::abs(double)
#include "fastjet/Error.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/SharedPtr.hh"
#include "fastjet/LimitedWarning.hh"
#include "fastjet/FunctionOfPseudoJet.hh"
#include "fastjet/ClusterSequenceStructure.hh"
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
// forward declarations
class ClusterSequenceStructure;
class DynamicNearestNeighbours;
/// @ingroup basic_classes
/// \class ClusterSequence
/// deals with clustering
class ClusterSequence {
public:
/// default constructor
ClusterSequence () : _deletes_self_when_unused(false) {}
// /// create a clustersequence starting from the supplied set
// /// of pseudojets and clustering them with the long-invariant
// /// kt algorithm (E-scheme recombination) with the supplied
// /// value for R.
// ///
// /// If strategy=DumbN3 a very stupid N^3 algorithm is used for the
// /// clustering; otherwise strategy = NlnN* uses cylinders algorithms
// /// with some number of pi coverage. If writeout_combinations=true a
// /// summary of the recombination sequence is written out
// template<class L> ClusterSequence (const std::vector<L> & pseudojets,
// const double & R = 1.0,
// const Strategy & strategy = Best,
// const bool & writeout_combinations = false);
/// create a clustersequence starting from the supplied set
/// of pseudojets and clustering them with jet definition specified
/// by jet_def (which also specifies the clustering strategy)
template<class L> ClusterSequence (
const std::vector<L> & pseudojets,
const JetDefinition & jet_def,
const bool & writeout_combinations = false);
/// copy constructor for a ClusterSequence
ClusterSequence (const ClusterSequence & cs) : _deletes_self_when_unused(false) {
transfer_from_sequence(cs);
}
// virtual ClusterSequence destructor, in case any derived class
// thinks of needing a destructor at some point
virtual ~ClusterSequence (); //{}
// NB: in the routines that follow, for extracting lists of jets, a
// list structure might be more efficient, if sometimes a little
// more awkward to use (at least for old fortran hands).
/// return a vector of all jets (in the sense of the inclusive
/// algorithm) with pt >= ptmin. Time taken should be of the order
/// of the number of jets returned.
std::vector<PseudoJet> inclusive_jets (const double & ptmin = 0.0) const;
/// return the number of jets (in the sense of the exclusive
/// algorithm) that would be obtained when running the algorithm
/// with the given dcut.
int n_exclusive_jets (const double & dcut) const;
/// return a vector of all jets (in the sense of the exclusive
/// algorithm) that would be obtained when running the algorithm
/// with the given dcut.
std::vector<PseudoJet> exclusive_jets (const double & dcut) const;
/// return a vector of all jets when the event is clustered (in the
/// exclusive sense) to exactly njets.
///
/// If there are fewer than njets particles in the ClusterSequence
/// an error is thrown
std::vector<PseudoJet> exclusive_jets (const int & njets) const;
/// return a vector of all jets when the event is clustered (in the
/// exclusive sense) to exactly njets.
///
/// If there are fewer than njets particles in the ClusterSequence
/// the function just returns however many particles there were.
std::vector<PseudoJet> exclusive_jets_up_to (const int & njets) const;
/// return the dmin corresponding to the recombination that went
/// from n+1 to n jets (sometimes known as d_{n n+1}). If the number
/// of particles in the event is <= njets, the function returns 0.
double exclusive_dmerge (const int & njets) const;
/// return the maximum of the dmin encountered during all recombinations
/// up to the one that led to an n-jet final state; identical to
/// exclusive_dmerge, except in cases where the dmin do not increase
/// monotonically.
double exclusive_dmerge_max (const int & njets) const;
/// return the ymin corresponding to the recombination that went from
/// n+1 to n jets (sometimes known as y_{n n+1}).
double exclusive_ymerge (int njets) const {return exclusive_dmerge(njets) / Q2();}
/// same as exclusive_dmerge_max, but normalised to squared total energy
double exclusive_ymerge_max (int njets) const {return exclusive_dmerge_max(njets)/Q2();}
/// the number of exclusive jets at the given ycut
int n_exclusive_jets_ycut (double ycut) const {return n_exclusive_jets(ycut*Q2());}
/// the exclusive jets obtained at the given ycut
std::vector<PseudoJet> exclusive_jets_ycut (double ycut) const {
int njets = n_exclusive_jets_ycut(ycut);
return exclusive_jets(njets);
}
//int n_exclusive_jets (const PseudoJet & jet, const double & dcut) const;
/// return a vector of all subjets of the current jet (in the sense
/// of the exclusive algorithm) that would be obtained when running
/// the algorithm with the given dcut.
///
/// Time taken is O(m ln m), where m is the number of subjets that
/// are found. If m gets to be of order of the total number of
/// constituents in the jet, this could be substantially slower than
/// just getting that list of constituents.
std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
const double & dcut) const;
/// return the size of exclusive_subjets(...); still n ln n with same
/// coefficient, but marginally more efficient than manually taking
/// exclusive_subjets.size()
int n_exclusive_subjets(const PseudoJet & jet,
const double & dcut) const;
/// return the list of subjets obtained by unclustering the supplied
/// jet down to nsub subjets. Throws an error if there are fewer than
/// nsub particles in the jet.
///
/// This requires nsub ln nsub time
std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
int nsub) const;
/// return the list of subjets obtained by unclustering the supplied
/// jet down to nsub subjets (or all constituents if there are fewer
/// than nsub).
///
/// This requires nsub ln nsub time
std::vector<PseudoJet> exclusive_subjets_up_to (const PseudoJet & jet,
int nsub) const;
/// return the dij that was present in the merging nsub+1 -> nsub
/// subjets inside this jet.
///
/// Returns 0 if there were nsub or fewer constituents in the jet.
double exclusive_subdmerge(const PseudoJet & jet, int nsub) const;
/// return the maximum dij that occurred in the whole event at the
/// stage that the nsub+1 -> nsub merge of subjets occurred inside
/// this jet.
///
/// Returns 0 if there were nsub or fewer constituents in the jet.
double exclusive_subdmerge_max(const PseudoJet & jet, int nsub) const;
//std::vector<PseudoJet> exclusive_jets (const PseudoJet & jet,
// const int & njets) const;
//double exclusive_dmerge (const PseudoJet & jet, const int & njets) const;
/// returns the sum of all energies in the event (relevant mainly for e+e-)
double Q() const {return _Qtot;}
/// return Q()^2
double Q2() const {return _Qtot*_Qtot;}
/// returns true iff the object is included in the jet.
///
/// NB: this is only sensible if the object is already registered
/// within the cluster sequence, so you cannot use it with an input
/// particle to the CS (since the particle won't have the history
/// index set properly).
///
/// For nice clustering structures it should run in O(ln(N)) time
/// but in worst cases (certain cone plugins) it can take O(n) time,
/// where n is the number of particles in the jet.
bool object_in_jet(const PseudoJet & object, const PseudoJet & jet) const;
/// if the jet has parents in the clustering, it returns true
/// and sets parent1 and parent2 equal to them.
///
/// if it has no parents it returns false and sets parent1 and
/// parent2 to zero
bool has_parents(const PseudoJet & jet, PseudoJet & parent1,
PseudoJet & parent2) const;
/// if the jet has a child then return true and give the child jet
/// otherwise return false and set the child to zero
bool has_child(const PseudoJet & jet, PseudoJet & child) const;
/// Version of has_child that sets a pointer to the child if the child
/// exists;
bool has_child(const PseudoJet & jet, const PseudoJet * & childp) const;
/// if this jet has a child (and so a partner) return true
/// and give the partner, otherwise return false and set the
/// partner to zero
bool has_partner(const PseudoJet & jet, PseudoJet & partner) const;
/// return a vector of the particles that make up jet
std::vector<PseudoJet> constituents (const PseudoJet & jet) const;
/// output the supplied vector of jets in a format that can be read
/// by an appropriate root script; the format is:
/// jet-n jet-px jet-py jet-pz jet-E
/// particle-n particle-rap particle-phi particle-pt
/// particle-n particle-rap particle-phi particle-pt
/// ...
/// #END
/// ... [i.e. above repeated]
void print_jets_for_root(const std::vector<PseudoJet> & jets,
std::ostream & ostr = std::cout) const;
/// print jets for root to the file labelled filename, with an
/// optional comment at the beginning
void print_jets_for_root(const std::vector<PseudoJet> & jets,
const std::string & filename,
const std::string & comment = "") const;
// Not yet. Perhaps in a future release.
// /// print out all inclusive jets with pt > ptmin
// virtual void print_jets (const double & ptmin=0.0) const;
/// add on to subjet_vector the constituents of jet (for internal use mainly)
void add_constituents (const PseudoJet & jet,
std::vector<PseudoJet> & subjet_vector) const;
/// return the enum value of the strategy used to cluster the event
inline Strategy strategy_used () const {return _strategy;}
/// return the name of the strategy used to cluster the event
std::string strategy_string () const {return strategy_string(_strategy);}
/// return the name of the strategy associated with the enum strategy_in
std::string strategy_string (Strategy strategy_in) const;
/// return a reference to the jet definition
const JetDefinition & jet_def() const {return _jet_def;}
/// by calling this routine you tell the ClusterSequence to delete
/// itself when all the Pseudojets associated with it have gone out
/// of scope.
///
/// At the time you call this, there must be at least one jet or
/// other object outside the CS that is associated with the CS
/// (e.g. the result of inclusive_jets()).
///
/// NB: after having made this call, the user is still allowed to
/// delete the CS or let it go out of scope. Jets associated with it
/// will then simply not be able to access their substructure after
/// that point.
void delete_self_when_unused();
/// return true if the object has been told to delete itself
/// when unused
bool will_delete_self_when_unused() const {return _deletes_self_when_unused;}
/// tell the ClusterSequence it's about to be self deleted (internal use only)
void signal_imminent_self_deletion() const;
/// returns the scale associated with a jet as required for this
/// clustering algorithm (kt^2 for the kt-algorithm, 1 for the
/// Cambridge algorithm). [May become virtual at some point]
double jet_scale_for_algorithm(const PseudoJet & jet) const;
///
//----- next follow functions designed specifically for plugins, which
// may only be called when plugin_activated() returns true
/// record the fact that there has been a recombination between
/// jets()[jet_i] and jets()[jet_k], with the specified dij, and
/// return the index (newjet_k) allocated to the new jet, whose
/// momentum is assumed to be the 4-vector sum of that of jet_i and
/// jet_j
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
int & newjet_k) {
assert(plugin_activated());
_do_ij_recombination_step(jet_i, jet_j, dij, newjet_k);
}
/// as for the simpler variant of plugin_record_ij_recombination,
/// except that the new jet is attributed the momentum and
/// user_index of newjet
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij,
const PseudoJet & newjet,
int & newjet_k);
/// record the fact that there has been a recombination between
/// jets()[jet_i] and the beam, with the specified diB; when looking
/// for inclusive jets, any iB recombination will returned to the user
/// as a jet.
void plugin_record_iB_recombination(int jet_i, double diB) {
assert(plugin_activated());
_do_iB_recombination_step(jet_i, diB);
}
/// @ingroup extra_info
/// \class Extras
/// base class to store extra information that plugins may provide
///
/// a class intended to serve as a base in case a plugin needs to
/// associate extra information with a ClusterSequence (see
/// SISConePlugin.* for an example).
class Extras {
public:
virtual ~Extras() {}
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";}
};
/// the plugin can associate some extra information with the
/// ClusterSequence object by calling this function
inline void plugin_associate_extras(std::auto_ptr<Extras> extras_in) {
//_extras = extras_in;
_extras.reset(extras_in.release());
}
/// returns true when the plugin is allowed to run the show.
inline bool plugin_activated() const {return _plugin_activated;}
/// returns a pointer to the extras object (may be null)
const Extras * extras() const {return _extras.get();}
/// allows a plugin to run a templated clustering (nearest-neighbour heuristic)
///
/// This has N^2 behaviour on "good" distance, but a worst case behaviour
/// of N^3 (and many algs trigger the worst case behaviour)
///
///
/// For more details on how this works, see GenBriefJet below
template<class GBJ> void plugin_simple_N2_cluster () {
assert(plugin_activated());
_simple_N2_cluster<GBJ>();
}
public:
//DEP /// set the default (static) jet finder across all current and future
//DEP /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
//DEP /// suppressed in a future release).
//DEP static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
//DEP /// same as above for backward compatibility
//DEP static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
/// \ingroup extra_info
/// \struct history_element
/// a single element in the clustering history
///
/// (see vector _history below).
struct history_element{
int parent1; /// index in _history where first parent of this jet
/// was created (InexistentParent if this jet is an
/// original particle)
int parent2; /// index in _history where second parent of this jet
/// was created (InexistentParent if this jet is an
/// original particle); BeamJet if this history entry
/// just labels the fact that the jet has recombined
/// with the beam)
int child; /// index in _history where the current jet is
/// recombined with another jet to form its child. It
/// is Invalid if this jet does not further
/// recombine.
int jetp_index; /// index in the _jets vector where we will find the
/// PseudoJet object corresponding to this jet
/// (i.e. the jet created at this entry of the
/// history). NB: if this element of the history
/// corresponds to a beam recombination, then
/// jetp_index=Invalid.
double dij; /// the distance corresponding to the recombination
/// at this stage of the clustering.
double max_dij_so_far; /// the largest recombination distance seen
/// so far in the clustering history.
};
enum JetType {Invalid=-3, InexistentParent = -2, BeamJet = -1};
/// allow the user to access the internally stored _jets() array,
/// which contains both the initial particles and the various
/// intermediate and final stages of recombination.
///
/// The first n_particles() entries are the original particles,
/// in the order in which they were supplied to the ClusterSequence
/// constructor. It can be useful to access them for example when
/// examining whether a given input object is part of a specific
/// jet, via the objects_in_jet(...) member function (which only takes
/// PseudoJets that are registered in the ClusterSequence).
///
/// One of the other (internal uses) is related to the fact
/// because we don't seem to be able to access protected elements of
/// the class for an object that is not "this" (at least in case where
/// "this" is of a slightly different kind from the object, both
/// derived from ClusterSequence).
const std::vector<PseudoJet> & jets() const;
/// allow the user to access the raw internal history.
///
/// This is present (as for jets()) in part so that protected
/// derived classes can access this information about other
/// ClusterSequences.
///
/// A user who wishes to follow the details of the ClusterSequence
/// can also make use of this information (and should consult the
/// history_element documentation for more information), but should
/// be aware that these internal structures may evolve in future
/// FastJet versions.
const std::vector<history_element> & history() const;
/// returns the number of particles that were provided to the
/// clustering algorithm (helps the user find their way around the
/// history and jets objects if they weren't paying attention
/// beforehand).
unsigned int n_particles() const;
/// returns a vector of size n_particles() which indicates, for
/// each of the initial particles (in the order in which they were
/// supplied), which of the supplied jets it belongs to; if it does
/// not belong to any of the supplied jets, the index is set to -1;
std::vector<int> particle_jet_indices(const std::vector<PseudoJet> &) const;
/// routine that returns an order in which to read the history
/// such that clusterings that lead to identical jet compositions
/// but different histories (because of degeneracies in the
/// clustering order) will have matching constituents for each
/// matching entry in the unique_history_order.
///
/// The order has the property that an entry's parents will always
/// appear prior to that entry itself.
///
/// Roughly speaking the order is such that we first provide all
/// steps that lead to the final jet containing particle 1; then we
/// have the steps that lead to reconstruction of the jet containing
/// the next-lowest-numbered unclustered particle, etc...
/// [see GPS CCN28-12 for more info -- of course a full explanation
/// here would be better...]
std::vector<int> unique_history_order() const;
/// return the set of particles that have not been clustered. For
/// kt and cam/aachen algorithms this should always be null, but for
/// cone type algorithms it can be non-null;
std::vector<PseudoJet> unclustered_particles() const;
/// Return the list of pseudojets in the ClusterSequence that do not
/// have children (and are not among the inclusive jets). They may
/// result from a clustering step or may be one of the pseudojets
/// returned by unclustered_particles().
std::vector<PseudoJet> childless_pseudojets() const;
/// returns true if the object (jet or particle) is contained by (ie
/// belongs to) this cluster sequence.
///
/// Tests performed: if thejet's interface is this cluster sequence
/// and its cluster history index is in a consistent range.
bool contains(const PseudoJet & object) const;
/// transfer the sequence contained in other_seq into our own;
/// any plugin "extras" contained in the from_seq will be lost
/// from there.
///
/// It also sets the ClusterSequence pointers of the PseudoJets in
/// the history to point to this ClusterSequence
///
/// When specified, the second argument is an action that will be
/// applied on every jets in the resulting ClusterSequence
void transfer_from_sequence(const ClusterSequence & from_seq,
const FunctionOfPseudoJet<PseudoJet> * action_on_jets = 0);
/// retrieve a shared pointer to the wrapper to this ClusterSequence
///
/// this may turn useful if you want to track when this
/// ClusterSequence goes out of scope
const SharedPtr<PseudoJetStructureBase> & structure_shared_ptr() const{
return _structure_shared_ptr;
}
/// the structure type associated with a jet belonging to a ClusterSequence
typedef ClusterSequenceStructure StructureType;
/// This is the function that is automatically called during
/// clustering to print the FastJet banner. Only the first call to
/// this function will result in the printout of the banner. Users
/// may wish to call this function themselves, during the
/// initialization phase of their program, in order to ensure that
/// the banner appears before other output. This call will not
/// affect 3rd-party banners, e.g. those from plugins.
static void print_banner();
/// \cond internal_doc
// [this line must be left as is to hide the doxygen comment]
/// A call to this function modifies the stream used to print
/// banners (by default cout). If a null pointer is passed, banner
/// printout is suppressed. This affects all banners, including
/// those from plugins.
///
/// Please note that if you distribute 3rd party code
/// that links with FastJet, that 3rd party code must not
/// use this call turn off the printing of FastJet banners
/// by default. This requirement reflects the spirit of
/// clause 2c of the GNU Public License (v2), under which
/// FastJet and its plugins are distributed.
static void set_fastjet_banner_stream(std::ostream * ostr) {_fastjet_banner_ostr = ostr;}
// [this line must be left as is to hide the doxygen comment]
/// \endcond
/// returns a pointer to the stream to be used to print banners
/// (cout by default). This function is used by plugins to determine
/// where to direct their banners. Plugins should properly handle
/// the case where the pointer is null.
static std::ostream * fastjet_banner_stream() {return _fastjet_banner_ostr;}
private:
/// \cond internal_doc
/// contains the actual stream to use for banners
static std::ostream * _fastjet_banner_ostr;
/// \endcond
protected:
//DEP static JetAlgorithm _default_jet_algorithm;
JetDefinition _jet_def;
/// transfer the vector<L> of input jets into our own vector<PseudoJet>
/// _jets (with some reserved space for future growth).
template<class L> void _transfer_input_jets(
const std::vector<L> & pseudojets);
/// This is what is called to do all the initialisation and
/// then run the clustering (may be called by various constructors).
/// It assumes _jets contains the momenta to be clustered.
void _initialise_and_run (const JetDefinition & jet_def,
const bool & writeout_combinations);
//// this performs the initialisation, minus the option-decanting
//// stage; for low multiplicity, initialising a few things in the
//// constructor, calling the decant_options_partial() and then this
//// is faster than going through _initialise_and_run.
void _initialise_and_run_no_decant();
//DEP /// This is an alternative routine for initialising and running the
//DEP /// clustering, provided for legacy purposes. The jet finder is that
//DEP /// specified in the static member _default_jet_algorithm.
//DEP void _initialise_and_run (const double & R,
//DEP const Strategy & strategy,
//DEP const bool & writeout_combinations);
/// fills in the various member variables with "decanted" options from
/// the jet_definition and writeout_combinations variables
void _decant_options(const JetDefinition & jet_def,
const bool & writeout_combinations);
/// assuming that the jet definition, writeout_combinations and
/// _structure_shared_ptr have been set (e.g. in an initialiser list
/// in the constructor), it handles the remaining decanting of
/// options.
void _decant_options_partial();
/// fill out the history (and jet cross refs) related to the initial
/// set of jets (assumed already to have been "transferred"),
/// without any clustering
void _fill_initial_history();
/// carry out the recombination between the jets numbered jet_i and
/// jet_j, at distance scale dij; return the index newjet_k of the
/// result of the recombination of i and j.
void _do_ij_recombination_step(const int & jet_i, const int & jet_j,
const double & dij, int & newjet_k);
/// carry out an recombination step in which _jets[jet_i] merges with
/// the beam,
void _do_iB_recombination_step(const int & jet_i, const double & diB);
/// every time a jet is added internally during clustering, this
/// should be called to set the jet's structure shared ptr to point
/// to the CS (and the count of internally associated objects is
/// also updated). This should not be called outside construction of
/// a CS object.
void _set_structure_shared_ptr(PseudoJet & j);
/// make sure that the CS's internal tally of the use count matches
/// that of the _structure_shared_ptr
void _update_structure_use_count();
/// This contains the physical PseudoJets; for each PseudoJet one
/// can find the corresponding position in the _history by looking
/// at _jets[i].cluster_hist_index().
std::vector<PseudoJet> _jets;
/// this vector will contain the branching history; for each stage,
/// _history[i].jetp_index indicates where to look in the _jets
/// vector to get the physical PseudoJet.
std::vector<history_element> _history;
/// set subhist to be a set pointers to history entries corresponding to the
/// subjets of this jet; one stops going working down through the
/// subjets either when
/// - there is no further to go
/// - one has found maxjet entries
/// - max_dij_so_far <= dcut
/// By setting maxjet=0 one can use just dcut; by setting dcut<0
/// one can use jet maxjet
void get_subhist_set(std::set<const history_element*> & subhist,
const PseudoJet & jet, double dcut, int maxjet) const;
bool _writeout_combinations;
int _initial_n;
double _Rparam, _R2, _invR2;
double _Qtot;
Strategy _strategy;
JetAlgorithm _jet_algorithm;
SharedPtr<PseudoJetStructureBase> _structure_shared_ptr; //< will actually be of type ClusterSequenceStructure
int _structure_use_count_after_construction; //< info of use when CS handles its own memory
/// if true then the CS will delete itself when the last external
/// object referring to it disappears. It is mutable so as to ensure
/// that signal_imminent_self_deletion() [const] can make relevant
/// changes.
mutable bool _deletes_self_when_unused;
private:
bool _plugin_activated;
//std::auto_ptr<Extras> _extras; // things the plugin might want to add
SharedPtr<Extras> _extras; // things the plugin might want to add
void _really_dumb_cluster ();
void _delaunay_cluster ();
//void _simple_N2_cluster ();
template<class BJ> void _simple_N2_cluster ();
void _tiled_N2_cluster ();
void _faster_tiled_N2_cluster ();
//
void _minheap_faster_tiled_N2_cluster();
// things needed specifically for Cambridge with Chan's 2D closest
// pairs method
void _CP2DChan_cluster();
void _CP2DChan_cluster_2pi2R ();
void _CP2DChan_cluster_2piMultD ();
void _CP2DChan_limited_cluster(double D);
void _do_Cambridge_inclusive_jets();
// NSqrtN method for C/A
void _fast_NsqrtN_cluster();
void _add_step_to_history(const int & step_number, const int & parent1,
const int & parent2, const int & jetp_index,
const double & dij);
/// internal routine associated with the construction of the unique
/// history order (following children in the tree)
void _extract_tree_children(int pos, std::valarray<bool> &,
const std::valarray<int> &, std::vector<int> &) const;
/// internal routine associated with the construction of the unique
/// history order (following parents in the tree)
void _extract_tree_parents (int pos, std::valarray<bool> &,
const std::valarray<int> &, std::vector<int> &) const;
// these will be useful shorthands in the Voronoi-based code
typedef std::pair<int,int> TwoVertices;
typedef std::pair<double,TwoVertices> DijEntry;
typedef std::multimap<double,TwoVertices> DistMap;
/// currently used only in the Voronoi based code
void _add_ktdistance_to_map(const int & ii,
DistMap & DijMap,
const DynamicNearestNeighbours * DNN);
/// will be set by default to be true for the first run
static bool _first_time;
/// record the number of warnings provided about the exclusive
/// algorithm -- so that we don't print it out more than a few
/// times.
static int _n_exclusive_warnings;
/// the limited warning member for notification of user that
/// their requested strategy has been overridden (usually because
/// they have R>2pi and not all strategies work then)
static LimitedWarning _changed_strategy_warning;
//----------------------------------------------------------------------
/// the fundamental structure which contains the minimal info about
/// a jet, as needed for our plain N^2 algorithm -- the idea is to
/// put all info that will be accessed N^2 times into an array of
/// BriefJets...
struct BriefJet {
double eta, phi, kt2, NN_dist;
BriefJet * NN;
int _jets_index;
};
/// structure analogous to BriefJet, but with the extra information
/// needed for dealing with tiles
class TiledJet {
public:
double eta, phi, kt2, NN_dist;
TiledJet * NN, *previous, * next;
int _jets_index, tile_index, diJ_posn;
// routines that are useful in the minheap version of tiled
// clustering ("misuse" the otherwise unused diJ_posn, so as
// to indicate whether jets need to have their minheap entries
// updated).
inline void label_minheap_update_needed() {diJ_posn = 1;}
inline void label_minheap_update_done() {diJ_posn = 0;}
inline bool minheap_update_needed() const {return diJ_posn==1;}
};
//-- some of the functions that follow are templates and will work
//as well for briefjet and tiled jets
/// set the kinematic and labelling info for jeta so that it corresponds
/// to _jets[_jets_index]
template <class J> void _bj_set_jetinfo( J * const jet,
const int _jets_index) const;
/// "remove" this jet, which implies updating links of neighbours and
/// perhaps modifying the tile structure
void _bj_remove_from_tiles( TiledJet * const jet) const;
/// return the distance between two BriefJet objects
template <class J> double _bj_dist(const J * const jeta,
const J * const jetb) const;
// return the diJ (multiplied by _R2) for this jet assuming its NN
// info is correct
template <class J> double _bj_diJ(const J * const jeta) const;
/// for testing purposes only: if in the range head--tail-1 there is a
/// a jet which corresponds to hist_index in the history, then
/// return a pointer to that jet; otherwise return tail.
template <class J> inline J * _bj_of_hindex(
const int hist_index,
J * const head, J * const tail)
const {
J * res;
for(res = head; res<tail; res++) {
if (_jets[res->_jets_index].cluster_hist_index() == hist_index) {break;}
}
return res;
}
//-- remaining functions are different in various cases, so we
// will use templates but are not sure if they're useful...
/// updates (only towards smaller distances) the NN for jeta without checking
/// whether in the process jeta itself might be a new NN of one of
/// the jets being scanned -- span the range head to tail-1 with
/// assumption that jeta is not contained in that range
template <class J> void _bj_set_NN_nocross(J * const jeta,
J * const head, const J * const tail) const;
/// reset the NN for jeta and DO check whether in the process jeta
/// itself might be a new NN of one of the jets being scanned --
/// span the range head to tail-1 with assumption that jeta is not
/// contained in that range
template <class J> void _bj_set_NN_crosscheck(J * const jeta,
J * const head, const J * const tail) const;
/// number of neighbours that a tile will have (rectangular geometry
/// gives 9 neighbours).
static const int n_tile_neighbours = 9;
//----------------------------------------------------------------------
/// The fundamental structures to be used for the tiled N^2 algorithm
/// (see CCN27-44 for some discussion of pattern of tiling)
struct Tile {
/// pointers to neighbouring tiles, including self
Tile * begin_tiles[n_tile_neighbours];
/// neighbouring tiles, excluding self
Tile ** surrounding_tiles;
/// half of neighbouring tiles, no self
Tile ** RH_tiles;
/// just beyond end of tiles
Tile ** end_tiles;
/// start of list of BriefJets contained in this tile
TiledJet * head;
/// sometimes useful to be able to tag a tile
bool tagged;
};
std::vector<Tile> _tiles;
double _tiles_eta_min, _tiles_eta_max;
double _tile_size_eta, _tile_size_phi;
int _n_tiles_phi,_tiles_ieta_min,_tiles_ieta_max;
// reasonably robust return of tile index given ieta and iphi, in particular
// it works even if iphi is negative
inline int _tile_index (int ieta, int iphi) const {
// note that (-1)%n = -1 so that we have to add _n_tiles_phi
// before performing modulo operation
return (ieta-_tiles_ieta_min)*_n_tiles_phi
+ (iphi+_n_tiles_phi) % _n_tiles_phi;
}
// routines for tiled case, including some overloads of the plain
// BriefJet cases
int _tile_index(const double & eta, const double & phi) const;
void _tj_set_jetinfo ( TiledJet * const jet, const int _jets_index);
void _bj_remove_from_tiles(TiledJet * const jet);
void _initialise_tiles();
void _print_tiles(TiledJet * briefjets ) const;
void _add_neighbours_to_tile_union(const int tile_index,
std::vector<int> & tile_union, int & n_near_tiles) const;
void _add_untagged_neighbours_to_tile_union(const int tile_index,
std::vector<int> & tile_union, int & n_near_tiles);
//----------------------------------------------------------------------
/// fundamental structure for e+e- clustering
struct EEBriefJet {
double NN_dist; // obligatorily present
double kt2; // obligatorily present == E^2 in general
EEBriefJet * NN; // must be present too
int _jets_index; // must also be present!
//...........................................................
double nx, ny, nz; // our internal storage for fast distance calcs
};
/// to help instantiation (fj 2.4.0; did not quite work on gcc 33 and os x 10.3?)
//void _dummy_N2_cluster_instantiation();
/// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
void _simple_N2_cluster_BriefJet();
/// to avoid issues with template instantiation (OS X 10.3, gcc 3.3)
void _simple_N2_cluster_EEBriefJet();
};
//**********************************************************************
//************** START OF INLINE MATERIAL ******************
//**********************************************************************
//----------------------------------------------------------------------
// Transfer the initial jets into our internal structure
template<class L> void ClusterSequence::_transfer_input_jets(
const std::vector<L> & pseudojets) {
// this will ensure that we can point to jets without difficulties
// arising.
_jets.reserve(pseudojets.size()*2);
// insert initial jets this way so that any type L that can be
// converted to a pseudojet will work fine (basically PseudoJet
// and any type that has [] subscript access to the momentum
// components, such as CLHEP HepLorentzVector).
for (unsigned int i = 0; i < pseudojets.size(); i++) {
_jets.push_back(pseudojets[i]);}
}
// //----------------------------------------------------------------------
// // initialise from some generic type... Has to be made available
// // here in order for it the template aspect of it to work...
// template<class L> ClusterSequence::ClusterSequence (
// const std::vector<L> & pseudojets,
// const double & R,
// const Strategy & strategy,
// const bool & writeout_combinations) {
//
// // transfer the initial jets (type L) into our own array
// _transfer_input_jets(pseudojets);
//
// // run the clustering
// _initialise_and_run(R,strategy,writeout_combinations);
// }
//----------------------------------------------------------------------
/// constructor of a jet-clustering sequence from a vector of
/// four-momenta, with the jet definition specified by jet_def
template<class L> ClusterSequence::ClusterSequence (
const std::vector<L> & pseudojets,
const JetDefinition & jet_def_in,
const bool & writeout_combinations) :
_jet_def(jet_def_in), _writeout_combinations(writeout_combinations),
_structure_shared_ptr(new ClusterSequenceStructure(this))
{
// transfer the initial jets (type L) into our own array
_transfer_input_jets(pseudojets);
// transfer the remaining options
_decant_options_partial();
// run the clustering
_initialise_and_run_no_decant();
}
inline const std::vector<PseudoJet> & ClusterSequence::jets () const {
return _jets;
}
inline const std::vector<ClusterSequence::history_element> & ClusterSequence::history () const {
return _history;
}
inline unsigned int ClusterSequence::n_particles() const {return _initial_n;}
//----------------------------------------------------------------------
template <class J> inline void ClusterSequence::_bj_set_jetinfo(
J * const jetA, const int _jets_index) const {
jetA->eta = _jets[_jets_index].rap();
jetA->phi = _jets[_jets_index].phi_02pi();
jetA->kt2 = jet_scale_for_algorithm(_jets[_jets_index]);
jetA->_jets_index = _jets_index;
// initialise NN info as well
jetA->NN_dist = _R2;
jetA->NN = NULL;
}
//----------------------------------------------------------------------
template <class J> inline double ClusterSequence::_bj_dist(
const J * const jetA, const J * const jetB) const {
double dphi = std::abs(jetA->phi - jetB->phi);
double deta = (jetA->eta - jetB->eta);
if (dphi > pi) {dphi = twopi - dphi;}
return dphi*dphi + deta*deta;
}
//----------------------------------------------------------------------
template <class J> inline double ClusterSequence::_bj_diJ(const J * const jet) const {
double kt2 = jet->kt2;
if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
return jet->NN_dist * kt2;
}
//----------------------------------------------------------------------
// set the NN for jet without checking whether in the process you might
// have discovered a new nearest neighbour for another jet
template <class J> inline void ClusterSequence::_bj_set_NN_nocross(
J * const jet, J * const head, const J * const tail) const {
double NN_dist = _R2;
J * NN = NULL;
if (head < jet) {
for (J * jetB = head; jetB != jet; jetB++) {
double dist = _bj_dist(jet,jetB);
if (dist < NN_dist) {
NN_dist = dist;
NN = jetB;
}
}
}
if (tail > jet) {
for (J * jetB = jet+1; jetB != tail; jetB++) {
double dist = _bj_dist(jet,jetB);
if (dist < NN_dist) {
NN_dist = dist;
NN = jetB;
}
}
}
jet->NN = NN;
jet->NN_dist = NN_dist;
}
//----------------------------------------------------------------------
template <class J> inline void ClusterSequence::_bj_set_NN_crosscheck(J * const jet,
J * const head, const J * const tail) const {
double NN_dist = _R2;
J * NN = NULL;
for (J * jetB = head; jetB != tail; jetB++) {
double dist = _bj_dist(jet,jetB);
if (dist < NN_dist) {
NN_dist = dist;
NN = jetB;
}
if (dist < jetB->NN_dist) {
jetB->NN_dist = dist;
jetB->NN = jet;
}
}
jet->NN = NN;
jet->NN_dist = NN_dist;
}
FASTJET_END_NAMESPACE
#endif // __FASTJET_CLUSTERSEQUENCE_HH__
|