FastJet  3.3.3
ClusterSequenceActiveAreaExplicitGhosts.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequenceActiveAreaExplicitGhosts.cc 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 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
32 #include<limits>
33 
34 using namespace std;
35 
36 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37 
38 // save some typing
39 typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
40 
41 
42 LimitedWarning ClustSeqActAreaEG::_warnings;
43 
44 //----------------------------------------------------------------------
45 ///
46 void ClustSeqActAreaEG::_add_ghosts (
47  const GhostedAreaSpec & ghost_spec) {
48 
49  // add the ghosts to the jets
50  ghost_spec.add_ghosts(_jets);
51 
52  // now add labelling...
53  for (unsigned i = _initial_hard_n; i < _jets.size(); i++) {
54  //_jets[i].set_user_index(1);
55  _is_pure_ghost.push_back(true);
56  }
57 
58  // and record some info from the ghost_spec
59  _ghost_area = ghost_spec.actual_ghost_area();
60  _n_ghosts = ghost_spec.n_ghosts();
61 }
62 
63 
64 //----------------------------------------------------------------------
65 // return the area of a jet
66 double ClustSeqActAreaEG::area (const PseudoJet & jet) const {
67  return _areas[jet.cluster_hist_index()];
68 }
69 
70 
71 //----------------------------------------------------------------------
72 // return the total area
73 double ClustSeqActAreaEG::total_area () const {
74  return _n_ghosts * _ghost_area;
75 }
76 
77 
78 //----------------------------------------------------------------------
79 // return the extended area of a jet
80 PseudoJet ClustSeqActAreaEG::area_4vector (const PseudoJet & jet) const {
81  return _area_4vectors[jet.cluster_hist_index()];
82 }
83 
84 //----------------------------------------------------------------------
85 bool ClustSeqActAreaEG::is_pure_ghost(const PseudoJet & jet) const
86 {
87  return _is_pure_ghost[jet.cluster_hist_index()];
88 }
89 
90 //----------------------------------------------------------------------
91 bool ClustSeqActAreaEG::is_pure_ghost(int hist_ix) const
92 {
93  return hist_ix >= 0 ? _is_pure_ghost[hist_ix] : false;
94 }
95 
96 //----------------------------------------------------------------------
97 double ClustSeqActAreaEG::empty_area(const Selector & selector) const {
98  // make sure that the selector applies jet by jet
99  if (! selector.applies_jet_by_jet()){
100  throw Error("ClusterSequenceActiveAreaExplicitGhosts: empty area can only be computed from selectors applying jet by jet");
101  }
102 
103  vector<PseudoJet> unclust = unclustered_particles();
104  double area_local = 0.0;
105  for (unsigned iu = 0; iu < unclust.size(); iu++) {
106  if (is_pure_ghost(unclust[iu]) && selector.pass(unclust[iu])) {
107  area_local += _ghost_area;
108  }
109  }
110  return area_local;
111 }
112 
113 //======================================================================
114 // sort out the areas
115 void ClustSeqActAreaEG::_post_process() {
116 
117  // first check for danger signals.
118  // Establish largest ghost transverse momentum
119  _max_ghost_perp2 = 0.0;
120  for (int i = 0; i < _initial_n; i++) {
121  if (_is_pure_ghost[i] && _jets[i].perp2() > _max_ghost_perp2)
122  _max_ghost_perp2 = _jets[i].perp2();
123  }
124 
125  // now find out if any of the particles are close to danger
126  double danger_ratio = numeric_limits<double>::epsilon();
127  danger_ratio = danger_ratio * danger_ratio;
128  _has_dangerous_particles = false;
129  for (int i = 0; i < _initial_n; i++) {
130  if (!_is_pure_ghost[i] &&
131  danger_ratio * _jets[i].perp2() <= _max_ghost_perp2) {
132  _has_dangerous_particles = true;
133  break;
134  }
135  }
136 
137  if (_has_dangerous_particles) _warnings.warn("ClusterSequenceActiveAreaExplicitGhosts: \n ghosts not sufficiently soft wrt some of the input particles\n a common cause is (unphysical?) input particles with pt=0 but finite rapidity");
138 
139  // sort out sizes
140  _areas.resize(_history.size());
141  _area_4vectors.resize(_history.size());
142  _is_pure_ghost.resize(_history.size());
143 
144 // copy(_jets.begin(), _jets.begin()+_initial_n, _area_4vectors.begin());
145 // for (int i = 0; i < _initial_n; i++) {
146 // if (_is_pure_ghost[i]) {
147 // _areas[i] = _ghost_area;
148 // // normalise pt to be _ghost_area (NB we make use of fact that
149 // // for initial particles, jet and clust_hist index are the same).
150 // _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
151 // } else {
152 // _areas[i] = 0;
153 // _area_4vectors[i].reset(0,0,0,0);
154 // }
155 // }
156 
157  // First set up areas for the initial particles (ghost=_ghost_area,
158  // real particles = 0); recall that _initial_n here is the number of
159  // particles including ghosts
160  for (int i = 0; i < _initial_n; i++) {
161  if (_is_pure_ghost[i]) {
162  _areas[i] = _ghost_area;
163  // normalise pt to be _ghost_area (NB we make use of fact that
164  // for initial particles, jet and clust_hist index are the same).
165  //_area_4vectors[i] = (_ghost_area/_jets[i].perp()) * _jets[i];
166 
167  // NB: we use reset_momentum here, to ensure that the area 4
168  // vectors do not acquire any structure (the structure would not
169  // be meaningful for an area, and it messes up the use count (->
170  // memory leaks if the user call delete_self_when_unused).
171  _area_4vectors[i].reset_momentum(_jets[i]);
172  _area_4vectors[i] *= (_ghost_area/_jets[i].perp());
173  } else {
174  _areas[i] = 0;
175  _area_4vectors[i] = PseudoJet(0.0,0.0,0.0,0.0);
176  }
177  }
178 
179  // next follow the branching through and set up the areas
180  // and ghost-nature at each step of the clustering (rather than
181  // each jet).
182  for (unsigned i = _initial_n; i < _history.size(); i++) {
183  if (_history[i].parent2 == BeamJet) {
184  _is_pure_ghost[i] = _is_pure_ghost[_history[i].parent1];
185  _areas[i] = _areas[_history[i].parent1];
186  _area_4vectors[i] = _area_4vectors[_history[i].parent1];
187  } else {
188  _is_pure_ghost[i] = _is_pure_ghost[_history[i].parent1] &&
189  _is_pure_ghost[_history[i].parent2] ;
190  _areas[i] = _areas[_history[i].parent1] +
191  _areas[_history[i].parent2] ;
192  _jet_def.recombiner()->recombine(_area_4vectors[_history[i].parent1],
193  _area_4vectors[_history[i].parent2],
194  _area_4vectors[i]);
195 // _area_4vectors[i] = _area_4vectors[_history[i].parent1] +
196 // _area_4vectors[_history[i].parent2] ;
197  }
198 
199  }
200 
201 }
202 
203 FASTJET_END_NAMESPACE
204 
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
Definition: Selector.hh:215
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Definition: Selector.hh:178
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition: PseudoJet.hh:784
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67