1 | #ifndef __SELECTOR_HH__ |
---|
2 | #define __SELECTOR_HH__ |
---|
3 | |
---|
4 | //STARTHEADER |
---|
5 | // $Id: Selector.hh 2687 2011-11-14 11:17:51Z soyez $ |
---|
6 | // |
---|
7 | // Copyright (c) 2009-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez |
---|
8 | // |
---|
9 | //---------------------------------------------------------------------- |
---|
10 | // This file is part of FastJet. |
---|
11 | // |
---|
12 | // FastJet is free software; you can redistribute it and/or modify |
---|
13 | // it under the terms of the GNU General Public License as published by |
---|
14 | // the Free Software Foundation; either version 2 of the License, or |
---|
15 | // (at your option) any later version. |
---|
16 | // |
---|
17 | // The algorithms that underlie FastJet have required considerable |
---|
18 | // development and are described in hep-ph/0512210. If you use |
---|
19 | // FastJet as part of work towards a scientific publication, please |
---|
20 | // include a citation to the FastJet paper. |
---|
21 | // |
---|
22 | // FastJet is distributed in the hope that it will be useful, |
---|
23 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
24 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
25 | // GNU General Public License for more details. |
---|
26 | // |
---|
27 | // You should have received a copy of the GNU General Public License |
---|
28 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>. |
---|
29 | //---------------------------------------------------------------------- |
---|
30 | //ENDHEADER |
---|
31 | |
---|
32 | #include "fastjet/PseudoJet.hh" |
---|
33 | #include "fastjet/RangeDefinition.hh" // for initialisation from a RangeDefinition |
---|
34 | #include <limits> |
---|
35 | #include <cmath> |
---|
36 | |
---|
37 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
38 | |
---|
39 | //---------------------------------------------------------------------- |
---|
40 | /// @ingroup selectors |
---|
41 | /// \class Selector |
---|
42 | /// Class that encodes information about cuts and other selection |
---|
43 | /// criteria that can be applied to PseudoJet(s). |
---|
44 | /// |
---|
45 | class Selector; |
---|
46 | //---------------------------------------------------------------------- |
---|
47 | |
---|
48 | /// @ingroup selectors |
---|
49 | /// \class SelectorWorker |
---|
50 | /// default selector worker is an abstract virtual base class |
---|
51 | /// |
---|
52 | /// The Selector class is only an interface, it is the SelectorWorker |
---|
53 | /// that really does the work. To implement various selectors, one |
---|
54 | /// thus has to overload this class. |
---|
55 | class SelectorWorker { |
---|
56 | public: |
---|
57 | //---------------------------------------------------------- |
---|
58 | // fundamental info |
---|
59 | //---------------------------------------------------------- |
---|
60 | /// default dtor |
---|
61 | virtual ~SelectorWorker() {} |
---|
62 | |
---|
63 | //---------------------------------------------------------- |
---|
64 | // basic operations for checking what gets selected |
---|
65 | //---------------------------------------------------------- |
---|
66 | |
---|
67 | /// returns true if a given object passes the selection criterion. |
---|
68 | /// This has to be overloaded by derived workers |
---|
69 | virtual bool pass(const PseudoJet & jet) const = 0; |
---|
70 | |
---|
71 | /// For each jet that does not pass the cuts, this routine sets the |
---|
72 | /// pointer to 0. |
---|
73 | /// |
---|
74 | /// It does not assume that the PseudoJet* passed as argumetn are not NULL |
---|
75 | virtual void terminator(std::vector<const PseudoJet *> & jets) const { |
---|
76 | for (unsigned i = 0; i < jets.size(); i++) { |
---|
77 | if (jets[i] && !pass(*jets[i])) jets[i] = NULL; |
---|
78 | } |
---|
79 | } |
---|
80 | |
---|
81 | /// returns true if this can be applied jet by jet |
---|
82 | virtual bool applies_jet_by_jet() const {return true;} |
---|
83 | |
---|
84 | /// returns a description of the worker |
---|
85 | virtual std::string description() const {return "missing description";} |
---|
86 | |
---|
87 | |
---|
88 | //---------------------------------------------------------- |
---|
89 | // operations for dealing with reference jets |
---|
90 | //---------------------------------------------------------- |
---|
91 | |
---|
92 | /// returns true if the worker is defined with respect to a reference jet |
---|
93 | virtual bool takes_reference() const { return false;} |
---|
94 | |
---|
95 | /// sets the reference jet for the selector |
---|
96 | /// NB: "reference" is commented to avoid unused-variable compiler warnings |
---|
97 | virtual void set_reference(const PseudoJet & /*reference*/){ |
---|
98 | throw Error("set_reference(...) cannot be used for a selector worker that does not take a reference"); |
---|
99 | } |
---|
100 | |
---|
101 | /// return a copy of the current object. |
---|
102 | /// |
---|
103 | /// This function is only called for objects that take a reference and need |
---|
104 | /// not be reimplemented otherwise. |
---|
105 | virtual SelectorWorker* copy(){ |
---|
106 | throw Error("this SelectorWorker has nothing to copy"); |
---|
107 | } |
---|
108 | |
---|
109 | //---------------------------------------------------------- |
---|
110 | // operations for area and extent |
---|
111 | //---------------------------------------------------------- |
---|
112 | |
---|
113 | /// returns the rapidity range for which it may return "true" |
---|
114 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const { |
---|
115 | rapmax = std::numeric_limits<double>::infinity(); |
---|
116 | rapmin = -rapmax; |
---|
117 | } |
---|
118 | |
---|
119 | /// check if it is a geometric selector (i.e. only puts constraints |
---|
120 | /// on rapidity and azimuthal angle) |
---|
121 | virtual bool is_geometric() const { return false;} |
---|
122 | |
---|
123 | /// check if it has a finite area |
---|
124 | virtual bool has_finite_area() const; |
---|
125 | |
---|
126 | /// check if it has an analytically computable area |
---|
127 | virtual bool has_known_area() const { return false;} |
---|
128 | |
---|
129 | /// if it has a computable area, return it |
---|
130 | virtual double known_area() const{ |
---|
131 | throw Error("this selector has no computable area"); |
---|
132 | } |
---|
133 | }; |
---|
134 | |
---|
135 | //---------------------------------------------------------------------- |
---|
136 | // class Selector |
---|
137 | // |
---|
138 | // Class that encodes information about cuts that |
---|
139 | class Selector{ |
---|
140 | public: |
---|
141 | /// default constructor produces a Selector whose action is undefined |
---|
142 | /// (any attempt to use it will lead to an error) |
---|
143 | Selector() {} |
---|
144 | |
---|
145 | /// constructor that causes the Selector to use the supplied worker |
---|
146 | /// |
---|
147 | /// Note that the Selector takes ownership of the pointer to the |
---|
148 | /// worker (and so will delete automatically when appropriate). |
---|
149 | Selector(SelectorWorker * worker_in) {_worker.reset(worker_in);} |
---|
150 | |
---|
151 | |
---|
152 | /// ctor from a RangeDefinition |
---|
153 | /// |
---|
154 | /// This is provided for backward compatibility and will be removed in |
---|
155 | /// a future major release of FastJet |
---|
156 | /// |
---|
157 | /// Watch out that the Selector will only hold a pointer to the |
---|
158 | /// range so the selector will crash if one tries to use it after |
---|
159 | /// the range has gone out of scope. We thus strongly advise against |
---|
160 | /// the direct use of this constructor. |
---|
161 | Selector(const RangeDefinition &range); |
---|
162 | |
---|
163 | /// dummy virtual dtor |
---|
164 | virtual ~Selector(){} |
---|
165 | |
---|
166 | /// return true if the jet passes the selection |
---|
167 | bool pass(const PseudoJet & jet) const { |
---|
168 | if (!validated_worker()->applies_jet_by_jet()) { |
---|
169 | throw Error("Cannot apply this selector to an individual jet"); |
---|
170 | } |
---|
171 | return _worker->pass(jet); |
---|
172 | } |
---|
173 | |
---|
174 | /// an operator way of knowing whether a given jet passes the selection or not |
---|
175 | bool operator()(const PseudoJet & jet) const { |
---|
176 | return pass(jet); |
---|
177 | } |
---|
178 | |
---|
179 | /// Return a count of the objects that pass the selection. |
---|
180 | /// |
---|
181 | /// This will often be more efficient that getting the vector of objects that |
---|
182 | /// passes and then evaluating the size of the vector |
---|
183 | unsigned int count(const std::vector<PseudoJet> & jets) const; |
---|
184 | |
---|
185 | /// sift the input jets into two vectors -- those that pass the selector |
---|
186 | /// and those that do not |
---|
187 | void sift(const std::vector<PseudoJet> & jets, |
---|
188 | std::vector<PseudoJet> & jets_that_pass, |
---|
189 | std::vector<PseudoJet> & jets_that_fail) const; |
---|
190 | |
---|
191 | /// returns true if this can be applied jet by jet |
---|
192 | bool applies_jet_by_jet() const { |
---|
193 | return validated_worker()->applies_jet_by_jet(); |
---|
194 | } |
---|
195 | |
---|
196 | /// returns a vector with the jets that pass the selection |
---|
197 | std::vector<PseudoJet> operator()(const std::vector<PseudoJet> & jets) const; |
---|
198 | |
---|
199 | /// For each jet that does not pass the cuts, this routine sets the |
---|
200 | /// pointer to 0. |
---|
201 | /// |
---|
202 | /// It is legitimate for some (or all) of the pointers that are |
---|
203 | /// passed to already be NULL. |
---|
204 | virtual void nullify_non_selected(std::vector<const PseudoJet *> & jets) const { |
---|
205 | validated_worker()->terminator(jets); |
---|
206 | } |
---|
207 | |
---|
208 | /// returns the rapidity range for which it may return "true" |
---|
209 | void get_rapidity_extent(double &rapmin, double &rapmax) const { |
---|
210 | return validated_worker()->get_rapidity_extent(rapmin, rapmax); |
---|
211 | } |
---|
212 | |
---|
213 | /// returns a textual description of the selector |
---|
214 | std::string description() const { |
---|
215 | return validated_worker()->description(); |
---|
216 | } |
---|
217 | |
---|
218 | /// returns true if it is a geometric selector (i.e. one that only puts |
---|
219 | /// constraints on rapidities and azimuthal angles) |
---|
220 | bool is_geometric() const{ |
---|
221 | return validated_worker()->is_geometric(); |
---|
222 | } |
---|
223 | |
---|
224 | /// returns true if it has a meaningful and finite area (i.e. the |
---|
225 | /// Selector has the property that is_geometric() returns true and |
---|
226 | /// the rapidity extent is finite). |
---|
227 | bool has_finite_area() const{ |
---|
228 | return validated_worker()->has_finite_area(); |
---|
229 | } |
---|
230 | |
---|
231 | /// returns the rapidity-phi area associated with the Selector |
---|
232 | /// (throws InvalidArea if the area does not make sense). |
---|
233 | /// |
---|
234 | /// If the result is not known analytically, the area will be |
---|
235 | /// estimated using a pseudo Monte Carlo method (as for jet areas), |
---|
236 | /// using the default ghost area from the GhostedAreaSpec class |
---|
237 | /// (0.01). The Monte Carlo estimate involves a time penalty |
---|
238 | /// proportional to the ratio of the rapidity extent of the Selector |
---|
239 | /// divided by the ghost area. |
---|
240 | double area() const; |
---|
241 | |
---|
242 | /// returns the rapidity-phi area associated with the Selector |
---|
243 | /// (throws InvalidArea if the area does not make sense). |
---|
244 | /// |
---|
245 | /// The behaviour is the as with the area() call, but with the |
---|
246 | /// ability to additionally specify the ghost area to be used in the |
---|
247 | /// case of a Monte Carlo area evaluation. |
---|
248 | /// |
---|
249 | double area(double ghost_area) const; |
---|
250 | |
---|
251 | /// returns a (reference to) the underlying worker's shared pointer |
---|
252 | const SharedPtr<SelectorWorker> & worker() const {return _worker;} |
---|
253 | |
---|
254 | /// returns a worker if there is a valid one, otherwise throws an InvalidWorker error |
---|
255 | const SelectorWorker* validated_worker() const { |
---|
256 | const SelectorWorker* worker_ptr = _worker.get(); |
---|
257 | if (worker_ptr == 0) throw InvalidWorker(); |
---|
258 | return worker_ptr; |
---|
259 | } |
---|
260 | |
---|
261 | /// returns true if this can be applied jet by jet |
---|
262 | bool takes_reference() const { |
---|
263 | return validated_worker()->takes_reference(); |
---|
264 | } |
---|
265 | |
---|
266 | /// set the reference jet for this Selector |
---|
267 | const Selector & set_reference(const PseudoJet &reference){ |
---|
268 | |
---|
269 | // if the worker does not take a reference jet, do nothing |
---|
270 | if (! validated_worker()->takes_reference()){ |
---|
271 | return *this; |
---|
272 | } |
---|
273 | |
---|
274 | // since this is a non-const operation, make sure we have a |
---|
275 | // correct behaviour with respect to shared workers |
---|
276 | _copy_worker_if_needed(); |
---|
277 | |
---|
278 | _worker->set_reference(reference); |
---|
279 | return *this; |
---|
280 | } |
---|
281 | |
---|
282 | /// class that gets throw when a Selector is applied despite it not |
---|
283 | /// having a valid underlying worker. |
---|
284 | class InvalidWorker : public Error { |
---|
285 | public: |
---|
286 | InvalidWorker() : Error("Attempt to use Selector with no valid underlying worker") {} |
---|
287 | }; |
---|
288 | |
---|
289 | /// class that gets throw when a Selector is applied despite it not |
---|
290 | /// having a valid underlying worker. |
---|
291 | class InvalidArea : public Error { |
---|
292 | public: |
---|
293 | InvalidArea() : Error("Attempt to obtain area from Selector for which this is not meaningful") {} |
---|
294 | }; |
---|
295 | |
---|
296 | // some operators (applying directly on a Selector) |
---|
297 | //---------------------------------------------------------------------- |
---|
298 | /// For 2 Selectors a and b, a &= b is eauivalent to a = a && b; |
---|
299 | Selector & operator &=(const Selector & b); |
---|
300 | |
---|
301 | /// For 2 Selectors a and b, a |= b is eauivalent to a = a || b; |
---|
302 | Selector & operator |=(const Selector & b); |
---|
303 | |
---|
304 | |
---|
305 | protected: |
---|
306 | /// Helper for copying selector workers if needed |
---|
307 | /// |
---|
308 | /// The following is needed if we want to modify a selectors that |
---|
309 | /// shares a worker with another selector. In that case, we need to |
---|
310 | /// get another copy of the worker to avoid interferences |
---|
311 | /// |
---|
312 | /// Note that any non-const operation has to call this to behave |
---|
313 | /// correctly w.r.t shared workers! |
---|
314 | void _copy_worker_if_needed(){ |
---|
315 | // do nothing if there's a sinlge user of the worker |
---|
316 | if (_worker.unique()) return; |
---|
317 | |
---|
318 | // call the worker's copy |
---|
319 | //std::cout << "will make a copy of " << description() << std::endl; |
---|
320 | _worker.reset(_worker->copy()); |
---|
321 | } |
---|
322 | |
---|
323 | private: |
---|
324 | SharedPtr<SelectorWorker> _worker; ///< the underlying worker |
---|
325 | }; |
---|
326 | |
---|
327 | |
---|
328 | //---------------------------------------------------------------------- |
---|
329 | // a list of specific selectors |
---|
330 | //---------------------------------------------------------------------- |
---|
331 | |
---|
332 | /// \addtogroup selectors |
---|
333 | /// @{ |
---|
334 | |
---|
335 | |
---|
336 | // fundamental selectors |
---|
337 | //---------------------------------------------------------------------- |
---|
338 | |
---|
339 | // "identity" selector that lets everything pass |
---|
340 | Selector SelectorIdentity(); |
---|
341 | |
---|
342 | // logical operations |
---|
343 | //---------------------------------------------------------------------- |
---|
344 | |
---|
345 | /// logical not applied on a selector |
---|
346 | /// |
---|
347 | /// This will keep objects that do not pass the 's' selector |
---|
348 | Selector operator!(const Selector & s); |
---|
349 | |
---|
350 | /// logical or between two selectors |
---|
351 | /// |
---|
352 | /// this will keep the objects that are selected by s1 or s2 |
---|
353 | Selector operator ||(const Selector & s1, const Selector & s2); |
---|
354 | |
---|
355 | |
---|
356 | /// logical and between two selectors |
---|
357 | /// |
---|
358 | /// this will keep the objects that are selected by both s1 and s2 |
---|
359 | /// |
---|
360 | /// watch out: for both s1 and s2, the selection is applied on the |
---|
361 | /// original list of objects. For successive applications of two |
---|
362 | /// selectors (convolution/multiplication) see the operator * |
---|
363 | Selector operator&&(const Selector & s1, const Selector & s2); |
---|
364 | |
---|
365 | /// successive application of 2 selectors |
---|
366 | /// |
---|
367 | /// Apply the selector s2, then the selector s1. |
---|
368 | /// |
---|
369 | /// watch out: the operator * acts like an operator product i.e. does |
---|
370 | /// not commute. The order of its arguments is therefore important. |
---|
371 | /// Whenever they commute (in particluar, when they apply jet by |
---|
372 | /// jet), this would have the same effect as the logical &&. |
---|
373 | Selector operator*(const Selector & s1, const Selector & s2); |
---|
374 | |
---|
375 | |
---|
376 | // selection with kinematic cuts |
---|
377 | //---------------------------------------------------------------------- |
---|
378 | Selector SelectorPtMin(double ptmin); ///< select objects with pt >= ptmin |
---|
379 | Selector SelectorPtMax(double ptmax); ///< select objects with pt <= ptmax |
---|
380 | Selector SelectorPtRange(double ptmin, double ptmax); ///< select objects with ptmin <= pt <= ptmax |
---|
381 | |
---|
382 | Selector SelectorEtMin(double Etmin); ///< select objects with Et >= Etmin |
---|
383 | Selector SelectorEtMax(double Etmax); ///< select objects with Et <= Etmax |
---|
384 | Selector SelectorEtRange(double Etmin, double Etmax); ///< select objects with Etmin <= Et <= Etmax |
---|
385 | |
---|
386 | Selector SelectorEMin(double Emin); ///< select objects with E >= Emin |
---|
387 | Selector SelectorEMax(double Emax); ///< select objects with E <= Emax |
---|
388 | Selector SelectorERange(double Emin, double Emax); ///< select objects with Emin <= E <= Emax |
---|
389 | |
---|
390 | Selector SelectorMassMin(double Mmin); ///< select objects with Mass >= Mmin |
---|
391 | Selector SelectorMassMax(double Mmax); ///< select objects with Mass <= Mmax |
---|
392 | Selector SelectorMassRange(double Mmin, double Mmax); ///< select objects with Mmin <= Mass <= Mmax |
---|
393 | |
---|
394 | Selector SelectorRapMin(double rapmin); ///< select objects with rap >= rapmin |
---|
395 | Selector SelectorRapMax(double rapmax); ///< select objects with rap <= rapmax |
---|
396 | Selector SelectorRapRange(double rapmin, double rapmax); ///< select objects with rapmin <= rap <= rapmax |
---|
397 | |
---|
398 | Selector SelectorAbsRapMin(double absrapmin); ///< select objects with |rap| >= absrapmin |
---|
399 | Selector SelectorAbsRapMax(double absrapmax); ///< select objects with |rap| <= absrapmax |
---|
400 | Selector SelectorAbsRapRange(double absrapmin, double absrapmax); ///< select objects with absrapmin <= |rap| <= absrapmax |
---|
401 | |
---|
402 | Selector SelectorEtaMin(double etamin); ///< select objects with eta >= etamin |
---|
403 | Selector SelectorEtaMax(double etamax); ///< select objects with eta <= etamax |
---|
404 | Selector SelectorEtaRange(double etamin, double etamax); ///< select objects with etamin <= eta <= etamax |
---|
405 | |
---|
406 | Selector SelectorAbsEtaMin(double absetamin); ///< select objects with |eta| >= absetamin |
---|
407 | Selector SelectorAbsEtaMax(double absetamax); ///< select objects with |eta| <= absetamax |
---|
408 | Selector SelectorAbsEtaRange(double absetamin, double absetamax); ///< select objects with absetamin <= |eta| <= absetamax |
---|
409 | |
---|
410 | Selector SelectorPhiRange(double phimin, double phimax); ///< select objects with phimin <= phi <= phimax |
---|
411 | |
---|
412 | /// select objects with rapmin <= rap <= rapmax && phimin <= phi <= phimax |
---|
413 | /// |
---|
414 | /// Note that this is essentially a combination of SelectorRapRange |
---|
415 | /// and SelectorPhiRange. We provide it as a Selector on its own in |
---|
416 | /// order to use the known area (which would otherwise be lost by the && |
---|
417 | /// operator) |
---|
418 | Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax); |
---|
419 | |
---|
420 | /// select the n hardest objects |
---|
421 | Selector SelectorNHardest(unsigned int n); |
---|
422 | |
---|
423 | |
---|
424 | // Selectors that take (require) a reference jet. |
---|
425 | //---------------------------------------------------------------------- |
---|
426 | |
---|
427 | /// select objets within a distance 'radius' from the location of the |
---|
428 | /// reference jet, set by Selector::set_reference(...) |
---|
429 | Selector SelectorCircle(const double & radius); |
---|
430 | |
---|
431 | /// select objets with distance from the reference jet is between 'radius_in' |
---|
432 | /// and 'radius_out'; the reference jet is set by Selector::set_reference(...) |
---|
433 | Selector SelectorDoughnut(const double & radius_in, const double & radius_out); |
---|
434 | |
---|
435 | /// select objets within a rapidity distance 'half_width' from the |
---|
436 | /// location of the reference jet, set by Selector::set_reference(...) |
---|
437 | Selector SelectorStrip(const double & half_width); |
---|
438 | |
---|
439 | /// select objets within rapidity distance 'half_rap_width' from the |
---|
440 | /// reference jet and azimuthal-angle distance within 'half_phi_width'; the |
---|
441 | /// reference jet is set by Selector::set_reference(...) |
---|
442 | Selector SelectorRectangle(const double & half_rap_width, const double & half_phi_width); |
---|
443 | |
---|
444 | |
---|
445 | /// select objects that carry at least a fraction "fraction" of the |
---|
446 | /// reference jet. The reference jet must have been set with |
---|
447 | /// Selector::set_reference(...) |
---|
448 | Selector SelectorPtFractionMin(double fraction); |
---|
449 | |
---|
450 | |
---|
451 | // additional (mostly helper) selectors |
---|
452 | //---------------------------------------------------------------------- |
---|
453 | |
---|
454 | /// select PseudoJet with 0 momentum |
---|
455 | Selector SelectorIsZero(); |
---|
456 | |
---|
457 | /// select objects that are (or are only made of) ghosts. |
---|
458 | /// PseudoJets for which has_area() are considered non-pure-ghost. |
---|
459 | Selector SelectorIsPureGhost(); |
---|
460 | |
---|
461 | /// @} |
---|
462 | |
---|
463 | FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh |
---|
464 | |
---|
465 | #endif // __SELECTOR_HH__ |
---|
466 | |
---|