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