Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-20 02:32:11

0001 //#define BROKENLINE_DEBUG
0002 //#define BL_DUMP_HITS
0003 
0004 #include <cstdint>
0005 
0006 #include <alpaka/alpaka.hpp>
0007 
0008 #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h"
0009 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0010 #include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforDevice.h"
0011 #include "RecoTracker/PixelTrackFitting/interface/alpaka/BrokenLine.h"
0012 
0013 #include "HelixFit.h"
0014 
0015 template <typename TrackerTraits>
0016 using Tuples = typename reco::TrackSoA<TrackerTraits>::HitContainer;
0017 template <typename TrackerTraits>
0018 using OutputSoAView = reco::TrackSoAView<TrackerTraits>;
0019 template <typename TrackerTraits>
0020 using TupleMultiplicity = caStructures::TupleMultiplicityT<TrackerTraits>;
0021 
0022 // #define BL_DUMP_HITS
0023 
0024 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0025   template <int N, typename TrackerTraits>
0026   class Kernel_BLFastFit {
0027   public:
0028     template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0029     ALPAKA_FN_ACC void operator()(TAcc const &acc,
0030                                   Tuples<TrackerTraits> const *__restrict__ foundNtuplets,
0031                                   TupleMultiplicity<TrackerTraits> const *__restrict__ tupleMultiplicity,
0032                                   TrackingRecHitSoAConstView<TrackerTraits> hh,
0033                                   pixelCPEforDevice::ParamsOnDeviceT<TrackerTraits> const *__restrict__ cpeParams,
0034                                   typename TrackerTraits::tindex_type *__restrict__ ptkids,
0035                                   double *__restrict__ phits,
0036                                   float *__restrict__ phits_ge,
0037                                   double *__restrict__ pfast_fit,
0038                                   uint32_t nHitsL,
0039                                   uint32_t nHitsH,
0040                                   int32_t offset) const {
0041       constexpr uint32_t hitsInFit = N;
0042       constexpr auto invalidTkId = std::numeric_limits<typename TrackerTraits::tindex_type>::max();
0043 
0044       ALPAKA_ASSERT_ACC(hitsInFit <= nHitsL);
0045       ALPAKA_ASSERT_ACC(nHitsL <= nHitsH);
0046       ALPAKA_ASSERT_ACC(phits);
0047       ALPAKA_ASSERT_ACC(pfast_fit);
0048       ALPAKA_ASSERT_ACC(foundNtuplets);
0049       ALPAKA_ASSERT_ACC(tupleMultiplicity);
0050 
0051       // look in bin for this hit multiplicity
0052       int totTK = tupleMultiplicity->end(nHitsH) - tupleMultiplicity->begin(nHitsL);
0053       ALPAKA_ASSERT_ACC(totTK <= int(tupleMultiplicity->size()));
0054       ALPAKA_ASSERT_ACC(totTK >= 0);
0055 
0056 #ifdef BROKENLINE_DEBUG
0057       const uint32_t threadIdx(alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0u]);
0058       if (cms::alpakatools::once_per_grid(acc)) {
0059         printf("%d total Ntuple\n", tupleMultiplicity->size());
0060         printf("%d Ntuple of size %d/%d for %d hits to fit\n", totTK, nHitsL, nHitsH, hitsInFit);
0061       }
0062 #endif
0063       const auto nt = riemannFit::maxNumberOfConcurrentFits;
0064       for (auto local_idx : cms::alpakatools::uniform_elements(acc, nt)) {
0065         auto tuple_idx = local_idx + offset;
0066         if ((int)tuple_idx >= totTK) {
0067           ptkids[local_idx] = invalidTkId;
0068           break;
0069         }
0070         // get it from the ntuple container (one to one to helix)
0071         auto tkid = *(tupleMultiplicity->begin(nHitsL) + tuple_idx);
0072         ALPAKA_ASSERT_ACC(static_cast<int>(tkid) < foundNtuplets->nOnes());
0073 
0074         ptkids[local_idx] = tkid;
0075 
0076         auto nHits = foundNtuplets->size(tkid);
0077 
0078         ALPAKA_ASSERT_ACC(nHits >= nHitsL);
0079         ALPAKA_ASSERT_ACC(nHits <= nHitsH);
0080 
0081         riemannFit::Map3xNd<N> hits(phits + local_idx);
0082         riemannFit::Map4d fast_fit(pfast_fit + local_idx);
0083         riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
0084 
0085 #ifdef BL_DUMP_HITS
0086         auto &&done = alpaka::declareSharedVar<int, __COUNTER__>(acc);
0087         done = 0;
0088         alpaka::syncBlockThreads(acc);
0089         bool dump =
0090             (foundNtuplets->size(tkid) == 5 && 0 == alpaka::atomicAdd(acc, &done, 1, alpaka::hierarchy::Blocks{}));
0091 #endif
0092 
0093         // Prepare data structure
0094         auto const *hitId = foundNtuplets->begin(tkid);
0095 
0096         // #define YERR_FROM_DC
0097 #ifdef YERR_FROM_DC
0098         // try to compute more precise error in y
0099         auto dx = hh[hitId[hitsInFit - 1]].xGlobal() - hh[hitId[0]].xGlobal();
0100         auto dy = hh[hitId[hitsInFit - 1]].yGlobal() - hh[hitId[0]].yGlobal();
0101         auto dz = hh[hitId[hitsInFit - 1]].zGlobal() - hh[hitId[0]].zGlobal();
0102         float ux, uy, uz;
0103 #endif
0104 
0105         float incr = std::max(1.f, float(nHits) / float(hitsInFit));
0106         float n = 0;
0107         for (uint32_t i = 0; i < hitsInFit; ++i) {
0108           int j = int(n + 0.5f);  // round
0109           if (hitsInFit - 1 == i)
0110             j = nHits - 1;  // force last hit to ensure max lever arm.
0111           ALPAKA_ASSERT_ACC(j < int(nHits));
0112           n += incr;
0113           auto hit = hitId[j];
0114           float ge[6];
0115 
0116 #ifdef YERR_FROM_DC
0117           auto const &dp = cpeParams->detParams(hh.detectorIndex(hit));
0118           auto status = hh[hit].chargeAndStatus().status;
0119           int qbin = CPEFastParametrisation::kGenErrorQBins - 1 - status.qBin;
0120           ALPAKA_ASSERT_ACC(qbin >= 0 && qbin < 5);
0121           bool nok = (status.isBigY | status.isOneY);
0122           // compute cotanbeta and use it to recompute error
0123           dp.frame.rotation().multiply(dx, dy, dz, ux, uy, uz);
0124           auto cb = std::abs(uy / uz);
0125           int bin =
0126               int(cb * (float(phase1PixelTopology::pixelThickess) / float(phase1PixelTopology::pixelPitchY)) * 8.f) - 4;
0127           int low_value = 0;
0128           int high_value = CPEFastParametrisation::kNumErrorBins - 1;
0129           // return estimated bin value truncated to [0, 15]
0130           bin = std::clamp(bin, low_value, high_value);
0131           float yerr = dp.sigmay[bin] * 1.e-4f;  // toCM
0132           yerr *= dp.yfact[qbin];                // inflate
0133           yerr *= yerr;
0134           yerr += dp.apeYY;
0135           yerr = nok ? hh[hit].yerrLocal() : yerr;
0136           dp.frame.toGlobal(hh[hit].xerrLocal(), 0, yerr, ge);
0137 #else
0138           cpeParams->detParams(hh[hit].detectorIndex()).frame.toGlobal(hh[hit].xerrLocal(), 0, hh[hit].yerrLocal(), ge);
0139 #endif
0140 
0141 #ifdef BL_DUMP_HITS
0142           bool dump = foundNtuplets->size(tkid) == 5;
0143           if (dump) {
0144             printf("Track id %d %d Hit %d on %d\nGlobal: hits.col(%d) << %f,%f,%f\n",
0145                    local_idx,
0146                    tkid,
0147                    hit,
0148                    hh[hit].detectorIndex(),
0149                    i,
0150                    hh[hit].xGlobal(),
0151                    hh[hit].yGlobal(),
0152                    hh[hit].zGlobal());
0153             printf("Error: hits_ge.col(%d) << %e,%e,%e,%e,%e,%e\n", i, ge[0], ge[1], ge[2], ge[3], ge[4], ge[5]);
0154           }
0155 #endif
0156 
0157           hits.col(i) << hh[hit].xGlobal(), hh[hit].yGlobal(), hh[hit].zGlobal();
0158           hits_ge.col(i) << ge[0], ge[1], ge[2], ge[3], ge[4], ge[5];
0159         }
0160         brokenline::fastFit(acc, hits, fast_fit);
0161 
0162         // no NaN here....
0163         ALPAKA_ASSERT_ACC(fast_fit(0) == fast_fit(0));
0164         ALPAKA_ASSERT_ACC(fast_fit(1) == fast_fit(1));
0165         ALPAKA_ASSERT_ACC(fast_fit(2) == fast_fit(2));
0166         ALPAKA_ASSERT_ACC(fast_fit(3) == fast_fit(3));
0167       }
0168     }
0169   };
0170 
0171   template <int N, typename TrackerTraits>
0172   struct Kernel_BLFit {
0173   public:
0174     template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0175     ALPAKA_FN_ACC void operator()(TAcc const &acc,
0176                                   TupleMultiplicity<TrackerTraits> const *__restrict__ tupleMultiplicity,
0177                                   double bField,
0178                                   OutputSoAView<TrackerTraits> results_view,
0179                                   typename TrackerTraits::tindex_type const *__restrict__ ptkids,
0180                                   double *__restrict__ phits,
0181                                   float *__restrict__ phits_ge,
0182                                   double *__restrict__ pfast_fit) const {
0183       ALPAKA_ASSERT_ACC(results_view.pt());
0184       ALPAKA_ASSERT_ACC(results_view.eta());
0185       ALPAKA_ASSERT_ACC(results_view.chi2());
0186       ALPAKA_ASSERT_ACC(pfast_fit);
0187       constexpr auto invalidTkId = std::numeric_limits<typename TrackerTraits::tindex_type>::max();
0188 
0189       // same as above...
0190       // look in bin for this hit multiplicity
0191       const auto nt = riemannFit::maxNumberOfConcurrentFits;
0192       for (auto local_idx : cms::alpakatools::uniform_elements(acc, nt)) {
0193         if (invalidTkId == ptkids[local_idx])
0194           break;
0195         auto tkid = ptkids[local_idx];
0196 
0197         ALPAKA_ASSERT_ACC(tkid < TrackerTraits::maxNumberOfTuples);
0198 
0199         riemannFit::Map3xNd<N> hits(phits + local_idx);
0200         riemannFit::Map4d fast_fit(pfast_fit + local_idx);
0201         riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
0202 
0203         brokenline::PreparedBrokenLineData<N> data;
0204 
0205         brokenline::karimaki_circle_fit circle;
0206         riemannFit::LineFit line;
0207 
0208         brokenline::prepareBrokenLineData(acc, hits, fast_fit, bField, data);
0209         brokenline::lineFit(acc, hits_ge, fast_fit, bField, data, line);
0210         brokenline::circleFit(acc, hits, hits_ge, fast_fit, bField, data, circle);
0211 
0212         TracksUtilities<TrackerTraits>::copyFromCircle(
0213             results_view, circle.par, circle.cov, line.par, line.cov, 1.f / float(bField), tkid);
0214         results_view[tkid].pt() = float(bField) / float(std::abs(circle.par(2)));
0215         results_view[tkid].eta() = alpaka::math::asinh(acc, line.par(0));
0216         results_view[tkid].chi2() = (circle.chi2 + line.chi2) / (2 * N - 5);
0217 
0218 #ifdef BROKENLINE_DEBUG
0219         if (!(circle.chi2 >= 0) || !(line.chi2 >= 0))
0220           printf("kernelBLFit failed! %f/%f\n", circle.chi2, line.chi2);
0221         printf("kernelBLFit size %d for %d hits circle.par(0,1,2): %d %f,%f,%f\n",
0222                N,
0223                N,
0224                tkid,
0225                circle.par(0),
0226                circle.par(1),
0227                circle.par(2));
0228         printf("kernelBLHits line.par(0,1): %d %f,%f\n", tkid, line.par(0), line.par(1));
0229         printf("kernelBLHits chi2 cov %f/%f  %e,%e,%e,%e,%e\n",
0230                circle.chi2,
0231                line.chi2,
0232                circle.cov(0, 0),
0233                circle.cov(1, 1),
0234                circle.cov(2, 2),
0235                line.cov(0, 0),
0236                line.cov(1, 1));
0237 #endif
0238       }
0239     }
0240   };
0241 
0242   template <typename TrackerTraits>
0243   void HelixFit<TrackerTraits>::launchBrokenLineKernels(
0244       const TrackingRecHitSoAConstView<TrackerTraits> &hv,
0245       pixelCPEforDevice::ParamsOnDeviceT<TrackerTraits> const *cpeParams,
0246       uint32_t hitsInFit,
0247       uint32_t maxNumberOfTuples,
0248       Queue &queue) {
0249     ALPAKA_ASSERT_ACC(tuples_);
0250 
0251     uint32_t blockSize = 64;
0252     uint32_t numberOfBlocks = cms::alpakatools::divide_up_by(maxNumberOfConcurrentFits_, blockSize);
0253     const WorkDiv1D workDivTriplets = cms::alpakatools::make_workdiv<Acc1D>(numberOfBlocks, blockSize);
0254     const WorkDiv1D workDivQuadsPenta = cms::alpakatools::make_workdiv<Acc1D>(numberOfBlocks / 4, blockSize);
0255 
0256     //  Fit internals
0257     auto tkidDevice =
0258         cms::alpakatools::make_device_buffer<typename TrackerTraits::tindex_type[]>(queue, maxNumberOfConcurrentFits_);
0259     auto hitsDevice = cms::alpakatools::make_device_buffer<double[]>(
0260         queue, maxNumberOfConcurrentFits_ * sizeof(riemannFit::Matrix3xNd<6>) / sizeof(double));
0261     auto hits_geDevice = cms::alpakatools::make_device_buffer<float[]>(
0262         queue, maxNumberOfConcurrentFits_ * sizeof(riemannFit::Matrix6xNf<6>) / sizeof(float));
0263     auto fast_fit_resultsDevice = cms::alpakatools::make_device_buffer<double[]>(
0264         queue, maxNumberOfConcurrentFits_ * sizeof(riemannFit::Vector4d) / sizeof(double));
0265 
0266     for (uint32_t offset = 0; offset < maxNumberOfTuples; offset += maxNumberOfConcurrentFits_) {
0267       // fit triplets
0268 
0269       alpaka::exec<Acc1D>(queue,
0270                           workDivTriplets,
0271                           Kernel_BLFastFit<3, TrackerTraits>{},
0272                           tuples_,
0273                           tupleMultiplicity_,
0274                           hv,
0275                           cpeParams,
0276                           tkidDevice.data(),
0277                           hitsDevice.data(),
0278                           hits_geDevice.data(),
0279                           fast_fit_resultsDevice.data(),
0280                           3,
0281                           3,
0282                           offset);
0283 
0284       alpaka::exec<Acc1D>(queue,
0285                           workDivTriplets,
0286                           Kernel_BLFit<3, TrackerTraits>{},
0287                           tupleMultiplicity_,
0288                           bField_,
0289                           outputSoa_,
0290                           tkidDevice.data(),
0291                           hitsDevice.data(),
0292                           hits_geDevice.data(),
0293                           fast_fit_resultsDevice.data());
0294 
0295       if (fitNas4_) {
0296         // fit all as 4
0297         riemannFit::rolling_fits<4, TrackerTraits::maxHitsOnTrack, 1>([this,
0298                                                                        &hv,
0299                                                                        &cpeParams,
0300                                                                        &tkidDevice,
0301                                                                        &hitsDevice,
0302                                                                        &hits_geDevice,
0303                                                                        &fast_fit_resultsDevice,
0304                                                                        &offset,
0305                                                                        &queue,
0306                                                                        &workDivQuadsPenta](auto i) {
0307           alpaka::exec<Acc1D>(queue,
0308                               workDivQuadsPenta,
0309                               Kernel_BLFastFit<4, TrackerTraits>{},
0310                               tuples_,
0311                               tupleMultiplicity_,
0312                               hv,
0313                               cpeParams,
0314                               tkidDevice.data(),
0315                               hitsDevice.data(),
0316                               hits_geDevice.data(),
0317                               fast_fit_resultsDevice.data(),
0318                               4,
0319                               4,
0320                               offset);
0321 
0322           alpaka::exec<Acc1D>(queue,
0323                               workDivQuadsPenta,
0324                               Kernel_BLFit<4, TrackerTraits>{},
0325                               tupleMultiplicity_,
0326                               bField_,
0327                               outputSoa_,
0328                               tkidDevice.data(),
0329                               hitsDevice.data(),
0330                               hits_geDevice.data(),
0331                               fast_fit_resultsDevice.data());
0332         });
0333 
0334       } else {
0335         riemannFit::rolling_fits<4, TrackerTraits::maxHitsOnTrackForFullFit, 1>([this,
0336                                                                                  &hv,
0337                                                                                  &cpeParams,
0338                                                                                  &tkidDevice,
0339                                                                                  &hitsDevice,
0340                                                                                  &hits_geDevice,
0341                                                                                  &fast_fit_resultsDevice,
0342                                                                                  &offset,
0343                                                                                  &queue,
0344                                                                                  &workDivQuadsPenta](auto i) {
0345           alpaka::exec<Acc1D>(queue,
0346                               workDivQuadsPenta,
0347                               Kernel_BLFastFit<i, TrackerTraits>{},
0348                               tuples_,
0349                               tupleMultiplicity_,
0350                               hv,
0351                               cpeParams,
0352                               tkidDevice.data(),
0353                               hitsDevice.data(),
0354                               hits_geDevice.data(),
0355                               fast_fit_resultsDevice.data(),
0356                               i,
0357                               i,
0358                               offset);
0359 
0360           alpaka::exec<Acc1D>(queue,
0361                               workDivQuadsPenta,
0362                               Kernel_BLFit<i, TrackerTraits>{},
0363                               tupleMultiplicity_,
0364                               bField_,
0365                               outputSoa_,
0366                               tkidDevice.data(),
0367                               hitsDevice.data(),
0368                               hits_geDevice.data(),
0369                               fast_fit_resultsDevice.data());
0370         });
0371 
0372         static_assert(TrackerTraits::maxHitsOnTrackForFullFit < TrackerTraits::maxHitsOnTrack);
0373 
0374         //Fit all the rest using the maximum from previous call
0375         alpaka::exec<Acc1D>(queue,
0376                             workDivQuadsPenta,
0377                             Kernel_BLFastFit<TrackerTraits::maxHitsOnTrackForFullFit, TrackerTraits>{},
0378                             tuples_,
0379                             tupleMultiplicity_,
0380                             hv,
0381                             cpeParams,
0382                             tkidDevice.data(),
0383                             hitsDevice.data(),
0384                             hits_geDevice.data(),
0385                             fast_fit_resultsDevice.data(),
0386                             TrackerTraits::maxHitsOnTrackForFullFit,
0387                             TrackerTraits::maxHitsOnTrack - 1,
0388                             offset);
0389 
0390         alpaka::exec<Acc1D>(queue,
0391                             workDivQuadsPenta,
0392                             Kernel_BLFit<TrackerTraits::maxHitsOnTrackForFullFit, TrackerTraits>{},
0393                             tupleMultiplicity_,
0394                             bField_,
0395                             outputSoa_,
0396                             tkidDevice.data(),
0397                             hitsDevice.data(),
0398                             hits_geDevice.data(),
0399                             fast_fit_resultsDevice.data());
0400       }
0401 
0402     }  // loop on concurrent fits
0403   }
0404 
0405   template class HelixFit<pixelTopology::Phase1>;
0406   template class HelixFit<pixelTopology::Phase2>;
0407   template class HelixFit<pixelTopology::HIonPhase1>;
0408 
0409 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE