interface/alpaka/TrackUtilities.h

0001 #ifndef DataFormats_TrackSoA_interface_alpaka_TrackUtilities_h
0002 #define DataFormats_TrackSoA_interface_alpaka_TrackUtilities_h
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <cstdint>
0007
0008 #include <alpaka/alpaka.hpp>
0009
0010 #include "DataFormats/TrackSoA/interface/TracksSoA.h"
0011 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0012
0013 // Methods that operate on View and ConstView of the TrackSoA, and cannot be class methods.
0014 template <typename TrackerTraits>
0015 struct TracksUtilities {
0016   using TrackSoAView = typename reco::TrackSoA<TrackerTraits>::template Layout<>::View;
0017   using TrackSoAConstView = typename reco::TrackSoA<TrackerTraits>::template Layout<>::ConstView;
0018   using hindex_type = typename reco::TrackSoA<TrackerTraits>::hindex_type;
0019
0020   // state at the beam spot: { phi, tip, 1/pt, cotan(theta), zip }
0021
0022   template <typename V3, typename M3, typename V2, typename M2>
0023   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE static constexpr void copyFromCircle(
0024       TrackSoAView &tracks, V3 const &cp, M3 const &ccov, V2 const &lp, M2 const &lcov, float b, int32_t i) {
0025     tracks[i].state() << cp.template cast<float>(), lp.template cast<float>();
0026
0027     tracks[i].state()(2) = tracks[i].state()(2) * b;
0028     auto cov = tracks[i].covariance();
0029     cov(0) = ccov(0, 0);
0030     cov(1) = ccov(0, 1);
0031     cov(2) = b * float(ccov(0, 2));
0032     cov(4) = cov(3) = 0;
0033     cov(5) = ccov(1, 1);
0034     cov(6) = b * float(ccov(1, 2));
0035     cov(8) = cov(7) = 0;
0036     cov(9) = b * b * float(ccov(2, 2));
0037     cov(11) = cov(10) = 0;
0038     cov(12) = lcov(0, 0);
0039     cov(13) = lcov(0, 1);
0040     cov(14) = lcov(1, 1);
0041   }
0042
0043   template <typename V5, typename M5>
0044   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE static constexpr void copyFromDense(TrackSoAView &tracks,
0045                                                                           V5 const &v,
0046                                                                           M5 const &cov,
0047                                                                           int32_t i) {
0048     tracks[i].state() = v.template cast<float>();
0049     for (int j = 0, ind = 0; j < 5; ++j)
0050       for (auto k = j; k < 5; ++k)
0051         tracks[i].covariance()(ind++) = cov(j, k);
0052   }
0053
0054   template <typename V5, typename M5>
0055   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE static constexpr void copyToDense(const TrackSoAConstView &tracks,
0056                                                                         V5 &v,
0057                                                                         M5 &cov,
0058                                                                         int32_t i) {
0059     v = tracks[i].state().template cast<typename V5::Scalar>();
0060     for (int j = 0, ind = 0; j < 5; ++j) {
0061       cov(j, j) = tracks[i].covariance()(ind++);
0062       for (auto k = j + 1; k < 5; ++k)
0063         cov(k, j) = cov(j, k) = tracks[i].covariance()(ind++);
0064     }
0065   }
0066
0067   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE static constexpr int computeNumberOfLayers(const TrackSoAConstView &tracks,
0068                                                                                  int32_t i) {
0069     auto pdet = tracks.detIndices().begin(i);
0070     int nl = 1;
0071     auto ol = pixelTopology::getLayer<TrackerTraits>(*pdet);
0072     for (; pdet < tracks.detIndices().end(i); ++pdet) {
0073       auto il = pixelTopology::getLayer<TrackerTraits>(*pdet);
0074       if (il != ol)
0075         ++nl;
0076       ol = il;
0077     }
0078     return nl;
0079   }
0080
0081   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE static constexpr int nHits(const TrackSoAConstView &tracks, int i) {
0082     return tracks.detIndices().size(i);
0083   }
0084 };
0085
0086 namespace pixelTrack {
0087
0088   template <typename TrackerTraits, typename Enable = void>
0089   struct QualityCutsT {};
0090
0091   template <typename TrackerTraits>
0092   struct QualityCutsT<TrackerTraits, pixelTopology::isPhase1Topology<TrackerTraits>> {
0093     using TrackSoAView = typename reco::TrackSoA<TrackerTraits>::template Layout<>::View;
0094     using TrackSoAConstView = typename reco::TrackSoA<TrackerTraits>::template Layout<>::ConstView;
0095     float chi2Coeff[4];
0096     float chi2MaxPt;  // GeV
0097     float chi2Scale;
0098
0099     struct Region {
0100       float maxTip;  // cm
0101       float minPt;   // GeV
0102       float maxZip;  // cm
0103     };
0104
0105     Region triplet;
0106     Region quadruplet;
0107
0108     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool isHP(const TrackSoAConstView &tracks, int nHits, int it) const {
0109       // impose "region cuts" based on the fit results (phi, Tip, pt, cotan(theta)), Zip)
0110       // default cuts:
0111       //   - for triplets:    |Tip| < 0.3 cm, pT > 0.5 GeV, |Zip| < 12.0 cm
0112       //   - for quadruplets: |Tip| < 0.5 cm, pT > 0.3 GeV, |Zip| < 12.0 cm
0113       // (see CAHitNtupletGeneratorGPU.cc)
0114       auto const &region = (nHits > 3) ? quadruplet : triplet;
0115       return (std::abs(reco::tip(tracks, it)) < region.maxTip) and (tracks.pt(it) > region.minPt) and
0116              (std::abs(reco::zip(tracks, it)) < region.maxZip);
0117     }
0118
0119     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool strictCut(const TrackSoAConstView &tracks, int it) const {
0120       auto roughLog = [](float x) {
0121         // max diff [0.5,12] at 1.25 0.16143
0122         // average diff  0.0662998
0123         union IF {
0124           uint32_t i;
0125           float f;
0126         };
0127         IF z;
0128         z.f = x;
0129         uint32_t lsb = 1 < 21;
0130         z.i += lsb;
0131         z.i >>= 21;
0132         auto f = z.i & 3;
0133         int ex = int(z.i >> 2) - 127;
0134
0135         // log2(1+0.25*f)
0136         // averaged over bins
0137         const float frac[4] = {0.160497f, 0.452172f, 0.694562f, 0.901964f};
0138         return float(ex) + frac[f];
0139       };
0140
0141       float pt = std::min<float>(tracks.pt(it), chi2MaxPt);
0142       float chi2Cut = chi2Scale * (chi2Coeff[0] + roughLog(pt) * chi2Coeff[1]);
0143       if (tracks.chi2(it) >= chi2Cut) {
0144 #ifdef NTUPLE_FIT_DEBUG
0145         printf("Bad chi2 %d pt %f eta %f chi2 %f\n", it, tracks.pt(it), tracks.eta(it), tracks.chi2(it));
0146 #endif
0147         return true;
0148       }
0149       return false;
0150     }
0151   };
0152
0153   template <typename TrackerTraits>
0154   struct QualityCutsT<TrackerTraits, pixelTopology::isPhase2Topology<TrackerTraits>> {
0155     using TrackSoAView = typename reco::TrackSoA<TrackerTraits>::template Layout<>::View;
0156     using TrackSoAConstView = typename reco::TrackSoA<TrackerTraits>::template Layout<>::ConstView;
0157
0158     float maxChi2;
0159     float minPt;
0160     float maxTip;
0161     float maxZip;
0162
0163     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool isHP(const TrackSoAConstView &tracks, int nHits, int it) const {
0164       return (std::abs(reco::tip(tracks, it)) < maxTip) and (tracks.pt(it) > minPt) and
0165              (std::abs(reco::zip(tracks, it)) < maxZip);
0166     }
0167     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool strictCut(const TrackSoAConstView &tracks, int it) const {
0168       return tracks.chi2(it) >= maxChi2;
0169     }
0170   };
0171
0172 }  // namespace pixelTrack
0173
0174 // TODO: Should those be placed in the ALPAKA_ACCELERATOR_NAMESPACE
0175 template struct TracksUtilities<pixelTopology::Phase1>;
0176 template struct TracksUtilities<pixelTopology::Phase2>;
0177
0178 #endif  // DataFormats_TrackSoA_interface_alpaka_TrackUtilities_h