FastJet 3.0.2
|
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__