Move LSTGeometry repo into ESProducer and standalone binary

RecoTracker/LST/plugins/alpaka/LSTModulesDevESProducer.cc

@@ -1,5 +1,6 @@
 // LST includes
 #include "RecoTracker/LSTCore/interface/alpaka/LST.h"
+#include "RecoTracker/LSTCore/interface/LSTGeometry/LSTGeometry.h"
 
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 
@@ -14,22 +15,25 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
   class LSTModulesDevESProducer : public ESProducer {
   private:
-    std::string ptCutLabel_;
+    double ptCut_;
+    edm::ESGetToken<lstgeometry::LSTGeometry, TrackerRecoGeometryRecord> lstGeoToken_;
 
   public:
     LSTModulesDevESProducer(edm::ParameterSet const& iConfig)
-        : ESProducer(iConfig), ptCutLabel_(iConfig.getParameter<std::string>("ptCutLabel")) {
-      setWhatProduced(this, ptCutLabel_);
+        : ESProducer(iConfig), ptCut_(iConfig.getParameter<double>("ptCut")) {
+      auto cc = setWhatProduced(this, "LSTModuleMaps");
+      lstGeoToken_ = cc.consumes<lstgeometry::LSTGeometry>(edm::ESInputTag("", "LSTGeometry"));
     }
 
     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
       edm::ParameterSetDescription desc;
-      desc.add<std::string>("ptCutLabel", "0.8");
+      desc.add<double>("ptCut", 0.8);
       descriptions.addWithDefaultLabel(desc);
     }
 
     std::unique_ptr<lst::LSTESData<DevHost>> produce(TrackerRecoGeometryRecord const& iRecord) {
-      return lst::loadAndFillESHost(ptCutLabel_);
+      const auto& lstg = iRecord.get(lstGeoToken_);
+      return lst::loadAndFillESHost(lstg);
     }
   };
 

RecoTracker/LST/plugins/alpaka/LSTProducer.cc

@@ -26,7 +26,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
     LSTProducer(edm::ParameterSet const& config)
         : EDProducer(config),
           lstInputToken_{consumes(config.getParameter<edm::InputTag>("lstInput"))},
-          lstESToken_{esConsumes(edm::ESInputTag("", config.getParameter<std::string>("ptCutLabel")))},
+          lstESToken_{esConsumes(edm::ESInputTag("", "LSTModuleMaps"))},
           verbose_(config.getParameter<bool>("verbose")),
           ptCut_(config.getParameter<double>("ptCut")),
           clustSizeCut_(static_cast<uint16_t>(config.getParameter<uint32_t>("clustSizeCut"))),
@@ -42,7 +42,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
 
       lst.run(iEvent.queue(),
               verbose_,
-              static_cast<float>(ptCut_),
+              ptCut_,
               clustSizeCut_,
               &lstESDeviceData,
               &lstInputDC,
@@ -60,7 +60,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE {
       desc.add<bool>("verbose", false);
       desc.add<double>("ptCut", 0.8);
       desc.add<uint32_t>("clustSizeCut", 16);
-      desc.add<std::string>("ptCutLabel", "0.8");
       desc.add<bool>("nopLSDupClean", false);
       desc.add<bool>("tcpLSTriplets", false);
       descriptions.addWithDefaultLabel(desc);

RecoTracker/LST/python/lstProducerTask_cff.py

@@ -4,4 +4,8 @@
 
 from RecoTracker.LST.lstModulesDevESProducer_cfi import lstModulesDevESProducer
 
-lstProducerTask = cms.Task(lstModulesDevESProducer, lstProducer)
+from RecoTracker.LST.lstInputProducer_cfi import lstInputProducer
+
+from RecoTracker.LSTCore.lstGeometryESProducer_cfi import lstGeometryESProducer
+
+lstProducerTask = cms.Task(lstGeometryESProducer, lstModulesDevESProducer, lstInputProducer, lstProducer)

RecoTracker/LSTCore/BuildFile.xml

@@ -5,6 +5,7 @@
 <use name="DataFormats/Common"/>
 <use name="DataFormats/TrackerRecHit2D"/>
 <use name="HeterogeneousCore/AlpakaInterface"/>
+<use name="eigen"/>
 <flags ALPAKA_BACKENDS="1"/>
 <export>
   <lib name="1"/>

RecoTracker/LSTCore/interface/EndcapGeometry.h

@@ -1,6 +1,8 @@
 #ifndef RecoTracker_LSTCore_interface_EndcapGeometry_h
 #define RecoTracker_LSTCore_interface_EndcapGeometry_h
 
+#include "RecoTracker/LSTCore/interface/LSTGeometry/LSTGeometry.h"
+
 #include <map>
 #include <string>
 #include <vector>
@@ -21,6 +23,8 @@ namespace lst {
     EndcapGeometry(std::string const& filename);
 
     void load(std::string const&);
+    void load(std::unordered_map<unsigned int, lstgeometry::SlopeData> const&,
+              std::unordered_map<unsigned int, lstgeometry::SensorInfo> const&);
     void fillGeoMapArraysExplicit();
     float getdxdy_slope(unsigned int detid) const;
   };

RecoTracker/LSTCore/interface/LSTESData.h

@@ -5,6 +5,7 @@
 #include "RecoTracker/LSTCore/interface/EndcapGeometryDevHostCollection.h"
 #include "RecoTracker/LSTCore/interface/ModulesHostCollection.h"
 #include "RecoTracker/LSTCore/interface/PixelMap.h"
+#include "RecoTracker/LSTCore/interface/LSTGeometry/LSTGeometry.h"
 
 #include "HeterogeneousCore/AlpakaInterface/interface/CopyToDevice.h"
 
@@ -41,6 +42,7 @@ namespace lst {
   };
 
   std::unique_ptr<LSTESData<alpaka_common::DevHost>> loadAndFillESHost(std::string& ptCutLabel);
+  std::unique_ptr<LSTESData<alpaka_common::DevHost>> loadAndFillESHost(lstgeometry::LSTGeometry const& lstg);
 
 }  // namespace lst
 

RecoTracker/LSTCore/interface/LSTGeometry/Centroid.h

@@ -0,0 +1,15 @@
+#ifndef RecoTracker_LSTCore_interface_LSTGeometry_Centroid_h
+#define RecoTracker_LSTCore_interface_LSTGeometry_Centroid_h
+
+namespace lstgeometry {
+
+  struct Centroid {
+    unsigned int moduleType;
+    double x;
+    double y;
+    double z;
+  };
+
+}  // namespace lstgeometry
+
+#endif

RecoTracker/LSTCore/interface/LSTGeometry/CentroidMethods.h

@@ -0,0 +1,80 @@
+#ifndef RecoTracker_LSTCore_interface_LSTGeometry_CentroidMethods_h
+#define RecoTracker_LSTCore_interface_LSTGeometry_CentroidMethods_h
+
+#include <stdexcept>
+#include <unordered_map>
+
+#include "Common.h"
+#include "RecoTracker/LSTCore/interface/LSTGeometry/Module.h"
+#include "RecoTracker/LSTCore/interface/LSTGeometry/Centroid.h"
+#include "RecoTracker/LSTCore/interface/LSTGeometry/SensorInfo.h"
+
+namespace lstgeometry {
+
+  unsigned int extractBits(unsigned int value, unsigned int start, unsigned int end) {
+    unsigned int mask = (1 << (end - start + 1)) - 1;
+    return (value >> start) & mask;
+  }
+
+  unsigned int firstDigit(unsigned int n) {
+    while (n >= 10) {
+      n /= 10;
+    }
+    return n;
+  }
+
+  // TODO: refactor to use Module class better
+  int parseModuleType(unsigned int detId) {
+    // Check if the first digit of detId is '3' for inner tracker
+    if (firstDigit(detId) == 3)
+      return -1;
+
+    unsigned int subdet = extractBits(detId, 25, 27);
+    unsigned int layer = subdet == Module::SubDet::Barrel ? extractBits(detId, 20, 22) : extractBits(detId, 18, 20);
+    unsigned int ring = subdet == Module::SubDet::Endcap ? extractBits(detId, 12, 15) : 0;
+
+    bool is_even_det_id = detId % 2 == 0;
+    if (subdet == Module::SubDet::Barrel) {
+      if (layer <= 3)
+        return is_even_det_id ? Module::ModuleType::PSS : Module::ModuleType::PSP;
+      else
+        return Module::ModuleType::TwoS;
+    } else if (subdet == Module::SubDet::Endcap) {
+      if (layer <= 2)
+        return is_even_det_id && ring <= 10 ? Module::ModuleType::PSS
+                                            : (ring <= 10 ? Module::ModuleType::PSP : Module::ModuleType::TwoS);
+      else
+        return is_even_det_id && ring <= 7 ? Module::ModuleType::PSS
+                                           : (ring <= 7 ? Module::ModuleType::PSP : Module::ModuleType::TwoS);
+    } else {
+      throw std::runtime_error("Invalid subdetector type");
+    }
+  }
+
+  std::unordered_map<unsigned int, Centroid> computeCentroids(
+      std::unordered_map<unsigned int, SensorInfo> const& sensors) {
+    std::unordered_map<unsigned int, Centroid> centroids;
+    for (auto const& [detId, sensor] : sensors) {
+      int moduleType = parseModuleType(detId);
+
+      // Remove sensors from inner tracker
+      if (moduleType == -1) {
+        continue;
+      }
+
+      // Convert from mm to cm
+      double z = sensor.sensorCenterZ_cm;
+      double rho = sensor.sensorCenterRho_cm;
+      double phi = sensor.phi_rad;
+      double x = rho * cos(phi);
+      double y = rho * sin(phi);
+
+      Centroid centroid{static_cast<unsigned int>(moduleType), x, y, z};
+      centroids[detId] = centroid;
+    }
+    return centroids;
+  }
+
+}  // namespace lstgeometry
+
+#endif

RecoTracker/LSTCore/interface/LSTGeometry/Common.h

@@ -0,0 +1,65 @@
+#ifndef RecoTracker_LSTCore_interface_LSTGeometry_Common_h
+#define RecoTracker_LSTCore_interface_LSTGeometry_Common_h
+
+#include <Eigen/Dense>
+#include <limits>
+
+namespace lstgeometry {
+
+  using MatrixD3x3 = Eigen::Matrix<double, 3, 3, Eigen::RowMajor>;
+  using MatrixD4x2 = Eigen::Matrix<double, 4, 2, Eigen::RowMajor>;
+  using MatrixD4x3 = Eigen::Matrix<double, 4, 3, Eigen::RowMajor>;
+  using MatrixD8x3 = Eigen::Matrix<double, 8, 3, Eigen::RowMajor>;
+  using MatrixDNx2 = Eigen::Matrix<double, Eigen::Dynamic, 2, Eigen::RowMajor>;
+  using MatrixDNx3 = Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor>;
+  using RowVectorD2 = Eigen::Matrix<double, 1, 2>;
+  using ColVectorD3 = Eigen::Matrix<double, 3, 1>;
+  using RowVectorD3 = Eigen::Matrix<double, 1, 3>;
+
+  // TODO: These should be moved to ../Common.h
+  constexpr double k2Rinv1GeVf = 0.00299792458;
+  constexpr double kB = 3.8112;
+
+  // For pixel maps
+  constexpr unsigned int kNEta = 25;
+  constexpr unsigned int kNPhi = 72;
+  constexpr unsigned int kNZ = 25;
+  constexpr std::array<double, 2> kPtBounds = {{2.0, 10'000.0}};
+
+  // This is defined as a constant in case the legacy value (123456789) needs to be used
+  double kDefaultSlope = std::numeric_limits<double>::infinity();
+
+  double degToRad(double degrees) { return degrees * (std::numbers::pi_v<double> / 180); }
+
+  double phi_mpi_pi(double phi) {
+    while (phi >= std::numbers::pi_v<double>)
+      phi -= 2 * std::numbers::pi_v<double>;
+    while (phi < -std::numbers::pi_v<double>)
+      phi += 2 * std::numbers::pi_v<double>;
+    return phi;
+  }
+
+  double roundAngle(double angle, double tol = 1e-3) {
+    const double pi = std::numbers::pi_v<double>;
+    if (std::fabs(angle) < tol) {
+      return 0.0;
+    } else if (std::fabs(angle - pi / 2) < tol) {
+      return pi / 2;
+    } else if (std::fabs(angle + pi / 2) < tol) {
+      return -pi / 2;
+    } else if (std::fabs(angle - pi) < tol || std::fabs(angle + pi) < tol) {
+      return -pi;
+    }
+    return angle;
+  }
+
+  double roundCoordinate(double coord, double tol = 1e-3) {
+    if (std::fabs(coord) < tol) {
+      return 0.0;
+    }
+    return coord;
+  }
+
+}  // namespace lstgeometry
+
+#endif