HGCalRecProducers/interface/HGCalLayerTiles.h

0001 // Authors: Felice Pantaleo - felice.pantaleo@cern.ch, Olivie Abigail Franklova - olivie.abigail.franklova@cern.ch
0002 // Date: 03/2019
0003
0004 #ifndef RecoLocalCalo_HGCalRecProducers_HGCalLayerTiles_h
0005 #define RecoLocalCalo_HGCalRecProducers_HGCalLayerTiles_h
0006
0007 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h"
0008 #include "RecoLocalCalo/HGCalRecProducers/interface/HFNoseTilesConstants.h"
0009 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesWrapper.h"
0010 #include "DataFormats/Math/interface/normalizedPhi.h"
0011
0012 #include <vector>
0013 #include <array>
0014 #include <cmath>
0015 #include <algorithm>
0016 #include <cassert>
0017
0018 template <typename T, typename WRAPPER>
0019 class HGCalLayerTilesT {
0020 public:
0021   typedef T type;
0022   /**
0023      * @brief fill the tile
0024      *
0025      * @param[in] dim1 represents x or eta
0026      * @param[in] dim2 represents y or phils
0027      *
0028     */
0029   void fill(const std::vector<float>& dim1, const std::vector<float>& dim2) {
0030     auto cellsSize = dim1.size();
0031     for (unsigned int i = 0; i < cellsSize; ++i) {
0032       auto idx = getGlobalBin(dim1[i], dim2[i]);
0033       tiles_[idx].push_back(i);
0034     }
0035   }
0036   /**
0037     * @brief compute bin for dim1 (x or eta)
0038     *
0039     * @param[in] dim for bining
0040     * @return computed bin
0041     */
0042   int getDim1Bin(float dim) const {
0043     constexpr float dimRange = T::maxDim1 - T::minDim1;
0044     static_assert(dimRange >= 0.);
0045     constexpr float r = T::nColumns / dimRange;
0046     int dimBin = (dim - T::minDim1) * r;
0047     dimBin = std::clamp(dimBin, 0, T::nColumns - 1);
0048     return dimBin;
0049   }
0050
0051   /**
0052     * @brief compute bin for dim2 (y or phi)
0053     *
0054     * @param[in] dim for bining
0055     * @return computed bin
0056     */
0057   int getDim2Bin(float dim2) const {
0058     if constexpr (std::is_same_v<WRAPPER, NoPhiWrapper>) {
0059       constexpr float dimRange = T::maxDim2 - T::minDim2;
0060       static_assert(dimRange >= 0.);
0061       constexpr float r = T::nRows / dimRange;
0062       int dimBin = (dim2 - T::minDim2) * r;
0063       dimBin = std::clamp(dimBin, 0, T::nRows - 1);
0064       return dimBin;
0065     } else {
0066       auto normPhi = normalizedPhi(dim2);
0067       constexpr float r = T::nRows * M_1_PI * 0.5f;
0068       int phiBin = (normPhi + M_PI) * r;
0069       return phiBin;
0070     }
0071   }
0072
0073   int mPiPhiBin = getDim2Bin(-M_PI);
0074   int pPiPhiBin = getDim2Bin(M_PI);
0075
0076   inline float distance2(float dim1Cell1, float dim2Cell1, float dim1Cell2, float dim2Cell2) const {  // distance squared
0077     float d1 = dim1Cell1 - dim1Cell2;
0078     float d2 = dim2Cell1 - dim2Cell2;
0079     if constexpr (std::is_same_v<WRAPPER, PhiWrapper>) {
0080       d2 = reco::deltaPhi(dim2Cell1, dim2Cell2);
0081     }
0082     return std::fmaf(d1, d1, d2 * d2);
0083   }
0084   int getGlobalBin(float dim1, float dim2) const { return getDim1Bin(dim1) + getDim2Bin(dim2) * T::nColumns; }
0085
0086   int getGlobalBinByBin(int dim1Bin, int dim2Bin) const { return dim1Bin + dim2Bin * T::nColumns; }
0087
0088   std::array<int, 4> searchBox(float dim1Min, float dim1Max, float dim2Min, float dim2Max) const {
0089     if constexpr (std::is_same_v<WRAPPER, PhiWrapper>) {
0090       if (dim1Max - dim1Min < 0) {
0091         return std::array<int, 4>({{0, 0, 0, 0}});
0092       }
0093     }
0094     int dim1BinMin = getDim1Bin(dim1Min);
0095     int dim1BinMax = getDim1Bin(dim1Max);
0096     int dim2BinMin = getDim2Bin(dim2Min);
0097     int dim2BinMax = getDim2Bin(dim2Max);
0098     if constexpr (std::is_same_v<WRAPPER, PhiWrapper>) {
0099       // If the search window cross the phi-bin boundary, add T::nPhiBins to the
0100       // MAx value. This guarantees that the caller can perform a valid doule
0101       // loop on eta and phi. It is the caller responsibility to perform a module
0102       // operation on the phiBin values returned by this function, to explore the
0103       // correct bins.
0104       if (dim2BinMax < dim2BinMin) {
0105         dim2BinMax += T::nRows;
0106       }
0107     }
0108     return std::array<int, 4>({{dim1BinMin, dim1BinMax, dim2BinMin, dim2BinMax}});
0109   }
0110
0111   void clear() {
0112     for (auto& t : tiles_)
0113       t.clear();
0114   }
0115
0116   const std::vector<int>& operator[](int globalBinId) const { return tiles_[globalBinId]; }
0117
0118 private:
0119   std::array<std::vector<int>, T::nTiles> tiles_;
0120 };
0121
0122 using HGCalSiliconLayerTiles = HGCalLayerTilesT<HGCalSiliconTilesConstants, NoPhiWrapper>;
0123 using HGCalScintillatorLayerTiles = HGCalLayerTilesT<HGCalScintillatorTilesConstants, PhiWrapper>;
0124 using HFNoseLayerTiles = HGCalLayerTilesT<HFNoseTilesConstants, NoPhiWrapper>;
0125 #endif