Make the CA, ntuplet fit, and track clustering modules configurable

RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGeneratorGPU.cc

@@ -34,7 +34,16 @@ CAHitQuadrupletGeneratorGPU::CAHitQuadrupletGeneratorGPU(const edm::ParameterSet
               cfg.getParameter<bool>("earlyFishbone"),
               cfg.getParameter<bool>("lateFishbone"),
               cfg.getParameter<bool>("idealConditions"),
-              cfg.getParameter<bool>("fillStatistics")),
+              cfg.getParameter<bool>("fillStatistics"),
+              cfg.getParameter<bool>("doClusterCut"),
+              cfg.getParameter<bool>("doZCut"),
+              cfg.getParameter<bool>("doPhiCut"),
+              cfg.getParameter<double>("ptmin"),
+              cfg.getParameter<double>("CAThetaCutBarrel"),
+              cfg.getParameter<double>("CAThetaCutForward"),
+              cfg.getParameter<double>("hardCurvCut"),
+              cfg.getParameter<double>("dcaCutInnerTriplet"),
+              cfg.getParameter<double>("dcaCutOuterTriplet")),
       fitter(cfg.getParameter<bool>("fit5as4")),
       caThetaCut(cfg.getParameter<double>("CAThetaCut")),
       caPhiCut(cfg.getParameter<double>("CAPhiCut")),
@@ -44,11 +53,23 @@ void CAHitQuadrupletGeneratorGPU::fillDescriptions(edm::ParameterSetDescription
   desc.add<double>("CAThetaCut", 0.00125);
   desc.add<double>("CAPhiCut", 10);
   desc.add<double>("CAHardPtCut", 0);
+  // 87 cm/GeV = 1/(3.8T * 0.3)
+  // take less than radius given by the hardPtCut and reject everything below
+  // auto hardCurvCut = 1.f/(0.35 * 87.f);
+  desc.add<double>("ptmin", 0.9f)->setComment("Cut on minimum pt");
+  desc.add<double>("CAThetaCutBarrel", 0.002f)->setComment("Cut on RZ alignement for Barrel");
+  desc.add<double>("CAThetaCutForward", 0.003f)->setComment("Cut on RZ alignment for Forward");
+  desc.add<double>("hardCurvCut", 1.f / (0.35 * 87.f))->setComment("Cut on minimum curvature");
+  desc.add<double>("dcaCutInnerTriplet", 0.15f)->setComment("Cut on origin radius when the inner hit is on BPix1");
+  desc.add<double>("dcaCutOuterTriplet", 0.25f)->setComment("Cut on origin radius when the outer hit is on BPix1");
   desc.add<bool>("earlyFishbone", false);
   desc.add<bool>("lateFishbone", true);
   desc.add<bool>("idealConditions", true), desc.add<bool>("fillStatistics", false),
       desc.add<unsigned int>("minHitsPerNtuplet", 4);
   desc.add<bool>("fit5as4", true);
+  desc.add<bool>("doClusterCut", true);
+  desc.add<bool>("doZCut", true);
+  desc.add<bool>("doPhiCut", true);
 }
 
 void CAHitQuadrupletGeneratorGPU::initEvent(edm::Event const &ev, edm::EventSetup const &es) {

RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGeneratorKernels.cu

@@ -163,15 +163,15 @@ __global__ void kernel_connect(AtomicPairCounter *apc1,
                                GPUCACell *cells,
                                uint32_t const *__restrict__ nCells,
                                CellNeighborsVector *cellNeighbors,
-                               GPUCACell::OuterHitOfCell const *__restrict__ isOuterHitOfCell) {
+                               GPUCACell::OuterHitOfCell const *__restrict__ isOuterHitOfCell,
+                               float hardCurvCut,
+                               float ptmin,
+                               float CAThetaCutBarrel,
+                               float CAThetaCutForward,
+                               float dcaCutInnerTriplet,
+                               float dcaCutOuterTriplet) {
   auto const &hh = *hhp;
 
-  // 87 cm/GeV = 1/(3.8T * 0.3)
-  // take less than radius given by the hardPtCut and reject everything below
-  // auto hardCurvCut = 1.f/(hardPtCut * 87.f);
-  constexpr auto hardCurvCut = 1.f / (0.35f * 87.f);  // FIXME VI tune
-  constexpr auto ptmin = 0.9f;                        // FIXME original "tune"
-
   auto cellIndex = threadIdx.y + blockIdx.y * blockDim.y;
   auto first = threadIdx.x;
   auto stride = blockDim.x;
@@ -193,7 +193,14 @@ __global__ void kernel_connect(AtomicPairCounter *apc1,
     auto otherCell = __ldg(vi + j);
     if (cells[otherCell].theDoubletId < 0)
       continue;
-    if (thisCell.check_alignment(hh, cells[otherCell], ptmin, hardCurvCut)) {
+    if (thisCell.check_alignment(hh,
+                                 cells[otherCell],
+                                 ptmin,
+                                 hardCurvCut,
+                                 CAThetaCutBarrel,
+                                 CAThetaCutForward,
+                                 dcaCutInnerTriplet,
+                                 dcaCutOuterTriplet)) {
       cells[otherCell].addOuterNeighbor(cellIndex, *cellNeighbors);
     }
   }
@@ -476,7 +483,13 @@ void CAHitQuadrupletGeneratorKernels::launchKernels(  // here goes algoparms....
       device_theCells_.get(),
       device_nCells_,
       device_theCellNeighbors_,
-      device_isOuterHitOfCell_.get());
+      device_isOuterHitOfCell_.get(),
+      hardCurvCut_,
+      ptmin_,
+      CAThetaCutBarrel_,
+      CAThetaCutForward_,
+      dcaCutInnerTriplet_,
+      dcaCutOuterTriplet_);
   cudaCheck(cudaGetLastError());
 
   kernel_find_ntuplets<<<numberOfBlocks, blockSize, 0, cudaStream>>>(hh.view(),
@@ -569,7 +582,10 @@ void CAHitQuadrupletGeneratorKernels::buildDoublets(HitsOnCPU const &hh, cuda::s
                                                                          device_theCellTracks_,
                                                                          hh.view(),
                                                                          device_isOuterHitOfCell_.get(),
-                                                                         idealConditions_);
+                                                                         idealConditions_,
+                                                                         doClusterCut_,
+                                                                         doZCut_,
+                                                                         doPhiCut_);
   cudaCheck(cudaGetLastError());
 }
 

RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGeneratorKernels.h

@@ -31,13 +31,34 @@ class CAHitQuadrupletGeneratorKernels {
   using HitToTuple = CAConstants::HitToTuple;
   using TupleMultiplicity = CAConstants::TupleMultiplicity;
 
-  CAHitQuadrupletGeneratorKernels(
-      uint32_t minHitsPerNtuplet, bool earlyFishbone, bool lateFishbone, bool idealConditions, bool doStats)
+  CAHitQuadrupletGeneratorKernels(uint32_t minHitsPerNtuplet,
+                                  bool earlyFishbone,
+                                  bool lateFishbone,
+                                  bool idealConditions,
+                                  bool doStats,
+                                  bool doClusterCut,
+                                  bool doZCut,
+                                  bool doPhiCut,
+                                  float ptmin,
+                                  float CAThetaCutBarrel,
+                                  float CAThetaCutForward,
+                                  float hardCurvCut,
+                                  float dcaCutInnerTriplet,
+                                  float dcaCutOuterTriplet)
       : minHitsPerNtuplet_(minHitsPerNtuplet),
         earlyFishbone_(earlyFishbone),
         lateFishbone_(lateFishbone),
         idealConditions_(idealConditions),
-        doStats_(doStats) {}
+        doStats_(doStats),
+        doClusterCut_(doClusterCut),
+        doZCut_(doZCut),
+        doPhiCut_(doPhiCut),
+        ptmin_(ptmin),
+        CAThetaCutBarrel_(CAThetaCutBarrel),
+        CAThetaCutForward_(CAThetaCutForward),
+        hardCurvCut_(hardCurvCut),
+        dcaCutInnerTriplet_(dcaCutInnerTriplet),
+        dcaCutOuterTriplet_(dcaCutOuterTriplet){};
   ~CAHitQuadrupletGeneratorKernels() { deallocateOnGPU(); }
 
   TupleMultiplicity const* tupleMultiplicity() const { return device_tupleMultiplicity_; }
@@ -77,6 +98,15 @@ class CAHitQuadrupletGeneratorKernels {
   const bool lateFishbone_;
   const bool idealConditions_;
   const bool doStats_;
+  const bool doClusterCut_;
+  const bool doZCut_;
+  const bool doPhiCut_;
+  const float ptmin_;
+  const float CAThetaCutBarrel_;
+  const float CAThetaCutForward_;
+  const float hardCurvCut_;
+  const float dcaCutInnerTriplet_;
+  const float dcaCutOuterTriplet_;
 };
 
 #endif  // RecoPixelVertexing_PixelTriplets_plugins_CAHitQuadrupletGeneratorKernels_h

RecoPixelVertexing/PixelTriplets/plugins/GPUCACell.h

@@ -86,8 +86,8 @@ class GPUCACell {
   __device__ __forceinline__ auto get_inner_iphi(Hits const& hh) const { return hh.iphi(theInnerHitId); }
   __device__ __forceinline__ auto get_outer_iphi(Hits const& hh) const { return hh.iphi(theOuterHitId); }
 
-  __device__ __forceinline__ float get_inner_detId(Hits const& hh) const { return hh.detectorIndex(theInnerHitId); }
-  __device__ __forceinline__ float get_outer_detId(Hits const& hh) const { return hh.detectorIndex(theOuterHitId); }
+  __device__ __forceinline__ float get_inner_detIndex(Hits const& hh) const { return hh.detectorIndex(theInnerHitId); }
+  __device__ __forceinline__ float get_outer_detIndex(Hits const& hh) const { return hh.detectorIndex(theOuterHitId); }
 
   constexpr unsigned int get_inner_hit_id() const { return theInnerHitId; }
   constexpr unsigned int get_outer_hit_id() const { return theOuterHitId; }
@@ -105,7 +105,16 @@ class GPUCACell {
   __device__ bool check_alignment(Hits const& hh,
                                   GPUCACell const& otherCell,
                                   const float ptmin,
-                                  const float hardCurvCut) const {
+                                  const float hardCurvCut,
+                                  const float CAThetaCutBarrel,
+                                  const float CAThetaCutForward,
+                                  const float dcaCutInnerTriplet,
+                                  const float dcaCutOuterTriplet) const {
+    // detIndex of the layerStart for the Phase1 Pixel Detector:
+    // [BPX1, BPX2, BPX3, BPX4,  FP1,  FP2,  FP3,  FN1,  FN2,  FN3, LAST_VALID]
+    // [   0,   96,  320,  672, 1184, 1296, 1408, 1520, 1632, 1744,       1856]
+    constexpr uint32_t last_bpix1_detIndex = 96;
+    constexpr uint32_t last_barrel_detIndex = 1184;
     auto ri = get_inner_r(hh);
     auto zi = get_inner_z(hh);
 
@@ -114,12 +123,21 @@ class GPUCACell {
 
     auto r1 = otherCell.get_inner_r(hh);
     auto z1 = otherCell.get_inner_z(hh);
-    auto isBarrel = otherCell.get_outer_detId(hh) < 1184;
-    bool aligned =
-        areAlignedRZ(r1, z1, ri, zi, ro, zo, ptmin, isBarrel ? 0.002f : 0.003f);  // 2.f*thetaCut); // FIXME tune cuts
+    auto isBarrel = otherCell.get_outer_detIndex(hh) < last_barrel_detIndex;
+    bool aligned = areAlignedRZ(r1,
+                                z1,
+                                ri,
+                                zi,
+                                ro,
+                                zo,
+                                ptmin,
+                                isBarrel ? CAThetaCutBarrel : CAThetaCutForward);  // 2.f*thetaCut); // FIXME tune cuts
     return (aligned &&
-            dcaCut(hh, otherCell, otherCell.get_inner_detId(hh) < 96 ? 0.15f : 0.25f, hardCurvCut));  // FIXME tune cuts
-    // region_origin_radius_plus_tolerance,  hardCurvCut));
+            dcaCut(hh,
+                   otherCell,
+                   otherCell.get_inner_detIndex(hh) < last_bpix1_detIndex ? dcaCutInnerTriplet : dcaCutOuterTriplet,
+                   hardCurvCut));  // FIXME tune cuts
+                                   // region_origin_radius_plus_tolerance,  hardCurvCut));
   }
 
   __device__ __forceinline__ static bool areAlignedRZ(
@@ -156,20 +174,27 @@ class GPUCACell {
   }
 
   __device__ inline bool hole(Hits const& hh, GPUCACell const& innerCell) const {
+    constexpr uint32_t max_ladder_bpx4 = 64;
+    constexpr float radius_even_ladder = 15.815f;
+    constexpr float radius_odd_ladder = 16.146f;
+    constexpr float ladder_length = 6.7f;
+    constexpr float ladder_tolerance = 0.2f;
+    constexpr float barrel_z_length = 26.f;
+    constexpr float forward_z_begin = 32.f;
     int p = get_outer_iphi(hh);
     if (p < 0)
       p += std::numeric_limits<unsigned short>::max();
-    p = (64 * p) / std::numeric_limits<unsigned short>::max();
+    p = (max_ladder_bpx4 * p) / std::numeric_limits<unsigned short>::max();
     p %= 2;
-    float r4 = p == 0 ? 15.815 : 16.146;  // later on from geom
+    float r4 = p == 0 ? radius_even_ladder : radius_odd_ladder;  // later on from geom
     auto ri = innerCell.get_inner_r(hh);
     auto zi = innerCell.get_inner_z(hh);
     auto ro = get_outer_r(hh);
     auto zo = get_outer_z(hh);
     auto z4 = std::abs(zi + (r4 - ri) * (zo - zi) / (ro - ri));
-    auto zm = z4 - 6.7 * int(z4 / 6.7);
-    auto h = zm < 0.2 || zm > 6.5;
-    return h || (z4 > 26 && z4 < 32.f);
+    auto z_in_ladder = z4 - ladder_length * int(z4 / ladder_length);
+    auto h = z_in_ladder < ladder_tolerance || z_in_ladder > (ladder_length - ladder_tolerance);
+    return h || (z4 > barrel_z_length && z4 < forward_z_begin);
   }
 
   // trying to free the track building process from hardcoded layers, leaving

RecoPixelVertexing/PixelTriplets/plugins/gpuFishbone.h

@@ -53,7 +53,7 @@ namespace gpuPixelDoublets {
       if (checkTrack && ci.tracks().empty())
         continue;
       cc[sg] = vc[ic];
-      d[sg] = ci.get_inner_detId(hh);
+      d[sg] = ci.get_inner_detIndex(hh);
       //      l[sg] = layer(d[sg]);
       x[sg] = ci.get_inner_x(hh) - xo;
       y[sg] = ci.get_inner_y(hh) - yo;

RecoPixelVertexing/PixelTriplets/plugins/gpuPixelDoublets.h

@@ -95,7 +95,10 @@ namespace gpuPixelDoublets {
                                                                     CellTracksVector* cellTracks,
                                                                     TrackingRecHit2DSOAView const* __restrict__ hhp,
                                                                     GPUCACell::OuterHitOfCell* isOuterHitOfCell,
-                                                                    bool ideal_cond) {
+                                                                    bool ideal_cond,
+                                                                    bool doClusterCut,
+                                                                    bool doZCut,
+                                                                    bool doPhiCut) {
     auto const& __restrict__ hh = *hhp;
     doubletsFromHisto(layerPairs,
                       nPairs,
@@ -109,7 +112,10 @@ namespace gpuPixelDoublets {
                       minz,
                       maxz,
                       maxr,
-                      ideal_cond);
+                      ideal_cond,
+                      doClusterCut,
+                      doZCut,
+                      doPhiCut);
   }
 
 }  // namespace gpuPixelDoublets

RecoPixelVertexing/PixelTriplets/plugins/gpuPixelDoubletsAlgos.h

@@ -15,11 +15,6 @@
 #include "CAConstants.h"
 #include "GPUCACell.h"
 
-// useful for benchmark
-// #define ONLY_PHICUT
-// #define NO_ZCUT
-// #define NO_CLSCUT
-
 namespace gpuPixelDoubletsAlgos {
 
   constexpr uint32_t MaxNumOfDoublets = CAConstants::maxNumberOfDoublets();  // not really relevant
@@ -43,14 +38,18 @@ namespace gpuPixelDoubletsAlgos {
                                                     float const* __restrict__ minz,
                                                     float const* __restrict__ maxz,
                                                     float const* __restrict__ maxr,
-                                                    bool ideal_cond) {
-#ifndef NO_CLSCUT
+                                                    bool ideal_cond,
+                                                    bool doClusterCut,
+                                                    bool doZCut,
+                                                    bool doPhiCut) {
     // ysize cuts (z in the barrel)  times 8
+    // these are used if doClusterCut is true
     constexpr int minYsizeB1 = 36;
     constexpr int minYsizeB2 = 28;
     constexpr int maxDYsize12 = 28;
     constexpr int maxDYsize = 20;
-#endif
+    int16_t mes;
+    bool isOuterLadder = ideal_cond;
 
     using Hist = TrackingRecHit2DSOAView::Hist;
 
@@ -107,29 +106,25 @@ namespace gpuPixelDoubletsAlgos {
       // found hit corresponding to our cuda thread, now do the job
       auto mez = hh.zGlobal(i);
 
-#ifndef NO_ZCUT
-      if (mez < minz[pairLayerId] || mez > maxz[pairLayerId])
+      if (doZCut && (mez < minz[pairLayerId] || mez > maxz[pairLayerId]))
         continue;
-#endif
 
-#ifndef NO_CLSCUT
-      auto mes = hh.clusterSizeY(i);
+      if (doClusterCut) {
+        auto mes = hh.clusterSizeY(i);
 
-      // if ideal treat inner ladder as outer
-      auto mi = hh.detectorIndex(i);
-      if (inner == 0)
-        assert(mi < 96);
-      const bool isOuterLadder =
-          ideal_cond ? true : 0 == (mi / 8) % 2;  // only for B1/B2/B3 B4 is opposite, FPIX:noclue...
-
-      if (inner == 0 && outer > 3 && isOuterLadder)  // B1 and F1
-        if (mes > 0 && mes < minYsizeB1)
-          continue;                 // only long cluster  (5*8)
-      if (inner == 1 && outer > 3)  // B2 and F1
-        if (mes > 0 && mes < minYsizeB2)
-          continue;
-#endif  // NO_CLSCUT
+        // if ideal treat inner ladder as outer
+        auto mi = hh.detectorIndex(i);
+        if (inner == 0)
+          assert(mi < 96);
+        isOuterLadder = ideal_cond ? true : 0 == (mi / 8) % 2;  // only for B1/B2/B3 B4 is opposite, FPIX:noclue...
 
+        if (inner == 0 && outer > 3 && isOuterLadder)  // B1 and F1
+          if (mes > 0 && mes < minYsizeB1)
+            continue;                 // only long cluster  (5*8)
+        if (inner == 1 && outer > 3)  // B2 and F1
+          if (mes > 0 && mes < minYsizeB2)
+            continue;
+      }
       auto mep = hh.iphi(i);
       auto mer = hh.rGlobal(i);
 
@@ -153,14 +148,12 @@ namespace gpuPixelDoubletsAlgos {
         return dr > maxr[pairLayerId] || dr < 0 || std::abs((mez * ro - mer * zo)) > z0cut * dr;
       };
 
-#ifndef NO_CLSCUT
       auto zsizeCut = [&](int j) {
         auto onlyBarrel = outer < 4;
         auto so = hh.clusterSizeY(j);
         auto dy = inner == 0 ? (isOuterLadder ? maxDYsize12 : 100) : maxDYsize;
         return onlyBarrel && mes > 0 && so > 0 && std::abs(so - mes) > dy;
       };
-#endif
 
       auto iphicut = phicuts[pairLayerId];
 
@@ -191,14 +184,12 @@ namespace gpuPixelDoubletsAlgos {
 
           if (std::min(std::abs(int16_t(hh.iphi(oi) - mep)), std::abs(int16_t(mep - hh.iphi(oi)))) > iphicut)
             continue;
-#ifndef ONLY_PHICUT
-#ifndef NO_CLSCUT
-          if (zsizeCut(oi))
-            continue;
-#endif
-          if (z0cutoff(oi) || ptcut(oi))
-            continue;
-#endif
+          if (doPhiCut) {
+            if (doClusterCut && zsizeCut(oi))
+              continue;
+            if (z0cutoff(oi) || ptcut(oi))
+              continue;
+          }
           auto ind = atomicAdd(nCells, 1);
           if (ind >= MaxNumOfDoublets) {
             atomicSub(nCells, 1);