FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ClusterSequenceActiveAreaExplicitGhosts.hh
1 //FJSTARTHEADER
2 // $Id: ClusterSequenceActiveAreaExplicitGhosts.hh 3433 2014-07-23 08:17:03Z 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_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;
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;
106 
107  /// true if a jet is made exclusively of ghosts
108  virtual bool is_pure_ghost(const PseudoJet & jet) const;
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 {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;
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_