FastJet  3.0.6
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__
mean_pt_strategies
enum providing a variety of tentative strategies for estimating the background (e.g.
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:764
bool has_dangerous_particles() const
returns true if there are any particles whose transverse momentum if so low that there's a risk of th...
virtual double area_error(const PseudoJet &jet) const
return the error (uncertainty) associated with the determination of the area of this jet; this base c...
Like ClusterSequence with computation of the active jet area.
virtual double area(const PseudoJet &jet) const
return the area associated with the given jet; this base class returns 0.
base class that sets interface for extensions of ClusterSequence that provide information about the a...
Like ClusterSequence with computation of the active jet area with the addition of explicit ghosts...
std::valarray< double > _average_area
child classes benefit from having these at their disposal
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
virtual PseudoJet area_4vector(const PseudoJet &jet) const
return a PseudoJet whose 4-vector is defined by the following integral
Parameters to configure the computation of jet areas using ghosts.
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:65
class that is intended to hold a full definition of the jet clusterer