Fork me on GitHub

source: git/external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.cc@ 82e07cf

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 82e07cf was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 7.2 KB
Line 
1//STARTHEADER
2// $Id$
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
32using namespace std;
33
34FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
35
36// save some typing
37typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
38
39
40LimitedWarning ClustSeqActAreaEG::_warnings;
41
42//----------------------------------------------------------------------
43///
44void 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
64double ClustSeqActAreaEG::area (const PseudoJet & jet) const {
65 return _areas[jet.cluster_hist_index()];
66}
67
68
69//----------------------------------------------------------------------
70// return the total area
71double ClustSeqActAreaEG::total_area () const {
72 return _n_ghosts * _ghost_area;
73}
74
75
76//----------------------------------------------------------------------
77// return the extended area of a jet
78PseudoJet ClustSeqActAreaEG::area_4vector (const PseudoJet & jet) const {
79 return _area_4vectors[jet.cluster_hist_index()];
80}
81
82//----------------------------------------------------------------------
83bool ClustSeqActAreaEG::is_pure_ghost(const PseudoJet & jet) const
84{
85 return _is_pure_ghost[jet.cluster_hist_index()];
86}
87
88//----------------------------------------------------------------------
89bool ClustSeqActAreaEG::is_pure_ghost(int hist_ix) const
90{
91 return hist_ix >= 0 ? _is_pure_ghost[hist_ix] : false;
92}
93
94//----------------------------------------------------------------------
95double 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
113void 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
201FASTJET_END_NAMESPACE
202
Note: See TracBrowser for help on using the repository browser.