Fork me on GitHub

source: svn/trunk/external/fastjet/tools/Pruner.cc@ 1397

Last change on this file since 1397 was 999, checked in by Pavel Demin, 12 years ago

add fastjet/tools

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 13.1 KB
Line 
1//STARTHEADER
2// $Id: Pruner.cc 999 2013-03-04 11:48:06Z pavel $
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/tools/Pruner.hh"
30#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
31#include "fastjet/Selector.hh"
32#include <cassert>
33#include <algorithm>
34#include <sstream>
35#include <typeinfo>
36
37using namespace std;
38
39
40FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41
42
43//----------------------------------------------------------------------
44// class Pruner
45//----------------------------------------------------------------------
46
47//----------------------------------------------------------------------
48// alternative (dynamic) ctor
49// \param jet_def the jet definition for the internal clustering
50// \param zcut_dyn dynamic pt-fraction cut in the pruning
51// \param Rcut_dyn dynamic angular distance cut in the pruning
52Pruner::Pruner(const JetDefinition &jet_def,
53 FunctionOfPseudoJet<double> *zcut_dyn,
54 FunctionOfPseudoJet<double> *Rcut_dyn)
55 : _jet_def(jet_def), _zcut(0), _Rcut_factor(0),
56 _zcut_dyn(zcut_dyn), _Rcut_dyn(Rcut_dyn), _get_recombiner_from_jet(false) {
57 assert(_zcut_dyn != 0 && _Rcut_dyn != 0);
58}
59
60//----------------------------------------------------------------------
61// action on a single jet
62PseudoJet Pruner::result(const PseudoJet &jet) const{
63 // pruning can only be applied to jets that have constituents
64 if (!jet.has_constituents()){
65 throw Error("Pruner: trying to apply the Pruner transformer to a jet that has no constituents");
66 }
67
68 // if the jet has area support and there are explicit ghosts, we can
69 // transfer that support to the internal re-clustering
70 bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
71
72 // build the pruning plugin
73 double Rcut = (_Rcut_dyn) ? (*_Rcut_dyn)(jet) : _Rcut_factor * 2.0*jet.m()/jet.perp();
74 double zcut = (_zcut_dyn) ? (*_zcut_dyn)(jet) : _zcut;
75 PruningPlugin * pruning_plugin;
76 // for some constructors, we get the recombiner from the
77 // input jet -- some acrobatics are needed (see plans for FJ3.1
78 // for a hopefully better solution).
79 if (_get_recombiner_from_jet) {
80 const JetDefinition::Recombiner * common_recombiner =
81 _get_common_recombiner(jet);
82 if (common_recombiner) {
83 JetDefinition jet_def = _jet_def;
84 if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
85 RecombinationScheme scheme =
86 static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
87 jet_def.set_recombination_scheme(scheme);
88 } else {
89 jet_def.set_recombiner(common_recombiner);
90 }
91 pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
92 } else {
93 // if there wasn't a common recombiner, we just use the default
94 // recombiner that was in _jet_def
95 pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
96 }
97 } else {
98 pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
99 }
100
101 // now recluster the constituents of the jet with that plugin
102 JetDefinition internal_jet_def(pruning_plugin);
103 // flag the plugin for automatic deletion _before_ we make
104 // copies (so that as long as the copies are also present
105 // it doesn't get deleted).
106 internal_jet_def.delete_plugin_when_unused();
107
108 ClusterSequence * cs;
109 if (do_areas){
110 vector<PseudoJet> particles, ghosts;
111 SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
112 // determine the ghost area from the 1st ghost (if none, any value
113 // will do, as the area will be 0 and subtraction will have
114 // no effect!)
115 double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
116 cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def,
117 ghosts, ghost_area);
118 } else {
119 cs = new ClusterSequence(jet.constituents(), internal_jet_def);
120 }
121
122 PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
123 PrunerStructure * s = new PrunerStructure(result_local);
124 result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
125
126 // make sure things remain persistent -- i.e. tell the jet definition
127 // and the cluster sequence that it is their responsibility to clean
128 // up memory once the "result" reaches the end of its life in the user's
129 // code. (The CS deletes itself when the result goes out of scope and
130 // that also triggers deletion of the plugin)
131 cs->delete_self_when_unused();
132
133 return result_local;
134}
135
136// check if the jet has explicit_ghosts (knowing that there is an
137// area support)
138bool Pruner::_check_explicit_ghosts(const PseudoJet &jet) const{
139 // if the jet comes from a Clustering check explicit ghosts in that
140 // clustering
141 if (jet.has_associated_cluster_sequence())
142 return jet.validated_csab()->has_explicit_ghosts();
143
144 // if the jet has pieces, recurse in the pieces
145 if (jet.has_pieces()){
146 vector<PseudoJet> pieces = jet.pieces();
147 for (unsigned int i=0;i<pieces.size(); i++)
148 if (!_check_explicit_ghosts(pieces[i])) return false;
149 // never returned false, so we're OK.
150 return true;
151 }
152
153 // return false for any other (unknown) structure
154 return false;
155}
156
157// see if there is a common recombiner among the pieces; if there
158// is return a pointer to it; otherwise, return NULL.
159//
160// NB: this way of doing things is not ideal, because quite some work
161// is needed to get a correct handling of the final recombiner
162// (e.g. default v. non-default). In future add
163// set_recombiner(jet_def) to JetDefinition, maybe also add
164// an invalid_scheme to the default recombiner and then
165// do all the work below directly with a JetDefinition directly
166// together with JD::has_same_recombiner(...)
167const JetDefinition::Recombiner * Pruner::_get_common_recombiner(const PseudoJet &jet) const{
168 if (jet.has_associated_cluster_sequence())
169 return jet.validated_cs()->jet_def().recombiner();
170
171 // if the jet has pieces, recurse in the pieces
172 if (jet.has_pieces()){
173 vector<PseudoJet> pieces = jet.pieces();
174 if (pieces.size() == 0) return 0;
175 const JetDefinition::Recombiner * reco = _get_common_recombiner(pieces[0]);
176 for (unsigned int i=1;i<pieces.size(); i++)
177 if (_get_common_recombiner(pieces[i]) != reco) return 0;
178 // never returned false, so we're OK.
179 return reco;
180 }
181
182 // return false for any other (unknown) structure
183 return 0;
184}
185
186
187// transformer description
188std::string Pruner::description() const{
189 ostringstream oss;
190 oss << "Pruner with jet_definition = (" << _jet_def.description() << ")";
191 if (_zcut_dyn) {
192 oss << ", dynamic zcut (" << _zcut_dyn->description() << ")"
193 << ", dynamic Rcut (" << _Rcut_dyn->description() << ")";
194 } else {
195 oss << ", zcut = " << _zcut
196 << ", Rcut_factor = " << _Rcut_factor;
197 }
198 return oss.str();
199}
200
201
202
203//----------------------------------------------------------------------
204// class PrunerStructure
205//----------------------------------------------------------------------
206
207// return the other jets that may have been found along with the
208// result of the pruning
209// The resulting vector is sorted in pt
210vector<PseudoJet> PrunerStructure::extra_jets() const{
211 return sorted_by_pt((!SelectorNHardest(1))(validated_cs()->inclusive_jets()));;
212}
213
214
215//----------------------------------------------------------------------
216// class PruningRecombiner
217//----------------------------------------------------------------------
218
219// decide whether to recombine things or not
220void PruningRecombiner::recombine(const PseudoJet &pa,
221 const PseudoJet &pb,
222 PseudoJet &pab) const{
223 PseudoJet p;
224 _recombiner->recombine(pa, pb, p);
225
226 // if the 2 particles are close enough, do the recombination
227 if (pa.squared_distance(pb)<=_Rcut2){
228 pab=p; return;
229 }
230
231 double pt2a = pa.perp2();
232 double pt2b = pb.perp2();
233
234 // check which is the softest
235 if (pt2a < pt2b){
236 if (pt2a<_zcut2*p.perp2()){
237 pab = pb; _rejected.push_back(pa.cluster_hist_index());
238 } else {
239 pab = p;
240 }
241 } else {
242 if (pt2b<_zcut2*p.perp2()) {
243 pab = pa; _rejected.push_back(pb.cluster_hist_index());
244 } else {
245 pab = p;
246 }
247 }
248}
249
250// description
251string PruningRecombiner::description() const{
252 ostringstream oss;
253 oss << "Pruning recombiner with zcut = " << sqrt(_zcut2)
254 << ", Rcut = " << sqrt(_Rcut2)
255 << ", and underlying recombiner = " << _recombiner->description();
256 return oss.str();
257}
258
259
260
261
262//----------------------------------------------------------------------
263// class PruningPlugin
264//----------------------------------------------------------------------
265// the actual clustering work for the plugin
266void PruningPlugin::run_clustering(ClusterSequence &input_cs) const{
267 // declare a pruning recombiner
268 PruningRecombiner pruning_recombiner(_zcut, _Rcut, _jet_def.recombiner());
269 JetDefinition jet_def = _jet_def;
270 jet_def.set_recombiner(&pruning_recombiner);
271
272 // cluster the particles using that recombiner
273 ClusterSequence internal_cs(input_cs.jets(), jet_def);
274 const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
275
276 // transfer the list of "childless" elements into a bool vector
277 vector<bool> kept(internal_hist.size(), true);
278 const vector<unsigned int> &pr_rej = pruning_recombiner.rejected();
279 for (unsigned int i=0;i<pr_rej.size(); i++) kept[pr_rej[i]]=false;
280
281 // browse the history, keeping only the elements that have not been
282 // vetoed.
283 //
284 // In the process we build a map for the history indices
285 vector<unsigned int> internal2input(internal_hist.size());
286 for (unsigned int i=0; i<input_cs.jets().size(); i++)
287 internal2input[i] = i;
288
289 for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
290 const ClusterSequence::history_element &he = internal_hist[i];
291
292 // deal with recombinations with the beam
293 if (he.parent2 == ClusterSequence::BeamJet){
294 int internal_jetp_index = internal_hist[he.parent1].jetp_index;
295 int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
296
297 int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
298
299 // cout << "Beam recomb for internal " << internal_hist_index
300 // << " (input jet index=" << input_jetp_index << endl;
301
302 input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
303 continue;
304 }
305
306 // now, deal with two-body recombinations
307 if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
308 internal2input[i]=internal2input[he.parent2];
309 // cout << "rejecting internal " << he.parent1
310 // << ", mapping internal " << i
311 // << " to internal " << he.parent2
312 // << " i.e. " << internal2input[i] << endl;
313 } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
314 internal2input[i]=internal2input[he.parent1];
315 // cout << "rejecting internal " << he.parent2
316 // << ", mapping internal " << i
317 // << " to internal " << he.parent1
318 // << " i.e. " << internal2input[i] << endl;
319 } else { // do the recombination
320 int new_index;
321 input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
322 input_cs.history()[internal2input[he.parent2]].jetp_index,
323 he.dij, internal_cs.jets()[he.jetp_index], new_index);
324 internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
325 // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
326 // << " and " << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
327 // << " into " << internal2input[i] << " (int: " << i << ")" << endl;
328 }
329 }
330}
331
332// returns the plugin description
333string PruningPlugin::description() const{
334 ostringstream oss;
335 oss << "Pruning plugin with jet_definition = (" << _jet_def.description()
336 <<"), zcut = " << _zcut
337 << ", Rcut = " << _Rcut;
338 return oss.str();
339}
340
341
342FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.