1 | //STARTHEADER |
---|
2 | // $Id: Selector.cc 859 2012-11-28 01:49:23Z 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 | |
---|
30 | #include <sstream> |
---|
31 | #include <algorithm> |
---|
32 | #include "fastjet/Selector.hh" |
---|
33 | #include "fastjet/GhostedAreaSpec.hh" // for area support |
---|
34 | |
---|
35 | using namespace std; |
---|
36 | |
---|
37 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
38 | |
---|
39 | //---------------------------------------------------------------------- |
---|
40 | // implementations of some of the more complex bits of Selector |
---|
41 | //---------------------------------------------------------------------- |
---|
42 | |
---|
43 | // implementation of the operator() acting on a vector of jets |
---|
44 | std::vector<PseudoJet> Selector::operator()(const std::vector<PseudoJet> & jets) const { |
---|
45 | std::vector<PseudoJet> result; |
---|
46 | const SelectorWorker * worker_local = validated_worker(); |
---|
47 | if (worker_local->applies_jet_by_jet()) { |
---|
48 | //if (false) { |
---|
49 | // for workers that apply jet by jet, this is more efficient |
---|
50 | for (std::vector<PseudoJet>::const_iterator jet = jets.begin(); |
---|
51 | jet != jets.end(); jet++) { |
---|
52 | if (worker_local->pass(*jet)) result.push_back(*jet); |
---|
53 | } |
---|
54 | } else { |
---|
55 | // for workers that can only be applied to entire vectors, |
---|
56 | // go through the following |
---|
57 | std::vector<const PseudoJet *> jetptrs(jets.size()); |
---|
58 | for (unsigned i = 0; i < jets.size(); i++) { |
---|
59 | jetptrs[i] = & jets[i]; |
---|
60 | } |
---|
61 | worker_local->terminator(jetptrs); |
---|
62 | for (unsigned i = 0; i < jetptrs.size(); i++) { |
---|
63 | if (jetptrs[i]) result.push_back(jets[i]); |
---|
64 | } |
---|
65 | } |
---|
66 | return result; |
---|
67 | } |
---|
68 | |
---|
69 | |
---|
70 | //---------------------------------------------------------------------- |
---|
71 | // count the number of jets that pass the cuts |
---|
72 | unsigned int Selector::count(const std::vector<PseudoJet> & jets) const { |
---|
73 | unsigned n = 0; |
---|
74 | const SelectorWorker * worker_local = validated_worker(); |
---|
75 | |
---|
76 | // separate strategies according to whether the worker applies jet by jet |
---|
77 | if (worker_local->applies_jet_by_jet()) { |
---|
78 | for (unsigned i = 0; i < jets.size(); i++) { |
---|
79 | if (worker_local->pass(jets[i])) n++; |
---|
80 | } |
---|
81 | } else { |
---|
82 | std::vector<const PseudoJet *> jetptrs(jets.size()); |
---|
83 | for (unsigned i = 0; i < jets.size(); i++) { |
---|
84 | jetptrs[i] = & jets[i]; |
---|
85 | } |
---|
86 | worker_local->terminator(jetptrs); |
---|
87 | for (unsigned i = 0; i < jetptrs.size(); i++) { |
---|
88 | if (jetptrs[i]) n++; |
---|
89 | } |
---|
90 | } |
---|
91 | |
---|
92 | return n; |
---|
93 | } |
---|
94 | |
---|
95 | |
---|
96 | //---------------------------------------------------------------------- |
---|
97 | // sift the input jets into two vectors -- those that pass the selector |
---|
98 | // and those that do not |
---|
99 | void Selector::sift(const std::vector<PseudoJet> & jets, |
---|
100 | std::vector<PseudoJet> & jets_that_pass, |
---|
101 | std::vector<PseudoJet> & jets_that_fail |
---|
102 | ) const { |
---|
103 | const SelectorWorker * worker_local = validated_worker(); |
---|
104 | |
---|
105 | jets_that_pass.clear(); |
---|
106 | jets_that_fail.clear(); |
---|
107 | |
---|
108 | // separate strategies according to whether the worker applies jet by jet |
---|
109 | if (worker_local->applies_jet_by_jet()) { |
---|
110 | for (unsigned i = 0; i < jets.size(); i++) { |
---|
111 | if (worker_local->pass(jets[i])) { |
---|
112 | jets_that_pass.push_back(jets[i]); |
---|
113 | } else { |
---|
114 | jets_that_fail.push_back(jets[i]); |
---|
115 | } |
---|
116 | } |
---|
117 | } else { |
---|
118 | std::vector<const PseudoJet *> jetptrs(jets.size()); |
---|
119 | for (unsigned i = 0; i < jets.size(); i++) { |
---|
120 | jetptrs[i] = & jets[i]; |
---|
121 | } |
---|
122 | worker_local->terminator(jetptrs); |
---|
123 | for (unsigned i = 0; i < jetptrs.size(); i++) { |
---|
124 | if (jetptrs[i]) { |
---|
125 | jets_that_pass.push_back(jets[i]); |
---|
126 | } else { |
---|
127 | jets_that_fail.push_back(jets[i]); |
---|
128 | } |
---|
129 | } |
---|
130 | } |
---|
131 | } |
---|
132 | |
---|
133 | // area using default ghost area |
---|
134 | double Selector::area() const{ |
---|
135 | return area(gas::def_ghost_area); |
---|
136 | } |
---|
137 | |
---|
138 | // implementation of the Selector's area function |
---|
139 | double Selector::area(double ghost_area) const{ |
---|
140 | if (! is_geometric()) throw InvalidArea(); |
---|
141 | |
---|
142 | // has area will already check we've got a valid worker |
---|
143 | if (_worker->has_known_area()) return _worker->known_area(); |
---|
144 | |
---|
145 | // generate a set of "ghosts" |
---|
146 | double rapmin, rapmax; |
---|
147 | get_rapidity_extent(rapmin, rapmax); |
---|
148 | GhostedAreaSpec ghost_spec(rapmin, rapmax, 1, ghost_area); |
---|
149 | std::vector<PseudoJet> ghosts; |
---|
150 | ghost_spec.add_ghosts(ghosts); |
---|
151 | |
---|
152 | // check what passes the selection |
---|
153 | return ghost_spec.ghost_area() * ((*this)(ghosts)).size(); |
---|
154 | } |
---|
155 | |
---|
156 | |
---|
157 | //---------------------------------------------------------------------- |
---|
158 | // implementations of some of the more complex bits of SelectorWorker |
---|
159 | //---------------------------------------------------------------------- |
---|
160 | // check if it has a finite area |
---|
161 | bool SelectorWorker::has_finite_area() const { |
---|
162 | if (! is_geometric()) return false; |
---|
163 | double rapmin, rapmax; |
---|
164 | get_rapidity_extent(rapmin, rapmax); |
---|
165 | return (rapmax != std::numeric_limits<double>::infinity()) |
---|
166 | && (-rapmin != std::numeric_limits<double>::infinity()); |
---|
167 | } |
---|
168 | |
---|
169 | |
---|
170 | |
---|
171 | //---------------------------------------------------------------------- |
---|
172 | // very basic set of selectors (at the moment just the identity!) |
---|
173 | //---------------------------------------------------------------------- |
---|
174 | |
---|
175 | //---------------------------------------------------------------------- |
---|
176 | /// helper for selecting the n hardest jets |
---|
177 | class SW_Identity : public SelectorWorker { |
---|
178 | public: |
---|
179 | /// ctor with specification of the number of objects to keep |
---|
180 | SW_Identity(){} |
---|
181 | |
---|
182 | /// just let everything pass |
---|
183 | virtual bool pass(const PseudoJet &) const { |
---|
184 | return true; |
---|
185 | } |
---|
186 | |
---|
187 | /// For each jet that does not pass the cuts, this routine sets the |
---|
188 | /// pointer to 0. |
---|
189 | virtual void terminator(vector<const PseudoJet *> &) const { |
---|
190 | // everything passes, hence nothing to nullify |
---|
191 | return; |
---|
192 | } |
---|
193 | |
---|
194 | /// returns a description of the worker |
---|
195 | virtual string description() const { return "Identity";} |
---|
196 | |
---|
197 | /// strictly speaking, this is geometric |
---|
198 | virtual bool is_geometric() const { return true;} |
---|
199 | }; |
---|
200 | |
---|
201 | |
---|
202 | // returns an "identity" selector that lets everything pass |
---|
203 | Selector SelectorIdentity() { |
---|
204 | return Selector(new SW_Identity); |
---|
205 | } |
---|
206 | |
---|
207 | |
---|
208 | //---------------------------------------------------------------------- |
---|
209 | // selector and workers for operators |
---|
210 | //---------------------------------------------------------------------- |
---|
211 | |
---|
212 | //---------------------------------------------------------------------- |
---|
213 | /// helper for combining selectors with a logical not |
---|
214 | class SW_Not : public SelectorWorker { |
---|
215 | public: |
---|
216 | /// ctor |
---|
217 | SW_Not(const Selector & s) : _s(s) {} |
---|
218 | |
---|
219 | /// return a copy of the current object |
---|
220 | virtual SelectorWorker* copy(){ return new SW_Not(*this);} |
---|
221 | |
---|
222 | /// returns true if a given object passes the selection criterium |
---|
223 | /// this has to be overloaded by derived workers |
---|
224 | virtual bool pass(const PseudoJet & jet) const { |
---|
225 | // make sure that the "pass" can be applied on both selectors |
---|
226 | if (!applies_jet_by_jet()) |
---|
227 | throw Error("Cannot apply this selector worker to an individual jet"); |
---|
228 | |
---|
229 | return ! _s.pass(jet); |
---|
230 | } |
---|
231 | |
---|
232 | /// returns true if this can be applied jet by jet |
---|
233 | virtual bool applies_jet_by_jet() const {return _s.applies_jet_by_jet();} |
---|
234 | |
---|
235 | /// select the jets in the list that pass both selectors |
---|
236 | virtual void terminator(vector<const PseudoJet *> & jets) const { |
---|
237 | // if we can apply the selector jet-by-jet, call the base selector |
---|
238 | // that does exactly that |
---|
239 | if (applies_jet_by_jet()){ |
---|
240 | SelectorWorker::terminator(jets); |
---|
241 | return; |
---|
242 | } |
---|
243 | |
---|
244 | // check the effect of the selector we want to negate |
---|
245 | vector<const PseudoJet *> s_jets = jets; |
---|
246 | _s.worker()->terminator(s_jets); |
---|
247 | |
---|
248 | // now apply the negation: all the jets that pass the base |
---|
249 | // selector (i.e. are not NULL) have to be set to NULL |
---|
250 | for (unsigned int i=0; i<s_jets.size(); i++){ |
---|
251 | if (s_jets[i]) jets[i] = NULL; |
---|
252 | } |
---|
253 | } |
---|
254 | |
---|
255 | /// returns a description of the worker |
---|
256 | virtual string description() const { |
---|
257 | ostringstream ostr; |
---|
258 | ostr << "!(" << _s.description() << ")"; |
---|
259 | return ostr.str(); |
---|
260 | } |
---|
261 | |
---|
262 | /// is geometric if the underlying selector is |
---|
263 | virtual bool is_geometric() const { return _s.is_geometric();} |
---|
264 | |
---|
265 | /// returns true if the worker can be set_referenced |
---|
266 | virtual bool takes_reference() const { return _s.takes_reference();} |
---|
267 | |
---|
268 | /// set the reference jet for this selector |
---|
269 | virtual void set_reference(const PseudoJet &ref) { _s.set_reference(ref);} |
---|
270 | |
---|
271 | protected: |
---|
272 | Selector _s; |
---|
273 | }; |
---|
274 | |
---|
275 | |
---|
276 | // logical not applied on a selector |
---|
277 | Selector operator!(const Selector & s) { |
---|
278 | return Selector(new SW_Not(s)); |
---|
279 | } |
---|
280 | |
---|
281 | |
---|
282 | //---------------------------------------------------------------------- |
---|
283 | /// Base class for binary operators |
---|
284 | class SW_BinaryOperator: public SelectorWorker { |
---|
285 | public: |
---|
286 | /// ctor |
---|
287 | SW_BinaryOperator(const Selector & s1, const Selector & s2) : _s1(s1), _s2(s2) { |
---|
288 | // stores info for more efficient access to the selector's properties |
---|
289 | |
---|
290 | // we can apply jet by jet only if this is the case for both sub-selectors |
---|
291 | _applies_jet_by_jet = _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet(); |
---|
292 | |
---|
293 | // the selector takes a reference if either of the sub-selectors does |
---|
294 | _takes_reference = _s1.takes_reference() || _s2.takes_reference(); |
---|
295 | |
---|
296 | // we have a well-defined area provided the two objects have one |
---|
297 | _is_geometric = _s1.is_geometric() && _s2.is_geometric(); |
---|
298 | } |
---|
299 | |
---|
300 | /// returns true if this can be applied jet by jet |
---|
301 | virtual bool applies_jet_by_jet() const {return _applies_jet_by_jet;} |
---|
302 | |
---|
303 | /// returns true if this takes a reference jet |
---|
304 | virtual bool takes_reference() const{ |
---|
305 | return _takes_reference; |
---|
306 | } |
---|
307 | |
---|
308 | /// sets the reference jet |
---|
309 | virtual void set_reference(const PseudoJet ¢re){ |
---|
310 | _s1.set_reference(centre); |
---|
311 | _s2.set_reference(centre); |
---|
312 | } |
---|
313 | |
---|
314 | /// check if it has a finite area |
---|
315 | virtual bool is_geometric() const { return _is_geometric;} |
---|
316 | |
---|
317 | protected: |
---|
318 | Selector _s1, _s2; |
---|
319 | bool _applies_jet_by_jet; |
---|
320 | bool _takes_reference; |
---|
321 | bool _is_geometric; |
---|
322 | }; |
---|
323 | |
---|
324 | |
---|
325 | |
---|
326 | //---------------------------------------------------------------------- |
---|
327 | /// helper for combining selectors with a logical and |
---|
328 | class SW_And: public SW_BinaryOperator { |
---|
329 | public: |
---|
330 | /// ctor |
---|
331 | SW_And(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2){} |
---|
332 | |
---|
333 | /// return a copy of this |
---|
334 | virtual SelectorWorker* copy(){ return new SW_And(*this);} |
---|
335 | |
---|
336 | /// returns true if a given object passes the selection criterium |
---|
337 | /// this has to be overloaded by derived workers |
---|
338 | virtual bool pass(const PseudoJet & jet) const { |
---|
339 | // make sure that the "pass" can be applied on both selectors |
---|
340 | if (!applies_jet_by_jet()) |
---|
341 | throw Error("Cannot apply this selector worker to an individual jet"); |
---|
342 | |
---|
343 | return _s1.pass(jet) && _s2.pass(jet); |
---|
344 | } |
---|
345 | |
---|
346 | /// select the jets in the list that pass both selectors |
---|
347 | virtual void terminator(vector<const PseudoJet *> & jets) const { |
---|
348 | // if we can apply the selector jet-by-jet, call the base selector |
---|
349 | // that does exactly that |
---|
350 | if (applies_jet_by_jet()){ |
---|
351 | SelectorWorker::terminator(jets); |
---|
352 | return; |
---|
353 | } |
---|
354 | |
---|
355 | // check the effect of the first selector |
---|
356 | vector<const PseudoJet *> s1_jets = jets; |
---|
357 | _s1.worker()->terminator(s1_jets); |
---|
358 | |
---|
359 | // apply the second |
---|
360 | _s2.worker()->terminator(jets); |
---|
361 | |
---|
362 | // terminate the jets that wiuld be terminated by _s1 |
---|
363 | for (unsigned int i=0; i<jets.size(); i++){ |
---|
364 | if (! s1_jets[i]) jets[i] = NULL; |
---|
365 | } |
---|
366 | } |
---|
367 | |
---|
368 | /// returns the rapidity range for which it may return "true" |
---|
369 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const { |
---|
370 | double s1min, s1max, s2min, s2max; |
---|
371 | _s1.get_rapidity_extent(s1min, s1max); |
---|
372 | _s2.get_rapidity_extent(s2min, s2max); |
---|
373 | rapmax = min(s1max, s2max); |
---|
374 | rapmin = max(s1min, s2min); |
---|
375 | } |
---|
376 | |
---|
377 | /// returns a description of the worker |
---|
378 | virtual string description() const { |
---|
379 | ostringstream ostr; |
---|
380 | ostr << "(" << _s1.description() << " && " << _s2.description() << ")"; |
---|
381 | return ostr.str(); |
---|
382 | } |
---|
383 | }; |
---|
384 | |
---|
385 | |
---|
386 | // logical and between two selectors |
---|
387 | Selector operator&&(const Selector & s1, const Selector & s2) { |
---|
388 | return Selector(new SW_And(s1,s2)); |
---|
389 | } |
---|
390 | |
---|
391 | |
---|
392 | |
---|
393 | //---------------------------------------------------------------------- |
---|
394 | /// helper for combining selectors with a logical or |
---|
395 | class SW_Or: public SW_BinaryOperator { |
---|
396 | public: |
---|
397 | /// ctor |
---|
398 | SW_Or(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2) {} |
---|
399 | |
---|
400 | /// return a copy of this |
---|
401 | virtual SelectorWorker* copy(){ return new SW_Or(*this);} |
---|
402 | |
---|
403 | /// returns true if a given object passes the selection criterium |
---|
404 | /// this has to be overloaded by derived workers |
---|
405 | virtual bool pass(const PseudoJet & jet) const { |
---|
406 | // make sure that the "pass" can be applied on both selectors |
---|
407 | if (!applies_jet_by_jet()) |
---|
408 | throw Error("Cannot apply this selector worker to an individual jet"); |
---|
409 | |
---|
410 | return _s1.pass(jet) || _s2.pass(jet); |
---|
411 | } |
---|
412 | |
---|
413 | /// returns true if this can be applied jet by jet |
---|
414 | virtual bool applies_jet_by_jet() const { |
---|
415 | // watch out, even though it's the "OR" selector, to be applied jet |
---|
416 | // by jet, both the base selectors need to be jet-by-jet-applicable, |
---|
417 | // so the use of a && in the line below |
---|
418 | return _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet(); |
---|
419 | } |
---|
420 | |
---|
421 | /// select the jets in the list that pass both selectors |
---|
422 | virtual void terminator(vector<const PseudoJet *> & jets) const { |
---|
423 | // if we can apply the selector jet-by-jet, call the base selector |
---|
424 | // that does exactly that |
---|
425 | if (applies_jet_by_jet()){ |
---|
426 | SelectorWorker::terminator(jets); |
---|
427 | return; |
---|
428 | } |
---|
429 | |
---|
430 | // check the effect of the first selector |
---|
431 | vector<const PseudoJet *> s1_jets = jets; |
---|
432 | _s1.worker()->terminator(s1_jets); |
---|
433 | |
---|
434 | // apply the second |
---|
435 | _s2.worker()->terminator(jets); |
---|
436 | |
---|
437 | // resurrect any jet that has been terminated by the second one |
---|
438 | // and not by the first one |
---|
439 | for (unsigned int i=0; i<jets.size(); i++){ |
---|
440 | if (s1_jets[i]) jets[i] = s1_jets[i]; |
---|
441 | } |
---|
442 | } |
---|
443 | |
---|
444 | /// returns a description of the worker |
---|
445 | virtual string description() const { |
---|
446 | ostringstream ostr; |
---|
447 | ostr << "(" << _s1.description() << " || " << _s2.description() << ")"; |
---|
448 | return ostr.str(); |
---|
449 | } |
---|
450 | |
---|
451 | /// returns the rapidity range for which it may return "true" |
---|
452 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const { |
---|
453 | double s1min, s1max, s2min, s2max; |
---|
454 | _s1.get_rapidity_extent(s1min, s1max); |
---|
455 | _s2.get_rapidity_extent(s2min, s2max); |
---|
456 | rapmax = max(s1max, s2max); |
---|
457 | rapmin = min(s1min, s2min); |
---|
458 | } |
---|
459 | }; |
---|
460 | |
---|
461 | |
---|
462 | // logical or between two selectors |
---|
463 | Selector operator ||(const Selector & s1, const Selector & s2) { |
---|
464 | return Selector(new SW_Or(s1,s2)); |
---|
465 | } |
---|
466 | |
---|
467 | //---------------------------------------------------------------------- |
---|
468 | /// helper for multiplying two selectors (in an operator-like way) |
---|
469 | class SW_Mult: public SW_And { |
---|
470 | public: |
---|
471 | /// ctor |
---|
472 | SW_Mult(const Selector & s1, const Selector & s2) : SW_And(s1,s2) {} |
---|
473 | |
---|
474 | /// return a copy of this |
---|
475 | virtual SelectorWorker* copy(){ return new SW_Mult(*this);} |
---|
476 | |
---|
477 | /// select the jets in the list that pass both selectors |
---|
478 | virtual void terminator(vector<const PseudoJet *> & jets) const { |
---|
479 | // if we can apply the selector jet-by-jet, call the base selector |
---|
480 | // that does exactly that |
---|
481 | if (applies_jet_by_jet()){ |
---|
482 | SelectorWorker::terminator(jets); |
---|
483 | return; |
---|
484 | } |
---|
485 | |
---|
486 | // first apply _s2 |
---|
487 | _s2.worker()->terminator(jets); |
---|
488 | |
---|
489 | // then apply _s1 |
---|
490 | _s1.worker()->terminator(jets); |
---|
491 | } |
---|
492 | |
---|
493 | /// returns a description of the worker |
---|
494 | virtual string description() const { |
---|
495 | ostringstream ostr; |
---|
496 | ostr << "(" << _s1.description() << " * " << _s2.description() << ")"; |
---|
497 | return ostr.str(); |
---|
498 | } |
---|
499 | }; |
---|
500 | |
---|
501 | |
---|
502 | // logical and between two selectors |
---|
503 | Selector operator*(const Selector & s1, const Selector & s2) { |
---|
504 | return Selector(new SW_Mult(s1,s2)); |
---|
505 | } |
---|
506 | |
---|
507 | |
---|
508 | //---------------------------------------------------------------------- |
---|
509 | // selector and workers for kinematic cuts |
---|
510 | //---------------------------------------------------------------------- |
---|
511 | |
---|
512 | //---------------------------------------------------------------------- |
---|
513 | // a series of basic classes that allow easy implementations of |
---|
514 | // min, max and ranges on a quantity-to-be-defined |
---|
515 | |
---|
516 | // generic holder for a quantity |
---|
517 | class QuantityBase{ |
---|
518 | public: |
---|
519 | QuantityBase(double q) : _q(q){} |
---|
520 | virtual ~QuantityBase(){} |
---|
521 | virtual double operator()(const PseudoJet & jet ) const =0; |
---|
522 | virtual string description() const =0; |
---|
523 | virtual bool is_geometric() const { return false;} |
---|
524 | virtual double comparison_value() const {return _q;} |
---|
525 | virtual double description_value() const {return comparison_value();} |
---|
526 | protected: |
---|
527 | double _q; |
---|
528 | }; |
---|
529 | |
---|
530 | // generic holder for a squared quantity |
---|
531 | class QuantitySquareBase : public QuantityBase{ |
---|
532 | public: |
---|
533 | QuantitySquareBase(double sqrtq) : QuantityBase(sqrtq*sqrtq), _sqrtq(sqrtq){} |
---|
534 | virtual double description_value() const {return _sqrtq;} |
---|
535 | protected: |
---|
536 | double _sqrtq; |
---|
537 | }; |
---|
538 | |
---|
539 | // generic_quantity >= minimum |
---|
540 | template<typename QuantityType> |
---|
541 | class SW_QuantityMin : public SelectorWorker{ |
---|
542 | public: |
---|
543 | /// detfault ctor (initialises the pt cut) |
---|
544 | SW_QuantityMin(double qmin) : _qmin(qmin) {} |
---|
545 | |
---|
546 | /// returns true is the given object passes the selection pt cut |
---|
547 | virtual bool pass(const PseudoJet & jet) const {return _qmin(jet) >= _qmin.comparison_value();} |
---|
548 | |
---|
549 | /// returns a description of the worker |
---|
550 | virtual string description() const { |
---|
551 | ostringstream ostr; |
---|
552 | ostr << _qmin.description() << " >= " << _qmin.description_value(); |
---|
553 | return ostr.str(); |
---|
554 | } |
---|
555 | |
---|
556 | virtual bool is_geometric() const { return _qmin.is_geometric();} |
---|
557 | |
---|
558 | protected: |
---|
559 | QuantityType _qmin; ///< the cut |
---|
560 | }; |
---|
561 | |
---|
562 | |
---|
563 | // generic_quantity <= maximum |
---|
564 | template<typename QuantityType> |
---|
565 | class SW_QuantityMax : public SelectorWorker { |
---|
566 | public: |
---|
567 | /// detfault ctor (initialises the pt cut) |
---|
568 | SW_QuantityMax(double qmax) : _qmax(qmax) {} |
---|
569 | |
---|
570 | /// returns true is the given object passes the selection pt cut |
---|
571 | virtual bool pass(const PseudoJet & jet) const {return _qmax(jet) <= _qmax.comparison_value();} |
---|
572 | |
---|
573 | /// returns a description of the worker |
---|
574 | virtual string description() const { |
---|
575 | ostringstream ostr; |
---|
576 | ostr << _qmax.description() << " <= " << _qmax.description_value(); |
---|
577 | return ostr.str(); |
---|
578 | } |
---|
579 | |
---|
580 | virtual bool is_geometric() const { return _qmax.is_geometric();} |
---|
581 | |
---|
582 | protected: |
---|
583 | QuantityType _qmax; ///< the cut |
---|
584 | }; |
---|
585 | |
---|
586 | |
---|
587 | // generic quantity in [minimum:maximum] |
---|
588 | template<typename QuantityType> |
---|
589 | class SW_QuantityRange : public SelectorWorker { |
---|
590 | public: |
---|
591 | /// detfault ctor (initialises the pt cut) |
---|
592 | SW_QuantityRange(double qmin, double qmax) : _qmin(qmin), _qmax(qmax) {} |
---|
593 | |
---|
594 | /// returns true is the given object passes the selection pt cut |
---|
595 | virtual bool pass(const PseudoJet & jet) const { |
---|
596 | double q = _qmin(jet); // we could identically use _qmax |
---|
597 | return (q >= _qmin.comparison_value()) && (q <= _qmax.comparison_value()); |
---|
598 | } |
---|
599 | |
---|
600 | /// returns a description of the worker |
---|
601 | virtual string description() const { |
---|
602 | ostringstream ostr; |
---|
603 | ostr << _qmin.description_value() << " <= " << _qmin.description() << " <= " << _qmax.description_value(); |
---|
604 | return ostr.str(); |
---|
605 | } |
---|
606 | |
---|
607 | virtual bool is_geometric() const { return _qmin.is_geometric();} |
---|
608 | |
---|
609 | protected: |
---|
610 | QuantityType _qmin; // the lower cut |
---|
611 | QuantityType _qmax; // the upper cut |
---|
612 | }; |
---|
613 | |
---|
614 | |
---|
615 | //---------------------------------------------------------------------- |
---|
616 | /// helper class for selecting on pt |
---|
617 | class QuantityPt2 : public QuantitySquareBase{ |
---|
618 | public: |
---|
619 | QuantityPt2(double pt) : QuantitySquareBase(pt){} |
---|
620 | virtual double operator()(const PseudoJet & jet ) const { return jet.perp2();} |
---|
621 | virtual string description() const {return "pt";} |
---|
622 | }; |
---|
623 | |
---|
624 | // returns a selector for a minimum pt |
---|
625 | Selector SelectorPtMin(double ptmin) { |
---|
626 | return Selector(new SW_QuantityMin<QuantityPt2>(ptmin)); |
---|
627 | } |
---|
628 | |
---|
629 | // returns a selector for a maximum pt |
---|
630 | Selector SelectorPtMax(double ptmax) { |
---|
631 | return Selector(new SW_QuantityMax<QuantityPt2>(ptmax)); |
---|
632 | } |
---|
633 | |
---|
634 | // returns a selector for a pt range |
---|
635 | Selector SelectorPtRange(double ptmin, double ptmax) { |
---|
636 | return Selector(new SW_QuantityRange<QuantityPt2>(ptmin, ptmax)); |
---|
637 | } |
---|
638 | |
---|
639 | |
---|
640 | //---------------------------------------------------------------------- |
---|
641 | /// helper class for selecting on transverse energy |
---|
642 | class QuantityEt2 : public QuantitySquareBase{ |
---|
643 | public: |
---|
644 | QuantityEt2(double Et) : QuantitySquareBase(Et){} |
---|
645 | virtual double operator()(const PseudoJet & jet ) const { return jet.Et2();} |
---|
646 | virtual string description() const {return "Et";} |
---|
647 | }; |
---|
648 | |
---|
649 | // returns a selector for a minimum Et |
---|
650 | Selector SelectorEtMin(double Etmin) { |
---|
651 | return Selector(new SW_QuantityMin<QuantityEt2>(Etmin)); |
---|
652 | } |
---|
653 | |
---|
654 | // returns a selector for a maximum Et |
---|
655 | Selector SelectorEtMax(double Etmax) { |
---|
656 | return Selector(new SW_QuantityMax<QuantityEt2>(Etmax)); |
---|
657 | } |
---|
658 | |
---|
659 | // returns a selector for a Et range |
---|
660 | Selector SelectorEtRange(double Etmin, double Etmax) { |
---|
661 | return Selector(new SW_QuantityRange<QuantityEt2>(Etmin, Etmax)); |
---|
662 | } |
---|
663 | |
---|
664 | |
---|
665 | //---------------------------------------------------------------------- |
---|
666 | /// helper class for selecting on energy |
---|
667 | class QuantityE : public QuantityBase{ |
---|
668 | public: |
---|
669 | QuantityE(double E) : QuantityBase(E){} |
---|
670 | virtual double operator()(const PseudoJet & jet ) const { return jet.E();} |
---|
671 | virtual string description() const {return "E";} |
---|
672 | }; |
---|
673 | |
---|
674 | // returns a selector for a minimum E |
---|
675 | Selector SelectorEMin(double Emin) { |
---|
676 | return Selector(new SW_QuantityMin<QuantityE>(Emin)); |
---|
677 | } |
---|
678 | |
---|
679 | // returns a selector for a maximum E |
---|
680 | Selector SelectorEMax(double Emax) { |
---|
681 | return Selector(new SW_QuantityMax<QuantityE>(Emax)); |
---|
682 | } |
---|
683 | |
---|
684 | // returns a selector for a E range |
---|
685 | Selector SelectorERange(double Emin, double Emax) { |
---|
686 | return Selector(new SW_QuantityRange<QuantityE>(Emin, Emax)); |
---|
687 | } |
---|
688 | |
---|
689 | |
---|
690 | //---------------------------------------------------------------------- |
---|
691 | /// helper class for selecting on mass |
---|
692 | class QuantityM2 : public QuantitySquareBase{ |
---|
693 | public: |
---|
694 | QuantityM2(double m) : QuantitySquareBase(m){} |
---|
695 | virtual double operator()(const PseudoJet & jet ) const { return jet.m2();} |
---|
696 | virtual string description() const {return "mass";} |
---|
697 | }; |
---|
698 | |
---|
699 | // returns a selector for a minimum mass |
---|
700 | Selector SelectorMassMin(double mmin) { |
---|
701 | return Selector(new SW_QuantityMin<QuantityM2>(mmin)); |
---|
702 | } |
---|
703 | |
---|
704 | // returns a selector for a maximum mass |
---|
705 | Selector SelectorMassMax(double mmax) { |
---|
706 | return Selector(new SW_QuantityMax<QuantityM2>(mmax)); |
---|
707 | } |
---|
708 | |
---|
709 | // returns a selector for a mass range |
---|
710 | Selector SelectorMassRange(double mmin, double mmax) { |
---|
711 | return Selector(new SW_QuantityRange<QuantityM2>(mmin, mmax)); |
---|
712 | } |
---|
713 | |
---|
714 | |
---|
715 | |
---|
716 | //---------------------------------------------------------------------- |
---|
717 | /// helper for selecting on rapidities: quantity |
---|
718 | class QuantityRap : public QuantityBase{ |
---|
719 | public: |
---|
720 | QuantityRap(double rap) : QuantityBase(rap){} |
---|
721 | virtual double operator()(const PseudoJet & jet ) const { return jet.rap();} |
---|
722 | virtual string description() const {return "rap";} |
---|
723 | virtual bool is_geometric() const { return true;} |
---|
724 | }; |
---|
725 | |
---|
726 | |
---|
727 | /// helper for selecting on rapidities: min |
---|
728 | class SW_RapMin : public SW_QuantityMin<QuantityRap>{ |
---|
729 | public: |
---|
730 | SW_RapMin(double rapmin) : SW_QuantityMin<QuantityRap>(rapmin){} |
---|
731 | virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{ |
---|
732 | rapmax = std::numeric_limits<double>::max(); |
---|
733 | rapmin = _qmin.comparison_value(); |
---|
734 | } |
---|
735 | }; |
---|
736 | |
---|
737 | /// helper for selecting on rapidities: max |
---|
738 | class SW_RapMax : public SW_QuantityMax<QuantityRap>{ |
---|
739 | public: |
---|
740 | SW_RapMax(double rapmax) : SW_QuantityMax<QuantityRap>(rapmax){} |
---|
741 | virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{ |
---|
742 | rapmax = _qmax.comparison_value(); |
---|
743 | rapmin = -std::numeric_limits<double>::max(); |
---|
744 | } |
---|
745 | }; |
---|
746 | |
---|
747 | /// helper for selecting on rapidities: range |
---|
748 | class SW_RapRange : public SW_QuantityRange<QuantityRap>{ |
---|
749 | public: |
---|
750 | SW_RapRange(double rapmin, double rapmax) : SW_QuantityRange<QuantityRap>(rapmin, rapmax){ |
---|
751 | assert(rapmin<=rapmax); |
---|
752 | } |
---|
753 | virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{ |
---|
754 | rapmax = _qmax.comparison_value(); |
---|
755 | rapmin = _qmin.comparison_value(); |
---|
756 | } |
---|
757 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
758 | virtual double known_area() const { |
---|
759 | return twopi * (_qmax.comparison_value()-_qmin.comparison_value()); |
---|
760 | } |
---|
761 | }; |
---|
762 | |
---|
763 | // returns a selector for a minimum rapidity |
---|
764 | Selector SelectorRapMin(double rapmin) { |
---|
765 | return Selector(new SW_RapMin(rapmin)); |
---|
766 | } |
---|
767 | |
---|
768 | // returns a selector for a maximum rapidity |
---|
769 | Selector SelectorRapMax(double rapmax) { |
---|
770 | return Selector(new SW_RapMax(rapmax)); |
---|
771 | } |
---|
772 | |
---|
773 | // returns a selector for a rapidity range |
---|
774 | Selector SelectorRapRange(double rapmin, double rapmax) { |
---|
775 | return Selector(new SW_RapRange(rapmin, rapmax)); |
---|
776 | } |
---|
777 | |
---|
778 | |
---|
779 | //---------------------------------------------------------------------- |
---|
780 | /// helper for selecting on |rapidities| |
---|
781 | class QuantityAbsRap : public QuantityBase{ |
---|
782 | public: |
---|
783 | QuantityAbsRap(double absrap) : QuantityBase(absrap){} |
---|
784 | virtual double operator()(const PseudoJet & jet ) const { return abs(jet.rap());} |
---|
785 | virtual string description() const {return "|rap|";} |
---|
786 | virtual bool is_geometric() const { return true;} |
---|
787 | }; |
---|
788 | |
---|
789 | |
---|
790 | /// helper for selecting on |rapidities|: max |
---|
791 | class SW_AbsRapMax : public SW_QuantityMax<QuantityAbsRap>{ |
---|
792 | public: |
---|
793 | SW_AbsRapMax(double absrapmax) : SW_QuantityMax<QuantityAbsRap>(absrapmax){} |
---|
794 | virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{ |
---|
795 | rapmax = _qmax.comparison_value(); |
---|
796 | rapmin = -_qmax.comparison_value(); |
---|
797 | } |
---|
798 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
799 | virtual double known_area() const { |
---|
800 | return twopi * 2 * _qmax.comparison_value(); |
---|
801 | } |
---|
802 | }; |
---|
803 | |
---|
804 | /// helper for selecting on |rapidities|: max |
---|
805 | class SW_AbsRapRange : public SW_QuantityRange<QuantityAbsRap>{ |
---|
806 | public: |
---|
807 | SW_AbsRapRange(double absrapmin, double absrapmax) : SW_QuantityRange<QuantityAbsRap>(absrapmin, absrapmax){} |
---|
808 | virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{ |
---|
809 | rapmax = _qmax.comparison_value(); |
---|
810 | rapmin = -_qmax.comparison_value(); |
---|
811 | } |
---|
812 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
813 | virtual double known_area() const { |
---|
814 | return twopi * 2 * (_qmax.comparison_value()-max(_qmin.comparison_value(),0.0)); // this should handle properly absrapmin<0 |
---|
815 | } |
---|
816 | }; |
---|
817 | |
---|
818 | // returns a selector for a minimum |rapidity| |
---|
819 | Selector SelectorAbsRapMin(double absrapmin) { |
---|
820 | return Selector(new SW_QuantityMin<QuantityAbsRap>(absrapmin)); |
---|
821 | } |
---|
822 | |
---|
823 | // returns a selector for a maximum |rapidity| |
---|
824 | Selector SelectorAbsRapMax(double absrapmax) { |
---|
825 | return Selector(new SW_AbsRapMax(absrapmax)); |
---|
826 | } |
---|
827 | |
---|
828 | // returns a selector for a |rapidity| range |
---|
829 | Selector SelectorAbsRapRange(double rapmin, double rapmax) { |
---|
830 | return Selector(new SW_AbsRapRange(rapmin, rapmax)); |
---|
831 | } |
---|
832 | |
---|
833 | |
---|
834 | //---------------------------------------------------------------------- |
---|
835 | /// helper for selecting on pseudo-rapidities |
---|
836 | class QuantityEta : public QuantityBase{ |
---|
837 | public: |
---|
838 | QuantityEta(double eta) : QuantityBase(eta){} |
---|
839 | virtual double operator()(const PseudoJet & jet ) const { return jet.eta();} |
---|
840 | virtual string description() const {return "eta";} |
---|
841 | // virtual bool is_geometric() const { return true;} // not strictly only y and phi-dependent |
---|
842 | }; |
---|
843 | |
---|
844 | // returns a selector for a pseudo-minimum rapidity |
---|
845 | Selector SelectorEtaMin(double etamin) { |
---|
846 | return Selector(new SW_QuantityMin<QuantityEta>(etamin)); |
---|
847 | } |
---|
848 | |
---|
849 | // returns a selector for a pseudo-maximum rapidity |
---|
850 | Selector SelectorEtaMax(double etamax) { |
---|
851 | return Selector(new SW_QuantityMax<QuantityEta>(etamax)); |
---|
852 | } |
---|
853 | |
---|
854 | // returns a selector for a pseudo-rapidity range |
---|
855 | Selector SelectorEtaRange(double etamin, double etamax) { |
---|
856 | return Selector(new SW_QuantityRange<QuantityEta>(etamin, etamax)); |
---|
857 | } |
---|
858 | |
---|
859 | |
---|
860 | //---------------------------------------------------------------------- |
---|
861 | /// helper for selecting on |pseudo-rapidities| |
---|
862 | class QuantityAbsEta : public QuantityBase{ |
---|
863 | public: |
---|
864 | QuantityAbsEta(double abseta) : QuantityBase(abseta){} |
---|
865 | virtual double operator()(const PseudoJet & jet ) const { return abs(jet.eta());} |
---|
866 | virtual string description() const {return "|eta|";} |
---|
867 | virtual bool is_geometric() const { return true;} |
---|
868 | }; |
---|
869 | |
---|
870 | // returns a selector for a minimum |pseudo-rapidity| |
---|
871 | Selector SelectorAbsEtaMin(double absetamin) { |
---|
872 | return Selector(new SW_QuantityMin<QuantityAbsEta>(absetamin)); |
---|
873 | } |
---|
874 | |
---|
875 | // returns a selector for a maximum |pseudo-rapidity| |
---|
876 | Selector SelectorAbsEtaMax(double absetamax) { |
---|
877 | return Selector(new SW_QuantityMax<QuantityAbsEta>(absetamax)); |
---|
878 | } |
---|
879 | |
---|
880 | // returns a selector for a |pseudo-rapidity| range |
---|
881 | Selector SelectorAbsEtaRange(double absetamin, double absetamax) { |
---|
882 | return Selector(new SW_QuantityRange<QuantityAbsEta>(absetamin, absetamax)); |
---|
883 | } |
---|
884 | |
---|
885 | |
---|
886 | //---------------------------------------------------------------------- |
---|
887 | /// helper for selecting on azimuthal angle |
---|
888 | /// |
---|
889 | /// Note that the bounds have to be specified as min<max |
---|
890 | /// phimin has to be > -2pi |
---|
891 | /// phimax has to be < 4pi |
---|
892 | class SW_PhiRange : public SelectorWorker { |
---|
893 | public: |
---|
894 | /// detfault ctor (initialises the pt cut) |
---|
895 | SW_PhiRange(double phimin, double phimax) : _phimin(phimin), _phimax(phimax){ |
---|
896 | assert(_phimin<_phimax); |
---|
897 | assert(_phimin>-twopi); |
---|
898 | assert(_phimax<2*twopi); |
---|
899 | |
---|
900 | _phispan = _phimax - _phimin; |
---|
901 | } |
---|
902 | |
---|
903 | /// returns true is the given object passes the selection pt cut |
---|
904 | virtual bool pass(const PseudoJet & jet) const { |
---|
905 | double dphi=jet.phi()-_phimin; |
---|
906 | if (dphi >= twopi) dphi -= twopi; |
---|
907 | if (dphi < 0) dphi += twopi; |
---|
908 | return (dphi <= _phispan); |
---|
909 | } |
---|
910 | |
---|
911 | /// returns a description of the worker |
---|
912 | virtual string description() const { |
---|
913 | ostringstream ostr; |
---|
914 | ostr << _phimin << " <= phi <= " << _phimax; |
---|
915 | return ostr.str(); |
---|
916 | } |
---|
917 | |
---|
918 | virtual bool is_geometric() const { return true;} |
---|
919 | |
---|
920 | protected: |
---|
921 | double _phimin; // the lower cut |
---|
922 | double _phimax; // the upper cut |
---|
923 | double _phispan; // the span of the range |
---|
924 | }; |
---|
925 | |
---|
926 | |
---|
927 | // returns a selector for a phi range |
---|
928 | Selector SelectorPhiRange(double phimin, double phimax) { |
---|
929 | return Selector(new SW_PhiRange(phimin, phimax)); |
---|
930 | } |
---|
931 | |
---|
932 | //---------------------------------------------------------------------- |
---|
933 | /// helper for selecting on both rapidity and azimuthal angle |
---|
934 | class SW_RapPhiRange : public SW_And{ |
---|
935 | public: |
---|
936 | SW_RapPhiRange(double rapmin, double rapmax, double phimin, double phimax) |
---|
937 | : SW_And(SelectorRapRange(rapmin, rapmax), SelectorPhiRange(phimin, phimax)){ |
---|
938 | _known_area = ((phimax-phimin > twopi) ? twopi : phimax-phimin) * (rapmax-rapmin); |
---|
939 | } |
---|
940 | |
---|
941 | /// if it has a computable area, return it |
---|
942 | virtual double known_area() const{ |
---|
943 | return _known_area; |
---|
944 | } |
---|
945 | |
---|
946 | protected: |
---|
947 | double _known_area; |
---|
948 | }; |
---|
949 | |
---|
950 | Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax) { |
---|
951 | return Selector(new SW_RapPhiRange(rapmin, rapmax, phimin, phimax)); |
---|
952 | } |
---|
953 | |
---|
954 | |
---|
955 | //---------------------------------------------------------------------- |
---|
956 | /// helper for selecting the n hardest jets |
---|
957 | class SW_NHardest : public SelectorWorker { |
---|
958 | public: |
---|
959 | /// ctor with specification of the number of objects to keep |
---|
960 | SW_NHardest(unsigned int n) : _n(n) {}; |
---|
961 | |
---|
962 | /// pass makes no sense here normally the parent selector will throw |
---|
963 | /// an error but for internal use in the SW, we'll throw one from |
---|
964 | /// here by security |
---|
965 | virtual bool pass(const PseudoJet &) const { |
---|
966 | if (!applies_jet_by_jet()) |
---|
967 | throw Error("Cannot apply this selector worker to an individual jet"); |
---|
968 | return false; |
---|
969 | } |
---|
970 | |
---|
971 | /// For each jet that does not pass the cuts, this routine sets the |
---|
972 | /// pointer to 0. |
---|
973 | virtual void terminator(vector<const PseudoJet *> & jets) const { |
---|
974 | // nothing to do if the size is too small |
---|
975 | if (jets.size() < _n) return; |
---|
976 | |
---|
977 | // do we want to first chech if things are already ordered before |
---|
978 | // going through the ordering process? For now, no. Maybe carry |
---|
979 | // out timing tests at some point to establish the optimal |
---|
980 | // strategy. |
---|
981 | |
---|
982 | vector<double> minus_pt2(jets.size()); |
---|
983 | vector<unsigned int> indices(jets.size()); |
---|
984 | |
---|
985 | for (unsigned int i=0; i<jets.size(); i++){ |
---|
986 | indices[i] = i; |
---|
987 | |
---|
988 | // we need to make sure that the object has not already been |
---|
989 | // nullified. Note that if we have less than _n jets, this |
---|
990 | // whole n-hardest selection will not have any effect. |
---|
991 | minus_pt2[i] = jets[i] ? -jets[i]->perp2() : 0.0; |
---|
992 | } |
---|
993 | |
---|
994 | IndexedSortHelper sort_helper(& minus_pt2); |
---|
995 | |
---|
996 | partial_sort(indices.begin(), indices.begin()+_n, indices.end(), sort_helper); |
---|
997 | |
---|
998 | for (unsigned int i=_n; i<jets.size(); i++) |
---|
999 | jets[indices[i]] = NULL; |
---|
1000 | } |
---|
1001 | |
---|
1002 | /// returns true if this can be applied jet by jet |
---|
1003 | virtual bool applies_jet_by_jet() const {return false;} |
---|
1004 | |
---|
1005 | /// returns a description of the worker |
---|
1006 | virtual string description() const { |
---|
1007 | ostringstream ostr; |
---|
1008 | ostr << _n << " hardest"; |
---|
1009 | return ostr.str(); |
---|
1010 | } |
---|
1011 | |
---|
1012 | protected: |
---|
1013 | unsigned int _n; |
---|
1014 | }; |
---|
1015 | |
---|
1016 | |
---|
1017 | // returns a selector for the n hardest jets |
---|
1018 | Selector SelectorNHardest(unsigned int n) { |
---|
1019 | return Selector(new SW_NHardest(n)); |
---|
1020 | } |
---|
1021 | |
---|
1022 | |
---|
1023 | |
---|
1024 | //---------------------------------------------------------------------- |
---|
1025 | // selector and workers for geometric ranges |
---|
1026 | //---------------------------------------------------------------------- |
---|
1027 | |
---|
1028 | //---------------------------------------------------------------------- |
---|
1029 | /// a generic class for objects that contain a position |
---|
1030 | class SW_WithReference : public SelectorWorker{ |
---|
1031 | public: |
---|
1032 | /// ctor |
---|
1033 | SW_WithReference() : _is_initialised(false){}; |
---|
1034 | |
---|
1035 | /// returns true if the worker takes a reference jet |
---|
1036 | virtual bool takes_reference() const { return true;} |
---|
1037 | |
---|
1038 | /// sets the reference jet |
---|
1039 | virtual void set_reference(const PseudoJet ¢re){ |
---|
1040 | _is_initialised = true; |
---|
1041 | _reference = centre; |
---|
1042 | } |
---|
1043 | |
---|
1044 | protected: |
---|
1045 | PseudoJet _reference; |
---|
1046 | bool _is_initialised; |
---|
1047 | }; |
---|
1048 | |
---|
1049 | //---------------------------------------------------------------------- |
---|
1050 | /// helper for selecting on objects within a distance 'radius' of a reference |
---|
1051 | class SW_Circle : public SW_WithReference { |
---|
1052 | public: |
---|
1053 | SW_Circle(const double &radius) : _radius2(radius*radius) {} |
---|
1054 | |
---|
1055 | /// return a copy of the current object |
---|
1056 | virtual SelectorWorker* copy(){ return new SW_Circle(*this);} |
---|
1057 | |
---|
1058 | /// returns true if a given object passes the selection criterium |
---|
1059 | /// this has to be overloaded by derived workers |
---|
1060 | virtual bool pass(const PseudoJet & jet) const { |
---|
1061 | // make sure the centre is initialised |
---|
1062 | if (! _is_initialised) |
---|
1063 | throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1064 | |
---|
1065 | return jet.squared_distance(_reference) <= _radius2; |
---|
1066 | } |
---|
1067 | |
---|
1068 | /// returns a description of the worker |
---|
1069 | virtual string description() const { |
---|
1070 | ostringstream ostr; |
---|
1071 | ostr << "distance from the centre <= " << sqrt(_radius2); |
---|
1072 | return ostr.str(); |
---|
1073 | } |
---|
1074 | |
---|
1075 | /// returns the rapidity range for which it may return "true" |
---|
1076 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{ |
---|
1077 | // make sure the centre is initialised |
---|
1078 | if (! _is_initialised) |
---|
1079 | throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1080 | |
---|
1081 | rapmax = _reference.rap()+sqrt(_radius2); |
---|
1082 | rapmin = _reference.rap()-sqrt(_radius2); |
---|
1083 | } |
---|
1084 | |
---|
1085 | virtual bool is_geometric() const { return true;} ///< implies a finite area |
---|
1086 | virtual bool has_finite_area() const { return true;} ///< regardless of the reference |
---|
1087 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
1088 | virtual double known_area() const { |
---|
1089 | return pi * _radius2; |
---|
1090 | } |
---|
1091 | |
---|
1092 | protected: |
---|
1093 | double _radius2; |
---|
1094 | }; |
---|
1095 | |
---|
1096 | |
---|
1097 | // select on objets within a distance 'radius' of a variable location |
---|
1098 | Selector SelectorCircle(const double & radius) { |
---|
1099 | return Selector(new SW_Circle(radius)); |
---|
1100 | } |
---|
1101 | |
---|
1102 | |
---|
1103 | //---------------------------------------------------------------------- |
---|
1104 | /// helper for selecting on objects with a distance to a reference |
---|
1105 | /// betwene 'radius_in' and 'radius_out' |
---|
1106 | class SW_Doughnut : public SW_WithReference { |
---|
1107 | public: |
---|
1108 | SW_Doughnut(const double &radius_in, const double &radius_out) |
---|
1109 | : _radius_in2(radius_in*radius_in), _radius_out2(radius_out*radius_out) {} |
---|
1110 | |
---|
1111 | /// return a copy of the current object |
---|
1112 | virtual SelectorWorker* copy(){ return new SW_Doughnut(*this);} |
---|
1113 | |
---|
1114 | /// returns true if a given object passes the selection criterium |
---|
1115 | /// this has to be overloaded by derived workers |
---|
1116 | virtual bool pass(const PseudoJet & jet) const { |
---|
1117 | // make sure the centre is initialised |
---|
1118 | if (! _is_initialised) |
---|
1119 | throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1120 | |
---|
1121 | double distance2 = jet.squared_distance(_reference); |
---|
1122 | |
---|
1123 | return (distance2 <= _radius_out2) && (distance2 >= _radius_in2); |
---|
1124 | } |
---|
1125 | |
---|
1126 | /// returns a description of the worker |
---|
1127 | virtual string description() const { |
---|
1128 | ostringstream ostr; |
---|
1129 | ostr << sqrt(_radius_in2) << " <= distance from the centre <= " << sqrt(_radius_out2); |
---|
1130 | return ostr.str(); |
---|
1131 | } |
---|
1132 | |
---|
1133 | /// returns the rapidity range for which it may return "true" |
---|
1134 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{ |
---|
1135 | // make sure the centre is initialised |
---|
1136 | if (! _is_initialised) |
---|
1137 | throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1138 | |
---|
1139 | rapmax = _reference.rap()+sqrt(_radius_out2); |
---|
1140 | rapmin = _reference.rap()-sqrt(_radius_out2); |
---|
1141 | } |
---|
1142 | |
---|
1143 | virtual bool is_geometric() const { return true;} ///< implies a finite area |
---|
1144 | virtual bool has_finite_area() const { return true;} ///< regardless of the reference |
---|
1145 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
1146 | virtual double known_area() const { |
---|
1147 | return pi * (_radius_out2-_radius_in2); |
---|
1148 | } |
---|
1149 | |
---|
1150 | protected: |
---|
1151 | double _radius_in2, _radius_out2; |
---|
1152 | }; |
---|
1153 | |
---|
1154 | |
---|
1155 | |
---|
1156 | // select on objets with distance from the centre is between 'radius_in' and 'radius_out' |
---|
1157 | Selector SelectorDoughnut(const double & radius_in, const double & radius_out) { |
---|
1158 | return Selector(new SW_Doughnut(radius_in, radius_out)); |
---|
1159 | } |
---|
1160 | |
---|
1161 | |
---|
1162 | //---------------------------------------------------------------------- |
---|
1163 | /// helper for selecting on objects with rapidity within a distance 'delta' of a reference |
---|
1164 | class SW_Strip : public SW_WithReference { |
---|
1165 | public: |
---|
1166 | SW_Strip(const double &delta) : _delta(delta) {} |
---|
1167 | |
---|
1168 | /// return a copy of the current object |
---|
1169 | virtual SelectorWorker* copy(){ return new SW_Strip(*this);} |
---|
1170 | |
---|
1171 | /// returns true if a given object passes the selection criterium |
---|
1172 | /// this has to be overloaded by derived workers |
---|
1173 | virtual bool pass(const PseudoJet & jet) const { |
---|
1174 | // make sure the centre is initialised |
---|
1175 | if (! _is_initialised) |
---|
1176 | throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1177 | |
---|
1178 | return abs(jet.rap()-_reference.rap()) <= _delta; |
---|
1179 | } |
---|
1180 | |
---|
1181 | /// returns a description of the worker |
---|
1182 | virtual string description() const { |
---|
1183 | ostringstream ostr; |
---|
1184 | ostr << "|rap - rap_reference| <= " << _delta; |
---|
1185 | return ostr.str(); |
---|
1186 | } |
---|
1187 | |
---|
1188 | /// returns the rapidity range for which it may return "true" |
---|
1189 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{ |
---|
1190 | // make sure the centre is initialised |
---|
1191 | if (! _is_initialised) |
---|
1192 | throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1193 | |
---|
1194 | rapmax = _reference.rap()+_delta; |
---|
1195 | rapmin = _reference.rap()-_delta; |
---|
1196 | } |
---|
1197 | |
---|
1198 | virtual bool is_geometric() const { return true;} ///< implies a finite area |
---|
1199 | virtual bool has_finite_area() const { return true;} ///< regardless of the reference |
---|
1200 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
1201 | virtual double known_area() const { |
---|
1202 | return twopi * 2 * _delta; |
---|
1203 | } |
---|
1204 | |
---|
1205 | protected: |
---|
1206 | double _delta; |
---|
1207 | }; |
---|
1208 | |
---|
1209 | |
---|
1210 | // select on objets within a distance 'radius' of a variable location |
---|
1211 | Selector SelectorStrip(const double & half_width) { |
---|
1212 | return Selector(new SW_Strip(half_width)); |
---|
1213 | } |
---|
1214 | |
---|
1215 | |
---|
1216 | //---------------------------------------------------------------------- |
---|
1217 | /// helper for selecting on objects with rapidity within a distance |
---|
1218 | /// 'delta_rap' of a reference and phi within a distanve delta_phi of |
---|
1219 | /// a reference |
---|
1220 | class SW_Rectangle : public SW_WithReference { |
---|
1221 | public: |
---|
1222 | SW_Rectangle(const double &delta_rap, const double &delta_phi) |
---|
1223 | : _delta_rap(delta_rap), _delta_phi(delta_phi) {} |
---|
1224 | |
---|
1225 | /// return a copy of the current object |
---|
1226 | virtual SelectorWorker* copy(){ return new SW_Rectangle(*this);} |
---|
1227 | |
---|
1228 | /// returns true if a given object passes the selection criterium |
---|
1229 | /// this has to be overloaded by derived workers |
---|
1230 | virtual bool pass(const PseudoJet & jet) const { |
---|
1231 | // make sure the centre is initialised |
---|
1232 | if (! _is_initialised) |
---|
1233 | throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1234 | |
---|
1235 | return (abs(jet.rap()-_reference.rap()) <= _delta_rap) && (abs(jet.delta_phi_to(_reference)) <= _delta_phi); |
---|
1236 | } |
---|
1237 | |
---|
1238 | /// returns a description of the worker |
---|
1239 | virtual string description() const { |
---|
1240 | ostringstream ostr; |
---|
1241 | ostr << "|rap - rap_reference| <= " << _delta_rap << " && |phi - phi_reference| <= " << _delta_phi ; |
---|
1242 | return ostr.str(); |
---|
1243 | } |
---|
1244 | |
---|
1245 | /// returns the rapidity range for which it may return "true" |
---|
1246 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{ |
---|
1247 | // make sure the centre is initialised |
---|
1248 | if (! _is_initialised) |
---|
1249 | throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1250 | |
---|
1251 | rapmax = _reference.rap()+_delta_rap; |
---|
1252 | rapmin = _reference.rap()-_delta_rap; |
---|
1253 | } |
---|
1254 | |
---|
1255 | virtual bool is_geometric() const { return true;} ///< implies a finite area |
---|
1256 | virtual bool has_finite_area() const { return true;} ///< regardless of the reference |
---|
1257 | virtual bool has_known_area() const { return true;} ///< the area is analytically known |
---|
1258 | virtual double known_area() const { |
---|
1259 | return 4 * _delta_rap * _delta_phi; |
---|
1260 | } |
---|
1261 | |
---|
1262 | protected: |
---|
1263 | double _delta_rap, _delta_phi; |
---|
1264 | }; |
---|
1265 | |
---|
1266 | |
---|
1267 | // select on objets within a distance 'radius' of a variable location |
---|
1268 | Selector SelectorRectangle(const double & half_rap_width, const double & half_phi_width) { |
---|
1269 | return Selector(new SW_Rectangle(half_rap_width, half_phi_width)); |
---|
1270 | } |
---|
1271 | |
---|
1272 | |
---|
1273 | //---------------------------------------------------------------------- |
---|
1274 | /// helper for selecting the jets that carry at least a given fraction |
---|
1275 | /// of the reference jet |
---|
1276 | class SW_PtFractionMin : public SW_WithReference { |
---|
1277 | public: |
---|
1278 | /// ctor with specification of the number of objects to keep |
---|
1279 | SW_PtFractionMin(double fraction) : _fraction2(fraction*fraction){} |
---|
1280 | |
---|
1281 | /// return a copy of the current object |
---|
1282 | virtual SelectorWorker* copy(){ return new SW_PtFractionMin(*this);} |
---|
1283 | |
---|
1284 | /// return true if the jet carries a large enough fraction of the reference. |
---|
1285 | /// Throw an error if the reference is not initialised. |
---|
1286 | virtual bool pass(const PseudoJet & jet) const { |
---|
1287 | // make sure the centre is initialised |
---|
1288 | if (! _is_initialised) |
---|
1289 | throw Error("To use a SelectorPtFractionMin (or any selector that requires a reference), you first have to call set_reference(...)"); |
---|
1290 | |
---|
1291 | // otherwise, just call that method on the jet |
---|
1292 | return (jet.perp2() >= _fraction2*_reference.perp2()); |
---|
1293 | } |
---|
1294 | |
---|
1295 | /// returns a description of the worker |
---|
1296 | virtual string description() const { |
---|
1297 | ostringstream ostr; |
---|
1298 | ostr << "pt >= " << sqrt(_fraction2) << "* pt_ref"; |
---|
1299 | return ostr.str(); |
---|
1300 | } |
---|
1301 | |
---|
1302 | protected: |
---|
1303 | double _fraction2; |
---|
1304 | }; |
---|
1305 | |
---|
1306 | |
---|
1307 | // select objects that carry at least a fraction "fraction" of the reference jet |
---|
1308 | // (Note that this selectir takes a reference) |
---|
1309 | Selector SelectorPtFractionMin(double fraction){ |
---|
1310 | return Selector(new SW_PtFractionMin(fraction)); |
---|
1311 | } |
---|
1312 | |
---|
1313 | |
---|
1314 | //---------------------------------------------------------------------- |
---|
1315 | // additional (mostly helper) selectors |
---|
1316 | //---------------------------------------------------------------------- |
---|
1317 | |
---|
1318 | //---------------------------------------------------------------------- |
---|
1319 | /// helper for selecting the 0-momentum jets |
---|
1320 | class SW_IsZero : public SelectorWorker { |
---|
1321 | public: |
---|
1322 | /// ctor |
---|
1323 | SW_IsZero(){} |
---|
1324 | |
---|
1325 | /// return true if the jet has zero momentum |
---|
1326 | virtual bool pass(const PseudoJet & jet) const { |
---|
1327 | return jet==0; |
---|
1328 | } |
---|
1329 | |
---|
1330 | /// rereturns a description of the worker |
---|
1331 | virtual string description() const { return "zero";} |
---|
1332 | }; |
---|
1333 | |
---|
1334 | |
---|
1335 | // select objects with zero momentum |
---|
1336 | Selector SelectorIsZero(){ |
---|
1337 | return Selector(new SW_IsZero()); |
---|
1338 | } |
---|
1339 | |
---|
1340 | |
---|
1341 | //---------------------------------------------------------------------- |
---|
1342 | /// helper for selecting the pure ghost |
---|
1343 | class SW_IsPureGhost : public SelectorWorker { |
---|
1344 | public: |
---|
1345 | /// ctor |
---|
1346 | SW_IsPureGhost(){} |
---|
1347 | |
---|
1348 | /// return true if the jet is a pure-ghost jet |
---|
1349 | virtual bool pass(const PseudoJet & jet) const { |
---|
1350 | // if the jet has no area support then it's certainly not a ghost |
---|
1351 | if (!jet.has_area()) return false; |
---|
1352 | |
---|
1353 | // otherwise, just call that method on the jet |
---|
1354 | return jet.is_pure_ghost(); |
---|
1355 | } |
---|
1356 | |
---|
1357 | /// rereturns a description of the worker |
---|
1358 | virtual string description() const { return "pure ghost";} |
---|
1359 | }; |
---|
1360 | |
---|
1361 | |
---|
1362 | // select objects that are (or are only made of) ghosts |
---|
1363 | Selector SelectorIsPureGhost(){ |
---|
1364 | return Selector(new SW_IsPureGhost()); |
---|
1365 | } |
---|
1366 | |
---|
1367 | |
---|
1368 | //---------------------------------------------------------------------- |
---|
1369 | // Selector and workers for obtaining a Selector from an old |
---|
1370 | // RangeDefinition |
---|
1371 | // |
---|
1372 | // This is mostly intended for backward compatibility and is likely to |
---|
1373 | // be removed in a future major release of FastJet |
---|
1374 | //---------------------------------------------------------------------- |
---|
1375 | |
---|
1376 | //---------------------------------------------------------------------- |
---|
1377 | /// helper for selecting on both rapidity and azimuthal angle |
---|
1378 | class SW_RangeDefinition : public SelectorWorker{ |
---|
1379 | public: |
---|
1380 | /// ctor from a RangeDefinition |
---|
1381 | SW_RangeDefinition(const RangeDefinition &range) : _range(&range){} |
---|
1382 | |
---|
1383 | /// transfer the selection creterium to the underlying RangeDefinition |
---|
1384 | virtual bool pass(const PseudoJet & jet) const { |
---|
1385 | return _range->is_in_range(jet); |
---|
1386 | } |
---|
1387 | |
---|
1388 | /// returns a description of the worker |
---|
1389 | virtual string description() const { |
---|
1390 | return _range->description(); |
---|
1391 | } |
---|
1392 | |
---|
1393 | /// returns the rapidity range for which it may return "true" |
---|
1394 | virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{ |
---|
1395 | _range->get_rap_limits(rapmin, rapmax); |
---|
1396 | } |
---|
1397 | |
---|
1398 | /// check if it has a finite area |
---|
1399 | virtual bool is_geometric() const { return true;} |
---|
1400 | |
---|
1401 | /// check if it has an analytically computable area |
---|
1402 | virtual bool has_known_area() const { return true;} |
---|
1403 | |
---|
1404 | /// if it has a computable area, return it |
---|
1405 | virtual double known_area() const{ |
---|
1406 | return _range->area(); |
---|
1407 | } |
---|
1408 | |
---|
1409 | protected: |
---|
1410 | const RangeDefinition *_range; |
---|
1411 | }; |
---|
1412 | |
---|
1413 | |
---|
1414 | // ctor from a RangeDefinition |
---|
1415 | //---------------------------------------------------------------------- |
---|
1416 | // |
---|
1417 | // This is provided for backward compatibility and will be removed in |
---|
1418 | // a future major release of FastJet |
---|
1419 | Selector::Selector(const RangeDefinition &range) { |
---|
1420 | _worker.reset(new SW_RangeDefinition(range)); |
---|
1421 | } |
---|
1422 | |
---|
1423 | |
---|
1424 | // operators applying directly on a Selector |
---|
1425 | //---------------------------------------------------------------------- |
---|
1426 | |
---|
1427 | // operator &= |
---|
1428 | // For 2 Selectors a and b, a &= b is eauivalent to a = a & b; |
---|
1429 | Selector & Selector::operator &=(const Selector & b){ |
---|
1430 | _worker.reset(new SW_And(*this, b)); |
---|
1431 | return *this; |
---|
1432 | } |
---|
1433 | |
---|
1434 | // operator &= |
---|
1435 | // For 2 Selectors a and b, a &= b is eauivalent to a = a & b; |
---|
1436 | Selector & Selector::operator |=(const Selector & b){ |
---|
1437 | _worker.reset(new SW_Or(*this, b)); |
---|
1438 | return *this; |
---|
1439 | } |
---|
1440 | |
---|
1441 | FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh |
---|