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     float eta(float eta, float phi) const { return grid().eta(grid_->find_cell(eta, phi)); }
0105     float phi(float eta, float phi) const { return grid().phi(grid_->find_cell(eta, phi)); }
0106
0107     const Grid &grid() const { return *grid_; }
0108
0109     unsigned int size() const { return data_.size(); }
0110
0111     float eta(int icell) const { return grid().eta(icell); }
0112     float phi(int icell) const { return grid().phi(icell); }
0113     int ieta(int icell) const { return grid().ieta(icell); }
0114     int iphi(int icell) const { return grid().iphi(icell); }
0115
0116     T &operator[](int icell) { return data_[icell]; }
0117     const T &operator[](int icell) const { return data_[icell]; }
0118
0119     const T &neigh(int icell, unsigned int idx) const {
0120       int ineigh = grid_->neighbour(icell, idx);
0121       return (ineigh < 0 ? empty_ : data_[ineigh]);
0122     }
0123
0124     GridData<T> &operator=(const GridData<T> &other) {
0125       assert(grid_ == other.grid_);
0126       data_ = other.data_;
0127       return *this;
0128     }
0129     GridData<T> &operator+=(const GridData<T> &other) {
0130       assert(grid_ == other.grid_);
0131       for (unsigned int i = 0, n = data_.size(); i < n; ++i) {
0132         data_[i] += other.data_[i];
0133       }
0134       return *this;
0135     }
0136
0137     // always defined
0138     void fill(const T &val) { std::fill(data_.begin(), data_.end(), val); }
0139     void zero() { fill(T()); }
0140
0141     // defined only if T has a 'clear' method
0142     void clear() {
0143       for (T &t : data_)
0144         t.clear();
0145     }
0146
0147   private:
0148     const Grid *grid_;
0149     std::vector<T> data_;
0150     const T empty_;
0151   };
0152   typedef GridData<float> EtaPhiCenterGrid;
0153   typedef GridData<float> EtGrid;
0154   typedef GridData<int> IndexGrid;
0155
0156   struct PreCluster {
0157     PreCluster() : ptLocalMax(0), ptOverNeighLocalMaxSum(0) {}
0158     float ptLocalMax;              // pt if it's a local max, zero otherwise
0159     float ptOverNeighLocalMaxSum;  // pt / (sum of ptLocalMax of neighbours); zero if no neighbours
0160     void clear() { ptLocalMax = ptOverNeighLocalMaxSum = 0; }
0161   };
0162   typedef GridData<PreCluster> PreClusterGrid;
0163
0164   struct Cluster {
0165     Cluster() : et(0), eta(0), phi(0) {}
0166     float et, eta, phi;
0167     std::vector<std::pair<int, float>> constituents;
0168     void clear() {
0169       et = eta = phi = 0;
0170       constituents.clear();
0171     }
0172   };
0173
0174   struct CombinedCluster : public Cluster {
0175     float ecal_et, hcal_et;
0176     float ecal_eta, ecal_phi;
0177     void clear() {
0178       Cluster::clear();
0179       ecal_et = hcal_et = 0;
0180       ecal_eta = ecal_phi = 0;
0181     }
0182   };
0183
0184   const Grid *getGrid(const std::string &type);
0185
0186   class SingleCaloClusterer {
0187   public:
0188     SingleCaloClusterer(const edm::ParameterSet &pset);
0189     ~SingleCaloClusterer();
0190     void clear();
0191     void add(const reco::Candidate &c) { add(c.pt(), c.eta(), c.phi()); }
0192     void add(float pt, float eta, float phi) {
0193       rawet_(eta, phi) += pt;
0194       if (preciseEtaPhi_) {
0195         float newet = rawet_(eta, phi);
0196         float prevw = (newet - pt) / newet;
0197         float nextw = pt / newet;
0198         eta_center_(eta, phi) = eta_center_(eta, phi) * prevw + eta * nextw;
0199         phi_center_(eta, phi) = phi_center_(eta, phi) * prevw + phi * nextw;
0200       }
0201     }
0202     void run();
0203
0204     /// possibly grow clusters by adding unclustered energy on the sides
0205     //  note: there can be some double-counting as the same unclustered energy can go into more clusters
0206     void grow();
0207
0208     const EtGrid &raw() const { return rawet_; }
0209     const EtaPhiCenterGrid &etaCenter() const { return eta_center_; }
0210     const EtaPhiCenterGrid &phiCenter() const { return phi_center_; }
0211     const IndexGrid &indexGrid() const { return clusterIndex_; }
0212     const std::vector<Cluster> &clusters() const { return clusters_; }
0213     const Cluster &cluster(int i) const {
0214       return (i == -1 || clusterIndex_[i] == -1) ? nullCluster_ : clusters_[clusterIndex_[i]];
0215     }
0216
0217     /// non-const access to the energy: be careful to use it only before 'run()'
0218     EtGrid &raw() { return rawet_; }
0219     EtaPhiCenterGrid &etaCenter() { return eta_center_; }
0220     EtaPhiCenterGrid &phiCenter() { return phi_center_; }
0221
0222     // for the moment, generic interface that takes a cluster and returns the corrected pt
0223     template <typename Corrector>
0224     void correct(const Corrector &corrector) {
0225       for (Cluster &c : clusters_) {
0226         c.et = corrector(c);
0227       }
0228     }
0229
0230     std::unique_ptr<l1t::PFClusterCollection> fetchCells(bool unclusteredOnly = false, float ptMin = 0.) const;
0231
0232     std::unique_ptr<l1t::PFClusterCollection> fetch(float ptMin = 0.) const;
0233     std::unique_ptr<l1t::PFClusterCollection> fetch(const edm::OrphanHandle<l1t::PFClusterCollection> &cells,
0234                                                     float ptMin = 0.) const;
0235
0236   private:
0237     enum class EnergyShareAlgo {
0238       Fractions, /* each local maximum neighbour takes a share proportional to its value */
0239       None,      /* each local maximum neighbour takes all the value (double counting!) */
0240       Greedy,    /* assing cell to the highest local maximum neighbour */
0241       Crude
0242     }; /* if there's more than one local maximum neighbour, they all take half of the value (no fp division) */
0243     const Grid *grid_;
0244     EtGrid rawet_, unclustered_;
0245     EtaPhiCenterGrid eta_center_;
0246     EtaPhiCenterGrid phi_center_;
0247     PreClusterGrid precluster_;
0248     IndexGrid clusterIndex_, cellKey_;
0249     bool preciseEtaPhi_;
0250     std::vector<double> etaBounds_;
0251     std::vector<double> phiBounds_;
0252     std::vector<unsigned int> maxClustersEtaPhi_;  //eta x phi
0253     std::vector<Cluster> clusters_;
0254     const Cluster nullCluster_;
0255     float zsEt_, seedEt_, minClusterEt_, minEtToGrow_;
0256     EnergyShareAlgo energyShareAlgo_;
0257     bool energyWeightedPosition_;  // do the energy-weighted cluster position instead of the cell center
0258     std::vector<int> neighborCells_;
0259   };
0260
0261   class SimpleCaloLinkerBase {
0262   public:
0263     SimpleCaloLinkerBase(const edm::ParameterSet &pset,
0264                          const SingleCaloClusterer &ecal,
0265                          const SingleCaloClusterer &hcal);
0266     virtual ~SimpleCaloLinkerBase();
0267     virtual void clear() { clearBase(); }
0268     virtual void run() = 0;
0269     void clearBase() {
0270       clusters_.clear();
0271       clusterIndex_.fill(-1);
0272     }
0273
0274     // for the moment, generic interface that takes a cluster and returns the corrected pt
0275     template <typename Corrector>
0276     void correct(const Corrector &corrector) {
0277       for (CombinedCluster &c : clusters_) {
0278         c.et = corrector(c);
0279       }
0280     }
0281
0282     std::unique_ptr<l1t::PFClusterCollection> fetch() const;
0283     std::unique_ptr<l1t::PFClusterCollection> fetch(const edm::OrphanHandle<l1t::PFClusterCollection> &ecal,
0284                                                     const edm::OrphanHandle<l1t::PFClusterCollection> &hcal) const;
0285
0286   protected:
0287     const Grid *grid_;
0288     const SingleCaloClusterer &ecal_, &hcal_;
0289     IndexGrid clusterIndex_;
0290     std::vector<CombinedCluster> clusters_;
0291     std::vector<double> etaBounds_;
0292     std::vector<double> phiBounds_;
0293     std::vector<unsigned int> maxClustersEtaPhi_;  //eta x phi
0294     float hoeCut_, minPhotonEt_, minHadronRawEt_, minHadronEt_;
0295     bool noEmInHGC_;
0296   };
0297
0298   class SimpleCaloLinker : public SimpleCaloLinkerBase {
0299   public:
0300     SimpleCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer &ecal, const SingleCaloClusterer &hcal);
0301     ~SimpleCaloLinker() override;
0302     void clear() override;
0303     void run() override;
0304
0305   protected:
0306     PreClusterGrid ecalToHCal_;
0307   };
0308   class FlatCaloLinker : public SimpleCaloLinkerBase {
0309   public:
0310     FlatCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer &ecal, const SingleCaloClusterer &hcal);
0311     ~FlatCaloLinker() override;
0312     void clear() override;
0313     void run() override;
0314
0315   protected:
0316     SingleCaloClusterer combClusterer_;
0317   };
0318
0319   class CombinedCaloLinker : public SimpleCaloLinkerBase {
0320   public:
0321     CombinedCaloLinker(const edm::ParameterSet &pset, const SingleCaloClusterer &ecal, const SingleCaloClusterer &hcal);
0322     ~CombinedCaloLinker() override;
0323     void clear() override;
0324     void run() override;
0325
0326   protected:
0327     SingleCaloClusterer combClusterer_;
0328   };
0329
0330   class GridSelector {
0331   public:
0332     GridSelector(std::vector<double> etaBounds, std::vector<double> phiBounds, std::vector<unsigned int> maxClusters);
0333     ~GridSelector() {}
0334     void fill(float pt, float eta, float phi, unsigned int index);
0335     std::vector<unsigned int> returnSorted();
0336
0337   private:
0338     const std::vector<double> etaBounds_;
0339     const std::vector<double> phiBounds_;
0340     const std::vector<unsigned int> maxClustersEtaPhi_;
0341     std::vector<std::vector<std::pair<float, unsigned int>>> regionPtIndices_;  //pt and index pairs in each region
0342   };
0343
0344   // makes a calo linker (pointer will be owned by the callee)
0345   std::unique_ptr<SimpleCaloLinkerBase> makeCaloLinker(const edm::ParameterSet &pset,
0346                                                        const SingleCaloClusterer &ecal,
0347                                                        const SingleCaloClusterer &hcal);
0348
0349 }  // namespace l1tpf_calo
0350
0351 #endif