+++ /dev/null
-//STARTHEADER
-// $Id: ClusterSequence.hh 1435 2009-02-12 21:11:04Z salam $
-//
-// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
-//
-//----------------------------------------------------------------------
-// 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, write to the Free Software
-// Foundation, Inc.:
-// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//----------------------------------------------------------------------
-//ENDHEADER
-
-
-//----------------------------------------------------------------------
-// here's where we put the main page for fastjet (as explained in the
-// Doxygen faq)
-//......................................................................
-/*! \mainpage FastJet code documentation
- *
- * These pages provide automatically generated documentation for the
- * FastJet package.
- *
- * For further information and normal documentation, see the main <a
- * href="http://www.lpthe.jussieu.fr/~salam/fastjet">FastJet</a> page.
- */
-//----------------------------------------------------------------------
-
-#ifndef __FASTJET_CLUSTERSEQUENCE_HH__
-#define __FASTJET_CLUSTERSEQUENCE_HH__
-
-#include<vector>
-#include<map>
-#include "fastjet/internal/DynamicNearestNeighbours.hh"
-#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"
-
-namespace fastjet { // defined in fastjet/internal/base.hh
-
-
-/// deals with clustering
-class ClusterSequence {
-
-
- public:
-
- /// default constructor
- ClusterSequence () {}
-
- /// 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);
-
- // 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.
- std::vector<PseudoJet> exclusive_jets (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}).
- 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_ymerge_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 n subjets (or all constituents if there are fewer
- /// than n).
- ///
- /// requires n ln n time
- std::vector<PseudoJet> exclusive_subjets (const PseudoJet & jet,
- int nsub) const;
-
- /// return the dij that was present in the merging nsub+1 -> nsub
- /// subjets inside this 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.
- 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 _Q;}
- /// return Q()^2
- double Q2() const {return _Q*_Q;}
-
- /// 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;}
- std::string strategy_string () const;
-
- /// return a reference to the jet definition
- const JetDefinition & jet_def() const {return _jet_def;}
-
- /// 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);
- }
-
- /// 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 associated 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;
- }
-
- /// 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>();
- }
-
- //----------------------------------------------------------------------
- /// class to help with a generic clustering sequence
- class GenBriefJet {
- public:
- /// function that initialises the GenBriefJet given a PseudoJet.
- ///
- /// In a derived class, this member has a responsability to call
- ///
- /// - set_scale_squared
- /// - set_geom_iB
- ///
- /// The clustering will be performed by finding the minimum of
- ///
- /// diB = scale_squared[i] * geom_iB * _invR2
- /// dij = min(scale_squared[i],scale_squared[j]) * geom_ij * _invR2
- ///
- virtual void init(const PseudoJet & jet) = 0;
-
- /// Returns the "geometric" part of distance between this jet
- /// and jet_j
- virtual double geom_ij(const GenBriefJet * jet_j) const = 0;
-
- void set_scale_squared(double scale_squared) {kt2 = scale_squared;}
- void set_geom_iB(double diB) {NN_dist = diB; NN = NULL;}
-
- public: // formally public: but users should think of it as private!
- double NN_dist; // dij
- double kt2; // squared scale
- GenBriefJet * NN; // pointer to nearest neighbour
- int _jets_index; // index of this jet
- };
-
-
-public:
- /// set the default (static) jet finder across all current and future
- /// ClusterSequence objects -- deprecated and obsolescent (i.e. may be
- /// suppressed in a future release).
- static void set_jet_algorithm (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
- /// same as above for backward compatibility
- static void set_jet_finder (JetAlgorithm jet_algorithm) {_default_jet_algorithm = jet_algorithm;}
-
-
- /// 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;
-
- /// transfer the sequence contained in other_seq into our own;
- /// any plugin "extras" contained in the from_seq will be lost
- /// from there.
- void transfer_from_sequence(ClusterSequence & from_seq);
-
-
-protected:
- static JetAlgorithm _default_jet_algorithm;
- JetDefinition _jet_def;
-
- /// returns true if the jet has a history index contained within
- /// the range of this CS
- bool _potentially_valid(const PseudoJet & jet) const {
- return jet.cluster_hist_index() >= 0
- && jet.cluster_hist_index() < int(_history.size());
- }
-
- /// 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 the routine that will do all the initialisation and
- /// then run the clustering (may be called by various constructors).
- /// It assumes _jets contains the momenta to be clustered.
- void _initialise_and_run (const JetDefinition & jet_def,
- const bool & writeout_combinations);
-
- /// This is an alternative routine for initialising and running the
- /// clustering, provided for legacy purposes. The jet finder is that
- /// specified in the static member _default_jet_algorithm.
- void _initialise_and_run (const double & R,
- const Strategy & strategy,
- 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);
-
- /// 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);
-
-
- /// 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 _Q;
- Strategy _strategy;
- JetAlgorithm _jet_algorithm;
-
-
- private:
-
- bool _plugin_activated;
- std::auto_ptr<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();
-
- 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);
-
- /// for making sure the user knows what it is they're running...
- void _print_banner();
- /// 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 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
- void _dummy_N2_cluster_instantiation();
-
-};
-
-
-//**********************************************************************
-//************** 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,
- const bool & writeout_combinations) {
-
- // transfer the initial jets (type L) into our own array
- _transfer_input_jets(pseudojets);
-
- // run the clustering
- _initialise_and_run(jet_def,writeout_combinations);
-}
-
-
-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 namespace
-
-#endif // __FASTJET_CLUSTERSEQUENCE_HH__
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff