Track/interface/PixelTrackUtilities.h

0001 #ifndef CUDADataFormats_Track_PixelTrackUtilities_h
0002 #define CUDADataFormats_Track_PixelTrackUtilities_h
0003
0004 #include <Eigen/Dense>
0005 #include <Eigen/Core>
0006 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0007 #include "HeterogeneousCore/CUDAUtilities/interface/HistoContainer.h"
0008 #include "DataFormats/SoATemplate/interface/SoALayout.h"
0009
0010 namespace pixelTrack {
0011
0012   enum class Quality : uint8_t { bad = 0, edup, dup, loose, strict, tight, highPurity, notQuality };
0013   constexpr uint32_t qualitySize{uint8_t(Quality::notQuality)};
0014   const std::string qualityName[qualitySize]{"bad", "edup", "dup", "loose", "strict", "tight", "highPurity"};
0015   inline Quality qualityByName(std::string const &name) {
0016     auto qp = std::find(qualityName, qualityName + qualitySize, name) - qualityName;
0017     return static_cast<Quality>(qp);
0018   }
0019
0020 }  // namespace pixelTrack
0021
0022 template <typename TrackerTraits>
0023 struct TrackSoA {
0024   static constexpr int32_t S = TrackerTraits::maxNumberOfTuples;
0025   static constexpr int32_t H = TrackerTraits::avgHitsPerTrack;
0026   // Aliases in order to not confuse the GENERATE_SOA_LAYOUT
0027   // macro with weird colons and angled brackets.
0028   using Vector5f = Eigen::Matrix<float, 5, 1>;
0029   using Vector15f = Eigen::Matrix<float, 15, 1>;
0030   using Quality = pixelTrack::Quality;
0031
0032   using hindex_type = uint32_t;
0033
0034   using HitContainer = cms::cuda::OneToManyAssoc<hindex_type, S + 1, H * S>;
0035
0036   GENERATE_SOA_LAYOUT(TrackSoALayout,
0037                       SOA_COLUMN(Quality, quality),
0038                       SOA_COLUMN(float, chi2),
0039                       SOA_COLUMN(int8_t, nLayers),
0040                       SOA_COLUMN(float, eta),
0041                       SOA_COLUMN(float, pt),
0042                       SOA_EIGEN_COLUMN(Vector5f, state),
0043                       SOA_EIGEN_COLUMN(Vector15f, covariance),
0044                       SOA_SCALAR(int, nTracks),
0045                       SOA_SCALAR(HitContainer, hitIndices),
0046                       SOA_SCALAR(HitContainer, detIndices))
0047 };
0048
0049 // Methods that operate on View and ConstView of the TrackSoA, and cannot be class methods.
0050
0051 template <typename TrackerTraits>
0052 struct TracksUtilities {
0053   using TrackSoAView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::View;
0054   using TrackSoAConstView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::ConstView;
0055   using hindex_type = typename TrackSoA<TrackerTraits>::hindex_type;
0056
0057   // State at the Beam spot
0058   // phi,tip,1/pt,cotan(theta),zip
0059   static __host__ __device__ inline float charge(const TrackSoAConstView &tracks, int32_t i) {
0060     return std::copysign(1.f, tracks[i].state()(2));
0061   }
0062
0063   static constexpr __host__ __device__ inline float phi(const TrackSoAConstView &tracks, int32_t i) {
0064     return tracks[i].state()(0);
0065   }
0066
0067   static constexpr __host__ __device__ inline float tip(const TrackSoAConstView &tracks, int32_t i) {
0068     return tracks[i].state()(1);
0069   }
0070
0071   static constexpr __host__ __device__ inline float zip(const TrackSoAConstView &tracks, int32_t i) {
0072     return tracks[i].state()(4);
0073   }
0074
0075   static constexpr __host__ __device__ inline bool isTriplet(const TrackSoAConstView &tracks, int i) {
0076     return tracks[i].nLayers() == 3;
0077   }
0078
0079   template <typename V3, typename M3, typename V2, typename M2>
0080   static constexpr __host__ __device__ inline void copyFromCircle(
0081       TrackSoAView &tracks, V3 const &cp, M3 const &ccov, V2 const &lp, M2 const &lcov, float b, int32_t i) {
0082     tracks[i].state() << cp.template cast<float>(), lp.template cast<float>();
0083
0084     tracks[i].state()(2) = tracks[i].state()(2) * b;
0085     auto cov = tracks[i].covariance();
0086     cov(0) = ccov(0, 0);
0087     cov(1) = ccov(0, 1);
0088     cov(2) = b * float(ccov(0, 2));
0089     cov(4) = cov(3) = 0;
0090     cov(5) = ccov(1, 1);
0091     cov(6) = b * float(ccov(1, 2));
0092     cov(8) = cov(7) = 0;
0093     cov(9) = b * b * float(ccov(2, 2));
0094     cov(11) = cov(10) = 0;
0095     cov(12) = lcov(0, 0);
0096     cov(13) = lcov(0, 1);
0097     cov(14) = lcov(1, 1);
0098   }
0099
0100   template <typename V5, typename M5>
0101   static constexpr __host__ __device__ inline void copyFromDense(TrackSoAView &tracks,
0102                                                                  V5 const &v,
0103                                                                  M5 const &cov,
0104                                                                  int32_t i) {
0105     tracks[i].state() = v.template cast<float>();
0106     for (int j = 0, ind = 0; j < 5; ++j)
0107       for (auto k = j; k < 5; ++k)
0108         tracks[i].covariance()(ind++) = cov(j, k);
0109   }
0110
0111   template <typename V5, typename M5>
0112   static constexpr __host__ __device__ inline void copyToDense(const TrackSoAConstView &tracks,
0113                                                                V5 &v,
0114                                                                M5 &cov,
0115                                                                int32_t i) {
0116     v = tracks[i].state().template cast<typename V5::Scalar>();
0117     for (int j = 0, ind = 0; j < 5; ++j) {
0118       cov(j, j) = tracks[i].covariance()(ind++);
0119       for (auto k = j + 1; k < 5; ++k)
0120         cov(k, j) = cov(j, k) = tracks[i].covariance()(ind++);
0121     }
0122   }
0123
0124   static constexpr __host__ __device__ inline int computeNumberOfLayers(const TrackSoAConstView &tracks, int32_t i) {
0125     auto pdet = tracks.detIndices().begin(i);
0126     int nl = 1;
0127     auto ol = pixelTopology::getLayer<TrackerTraits>(*pdet);
0128     for (; pdet < tracks.detIndices().end(i); ++pdet) {
0129       auto il = pixelTopology::getLayer<TrackerTraits>(*pdet);
0130       if (il != ol)
0131         ++nl;
0132       ol = il;
0133     }
0134     return nl;
0135   }
0136
0137   static constexpr __host__ __device__ inline int nHits(const TrackSoAConstView &tracks, int i) {
0138     return tracks.detIndices().size(i);
0139   }
0140 };
0141
0142 namespace pixelTrack {
0143
0144   template <typename TrackerTraits, typename Enable = void>
0145   struct QualityCutsT {};
0146
0147   template <typename TrackerTraits>
0148   struct QualityCutsT<TrackerTraits, pixelTopology::isPhase1Topology<TrackerTraits>> {
0149     using TrackSoAView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::View;
0150     using TrackSoAConstView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::ConstView;
0151     using tracksHelper = TracksUtilities<TrackerTraits>;
0152     // chi2 cut = chi2Scale * (chi2Coeff[0] + pT/GeV * (chi2Coeff[1] + pT/GeV * (chi2Coeff[2] + pT/GeV * chi2Coeff[3])))
0153     float chi2Coeff[4];
0154     float chi2MaxPt;  // GeV
0155     float chi2Scale;
0156
0157     struct Region {
0158       float maxTip;  // cm
0159       float minPt;   // GeV
0160       float maxZip;  // cm
0161     };
0162
0163     Region triplet;
0164     Region quadruplet;
0165
0166     __device__ __forceinline__ bool isHP(const TrackSoAConstView &tracks, int nHits, int it) const {
0167       // impose "region cuts" based on the fit results (phi, Tip, pt, cotan(theta)), Zip)
0168       // default cuts:
0169       //   - for triplets:    |Tip| < 0.3 cm, pT > 0.5 GeV, |Zip| < 12.0 cm
0170       //   - for quadruplets: |Tip| < 0.5 cm, pT > 0.3 GeV, |Zip| < 12.0 cm
0171       // (see CAHitNtupletGeneratorGPU.cc)
0172       auto const &region = (nHits > 3) ? quadruplet : triplet;
0173       return (std::abs(tracksHelper::tip(tracks, it)) < region.maxTip) and (tracks.pt(it) > region.minPt) and
0174              (std::abs(tracksHelper::zip(tracks, it)) < region.maxZip);
0175     }
0176
0177     __device__ __forceinline__ bool strictCut(const TrackSoAConstView &tracks, int it) const {
0178       auto roughLog = [](float x) {
0179         // max diff [0.5,12] at 1.25 0.16143
0180         // average diff  0.0662998
0181         union IF {
0182           uint32_t i;
0183           float f;
0184         };
0185         IF z;
0186         z.f = x;
0187         uint32_t lsb = 1 < 21;
0188         z.i += lsb;
0189         z.i >>= 21;
0190         auto f = z.i & 3;
0191         int ex = int(z.i >> 2) - 127;
0192
0193         // log2(1+0.25*f)
0194         // averaged over bins
0195         const float frac[4] = {0.160497f, 0.452172f, 0.694562f, 0.901964f};
0196         return float(ex) + frac[f];
0197       };
0198
0199       float pt = std::min<float>(tracks.pt(it), chi2MaxPt);
0200       float chi2Cut = chi2Scale * (chi2Coeff[0] + roughLog(pt) * chi2Coeff[1]);
0201       if (tracks.chi2(it) >= chi2Cut) {
0202 #ifdef NTUPLE_FIT_DEBUG
0203         printf("Bad chi2 %d pt %f eta %f chi2 %f\n", it, tracks.pt(it), tracks.eta(it), tracks.chi2(it));
0204 #endif
0205         return true;
0206       }
0207       return false;
0208     }
0209   };
0210
0211   template <typename TrackerTraits>
0212   struct QualityCutsT<TrackerTraits, pixelTopology::isPhase2Topology<TrackerTraits>> {
0213     using TrackSoAView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::View;
0214     using TrackSoAConstView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::ConstView;
0215     using tracksHelper = TracksUtilities<TrackerTraits>;
0216
0217     float maxChi2;
0218     float minPt;
0219     float maxTip;
0220     float maxZip;
0221
0222     __device__ __forceinline__ bool isHP(const TrackSoAConstView &tracks, int nHits, int it) const {
0223       return (std::abs(tracksHelper::tip(tracks, it)) < maxTip) and (tracks.pt(it) > minPt) and
0224              (std::abs(tracksHelper::zip(tracks, it)) < maxZip);
0225     }
0226     __device__ __forceinline__ bool strictCut(const TrackSoAConstView &tracks, int it) const {
0227       return tracks.chi2(it) >= maxChi2;
0228     }
0229   };
0230
0231 }  // namespace pixelTrack
0232
0233 template <typename TrackerTraits>
0234 using TrackLayout = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>;
0235 template <typename TrackerTraits>
0236 using TrackSoAView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::View;
0237 template <typename TrackerTraits>
0238 using TrackSoAConstView = typename TrackSoA<TrackerTraits>::template TrackSoALayout<>::ConstView;
0239
0240 template struct TracksUtilities<pixelTopology::Phase1>;
0241 template struct TracksUtilities<pixelTopology::Phase2>;
0242
0243 #endif