Fork me on GitHub

source: git/external/fastjet/ClusterSequenceActiveAreaExplicitGhosts.cc@ 601a13b

3.4.3pre09
Last change on this file since 601a13b was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 7.4 KB
Line 
1//FJSTARTHEADER
2// $Id: ClusterSequenceActiveAreaExplicitGhosts.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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
34using namespace std;
35
36FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37
38// save some typing
39typedef ClusterSequenceActiveAreaExplicitGhosts ClustSeqActAreaEG;
40
41
42LimitedWarning ClustSeqActAreaEG::_warnings;
43
44//----------------------------------------------------------------------
45///
46void 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
66double ClustSeqActAreaEG::area (const PseudoJet & jet) const {
67 return _areas[jet.cluster_hist_index()];
68}
69
70
71//----------------------------------------------------------------------
72// return the total area
73double ClustSeqActAreaEG::total_area () const {
74 return _n_ghosts * _ghost_area;
75}
76
77
78//----------------------------------------------------------------------
79// return the extended area of a jet
80PseudoJet ClustSeqActAreaEG::area_4vector (const PseudoJet & jet) const {
81 return _area_4vectors[jet.cluster_hist_index()];
82}
83
84//----------------------------------------------------------------------
85bool ClustSeqActAreaEG::is_pure_ghost(const PseudoJet & jet) const
86{
87 return _is_pure_ghost[jet.cluster_hist_index()];
88}
89
90//----------------------------------------------------------------------
91bool ClustSeqActAreaEG::is_pure_ghost(int hist_ix) const
92{
93 return hist_ix >= 0 ? _is_pure_ghost[hist_ix] : false;
94}
95
96//----------------------------------------------------------------------
97double 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
115void 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
203FASTJET_END_NAMESPACE
204
Note: See TracBrowser for help on using the repository browser.