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