1 | //STARTHEADER
|
---|
2 | // $Id: Filter.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/Filter.hh"
|
---|
30 | #include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
|
---|
31 | #include <cassert>
|
---|
32 | #include <algorithm>
|
---|
33 | #include <sstream>
|
---|
34 | #include <typeinfo>
|
---|
35 |
|
---|
36 | using namespace std;
|
---|
37 |
|
---|
38 |
|
---|
39 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
---|
40 |
|
---|
41 | //----------------------------------------------------------------------
|
---|
42 | // Filter class implementation
|
---|
43 | //----------------------------------------------------------------------
|
---|
44 |
|
---|
45 | // class description
|
---|
46 | string Filter::description() const {
|
---|
47 | ostringstream ostr;
|
---|
48 | ostr << "Filter with subjet_def = ";
|
---|
49 | if (_Rfiltfunc) {
|
---|
50 | ostr << "Cambridge/Aachen algorithm with dynamic Rfilt"
|
---|
51 | << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
|
---|
52 | } else if (_Rfilt > 0) {
|
---|
53 | ostr << "Cambridge/Aachen algorithm with Rfilt = "
|
---|
54 | << _Rfilt
|
---|
55 | << " (recomb. scheme deduced from jet, or E-scheme if not unique)";
|
---|
56 | } else {
|
---|
57 | ostr << _subjet_def.description();
|
---|
58 | }
|
---|
59 | ostr<< ", selection " << _selector.description();
|
---|
60 | if (_subtractor) {
|
---|
61 | ostr << ", subtractor: " << _subtractor->description();
|
---|
62 | } else if (_rho != 0) {
|
---|
63 | ostr << ", subtracting with rho = " << _rho;
|
---|
64 | }
|
---|
65 | return ostr.str();
|
---|
66 | }
|
---|
67 |
|
---|
68 |
|
---|
69 | // core functions
|
---|
70 | //----------------------------------------------------------------------
|
---|
71 |
|
---|
72 | // return a vector of subjets, which are the ones that would be kept
|
---|
73 | // by the filtering
|
---|
74 | PseudoJet Filter::result(const PseudoJet &jet) const {
|
---|
75 | // start by getting the list of subjets (including a list of sanity
|
---|
76 | // checks)
|
---|
77 | // NB: subjets is empty to begin with (see the comment for
|
---|
78 | // _set_filtered_elements_cafilt)
|
---|
79 | vector<PseudoJet> subjets;
|
---|
80 | JetDefinition subjet_def;
|
---|
81 | bool discard_area;
|
---|
82 | // NB: on return, subjet_def is set to the jet definition actually
|
---|
83 | // used (so that we can make use of its recombination scheme
|
---|
84 | // when joining the jets to be kept).
|
---|
85 | _set_filtered_elements(jet, subjets, subjet_def, discard_area);
|
---|
86 |
|
---|
87 | // now build the vector of kept and rejected subjets
|
---|
88 | vector<PseudoJet> kept, rejected;
|
---|
89 | // Note that in the following line we make a copy of the _selector
|
---|
90 | // to avoid issues with needing a mutable _selector
|
---|
91 | Selector selector_copy = _selector;
|
---|
92 | if (selector_copy.takes_reference()) selector_copy.set_reference(jet);
|
---|
93 | selector_copy.sift(subjets, kept, rejected);
|
---|
94 |
|
---|
95 | // gather the info under the form of a PseudoJet
|
---|
96 | return _finalise(jet, kept, rejected, subjet_def, discard_area);
|
---|
97 | }
|
---|
98 |
|
---|
99 |
|
---|
100 | // sets filtered_elements to be all the subjets on which filtering will work
|
---|
101 | void Filter::_set_filtered_elements(const PseudoJet & jet,
|
---|
102 | vector<PseudoJet> & filtered_elements,
|
---|
103 | JetDefinition & subjet_def,
|
---|
104 | bool & discard_area) const {
|
---|
105 | // sanity checks
|
---|
106 | //-------------------------------------------------------------------
|
---|
107 | // make sure that the jet has constituents
|
---|
108 | if (! jet.has_constituents())
|
---|
109 | throw Error("Filter can only be applied on jets having constituents");
|
---|
110 |
|
---|
111 | // for a whole variety of checks, we shall need the "recursive"
|
---|
112 | // pieces of the jet (the jet itself or recursing down to its most
|
---|
113 | // fundamental pieces).
|
---|
114 | // So we do compute these once and for all
|
---|
115 | vector<PseudoJet> all_pieces; //.clear();
|
---|
116 | if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0))
|
---|
117 | throw Error("Attempt to filter a jet that has no associated ClusterSequence or is not a superposition of jets associated with a ClusterSequence");
|
---|
118 |
|
---|
119 | // if the filter uses subtraction, make sure we have a CS that supports area and has
|
---|
120 | // explicit ghosts
|
---|
121 | if (_uses_subtraction()) {
|
---|
122 | if (!jet.has_area())
|
---|
123 | throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without area info for the original jet");
|
---|
124 |
|
---|
125 | if (!_check_explicit_ghosts(all_pieces))
|
---|
126 | throw Error("Attempt to filter and subtract (non-zero rho or subtractor) without explicit ghosts");
|
---|
127 | }
|
---|
128 |
|
---|
129 | // if we're dealing with a dynamic determination of the filtering
|
---|
130 | // radius, do it now
|
---|
131 | if ((_Rfilt>=0) || (_Rfiltfunc)){
|
---|
132 | double Rfilt = (_Rfiltfunc) ? (*_Rfiltfunc)(jet) : _Rfilt;
|
---|
133 | const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
|
---|
134 | if (common_recombiner) {
|
---|
135 | if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
|
---|
136 | RecombinationScheme scheme =
|
---|
137 | static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
|
---|
138 | subjet_def = JetDefinition(cambridge_algorithm, Rfilt, scheme);
|
---|
139 | } else {
|
---|
140 | subjet_def = JetDefinition(cambridge_algorithm, Rfilt, common_recombiner);
|
---|
141 | }
|
---|
142 | } else {
|
---|
143 | subjet_def = JetDefinition(cambridge_algorithm, Rfilt);
|
---|
144 | }
|
---|
145 | } else {
|
---|
146 | subjet_def = _subjet_def;
|
---|
147 | }
|
---|
148 |
|
---|
149 | // get the jet definition to be use and whether we can apply our
|
---|
150 | // simplified C/A+C/A filter
|
---|
151 | //
|
---|
152 | // we apply C/A clustering iff
|
---|
153 | // - the request subjet_def is C/A
|
---|
154 | // - the jet is either directly coming from C/A or if it is a
|
---|
155 | // superposition of C/A jets
|
---|
156 | // - the pieces agree with the recombination scheme of subjet_def
|
---|
157 | //------------------------------------------------------------------
|
---|
158 | bool simple_cafilt = _check_ca(all_pieces);
|
---|
159 |
|
---|
160 | // extract the subjets
|
---|
161 | //-------------------------------------------------------------------
|
---|
162 | discard_area = false;
|
---|
163 | if (simple_cafilt){
|
---|
164 | // first make sure that 'filtered_elements' is empty
|
---|
165 | filtered_elements.clear();
|
---|
166 | _set_filtered_elements_cafilt(jet, filtered_elements, subjet_def.R());
|
---|
167 | // in the following case, areas can be erroneous and will be discarded
|
---|
168 | discard_area = (!_uses_subtraction()) && (jet.has_area()) && (!_check_explicit_ghosts(all_pieces));
|
---|
169 | } else {
|
---|
170 | // here we'll simply recluster the jets.
|
---|
171 | //
|
---|
172 | // To include an area support we need
|
---|
173 | // - the jet to have an area
|
---|
174 | // - subtraction requested or explicit ghosts
|
---|
175 | bool do_areas = (jet.has_area()) && ((_uses_subtraction()) || (_check_explicit_ghosts(all_pieces)));
|
---|
176 | _set_filtered_elements_generic(jet, filtered_elements, subjet_def, do_areas);
|
---|
177 | }
|
---|
178 |
|
---|
179 | // order the filtered elements in pt
|
---|
180 | filtered_elements = sorted_by_pt(filtered_elements);
|
---|
181 | }
|
---|
182 |
|
---|
183 | // set the filtered elements in the simple case of C/A+C/A
|
---|
184 | //
|
---|
185 | // WATCH OUT: this could be recursively called, so filtered elements
|
---|
186 | // of 'jet' are APPENDED to 'filtered_elements'
|
---|
187 | void Filter::_set_filtered_elements_cafilt(const PseudoJet & jet,
|
---|
188 | vector<PseudoJet> & filtered_elements,
|
---|
189 | double Rfilt) const{
|
---|
190 | // we know that the jet is either a C/A jet or a superposition of
|
---|
191 | // such pieces
|
---|
192 | if (jet.has_associated_cluster_sequence()){
|
---|
193 | // just extract the exclusive subjets of 'jet'
|
---|
194 | const ClusterSequence *cs = jet.associated_cluster_sequence();
|
---|
195 | vector<PseudoJet> local_fe;
|
---|
196 |
|
---|
197 | double dcut = Rfilt / cs->jet_def().R();
|
---|
198 | if (dcut>=1.0){
|
---|
199 | local_fe.push_back(jet);
|
---|
200 | } else {
|
---|
201 | dcut *= dcut;
|
---|
202 | local_fe = jet.exclusive_subjets(dcut);
|
---|
203 | }
|
---|
204 |
|
---|
205 | // subtract the jets if needed
|
---|
206 | // Note that this one would work on pieces!!
|
---|
207 | //-----------------------------------------------------------------
|
---|
208 | if (_uses_subtraction()){
|
---|
209 | const ClusterSequenceAreaBase * csab = jet.validated_csab();
|
---|
210 | for (unsigned int i=0;i<local_fe.size();i++) {
|
---|
211 | if (_subtractor) {
|
---|
212 | local_fe[i] = (*_subtractor)(local_fe[i]);
|
---|
213 | } else {
|
---|
214 | local_fe[i] = csab->subtracted_jet(local_fe[i], _rho);
|
---|
215 | }
|
---|
216 | }
|
---|
217 | }
|
---|
218 |
|
---|
219 | copy(local_fe.begin(), local_fe.end(), back_inserter(filtered_elements));
|
---|
220 | return;
|
---|
221 | }
|
---|
222 |
|
---|
223 | // just recurse into the pieces
|
---|
224 | const vector<PseudoJet> & pieces = jet.pieces();
|
---|
225 | for (vector<PseudoJet>::const_iterator it = pieces.begin();
|
---|
226 | it!=pieces.end(); it++)
|
---|
227 | _set_filtered_elements_cafilt(*it, filtered_elements, Rfilt);
|
---|
228 | }
|
---|
229 |
|
---|
230 |
|
---|
231 | // set the filtered elements in the generic re-clustering case (wo
|
---|
232 | // subtraction)
|
---|
233 | void Filter::_set_filtered_elements_generic(const PseudoJet & jet,
|
---|
234 | vector<PseudoJet> & filtered_elements,
|
---|
235 | const JetDefinition & subjet_def,
|
---|
236 | bool do_areas) const{
|
---|
237 | // create a new, internal, ClusterSequence from the jet constituents
|
---|
238 | // get the subjets directly from there
|
---|
239 | //
|
---|
240 | // If the jet has area support then we separate the ghosts from the
|
---|
241 | // "regular" particles so the subjets will also have area
|
---|
242 | // support. Note that we do this regardless of whether rho is zero
|
---|
243 | // or not.
|
---|
244 | //
|
---|
245 | // Note that to be able to separate the ghosts, one needs explicit
|
---|
246 | // ghosts!!
|
---|
247 | // ---------------------------------------------------------------
|
---|
248 | if (do_areas){
|
---|
249 | vector<PseudoJet> all_constituents = jet.constituents();
|
---|
250 | vector<PseudoJet> regular_constituents, ghosts;
|
---|
251 |
|
---|
252 | for (vector<PseudoJet>::iterator it = all_constituents.begin();
|
---|
253 | it != all_constituents.end(); it++){
|
---|
254 | if (it->is_pure_ghost())
|
---|
255 | ghosts.push_back(*it);
|
---|
256 | else
|
---|
257 | regular_constituents.push_back(*it);
|
---|
258 | }
|
---|
259 |
|
---|
260 | // figure out the ghost area from the 1st ghost (if none, any value
|
---|
261 | // would probably do as the area will be 0 and subtraction will have
|
---|
262 | // no effect!)
|
---|
263 | double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
|
---|
264 | ClusterSequenceActiveAreaExplicitGhosts * csa
|
---|
265 | = new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents,
|
---|
266 | subjet_def,
|
---|
267 | ghosts, ghost_area);
|
---|
268 |
|
---|
269 | // get the subjets: we use the subtracted or unsubtracted ones
|
---|
270 | // depending on rho or _subtractor being non-zero
|
---|
271 | if (_uses_subtraction()) {
|
---|
272 | if (_subtractor) {
|
---|
273 | filtered_elements = (*_subtractor)(csa->inclusive_jets());
|
---|
274 | } else {
|
---|
275 | filtered_elements = csa->subtracted_jets(_rho);
|
---|
276 | }
|
---|
277 | } else {
|
---|
278 | filtered_elements = csa->inclusive_jets();
|
---|
279 | }
|
---|
280 |
|
---|
281 | // allow the cs to be deleted when it's no longer used
|
---|
282 | csa->delete_self_when_unused();
|
---|
283 | } else {
|
---|
284 | ClusterSequence * cs = new ClusterSequence(jet.constituents(), subjet_def);
|
---|
285 | filtered_elements = cs->inclusive_jets();
|
---|
286 | // allow the cs to be deleted when it's no longer used
|
---|
287 | cs->delete_self_when_unused();
|
---|
288 | }
|
---|
289 | }
|
---|
290 |
|
---|
291 |
|
---|
292 | // gather the information about what is kept and rejected under the
|
---|
293 | // form of a PseudoJet with a special ClusterSequenceInfo
|
---|
294 | PseudoJet Filter::_finalise(const PseudoJet & /*jet*/,
|
---|
295 | vector<PseudoJet> & kept,
|
---|
296 | vector<PseudoJet> & rejected,
|
---|
297 | const JetDefinition & subjet_def,
|
---|
298 | const bool discard_area) const {
|
---|
299 | // figure out which recombiner to use
|
---|
300 | const JetDefinition::Recombiner &rec = *(subjet_def.recombiner());
|
---|
301 |
|
---|
302 | // create an appropriate structure and transfer the info to it
|
---|
303 | PseudoJet filtered_jet = join<StructureType>(kept, rec);
|
---|
304 | StructureType *fs = (StructureType*) filtered_jet.structure_non_const_ptr();
|
---|
305 | // fs->_original_jet = jet;
|
---|
306 | fs->_rejected = rejected;
|
---|
307 |
|
---|
308 | if (discard_area){
|
---|
309 | // safety check: make sure there is an area to discard!!!
|
---|
310 | assert(fs->_area_4vector_ptr);
|
---|
311 | delete fs->_area_4vector_ptr;
|
---|
312 | fs->_area_4vector_ptr=0;
|
---|
313 | }
|
---|
314 |
|
---|
315 | return filtered_jet;
|
---|
316 | }
|
---|
317 |
|
---|
318 | // various checks
|
---|
319 | //----------------------------------------------------------------------
|
---|
320 |
|
---|
321 | // get the pieces down to the fundamental pieces
|
---|
322 | //
|
---|
323 | // Note that this just checks that there is an associated CS to the
|
---|
324 | // fundamental pieces, not that it is still valid
|
---|
325 | bool Filter::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
|
---|
326 | if (jet.has_associated_cluster_sequence()){
|
---|
327 | all_pieces.push_back(jet);
|
---|
328 | return true;
|
---|
329 | }
|
---|
330 |
|
---|
331 | if (jet.has_pieces()){
|
---|
332 | const vector<PseudoJet> pieces = jet.pieces();
|
---|
333 | for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
|
---|
334 | if (!_get_all_pieces(*it, all_pieces)) return false;
|
---|
335 | return true;
|
---|
336 | }
|
---|
337 |
|
---|
338 | return false;
|
---|
339 | }
|
---|
340 |
|
---|
341 | // get the common recombiner to all pieces (NULL if none)
|
---|
342 | //
|
---|
343 | // Note that if the jet has an associated cluster sequence that is no
|
---|
344 | // longer valid, an error will be thrown (needed since it could be the
|
---|
345 | // 1st check called after the enumeration of the pieces)
|
---|
346 | const JetDefinition::Recombiner* Filter::_get_common_recombiner(const vector<PseudoJet> &all_pieces) const{
|
---|
347 | const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
|
---|
348 | for (unsigned int i=1; i<all_pieces.size(); i++)
|
---|
349 | if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)) return NULL;
|
---|
350 |
|
---|
351 | return jd_ref.recombiner();
|
---|
352 | }
|
---|
353 |
|
---|
354 | // check if the jet (or all its pieces) have explicit ghosts
|
---|
355 | // (assuming the jet has area support
|
---|
356 | //
|
---|
357 | // Note that if the jet has an associated cluster sequence that is no
|
---|
358 | // longer valid, an error will be thrown (needed since it could be the
|
---|
359 | // 1st check called after the enumeration of the pieces)
|
---|
360 | bool Filter::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
|
---|
361 | for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
|
---|
362 | if (! it->validated_csab()->has_explicit_ghosts()) return false;
|
---|
363 | return true;
|
---|
364 | }
|
---|
365 |
|
---|
366 | // check if one can apply the simplification for C/A subjets
|
---|
367 | //
|
---|
368 | // This includes:
|
---|
369 | // - the subjet definition asks for C/A subjets
|
---|
370 | // - all the pieces share the same CS
|
---|
371 | // - that CS is C/A with the same recombiner as the subjet def
|
---|
372 | // - the filtering radius is not larger than any of the pairwise
|
---|
373 | // distance between the pieces
|
---|
374 | //
|
---|
375 | // Note that if the jet has an associated cluster sequence that is no
|
---|
376 | // longer valid, an error will be thrown (needed since it could be the
|
---|
377 | // 1st check called after the enumeration of the pieces)
|
---|
378 | bool Filter::_check_ca(const vector<PseudoJet> &all_pieces) const{
|
---|
379 | if (_subjet_def.jet_algorithm() != cambridge_algorithm) return false;
|
---|
380 |
|
---|
381 | // check that the 1st of all the pieces (we're sure there is at
|
---|
382 | // least one) is coming from a C/A clustering. Then check that all
|
---|
383 | // the following pieces share the same ClusterSequence
|
---|
384 | const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
|
---|
385 | if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
|
---|
386 | for (unsigned int i=1; i<all_pieces.size(); i++)
|
---|
387 | if (all_pieces[i].validated_cs() != cs_ref) return false;
|
---|
388 |
|
---|
389 | // check that the 1st peice has the same recombiner as the one used
|
---|
390 | // for the subjet clustering
|
---|
391 | // Note that since they share the same CS, checking the 2st one is enough
|
---|
392 | if (!cs_ref->jet_def().has_same_recombiner(_subjet_def)) return false;
|
---|
393 |
|
---|
394 | // we also have to make sure that the filtering radius is not larger
|
---|
395 | // than any of the inter-pieces distance
|
---|
396 | double Rfilt2 = _subjet_def.R();
|
---|
397 | Rfilt2 *= Rfilt2;
|
---|
398 | for (unsigned int i=0; i<all_pieces.size()-1; i++){
|
---|
399 | for (unsigned int j=i+1; j<all_pieces.size(); j++){
|
---|
400 | if (all_pieces[i].squared_distance(all_pieces[j]) < Rfilt2) return false;
|
---|
401 | }
|
---|
402 | }
|
---|
403 |
|
---|
404 | return true;
|
---|
405 | }
|
---|
406 |
|
---|
407 | //----------------------------------------------------------------------
|
---|
408 | // FilterInterface implementation
|
---|
409 | //----------------------------------------------------------------------
|
---|
410 |
|
---|
411 |
|
---|
412 | FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
|
---|