Merge pull request #37555 from rovere/EtaPhiSearchInTile

Eta phi search in tile

DataFormats/HGCalReco/interface/TICLLayerTile.h

@@ -38,6 +38,14 @@ class TICLLayerTileT {
     int etaBinMax = etaBin(etaMax);
     int phiBinMin = phiBin(phiMin);
     int phiBinMax = phiBin(phiMax);
+    // If the search window cross the phi-bin boundary, add T::nPhiBins to the
+    // MAx value. This guarantees that the caller can perform a valid doule
+    // loop on eta and phi. It is the caller responsibility to perform a module
+    // operation on the phiBin values returned by this function, to explore the
+    // correct bins.
+    if (phiBinMax < phiBinMin) {
+      phiBinMax += T::nPhiBins;
+    }
     return std::array<int, 4>({{etaBinMin, etaBinMax, phiBinMin, phiBinMax}});
   }
 

DataFormats/HGCalReco/test/BuildFile.xml

@@ -0,0 +1,3 @@
+<bin   name="EtaPhiSearchInTile" file="EtaPhiSearchInTile_t.cpp">
+  <use name="DataFormats/HGCalReco" />
+</bin>

DataFormats/HGCalReco/test/EtaPhiSearchInTile_t.cpp

@@ -0,0 +1,60 @@
+#include <cassert>
+#include <cmath>
+#include <iostream>
+
+#include "DataFormats/HGCalReco/interface/TICLLayerTile.h"
+
+using namespace ticl;
+
+void runTest(TICLLayerTile const& t, int expected, float etaMin, float etaMax, float phiMin, float phiMax) {
+  auto limits = t.searchBoxEtaPhi(etaMin, etaMax, phiMin, phiMax);
+
+  std::cout << "Limits are: " << limits[0] << " " << limits[1] << " " << limits[2] << " " << limits[3] << std::endl;
+  assert(limits[0] <= limits[1]);
+  assert(limits[2] <= limits[3]);
+
+  int entries = 0;
+  for (int e = limits[0]; e <= limits[1]; ++e) {
+    for (int p = limits[2]; p <= limits[3]; ++p) {
+      int phi = (p % TICLLayerTile::type::nPhiBins);
+      auto global_bin = t.globalBin(e, phi);
+      entries += t[global_bin].size();
+    }
+  }
+
+  std::cout << "Found " << entries << " entries, expected " << expected << std::endl;
+  assert(entries == expected);
+}
+
+int main(int argc, char* argv[]) {
+  auto constexpr phiBins = TICLLayerTile::type::nPhiBins;
+  auto constexpr phi_bin_width = 2. * M_PI / phiBins;
+  auto constexpr phi_transition_left = M_PI - phi_bin_width;
+  auto constexpr phi_transition_right = M_PI + phi_bin_width;
+  auto constexpr phi_transition_right2 = -M_PI + 3. * phi_bin_width;
+  unsigned int constexpr entries_left = 11;
+  unsigned int constexpr entries_right = 7;
+  float constexpr eta = 2.0;
+
+  TICLLayerTile t, t2;
+  std::cout << "Testing a Tile with " << phiBins << " bins with binwidth: " << phi_bin_width << " at bin transition"
+            << std::endl;
+  std::cout << "Filling left-pi bin: " << t.phiBin(phi_transition_left) << std::endl;
+  std::cout << "Filling right-pi bin: " << t.phiBin(phi_transition_right) << std::endl;
+  std::cout << "Filling right2-pi bin: " << t.phiBin(phi_transition_right2) << std::endl;
+
+  for (unsigned int i = 0; i < entries_left; ++i) {
+    t.fill(eta, phi_transition_left, i);
+    t2.fill(eta, phi_transition_left, i);
+  }
+
+  for (unsigned int i = 0; i < entries_right; ++i) {
+    t.fill(eta, phi_transition_right, i);
+    t2.fill(eta, phi_transition_right2, i);
+  }
+
+  runTest(t, entries_left + entries_right, 1.95, 2.05, phi_transition_left, phi_transition_right);
+  runTest(t2, entries_left + entries_right, 1.95, 2.05, phi_transition_left, phi_transition_right2);
+
+  return 0;
+}

RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h

@@ -6,6 +6,7 @@
 
 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h"
 #include "RecoLocalCalo/HGCalRecProducers/interface/HFNoseTilesConstants.h"
+#include "DataFormats/Math/interface/normalizedPhi.h"
 
 #include <vector>
 #include <array>
@@ -16,6 +17,7 @@
 template <typename T>
 class HGCalLayerTilesT {
 public:
+  typedef T type;
   void fill(const std::vector<float>& x,
             const std::vector<float>& y,
             const std::vector<float>& eta,
@@ -26,14 +28,6 @@ class HGCalLayerTilesT {
       tiles_[getGlobalBin(x[i], y[i])].push_back(i);
       if (!isSi[i]) {
         tiles_[getGlobalBinEtaPhi(eta[i], phi[i])].push_back(i);
-        // Copy cells in phi=[-3.15,-3.] to the last bin
-        if (getPhiBin(phi[i]) == mPiPhiBin) {
-          tiles_[getGlobalBinEtaPhi(eta[i], phi[i] + 2 * M_PI)].push_back(i);
-        }
-        // Copy cells in phi=[3.,3.15] to the first bin
-        if (getPhiBin(phi[i]) == pPiPhiBin) {
-          tiles_[getGlobalBinEtaPhi(eta[i], phi[i] - 2 * M_PI)].push_back(i);
-        }
       }
     }
   }
@@ -66,11 +60,9 @@ class HGCalLayerTilesT {
   }
 
   int getPhiBin(float phi) const {
-    constexpr float phiRange = T::maxPhi - T::minPhi;
-    static_assert(phiRange >= 0.);
-    constexpr float r = T::nRowsPhi / phiRange;
-    int phiBin = (phi - T::minPhi) * r;
-    phiBin = std::clamp(phiBin, 0, T::nRowsPhi - 1);
+    auto normPhi = normalizedPhi(phi);
+    constexpr float r = T::nRowsPhi * M_1_PI * 0.5f;
+    int phiBin = (normPhi + M_PI) * r;
     return phiBin;
   }
 
@@ -102,6 +94,15 @@ class HGCalLayerTilesT {
     int etaBinMax = getEtaBin(etaMax);
     int phiBinMin = getPhiBin(phiMin);
     int phiBinMax = getPhiBin(phiMax);
+    // If the search window cross the phi-bin boundary, add T::nPhiBins to the
+    // MAx value. This guarantees that the caller can perform a valid doule
+    // loop on eta and phi. It is the caller responsibility to perform a module
+    // operation on the phiBin values returned by this function, to explore the
+    // correct bins.
+    if (phiBinMax < phiBinMin) {
+      phiBinMax += T::nRowsPhi;
+    }
+
     return std::array<int, 4>({{etaBinMin, etaBinMax, phiBinMin, phiBinMax}});
   }
 

RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h

@@ -5,6 +5,7 @@
 #define RecoLocalCalo_HGCalRecProducer_HGCalTilesConstants_h
 
 #include "DataFormats/Math/interface/constexpr_cmath.h"
+#include <cmath>
 #include <cstdint>
 #include <array>
 
@@ -19,12 +20,8 @@ struct HGCalTilesConstants {
   static constexpr float tileSizeEtaPhi = 0.15f;
   static constexpr float minEta = -3.f;
   static constexpr float maxEta = 3.f;
-  //To properly construct search box for cells in phi=[-3.15,-3.] and [3.,3.15], cells in phi=[3.,3.15]
-  //are copied to the first bin and cells in phi=[-3.15,-3.] are copied to the last bin
-  static constexpr float minPhi = -3.3f;
-  static constexpr float maxPhi = 3.3f;
   static constexpr int nColumnsEta = reco::ceil((maxEta - minEta) / tileSizeEtaPhi);
-  static constexpr int nRowsPhi = reco::ceil((maxPhi - minPhi) / tileSizeEtaPhi);
+  static constexpr int nRowsPhi = reco::ceil(2. * M_PI / tileSizeEtaPhi);
   static constexpr int nTiles = nColumns * nRows + nColumnsEta * nRowsPhi;
 };
 

RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc

@@ -285,7 +285,8 @@ void HGCalCLUEAlgoT<T>::calculateLocalDensity(const T& lt, const unsigned int la
 
       for (int etaBin = search_box[0]; etaBin < search_box[1] + 1; ++etaBin) {
         for (int phiBin = search_box[2]; phiBin < search_box[3] + 1; ++phiBin) {
-          int binId = lt.getGlobalBinByBinEtaPhi(etaBin, phiBin);
+          int phi = (phiBin % T::type::nRowsPhi);
+          int binId = lt.getGlobalBinByBinEtaPhi(etaBin, phi);
           size_t binSize = lt[binId].size();
 
           for (unsigned int j = 0; j < binSize; j++) {
@@ -416,7 +417,8 @@ void HGCalCLUEAlgoT<T>::calculateDistanceToHigher(const T& lt,
       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
           // get the id of this bin
-          size_t binId = lt.getGlobalBinByBinEtaPhi(xBin, yBin);
+          int phi = (yBin % T::type::nRowsPhi);
+          size_t binId = lt.getGlobalBinByBinEtaPhi(xBin, phi);
           // get the size of this bin
           size_t binSize = lt[binId].size();
 

RecoLocalCalo/HGCalRecProducers/test/BuildFile.xml

@@ -0,0 +1,3 @@
+<bin   name="EtaPhiSearchInTileLC" file="EtaPhiSearchInTile_t.cpp">
+  <use name="RecoLocalCalo/HGCalRecProducers" />
+</bin>

RecoLocalCalo/HGCalRecProducers/test/EtaPhiSearchInTile_t.cpp

@@ -0,0 +1,71 @@
+#include <cassert>
+#include <cmath>
+#include <iostream>
+
+#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h"
+
+void runTest(HGCalLayerTiles const& t, int expected, float etaMin, float etaMax, float phiMin, float phiMax) {
+  auto limits = t.searchBoxEtaPhi(etaMin, etaMax, phiMin, phiMax);
+
+  std::cout << "Limits are: " << limits[0] << " " << limits[1] << " " << limits[2] << " " << limits[3] << std::endl;
+  assert(limits[0] <= limits[1]);
+  assert(limits[2] <= limits[3]);
+
+  int entries = 0;
+  for (int e = limits[0]; e <= limits[1]; ++e) {
+    for (int p = limits[2]; p <= limits[3]; ++p) {
+      int phi = (p % HGCalLayerTiles::type::nRowsPhi);
+      auto global_bin = t.getGlobalBinByBinEtaPhi(e, phi);
+      entries += t[global_bin].size();
+    }
+  }
+
+  std::cout << "Found " << entries << " entries, expected " << expected << std::endl;
+  assert(entries == expected);
+}
+
+int main(int argc, char* argv[]) {
+  auto constexpr phiBins = HGCalLayerTiles::type::nRowsPhi;
+  auto constexpr phi_bin_width = 2. * M_PI / phiBins;
+  auto constexpr phi_transition_left = M_PI - 3. * phi_bin_width;
+  auto constexpr phi_transition_right = M_PI + 3. * phi_bin_width;
+  auto constexpr phi_transition_right2 = -M_PI + 5. * phi_bin_width;
+  unsigned int constexpr entries_left = 11;
+  unsigned int constexpr entries_right = 7;
+  float constexpr eta = 2.0;
+
+  // Force filling using eta/phi vectors
+  std::vector<bool> isSilicon(entries_left + entries_right, false);
+  std::vector<float> dummy(entries_left + entries_right, 0.);
+  std::vector<float> etas(entries_left + entries_right, eta);
+  std::vector<float> phis_left(entries_left, phi_transition_left);
+  std::vector<float> phis_right(entries_right, phi_transition_right);
+  std::vector<float> phis_right2(entries_right, phi_transition_right2);
+  std::vector<float> phis;
+  std::vector<float> phis2;
+  phis.reserve(entries_left + entries_right);
+  phis.insert(phis.end(), phis_left.begin(), phis_left.end());
+  phis.insert(phis.end(), phis_right.begin(), phis_right.end());
+  phis2.reserve(entries_left + entries_right);
+  phis2.insert(phis2.end(), phis_left.begin(), phis_left.end());
+  phis2.insert(phis2.end(), phis_right2.begin(), phis_right2.end());
+
+  HGCalLayerTiles t, t2;
+  t.fill(dummy, dummy, etas, phis, isSilicon);
+  t2.fill(dummy, dummy, etas, phis2, isSilicon);
+
+  std::cout << "Testing a Tile with " << phiBins << " bins with binwidth: " << phi_bin_width << " at pi transition"
+            << std::endl;
+  std::cout << "-M_PI bin: " << t.mPiPhiBin << " M_PI bin: " << t.pPiPhiBin << std::endl;
+  std::cout << "Filling phi value: " << phi_transition_left << " at left-pi bin: " << t.getPhiBin(phi_transition_left)
+            << std::endl;
+  std::cout << "Filling phi value: " << phi_transition_right
+            << " at right-pi bin: " << t.getPhiBin(phi_transition_right) << std::endl;
+  std::cout << "Filling phi value: " << phi_transition_right2
+            << " at right-pi bin: " << t.getPhiBin(phi_transition_right2) << std::endl;
+
+  runTest(t, entries_right + entries_left, 1.95, 2.05, phi_transition_left, phi_transition_right);
+  runTest(t2, entries_right + entries_left, 1.95, 2.05, phi_transition_left, phi_transition_right2);
+
+  return 0;
+}