fastjet/contrib/EventStorage.hh

0001 // JetsWithoutJets Package
0002 // Questions/Comments? danbert@mit.edu jthaler@mit.edu
0003 //
0004 // Copyright (c) 2013
0005 // Daniele Bertolini and Jesse Thaler
0006 //
0007 // $Id: EventStorage.hh 554 2014-02-21 19:02:08Z danbert $
0008 //----------------------------------------------------------------------
0009 // This file is part of FastJet contrib.
0010 //
0011 // It is free software; you can redistribute it and/or modify it under
0012 // the terms of the GNU General Public License as published by the
0013 // Free Software Foundation; either version 2 of the License, or (at
0014 // your option) any later version.
0015 //
0016 // It is distributed in the hope that it will be useful, but WITHOUT
0017 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
0018 // or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
0019 // License for more details.
0020 //
0021 // You should have received a copy of the GNU General Public License
0022 // along with this code. If not, see <http://www.gnu.org/licenses/>.
0023 //----------------------------------------------------------------------
0024
0025 #ifndef __FASTJET_CONTRIB_EVENTSTORAGE_HH__
0026 #define __FASTJET_CONTRIB_EVENTSTORAGE_HH__
0027
0028 #include "fastjet/PseudoJet.hh"
0029 #include <sstream>
0030
0031 FASTJET_BEGIN_NAMESPACE
0032
0033 namespace jwj {
0034
0035
0036 //////
0037 //
0038 // ParticleStorage
0039 //
0040 //////
0041
0042
0043 //----------------------------------------------------------------------
0044 // ParticleStorage
0045 // ParticleStorage stores information about single particles.
0046 // By default it stores the full pseudojet, rapidity, phi, pt, mass, px and py.
0047 // This information is used multiple times, so it is convenient to cache it.
0048 // By design, we are storing redundant information (i.e. the full pseudo jet as well as rap/phi, px/py, etc.)
0049 // because we don't want to be ever rely on the internal way the PseudoJet stores/computes this information
0050 // In addition, information about the neighborhood of the particle (i.e. within Rjet/Rsub circles),
0051 // such as pT_in_Rjet, pT_in_Rsub, mass_in_Rjet, indices identifying neighbors, can be stored.
0052 // Storing of this additional information is done externally by the EventStorage class.
0053 // Note that the EventStorage class uses data already cached in ParticleStorage to calculate information about neighborhoods.
0054
0055
0056 class ParticleStorage {
0057
0058 private:
0059
0060    PseudoJet _pj;
0061    double _rap,_phi,_pt,_m,_px,_py,_pt_in_Rjet,_pt_in_Rsub,_m_in_Rjet,_weight;
0062    bool _includeParticle;
0063    std::vector<unsigned int> _neighbors;
0064
0065 public:
0066
0067    ParticleStorage(){}
0068
0069    // by default it stores the full pseudojet, rapidity, phi, pt, mass, px and py (need px and py to calculate vector pt) of the particle
0070    // This information is typically set by EventStorage::_establishBasicStorage().
0071    // Also, _pt_in_Rjet,_pt_in_Rsub,_m_in_Rjet,_weight,_includeParticle are initialized to default values. They can be reset through
0072    // corresponding set functions (typically done by EventStorage::_establishDerivedStorage()).
0073    ParticleStorage(const PseudoJet & myParticle): _pj(myParticle),_rap(myParticle.rap()),_phi(myParticle.phi()),_pt(myParticle.pt()),_m(myParticle.m()),_px(myParticle.px()), _py(myParticle.py()),_pt_in_Rjet(0.0),_pt_in_Rsub(0.0),_m_in_Rjet(0.0),_weight(0.0),_includeParticle(false) {}
0074    ~ParticleStorage(){}
0075
0076    // return stored information
0077    const PseudoJet & pseudoJet() const {return _pj;}
0078    const double & rap() const {return _rap;}
0079    const double & phi() const {return _phi;}
0080    const double & pt() const {return _pt;}
0081    const double & m() const {return _m;}
0082    const double & px() const {return _px;}
0083    const double & py() const {return _py;}
0084
0085    // deltaRsq returns the squared distance in rapidity-azimuth plane between this particle and other
0086    double deltaRsq(const ParticleStorage & other) const;
0087
0088    // The following functions set and call information about the neighborhood of the particle.
0089    // This information will typically be set by EventStorage::_establishDerivedStorage().
0090
0091    // scalar sum pt in a cone of radius Rjet around this particle
0092    void set_pt_in_Rjet(const double & pt_in_Rjet) {_pt_in_Rjet=pt_in_Rjet;}
0093    const double & pt_in_Rjet() const {return _pt_in_Rjet;}
0094
0095    // scalar sum pt in a cone of radius Rsub around this particle
0096    void set_pt_in_Rsub(const double & pt_in_Rsub) {_pt_in_Rsub=pt_in_Rsub;}
0097    const double & pt_in_Rsub() const {return _pt_in_Rsub;}
0098
0099    // mass in a cone of radius Rjet around this particle
0100    void set_m_in_Rjet(const double & m_in_Rjet) {_m_in_Rjet=m_in_Rjet;}
0101    const double & m_in_Rjet() const {return _m_in_Rjet;}
0102
0103    // weight=pt/pt_in_Rjet for this particle
0104    void set_weight(const double & weight) {_weight=weight;}
0105    const double & weight() const {return _weight;}
0106
0107    // true if pt_in_Rjet >= ptcut
0108    void set_includeParticle(const bool & includeParticle) {_includeParticle=includeParticle;}
0109    const bool & includeParticle() const {return _includeParticle;}
0110
0111    // vector of indices of particles within a distance Rjet from this particle
0112    void set_neighbors(const std::vector<unsigned int> & neighbors) {_neighbors=neighbors;}
0113    const std::vector<unsigned int> & neighbors() const {return _neighbors;}
0114
0115 };
0116
0117
0118 //////
0119 //
0120 // Local Storage (to reduce computation time)
0121 //
0122 //////
0123
0124 //----------------------------------------------------------------------
0125 // LocalStorage
0126 // This class divides an event into (overlapping) regions of size 2R x 2R
0127 // and caches the partition for later use.  Partitions that are below
0128 // the ptcut can be ignored entirely.  This saves some computational time by
0129 // avoiding having to calculate irrelevant distances.  For speed purposes, this
0130 // class uses ParticleStorage information (instead of PseudoJets)
0131
0132 class LocalStorage {
0133
0134 private:
0135
0136    double _Rjet;
0137    double _ptcut;
0138
0139 public:
0140
0141    LocalStorage() : _Rjet(0.0) {
0142       // at the moment, this value is hard coded
0143       _rapmax = 10.0;
0144    }
0145
0146    // make the local storage
0147    void establishStorage(const std::vector<ParticleStorage> & myParticles, double Rjet, double ptcut);
0148
0149    // give the 2R x 2R region relevant for a given particle
0150    const std::vector<unsigned int> & getStorageFor(const ParticleStorage & myParticle) const;
0151
0152    // give whether the 2R x 2R region has enough pT to merit further examination
0153    bool aboveCutFor(const ParticleStorage & myParticle);
0154
0155 private:
0156
0157    // storage for partitions of particles
0158    std::vector<std::vector<std::vector<unsigned int> > > _regionStorage;
0159
0160    // bool for whether region needs to be considered
0161    std::vector<std::vector<bool> > _aboveCutBool;
0162
0163    double _rapmax;
0164    int _maxRapIndex;
0165    double _rapSpread;
0166    int _maxPhiIndex;
0167    double _phiSpread;
0168
0169    // helper functions to determine which region a pseudo jet is in
0170    int getRapIndex(const ParticleStorage & myParticle) const;
0171    int getPhiIndex(const ParticleStorage & myParticle) const;
0172
0173    // helper function to determine sum scalar pt of a subset of particles identifed by myIds
0174    double getSumPt(const std::vector<ParticleStorage> & Particles, const std::vector<unsigned int> myIds) const;
0175
0176 };
0177
0178
0179 //////
0180 //
0181 // EventStorage
0182 //
0183 //////
0184
0185 //----------------------------------------------------------------------
0186 // EventStorage
0187 // This class is a wrapper for a vectors of ParticleStorages.
0188 // It uses default information stored in ParticleStorage to calculate additional infomation
0189 // such pt_in_Rjet,pt_in_Rsub etc and it stores this information back to ParticleStorage.
0190 // Can choose if want to use LocalStorage to reduce computation time (on by default)
0191
0192 class EventStorage {
0193
0194 private:
0195
0196    double _Rjet,_ptcut,_Rsub,_fcut;
0197    std::vector<unsigned int> _ids;
0198    std::vector<ParticleStorage> _storage;
0199    bool _useLocalStorage;
0200    bool _storeNeighbors, _storeMass;
0201
0202 public:
0203
0204    EventStorage(){}
0205    // standard constructor, by default it uses LocalStorage and for each particle it stores particles within a distance Rjet (neighbors) and it does
0206    // not store mass of the jet made of neighbors
0207    EventStorage(double Rjet, double ptcut, double Rsub, double fcut, bool useLocalStorage=true, bool storeNeighbors=true, bool storeMass=false):
0208    _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rsub), _fcut(fcut), _useLocalStorage(useLocalStorage), _storeNeighbors(storeNeighbors), _storeMass(storeMass) {}
0209
0210    // constructor with Rsub and fcut automatically initialized to Rjet and 1.0 respectively.
0211    // By default it uses LocalStorage and for each particle it stores particles within a distance Rjet (neighbors) and it does
0212    // not store mass of the jet made of neighbors.
0213    // Use this constructor if don't need to do trimming.
0214    EventStorage(double Rjet, double ptcut, bool useLocalStorage=true, bool storeNeighbors=true, bool storeMass=false):
0215    _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rjet), _fcut(1.0), _useLocalStorage(useLocalStorage), _storeNeighbors(storeNeighbors), _storeMass(storeMass) {}
0216    ~EventStorage() {}
0217
0218    // establish the storage
0219    void establishStorage(const std::vector<PseudoJet> & particles){
0220       _establishBasicStorage(particles);
0221       _establishDerivedStorage();
0222    }
0223
0224    // get general info about the storage
0225    unsigned int size() const {return _storage.size();}
0226    const double & Rjet() const {return _Rjet;}
0227    const double & ptcut() const {return _ptcut;}
0228    const double & Rsub() const {return _Rsub;}
0229    const double & fcut() const {return _fcut;}
0230    const bool & storeNeighbors() const {return _storeNeighbors;}
0231    const bool & storeMass() const {return _storeMass;}
0232
0233    // access individual ParticleStorage through [] operator
0234    ParticleStorage operator[](const unsigned int i) const {return _storage[i];}
0235
0236    // convert stored indices of neighbors (particles within a distance Rjet) into a vector<PseudoJet>
0237    std::vector<PseudoJet> particles_near_to(const unsigned int id) const {
0238       std::vector<unsigned int> neighbors=_storage[id].neighbors();
0239       std::vector<PseudoJet> answer;
0240       for(unsigned int i=0; i<neighbors.size(); i++) answer.push_back(_storage[neighbors[i]].pseudoJet());
0241       return answer;
0242    }
0243
0244    // Description
0245    std::string parameterString() const {
0246       std::stringstream stream;
0247       stream << "R_jet=" << _Rjet << ", pT_cut=" << _ptcut << ", R_sub=" << _Rsub <<", fcut=" << _fcut;
0248       return stream.str();
0249    }
0250    std::string description() const {
0251       return "Event Storage with "+parameterString();
0252    }
0253
0254 private:
0255
0256    void _establishBasicStorage(const std::vector<PseudoJet> & particles);
0257    void _establishDerivedStorage();
0258    void _get_local_info(const unsigned int id, const std::vector<unsigned int>* myLocalRegion, double & pt_in_Rjet, double & pt_in_Rsub, double & m_in_Rjet, std::vector<unsigned int> & neighbors) const;
0259
0260 };
0261
0262 } // namespace jwj
0263
0264 FASTJET_END_NAMESPACE
0265
0266 #endif