FastJet 3.0.2
ClusterSequenceActiveArea.hh
00001 //STARTHEADER
00002 // $Id: ClusterSequenceActiveArea.hh 2687 2011-11-14 11:17:51Z soyez $
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 #ifndef __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
00030 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
00031 
00032 
00033 #include "fastjet/PseudoJet.hh"
00034 #include "fastjet/ClusterSequenceAreaBase.hh"
00035 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
00036 #include<iostream>
00037 #include<vector>
00038 
00039 //------------ backwards compatibility with version 2.1 -------------
00040 // for backwards compatibility make ActiveAreaSpec name available
00041 //#include "fastjet/ActiveAreaSpec.hh"
00042 //#include "fastjet/ClusterSequenceWithArea.hh"
00043 //--------------------------------------------------------------------
00044 
00045 
00046 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00047 
00048 //using namespace std;
00049 
00050 /// @ingroup sec_area_classes
00051 /// \class ClusterSequenceActiveArea
00052 /// Like ClusterSequence with computation of the active jet area
00053 ///
00054 /// Class that behaves essentially like ClusterSequence except
00055 /// that it also provides access to the area of a jet (which
00056 /// will be a random quantity... Figure out what to do about seeds 
00057 /// later...)
00058 ///
00059 /// This class should not be used directly. Rather use
00060 /// ClusterSequenceArea with the appropriate AreaDefinition
00061 class ClusterSequenceActiveArea : public ClusterSequenceAreaBase {
00062 public:
00063 
00064   /// default constructor
00065   ClusterSequenceActiveArea() {}
00066 
00067   /// constructor based on JetDefinition and GhostedAreaSpec
00068   template<class L> ClusterSequenceActiveArea
00069          (const std::vector<L> & pseudojets, 
00070           const JetDefinition & jet_def_in,
00071           const GhostedAreaSpec & ghost_spec,
00072           const bool & writeout_combinations = false) ;
00073 
00074   virtual double area (const PseudoJet & jet) const {
00075                              return _average_area[jet.cluster_hist_index()];};
00076   virtual double area_error (const PseudoJet & jet) const {
00077                              return _average_area2[jet.cluster_hist_index()];};
00078 
00079   virtual PseudoJet area_4vector (const PseudoJet & jet) const {
00080                     return _average_area_4vector[jet.cluster_hist_index()];};
00081 
00082   /// enum providing a variety of tentative strategies for estimating
00083   /// the background (e.g. non-jet) activity in a highly populated event; the
00084   /// one that has been most extensively tested is median.
00085   enum mean_pt_strategies{median=0, non_ghost_median, pttot_over_areatot, 
00086                           pttot_over_areatot_cut, mean_ratio_cut, play,
00087                           median_4vector};
00088 
00089   /// return the transverse momentum per unit area according to one
00090   /// of the above strategies; for some strategies (those with "cut"
00091   /// in their name) the parameter "range" allows one to exclude a
00092   /// subset of the jets for the background estimation, those that
00093   /// have pt/area > median(pt/area)*range.
00094   ///
00095   /// NB: This call is OBSOLETE; use media_pt_per_unit_area from the
00096   //      ClusterSequenceAreaBase class instead
00097   double pt_per_unit_area(mean_pt_strategies strat=median, 
00098                           double range=2.0 ) const;
00099 
00100   // following code removed -- now dealt with by AreaBase class (and
00101   // this definition here conflicts with it).
00102 //   /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range
00103 //   /// abs(y)<raprange (for negative raprange, it defaults to
00104 //   /// _safe_rap_for_area).
00105 //   void parabolic_pt_per_unit_area(double & a,double & b, double raprange=-1.0,
00106 //                                double exclude_above=-1.0, 
00107 //                                bool use_area_4vector=false ) const;
00108 // 
00109   /// rewrite the empty area from the parent class, so as to use
00110   /// all info at our disposal
00111   /// return the total area, corresponding to a given Selector, that
00112   /// consists of ghost jets or unclustered ghosts
00113   ///
00114   /// The selector passed as an argument needs to apply jet by jet.
00115   virtual double empty_area(const Selector & selector) const;
00116 
00117   /// return the true number of empty jets (replaces
00118   /// ClusterSequenceAreaBase::n_empty_jets(...))
00119   virtual double n_empty_jets(const Selector & selector) const;
00120 
00121 protected:
00122   void _resize_and_zero_AA ();
00123   void _initialise_AA(const JetDefinition & jet_def,
00124                       const GhostedAreaSpec & ghost_spec,
00125                       const bool & writeout_combinations,
00126                       bool & continue_running);
00127 
00128   void _run_AA(const GhostedAreaSpec & ghost_spec);
00129 
00130   void _postprocess_AA(const GhostedAreaSpec & ghost_spec);
00131 
00132   /// does the initialisation and running specific to the active
00133   /// areas class
00134   void _initialise_and_run_AA (const JetDefinition & jet_def,
00135                                const GhostedAreaSpec & ghost_spec,
00136                                const bool & writeout_combinations = false);
00137 
00138   /// transfer the history (and jet-momenta) from clust_seq to our
00139   /// own internal structure while removing ghosts
00140   void _transfer_ghost_free_history(
00141            const ClusterSequenceActiveAreaExplicitGhosts & clust_seq);
00142 
00143 
00144   /// transfer areas from the ClusterSequenceActiveAreaExplicitGhosts
00145   /// object into our internal area bookkeeping...
00146   void _transfer_areas(const std::vector<int> & unique_hist_order, 
00147                        const ClusterSequenceActiveAreaExplicitGhosts & );
00148 
00149   /// child classes benefit from having these at their disposal
00150   std::valarray<double> _average_area, _average_area2;
00151   std::valarray<PseudoJet> _average_area_4vector;
00152 
00153   /// returns true if there are any particles whose transverse momentum
00154   /// if so low that there's a risk of the ghosts having modified the
00155   /// clustering sequence
00156   bool has_dangerous_particles() const {return _has_dangerous_particles;}
00157 
00158 private:
00159 
00160 
00161   double           _non_jet_area, _non_jet_area2, _non_jet_number;
00162 
00163   double _maxrap_for_area; // max rap where we put ghosts
00164   double _safe_rap_for_area; // max rap where we trust jet areas
00165 
00166   bool   _has_dangerous_particles; 
00167 
00168 
00169   /// routine for extracting the tree in an order that will be independent
00170   /// of any degeneracies in the recombination sequence that don't
00171   /// affect the composition of the final jets
00172   void _extract_tree(std::vector<int> &) const;
00173   /// do the part of the extraction associated with pos, working
00174   /// through its children and their parents
00175   void _extract_tree_children(int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const;
00176   /// do the part of the extraction associated with the parents of pos.
00177   void _extract_tree_parents (int pos, std::valarray<bool> &, const std::valarray<int> &,  std::vector<int> &) const;
00178 
00179   /// check if two jets have the same momentum to within the
00180   /// tolerance (and if pt's are not the same we're forgiving and
00181   /// look to see if the energy is the same)
00182   void _throw_unless_jets_have_same_perp_or_E(const PseudoJet & jet, 
00183                                               const PseudoJet & refjet, 
00184                                               double tolerance,
00185              const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
00186                                               ) const;
00187 
00188   /// since we are playing nasty games with seeds, we should warn
00189   /// the user a few times
00190   //static int _n_seed_warnings;
00191   //const static int _max_seed_warnings = 10;
00192 
00193   // record the number of repeats
00194   int _ghost_spec_repeat;
00195 
00196   /// a class for our internal storage of ghost jets
00197   class GhostJet : public PseudoJet {
00198   public:
00199     GhostJet(const PseudoJet & j, double a) : PseudoJet(j), area(a){}
00200     double area;
00201   };
00202 
00203   std::vector<GhostJet> _ghost_jets;
00204   std::vector<GhostJet> _unclustered_ghosts;
00205 };
00206 
00207 
00208 
00209 
00210 template<class L> ClusterSequenceActiveArea::ClusterSequenceActiveArea 
00211 (const std::vector<L> & pseudojets, 
00212  const JetDefinition & jet_def_in,
00213  const GhostedAreaSpec & ghost_spec,
00214  const bool & writeout_combinations) {
00215 
00216   // transfer the initial jets (type L) into our own array
00217   _transfer_input_jets(pseudojets);
00218 
00219   // run the clustering for active areas
00220   _initialise_and_run_AA(jet_def_in, ghost_spec, writeout_combinations);
00221 
00222 }
00223 
00224 
00225   
00226 FASTJET_END_NAMESPACE
00227 
00228 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends