updates to handle pLS with OT hits

HLTrigger/Configuration/python/HLT_75e33/modules/hltInitialStepSeeds_cfi.py

@@ -14,14 +14,13 @@
     TTRHBuilder = cms.string('WithTrackAngle'),
     usePV = cms.bool(False),
     includeFourthHit = cms.bool(True),
-    removeOTRechits = cms.bool(True),
+    removeOTRechits = cms.bool(False),
     produceComplement = cms.bool(False)
 )
 
 from Configuration.ProcessModifiers.hltPhase2LegacyTracking_cff import hltPhase2LegacyTracking
 hltPhase2LegacyTracking.toModify(hltInitialStepSeeds,
-    includeFourthHit = False,
-    removeOTRechits = False
+    includeFourthHit = False
 )
 
 from Configuration.ProcessModifiers.phase2_hlt_vertexTrimming_cff import phase2_hlt_vertexTrimming

HLTrigger/Configuration/python/HLT_75e33/modules/hltInitialStepTrackCandidates_cfi.py

@@ -64,6 +64,8 @@
     lstPixelSeeds = cms.InputTag('hltInputLST'),
     includeT5s = cms.bool(True),
     includeNonpLSTSs = cms.bool(False),
+    dropOTHitsPurePLS = cms.bool(True),
+    maxITHitsToDropOTHitsPurePLS = cms.int32(3),
     produceSeeds = cms.bool(False),
     produceTrackCandidates = cms.bool(True),
     propagatorAlong = cms.ESInputTag('', 'PropagatorWithMaterial'),

HLTrigger/Configuration/python/HLT_75e33/modules/hltInitialStepTrajectorySeedsLST_cfi.py

@@ -6,6 +6,8 @@
     lstPixelSeeds = cms.InputTag('hltInputLST'),
     includeT5s = cms.bool(True),
     includeNonpLSTSs = cms.bool(True),
+    dropOTHitsPurePLS = cms.bool(True),
+    maxITHitsToDropOTHitsPurePLS = cms.int32(3),
     produceSeeds = cms.bool(True),
     produceTrackCandidates = cms.bool(False),
     propagatorAlong = cms.ESInputTag('', 'PropagatorWithMaterial'),

RecoTracker/LST/plugins/LSTOutputConverter.cc

@@ -25,6 +25,8 @@
 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
 
+#include <ranges>
+
 class LSTOutputConverter : public edm::stream::EDProducer<> {
 public:
   explicit LSTOutputConverter(edm::ParameterSet const& iConfig);
@@ -40,6 +42,8 @@ class LSTOutputConverter : public edm::stream::EDProducer<> {
   const edm::EDGetTokenT<TrajectorySeedCollection> lstPixelSeedToken_;
   const bool includeT5s_;
   const bool includeNonpLSTSs_;
+  const bool dropOTHitsPurePLS_;
+  const int maxITHitsToDropOTHitsPurePLS_;
   const bool produceSeeds_;
   const bool produceTrackCandidates_;
   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken_;
@@ -68,6 +72,8 @@ LSTOutputConverter::LSTOutputConverter(edm::ParameterSet const& iConfig)
       lstPixelSeedToken_{consumes(iConfig.getParameter<edm::InputTag>("lstPixelSeeds"))},
       includeT5s_(iConfig.getParameter<bool>("includeT5s")),
       includeNonpLSTSs_(iConfig.getParameter<bool>("includeNonpLSTSs")),
+      dropOTHitsPurePLS_(iConfig.getParameter<bool>("dropOTHitsPurePLS")),
+      maxITHitsToDropOTHitsPurePLS_(iConfig.getParameter<int>("maxITHitsToDropOTHitsPurePLS")),
       produceSeeds_(iConfig.getParameter<bool>("produceSeeds")),
       produceTrackCandidates_(iConfig.getParameter<bool>("produceTrackCandidates")),
       mfToken_(esConsumes()),
@@ -110,6 +116,8 @@ void LSTOutputConverter::fillDescriptions(edm::ConfigurationDescriptions& descri
   desc.add<edm::InputTag>("lstPixelSeeds", edm::InputTag("lstInputProducer"));
   desc.add<bool>("includeT5s", true);
   desc.add<bool>("includeNonpLSTSs", false);
+  desc.add<bool>("dropOTHitsPurePLS", false);
+  desc.add<int>("maxITHitsToDropOTHitsPurePLS", 3);
   desc.add<bool>("produceSeeds", true);
   desc.add<bool>("produceTrackCandidates", true);
   desc.add("propagatorAlong", edm::ESInputTag{"", "PropagatorWithMaterial"});
@@ -167,16 +175,21 @@ void LSTOutputConverter::produce(edm::Event& iEvent, const edm::EventSetup& iSet
   for (unsigned int i = 0; i < nTrackCandidates; i++) {
     auto iType = lstOutput_view.trackCandidateType()[i];
     bool const isT5orT4 = (iType == lst::LSTObjType::T5 || iType == lst::LSTObjType::T4);
-    LogDebug("LSTOutputConverter") << " cand " << i << " " << iType << " " << lstOutput_view.pixelSeedIndex()[i];
+    const auto iSeed = lstOutput_view.pixelSeedIndex()[i];
+    const bool dropHitsOTpL =
+        iType == lst::LSTObjType::pLS && dropOTHitsPurePLS_ &&
+        std::prev(pixelSeeds[iSeed].recHits().end())->geographicalId().subdetId() > PixelSubdetector::PixelEndcap &&
+        std::ranges::count_if(pixelSeeds[iSeed].recHits(), [](const auto& h) {
+          return h.geographicalId().subdetId() <= PixelSubdetector::PixelEndcap;
+        }) <= maxITHitsToDropOTHitsPurePLS_;
+    LogDebug("LSTOutputConverter") << " cand " << i << " " << iType << " " << iSeed;
     TrajectorySeed seed;
     edm::RefToBase<TrajectorySeed> seedRef;
-    if (!isT5orT4) {
-      seed = pixelSeeds[lstOutput_view.pixelSeedIndex()[i]];
-      seedRef = {pixelSeedsRBP, lstOutput_view.pixelSeedIndex()[i]};
-    }
-
     edm::OwnVector<TrackingRecHit> recHits;
-    if (!isT5orT4) {
+    if (!isT5orT4 && !dropHitsOTpL) {
+      seed = pixelSeeds[iSeed];
+      seedRef = {pixelSeedsRBP, iSeed};
+
       for (auto const& hit : seed.recHits())
         recHits.push_back(hit.clone());
     }
@@ -189,7 +202,14 @@ void LSTOutputConverter::produce(edm::Event& iEvent, const edm::EventSetup& iSet
           unsigned int hitIdx = lstOutput_view.hitIndices()[i][layerSlot][hitSlot];
           if (hitIdx == lst::kTCEmptyHitIdx)
             continue;
-          recHits.push_back(OTHits[hitIdx]->clone());
+          bool hitOK = true;
+          for (auto const& hit : recHits)
+            if (hit.sharesInput(OTHits[hitIdx], TrackingRecHit::all)) {
+              hitOK = false;
+              break;
+            }
+          if (hitOK)
+            recHits.push_back(OTHits[hitIdx]->clone());
         }
       }
 
@@ -212,9 +232,7 @@ void LSTOutputConverter::produce(edm::Event& iEvent, const edm::EventSetup& iSet
           }
         }
       });
-    }
 
-    if (iType != lst::LSTObjType::pLS) {
       // For T5/T4: makeSeed is needed whenever seeds or TCs are produced, since the resulting
       // seed is the only source of initial state for T5/T4 track candidates.
       // For other pT objects: makeSeed is only needed for seed output.
@@ -264,11 +282,27 @@ void LSTOutputConverter::produce(edm::Event& iEvent, const edm::EventSetup& iSet
         LogDebug("LSTOutputConverter") << "Created a seed with " << trajectorySeed.nHits() << " " << ss.detId() << " "
                                        << ss.pt() << " " << ss.parameters().vector() << " " << ss.error(0);
       }
-    } else {
-      if (produceSeeds_) {
-        outputTS.emplace_back(seed);
-        outputpLSTS.emplace_back(seed);
+    } else if (produceSeeds_) {  // if (iType == lst::LSTObjType::pLS)
+      if (dropHitsOTpL) {        // true if need to drop OT hits
+        using Hit = SeedingHitSet::ConstRecHitPointer;
+        std::vector<Hit> hitsForSeed;
+        hitsForSeed.reserve(std::ranges::size(pixelSeeds[iSeed].recHits()));
+        for (auto const& hit : pixelSeeds[iSeed].recHits()) {
+          if (hit.geographicalId().subdetId() > PixelSubdetector::PixelEndcap)
+            continue;
+          hitsForSeed.emplace_back(dynamic_cast<Hit>(&hit));
+          recHits.push_back(hit.clone());
+        }
+        GlobalTrackingRegion region;
+        seedCreator_->init(region, iSetup, nullptr);
+        seedCreator_->makeSeed(seeds, hitsForSeed);
+        if (seeds.empty())
+          edm::LogInfo("LSTOutputConverter") << "failed to convert a pLS object to a seed" << i << iSeed;
+        seed = seeds[0];
+        seedRef = edm::RefToBase<TrajectorySeed>(edm::Ref(outputTSRP, outputTS.size()));
       }
+      outputTS.emplace_back(seed);
+      outputpLSTS.emplace_back(seed);
     }
 
     if (!produceTrackCandidates_)

RecoTracker/LSTCore/interface/Common.h

@@ -24,6 +24,9 @@ namespace lst {
 
   constexpr unsigned int kPixelModuleId = 1;
 
+  constexpr unsigned int kMaxPLSHitsInHitsSoA = 4;
+  constexpr unsigned int kMaxPLSHitBitsInHitsSoA = 8;
+
   constexpr unsigned int max_blocks = 80;
   constexpr unsigned int max_connected_modules = 40;
 

RecoTracker/LSTCore/interface/LSTInputSoA.h

@@ -26,7 +26,9 @@ namespace lst {
                       SOA_SCALAR(unsigned int, nHitsOT))
 
   GENERATE_SOA_LAYOUT(PixelSeedsSoALayout,
-                      SOA_COLUMN(Params_pLS::ArrayUxHits, hitIndices),
+                      SOA_COLUMN(unsigned int, firstHit),
+                      SOA_COLUMN(uint8_t, nHits),
+                      SOA_COLUMN(uint8_t, hitDetBits),  // IT 0, OT 1 up to min(nHits, 8)
                       SOA_COLUMN(float, deltaPhi),
                       SOA_COLUMN(unsigned int, seedIdx),
                       SOA_COLUMN(int, charge),

RecoTracker/LSTCore/interface/LSTPrepareInput.h

@@ -55,7 +55,9 @@ namespace lst {
     std::vector<uint16_t> hitClustSize;
     std::vector<unsigned int> hitId;
     std::vector<unsigned int> hitIdxs;
-    std::vector<Params_pLS::ArrayUxHits> hitIndices_vec;
+    std::vector<unsigned int> firstHit_vec;
+    std::vector<uint8_t> nHits_vec;
+    std::vector<uint8_t> hitDetBits_vec;
     std::vector<float> deltaPhi_vec;
     std::vector<float> ptIn_vec;
     std::vector<float> ptErr_vec;
@@ -78,7 +80,9 @@ namespace lst {
     px_vec.reserve(n_see);
     py_vec.reserve(n_see);
     pz_vec.reserve(n_see);
-    hitIndices_vec.reserve(n_see);
+    firstHit_vec.reserve(n_see);
+    nHits_vec.reserve(n_see);
+    hitDetBits_vec.reserve(n_see);
     ptIn_vec.reserve(n_see);
     ptErr_vec.reserve(n_see);
     etaErr_vec.reserve(n_see);
@@ -96,6 +100,7 @@ namespace lst {
     hitIdxs.resize(hit_size);
 
     std::iota(hitIdxs.begin(), hitIdxs.end(), 0);
+    int nHitsOT = ph2_x.size();
 
     for (size_t iSeed = 0; iSeed < n_see; iSeed++) {
       // Only needed for standalone
@@ -133,20 +138,31 @@ namespace lst {
         } else
           continue;
 
-        unsigned int hitIdx0 = hit_size + count;
-        count++;
-        unsigned int hitIdx1 = hit_size + count;
-        count++;
-        unsigned int hitIdx2 = hit_size + count;
-        count++;
-        unsigned int hitIdx3;
-        if (see_hitIdx[iSeed].size() <= 3)
-          hitIdx3 = hitIdx2;
-        else {
-          hitIdx3 = hit_size + count;
-          count++;
+        firstHit_vec.push_back(hit_size + count);
+        nHits_vec.push_back(see_hitIdx[iSeed].size());
+
+        const unsigned int nHitsToSoA =
+            std::min(static_cast<unsigned int>(see_hitIdx[iSeed].size()), kMaxPLSHitsInHitsSoA);
+        assert(nHitsToSoA >= 3);  // need to rework the logic for less than 3 hits
+        count += nHitsToSoA;
+
+        auto const& hTypes = see_hitType[iSeed];
+        auto constexpr intPixel = static_cast<int>(HitType::Pixel);
+        auto const& hIdxs = see_hitIdx[iSeed];
+        for (unsigned int iSH = 0; iSH < nHitsToSoA; iSH++) {
+          auto iH = iSH + 1 == nHitsToSoA ? see_hitIdx[iSeed].size() - 1 : iSH;  // include the last
+          hitId.push_back(hTypes[iH] == intPixel ? kPixelModuleId : ph2_detId[hIdxs[iH]]);
+          hitClustSize.push_back(hTypes[iH] == intPixel ? 1 : ph2_clustSize[hIdxs[iH]]);
+        }
+        uint8_t hitDetBits = 0;
+        uint8_t nToBits = std::min(kMaxPLSHitBitsInHitsSoA, static_cast<unsigned int>(see_hitIdx[iSeed].size()));
+        for (int iSH = 0; iSH < nToBits; iSH++) {
+          auto iH = iSH + 1 == nToBits ? see_hitIdx[iSeed].size() - 1 : iSH;  // include the last
+          hitDetBits |= (hTypes[iH] != intPixel) << iSH;
         }
+        hitDetBits_vec.push_back(hitDetBits);
 
+        // eventually these trk[XYZ] should be moved to the PixelSeeds SoA
         trkX.push_back(r3PCA.x());
         trkY.push_back(r3PCA.y());
         trkZ.push_back(r3PCA.z());
@@ -158,25 +174,17 @@ namespace lst {
         trkX.push_back(r3LH.x());
         trkY.push_back(r3LH.y());
         trkZ.push_back(r3LH.z());
-        auto const& hTypes = see_hitType[iSeed];
-        auto constexpr intPixel = static_cast<int>(HitType::Pixel);
-        auto const& hIdxs = see_hitIdx[iSeed];
-        for (int iSH = 0; iSH < 3; iSH++) {
-          hitId.push_back(hTypes[iSH] == intPixel ? kPixelModuleId : ph2_detId[hIdxs[iSH]]);
-          hitClustSize.push_back(hTypes[iSH] == intPixel ? 1 : ph2_clustSize[hIdxs[iSH]]);
-        }
-        if (see_hitIdx[iSeed].size() > 3) {
+        for (unsigned int iH = 3; iH < nHitsToSoA; iH++) {
           trkX.push_back(r3LH.x());
           trkY.push_back(see_dxy[iSeed]);
           trkZ.push_back(see_dz[iSeed]);
-          hitId.push_back(hTypes[3] == intPixel ? kPixelModuleId : ph2_detId[hIdxs[3]]);
-          hitClustSize.push_back(hTypes[3] == intPixel ? 1 : ph2_clustSize[hIdxs[3]]);
         }
+        assert(trkX.size() == count);
+
         px_vec.push_back(px);
         py_vec.push_back(py);
         pz_vec.push_back(pz);
 
-        hitIndices_vec.push_back({{hitIdx0, hitIdx1, hitIdx2, hitIdx3}});
         ptIn_vec.push_back(ptIn);
         ptErr_vec.push_back(ptErr);
         etaErr_vec.push_back(etaErr);
@@ -193,7 +201,7 @@ namespace lst {
         char isQuad = false;
         if (see_hitIdx[iSeed].size() > 3) {
           isQuad = true;
-          hitIdxs.push_back(see_hitIdx[iSeed][3]);
+          hitIdxs.push_back(see_hitIdx[iSeed].back());
         }
         float neta = 25.;
         float nphi = 72.;
@@ -209,7 +217,6 @@ namespace lst {
     }
 
     // Build the SoAs
-    int nHitsOT = ph2_x.size();
     int nHitsIT = trkX.size();
     int nPixelSeeds = ptIn_vec.size();
     if (static_cast<unsigned int>(nPixelSeeds) > n_max_pixel_segments_per_module) {
@@ -241,7 +248,9 @@ namespace lst {
     std::copy_n(hitIdxs.data(), nHitsIT + nHitsOT, hits.idxs().data());
 
     auto pixelSeeds = lstInputHC.view().pixelSeeds();
-    std::copy_n(hitIndices_vec.data(), nPixelSeeds, pixelSeeds.hitIndices().data());
+    std::copy_n(firstHit_vec.data(), nPixelSeeds, pixelSeeds.firstHit().data());
+    std::copy_n(nHits_vec.data(), nPixelSeeds, pixelSeeds.nHits().data());
+    std::copy_n(hitDetBits_vec.data(), nPixelSeeds, pixelSeeds.hitDetBits().data());
     std::copy_n(deltaPhi_vec.data(), nPixelSeeds, pixelSeeds.deltaPhi().data());
     std::copy_n(ptIn_vec.data(), nPixelSeeds, pixelSeeds.ptIn().data());
     std::copy_n(ptErr_vec.data(), nPixelSeeds, pixelSeeds.ptErr().data());

RecoTracker/LSTCore/interface/PixelSegmentsSoA.h

@@ -9,7 +9,7 @@
 namespace lst {
 
   GENERATE_SOA_LAYOUT(PixelSegmentsSoALayout,
-                      SOA_COLUMN(Params_pLS::ArrayUxHits, pLSHitsIdxs),
+                      SOA_COLUMN(Params_pLS::ArrayUxHits, pLSHitsIdxs),  // packed OT/IT<<32 + index
                       SOA_COLUMN(Params_pLS::ArrayFxEmbed, plsEmbed),
                       SOA_COLUMN(char, isDup),
                       SOA_COLUMN(bool, partOfPT5),

RecoTracker/LSTCore/interface/alpaka/Common.h

@@ -47,6 +47,7 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
   // To be updated with std::numeric_limits<float>::infinity() in the code and data files
   HOST_DEVICE_CONSTANT float kVerticalModuleSlope = 123456789.0;
   HOST_DEVICE_CONSTANT int kLogicalOTLayers = 11;  // logical OT layers are 1..11
+  HOST_DEVICE_CONSTANT auto kMaxPLSHitBitsInHitsSoA = ::lst::kMaxPLSHitBitsInHitsSoA;
 
   HOST_DEVICE_CONSTANT float kMiniDeltaTilted[3] = {0.26f, 0.26f, 0.26f};
   HOST_DEVICE_CONSTANT float kMiniDeltaFlat[6] = {0.26f, 0.16f, 0.16f, 0.18f, 0.18f, 0.18f};

RecoTracker/LSTCore/src/alpaka/Hit.h

@@ -16,6 +16,11 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
     return cms::alpakatools::deltaPhi(acc, x1, y1, x2 - x1, y2 - y1);
   }
 
+  ALPAKA_FN_ACC ALPAKA_FN_INLINE unsigned int packedHitIdx(unsigned int ih, HitsBaseConst hitsBase) {
+    constexpr int kOTBit = 1 << 31;
+    return hitsBase.idxs()[ih] | (hitsBase.detid()[ih] == kPixelModuleId ? 0 : kOTBit);
+  }
+
   struct ModuleRangesKernel {
     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
                                   ModulesConst modules,

RecoTracker/LSTCore/src/alpaka/NeuralNetwork.h

@@ -155,7 +155,6 @@ namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
       constexpr unsigned int kInputFeatures = 7;
       constexpr unsigned int kHiddenFeatures = 32;
       constexpr unsigned int kOutputFeatures = 1;
-
       float x[kInputFeatures] = {
           alpaka::math::log10(acc, rPhiChiSquared),
           alpaka::math::log10(acc, tripletRadius),