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