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_
fastjet::ClusterSequenceActiveAreaExplicitGhosts::max_ghost_perp2
double max_ghost_perp2() const
returns the largest squared transverse momentum among all ghosts
Definition: ClusterSequenceActiveAreaExplicitGhosts.hh:127
fastjet::ClusterSequenceActiveAreaExplicitGhosts
Definition: ClusterSequenceActiveAreaExplicitGhosts.hh:55
fastjet::JetDefinition
Definition: JetDefinition.hh:173
fastjet::GhostedAreaSpec
Definition: GhostedAreaSpec.hh:62
fastjet::PseudoJet
Definition: PseudoJet.hh:65
fastjet::ClusterSequenceAreaBase
Definition: ClusterSequenceAreaBase.hh:45
fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles
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...
Definition: ClusterSequenceActiveAreaExplicitGhosts.hh:132
fastjet::Selector
Definition: Selector.hh:141
fastjet::LimitedWarning
Definition: LimitedWarning.hh:45
fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_explicit_ghosts
virtual bool has_explicit_ghosts() const
this class does have explicit ghosts
Definition: ClusterSequenceActiveAreaExplicitGhosts.hh:114