T4 Creation + DNN Draft PR

RecoTracker/LSTCore/interface/Constants.h

@@ -45,7 +45,7 @@ namespace lst {
   constexpr unsigned int n_max_pixel_quintuplets = 15000;
 
   constexpr unsigned int n_max_pixel_track_candidates = 30000;
-  constexpr unsigned int n_max_nonpixel_track_candidates = 1000;
+  constexpr unsigned int n_max_nonpixel_track_candidates = 11000;
 
   constexpr unsigned int size_superbins = 45000;
 
@@ -60,6 +60,9 @@ namespace lst {
   struct Params_T3 {
     static constexpr int kLayers = 3, kHits = 6;
   };
+  struct Params_T4 {
+    static constexpr int kLayers = 4, kHits = 8;
+  };
   struct Params_pT3 {
     static constexpr int kLayers = 5, kHits = 10;
   };

RecoTracker/LSTCore/interface/alpaka/Constants.h

@@ -60,7 +60,7 @@ namespace lst {
   // The constants below are usually used in functions like alpaka::math::min(),
   // expecting a reference (T const&) in the arguments. Hence,
   // ALPAKA_STATIC_ACC_MEM_GLOBAL needs to be used in addition to constexpr.
-
+  ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPi = float(M_PI);
   // 15 MeV constant from the approximate Bethe-Bloch formula
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMulsInGeV = 0.015;
   ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMiniMulsPtScaleBarrel[6] = {
@@ -104,6 +104,18 @@ namespace lst {
         {0.268, 0.3172, 0.3907, 0.4653, 0.4364, 0.4696, 0.6018, 0.6487, 0.7401, 0.7146},
         {0.1654, 0.2385, 0.2935, 0.3534, 0.2455, 0.1748, 0.1565, 0.1811, 0.3435, 0.1784}};
   }  // namespace t5dnn
+  namespace t4dnn {
+    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.0623, 0.0897, 0.1367, 0.2101, 0.332, 0.4381, 0.6449, 0.7911, 0.8863, 0.9467}, //98% retention
+        {0.0041, 0.0096, 0.0006, 0.0053, 0.0179, 0.0382, 0.0852, 0.1545, 0.5131, 0.3933}};
+  }  // namespace t4dnn
 
 }  //namespace lst
 #endif

RecoTracker/LSTCore/src/alpaka/Event.h

@@ -14,6 +14,7 @@
 #include "PixelQuintuplet.h"
 #include "PixelTriplet.h"
 #include "TrackCandidate.h"
+#include "Quadruplet.h"
 
 #include "HeterogeneousCore/AlpakaInterface/interface/host.h"
 
@@ -45,6 +46,8 @@ namespace lst {
     std::array<unsigned int, 5> n_trackCandidates_by_layer_endcap_;
     std::array<unsigned int, 6> n_quintuplets_by_layer_barrel_;
     std::array<unsigned int, 5> n_quintuplets_by_layer_endcap_;
+    std::array<unsigned int, 6> n_quadruplets_by_layer_barrel_;
+    std::array<unsigned int, 5> n_quadruplets_by_layer_endcap_;
 
     //Device stuff
     unsigned int nTotalSegments;
@@ -66,6 +69,8 @@ namespace lst {
     PixelTripletsBuffer<Device>* pixelTripletsBuffers;
     PixelQuintuplets* pixelQuintupletsInGPU;
     PixelQuintupletsBuffer<Device>* pixelQuintupletsBuffers;
+    Quadruplets* quadrupletsInGPU;
+    QuadrupletsBuffer<Device>* quadrupletsBuffers;
 
     //CPU interface stuff
     ObjectRangesBuffer<DevHost>* rangesInCPU;
@@ -78,6 +83,7 @@ namespace lst {
     QuintupletsBuffer<DevHost>* quintupletsInCPU;
     PixelTripletsBuffer<DevHost>* pixelTripletsInCPU;
     PixelQuintupletsBuffer<DevHost>* pixelQuintupletsInCPU;
+    QuadrupletsBuffer<DevHost>* quadrupletsInCPU;
 
     void init(bool verbose);
 
@@ -145,6 +151,7 @@ namespace lst {
     void addSegmentsToEventExplicit();
     void addTripletsToEventExplicit();
     void addQuintupletsToEventExplicit();
+    void addQuadrupletsToEventExplicit();
     void resetObjectsInModule();
 
     void createMiniDoublets();
@@ -158,6 +165,7 @@ namespace lst {
     void createPixelTriplets();
     void createPixelQuintuplets();
     void pixelLineSegmentCleaning(bool no_pls_dupclean);
+    void createQuadruplets();
 
     unsigned int getNumberOfHits();
     unsigned int getNumberOfHitsByLayer(unsigned int layer);
@@ -185,6 +193,7 @@ namespace lst {
     int getNumberOfPT3TrackCandidates();
     int getNumberOfT5TrackCandidates();
     int getNumberOfPLSTrackCandidates();
+    int getNumberOfT4TrackCandidates();
 
     unsigned int getNumberOfQuintuplets();
     unsigned int getNumberOfQuintupletsByLayer(unsigned int layer);
@@ -194,6 +203,11 @@ namespace lst {
     int getNumberOfPixelTriplets();
     int getNumberOfPixelQuintuplets();
 
+    unsigned int getNumberOfQuadruplets();
+    unsigned int getNumberOfQuadrupletsByLayer(unsigned int layer);
+    unsigned int getNumberOfQuadrupletsByLayerBarrel(unsigned int layer);
+    unsigned int getNumberOfQuadrupletsByLayerEndcap(unsigned int layer);
+
     ObjectRangesBuffer<DevHost>* getRanges();
     HitsBuffer<DevHost>* getHits();
     HitsBuffer<DevHost>* getHitsInCMSSW();
@@ -206,6 +220,7 @@ namespace lst {
     PixelTripletsBuffer<DevHost>* getPixelTriplets();
     PixelQuintupletsBuffer<DevHost>* getPixelQuintuplets();
     ModulesBuffer<DevHost>* getModules(bool isFull = false);
+    QuadrupletsBuffer<DevHost>* getQuadruplets();
   };
 
 }  // namespace lst

RecoTracker/LSTCore/src/alpaka/Kernels.h

@@ -12,6 +12,7 @@
 #include "Quintuplet.h"
 #include "PixelQuintuplet.h"
 #include "PixelTriplet.h"
+#include "Quadruplet.h"
 
 namespace lst {
   ALPAKA_FN_ACC ALPAKA_FN_INLINE void rmQuintupletFromMemory(lst::Quintuplets& quintupletsInGPU,
@@ -36,6 +37,11 @@ namespace lst {
     segmentsInGPU.isDup[pixelSegmentArrayIndex] |= 1 + secondpass;
   };
 
+  ALPAKA_FN_ACC ALPAKA_FN_INLINE void rmQuadrupletFromMemory(lst::Quadruplets& quadrupletsInGPU,
+                                                             unsigned int quadrupletIndex) {
+    quadrupletsInGPU.isDup[quadrupletIndex] = true;
+  };
+
   ALPAKA_FN_ACC ALPAKA_FN_INLINE int checkHitsT5(unsigned int ix,
                                                  unsigned int jx,
                                                  lst::Quintuplets const& quintupletsInGPU) {
@@ -142,6 +148,34 @@ namespace lst {
     matched[1] = nMatched;
   };
 
+  ALPAKA_FN_ACC ALPAKA_FN_INLINE int checkHitsT4(unsigned int ix,
+                                                 unsigned int jx,
+                                                 lst::Quadruplets const& quadrupletsInGPU) {
+    unsigned int hits1[Params_T4::kHits];
+    unsigned int hits2[Params_T4::kHits];
+
+    for (int i = 0; i < Params_T4::kHits; i++) {
+      hits1[i] = quadrupletsInGPU.hitIndices[Params_T4::kHits * ix + i];
+      hits2[i] = quadrupletsInGPU.hitIndices[Params_T4::kHits * jx + i];
+    }
+
+    int nMatched = 0;
+    for (int i = 0; i < Params_T4::kHits; i++) {
+      bool matched = false;
+      for (int j = 0; j < Params_T4::kHits; j++) {
+        if (hits1[i] == hits2[j]) {
+          matched = true;
+          break;
+        }
+      }
+      if (matched) {
+        nMatched++;
+      }
+    }
+    return nMatched;
+  };
+
+
   struct removeDupQuintupletsInGPUAfterBuild {
     template <typename TAcc>
     ALPAKA_FN_ACC void operator()(TAcc const& acc,
@@ -274,6 +308,107 @@ namespace lst {
     }
   };
 
+  struct removeDupQuadrupletsInGPUAfterBuild {
+    template <typename TAcc>
+    ALPAKA_FN_ACC void operator()(TAcc const& acc,
+                                  lst::Modules modulesInGPU,
+                                  lst::Quadruplets quadrupletsInGPU,
+                                  lst::ObjectRanges rangesInGPU) const {
+      auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
+      auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
+
+      for (unsigned int lowmod = globalThreadIdx[0]; lowmod < *modulesInGPU.nLowerModules;
+            lowmod += gridThreadExtent[0]) {
+        unsigned int nQuadruplets_lowmod = quadrupletsInGPU.nQuadruplets[lowmod];
+        int quadrupletModuleIndices_lowmod = rangesInGPU.quadrupletModuleIndices[lowmod];
+
+        for (unsigned int ix1 = globalThreadIdx[1]; ix1 < nQuadruplets_lowmod; ix1 += gridThreadExtent[1]) {
+          unsigned int ix = quadrupletModuleIndices_lowmod + ix1;
+          float eta1 = __H2F(quadrupletsInGPU.eta[ix]);
+          float phi1 = __H2F(quadrupletsInGPU.phi[ix]);
+
+          for (unsigned int jx1 = globalThreadIdx[2] + ix1 + 1; jx1 < nQuadruplets_lowmod; jx1 += gridThreadExtent[2]) {
+            unsigned int jx = quadrupletModuleIndices_lowmod + jx1;
+
+            float eta2 = __H2F(quadrupletsInGPU.eta[jx]);
+            float phi2 = __H2F(quadrupletsInGPU.phi[jx]);
+            float dEta = alpaka::math::abs(acc, eta1 - eta2);
+            float dPhi = lst::calculate_dPhi(phi1, phi2);
+
+            if (dEta > 0.1f)
+              continue;
+
+            if (alpaka::math::abs(acc, dPhi) > 0.1f)
+              continue;
+
+            int nMatched = checkHitsT4(ix, jx, quadrupletsInGPU);
+            const int minNHitsForDup_T4 = 6;
+            if (nMatched >= minNHitsForDup_T4) {
+              rmQuadrupletFromMemory(quadrupletsInGPU, (ix > jx ? ix : jx));
+            }
+          }
+        }
+      }
+      // add in second pass of dup removal
+      for (unsigned int lowmodIdx1 = globalThreadIdx[1]; lowmodIdx1 < *(rangesInGPU.nEligibleT4Modules);
+           lowmodIdx1 += gridThreadExtent[1]) {
+        uint16_t lowmod1 = rangesInGPU.indicesOfEligibleT4Modules[lowmodIdx1];
+        unsigned int nQuadruplets_lowmod1 = quadrupletsInGPU.nQuadruplets[lowmod1];
+        if (nQuadruplets_lowmod1 == 0)
+          continue;
+
+        unsigned int quadrupletModuleIndices_lowmod1 = rangesInGPU.quadrupletModuleIndices[lowmod1];
+
+        for (unsigned int lowmodIdx2 = globalThreadIdx[2] + lowmodIdx1; lowmodIdx2 < *(rangesInGPU.nEligibleT4Modules);
+             lowmodIdx2 += gridThreadExtent[2]) {
+          uint16_t lowmod2 = rangesInGPU.indicesOfEligibleT4Modules[lowmodIdx2];
+          unsigned int nQuadruplets_lowmod2 = quadrupletsInGPU.nQuadruplets[lowmod2];
+          if (nQuadruplets_lowmod2 == 0)
+            continue;
+
+          unsigned int quadrupletModuleIndices_lowmod2 = rangesInGPU.quadrupletModuleIndices[lowmod2];
+
+          for (unsigned int ix1 = 0; ix1 < nQuadruplets_lowmod1; ix1 += 1) {
+            unsigned int ix = quadrupletModuleIndices_lowmod1 + ix1;
+            if ((quadrupletsInGPU.isDup[ix] & 1))
+              continue;
+
+            for (unsigned int jx1 = 0; jx1 < nQuadruplets_lowmod2; jx1++) {
+              unsigned int jx = quadrupletModuleIndices_lowmod2 + jx1;
+              if (ix == jx)
+                continue;
+
+              if ((quadrupletsInGPU.isDup[jx] & 1))
+                continue;
+
+              float eta1 = __H2F(quadrupletsInGPU.eta[ix]);
+              float phi1 = __H2F(quadrupletsInGPU.phi[ix]);
+
+              float eta2 = __H2F(quadrupletsInGPU.eta[jx]);
+              float phi2 = __H2F(quadrupletsInGPU.phi[jx]);
+
+              float dEta = alpaka::math::abs(acc, eta1 - eta2);
+              float dPhi = lst::calculate_dPhi(phi1, phi2);
+
+              if (dEta > 0.1f)
+                continue;
+
+              if (alpaka::math::abs(acc, dPhi) > 0.1f)
+                continue;
+
+              float dR2 = dEta * dEta + dPhi * dPhi;
+              int nMatched = checkHitsT4(ix, jx, quadrupletsInGPU);
+              const int minNHitsForDup_T4 = 6;
+              if (dR2 < 0.001f || nMatched >= minNHitsForDup_T4) {
+                  rmQuadrupletFromMemory(quadrupletsInGPU, (ix < jx ? ix : jx)); // just remove one randomly for now
+              }
+            }
+          }
+        }
+      }
+    }
+  };
+
   struct removeDupPixelTripletsInGPUFromMap {
     template <typename TAcc>
     ALPAKA_FN_ACC void operator()(TAcc const& acc, lst::PixelTriplets pixelTripletsInGPU) const {

RecoTracker/LSTCore/src/alpaka/LST.dev.cc

@@ -418,6 +418,22 @@ void lst::LST<Acc3D>::run(Queue& queue,
   event.createPixelTriplets();
   if (verbose)
     printf("# of Pixel T3s produced: %d\n", event.getNumberOfPixelTriplets());
+
+  event.createQuadruplets();
+  if (verbose) {
+    printf("# of Quadruplets produced: %d\n", event.getNumberOfQuadruplets());
+    printf("# of Quadruplets produced layer 1-2-3-4: %d\n", event.getNumberOfQuadrupletsByLayerBarrel(0));
+    printf("# of Quadruplets produced layer 2: %d\n", event.getNumberOfQuadrupletsByLayerBarrel(1));
+    printf("# of Quadruplets produced layer 3: %d\n", event.getNumberOfQuadrupletsByLayerBarrel(2));
+    printf("# of Quadruplets produced layer 4: %d\n", event.getNumberOfQuadrupletsByLayerBarrel(3));
+    printf("# of Quadruplets produced layer 5: %d\n", event.getNumberOfQuadrupletsByLayerBarrel(4));
+    printf("# of Quadruplets produced layer 6: %d\n", event.getNumberOfQuadrupletsByLayerBarrel(5));
+    printf("# of Quadruplets produced endcap layer 1: %d\n", event.getNumberOfQuadrupletsByLayerEndcap(0));
+    printf("# of Quadruplets produced endcap layer 2: %d\n", event.getNumberOfQuadrupletsByLayerEndcap(1));
+    printf("# of Quadruplets produced endcap layer 3: %d\n", event.getNumberOfQuadrupletsByLayerEndcap(2));
+    printf("# of Quadruplets produced endcap layer 4: %d\n", event.getNumberOfQuadrupletsByLayerEndcap(3));
+    printf("# of Quadruplets produced endcap layer 5: %d\n", event.getNumberOfQuadrupletsByLayerEndcap(4));
+  }
 
   event.createTrackCandidates(no_pls_dupclean, tc_pls_triplets);
   if (verbose) {
@@ -427,6 +443,7 @@ void lst::LST<Acc3D>::run(Queue& queue,
     printf("        # of pT3 TrackCandidates produced: %d\n", event.getNumberOfPT3TrackCandidates());
     printf("        # of pLS TrackCandidates produced: %d\n", event.getNumberOfPLSTrackCandidates());
     printf("        # of T5 TrackCandidates produced: %d\n", event.getNumberOfT5TrackCandidates());
+    printf("        # of T4 TrackCandidates produced: %d\n", event.getNumberOfT4TrackCandidates());
   }
 
   getOutput(event);