FastJet 3.0.2
Selector.cc
00001 //STARTHEADER
00002 // $Id: Selector.cc 2687 2011-11-14 11:17:51Z soyez $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 #include <sstream>
00031 #include <algorithm>
00032 #include "fastjet/Selector.hh"
00033 #include "fastjet/GhostedAreaSpec.hh"  // for area support
00034 
00035 using namespace std;
00036 
00037 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00038 
00039 //----------------------------------------------------------------------
00040 // implementations of some of the more complex bits of Selector
00041 //----------------------------------------------------------------------
00042 
00043 // implementation of the operator() acting on a vector of jets
00044 std::vector<PseudoJet> Selector::operator()(const std::vector<PseudoJet> & jets) const {
00045   std::vector<PseudoJet> result;
00046   const SelectorWorker * worker_local = validated_worker();
00047   if (worker_local->applies_jet_by_jet()) {
00048     //if (false) {
00049     // for workers that apply jet by jet, this is more efficient
00050     for (std::vector<PseudoJet>::const_iterator jet = jets.begin(); 
00051          jet != jets.end(); jet++) {
00052       if (worker_local->pass(*jet)) result.push_back(*jet);
00053     }
00054   } else {
00055     // for workers that can only be applied to entire vectors,
00056     // go through the following
00057     std::vector<const PseudoJet *> jetptrs(jets.size());
00058     for (unsigned i = 0; i < jets.size(); i++) {
00059       jetptrs[i] = & jets[i];
00060     }
00061     worker_local->terminator(jetptrs);
00062     for (unsigned i = 0; i < jetptrs.size(); i++) {
00063       if (jetptrs[i]) result.push_back(jets[i]);
00064     }
00065   }
00066   return result;
00067 }
00068 
00069 
00070 //----------------------------------------------------------------------
00071 // count the number of jets that pass the cuts
00072 unsigned int Selector::count(const std::vector<PseudoJet> & jets) const {
00073   unsigned n = 0;
00074   const SelectorWorker * worker_local = validated_worker();
00075   
00076   // separate strategies according to whether the worker applies jet by jet
00077   if (worker_local->applies_jet_by_jet()) {
00078     for (unsigned i = 0; i < jets.size(); i++) {
00079       if (worker_local->pass(jets[i])) n++;
00080     }
00081   } else {
00082     std::vector<const PseudoJet *> jetptrs(jets.size());
00083     for (unsigned i = 0; i < jets.size(); i++) {
00084       jetptrs[i] = & jets[i];
00085     }
00086     worker_local->terminator(jetptrs);
00087     for (unsigned i = 0; i < jetptrs.size(); i++) {
00088       if (jetptrs[i]) n++;
00089     }
00090   }
00091 
00092   return n;
00093 }
00094 
00095 
00096 //----------------------------------------------------------------------
00097 // sift the input jets into two vectors -- those that pass the selector
00098 // and those that do not
00099 void Selector::sift(const std::vector<PseudoJet> & jets,
00100                     std::vector<PseudoJet> & jets_that_pass,
00101                     std::vector<PseudoJet> & jets_that_fail
00102                     ) const {
00103   const SelectorWorker * worker_local = validated_worker();
00104   
00105   jets_that_pass.clear();
00106   jets_that_fail.clear();
00107   
00108   // separate strategies according to whether the worker applies jet by jet
00109   if (worker_local->applies_jet_by_jet()) {
00110     for (unsigned i = 0; i < jets.size(); i++) {
00111       if (worker_local->pass(jets[i])) {
00112         jets_that_pass.push_back(jets[i]);
00113       } else {
00114         jets_that_fail.push_back(jets[i]);
00115       }
00116     }
00117   } else {
00118     std::vector<const PseudoJet *> jetptrs(jets.size());
00119     for (unsigned i = 0; i < jets.size(); i++) {
00120       jetptrs[i] = & jets[i];
00121     }
00122     worker_local->terminator(jetptrs);
00123     for (unsigned i = 0; i < jetptrs.size(); i++) {
00124       if (jetptrs[i]) {
00125         jets_that_pass.push_back(jets[i]);
00126       } else {
00127         jets_that_fail.push_back(jets[i]);
00128       }
00129     }
00130   }
00131 }
00132 
00133 // area using default ghost area
00134 double Selector::area() const{
00135   return area(gas::def_ghost_area);
00136 }
00137 
00138 // implementation of the Selector's area function
00139 double Selector::area(double ghost_area) const{
00140   if (! is_geometric()) throw InvalidArea();
00141   
00142   // has area will already check we've got a valid worker
00143   if (_worker->has_known_area()) return _worker->known_area();
00144   
00145   // generate a set of "ghosts"
00146   double rapmin, rapmax;
00147   get_rapidity_extent(rapmin, rapmax);
00148   GhostedAreaSpec ghost_spec(rapmin, rapmax, 1, ghost_area);
00149   std::vector<PseudoJet> ghosts;
00150   ghost_spec.add_ghosts(ghosts);
00151   
00152   // check what passes the selection
00153   return ghost_spec.ghost_area() * ((*this)(ghosts)).size();
00154 }
00155 
00156 
00157 //----------------------------------------------------------------------
00158 // implementations of some of the more complex bits of SelectorWorker
00159 //----------------------------------------------------------------------
00160 // check if it has a finite area
00161 bool SelectorWorker::has_finite_area() const { 
00162   if (! is_geometric()) return false;
00163   double rapmin, rapmax;
00164   get_rapidity_extent(rapmin, rapmax);
00165   return (rapmax != std::numeric_limits<double>::infinity())
00166     &&  (-rapmin != std::numeric_limits<double>::infinity());
00167 }
00168 
00169 
00170 
00171 //----------------------------------------------------------------------
00172 // very basic set of selectors (at the moment just the identity!)
00173 //----------------------------------------------------------------------
00174 
00175 //----------------------------------------------------------------------
00176 /// helper for selecting the n hardest jets
00177 class SW_Identity : public SelectorWorker {
00178 public:
00179   /// ctor with specification of the number of objects to keep
00180   SW_Identity(){}
00181 
00182   /// just let everything pass
00183   virtual bool pass(const PseudoJet &) const {
00184     return true;
00185   }
00186 
00187   /// For each jet that does not pass the cuts, this routine sets the 
00188   /// pointer to 0. 
00189   virtual void terminator(vector<const PseudoJet *> &) const {
00190     // everything passes, hence nothing to nullify
00191     return;
00192   }
00193   
00194   /// returns a description of the worker
00195   virtual string description() const { return "Identity";}
00196 
00197   /// strictly speaking, this is geometric
00198   virtual bool is_geometric() const { return true;}
00199 };
00200 
00201 
00202 // returns an "identity" selector that lets everything pass
00203 Selector SelectorIdentity() {
00204   return Selector(new SW_Identity);
00205 }
00206 
00207 
00208 //----------------------------------------------------------------------
00209 // selector and workers for operators
00210 //----------------------------------------------------------------------
00211 
00212 //----------------------------------------------------------------------
00213 /// helper for combining selectors with a logical not
00214 class SW_Not : public SelectorWorker {
00215 public:
00216   /// ctor
00217   SW_Not(const Selector & s) : _s(s) {}
00218 
00219   /// return a copy of the current object
00220   virtual SelectorWorker* copy(){ return new SW_Not(*this);}
00221 
00222   /// returns true if a given object passes the selection criterium
00223   /// this has to be overloaded by derived workers
00224   virtual bool pass(const PseudoJet & jet) const {
00225     // make sure that the "pass" can be applied on both selectors
00226     if (!applies_jet_by_jet())
00227       throw Error("Cannot apply this selector worker to an individual jet");
00228     
00229     return ! _s.pass(jet);
00230   } 
00231 
00232   /// returns true if this can be applied jet by jet
00233   virtual bool applies_jet_by_jet() const {return _s.applies_jet_by_jet();}
00234 
00235   /// select the jets in the list that pass both selectors
00236   virtual void terminator(vector<const PseudoJet *> & jets) const {
00237     // if we can apply the selector jet-by-jet, call the base selector
00238     // that does exactly that
00239     if (applies_jet_by_jet()){
00240       SelectorWorker::terminator(jets);
00241       return;
00242     }
00243 
00244     // check the effect of the selector we want to negate
00245     vector<const PseudoJet *> s_jets = jets;
00246     _s.worker()->terminator(s_jets);
00247 
00248     // now apply the negation: all the jets that pass the base
00249     // selector (i.e. are not NULL) have to be set to NULL
00250     for (unsigned int i=0; i<s_jets.size(); i++){
00251       if (s_jets[i]) jets[i] = NULL;
00252     }
00253   }
00254 
00255   /// returns a description of the worker
00256   virtual string description() const {
00257     ostringstream ostr;
00258     ostr << "!(" << _s.description() << ")";
00259     return ostr.str();
00260   }
00261 
00262   /// is geometric if the underlying selector is
00263   virtual bool is_geometric() const { return _s.is_geometric();}
00264 
00265   /// returns true if the worker can be set_referenced
00266   virtual bool takes_reference() const { return _s.takes_reference();}
00267 
00268   /// set the reference jet for this selector
00269   virtual void set_reference(const PseudoJet &ref) { _s.set_reference(ref);}
00270 
00271 protected:
00272   Selector _s;
00273 };
00274 
00275 
00276 // logical not applied on a selector
00277 Selector operator!(const Selector & s) {
00278   return Selector(new SW_Not(s));
00279 }
00280 
00281 
00282 //----------------------------------------------------------------------
00283 /// Base class for binary operators
00284 class SW_BinaryOperator: public SelectorWorker {
00285 public:
00286   /// ctor
00287   SW_BinaryOperator(const Selector & s1, const Selector & s2) : _s1(s1), _s2(s2) {
00288     // stores info for more efficient access to the selector's properties
00289 
00290     // we can apply jet by jet only if this is the case for both sub-selectors
00291     _applies_jet_by_jet = _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
00292 
00293     // the selector takes a reference if either of the sub-selectors does
00294     _takes_reference = _s1.takes_reference() || _s2.takes_reference();
00295 
00296     // we have a well-defined area provided the two objects have one
00297     _is_geometric = _s1.is_geometric() && _s2.is_geometric();
00298   }
00299 
00300   /// returns true if this can be applied jet by jet
00301   virtual bool applies_jet_by_jet() const {return _applies_jet_by_jet;}
00302 
00303   /// returns true if this takes a reference jet
00304   virtual bool takes_reference() const{ 
00305     return _takes_reference;
00306   }
00307 
00308   /// sets the reference jet
00309   virtual void set_reference(const PseudoJet &centre){
00310     _s1.set_reference(centre);
00311     _s2.set_reference(centre);
00312   }
00313 
00314   /// check if it has a finite area
00315   virtual bool is_geometric() const { return _is_geometric;} 
00316 
00317 protected:
00318   Selector _s1, _s2;
00319   bool _applies_jet_by_jet;
00320   bool _takes_reference;
00321   bool _is_geometric;
00322 };
00323 
00324 
00325 
00326 //----------------------------------------------------------------------
00327 /// helper for combining selectors with a logical and
00328 class SW_And: public SW_BinaryOperator {
00329 public:
00330   /// ctor
00331   SW_And(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2){}
00332 
00333   /// return a copy of this
00334   virtual SelectorWorker* copy(){ return new SW_And(*this);}
00335 
00336   /// returns true if a given object passes the selection criterium
00337   /// this has to be overloaded by derived workers
00338   virtual bool pass(const PseudoJet & jet) const {
00339     // make sure that the "pass" can be applied on both selectors
00340     if (!applies_jet_by_jet())
00341       throw Error("Cannot apply this selector worker to an individual jet");
00342     
00343     return _s1.pass(jet) && _s2.pass(jet);
00344   }
00345 
00346   /// select the jets in the list that pass both selectors
00347   virtual void terminator(vector<const PseudoJet *> & jets) const {
00348     // if we can apply the selector jet-by-jet, call the base selector
00349     // that does exactly that
00350     if (applies_jet_by_jet()){
00351       SelectorWorker::terminator(jets);
00352       return;
00353     }
00354 
00355     // check the effect of the first selector
00356     vector<const PseudoJet *> s1_jets = jets;
00357     _s1.worker()->terminator(s1_jets);
00358 
00359     // apply the second
00360     _s2.worker()->terminator(jets);
00361 
00362     // terminate the jets that wiuld be terminated by _s1
00363     for (unsigned int i=0; i<jets.size(); i++){
00364       if (! s1_jets[i]) jets[i] = NULL;
00365     }
00366   }
00367 
00368   /// returns the rapidity range for which it may return "true"
00369   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
00370     double s1min, s1max, s2min, s2max;
00371     _s1.get_rapidity_extent(s1min, s1max);
00372     _s2.get_rapidity_extent(s2min, s2max);
00373     rapmax = min(s1max, s2max);
00374     rapmin = max(s1min, s2min);
00375   }
00376 
00377   /// returns a description of the worker
00378   virtual string description() const {
00379     ostringstream ostr;
00380     ostr << "(" << _s1.description() << " && " << _s2.description() << ")";
00381     return ostr.str();
00382   }
00383 };
00384 
00385 
00386 // logical and between two selectors
00387 Selector operator&&(const Selector & s1, const Selector & s2) {
00388   return Selector(new SW_And(s1,s2));
00389 }
00390 
00391 
00392 
00393 //----------------------------------------------------------------------
00394 /// helper for combining selectors with a logical or
00395 class SW_Or: public SW_BinaryOperator {
00396 public:
00397   /// ctor
00398   SW_Or(const Selector & s1, const Selector & s2) : SW_BinaryOperator(s1,s2) {}
00399 
00400   /// return a copy of this
00401   virtual SelectorWorker* copy(){ return new SW_Or(*this);}
00402 
00403   /// returns true if a given object passes the selection criterium
00404   /// this has to be overloaded by derived workers
00405   virtual bool pass(const PseudoJet & jet) const {
00406     // make sure that the "pass" can be applied on both selectors
00407     if (!applies_jet_by_jet())
00408       throw Error("Cannot apply this selector worker to an individual jet");
00409     
00410     return _s1.pass(jet) || _s2.pass(jet);
00411   }
00412 
00413   /// returns true if this can be applied jet by jet
00414   virtual bool applies_jet_by_jet() const {
00415     // watch out, even though it's the "OR" selector, to be applied jet
00416     // by jet, both the base selectors need to be jet-by-jet-applicable,
00417     // so the use of a && in the line below
00418     return _s1.applies_jet_by_jet() && _s2.applies_jet_by_jet();
00419   }
00420 
00421   /// select the jets in the list that pass both selectors
00422   virtual void terminator(vector<const PseudoJet *> & jets) const {
00423     // if we can apply the selector jet-by-jet, call the base selector
00424     // that does exactly that
00425     if (applies_jet_by_jet()){
00426       SelectorWorker::terminator(jets);
00427       return;
00428     }
00429 
00430     // check the effect of the first selector
00431     vector<const PseudoJet *> s1_jets = jets;
00432     _s1.worker()->terminator(s1_jets);
00433 
00434     // apply the second
00435     _s2.worker()->terminator(jets);
00436 
00437     // resurrect any jet that has been terminated by the second one
00438     // and not by the first one
00439     for (unsigned int i=0; i<jets.size(); i++){
00440       if (s1_jets[i]) jets[i] = s1_jets[i];
00441     }
00442   }
00443 
00444   /// returns a description of the worker
00445   virtual string description() const {
00446     ostringstream ostr;
00447     ostr << "(" << _s1.description() << " || " << _s2.description() << ")";
00448     return ostr.str();
00449   }
00450 
00451   /// returns the rapidity range for which it may return "true"
00452   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const {
00453     double s1min, s1max, s2min, s2max;
00454     _s1.get_rapidity_extent(s1min, s1max);
00455     _s2.get_rapidity_extent(s2min, s2max);
00456     rapmax = max(s1max, s2max);
00457     rapmin = min(s1min, s2min);
00458   }
00459 };
00460 
00461 
00462 // logical or between two selectors
00463 Selector operator ||(const Selector & s1, const Selector & s2) {
00464   return Selector(new SW_Or(s1,s2));
00465 }
00466 
00467 //----------------------------------------------------------------------
00468 /// helper for multiplying two selectors (in an operator-like way)
00469 class SW_Mult: public SW_And {
00470 public:
00471   /// ctor
00472   SW_Mult(const Selector & s1, const Selector & s2) : SW_And(s1,s2) {}
00473 
00474   /// return a copy of this
00475   virtual SelectorWorker* copy(){ return new SW_Mult(*this);}
00476 
00477   /// select the jets in the list that pass both selectors
00478   virtual void terminator(vector<const PseudoJet *> & jets) const {
00479     // if we can apply the selector jet-by-jet, call the base selector
00480     // that does exactly that
00481     if (applies_jet_by_jet()){
00482       SelectorWorker::terminator(jets);
00483       return;
00484     }
00485 
00486     // first apply _s2
00487     _s2.worker()->terminator(jets);
00488 
00489     // then apply _s1
00490     _s1.worker()->terminator(jets);
00491   }
00492 
00493   /// returns a description of the worker
00494   virtual string description() const {
00495     ostringstream ostr;
00496     ostr << "(" << _s1.description() << " * " << _s2.description() << ")";
00497     return ostr.str();
00498   }
00499 };
00500 
00501 
00502 // logical and between two selectors
00503 Selector operator*(const Selector & s1, const Selector & s2) {
00504   return Selector(new SW_Mult(s1,s2));
00505 }
00506 
00507 
00508 //----------------------------------------------------------------------
00509 // selector and workers for kinematic cuts
00510 //----------------------------------------------------------------------
00511 
00512 //----------------------------------------------------------------------
00513 // a series of basic classes that allow easy implementations of
00514 // min, max and ranges on a quantity-to-be-defined
00515 
00516 // generic holder for a quantity
00517 class QuantityBase{
00518 public:
00519   QuantityBase(double q) : _q(q){}
00520   virtual ~QuantityBase(){}
00521   virtual double operator()(const PseudoJet & jet ) const =0;
00522   virtual string description() const =0;
00523   virtual bool is_geometric() const { return false;}
00524   virtual double comparison_value() const {return _q;}
00525   virtual double description_value() const {return comparison_value();}
00526 protected:
00527   double _q;
00528 };  
00529 
00530 // generic holder for a squared quantity
00531 class QuantitySquareBase : public QuantityBase{
00532 public:
00533   QuantitySquareBase(double sqrtq) : QuantityBase(sqrtq*sqrtq), _sqrtq(sqrtq){}
00534   virtual double description_value() const {return _sqrtq;}
00535 protected:
00536   double _sqrtq;
00537 };  
00538 
00539 // generic_quantity >= minimum
00540 template<typename QuantityType>
00541 class SW_QuantityMin : public SelectorWorker{
00542 public:
00543   /// detfault ctor (initialises the pt cut)
00544   SW_QuantityMin(double qmin) : _qmin(qmin) {}
00545 
00546   /// returns true is the given object passes the selection pt cut
00547   virtual bool pass(const PseudoJet & jet) const {return _qmin(jet) >= _qmin.comparison_value();}
00548 
00549   /// returns a description of the worker
00550   virtual string description() const {
00551     ostringstream ostr;
00552     ostr << _qmin.description() << " >= " << _qmin.description_value();
00553     return ostr.str();
00554   }
00555 
00556   virtual bool is_geometric() const { return _qmin.is_geometric();}
00557 
00558 protected:
00559   QuantityType _qmin;     ///< the cut
00560 };
00561 
00562 
00563 // generic_quantity <= maximum
00564 template<typename QuantityType>
00565 class SW_QuantityMax : public SelectorWorker {
00566 public:
00567   /// detfault ctor (initialises the pt cut)
00568   SW_QuantityMax(double qmax) : _qmax(qmax) {}
00569 
00570   /// returns true is the given object passes the selection pt cut
00571   virtual bool pass(const PseudoJet & jet) const {return _qmax(jet) <= _qmax.comparison_value();}
00572 
00573   /// returns a description of the worker
00574   virtual string description() const {
00575     ostringstream ostr;
00576     ostr << _qmax.description() << " <= " << _qmax.description_value();
00577     return ostr.str();
00578   }
00579 
00580   virtual bool is_geometric() const { return _qmax.is_geometric();}
00581 
00582 protected:
00583   QuantityType _qmax;   ///< the cut
00584 };
00585 
00586 
00587 // generic quantity in [minimum:maximum]
00588 template<typename QuantityType>
00589 class SW_QuantityRange : public SelectorWorker {
00590 public:
00591   /// detfault ctor (initialises the pt cut)
00592   SW_QuantityRange(double qmin, double qmax) : _qmin(qmin), _qmax(qmax) {}
00593 
00594   /// returns true is the given object passes the selection pt cut
00595   virtual bool pass(const PseudoJet & jet) const {
00596     double q = _qmin(jet); // we could identically use _qmax
00597     return (q >= _qmin.comparison_value()) && (q <= _qmax.comparison_value());
00598   }
00599 
00600   /// returns a description of the worker
00601   virtual string description() const {
00602     ostringstream ostr;
00603     ostr << _qmin.description_value() << " <= " << _qmin.description() << " <= " << _qmax.description_value();
00604     return ostr.str();
00605   }
00606 
00607   virtual bool is_geometric() const { return _qmin.is_geometric();}
00608 
00609 protected:
00610   QuantityType _qmin;   // the lower cut 
00611   QuantityType _qmax;   // the upper cut
00612 };
00613 
00614 
00615 //----------------------------------------------------------------------
00616 /// helper class for selecting on pt
00617 class QuantityPt2 : public QuantitySquareBase{
00618 public:
00619   QuantityPt2(double pt) : QuantitySquareBase(pt){}
00620   virtual double operator()(const PseudoJet & jet ) const { return jet.perp2();}
00621   virtual string description() const {return "pt";}
00622 };  
00623 
00624 // returns a selector for a minimum pt
00625 Selector SelectorPtMin(double ptmin) {
00626   return Selector(new SW_QuantityMin<QuantityPt2>(ptmin));
00627 }
00628 
00629 // returns a selector for a maximum pt
00630 Selector SelectorPtMax(double ptmax) {
00631   return Selector(new SW_QuantityMax<QuantityPt2>(ptmax));
00632 }
00633 
00634 // returns a selector for a pt range
00635 Selector SelectorPtRange(double ptmin, double ptmax) {
00636   return Selector(new SW_QuantityRange<QuantityPt2>(ptmin, ptmax));
00637 }
00638 
00639 
00640 //----------------------------------------------------------------------
00641 /// helper class for selecting on transverse energy
00642 class QuantityEt2 : public QuantitySquareBase{
00643 public:
00644   QuantityEt2(double Et) : QuantitySquareBase(Et){}
00645   virtual double operator()(const PseudoJet & jet ) const { return jet.Et2();}
00646   virtual string description() const {return "Et";}
00647 };  
00648 
00649 // returns a selector for a minimum Et
00650 Selector SelectorEtMin(double Etmin) {
00651   return Selector(new SW_QuantityMin<QuantityEt2>(Etmin));
00652 }
00653 
00654 // returns a selector for a maximum Et
00655 Selector SelectorEtMax(double Etmax) {
00656   return Selector(new SW_QuantityMax<QuantityEt2>(Etmax));
00657 }
00658 
00659 // returns a selector for a Et range
00660 Selector SelectorEtRange(double Etmin, double Etmax) {
00661   return Selector(new SW_QuantityRange<QuantityEt2>(Etmin, Etmax));
00662 }
00663 
00664 
00665 //----------------------------------------------------------------------
00666 /// helper class for selecting on energy
00667 class QuantityE : public QuantityBase{
00668 public:
00669   QuantityE(double E) : QuantityBase(E){}
00670   virtual double operator()(const PseudoJet & jet ) const { return jet.E();}
00671   virtual string description() const {return "E";}
00672 };  
00673 
00674 // returns a selector for a minimum E
00675 Selector SelectorEMin(double Emin) {
00676   return Selector(new SW_QuantityMin<QuantityE>(Emin));
00677 }
00678 
00679 // returns a selector for a maximum E
00680 Selector SelectorEMax(double Emax) {
00681   return Selector(new SW_QuantityMax<QuantityE>(Emax));
00682 }
00683 
00684 // returns a selector for a E range
00685 Selector SelectorERange(double Emin, double Emax) {
00686   return Selector(new SW_QuantityRange<QuantityE>(Emin, Emax));
00687 }
00688 
00689 
00690 //----------------------------------------------------------------------
00691 /// helper class for selecting on mass
00692 class QuantityM2 : public QuantitySquareBase{
00693 public:
00694   QuantityM2(double m) : QuantitySquareBase(m){}
00695   virtual double operator()(const PseudoJet & jet ) const { return jet.m2();}
00696   virtual string description() const {return "mass";}
00697 };  
00698 
00699 // returns a selector for a minimum mass
00700 Selector SelectorMassMin(double mmin) {
00701   return Selector(new SW_QuantityMin<QuantityM2>(mmin));
00702 }
00703 
00704 // returns a selector for a maximum mass
00705 Selector SelectorMassMax(double mmax) {
00706   return Selector(new SW_QuantityMax<QuantityM2>(mmax));
00707 }
00708 
00709 // returns a selector for a mass range
00710 Selector SelectorMassRange(double mmin, double mmax) {
00711   return Selector(new SW_QuantityRange<QuantityM2>(mmin, mmax));
00712 }
00713 
00714 
00715 
00716 //----------------------------------------------------------------------
00717 /// helper for selecting on rapidities: quantity
00718 class QuantityRap : public QuantityBase{
00719 public:
00720   QuantityRap(double rap) : QuantityBase(rap){}
00721   virtual double operator()(const PseudoJet & jet ) const { return jet.rap();}
00722   virtual string description() const {return "rap";}
00723   virtual bool is_geometric() const { return true;}
00724 };  
00725 
00726 
00727 /// helper for selecting on rapidities: min
00728 class SW_RapMin : public SW_QuantityMin<QuantityRap>{
00729 public:
00730   SW_RapMin(double rapmin) : SW_QuantityMin<QuantityRap>(rapmin){}
00731   virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
00732     rapmax = std::numeric_limits<double>::max();     
00733     rapmin = _qmin.comparison_value();
00734   }
00735 };
00736 
00737 /// helper for selecting on rapidities: max
00738 class SW_RapMax : public SW_QuantityMax<QuantityRap>{
00739 public:
00740   SW_RapMax(double rapmax) : SW_QuantityMax<QuantityRap>(rapmax){}
00741   virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
00742     rapmax = _qmax.comparison_value(); 
00743     rapmin = -std::numeric_limits<double>::max();
00744   }
00745 };
00746 
00747 /// helper for selecting on rapidities: range
00748 class SW_RapRange : public SW_QuantityRange<QuantityRap>{
00749 public:
00750   SW_RapRange(double rapmin, double rapmax) : SW_QuantityRange<QuantityRap>(rapmin, rapmax){
00751     assert(rapmin<=rapmax);
00752   }
00753   virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
00754     rapmax = _qmax.comparison_value();      
00755     rapmin = _qmin.comparison_value(); 
00756   }
00757   virtual bool has_known_area() const { return true;} ///< the area is analytically known
00758   virtual double known_area() const { 
00759     return twopi * (_qmax.comparison_value()-_qmin.comparison_value());
00760   }
00761 };
00762 
00763 // returns a selector for a minimum rapidity
00764 Selector SelectorRapMin(double rapmin) {
00765   return Selector(new SW_RapMin(rapmin));
00766 }
00767 
00768 // returns a selector for a maximum rapidity
00769 Selector SelectorRapMax(double rapmax) {
00770   return Selector(new SW_RapMax(rapmax));
00771 }
00772 
00773 // returns a selector for a rapidity range
00774 Selector SelectorRapRange(double rapmin, double rapmax) {
00775   return Selector(new SW_RapRange(rapmin, rapmax));
00776 }
00777 
00778 
00779 //----------------------------------------------------------------------
00780 /// helper for selecting on |rapidities|
00781 class QuantityAbsRap : public QuantityBase{
00782 public:
00783   QuantityAbsRap(double absrap) : QuantityBase(absrap){}
00784   virtual double operator()(const PseudoJet & jet ) const { return abs(jet.rap());}
00785   virtual string description() const {return "|rap|";}
00786   virtual bool is_geometric() const { return true;}
00787 };  
00788 
00789 
00790 /// helper for selecting on |rapidities|: max
00791 class SW_AbsRapMax : public SW_QuantityMax<QuantityAbsRap>{
00792 public:
00793   SW_AbsRapMax(double absrapmax) : SW_QuantityMax<QuantityAbsRap>(absrapmax){}
00794   virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
00795     rapmax =  _qmax.comparison_value(); 
00796     rapmin = -_qmax.comparison_value();
00797   }
00798   virtual bool has_known_area() const { return true;}   ///< the area is analytically known
00799   virtual double known_area() const { 
00800     return twopi * 2 * _qmax.comparison_value();
00801   }
00802 };
00803 
00804 /// helper for selecting on |rapidities|: max
00805 class SW_AbsRapRange : public SW_QuantityRange<QuantityAbsRap>{
00806 public:
00807   SW_AbsRapRange(double absrapmin, double absrapmax) : SW_QuantityRange<QuantityAbsRap>(absrapmin, absrapmax){}
00808   virtual void get_rapidity_extent(double &rapmin, double & rapmax) const{
00809     rapmax =  _qmax.comparison_value(); 
00810     rapmin = -_qmax.comparison_value();
00811   }
00812   virtual bool has_known_area() const { return true;} ///< the area is analytically known
00813   virtual double known_area() const { 
00814     return twopi * 2 * (_qmax.comparison_value()-max(_qmin.comparison_value(),0.0)); // this should handle properly absrapmin<0
00815   }
00816 };
00817 
00818 // returns a selector for a minimum |rapidity|
00819 Selector SelectorAbsRapMin(double absrapmin) {
00820   return Selector(new SW_QuantityMin<QuantityAbsRap>(absrapmin));
00821 }
00822 
00823 // returns a selector for a maximum |rapidity|
00824 Selector SelectorAbsRapMax(double absrapmax) {
00825   return Selector(new SW_AbsRapMax(absrapmax));
00826 }
00827 
00828 // returns a selector for a |rapidity| range
00829 Selector SelectorAbsRapRange(double rapmin, double rapmax) {
00830   return Selector(new SW_AbsRapRange(rapmin, rapmax));
00831 }
00832 
00833 
00834 //----------------------------------------------------------------------
00835 /// helper for selecting on pseudo-rapidities
00836 class QuantityEta : public QuantityBase{
00837 public:
00838   QuantityEta(double eta) : QuantityBase(eta){}
00839   virtual double operator()(const PseudoJet & jet ) const { return jet.eta();}
00840   virtual string description() const {return "eta";}
00841   // virtual bool is_geometric() const { return true;} // not strictly only y and phi-dependent
00842 };  
00843 
00844 // returns a selector for a pseudo-minimum rapidity
00845 Selector SelectorEtaMin(double etamin) {
00846   return Selector(new SW_QuantityMin<QuantityEta>(etamin));
00847 }
00848 
00849 // returns a selector for a pseudo-maximum rapidity
00850 Selector SelectorEtaMax(double etamax) {
00851   return Selector(new SW_QuantityMax<QuantityEta>(etamax));
00852 }
00853 
00854 // returns a selector for a pseudo-rapidity range
00855 Selector SelectorEtaRange(double etamin, double etamax) {
00856   return Selector(new SW_QuantityRange<QuantityEta>(etamin, etamax));
00857 }
00858 
00859 
00860 //----------------------------------------------------------------------
00861 /// helper for selecting on |pseudo-rapidities|
00862 class QuantityAbsEta : public QuantityBase{
00863 public:
00864   QuantityAbsEta(double abseta) : QuantityBase(abseta){}
00865   virtual double operator()(const PseudoJet & jet ) const { return abs(jet.eta());}
00866   virtual string description() const {return "|eta|";}
00867   virtual bool is_geometric() const { return true;}
00868 };  
00869 
00870 // returns a selector for a minimum |pseudo-rapidity|
00871 Selector SelectorAbsEtaMin(double absetamin) {
00872   return Selector(new SW_QuantityMin<QuantityAbsEta>(absetamin));
00873 }
00874 
00875 // returns a selector for a maximum |pseudo-rapidity|
00876 Selector SelectorAbsEtaMax(double absetamax) {
00877   return Selector(new SW_QuantityMax<QuantityAbsEta>(absetamax));
00878 }
00879 
00880 // returns a selector for a |pseudo-rapidity| range
00881 Selector SelectorAbsEtaRange(double absetamin, double absetamax) {
00882   return Selector(new SW_QuantityRange<QuantityAbsEta>(absetamin, absetamax));
00883 }
00884 
00885 
00886 //----------------------------------------------------------------------
00887 /// helper for selecting on azimuthal angle
00888 ///
00889 /// Note that the bounds have to be specified as min<max
00890 /// phimin has to be > -2pi
00891 /// phimax has to be <  4pi
00892 class SW_PhiRange : public SelectorWorker {
00893 public:
00894   /// detfault ctor (initialises the pt cut)
00895   SW_PhiRange(double phimin, double phimax) : _phimin(phimin), _phimax(phimax){
00896     assert(_phimin<_phimax);
00897     assert(_phimin>-twopi);
00898     assert(_phimax<2*twopi);
00899 
00900     _phispan = _phimax - _phimin;
00901   }
00902 
00903   /// returns true is the given object passes the selection pt cut
00904   virtual bool pass(const PseudoJet & jet) const {
00905     double dphi=jet.phi()-_phimin;
00906     if (dphi >= twopi) dphi -= twopi;
00907     if (dphi < 0)      dphi += twopi;
00908     return (dphi <= _phispan);
00909   }
00910 
00911   /// returns a description of the worker
00912   virtual string description() const {
00913     ostringstream ostr;
00914     ostr << _phimin << " <= phi <= " << _phimax;
00915     return ostr.str();
00916   }
00917 
00918   virtual bool is_geometric() const { return true;}
00919 
00920 protected:
00921   double _phimin;   // the lower cut 
00922   double _phimax;   // the upper cut
00923   double _phispan;  // the span of the range
00924 };
00925 
00926 
00927 // returns a selector for a phi range
00928 Selector SelectorPhiRange(double phimin, double phimax) {
00929   return Selector(new SW_PhiRange(phimin, phimax));
00930 }
00931 
00932 //----------------------------------------------------------------------
00933 /// helper for selecting on both rapidity and azimuthal angle
00934 class SW_RapPhiRange : public SW_And{
00935 public:
00936   SW_RapPhiRange(double rapmin, double rapmax, double phimin, double phimax)
00937     : SW_And(SelectorRapRange(rapmin, rapmax), SelectorPhiRange(phimin, phimax)){
00938     _known_area = ((phimax-phimin > twopi) ? twopi : phimax-phimin) * (rapmax-rapmin);
00939   }
00940 
00941   /// if it has a computable area, return it
00942   virtual double known_area() const{
00943     return _known_area;
00944   }
00945 
00946 protected:
00947   double _known_area;
00948 };
00949 
00950 Selector SelectorRapPhiRange(double rapmin, double rapmax, double phimin, double phimax) {
00951   return Selector(new SW_RapPhiRange(rapmin, rapmax, phimin, phimax));
00952 }
00953 
00954 
00955 //----------------------------------------------------------------------
00956 /// helper for selecting the n hardest jets
00957 class SW_NHardest : public SelectorWorker {
00958 public:
00959   /// ctor with specification of the number of objects to keep
00960   SW_NHardest(unsigned int n) : _n(n) {};
00961 
00962   /// pass makes no sense here normally the parent selector will throw
00963   /// an error but for internal use in the SW, we'll throw one from
00964   /// here by security
00965   virtual bool pass(const PseudoJet &) const {
00966     if (!applies_jet_by_jet())
00967       throw Error("Cannot apply this selector worker to an individual jet");
00968     return false;
00969   }
00970 
00971   /// For each jet that does not pass the cuts, this routine sets the 
00972   /// pointer to 0. 
00973   virtual void terminator(vector<const PseudoJet *> & jets) const {
00974     // nothing to do if the size is too small
00975     if (jets.size() < _n) return;
00976 
00977     // do we want to first chech if things are already ordered before
00978     // going through the ordering process? For now, no. Maybe carry
00979     // out timing tests at some point to establish the optimal
00980     // strategy.
00981 
00982     vector<double> minus_pt2(jets.size());
00983     vector<unsigned int> indices(jets.size());
00984 
00985     for (unsigned int i=0; i<jets.size(); i++){
00986       indices[i] = i;
00987 
00988       // we need to make sure that the object has not already been
00989       // nullified.  Note that if we have less than _n jets, this
00990       // whole n-hardest selection will not have any effect.
00991       minus_pt2[i] = jets[i] ? -jets[i]->perp2() : 0.0;
00992     }
00993     
00994     IndexedSortHelper sort_helper(& minus_pt2);
00995     
00996     partial_sort(indices.begin(), indices.begin()+_n, indices.end(), sort_helper);
00997     
00998     for (unsigned int i=_n; i<jets.size(); i++)
00999       jets[indices[i]] = NULL;
01000   }
01001   
01002   /// returns true if this can be applied jet by jet
01003   virtual bool applies_jet_by_jet() const {return false;}
01004   
01005   /// returns a description of the worker
01006   virtual string description() const {
01007     ostringstream ostr;
01008     ostr << _n << " hardest";
01009     return ostr.str();
01010   }
01011   
01012 protected:
01013   unsigned int _n;
01014 };
01015 
01016 
01017 // returns a selector for the n hardest jets
01018 Selector SelectorNHardest(unsigned int n) {
01019   return Selector(new SW_NHardest(n));
01020 }
01021 
01022 
01023 
01024 //----------------------------------------------------------------------
01025 // selector and workers for geometric ranges
01026 //----------------------------------------------------------------------
01027 
01028 //----------------------------------------------------------------------
01029 /// a generic class for objects that contain a position
01030 class SW_WithReference : public SelectorWorker{
01031 public:
01032   /// ctor
01033   SW_WithReference() : _is_initialised(false){};
01034 
01035   /// returns true if the worker takes a reference jet
01036   virtual bool takes_reference() const { return true;}
01037 
01038   /// sets the reference jet
01039   virtual void set_reference(const PseudoJet &centre){
01040     _is_initialised = true;
01041     _reference = centre;
01042   }
01043 
01044 protected:
01045   PseudoJet _reference;
01046   bool _is_initialised;
01047 };
01048 
01049 //----------------------------------------------------------------------
01050 /// helper for selecting on objects within a distance 'radius' of a reference
01051 class SW_Circle : public SW_WithReference {
01052 public:
01053   SW_Circle(const double &radius) : _radius2(radius*radius) {}
01054 
01055   /// return a copy of the current object
01056   virtual SelectorWorker* copy(){ return new SW_Circle(*this);}
01057 
01058   /// returns true if a given object passes the selection criterium
01059   /// this has to be overloaded by derived workers
01060   virtual bool pass(const PseudoJet & jet) const {
01061     // make sure the centre is initialised
01062     if (! _is_initialised)
01063       throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
01064     
01065     return jet.squared_distance(_reference) <= _radius2;
01066   } 
01067 
01068   /// returns a description of the worker
01069   virtual string description() const {
01070     ostringstream ostr;
01071     ostr << "distance from the centre <= " << sqrt(_radius2);
01072     return ostr.str();
01073   }
01074 
01075   /// returns the rapidity range for which it may return "true"
01076   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
01077     // make sure the centre is initialised
01078     if (! _is_initialised)
01079       throw Error("To use a SelectorCircle (or any selector that requires a reference), you first have to call set_reference(...)");
01080     
01081     rapmax = _reference.rap()+sqrt(_radius2);
01082     rapmin = _reference.rap()-sqrt(_radius2);
01083   }
01084 
01085   virtual bool is_geometric() const { return true;}    ///< implies a finite area
01086   virtual bool has_finite_area() const { return true;} ///< regardless of the reference 
01087   virtual bool has_known_area() const { return true;}  ///< the area is analytically known
01088   virtual double known_area() const { 
01089     return pi * _radius2;
01090   }
01091 
01092 protected:
01093   double _radius2;
01094 };
01095 
01096 
01097 // select on objets within a distance 'radius' of a variable location
01098 Selector SelectorCircle(const double & radius) {
01099   return Selector(new SW_Circle(radius));
01100 }
01101 
01102 
01103 //----------------------------------------------------------------------
01104 /// helper for selecting on objects with a distance to a reference
01105 /// betwene 'radius_in' and 'radius_out'
01106 class SW_Doughnut : public SW_WithReference {
01107 public:
01108   SW_Doughnut(const double &radius_in, const double &radius_out)
01109     : _radius_in2(radius_in*radius_in), _radius_out2(radius_out*radius_out) {}
01110 
01111   /// return a copy of the current object
01112   virtual SelectorWorker* copy(){ return new SW_Doughnut(*this);}
01113 
01114   /// returns true if a given object passes the selection criterium
01115   /// this has to be overloaded by derived workers
01116   virtual bool pass(const PseudoJet & jet) const {
01117     // make sure the centre is initialised
01118     if (! _is_initialised)
01119       throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
01120 
01121     double distance2 = jet.squared_distance(_reference);
01122 
01123     return (distance2 <= _radius_out2) && (distance2 >= _radius_in2);
01124   } 
01125 
01126   /// returns a description of the worker
01127   virtual string description() const {
01128     ostringstream ostr;
01129     ostr << sqrt(_radius_in2) << " <= distance from the centre <= " << sqrt(_radius_out2);
01130     return ostr.str();
01131   }
01132 
01133   /// returns the rapidity range for which it may return "true"
01134   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
01135     // make sure the centre is initialised
01136     if (! _is_initialised)
01137       throw Error("To use a SelectorDoughnut (or any selector that requires a reference), you first have to call set_reference(...)");
01138 
01139     rapmax = _reference.rap()+sqrt(_radius_out2);
01140     rapmin = _reference.rap()-sqrt(_radius_out2);
01141   }
01142 
01143   virtual bool is_geometric() const { return true;}    ///< implies a finite area
01144   virtual bool has_finite_area() const { return true;} ///< regardless of the reference 
01145   virtual bool has_known_area() const { return true;}  ///< the area is analytically known
01146   virtual double known_area() const { 
01147     return pi * (_radius_out2-_radius_in2);
01148   }
01149 
01150 protected:
01151   double _radius_in2, _radius_out2;
01152 };
01153 
01154 
01155 
01156 // select on objets with distance from the centre is between 'radius_in' and 'radius_out' 
01157 Selector SelectorDoughnut(const double & radius_in, const double & radius_out) {
01158   return Selector(new SW_Doughnut(radius_in, radius_out));
01159 }
01160 
01161 
01162 //----------------------------------------------------------------------
01163 /// helper for selecting on objects with rapidity within a distance 'delta' of a reference
01164 class SW_Strip : public SW_WithReference {
01165 public:
01166   SW_Strip(const double &delta) : _delta(delta) {}
01167 
01168   /// return a copy of the current object
01169   virtual SelectorWorker* copy(){ return new SW_Strip(*this);}
01170 
01171   /// returns true if a given object passes the selection criterium
01172   /// this has to be overloaded by derived workers
01173   virtual bool pass(const PseudoJet & jet) const {
01174     // make sure the centre is initialised
01175     if (! _is_initialised)
01176       throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
01177     
01178     return abs(jet.rap()-_reference.rap()) <= _delta;
01179   } 
01180 
01181   /// returns a description of the worker
01182   virtual string description() const {
01183     ostringstream ostr;
01184     ostr << "|rap - rap_reference| <= " << _delta;
01185     return ostr.str();
01186   }
01187 
01188   /// returns the rapidity range for which it may return "true"
01189   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
01190     // make sure the centre is initialised
01191     if (! _is_initialised)
01192       throw Error("To use a SelectorStrip (or any selector that requires a reference), you first have to call set_reference(...)");
01193     
01194     rapmax = _reference.rap()+_delta;
01195     rapmin = _reference.rap()-_delta;
01196   }
01197 
01198   virtual bool is_geometric() const { return true;}    ///< implies a finite area
01199   virtual bool has_finite_area() const { return true;} ///< regardless of the reference 
01200   virtual bool has_known_area() const { return true;}  ///< the area is analytically known
01201   virtual double known_area() const { 
01202     return twopi * 2 * _delta;
01203   }
01204 
01205 protected:
01206   double _delta;
01207 };
01208 
01209 
01210 // select on objets within a distance 'radius' of a variable location
01211 Selector SelectorStrip(const double & half_width) {
01212   return Selector(new SW_Strip(half_width));
01213 }
01214 
01215 
01216 //----------------------------------------------------------------------
01217 /// helper for selecting on objects with rapidity within a distance
01218 /// 'delta_rap' of a reference and phi within a distanve delta_phi of
01219 /// a reference
01220 class SW_Rectangle : public SW_WithReference {
01221 public:
01222   SW_Rectangle(const double &delta_rap, const double &delta_phi)
01223     : _delta_rap(delta_rap),  _delta_phi(delta_phi) {}
01224 
01225   /// return a copy of the current object
01226   virtual SelectorWorker* copy(){ return new SW_Rectangle(*this);}
01227 
01228   /// returns true if a given object passes the selection criterium
01229   /// this has to be overloaded by derived workers
01230   virtual bool pass(const PseudoJet & jet) const {
01231     // make sure the centre is initialised
01232     if (! _is_initialised)
01233       throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
01234 
01235     return (abs(jet.rap()-_reference.rap()) <= _delta_rap) && (abs(jet.delta_phi_to(_reference)) <= _delta_phi);
01236   } 
01237 
01238   /// returns a description of the worker
01239   virtual string description() const {
01240     ostringstream ostr;
01241     ostr << "|rap - rap_reference| <= " << _delta_rap << " && |phi - phi_reference| <= " << _delta_phi ;
01242     return ostr.str();
01243   }
01244 
01245   /// returns the rapidity range for which it may return "true"
01246   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
01247     // make sure the centre is initialised
01248     if (! _is_initialised)
01249       throw Error("To use a SelectorRectangle (or any selector that requires a reference), you first have to call set_reference(...)");
01250 
01251     rapmax = _reference.rap()+_delta_rap;
01252     rapmin = _reference.rap()-_delta_rap;
01253   }
01254 
01255   virtual bool is_geometric() const { return true;}    ///< implies a finite area
01256   virtual bool has_finite_area() const { return true;} ///< regardless of the reference 
01257   virtual bool has_known_area() const { return true;}  ///< the area is analytically known
01258   virtual double known_area() const { 
01259     return 4 * _delta_rap * _delta_phi;
01260   }
01261 
01262 protected:
01263   double _delta_rap, _delta_phi;
01264 };
01265 
01266 
01267 // select on objets within a distance 'radius' of a variable location
01268 Selector SelectorRectangle(const double & half_rap_width, const double & half_phi_width) {
01269   return Selector(new SW_Rectangle(half_rap_width, half_phi_width));
01270 }
01271 
01272 
01273 //----------------------------------------------------------------------
01274 /// helper for selecting the jets that carry at least a given fraction
01275 /// of the reference jet
01276 class SW_PtFractionMin : public SW_WithReference {
01277 public:
01278   /// ctor with specification of the number of objects to keep
01279   SW_PtFractionMin(double fraction) : _fraction2(fraction*fraction){}
01280 
01281   /// return a copy of the current object
01282   virtual SelectorWorker* copy(){ return new SW_PtFractionMin(*this);}
01283 
01284   /// return true if the jet carries a large enough fraction of the reference.
01285   /// Throw an error if the reference is not initialised.
01286   virtual bool pass(const PseudoJet & jet) const {
01287     // make sure the centre is initialised
01288     if (! _is_initialised)
01289       throw Error("To use a SelectorPtFractionMin (or any selector that requires a reference), you first have to call set_reference(...)");
01290 
01291     // otherwise, just call that method on the jet
01292     return (jet.perp2() >= _fraction2*_reference.perp2());
01293   }
01294   
01295   /// returns a description of the worker
01296   virtual string description() const {
01297     ostringstream ostr;
01298     ostr << "pt >= " << sqrt(_fraction2) << "* pt_ref";
01299     return ostr.str();
01300   }
01301 
01302 protected:
01303   double _fraction2;
01304 };
01305 
01306 
01307 // select objects that carry at least a fraction "fraction" of the reference jet
01308 // (Note that this selectir takes a reference)
01309 Selector SelectorPtFractionMin(double fraction){
01310   return Selector(new SW_PtFractionMin(fraction));
01311 }
01312 
01313 
01314 //----------------------------------------------------------------------
01315 // additional (mostly helper) selectors
01316 //----------------------------------------------------------------------
01317 
01318 //----------------------------------------------------------------------
01319 /// helper for selecting the 0-momentum jets
01320 class SW_IsZero : public SelectorWorker {
01321 public:
01322   /// ctor
01323   SW_IsZero(){}
01324 
01325   /// return true if the jet has zero momentum
01326   virtual bool pass(const PseudoJet & jet) const {
01327     return jet==0;
01328   }
01329   
01330   /// rereturns a description of the worker
01331   virtual string description() const { return "zero";}
01332 };
01333 
01334 
01335 // select objects with zero momentum
01336 Selector SelectorIsZero(){
01337   return Selector(new SW_IsZero());
01338 }
01339 
01340 
01341 //----------------------------------------------------------------------
01342 /// helper for selecting the pure ghost
01343 class SW_IsPureGhost : public SelectorWorker {
01344 public:
01345   /// ctor
01346   SW_IsPureGhost(){}
01347 
01348   /// return true if the jet is a pure-ghost jet
01349   virtual bool pass(const PseudoJet & jet) const {
01350     // if the jet has no area support then it's certainly not a ghost
01351     if (!jet.has_area()) return false;
01352 
01353     // otherwise, just call that method on the jet
01354     return jet.is_pure_ghost();
01355   }
01356   
01357   /// rereturns a description of the worker
01358   virtual string description() const { return "pure ghost";}
01359 };
01360 
01361 
01362 // select objects that are (or are only made of) ghosts
01363 Selector SelectorIsPureGhost(){
01364   return Selector(new SW_IsPureGhost());
01365 }
01366 
01367 
01368 //----------------------------------------------------------------------
01369 // Selector and workers for obtaining a Selector from an old
01370 // RangeDefinition
01371 //
01372 // This is mostly intended for backward compatibility and is likely to
01373 // be removed in a future major release of FastJet
01374 //----------------------------------------------------------------------
01375 
01376 //----------------------------------------------------------------------
01377 /// helper for selecting on both rapidity and azimuthal angle
01378 class SW_RangeDefinition : public SelectorWorker{
01379 public:
01380   /// ctor from a RangeDefinition
01381   SW_RangeDefinition(const RangeDefinition &range) : _range(&range){}
01382 
01383   /// transfer the selection creterium to the underlying RangeDefinition
01384   virtual bool pass(const PseudoJet & jet) const {
01385     return _range->is_in_range(jet);
01386   } 
01387 
01388   /// returns a description of the worker
01389   virtual string description() const {
01390     return _range->description();
01391   }
01392 
01393   /// returns the rapidity range for which it may return "true"
01394   virtual void get_rapidity_extent(double & rapmin, double & rapmax) const{
01395     _range->get_rap_limits(rapmin, rapmax);
01396   }
01397 
01398   /// check if it has a finite area
01399   virtual bool is_geometric() const { return true;}
01400 
01401   /// check if it has an analytically computable area
01402   virtual bool has_known_area() const { return true;}
01403   
01404   /// if it has a computable area, return it
01405   virtual double known_area() const{
01406     return _range->area();
01407   }
01408 
01409 protected:
01410   const RangeDefinition *_range;
01411 };
01412 
01413 
01414 // ctor from a RangeDefinition
01415 //----------------------------------------------------------------------
01416 //
01417 // This is provided for backward compatibility and will be removed in
01418 // a future major release of FastJet
01419 Selector::Selector(const RangeDefinition &range) {
01420   _worker.reset(new SW_RangeDefinition(range));
01421 }
01422 
01423 
01424 // operators applying directly on a Selector
01425 //----------------------------------------------------------------------
01426 
01427 // operator &=
01428 // For 2 Selectors a and b, a &= b is eauivalent to a = a & b;
01429 Selector & Selector::operator &=(const Selector & b){
01430   _worker.reset(new SW_And(*this, b));
01431   return *this;
01432 }
01433 
01434 // operator &=
01435 // For 2 Selectors a and b, a &= b is eauivalent to a = a & b;
01436 Selector & Selector::operator |=(const Selector & b){
01437   _worker.reset(new SW_Or(*this, b));
01438   return *this;
01439 }
01440 
01441 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends