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