FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Selector.cc
1 //STARTHEADER
2 // $Id: Selector.cc 2687 2011-11-14 11:17:51Z soyez $
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
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 &centre){
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 &centre){
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)
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
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
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