Merge pull request #46857 from SegmentLinking/CMSSW_14_1_0_pre3_LST_X…

…_LSTCore_realfiles_batch1_devel_rebased LST algorithm developments accumulated during integration PR
cms-sw · Dec 4, 2024 · 9628e8f · 9628e8f
2 parents 15e4fe7 + 8c9753d
commit 9628e8f
Show file tree

Hide file tree

Showing 35 changed files with 3,514 additions and 2,616 deletions.
diff --git a/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc b/RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc
@@ -13,16 +13,23 @@
 namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
   class LSTModulesDevESProducer : public ESProducer {
+  private:
+    std::string ptCutLabel_;
+
   public:
-    LSTModulesDevESProducer(edm::ParameterSet const& iConfig) : ESProducer(iConfig) { setWhatProduced(this); }
+    LSTModulesDevESProducer(edm::ParameterSet const& iConfig)
+        : ESProducer(iConfig), ptCutLabel_(iConfig.getParameter<std::string>("ptCutLabel")) {
+      setWhatProduced(this, ptCutLabel_);
+    }
 
     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
       edm::ParameterSetDescription desc;
+      desc.add<std::string>("ptCutLabel", "0.8");
       descriptions.addWithDefaultLabel(desc);
     }
 
     std::unique_ptr<lst::LSTESData<DevHost>> produce(TrackerRecoGeometryRecord const& iRecord) {
-      return lst::loadAndFillESHost();
+      return lst::loadAndFillESHost(ptCutLabel_);
     }
   };
 

diff --git a/RecoTracker/LST/plugins/alpaka/LSTProducer.cc b/RecoTracker/LST/plugins/alpaka/LSTProducer.cc
@@ -28,8 +28,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     LSTProducer(edm::ParameterSet const& config)
         : lstPixelSeedInputToken_{consumes(config.getParameter<edm::InputTag>("pixelSeedInput"))},
           lstPhase2OTHitsInputToken_{consumes(config.getParameter<edm::InputTag>("phase2OTHitsInput"))},
-          lstESToken_{esConsumes()},
+          lstESToken_{esConsumes(edm::ESInputTag("", config.getParameter<std::string>("ptCutLabel")))},
           verbose_(config.getParameter<bool>("verbose")),
+          ptCut_(config.getParameter<double>("ptCut")),
           nopLSDupClean_(config.getParameter<bool>("nopLSDupClean")),
           tcpLSTriplets_(config.getParameter<bool>("tcpLSTriplets")),
           lstOutputToken_{produces()} {}
@@ -43,6 +44,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
       lst_.run(event.queue(),
                verbose_,
+               static_cast<float>(ptCut_),
                &lstESDeviceData,
                pixelSeeds.px(),
                pixelSeeds.py(),
@@ -78,6 +80,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
       desc.add<edm::InputTag>("pixelSeedInput", edm::InputTag{"lstPixelSeedInputProducer"});
       desc.add<edm::InputTag>("phase2OTHitsInput", edm::InputTag{"lstPhase2OTHitsInputProducer"});
       desc.add<bool>("verbose", false);
+      desc.add<double>("ptCut", 0.8);
+      desc.add<std::string>("ptCutLabel", "0.8");
       desc.add<bool>("nopLSDupClean", false);
       desc.add<bool>("tcpLSTriplets", false);
       descriptions.addWithDefaultLabel(desc);
@@ -87,7 +91,10 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     edm::EDGetTokenT<LSTPixelSeedInput> lstPixelSeedInputToken_;
     edm::EDGetTokenT<LSTPhase2OTHitsInput> lstPhase2OTHitsInputToken_;
     device::ESGetToken<lst::LSTESData<Device>, TrackerRecoGeometryRecord> lstESToken_;
-    const bool verbose_, nopLSDupClean_, tcpLSTriplets_;
+    const bool verbose_;
+    const double ptCut_;
+    const bool nopLSDupClean_;
+    const bool tcpLSTriplets_;
     edm::EDPutTokenT<LSTOutput> lstOutputToken_;
 
     lst::LST lst_;

diff --git a/RecoTracker/LSTCore/BuildFile.xml b/RecoTracker/LSTCore/BuildFile.xml
@@ -2,7 +2,6 @@
 <use name="boost_header"/>
 <use name="root"/>
 <use name="HeterogeneousCore/AlpakaInterface"/>
-<flags CXXFLAGS="-DPT_CUT=0.8"/>
 <use name="RecoTracker/LSTCore" for="alpaka"/>
 <flags ALPAKA_BACKENDS="1"/>
 <export>

diff --git a/RecoTracker/LSTCore/interface/Common.h b/RecoTracker/LSTCore/interface/Common.h
@@ -20,11 +20,6 @@ namespace lst {
   // Named types for LST objects
   enum LSTObjType { T5 = 4, pT3 = 5, pT5 = 7, pLS = 8 };
 
-// If a compile time flag does not define PT_CUT, default to 0.8 (GeV)
-#ifndef PT_CUT
-  constexpr float PT_CUT = 0.8f;
-#endif
-
   constexpr unsigned int max_blocks = 80;
   constexpr unsigned int max_connected_modules = 40;
 

diff --git a/RecoTracker/LSTCore/interface/LSTESData.h b/RecoTracker/LSTCore/interface/LSTESData.h
@@ -40,7 +40,7 @@ namespace lst {
           pixelMapping(pixelMappingIn) {}
   };
 
-  std::unique_ptr<LSTESData<alpaka_common::DevHost>> loadAndFillESHost();
+  std::unique_ptr<LSTESData<alpaka_common::DevHost>> loadAndFillESHost(std::string& ptCutLabel);
 
 }  // namespace lst
 

diff --git a/RecoTracker/LSTCore/interface/QuintupletsSoA.h b/RecoTracker/LSTCore/interface/QuintupletsSoA.h
@@ -28,7 +28,9 @@ namespace lst {
                       SOA_COLUMN(float, regressionF),
                       SOA_COLUMN(float, rzChiSquared),  // r-z only chi2
                       SOA_COLUMN(float, chiSquared),
-                      SOA_COLUMN(float, nonAnchorChiSquared));
+                      SOA_COLUMN(float, nonAnchorChiSquared),
+                      SOA_COLUMN(float, dBeta1),
+                      SOA_COLUMN(float, dBeta2));
 
   using QuintupletsSoA = QuintupletsSoALayout<>;
   using Quintuplets = QuintupletsSoA::View;

diff --git a/RecoTracker/LSTCore/interface/alpaka/Common.h b/RecoTracker/LSTCore/interface/alpaka/Common.h
@@ -49,7 +49,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float k2Rinv1GeVf = (2.99792458e-3 * 3.8) / 2;
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kR1GeVf = 1. / (2.99792458e-3 * 3.8);
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kSinAlphaMax = 0.95;
-  ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float ptCut = PT_CUT;
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kDeltaZLum = 15.0;
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPixelPSZpitch = 0.15;
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kStripPSZpitch = 2.4;
@@ -60,22 +59,27 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
   // To be updated with std::numeric_limits<float>::infinity() in the code and data files
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kVerticalModuleSlope = 123456789.0;
 
-  namespace t5dnn {
-
-    // Working points matching LST fake rate (43.9%) or signal acceptance (82.0%)
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kLSTWp1 = 0.3418833f;  // 94.0% TPR, 43.9% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kLSTWp2 = 0.6177366f;  // 82.0% TPR, 20.0% FPR
-    // Other working points
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp70 = 0.7776195f;    // 70.0% TPR, 10.0% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp75 = 0.7181118f;    // 75.0% TPR, 13.5% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp80 = 0.6492643f;    // 80.0% TPR, 17.9% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp85 = 0.5655319f;    // 85.0% TPR, 23.8% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp90 = 0.4592205f;    // 90.0% TPR, 32.6% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp95 = 0.3073708f;    // 95.0% TPR, 47.7% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp97p5 = 0.2001348f;  // 97.5% TPR, 61.2% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp99 = 0.1120605f;    // 99.0% TPR, 75.9% FPR
-    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp99p9 = 0.0218196f;  // 99.9% TPR, 95.4% FPR
+  ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMiniDeltaTilted[3] = {0.26f, 0.26f, 0.26f};
+  ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMiniDeltaFlat[6] = {0.26f, 0.16f, 0.16f, 0.18f, 0.18f, 0.18f};
+  ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMiniDeltaLooseTilted[3] = {0.4f, 0.4f, 0.4f};
+  ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMiniDeltaEndcap[5][15] = {
+      {0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
+      {0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
+      {0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
+      {0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
+      {0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f}};
 
+  namespace t5dnn {
+    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kZ_max = 267.2349854f;
+    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kR_max = 110.1099396f;
+    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kEta_norm = 2.5f;
+    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPhi_norm = kPi;
+    // pt, eta binned
+    constexpr unsigned int kPtBins = 2;
+    constexpr unsigned int kEtaBins = 10;
+    ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWp[kPtBins][kEtaBins] = {
+        {0.4493, 0.4939, 0.5715, 0.6488, 0.5709, 0.5938, 0.7164, 0.7565, 0.8103, 0.8593},
+        {0.4488, 0.4448, 0.5067, 0.5929, 0.4836, 0.4112, 0.4968, 0.4403, 0.5597, 0.5067}};
   }  // namespace t5dnn
 
 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst

diff --git a/RecoTracker/LSTCore/interface/alpaka/LST.h b/RecoTracker/LSTCore/interface/alpaka/LST.h
@@ -17,6 +17,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
 
     void run(Queue& queue,
              bool verbose,
+             const float ptCut,
              LSTESData<Device> const* deviceESData,
              std::vector<float> const& see_px,
              std::vector<float> const& see_py,
@@ -63,7 +64,8 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
                       std::vector<unsigned int> const& ph2_detId,
                       std::vector<float> const& ph2_x,
                       std::vector<float> const& ph2_y,
-                      std::vector<float> const& ph2_z);
+                      std::vector<float> const& ph2_z,
+                      const float ptCut);
 
     void getOutput(LSTEvent& event);
 

diff --git a/RecoTracker/LSTCore/src/LSTESData.cc b/RecoTracker/LSTCore/src/LSTESData.cc
@@ -10,25 +10,25 @@
 
 namespace {
   std::string geometryDataDir() {
-    const char* path_lst_base = std::getenv("LST_BASE");
+    std::string path_str, path;
     const char* path_tracklooperdir = std::getenv("TRACKLOOPERDIR");
-    std::string path_str;
-    if (path_lst_base != nullptr) {
-      path_str = path_lst_base;
-    } else if (path_tracklooperdir != nullptr) {
+    std::stringstream search_path;
+    search_path << std::getenv("CMSSW_SEARCH_PATH");
+
+    while (std::getline(search_path, path, ':')) {
+      if (std::filesystem::exists(path + "/RecoTracker/LSTCore/data")) {
+        path_str = path;
+        break;
+      }
+    }
+
+    if (path_str.empty()) {
       path_str = path_tracklooperdir;
-      path_str += "/../";
+      path_str += "/..";
     } else {
-      std::stringstream search_path(std::getenv("CMSSW_SEARCH_PATH"));
-      std::string path;
-      while (std::getline(search_path, path, ':')) {
-        if (std::filesystem::exists(path + "/RecoTracker/LSTCore/data")) {
-          path_str = path;
-          break;
-        }
-      }
       path_str += "/RecoTracker/LSTCore";
     }
+
     return path_str;
   }
 
@@ -43,21 +43,22 @@ namespace {
   void loadMapsHost(lst::MapPLStoLayer& pLStoLayer,
                     lst::EndcapGeometry& endcapGeometry,
                     lst::TiltedGeometry& tiltedGeometry,
-                    lst::ModuleConnectionMap& moduleConnectionMap) {
+                    lst::ModuleConnectionMap& moduleConnectionMap,
+                    std::string& ptCutLabel) {
     // Module orientation information (DrDz or phi angles)
-    auto endcap_geom =
-        get_absolute_path_after_check_file_exists(geometryDataDir() + "/data/OT800_IT615_pt0.8/endcap_orientation.bin");
-    auto tilted_geom = get_absolute_path_after_check_file_exists(
-        geometryDataDir() + "/data/OT800_IT615_pt0.8/tilted_barrel_orientation.bin");
+    auto endcap_geom = get_absolute_path_after_check_file_exists(geometryDataDir() + "/data/OT800_IT615_pt" +
+                                                                 ptCutLabel + "/endcap_orientation.bin");
+    auto tilted_geom = get_absolute_path_after_check_file_exists(geometryDataDir() + "/data/OT800_IT615_pt" +
+                                                                 ptCutLabel + "/tilted_barrel_orientation.bin");
     // Module connection map (for line segment building)
-    auto mappath = get_absolute_path_after_check_file_exists(
-        geometryDataDir() + "/data/OT800_IT615_pt0.8/module_connection_tracing_merged.bin");
+    auto mappath = get_absolute_path_after_check_file_exists(geometryDataDir() + "/data/OT800_IT615_pt" + ptCutLabel +
+                                                             "/module_connection_tracing_merged.bin");
 
     endcapGeometry.load(endcap_geom);
     tiltedGeometry.load(tilted_geom);
     moduleConnectionMap.load(mappath);
 
-    auto pLSMapDir = geometryDataDir() + "/data/OT800_IT615_pt0.8/pixelmap/pLS_map";
+    auto pLSMapDir = geometryDataDir() + "/data/OT800_IT615_pt" + ptCutLabel + "/pixelmap/pLS_map";
     const std::array<std::string, 4> connects{
         {"_layer1_subdet5", "_layer2_subdet5", "_layer1_subdet4", "_layer2_subdet4"}};
     std::string path;
@@ -78,7 +79,7 @@ namespace {
   }
 }  // namespace
 
-std::unique_ptr<lst::LSTESData<alpaka_common::DevHost>> lst::loadAndFillESHost() {
+std::unique_ptr<lst::LSTESData<alpaka_common::DevHost>> lst::loadAndFillESHost(std::string& ptCutLabel) {
   uint16_t nModules;
   uint16_t nLowerModules;
   unsigned int nPixels;
@@ -87,7 +88,7 @@ std::unique_ptr<lst::LSTESData<alpaka_common::DevHost>> lst::loadAndFillESHost()
   TiltedGeometry tiltedGeometry;
   PixelMap pixelMapping;
   ModuleConnectionMap moduleConnectionMap;
-  ::loadMapsHost(pLStoLayer, endcapGeometry, tiltedGeometry, moduleConnectionMap);
+  ::loadMapsHost(pLStoLayer, endcapGeometry, tiltedGeometry, moduleConnectionMap, ptCutLabel);
 
   auto endcapGeometryDev =
       std::make_shared<EndcapGeometryDevHostCollection>(endcapGeometry.nEndCapMap, cms::alpakatools::host());
@@ -98,8 +99,8 @@ std::unique_ptr<lst::LSTESData<alpaka_common::DevHost>> lst::loadAndFillESHost()
               endcapGeometry.geoMapPhi_buf.data(),
               endcapGeometry.nEndCapMap * sizeof(float));
 
-  auto path =
-      get_absolute_path_after_check_file_exists(geometryDataDir() + "/data/OT800_IT615_pt0.8/sensor_centroids.bin");
+  auto path = get_absolute_path_after_check_file_exists(geometryDataDir() + "/data/OT800_IT615_pt" + ptCutLabel +
+                                                        "/sensor_centroids.bin");
   auto modulesBuffers = lst::loadModulesFromFile(pLStoLayer,
                                                  path.c_str(),
                                                  nModules,

diff --git a/RecoTracker/LSTCore/src/alpaka/Kernels.h b/RecoTracker/LSTCore/src/alpaka/Kernels.h
@@ -217,15 +217,22 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
 
           for (unsigned int ix1 = 0; ix1 < nQuintuplets_lowmod1; ix1 += 1) {
             unsigned int ix = quintupletModuleIndices_lowmod1 + ix1;
-            if (quintuplets.partOfPT5()[ix] || (quintuplets.isDup()[ix] & 1))
+            if (quintuplets.isDup()[ix] & 1)
               continue;
 
+            bool isPT5_ix = quintuplets.partOfPT5()[ix];
+
             for (unsigned int jx1 = 0; jx1 < nQuintuplets_lowmod2; jx1++) {
               unsigned int jx = quintupletModuleIndices_lowmod2 + jx1;
               if (ix == jx)
                 continue;
 
-              if (quintuplets.partOfPT5()[jx] || (quintuplets.isDup()[jx] & 1))
+              if (quintuplets.isDup()[jx] & 1)
+                continue;
+
+              bool isPT5_jx = quintuplets.partOfPT5()[jx];
+
+              if (isPT5_ix && isPT5_jx)
                 continue;
 
               float eta1 = __H2F(quintuplets.eta()[ix]);
@@ -249,9 +256,9 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
               int nMatched = checkHitsT5(ix, jx, quintuplets);
               const int minNHitsForDup_T5 = 5;
               if (dR2 < 0.001f || nMatched >= minNHitsForDup_T5) {
-                if (score_rphisum1 > score_rphisum2) {
+                if (isPT5_jx || score_rphisum1 > score_rphisum2) {
                   rmQuintupletFromMemory(quintuplets, ix, true);
-                } else if (score_rphisum1 < score_rphisum2) {
+                } else if (isPT5_ix || score_rphisum1 < score_rphisum2) {
                   rmQuintupletFromMemory(quintuplets, jx, true);
                 } else {
                   rmQuintupletFromMemory(quintuplets, (ix < jx ? ix : jx), true);

diff --git a/RecoTracker/LSTCore/src/alpaka/LST.cc b/RecoTracker/LSTCore/src/alpaka/LST.cc
@@ -75,7 +75,8 @@ void LST::prepareInput(std::vector<float> const& see_px,
                        std::vector<unsigned int> const& ph2_detId,
                        std::vector<float> const& ph2_x,
                        std::vector<float> const& ph2_y,
-                       std::vector<float> const& ph2_z) {
+                       std::vector<float> const& ph2_z,
+                       float const ptCut) {
   in_trkX_.clear();
   in_trkY_.clear();
   in_trkZ_.clear();
@@ -132,7 +133,7 @@ void LST::prepareInput(std::vector<float> const& see_px,
     float eta = p3LH.eta();
     float ptErr = see_ptErr[iSeed];
 
-    if ((ptIn > 0.8 - 2 * ptErr)) {
+    if ((ptIn > ptCut - 2 * ptErr)) {
       XYZVector r3LH(see_stateTrajGlbX[iSeed], see_stateTrajGlbY[iSeed], see_stateTrajGlbZ[iSeed]);
       XYZVector p3PCA(see_px[iSeed], see_py[iSeed], see_pz[iSeed]);
       XYZVector r3PCA(calculateR3FromPCA(p3PCA, see_dxy[iSeed], see_dz[iSeed]));
@@ -149,7 +150,7 @@ void LST::prepareInput(std::vector<float> const& see_px,
 
       if (ptIn >= 2.0)
         pixtype = PixelType::kHighPt;
-      else if (ptIn >= (0.8 - 2 * ptErr) and ptIn < 2.0) {
+      else if (ptIn >= (ptCut - 2 * ptErr) and ptIn < 2.0) {
         if (pixelSegmentDeltaPhiChange >= 0)
           pixtype = PixelType::kLowPtPosCurv;
         else
@@ -255,6 +256,7 @@ void LST::getOutput(LSTEvent& event) {
 
 void LST::run(Queue& queue,
               bool verbose,
+              float const ptCut,
               LSTESData<Device> const* deviceESData,
               std::vector<float> const& see_px,
               std::vector<float> const& see_py,
@@ -277,7 +279,7 @@ void LST::run(Queue& queue,
               std::vector<float> const& ph2_z,
               bool no_pls_dupclean,
               bool tc_pls_triplets) {
-  auto event = LSTEvent(verbose, queue, deviceESData);
+  auto event = LSTEvent(verbose, ptCut, queue, deviceESData);
   prepareInput(see_px,
                see_py,
                see_pz,
@@ -296,7 +298,8 @@ void LST::run(Queue& queue,
                ph2_detId,
                ph2_x,
                ph2_y,
-               ph2_z);
+               ph2_z,
+               ptCut);
 
   event.addHitToEvent(in_trkX_, in_trkY_, in_trkZ_, in_hitId_, in_hitIdxs_);
   event.addPixelSegmentToEvent(in_hitIndices_vec0_,