Patatrack integration - Pixel track reconstruction (10/N)

CUDADataFormats/Track/BuildFile.xml

@@ -0,0 +1,9 @@
+<use name="cuda"/>
+<use name="rootcore"/>
+<use name="CUDADataFormats/Common"/>
+<use name="DataFormats/Common"/>
+<use name="HeterogeneousCore/CUDAUtilities"/>
+<use name="eigen"/>
+<export>
+    <lib name="1"/>
+</export>

CUDADataFormats/Track/interface/PixelTrackHeterogeneous.h

@@ -0,0 +1,9 @@
+#ifndef CUDADataFormats_Track_PixelTrackHeterogeneous_h
+#define CUDADataFormats_Track_PixelTrackHeterogeneous_h
+
+#include "CUDADataFormats/Common/interface/HeterogeneousSoA.h"
+#include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h"
+
+using PixelTrackHeterogeneous = HeterogeneousSoA<pixelTrack::TrackSoA>;
+
+#endif  // #ifndef CUDADataFormats_Track_PixelTrackHeterogeneous_h

CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h

@@ -0,0 +1,73 @@
+#ifndef CUDADataFormats_Track_TrackHeterogeneousT_H
+#define CUDADataFormats_Track_TrackHeterogeneousT_H
+
+#include "CUDADataFormats/Track/interface/TrajectoryStateSoAT.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/HistoContainer.h"
+
+#include "CUDADataFormats/Common/interface/HeterogeneousSoA.h"
+
+namespace pixelTrack {
+  enum class Quality : uint8_t { bad = 0, dup, loose, strict, tight, highPurity };
+}
+
+template <int32_t S>
+class TrackSoAHeterogeneousT {
+public:
+  static constexpr int32_t stride() { return S; }
+
+  using Quality = pixelTrack::Quality;
+  using hindex_type = uint32_t;
+  using HitContainer = cms::cuda::OneToManyAssoc<hindex_type, S, 5 * S>;
+
+  // Always check quality is at least loose!
+  // CUDA does not support enums  in __lgc ...
+private:
+  eigenSoA::ScalarSoA<uint8_t, S> quality_;
+
+public:
+  constexpr Quality quality(int32_t i) const { return (Quality)(quality_(i)); }
+  constexpr Quality &quality(int32_t i) { return (Quality &)(quality_(i)); }
+  constexpr Quality const *qualityData() const { return (Quality const *)(quality_.data()); }
+  constexpr Quality *qualityData() { return (Quality *)(quality_.data()); }
+
+  // this is chi2/ndof as not necessarely all hits are used in the fit
+  eigenSoA::ScalarSoA<float, S> chi2;
+
+  constexpr int nHits(int i) const { return detIndices.size(i); }
+
+  // State at the Beam spot
+  // phi,tip,1/pt,cotan(theta),zip
+  TrajectoryStateSoAT<S> stateAtBS;
+  eigenSoA::ScalarSoA<float, S> eta;
+  eigenSoA::ScalarSoA<float, S> pt;
+  constexpr float charge(int32_t i) const { return std::copysign(1.f, stateAtBS.state(i)(2)); }
+  constexpr float phi(int32_t i) const { return stateAtBS.state(i)(0); }
+  constexpr float tip(int32_t i) const { return stateAtBS.state(i)(1); }
+  constexpr float zip(int32_t i) const { return stateAtBS.state(i)(4); }
+
+  // state at the detector of the outermost hit
+  // representation to be decided...
+  // not yet filled on GPU
+  // TrajectoryStateSoA<S> stateAtOuterDet;
+
+  HitContainer hitIndices;
+  HitContainer detIndices;
+};
+
+namespace pixelTrack {
+
+#ifdef GPU_SMALL_EVENTS
+  // kept for testing and debugging
+  constexpr uint32_t maxNumber() { return 2 * 1024; }
+#else
+  // tested on MC events with 55-75 pileup events
+  constexpr uint32_t maxNumber() { return 32 * 1024; }
+#endif
+
+  using TrackSoA = TrackSoAHeterogeneousT<maxNumber()>;
+  using TrajectoryState = TrajectoryStateSoAT<maxNumber()>;
+  using HitContainer = TrackSoA::HitContainer;
+
+}  // namespace pixelTrack
+
+#endif  // CUDADataFormats_Track_TrackHeterogeneousT_H

CUDADataFormats/Track/interface/TrajectoryStateSoAT.h

@@ -0,0 +1,59 @@
+#ifndef CUDADataFormats_Track_TrajectoryStateSOAT_H
+#define CUDADataFormats_Track_TrajectoryStateSOAT_H
+
+#include <Eigen/Dense>
+#include "HeterogeneousCore/CUDAUtilities/interface/eigenSoA.h"
+
+template <int32_t S>
+struct TrajectoryStateSoAT {
+  using Vector5f = Eigen::Matrix<float, 5, 1>;
+  using Vector15f = Eigen::Matrix<float, 15, 1>;
+
+  using Vector5d = Eigen::Matrix<double, 5, 1>;
+  using Matrix5d = Eigen::Matrix<double, 5, 5>;
+
+  static constexpr int32_t stride() { return S; }
+
+  eigenSoA::MatrixSoA<Vector5f, S> state;
+  eigenSoA::MatrixSoA<Vector15f, S> covariance;
+
+  template <typename V3, typename M3, typename V2, typename M2>
+  __host__ __device__ inline void copyFromCircle(
+      V3 const& cp, M3 const& ccov, V2 const& lp, M2 const& lcov, float b, int32_t i) {
+    state(i) << cp.template cast<float>(), lp.template cast<float>();
+    state(i)(2) *= b;
+    auto cov = covariance(i);
+    cov(0) = ccov(0, 0);
+    cov(1) = ccov(0, 1);
+    cov(2) = b * float(ccov(0, 2));
+    cov(4) = cov(3) = 0;
+    cov(5) = ccov(1, 1);
+    cov(6) = b * float(ccov(1, 2));
+    cov(8) = cov(7) = 0;
+    cov(9) = b * b * float(ccov(2, 2));
+    cov(11) = cov(10) = 0;
+    cov(12) = lcov(0, 0);
+    cov(13) = lcov(0, 1);
+    cov(14) = lcov(1, 1);
+  }
+
+  template <typename V5, typename M5>
+  __host__ __device__ inline void copyFromDense(V5 const& v, M5 const& cov, int32_t i) {
+    state(i) = v.template cast<float>();
+    for (int j = 0, ind = 0; j < 5; ++j)
+      for (auto k = j; k < 5; ++k)
+        covariance(i)(ind++) = cov(j, k);
+  }
+
+  template <typename V5, typename M5>
+  __host__ __device__ inline void copyToDense(V5& v, M5& cov, int32_t i) const {
+    v = state(i).template cast<typename V5::Scalar>();
+    for (int j = 0, ind = 0; j < 5; ++j) {
+      cov(j, j) = covariance(i)(ind++);
+      for (auto k = j + 1; k < 5; ++k)
+        cov(k, j) = cov(j, k) = covariance(i)(ind++);
+    }
+  }
+};
+
+#endif  // CUDADataFormats_Track_TrajectoryStateSOAT_H

CUDADataFormats/Track/src/classes.h

@@ -0,0 +1,9 @@
+#ifndef CUDADataFormats_Track_src_classes_h
+#define CUDADataFormats_Track_src_classes_h
+
+#include "CUDADataFormats/Common/interface/Product.h"
+#include "CUDADataFormats/Common/interface/HostProduct.h"
+#include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h"
+#include "DataFormats/Common/interface/Wrapper.h"
+
+#endif  // CUDADataFormats_Track_src_classes_h

CUDADataFormats/Track/src/classes_def.xml

@@ -0,0 +1,6 @@
+<lcgdict>
+  <class name="cms::cuda::Product<HeterogeneousSoA<pixelTrack::TrackSoA>>" persistent="false"/>
+  <class name="edm::Wrapper<cms::cuda::Product<HeterogeneousSoA<pixelTrack::TrackSoA>>>" persistent="false"/>
+  <class name="HeterogeneousSoA<pixelTrack::TrackSoA>" persistent="false"/>
+  <class name="edm::Wrapper<HeterogeneousSoA<pixelTrack::TrackSoA>>" persistent="false"/>
+</lcgdict>

CUDADataFormats/Track/test/BuildFile.xml

@@ -0,0 +1,13 @@
+<use name="HeterogeneousCore/CUDAUtilities"/>
+
+<bin file="TrajectoryStateSOA_t.cpp" name="cpuTrajectoryStateSOA_t">
+  <use name="eigen"/>
+  <flags CXXFLAGS="-g -DGPU_DEBUG"/>
+</bin>
+
+<bin file="TrajectoryStateSOA_t.cu" name="gpuTrajectoryStateSOA_t">
+  <use name="eigen"/>
+  <flags CUDA_FLAGS="-g -DGPU_DEBUG"/>
+  <flags CXXFLAGS="-g -DGPU_DEBUG"/>
+</bin>
+

CUDADataFormats/Track/test/TrajectoryStateSOA_t.cpp

@@ -0,0 +1 @@
+#include "TrajectoryStateSOA_t.h"

CUDADataFormats/Track/test/TrajectoryStateSOA_t.cu

@@ -0,0 +1 @@
+#include "TrajectoryStateSOA_t.h"

CUDADataFormats/Track/test/TrajectoryStateSOA_t.h

@@ -0,0 +1,75 @@
+#include "CUDADataFormats/Track/interface/TrajectoryStateSoAT.h"
+
+using Vector5d = Eigen::Matrix<double, 5, 1>;
+using Matrix5d = Eigen::Matrix<double, 5, 5>;
+
+__host__ __device__ Matrix5d loadCov(Vector5d const& e) {
+  Matrix5d cov;
+  for (int i = 0; i < 5; ++i)
+    cov(i, i) = e(i) * e(i);
+  for (int i = 0; i < 5; ++i) {
+    for (int j = 0; j < i; ++j) {
+      double v = 0.3 * std::sqrt(cov(i, i) * cov(j, j));  // this makes the matrix pos defined
+      cov(i, j) = (i + j) % 2 ? -0.4 * v : 0.1 * v;
+      cov(j, i) = cov(i, j);
+    }
+  }
+  return cov;
+}
+
+using TS = TrajectoryStateSoAT<128>;
+
+__global__ void testTSSoA(TS* pts, int n) {
+  assert(n <= 128);
+
+  Vector5d par0;
+  par0 << 0.2, 0.1, 3.5, 0.8, 0.1;
+  Vector5d e0;
+  e0 << 0.01, 0.01, 0.035, -0.03, -0.01;
+  auto cov0 = loadCov(e0);
+
+  TS& ts = *pts;
+
+  int first = threadIdx.x + blockIdx.x * blockDim.x;
+
+  for (int i = first; i < n; i += blockDim.x * gridDim.x) {
+    ts.copyFromDense(par0, cov0, i);
+    Vector5d par1;
+    Matrix5d cov1;
+    ts.copyToDense(par1, cov1, i);
+    Vector5d delV = par1 - par0;
+    Matrix5d delM = cov1 - cov0;
+    for (int j = 0; j < 5; ++j) {
+      assert(std::abs(delV(j)) < 1.e-5);
+      for (auto k = j; k < 5; ++k) {
+        assert(cov0(k, j) == cov0(j, k));
+        assert(cov1(k, j) == cov1(j, k));
+        assert(std::abs(delM(k, j)) < 1.e-5);
+      }
+    }
+  }
+}
+
+#ifdef __CUDACC__
+#include "HeterogeneousCore/CUDAUtilities/interface/requireDevices.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
+#endif
+
+int main() {
+#ifdef __CUDACC__
+  cms::cudatest::requireDevices();
+#endif
+
+  TS ts;
+
+#ifdef __CUDACC__
+  TS* ts_d;
+  cudaCheck(cudaMalloc(&ts_d, sizeof(TS)));
+  testTSSoA<<<1, 64>>>(ts_d, 128);
+  cudaCheck(cudaGetLastError());
+  cudaCheck(cudaMemcpy(&ts, ts_d, sizeof(TS), cudaMemcpyDefault));
+  cudaCheck(cudaDeviceSynchronize());
+#else
+  testTSSoA(&ts, 128);
+#endif
+}

DQM/TrackingMonitorClient/python/pixelTrackingEffFromHitPattern_cff.py

@@ -21,7 +21,10 @@ def _layers(suffix, quant, histoPostfix):
     ]
 
 pixelTrackingEffFromHitPattern = DQMEDHarvester("DQMGenericClient",
-    subDirs = cms.untracked.vstring("Tracking/PixelTrackParameters/HitEffFromHitPattern*"),
+    subDirs = cms.untracked.vstring("Tracking/PixelTrackParameters/pixelTracks/HitEffFromHitPattern*",
+                                    "Tracking/PixelTrackParameters/dzPV0p1/HitEffFromHitPattern*",
+                                    "Tracking/PixelTrackParameters/pt_0to1/HitEffFromHitPattern*",
+                                    "Tracking/PixelTrackParameters/pt_1/HitEffFromHitPattern*"),
     efficiency = cms.vstring(
         _layers("PU", "GoodNumVertices", "") +
         _layers("BX", "BX", "VsBX") +

DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py

@@ -1,23 +1,77 @@
 import FWCore.ParameterSet.Config as cms
 
 import DQM.TrackingMonitor.TrackerCollisionTrackingMonitor_cfi
-pixelTracksMonitoring = DQM.TrackingMonitor.TrackerCollisionTrackingMonitor_cfi.TrackerCollisionTrackMon.clone()
-pixelTracksMonitoring.FolderName                = 'Tracking/PixelTrackParameters'
-pixelTracksMonitoring.TrackProducer             = 'pixelTracks'
-pixelTracksMonitoring.allTrackProducer          = 'pixelTracks'
-pixelTracksMonitoring.beamSpot                  = 'offlineBeamSpot'
-pixelTracksMonitoring.primaryVertex             = 'pixelVertices'
-pixelTracksMonitoring.pvNDOF                    = 1
-pixelTracksMonitoring.doAllPlots                = True
-pixelTracksMonitoring.doLumiAnalysis            = True
-pixelTracksMonitoring.doProfilesVsLS            = True
-pixelTracksMonitoring.doDCAPlots                = True
-pixelTracksMonitoring.doProfilesVsLS            = True
-pixelTracksMonitoring.doPlotsVsGoodPVtx         = True
-pixelTracksMonitoring.doEffFromHitPatternVsPU   = False
-pixelTracksMonitoring.doEffFromHitPatternVsBX   = False
-pixelTracksMonitoring.doEffFromHitPatternVsLUMI = False
-pixelTracksMonitoring.doPlotsVsGoodPVtx         = True
-pixelTracksMonitoring.doPlotsVsLUMI             = True
-pixelTracksMonitoring.doPlotsVsBX               = True
+pixelTracksMonitor = DQM.TrackingMonitor.TrackerCollisionTrackingMonitor_cfi.TrackerCollisionTrackMon.clone()
+pixelTracksMonitor.FolderName                = 'Tracking/PixelTrackParameters/pixelTracks'
+pixelTracksMonitor.TrackProducer             = 'pixelTracks'
+pixelTracksMonitor.allTrackProducer          = 'pixelTracks'
+pixelTracksMonitor.beamSpot                  = 'offlineBeamSpot'
+pixelTracksMonitor.primaryVertex             = 'pixelVertices'
+pixelTracksMonitor.pvNDOF                    = 1
+pixelTracksMonitor.doAllPlots                = True
+pixelTracksMonitor.doLumiAnalysis            = True
+pixelTracksMonitor.doProfilesVsLS            = True
+pixelTracksMonitor.doDCAPlots                = True
+pixelTracksMonitor.doProfilesVsLS            = True
+pixelTracksMonitor.doPlotsVsGoodPVtx         = True
+pixelTracksMonitor.doEffFromHitPatternVsPU   = False
+pixelTracksMonitor.doEffFromHitPatternVsBX   = False
+pixelTracksMonitor.doEffFromHitPatternVsLUMI = False
+pixelTracksMonitor.doPlotsVsGoodPVtx         = True
+pixelTracksMonitor.doPlotsVsLUMI             = True
+pixelTracksMonitor.doPlotsVsBX               = True
 
+_trackSelector = cms.EDFilter('TrackSelector',
+    src = cms.InputTag('pixelTracks'),
+    cut = cms.string("")
+)
+
+pixelTracksPt0to1 = _trackSelector.clone(cut = "pt >= 0 & pt < 1 ")
+pixelTracksPt1 = _trackSelector.clone(cut = "pt >= 1 ")
+from DQM.TrackingMonitorSource.TrackCollections2monitor_cff import highPurityPV0p1 as _highPurityPV0p1
+pixelTracksPV0p1 = _highPurityPV0p1.clone(
+    src = "pixelTracks",
+    quality = "",
+    vertexTag = "goodPixelVertices"
+)
+
+pixelTracksMonitorPt0to1 = pixelTracksMonitor.clone(
+    TrackProducer = "pixelTracksPt0to1",
+    FolderName = "Tracking/PixelTrackParameters/pt_0to1"
+)
+pixelTracksMonitorPt1 = pixelTracksMonitor.clone(
+    TrackProducer = "pixelTracksPt1",
+    FolderName = "Tracking/PixelTrackParameters/pt_1"
+)
+pixelTracksMonitorPV0p1 = pixelTracksMonitor.clone(
+    TrackProducer = "pixelTracksPV0p1",
+    FolderName = "Tracking/PixelTrackParameters/dzPV0p1"
+)
+
+
+from CommonTools.ParticleFlow.goodOfflinePrimaryVertices_cfi import goodOfflinePrimaryVertices as _goodOfflinePrimaryVertices
+goodPixelVertices = _goodOfflinePrimaryVertices.clone(
+    src = "pixelVertices",
+)
+
+from DQM.TrackingMonitor.primaryVertexResolution_cfi import primaryVertexResolution as _primaryVertexResolution
+pixelVertexResolution = _primaryVertexResolution.clone(
+    vertexSrc = "goodPixelVertices",
+    rootFolder = "OfflinePixelPV/Resolution",
+)
+
+pixelTracksMonitoringTask = cms.Task(
+    goodPixelVertices,
+    pixelTracksPt0to1,
+    pixelTracksPt1,
+    pixelTracksPV0p1,
+)
+
+pixelTracksMonitoring = cms.Sequence(
+    pixelTracksMonitor +
+    pixelTracksMonitorPt0to1 +
+    pixelTracksMonitorPt1 +
+    pixelTracksMonitorPV0p1 +
+    pixelVertexResolution,
+    pixelTracksMonitoringTask
+)