FastJet  3.0.6
ClusterSequenceActiveAreaExplicitGhosts.hh
1 //STARTHEADER
2 // $Id: ClusterSequenceActiveAreaExplicitGhosts.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_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
30 #define __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
31 
32 #include "fastjet/PseudoJet.hh"
33 #include "fastjet/ClusterSequenceAreaBase.hh"
34 #include "fastjet/GhostedAreaSpec.hh"
35 #include "fastjet/LimitedWarning.hh"
36 #include<iostream>
37 #include<vector>
38 #include <cstdio>
39 
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41 
42 //======================================================================
43 /// @ingroup sec_area_classes
44 /// \class ClusterSequenceActiveAreaExplicitGhosts
45 /// Like ClusterSequence with computation of the active jet area with the
46 /// addition of explicit ghosts
47 ///
48 /// Class that behaves essentially like ClusterSequence except
49 /// that it also provides access to the area of a jet (which
50 /// will be a random quantity... Figure out what to do about seeds
51 /// later...)
52 ///
53 /// This class should not be used directly. Rather use
54 /// ClusterSequenceArea with the appropriate AreaDefinition
57 public:
58  /// constructor using a GhostedAreaSpec to specify how the area is
59  /// to be measured
61  (const std::vector<L> & pseudojets,
62  const JetDefinition & jet_def_in,
63  const GhostedAreaSpec & ghost_spec,
64  const bool & writeout_combinations = false)
66  std::vector<L> * ghosts = NULL;
67  _initialise(pseudojets,jet_def_in,&ghost_spec,ghosts,0.0,
68  writeout_combinations); }
69 
71  (const std::vector<L> & pseudojets,
72  const JetDefinition & jet_def_in,
73  const std::vector<L> & ghosts,
74  double ghost_area,
75  const bool & writeout_combinations = false)
77  const GhostedAreaSpec * ghost_spec = NULL;
78  _initialise(pseudojets,jet_def_in,ghost_spec,&ghosts,ghost_area,
79  writeout_combinations); }
80 
81 
82  /// does the actual work of initialisation
83  template<class L> void _initialise
84  (const std::vector<L> & pseudojets,
85  const JetDefinition & jet_def_in,
86  const GhostedAreaSpec * ghost_spec,
87  const std::vector<L> * ghosts,
88  double ghost_area,
89  const bool & writeout_combinations);
90 
91  //vector<PseudoJet> constituents (const PseudoJet & jet) const;
92 
93  /// returns the number of hard particles (i.e. those supplied by the user).
94  unsigned int n_hard_particles() const;
95 
96  /// returns the area of a jet
97  virtual double area (const PseudoJet & jet) const;
98 
99  /// returns a four vector corresponding to the sum (E-scheme) of the
100  /// ghost four-vectors composing the jet area, normalised such that
101  /// for a small contiguous area the p_t of the extended_area jet is
102  /// equal to area of the jet.
103  virtual PseudoJet area_4vector (const PseudoJet & jet) const;
104 
105  /// true if a jet is made exclusively of ghosts
106  virtual bool is_pure_ghost(const PseudoJet & jet) const;
107 
108  /// true if the entry in the history index corresponds to a
109  /// ghost; if hist_ix does not correspond to an actual particle
110  /// (i.e. hist_ix < 0), then the result is false.
111  bool is_pure_ghost(int history_index) const;
112 
113  /// this class does have explicit ghosts
114  virtual bool has_explicit_ghosts() const {return true;}
115 
116  /// return the total area, corresponding to a given Selector, that
117  /// consists of unclustered ghosts
118  ///
119  /// The selector needs to apply jet by jet
120  virtual double empty_area(const Selector & selector) const;
121 
122  /// returns the total area under study
123  double total_area () const;
124 
125  /// returns the largest squared transverse momentum among
126  /// all ghosts
127  double max_ghost_perp2() const {return _max_ghost_perp2;}
128 
129  /// returns true if there are any particles whose transverse momentum
130  /// if so low that there's a risk of the ghosts having modified the
131  /// clustering sequence
132  bool has_dangerous_particles() const {return _has_dangerous_particles;}
133 
134  /// get the area of the ghosts
135  //double ghost_area() const{return _ghost_area;}
136 
137 private:
138 
139  int _n_ghosts;
140  double _ghost_area;
141  std::vector<bool> _is_pure_ghost;
142  std::vector<double> _areas;
143  std::vector<PseudoJet> _area_4vectors;
144 
145  // things related to checks for dangerous particles
146  double _max_ghost_perp2;
147  bool _has_dangerous_particles;
148  static LimitedWarning _warnings;
149 
150  //static int _n_warn_dangerous_particles;
151  //static const int _max_warn_dangerous_particles = 5;
152 
153 
154  unsigned int _initial_hard_n;
155 
156  /// adds the "ghost" momenta, which will be used to estimate
157  /// the jet area
158  void _add_ghosts(const GhostedAreaSpec & ghost_spec);
159 
160  /// another way of adding ghosts
161  template<class L> void _add_ghosts (
162  const std::vector<L> & ghosts,
163  double ghost_area);
164 
165  /// routine to be called after the processing is done so as to
166  /// establish summary information on all the jets (areas, whether
167  /// pure ghost, etc.)
168  void _post_process();
169 
170 };
171 
172 
173 //----------------------------------------------------------------------
174 // initialise from some generic type... Has to be made available
175 // here in order for the template aspect of it to work...
176 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_initialise
177  (const std::vector<L> & pseudojets,
178  const JetDefinition & jet_def_in,
179  const GhostedAreaSpec * ghost_spec,
180  const std::vector<L> * ghosts,
181  double ghost_area,
182  const bool & writeout_combinations) {
183  // don't reserve space yet -- will be done below
184 
185  // insert initial jets this way so that any type L that can be
186  // converted to a pseudojet will work fine (basically PseudoJet
187  // and any type that has [] subscript access to the momentum
188  // components, such as CLHEP HepLorentzVector).
189  for (unsigned int i = 0; i < pseudojets.size(); i++) {
190  PseudoJet mom(pseudojets[i]);
191  //mom.set_user_index(0); // for user's particles (user index now lost...)
192  _jets.push_back(mom);
193  _is_pure_ghost.push_back(false);
194  }
195 
196  _initial_hard_n = _jets.size();
197 
198  if (ghost_spec != NULL) {
199  //std::cout << "about to reserve " << (_jets.size()+ghost_spec->n_ghosts())*2 << std::endl;
200  _jets.reserve((_jets.size()+ghost_spec->n_ghosts()));
201  _add_ghosts(*ghost_spec);
202  } else {
203  _jets.reserve(_jets.size()+ghosts->size());
204  _add_ghosts(*ghosts, ghost_area);
205  }
206 
207  if (writeout_combinations) {
208  std::cout << "# Printing particles including ghosts\n";
209  for (unsigned j = 0; j < _jets.size(); j++) {
210  printf("%5u %20.13f %20.13f %20.13e\n",
211  j,_jets[j].rap(),_jets[j].phi_02pi(),_jets[j].kt2());
212  }
213  std::cout << "# Finished printing particles including ghosts\n";
214  }
215 
216  // this will ensure that we can still point to jets without
217  // difficulties arising!
218  //std::cout << _jets.size() << " " << _jets.size()*2 << " " << _jets.max_size() << std::endl;
219  _jets.reserve(_jets.size()*2); //GPS tmp removed
220 
221  // run the clustering
222  _initialise_and_run(jet_def_in,writeout_combinations);
223 
224  // set up all other information
225  _post_process();
226 }
227 
228 
229 inline unsigned int ClusterSequenceActiveAreaExplicitGhosts::n_hard_particles() const {return _initial_hard_n;}
230 
231 
232 //----------------------------------------------------------------------
233 /// add an explicitly specified bunch of ghosts
234 template<class L> void ClusterSequenceActiveAreaExplicitGhosts::_add_ghosts (
235  const std::vector<L> & ghosts,
236  double ghost_area) {
237 
238 
239  for (unsigned i = 0; i < ghosts.size(); i++) {
240  _is_pure_ghost.push_back(true);
241  _jets.push_back(ghosts[i]);
242  }
243  // and record some info about ghosts
244  _ghost_area = ghost_area;
245  _n_ghosts = ghosts.size();
246 }
247 
248 
249 FASTJET_END_NAMESPACE
250 
251 #endif // __FASTJET_CLUSTERSEQUENCEACTIVEAREAEXPLICITGHOSTS_HH_
virtual bool has_explicit_ghosts() const
this class does have explicit ghosts
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
double max_ghost_perp2() const
returns the largest squared transverse momentum among all ghosts
base class that sets interface for extensions of ClusterSequence that provide information about the a...
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...
Like ClusterSequence with computation of the active jet area with the addition of explicit ghosts...
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:141
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