Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-17 02:47:07

0001 //
0002 // Author: Felice Pantaleo, CERN
0003 //
0004 
0005 #include <cstdint>
0006 
0007 #include <cuda_runtime.h>
0008 
0009 #include "CUDADataFormats/TrackingRecHit/interface/TrackingRecHit2DHeterogeneous.h"
0010 #include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
0011 #include "HeterogeneousCore/CUDAUtilities/interface/cuda_assert.h"
0012 #include "RecoLocalTracker/SiPixelRecHits/interface/pixelCPEforGPU.h"
0013 #include "RecoPixelVertexing/PixelTrackFitting/interface/RiemannFit.h"
0014 
0015 #include "HelixFitOnGPU.h"
0016 
0017 using HitsOnGPU = TrackingRecHit2DSOAView;
0018 using Tuples = pixelTrack::HitContainer;
0019 using OutputSoA = pixelTrack::TrackSoA;
0020 
0021 template <int N>
0022 __global__ void kernel_FastFit(Tuples const *__restrict__ foundNtuplets,
0023                                caConstants::TupleMultiplicity const *__restrict__ tupleMultiplicity,
0024                                uint32_t nHits,
0025                                HitsOnGPU const *__restrict__ hhp,
0026                                double *__restrict__ phits,
0027                                float *__restrict__ phits_ge,
0028                                double *__restrict__ pfast_fit,
0029                                uint32_t offset) {
0030   constexpr uint32_t hitsInFit = N;
0031 
0032   assert(hitsInFit <= nHits);
0033 
0034   assert(pfast_fit);
0035   assert(foundNtuplets);
0036   assert(tupleMultiplicity);
0037 
0038   // look in bin for this hit multiplicity
0039   auto local_start = blockIdx.x * blockDim.x + threadIdx.x;
0040 
0041 #ifdef RIEMANN_DEBUG
0042   if (0 == local_start)
0043     printf("%d Ntuple of size %d for %d hits to fit\n", tupleMultiplicity->size(nHits), nHits, hitsInFit);
0044 #endif
0045 
0046   for (int local_idx = local_start, nt = riemannFit::maxNumberOfConcurrentFits; local_idx < nt;
0047        local_idx += gridDim.x * blockDim.x) {
0048     auto tuple_idx = local_idx + offset;
0049     if (tuple_idx >= tupleMultiplicity->size(nHits))
0050       break;
0051 
0052     // get it from the ntuple container (one to one to helix)
0053     auto tkid = *(tupleMultiplicity->begin(nHits) + tuple_idx);
0054     assert(tkid < foundNtuplets->nOnes());
0055 
0056     assert(foundNtuplets->size(tkid) == nHits);
0057 
0058     riemannFit::Map3xNd<N> hits(phits + local_idx);
0059     riemannFit::Map4d fast_fit(pfast_fit + local_idx);
0060     riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
0061 
0062     // Prepare data structure
0063     auto const *hitId = foundNtuplets->begin(tkid);
0064     for (unsigned int i = 0; i < hitsInFit; ++i) {
0065       auto hit = hitId[i];
0066       // printf("Hit global: %f,%f,%f\n", hhp->xg_d[hit],hhp->yg_d[hit],hhp->zg_d[hit]);
0067       float ge[6];
0068       hhp->cpeParams()
0069           .detParams(hhp->detectorIndex(hit))
0070           .frame.toGlobal(hhp->xerrLocal(hit), 0, hhp->yerrLocal(hit), ge);
0071       // printf("Error: %d: %f,%f,%f,%f,%f,%f\n",hhp->detInd_d[hit],ge[0],ge[1],ge[2],ge[3],ge[4],ge[5]);
0072 
0073       hits.col(i) << hhp->xGlobal(hit), hhp->yGlobal(hit), hhp->zGlobal(hit);
0074       hits_ge.col(i) << ge[0], ge[1], ge[2], ge[3], ge[4], ge[5];
0075     }
0076     riemannFit::fastFit(hits, fast_fit);
0077 
0078     // no NaN here....
0079     assert(fast_fit(0) == fast_fit(0));
0080     assert(fast_fit(1) == fast_fit(1));
0081     assert(fast_fit(2) == fast_fit(2));
0082     assert(fast_fit(3) == fast_fit(3));
0083   }
0084 }
0085 
0086 template <int N>
0087 __global__ void kernel_CircleFit(caConstants::TupleMultiplicity const *__restrict__ tupleMultiplicity,
0088                                  uint32_t nHits,
0089                                  double bField,
0090                                  double *__restrict__ phits,
0091                                  float *__restrict__ phits_ge,
0092                                  double *__restrict__ pfast_fit_input,
0093                                  riemannFit::CircleFit *circle_fit,
0094                                  uint32_t offset) {
0095   assert(circle_fit);
0096   assert(N <= nHits);
0097 
0098   // same as above...
0099 
0100   // look in bin for this hit multiplicity
0101   auto local_start = blockIdx.x * blockDim.x + threadIdx.x;
0102   for (int local_idx = local_start, nt = riemannFit::maxNumberOfConcurrentFits; local_idx < nt;
0103        local_idx += gridDim.x * blockDim.x) {
0104     auto tuple_idx = local_idx + offset;
0105     if (tuple_idx >= tupleMultiplicity->size(nHits))
0106       break;
0107 
0108     riemannFit::Map3xNd<N> hits(phits + local_idx);
0109     riemannFit::Map4d fast_fit(pfast_fit_input + local_idx);
0110     riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
0111 
0112     riemannFit::VectorNd<N> rad = (hits.block(0, 0, 2, N).colwise().norm());
0113 
0114     riemannFit::Matrix2Nd<N> hits_cov = riemannFit::Matrix2Nd<N>::Zero();
0115     riemannFit::loadCovariance2D(hits_ge, hits_cov);
0116 
0117     circle_fit[local_idx] = riemannFit::circleFit(hits.block(0, 0, 2, N), hits_cov, fast_fit, rad, bField, true);
0118 
0119 #ifdef RIEMANN_DEBUG
0120 //    auto tkid = *(tupleMultiplicity->begin(nHits) + tuple_idx);
0121 //  printf("kernelCircleFit circle.par(0,1,2): %d %f,%f,%f\n", tkid,
0122 //         circle_fit[local_idx].par(0), circle_fit[local_idx].par(1), circle_fit[local_idx].par(2));
0123 #endif
0124   }
0125 }
0126 
0127 template <int N>
0128 __global__ void kernel_LineFit(caConstants::TupleMultiplicity const *__restrict__ tupleMultiplicity,
0129                                uint32_t nHits,
0130                                double bField,
0131                                OutputSoA *results,
0132                                double *__restrict__ phits,
0133                                float *__restrict__ phits_ge,
0134                                double *__restrict__ pfast_fit_input,
0135                                riemannFit::CircleFit *__restrict__ circle_fit,
0136                                uint32_t offset) {
0137   assert(results);
0138   assert(circle_fit);
0139   assert(N <= nHits);
0140 
0141   // same as above...
0142 
0143   // look in bin for this hit multiplicity
0144   auto local_start = (blockIdx.x * blockDim.x + threadIdx.x);
0145   for (int local_idx = local_start, nt = riemannFit::maxNumberOfConcurrentFits; local_idx < nt;
0146        local_idx += gridDim.x * blockDim.x) {
0147     auto tuple_idx = local_idx + offset;
0148     if (tuple_idx >= tupleMultiplicity->size(nHits))
0149       break;
0150 
0151     // get it for the ntuple container (one to one to helix)
0152     auto tkid = *(tupleMultiplicity->begin(nHits) + tuple_idx);
0153 
0154     riemannFit::Map3xNd<N> hits(phits + local_idx);
0155     riemannFit::Map4d fast_fit(pfast_fit_input + local_idx);
0156     riemannFit::Map6xNf<N> hits_ge(phits_ge + local_idx);
0157 
0158     auto const &line_fit = riemannFit::lineFit(hits, hits_ge, circle_fit[local_idx], fast_fit, bField, true);
0159 
0160     riemannFit::fromCircleToPerigee(circle_fit[local_idx]);
0161 
0162     results->stateAtBS.copyFromCircle(
0163         circle_fit[local_idx].par, circle_fit[local_idx].cov, line_fit.par, line_fit.cov, 1.f / float(bField), tkid);
0164     results->pt(tkid) = bField / std::abs(circle_fit[local_idx].par(2));
0165     results->eta(tkid) = asinhf(line_fit.par(0));
0166     results->chi2(tkid) = (circle_fit[local_idx].chi2 + line_fit.chi2) / (2 * N - 5);
0167 
0168 #ifdef RIEMANN_DEBUG
0169     printf("kernelLineFit size %d for %d hits circle.par(0,1,2): %d %f,%f,%f\n",
0170            N,
0171            nHits,
0172            tkid,
0173            circle_fit[local_idx].par(0),
0174            circle_fit[local_idx].par(1),
0175            circle_fit[local_idx].par(2));
0176     printf("kernelLineFit line.par(0,1): %d %f,%f\n", tkid, line_fit.par(0), line_fit.par(1));
0177     printf("kernelLineFit chi2 cov %f/%f %e,%e,%e,%e,%e\n",
0178            circle_fit[local_idx].chi2,
0179            line_fit.chi2,
0180            circle_fit[local_idx].cov(0, 0),
0181            circle_fit[local_idx].cov(1, 1),
0182            circle_fit[local_idx].cov(2, 2),
0183            line_fit.cov(0, 0),
0184            line_fit.cov(1, 1));
0185 #endif
0186   }
0187 }