Phase2L1ParticleFlow/interface/CaloClusterer.h

0001 #ifndef L1Trigger_Phase2L1ParticleFlow_CALOCLUSTERER_H
0002 #define L1Trigger_Phase2L1ParticleFlow_CALOCLUSTERER_H
0003 /**
0004  * Classes for calorimetric re-clustering
0005  * */
0006
0007 // fwd declarations
0008 namespace edm {
0009   class ParameterSet;
0010 }
0011
0012 // real includes
0013 #include <cstdint>
0014 #include <cmath>
0015 #include <vector>
0016 #include <array>
0017 #include <algorithm>
0018 #include "DataFormats/L1TParticleFlow/interface/PFCluster.h"
0019 #include "DataFormats/Common/interface/OrphanHandle.h"
0020
0021 namespace l1tpf_calo {
0022   class Grid {
0023   public:
0024     virtual ~Grid() {}
0025     unsigned int size() const { return ncells_; }
0026     virtual int find_cell(float eta, float phi) const = 0;
0027     int neighbour(int icell, unsigned int idx) const { return neighbours_[icell][idx]; }
0028     float eta(int icell) const { return eta_[icell]; }
0029     float phi(int icell) const { return phi_[icell]; }
0030     float etaWidth(int icell) const { return etaWidth_[icell]; }
0031     float phiWidth(int icell) const { return phiWidth_[icell]; }
0032     int ieta(int icell) const { return ieta_[icell]; }
0033     int iphi(int icell) const { return iphi_[icell]; }
0034
0035   protected:
0036     Grid(unsigned int size)
0037         : ncells_(size),
0038           eta_(size),
0039           etaWidth_(size),
0040           phi_(size),
0041           phiWidth_(size),
0042           ieta_(size),
0043           iphi_(size),
0044           neighbours_(size) {}
0045     unsigned int ncells_;
0046     std::vector<float> eta_, etaWidth_, phi_, phiWidth_;
0047     std::vector<int> ieta_, iphi_;
0048     std::vector<std::array<int, 8>> neighbours_;  // indices of the neigbours, -1 = none
0049   };
0050
0051   class Phase1GridBase : public Grid {
0052   public:
0053     Phase1GridBase(int nEta, int nPhi, int ietaCoarse, int ietaVeryCoarse, const float *towerEtas);
0054
0055     int find_cell(float eta, float phi) const override;
0056     int ifind_cell(int ieta, int iphi) const { return cell_map_[(ieta + nEta_) + 2 * nEta_ * (iphi - 1)]; }
0057
0058   protected:
0059     const int nEta_, nPhi_, ietaCoarse_, ietaVeryCoarse_;
0060     const float *towerEtas_;
0061     std::vector<int> cell_map_;
0062     // valid ieta, iphi (does not check for outside bounds, only for non-existence of ieta=0, iphi=0, and coarser towers at high eta)
0063     bool valid_ieta_iphi(int ieta, int iphi) const {
0064       if (ieta == 0 || iphi == 0)
0065         return false;
0066       if (std::abs(ieta) >= ietaVeryCoarse_ && (iphi % 4 != 1))
0067         return false;
0068       if (std::abs(ieta) >= ietaCoarse_ && (iphi % 2 != 1))
0069         return false;
0070       return true;
0071     }
0072     // move by +/-1 around a cell; return icell or -1 if not available
0073     int imove(int ieta, int iphi, int deta, int dphi);
0074   };
0075
0076   class Phase1Grid : public Phase1GridBase {
0077   public:
0078     Phase1Grid()
0079         : Phase1GridBase(phase1_nEta_, phase1_nPhi_, phase1_ietaCoarse_, phase1_ietaVeryCoarse_, phase1_towerEtas_) {}
0080
0081   protected:
0082     static const int phase1_nEta_ = 41, phase1_nPhi_ = 72, phase1_ietaCoarse_ = 29, phase1_ietaVeryCoarse_ = 40;
0083     static const float phase1_towerEtas_[phase1_nEta_];
0084   };
0085   class Phase2Grid : public Phase1GridBase {
0086   public:
0087     Phase2Grid()
0088         : Phase1GridBase(phase2_nEta_, phase2_nPhi_, phase2_ietaCoarse_, phase2_ietaVeryCoarse_, phase2_towerEtas_) {}
0089
0090   protected:
0091     static const int phase2_nEta_ = 48, phase2_nPhi_ = 72, phase2_ietaCoarse_ = 36, phase2_ietaVeryCoarse_ = 47;
0092     static const float phase2_towerEtas_[phase2_nEta_];
0093   };
0094
0095   template <typename T>
0096   class GridData {
0097   public:
0098     GridData() : grid_(nullptr), data_(), empty_() {}
0099     GridData(const Grid &grid) : grid_(&grid), data_(grid.size()), empty_() {}
0100
0101     T &operator()(float eta, float phi) { return data_[grid_->find_cell(eta, phi)]; }
0102     const T &operator()(float eta, float phi) const { return data_[grid_->find_cell(eta, phi)]; }
0103
0104     const Grid &grid() const { return *grid_; }
0105
0106     unsigned int size() const { return data_.size(); }
0107
0108     float eta(int icell) const { return grid().eta(icell); }
0109     float phi(int icell) const { return grid().phi(icell); }
0110     int ieta(int icell) const { return grid().ieta(icell); }
0111     int iphi(int icell) const { return grid().iphi(icell); }
0112
0113     T &operator[](int icell) { return data_[icell]; }
0114     const T &operator[](int icell) const { return data_[icell]; }
0115
0116     const T &neigh(int icell, unsigned int idx) const {
0117       int ineigh = grid_->neighbour(icell, idx);
0118       return (ineigh < 0 ? empty_ : data_[ineigh]);
0119     }
0120
0121     GridData<T> &operator=(const GridData<T> &other) {
0122       assert(grid_ == other.grid_);
0123       data_ = other.data_;
0124       return *this;
0125     }
0126     GridData<T> &operator+=(const GridData<T> &other) {
0127       assert(grid_ == other.grid_);
0128       for (unsigned int i = 0, n = data_.size(); i < n; ++i) {
0129         data_[i] += other.data_[i];
0130       }
0131       return *this;
0132     }
0133
0134     // always defined
0135     void fill(const T &val) { std::fill(data_.begin(), data_.end(), val); }
0136     void zero() { fill(T()); }
0137
0138     // defined only if T has a 'clear' method
0139     void clear() {
0140       for (T &t : data_)
0141         t.clear();
0142     }
0143
0144   private:
0145     const Grid *grid_;
0146     std::vector<T> data_;
0147     const T empty_;
0148   };
0149   typedef GridData<float> EtGrid;
0150   typedef GridData<int> IndexGrid;
0151
0152   struct PreCluster {
0153     PreCluster() : ptLocalMax(0), ptOverNeighLocalMaxSum(0) {}
0154     float ptLocalMax;              // pt if it's a local max, zero otherwise
0155     float ptOverNeighLocalMaxSum;  // pt / (sum of ptLocalMax of neighbours); zero if no neighbours
0156     void clear() { ptLocalMax = ptOverNeighLocalMaxSum = 0; }
0157   };
0158   typedef GridData<PreCluster> PreClusterGrid;
0159
0160   struct Cluster {
0161     Cluster() : et(0), eta(0), phi(0) {}
0162     float et, eta, phi;
0163     std::vector<std::pair<int, float>> constituents;
0164     void clear() {
0165       et = eta = phi = 0;
0166       constituents.clear();
0167     }
0168   };
0169
0170   struct CombinedCluster : public Cluster {
0171     float ecal_et, hcal_et;
0172     void clear() {
0173       Cluster::clear();
0174       ecal_et = hcal_et = 0;
0175     }
0176   };
0177
0178   const Grid *getGrid(const std::string &type);
0179
0180   class SingleCaloClusterer {
0181   public:
0182     SingleCaloClusterer(const edm::ParameterSet &pset);
0183     ~SingleCaloClusterer();
0184     void clear();
0185     void add(const reco::Candidate &c) { add(c.pt(), c.eta(), c.phi()); }
0186     void add(float pt, float eta, float phi) { rawet_(eta, phi) += pt; }
0187     void run();
0188
0189     /// possibly grow clusters by adding unclustered energy on the sides
0190     //  note: there can be some double-counting as the same unclustered energy can go into more clusters
0191     void grow();
0192
0193     const EtGrid &raw() const { return rawet_; }
0194     const IndexGrid &indexGrid() const { return clusterIndex_; }
0195     const std::vector<Cluster> &clusters() const { return clusters_; }
0196     const Cluster &cluster(int i) const {
0197       return (i == -1 || clusterIndex_[i] == -1) ? nullCluster_ : clusters_[clusterIndex_[i]];
0198     }
0199
0200     /// non-const access to the energy: be careful to use it only before 'run()'
0201     EtGrid &raw() { return rawet_; }
0202
0203     // for the moment, generic interface that takes a cluster and returns the corrected pt
0204     template <typename Corrector>
0205     void correct(const Corrector &corrector) {
0206       for (Cluster &c : clusters_) {
0207         c.et = corrector(c);
0208       }
0209     }
0210
0211     std::unique_ptr<l1t::PFClusterCollection> fetchCells(bool unclusteredOnly = false, float ptMin = 0.) const;
0212
0213     std::unique_ptr<l1t::PFClusterCollection> fetch(float ptMin = 0.) const;
0214     std::unique_ptr<l1t::PFClusterCollection> fetch(const edm::OrphanHandle<l1t::PFClusterCollection> &cells,
0215                                                     float ptMin = 0.) const;
0216
0217   private:
0218     enum class EnergyShareAlgo {
0219       Fractions, /* each local maximum neighbour takes a share proportional to its value */
0220       None,      /* each local maximum neighbour takes all the value (double counting!) */
0221       Greedy,    /* assing cell to the highest local maximum neighbour */
0222       Crude
0223     }; /* if there's more than one local maximum neighbour, they all take half of the value (no fp division) */
0224     const Grid *grid_;
0225     EtGrid rawet_, unclustered_;
0226     PreClusterGrid precluster_;
0227     IndexGrid clusterIndex_, cellKey_;
0228     std::vector<Cluster> clusters_;
0229     const Cluster nullCluster_;
0230     float zsEt_, seedEt_, minClusterEt_, minEtToGrow_;
0231     EnergyShareAlgo energyShareAlgo_;
0232     bool energyWeightedPosition_;  // do the energy-weighted cluster position instead of the cell center
0233   };
0234
0235   class SimpleCaloLinkerBase {
0236   public:
0237     SimpleCaloLinkerBase(const edm::ParameterSet &pset,
0238                          const SingleCaloClusterer &ecal,
0239                          const SingleCaloClusterer &hcal);
0240     virtual ~SimpleCaloLinkerBase();
0241     virtual void clear() { clearBase(); }
0242     virtual void run() = 0;
0243     void clearBase() {
0244       clusters_.clear();
0245       clusterIndex_.fill(-1);
0246     }
0247
0248     // for the moment, generic interface that takes a cluster and returns the corrected pt
0249     template <typename Corrector>
0250     void correct(const Corrector &corrector) {
0251       for (CombinedCluster &c : clusters_) {
0252         c.et = corrector(c);
0253       }
0254     }
0255
0256     std::unique_ptr<l1t::PFClusterCollection> fetch() const;
0257     std::unique_ptr<l1t::PFClusterCollection> fetch(const edm::OrphanHandle<l1t::PFClusterCollection> &ecal,
0258                                                     const edm::OrphanHandle<l1t::PFClusterCollection> &hcal) const;
0259
0260   protected:
0261     const Grid *grid_;
0262     const SingleCaloClusterer &ecal_, &hcal_;
0263     IndexGrid clusterIndex_;
0264     std::vector<CombinedCluster> clusters_;
0265     float hoeCut_, minPhotonEt_, minHadronRawEt_, minHadronEt_;
0266     bool noEmInHGC_;
0267   };
0268
0269   class SimpleCaloLinker : public SimpleCaloLinkerBase {
0270   public:
0271     SimpleCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer &ecal, const SingleCaloClusterer &hcal);
0272     ~SimpleCaloLinker() override;
0273     void clear() override;
0274     void run() override;
0275
0276   protected:
0277     PreClusterGrid ecalToHCal_;
0278   };
0279   class FlatCaloLinker : public SimpleCaloLinkerBase {
0280   public:
0281     FlatCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer &ecal, const SingleCaloClusterer &hcal);
0282     ~FlatCaloLinker() override;
0283     void clear() override;
0284     void run() override;
0285
0286   protected:
0287     SingleCaloClusterer combClusterer_;
0288   };
0289
0290   // makes a calo linker (pointer will be owned by the callee)
0291   std::unique_ptr<SimpleCaloLinkerBase> makeCaloLinker(const edm::ParameterSet &pset,
0292                                                        const SingleCaloClusterer &ecal,
0293                                                        const SingleCaloClusterer &hcal);
0294
0295 }  // namespace l1tpf_calo
0296
0297 #endif