FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ClusterSequenceActiveArea.hh
1 //STARTHEADER
2 // $Id: ClusterSequenceActiveArea.hh 2687 2011-11-14 11:17:51Z soyez $
3 //
4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #ifndef __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
30 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__
31 
32 
33 #include "fastjet/PseudoJet.hh"
34 #include "fastjet/ClusterSequenceAreaBase.hh"
35 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
36 #include<iostream>
37 #include<vector>
38 
39 //------------ backwards compatibility with version 2.1 -------------
40 // for backwards compatibility make ActiveAreaSpec name available
41 //#include "fastjet/ActiveAreaSpec.hh"
42 //#include "fastjet/ClusterSequenceWithArea.hh"
43 //--------------------------------------------------------------------
44 
45 
46 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
47 
48 //using namespace std;
49 
50 /// @ingroup sec_area_classes
51 /// \class ClusterSequenceActiveArea
52 /// Like ClusterSequence with computation of the active jet area
53 ///
54 /// Class that behaves essentially like ClusterSequence except
55 /// that it also provides access to the area of a jet (which
56 /// will be a random quantity... Figure out what to do about seeds
57 /// later...)
58 ///
59 /// This class should not be used directly. Rather use
60 /// ClusterSequenceArea with the appropriate AreaDefinition
62 public:
63 
64  /// default constructor
66 
67  /// constructor based on JetDefinition and GhostedAreaSpec
68  template<class L> ClusterSequenceActiveArea
69  (const std::vector<L> & pseudojets,
70  const JetDefinition & jet_def_in,
71  const GhostedAreaSpec & ghost_spec,
72  const bool & writeout_combinations = false) ;
73 
74  virtual double area (const PseudoJet & jet) const {
75  return _average_area[jet.cluster_hist_index()];};
76  virtual double area_error (const PseudoJet & jet) const {
77  return _average_area2[jet.cluster_hist_index()];};
78 
79  virtual PseudoJet area_4vector (const PseudoJet & jet) const {
80  return _average_area_4vector[jet.cluster_hist_index()];};
81 
82  /// enum providing a variety of tentative strategies for estimating
83  /// the background (e.g. non-jet) activity in a highly populated event; the
84  /// one that has been most extensively tested is median.
85  enum mean_pt_strategies{median=0, non_ghost_median, pttot_over_areatot,
86  pttot_over_areatot_cut, mean_ratio_cut, play,
87  median_4vector};
88 
89  /// return the transverse momentum per unit area according to one
90  /// of the above strategies; for some strategies (those with "cut"
91  /// in their name) the parameter "range" allows one to exclude a
92  /// subset of the jets for the background estimation, those that
93  /// have pt/area > median(pt/area)*range.
94  ///
95  /// NB: This call is OBSOLETE; use media_pt_per_unit_area from the
96  // ClusterSequenceAreaBase class instead
97  double pt_per_unit_area(mean_pt_strategies strat=median,
98  double range=2.0 ) const;
99 
100  // following code removed -- now dealt with by AreaBase class (and
101  // this definition here conflicts with it).
102 // /// fits a form pt_per_unit_area(y) = a + b*y^2 in the range
103 // /// abs(y)<raprange (for negative raprange, it defaults to
104 // /// _safe_rap_for_area).
105 // void parabolic_pt_per_unit_area(double & a,double & b, double raprange=-1.0,
106 // double exclude_above=-1.0,
107 // bool use_area_4vector=false ) const;
108 //
109  /// rewrite the empty area from the parent class, so as to use
110  /// all info at our disposal
111  /// return the total area, corresponding to a given Selector, that
112  /// consists of ghost jets or unclustered ghosts
113  ///
114  /// The selector passed as an argument needs to apply jet by jet.
115  virtual double empty_area(const Selector & selector) const;
116 
117  /// return the true number of empty jets (replaces
118  /// ClusterSequenceAreaBase::n_empty_jets(...))
119  virtual double n_empty_jets(const Selector & selector) const;
120 
121 protected:
122  void _resize_and_zero_AA ();
123  void _initialise_AA(const JetDefinition & jet_def,
124  const GhostedAreaSpec & ghost_spec,
125  const bool & writeout_combinations,
126  bool & continue_running);
127 
128  void _run_AA(const GhostedAreaSpec & ghost_spec);
129 
130  void _postprocess_AA(const GhostedAreaSpec & ghost_spec);
131 
132  /// does the initialisation and running specific to the active
133  /// areas class
134  void _initialise_and_run_AA (const JetDefinition & jet_def,
135  const GhostedAreaSpec & ghost_spec,
136  const bool & writeout_combinations = false);
137 
138  /// transfer the history (and jet-momenta) from clust_seq to our
139  /// own internal structure while removing ghosts
140  void _transfer_ghost_free_history(
141  const ClusterSequenceActiveAreaExplicitGhosts & clust_seq);
142 
143 
144  /// transfer areas from the ClusterSequenceActiveAreaExplicitGhosts
145  /// object into our internal area bookkeeping...
146  void _transfer_areas(const std::vector<int> & unique_hist_order,
148 
149  /// child classes benefit from having these at their disposal
150  std::valarray<double> _average_area, _average_area2;
151  std::valarray<PseudoJet> _average_area_4vector;
152 
153  /// returns true if there are any particles whose transverse momentum
154  /// if so low that there's a risk of the ghosts having modified the
155  /// clustering sequence
156  bool has_dangerous_particles() const {return _has_dangerous_particles;}
157 
158 private:
159 
160 
161  double _non_jet_area, _non_jet_area2, _non_jet_number;
162 
163  double _maxrap_for_area; // max rap where we put ghosts
164  double _safe_rap_for_area; // max rap where we trust jet areas
165 
166  bool _has_dangerous_particles;
167 
168 
169  /// routine for extracting the tree in an order that will be independent
170  /// of any degeneracies in the recombination sequence that don't
171  /// affect the composition of the final jets
172  void _extract_tree(std::vector<int> &) const;
173  /// do the part of the extraction associated with pos, working
174  /// through its children and their parents
175  void _extract_tree_children(int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const;
176  /// do the part of the extraction associated with the parents of pos.
177  void _extract_tree_parents (int pos, std::valarray<bool> &, const std::valarray<int> &, std::vector<int> &) const;
178 
179  /// check if two jets have the same momentum to within the
180  /// tolerance (and if pt's are not the same we're forgiving and
181  /// look to see if the energy is the same)
182  void _throw_unless_jets_have_same_perp_or_E(const PseudoJet & jet,
183  const PseudoJet & refjet,
184  double tolerance,
185  const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq
186  ) const;
187 
188  /// since we are playing nasty games with seeds, we should warn
189  /// the user a few times
190  //static int _n_seed_warnings;
191  //const static int _max_seed_warnings = 10;
192 
193  // record the number of repeats
194  int _ghost_spec_repeat;
195 
196  /// a class for our internal storage of ghost jets
197  class GhostJet : public PseudoJet {
198  public:
199  GhostJet(const PseudoJet & j, double a) : PseudoJet(j), area(a){}
200  double area;
201  };
202 
203  std::vector<GhostJet> _ghost_jets;
204  std::vector<GhostJet> _unclustered_ghosts;
205 };
206 
207 
208 
209 
210 template<class L> ClusterSequenceActiveArea::ClusterSequenceActiveArea
211 (const std::vector<L> & pseudojets,
212  const JetDefinition & jet_def_in,
213  const GhostedAreaSpec & ghost_spec,
214  const bool & writeout_combinations) {
215 
216  // transfer the initial jets (type L) into our own array
217  _transfer_input_jets(pseudojets);
218 
219  // run the clustering for active areas
220  _initialise_and_run_AA(jet_def_in, ghost_spec, writeout_combinations);
221 
222 }
223 
224 
225 
226 FASTJET_END_NAMESPACE
227 
228 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREA_HH__