Fork me on GitHub

source: git/external/fastjet/contribs/RecursiveTools/Recluster.cc@ 2ad823e

Timing
Last change on this file since 2ad823e was 1f1f858, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago

added RecursiveTools package from fastjet contribs

  • Property mode set to 100644
File size: 15.4 KB
Line 
1// $Id: Recluster.cc 699 2014-07-07 09:58:12Z gsoyez $
2//
3// Copyright (c) 2014-, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
4//
5//----------------------------------------------------------------------
6// This file is part of FastJet contrib.
7//
8// It is free software; you can redistribute it and/or modify it under
9// the terms of the GNU General Public License as published by the
10// Free Software Foundation; either version 2 of the License, or (at
11// your option) any later version.
12//
13// It is distributed in the hope that it will be useful, but WITHOUT
14// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16// License for more details.
17//
18// You should have received a copy of the GNU General Public License
19// along with this code. If not, see <http://www.gnu.org/licenses/>.
20//----------------------------------------------------------------------
21
22#include "Recluster.hh"
23#include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
24#include <sstream>
25#include <typeinfo>
26
27using namespace std;
28
29// Comments:
30//
31// - If the jet comes from a C/A clustering (or is a composite jet
32// made of C/A clusterings) and we recluster it with a C/A
33// algorithm, we just use exclusive jets instead of doing the
34// clustering explicitly. In this specific case, we need to make
35// sure that all the pieces share the same cluster sequence.
36//
37// - If the recombiner has to be taken from the original jet and that
38// jet is composite, we need to check that all the pieces were
39// obtained with the same recombiner.
40//
41// TODO:
42//
43// - check this actually works!!!
44
45FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
46
47namespace contrib{
48
49LimitedWarning Recluster::_explicit_ghost_warning;
50
51// class description
52string Recluster::description() const {
53 ostringstream ostr;
54 ostr << "Recluster with subjet_def = ";
55 if (_use_full_def) {
56 ostr << _subjet_def.description();
57 } else {
58 if (_subjet_alg == kt_algorithm) {
59 ostr << "Longitudinally invariant kt algorithm with R = " << _subjet_radius;
60 } else if (_subjet_alg == cambridge_algorithm) {
61 ostr << "Longitudinally invariant Cambridge/Aachen algorithm with R = " << _subjet_radius;
62 } else if (_subjet_alg == antikt_algorithm) {
63 ostr << "Longitudinally invariant anti-kt algorithm with R = " << _subjet_radius;
64 } else if (_subjet_alg == genkt_algorithm) {
65 ostr << "Longitudinally invariant generalised kt algorithm with R = " << _subjet_radius
66 << ", p = " << _subjet_extra;
67 } else if (_subjet_alg == cambridge_for_passive_algorithm) {
68 ostr << "Longitudinally invariant Cambridge/Aachen algorithm with R = " << _subjet_radius
69 << " and a special hack whereby particles with kt < "
70 << _subjet_extra << "are treated as passive ghosts";
71 } else if (_subjet_alg == ee_kt_algorithm) {
72 ostr << "e+e- kt (Durham) algorithm";
73 } else if (_subjet_alg == ee_genkt_algorithm) {
74 ostr << "e+e- generalised kt algorithm with R = " << _subjet_radius
75 << ", p = " << _subjet_extra;
76 } else if (_subjet_alg == undefined_jet_algorithm) {
77 ostr << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
78 } else {
79 ostr << "unrecognized jet_algorithm";
80 }
81 ostr << ", a recombiner obtained from the jet being reclustered";
82 }
83
84 if (_single)
85 ostr << " and keeping the hardest subjet";
86 else
87 ostr << " and joining all subjets in a composite jet";
88
89 return ostr.str();
90}
91
92// return a vector of subjets, which are the ones that would be kept
93// by the filtering
94PseudoJet Recluster::result(const PseudoJet &jet) const {
95 // generic sanity checks
96 //-------------------------------------------------------------------
97 // make sure that the jet has constituents
98 if (! jet.has_constituents())
99 throw Error("Filter can only be applied on jets having constituents");
100
101 // tests particular to certain configurations
102 //-------------------------------------------------------------------
103
104 // for a whole variety of checks, we shall need the "recursive"
105 // pieces of the jet (the jet itself or recursing down to its most
106 // fundamental pieces). So we do compute these once and for all.
107 //
108 // Note that the pieces are always needed (either for C/A or for the
109 // area checks)
110 vector<PseudoJet> all_pieces; //.clear();
111 if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0)){
112 throw Error("Recluster: failed to retrieve all the pieces composing the jet.");
113 }
114
115 // decide which jet definition to use
116 //-------------------------------------------------------------------
117 JetDefinition subjet_def;
118 if (_use_full_def){
119 subjet_def = _subjet_def;
120 } else {
121 _build_jet_def_with_recombiner(all_pieces, subjet_def);
122 }
123
124
125 // the vector that will ultimately hold the subjets
126 vector<PseudoJet> subjets;
127
128 // check if we can apply the simplification for C/A jets reclustered
129 // with C/A
130 //
131 // we apply C/A clustering iff
132 // - the request subjet_def is C/A
133 // - the jet is either directly coming from C/A or if it is a
134 // superposition of C/A jets from the same cluster sequence
135 // - the pieces agree with the recombination scheme of subjet_def
136 //
137 // Note that in this case area support will be automatically
138 // inherted so we can only worry about this later
139 //-------------------------------------------------------------------
140 if (_check_ca(all_pieces, subjet_def)){
141 _recluster_cafilt(all_pieces, subjets, subjet_def.R());
142 subjets = sorted_by_pt(subjets);
143 return _single
144 ? subjets[0]
145 : join(subjets, *(subjet_def.recombiner()));
146 }
147
148 // decide if area support has to be kept
149 //-------------------------------------------------------------------
150 bool include_area_support = jet.has_area();
151 if ((include_area_support) && (!_check_explicit_ghosts(all_pieces))){
152 _explicit_ghost_warning.warn("Recluster: the original cluster sequence is lacking explicit ghosts; area support will no longer be available after re-clustering");
153 include_area_support = false;
154 }
155
156 // extract the subjets
157 //-------------------------------------------------------------------
158 _recluster_generic(jet, subjets, subjet_def, include_area_support);
159 subjets = sorted_by_pt(subjets);
160
161 return _single
162 ? subjets[0]
163 : join(subjets, *(subjet_def.recombiner()));
164}
165
166//----------------------------------------------------------------------
167// the parts that really do the reclustering
168//----------------------------------------------------------------------
169
170// get the subjets in the simple case of C/A+C/A
171void Recluster::_recluster_cafilt(const vector<PseudoJet> & all_pieces,
172 vector<PseudoJet> & subjets,
173 double Rfilt) const{
174 subjets.clear();
175
176 // each individual piece should have a C/A cluster sequence
177 // associated to it. So a simple loop would do the job
178 for (vector<PseudoJet>::const_iterator piece_it = all_pieces.begin();
179 piece_it!=all_pieces.end(); piece_it++){
180 // just extract the exclusive subjets of 'jet'
181 const ClusterSequence *cs = piece_it->associated_cluster_sequence();
182 vector<PseudoJet> local_subjets;
183
184 double dcut = Rfilt / cs->jet_def().R();
185 if (dcut>=1.0){
186 local_subjets.push_back(*piece_it);
187 } else {
188 local_subjets = piece_it->exclusive_subjets(dcut*dcut);
189 }
190
191 copy(local_subjets.begin(), local_subjets.end(), back_inserter(subjets));
192 }
193}
194
195
196// set the filtered elements in the generic re-clustering case (w/o
197// subtraction)
198void Recluster::_recluster_generic(const PseudoJet & jet,
199 vector<PseudoJet> & subjets,
200 const JetDefinition & subjet_def,
201 bool do_areas) const{
202 // create a new, internal, ClusterSequence from the jet constituents
203 // get the subjets directly from there
204 //
205 // If the jet has area support then we separate the ghosts from the
206 // "regular" particles so the subjets will also have area
207 // support. Note that we do this regardless of whether rho is zero
208 // or not.
209 //
210 // Note that to be able to separate the ghosts, one needs explicit
211 // ghosts!!
212 // ---------------------------------------------------------------
213 if (do_areas){
214 vector<PseudoJet> all_constituents = jet.constituents();
215 vector<PseudoJet> regular_constituents, ghosts;
216
217 for (vector<PseudoJet>::iterator it = all_constituents.begin();
218 it != all_constituents.end(); it++){
219 if (it->is_pure_ghost())
220 ghosts.push_back(*it);
221 else
222 regular_constituents.push_back(*it);
223 }
224
225 // figure out the ghost area from the 1st ghost (if none, any value
226 // would probably do as the area will be 0 and subtraction will have
227 // no effect!)
228 double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
229 ClusterSequenceActiveAreaExplicitGhosts * csa
230 = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents,
231 subjet_def,
232 ghosts, ghost_area);
233
234 subjets = csa->inclusive_jets();
235
236 // allow the cs to be deleted when it's no longer used
237 //
238 // Note that there is at least one constituent in the jet so there
239 // is in principle at least one subjet But one may have used a
240 // nasty recombiner that left an empty set of subjets, so we'd
241 // rather play it safe
242 if (subjets.size())
243 csa->delete_self_when_unused();
244 else
245 delete csa;
246 } else {
247 ClusterSequence * cs = new ClusterSequence(jet.constituents(), subjet_def);
248 subjets = cs->inclusive_jets();
249 // allow the cs to be deleted when it's no longer used (again, we
250 // add an extra safety check)
251 if (subjets.size())
252 cs->delete_self_when_unused();
253 else
254 delete cs;
255 }
256}
257
258
259//----------------------------------------------------------------------
260// various checks and internal constructs
261//----------------------------------------------------------------------
262
263// fundamental info for CompositeJets
264//----------------------------------------------------------------------
265
266// get the pieces down to the fundamental pieces
267//
268// Note that this just checks that there is an associated CS to the
269// fundamental pieces, not that it is still valid
270bool Recluster::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
271 if (jet.has_associated_cluster_sequence()){
272 all_pieces.push_back(jet);
273 return true;
274 }
275
276 if (jet.has_pieces()){
277 const vector<PseudoJet> pieces = jet.pieces();
278 for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
279 if (!_get_all_pieces(*it, all_pieces)) return false;
280 return true;
281 }
282
283 return false;
284}
285
286// treatment of recombiners
287//----------------------------------------------------------------------
288// get the common recombiner to all pieces (NULL if none)
289//
290// Note that if the jet has an associated cluster sequence that is no
291// longer valid, an error will be thrown (needed since it could be the
292// 1st check called after the enumeration of the pieces)
293const JetDefinition::Recombiner* Recluster::_get_common_recombiner(const vector<PseudoJet> &all_pieces) const{
294 const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
295 for (unsigned int i=1; i<all_pieces.size(); i++)
296 if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)) return NULL;
297
298 return jd_ref.recombiner();
299}
300
301void Recluster::_build_jet_def_with_recombiner(const vector<PseudoJet> &all_pieces,
302 JetDefinition &subjet_def) const{
303 // the recombiner has to be guessed from the pieces
304 const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
305 if (common_recombiner) {
306 if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
307 RecombinationScheme scheme =
308 static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
309 if (_has_subjet_extra)
310 subjet_def = JetDefinition(_subjet_alg, _subjet_radius, _subjet_extra, scheme);
311 else if (_has_subjet_radius)
312 subjet_def = JetDefinition(_subjet_alg, _subjet_radius, scheme);
313 else
314 subjet_def = JetDefinition(_subjet_alg, scheme);
315 } else {
316 if (_has_subjet_extra)
317 subjet_def = JetDefinition(_subjet_alg, _subjet_radius, _subjet_extra, common_recombiner);
318 else if (_has_subjet_radius)
319 subjet_def = JetDefinition(_subjet_alg, _subjet_radius, common_recombiner);
320 else
321 subjet_def = JetDefinition(_subjet_alg, common_recombiner);
322 }
323 } else {
324 throw Error("Recluster: requested to guess the recombination scheme (or recombiner) from the original jet but an inconsistency was found between the pieces constituing that jet.");
325 }
326}
327
328// area support
329//----------------------------------------------------------------------
330
331// check if the jet (or all its pieces) have explicit ghosts
332// (assuming the jet has area support).
333//
334// Note that if the jet has an associated cluster sequence that is no
335// longer valid, an error will be thrown (needed since it could be the
336// 1st check called after the enumeration of the pieces)
337bool Recluster::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
338 for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
339 if (! it->validated_csab()->has_explicit_ghosts()) return false;
340 return true;
341}
342
343// C/A specific tests
344//----------------------------------------------------------------------
345
346// check if one can apply the simplification for C/A subjets
347//
348// This includes:
349// - the subjet definition asks for C/A subjets
350// - all the pieces share the same CS
351// - that CS is C/A with the same recombiner as the subjet def
352// - the filtering radius is not larger than any of the pairwise
353// distance between the pieces
354//
355// Note that if the jet has an associated cluster sequence that is no
356// longer valid, an error will be thrown (needed since it could be the
357// 1st check called after the enumeration of the pieces)
358bool Recluster::_check_ca(const vector<PseudoJet> &all_pieces,
359 const JetDefinition &subjet_def) const{
360 if (subjet_def.jet_algorithm() != cambridge_algorithm) return false;
361
362 // check that the 1st of all the pieces (we're sure there is at
363 // least one) is coming from a C/A clustering. Then check that all
364 // the following pieces share the same ClusterSequence
365 const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
366 if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
367 for (unsigned int i=1; i<all_pieces.size(); i++)
368 if (all_pieces[i].validated_cs() != cs_ref) return false;
369
370 // check that the 1st peice has the same recombiner as the one used
371 // for the subjet clustering
372 // Note that since they share the same CS, checking the 1st one is enough
373 if (!cs_ref->jet_def().has_same_recombiner(subjet_def)) return false;
374
375 // we also have to make sure that the reclustering radius is not larger
376 // than any of the inter-pieces distance
377 double Rsub2 = subjet_def.R();
378 Rsub2 *= Rsub2;
379 for (unsigned int i=0; i<all_pieces.size()-1; i++){
380 for (unsigned int j=i+1; j<all_pieces.size(); j++){
381 if (all_pieces[i].squared_distance(all_pieces[j]) < Rsub2) return false;
382 }
383 }
384
385 return true;
386}
387
388} // contrib namespace
389
390FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.