Rivet  1.8.0
FastJets.hh
00001 // -*- C++ -*-
00002 #ifndef RIVET_FastJets_HH
00003 #define RIVET_FastJets_HH
00004 
00005 #include "Rivet/Projection.hh"
00006 #include "Rivet/Projections/JetAlg.hh"
00007 #include "Rivet/Projections/FinalState.hh"
00008 #include "Rivet/Particle.hh"
00009 #include "Rivet/Jet.hh"
00010 
00011 #include "fastjet/JetDefinition.hh"
00012 #include "fastjet/AreaDefinition.hh"
00013 #include "fastjet/ClusterSequence.hh"
00014 #include "fastjet/ClusterSequenceArea.hh"
00015 #include "fastjet/PseudoJet.hh"
00016 
00017 #include "fastjet/SISConePlugin.hh"
00018 #include "fastjet/CMSIterativeConePlugin.hh"
00019 #include "fastjet/D0RunIIConePlugin.hh"
00020 #include "fastjet/TrackJetPlugin.hh"
00021 #include "fastjet/JadePlugin.hh"
00022 //#include "fastjet/PxConePlugin.hh"
00023 
00024 namespace Rivet {
00025 
00026 
00028   inline Vector3 momentum3(const fastjet::PseudoJet& pj) {
00029     return Vector3(pj.px(), pj.py(), pj.pz());
00030   }
00031 
00033   inline FourMomentum momentum(const fastjet::PseudoJet& pj) {
00034     return FourMomentum(pj.E(), pj.px(), pj.py(), pj.pz());
00035   }
00036 
00037 
00039   typedef vector<fastjet::PseudoJet> PseudoJets;
00040 
00041 
00042 
00044 
00045 
00046 
00048   class FastJets : public JetAlg {
00049 
00050   public:
00052     enum JetAlgName { KT, CAM, SISCONE, ANTIKT,
00053                       PXCONE,
00054                       ATLASCONE, CMSCONE,
00055                       CDFJETCLU, CDFMIDPOINT, D0ILCONE,
00056                       JADE, DURHAM, TRACKJET };
00057 
00058 
00060 
00061 
00066     FastJets(const FinalState& fsp, JetAlgName alg,
00067              double rparameter, double seed_threshold=1.0)
00068       : JetAlg(fsp), _adef(0) { _init1(alg, rparameter, seed_threshold); }
00069 
00071     FastJets(const FinalState& fsp, fastjet::JetAlgorithm type,
00072              fastjet::RecombinationScheme recom, double rparameter)
00073       : JetAlg(fsp), _adef(0) { _init2(type, recom, rparameter); }
00074 
00076     FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin* plugin)
00077       : JetAlg(fsp), _adef(0) { _init3(plugin); }
00079     FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin& plugin)
00080       : JetAlg(fsp), _adef(0) { _init3(&plugin); }
00081 
00082 
00084     FastJets(JetAlgName alg, double rparameter, double seed_threshold=1.0)
00085       : _adef(0) { _init1(alg, rparameter, seed_threshold); }
00087     FastJets(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter)
00088       : _adef(0) { _init2(type, recom, rparameter); }
00090     FastJets(fastjet::JetDefinition::Plugin* plugin)
00091       : _adef(0) { _init3(plugin); }
00093     FastJets(fastjet::JetDefinition::Plugin& plugin)
00094       : _adef(0) { _init3(&plugin); }
00095 
00096 
00098     virtual const Projection* clone() const {
00099       return new FastJets(*this);
00100     }
00101 
00103 
00104 
00105   public:
00106 
00108     void reset();
00109 
00113     void useJetArea(fastjet::AreaDefinition* adef) {
00114       _adef = adef;
00115     }
00116 
00118     size_t numJets(double ptmin = 0.0) const;
00119 
00121     size_t size() const {
00122       return numJets();
00123     }
00124 
00126     Jets _jets(double ptmin = 0.0) const;
00127 
00129     PseudoJets pseudoJets(double ptmin = 0.0) const;
00130 
00132     PseudoJets pseudoJetsByPt(double ptmin = 0.0) const {
00133       return sorted_by_pt(pseudoJets(ptmin));
00134     }
00135 
00137     PseudoJets pseudoJetsByE(double ptmin = 0.0) const {
00138       return sorted_by_E(pseudoJets(ptmin));
00139     }
00140 
00142     PseudoJets pseudoJetsByRapidity(double ptmin = 0.0) const {
00143       return sorted_by_rapidity(pseudoJets(ptmin));
00144     }
00145 
00147     const fastjet::ClusterSequence* clusterSeq() const {
00148       return _cseq.get();
00149     }
00150 
00152     const fastjet::ClusterSequenceArea* clusterSeqArea() const {
00154       if (_adef == 0) return (fastjet::ClusterSequenceArea*) 0;
00155       return dynamic_cast<fastjet::ClusterSequenceArea*>(_cseq.get());
00156     }
00157 
00159     const fastjet::JetDefinition& jetDef() const {
00160       return _jdef;
00161     }
00162 
00164     const fastjet::AreaDefinition* areaDef() const {
00165       return _adef;
00166     }
00167 
00169     vector<double> ySubJet(const fastjet::PseudoJet& jet) const;
00170 
00173     fastjet::PseudoJet splitJet(fastjet::PseudoJet jet, double& last_R) const;
00174 
00177     fastjet::PseudoJet filterJet(fastjet::PseudoJet jet, double& stingy_R, const double def_R) const;
00178 
00179   private:
00180 
00181     Jets _pseudojetsToJets(const PseudoJets& pjets) const;
00182 
00184     void _init1(JetAlgName alg, double rparameter, double seed_threshold);
00185     void _init2(fastjet::JetAlgorithm type,
00186                 fastjet::RecombinationScheme recom, double rparameter);
00187     void _init3(fastjet::JetDefinition::Plugin* plugin);
00188 
00189   protected:
00190 
00192     void project(const Event& e);
00193 
00195     int compare(const Projection& p) const;
00196 
00197   public:
00198 
00200     void calc(const ParticleVector& ps);
00201 
00202 
00203   private:
00204 
00206     fastjet::JetDefinition _jdef;
00207 
00209     fastjet::AreaDefinition* _adef;
00210 
00212     shared_ptr<fastjet::ClusterSequence> _cseq;
00213 
00215     shared_ptr<fastjet::JetDefinition::Plugin> _plugin;
00216 
00218     mutable map<int, vector<double> > _yscales;
00219 
00221     //set<Particle, ParticleBase::byPTAscending> _particles;
00222     map<int, Particle> _particles;
00223 
00224   };
00225 
00226 }
00227 
00228 #endif