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/ATLASConePlugin.hh"
00019 #include "fastjet/CMSIterativeConePlugin.hh"
00020 #include "fastjet/CDFJetCluPlugin.hh"
00021 #include "fastjet/CDFMidPointPlugin.hh"
00022 #include "fastjet/D0RunIIConePlugin.hh"
00023 #include "fastjet/TrackJetPlugin.hh"
00024 #include "fastjet/JadePlugin.hh"
00025 //#include "fastjet/PxConePlugin.hh"
00026 
00027 namespace Rivet {
00028 
00029 
00031   inline Vector3 momentum3(const fastjet::PseudoJet& pj) {
00032     return Vector3(pj.px(), pj.py(), pj.pz());
00033   }
00034 
00036   inline FourMomentum momentum(const fastjet::PseudoJet& pj) {
00037     return FourMomentum(pj.E(), pj.px(), pj.py(), pj.pz());
00038   }
00039 
00040 
00042   typedef vector<fastjet::PseudoJet> PseudoJets;
00043 
00044 
00045 
00047 
00048 
00049 
00051   class FastJets : public JetAlg {
00052 
00053   public:
00055     enum JetAlgName { KT, CAM, SISCONE, ANTIKT,
00056                       PXCONE,
00057                       ATLASCONE, CMSCONE,
00058                       CDFJETCLU, CDFMIDPOINT, D0ILCONE,
00059                       JADE, DURHAM, TRACKJET };
00060 
00061 
00063 
00064 
00069     FastJets(const FinalState& fsp, JetAlgName alg,
00070              double rparameter, double seed_threshold=1.0)
00071       : JetAlg(fsp), _adef(0) { _init1(alg, rparameter, seed_threshold); }
00072 
00074     FastJets(const FinalState& fsp, fastjet::JetAlgorithm type,
00075              fastjet::RecombinationScheme recom, double rparameter)
00076       : JetAlg(fsp), _adef(0) { _init2(type, recom, rparameter); }
00077 
00079     FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin* plugin)
00080       : JetAlg(fsp), _adef(0) { _init3(plugin); }
00082     FastJets(const FinalState& fsp, fastjet::JetDefinition::Plugin& plugin)
00083       : JetAlg(fsp), _adef(0) { _init3(&plugin); }
00084 
00085 
00087     FastJets(JetAlgName alg, double rparameter, double seed_threshold=1.0)
00088       : _adef(0) { _init1(alg, rparameter, seed_threshold); }
00090     FastJets(fastjet::JetAlgorithm type, fastjet::RecombinationScheme recom, double rparameter)
00091       : _adef(0) { _init2(type, recom, rparameter); }
00093     FastJets(fastjet::JetDefinition::Plugin* plugin)
00094       : _adef(0) { _init3(plugin); }
00096     FastJets(fastjet::JetDefinition::Plugin& plugin)
00097       : _adef(0) { _init3(&plugin); }
00098 
00099 
00101     virtual const Projection* clone() const {
00102       return new FastJets(*this);
00103     }
00104 
00106 
00107 
00108   public:
00109 
00111     void reset();
00112 
00116     void useJetArea(fastjet::AreaDefinition* adef) {
00117       _adef = adef;
00118     }
00119 
00121     size_t numJets(double ptmin = 0.0) const;
00122 
00124     size_t size() const {
00125       return numJets();
00126     }
00127 
00129     Jets _jets(double ptmin = 0.0) const;
00130 
00132     PseudoJets pseudoJets(double ptmin = 0.0) const;
00133 
00135     PseudoJets pseudoJetsByPt(double ptmin = 0.0) const {
00136       return sorted_by_pt(pseudoJets(ptmin));
00137     }
00138 
00140     PseudoJets pseudoJetsByE(double ptmin = 0.0) const {
00141       return sorted_by_E(pseudoJets(ptmin));
00142     }
00143 
00145     PseudoJets pseudoJetsByRapidity(double ptmin = 0.0) const {
00146       return sorted_by_rapidity(pseudoJets(ptmin));
00147     }
00148 
00150     const fastjet::ClusterSequence* clusterSeq() const {
00151       return _cseq.get();
00152     }
00153 
00155     const fastjet::ClusterSequenceArea* clusterSeqArea() const {
00157       if (_adef == 0) return (fastjet::ClusterSequenceArea*) 0;
00158       return dynamic_cast<fastjet::ClusterSequenceArea*>(_cseq.get());
00159     }
00160 
00162     const fastjet::JetDefinition& jetDef() const {
00163       return _jdef;
00164     }
00165 
00167     const fastjet::AreaDefinition* areaDef() const {
00168       return _adef;
00169     }
00170 
00172     vector<double> ySubJet(const fastjet::PseudoJet& jet) const;
00173 
00176     fastjet::PseudoJet splitJet(fastjet::PseudoJet jet, double& last_R) const;
00177 
00180     fastjet::PseudoJet filterJet(fastjet::PseudoJet jet, double& stingy_R, const double def_R) const;
00181 
00182   private:
00183 
00184     Jets _pseudojetsToJets(const PseudoJets& pjets) const;
00185 
00187     void _init1(JetAlgName alg, double rparameter, double seed_threshold);
00188     void _init2(fastjet::JetAlgorithm type,
00189                 fastjet::RecombinationScheme recom, double rparameter);
00190     void _init3(fastjet::JetDefinition::Plugin* plugin);
00191 
00192   protected:
00193 
00195     void project(const Event& e);
00196 
00198     int compare(const Projection& p) const;
00199 
00200   public:
00201 
00203     void calc(const ParticleVector& ps);
00204 
00205 
00206   private:
00207 
00209     fastjet::JetDefinition _jdef;
00210 
00212     fastjet::AreaDefinition* _adef;
00213 
00215     shared_ptr<fastjet::ClusterSequence> _cseq;
00216 
00218     shared_ptr<fastjet::JetDefinition::Plugin> _plugin;
00219 
00221     mutable map<int, vector<double> > _yscales;
00222 
00224     //set<Particle, ParticleBase::byPTAscending> _particles;
00225     map<int, Particle> _particles;
00226 
00227   };
00228 
00229 }
00230 
00231 #endif