Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/src/ClusterSequenceActiveAreaExplicitGhosts.cc@ 138

Last change on this file since 138 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

File size: 6.2 KB
Line 
1//STARTHEADER
2// $Id: ClusterSequenceActiveAreaExplicitGhosts.cc,v 1.1 2008-11-06 14:32:14 ovyn Exp $
3//
4// Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
26// Foundation, Inc.:
27// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28//----------------------------------------------------------------------
29//ENDHEADER
30
31#include "../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 & area_spec) {
48
49 // add the ghosts to the jets
50 area_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 area_spec
59 _ghost_area = area_spec.actual_ghost_area();
60 _n_ghosts = area_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 RangeDefinition & range) const {
98 vector<PseudoJet> unclust = unclustered_particles();
99 double area = 0.0;
100 for (unsigned iu = 0; iu < unclust.size(); iu++) {
101 if (is_pure_ghost(unclust[iu]) && range.is_in_range(unclust[iu])) {
102 area += _ghost_area;
103 }
104 }
105 return area;
106}
107
108//======================================================================
109// sort out the areas
110void ClustSeqActAreaEG::_post_process() {
111
112 // first check for danger signals.
113 // Establish largest ghost transverse momentum
114 _max_ghost_perp2 = 0.0;
115 for (int i = 0; i < _initial_n; i++) {
116 if (_is_pure_ghost[i] && _jets[i].perp2() > _max_ghost_perp2)
117 _max_ghost_perp2 = _jets[i].perp2();
118 }
119
120 // now find out if any of the particles are close to danger
121 double danger_ratio = numeric_limits<double>::epsilon();
122 danger_ratio = danger_ratio * danger_ratio;
123 _has_dangerous_particles = false;
124 for (int i = 0; i < _initial_n; i++) {
125 if (!_is_pure_ghost[i] &&
126 danger_ratio * _jets[i].perp2() <= _max_ghost_perp2) {
127 _has_dangerous_particles = true;
128 break;
129 }
130 }
131
132 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");
133
134 // sort out sizes
135 _areas.resize(_history.size());
136 _area_4vectors.resize(_history.size());
137 _is_pure_ghost.resize(_history.size());
138
139 // First set up areas for the initial particles (ghost=_ghost_area,
140 // real particles = 0); recall that _initial_n here is the number of
141 // particles including ghosts
142 for (int i = 0; i < _initial_n; i++) {
143 if (_is_pure_ghost[i]) {
144 _areas[i] = _ghost_area;
145 // normalise pt to be _ghost_area (NB we make use of fact that
146 // for initial particles, jet and clust_hist index are the same).
147 _area_4vectors[i] = (_ghost_area/_jets[i].perp()) * _jets[i];
148 } else {
149 _areas[i] = 0;
150 _area_4vectors[i] = PseudoJet(0.0,0.0,0.0,0.0);
151 }
152 }
153
154 // next follow the branching through and set up the areas
155 // and ghost-nature at each step of the clustering (rather than
156 // each jet).
157 for (unsigned i = _initial_n; i < _history.size(); i++) {
158 if (_history[i].parent2 == BeamJet) {
159 _is_pure_ghost[i] = _is_pure_ghost[_history[i].parent1];
160 _areas[i] = _areas[_history[i].parent1];
161 _area_4vectors[i] = _area_4vectors[_history[i].parent1];
162 } else {
163 _is_pure_ghost[i] = _is_pure_ghost[_history[i].parent1] &&
164 _is_pure_ghost[_history[i].parent2] ;
165 _areas[i] = _areas[_history[i].parent1] +
166 _areas[_history[i].parent2] ;
167 _jet_def.recombiner()->recombine(_area_4vectors[_history[i].parent1],
168 _area_4vectors[_history[i].parent2],
169 _area_4vectors[i]);
170// _area_4vectors[i] = _area_4vectors[_history[i].parent1] +
171// _area_4vectors[_history[i].parent2] ;
172 }
173
174 }
175
176}
177
178FASTJET_END_NAMESPACE
179
Note: See TracBrowser for help on using the repository browser.