PixelSeeding/plugins/CAHitNtupletGeneratorOnGPU.cc

0001 //
0002 // Original Author: Felice Pantaleo, CERN
0003 //
0004
0005 // #define GPU_DEBUG
0006 // #define DUMP_GPU_TK_TUPLES
0007
0008 #include <array>
0009 #include <cassert>
0010 #include <functional>
0011 #include <vector>
0012
0013 #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousDevice.h"
0014 #include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousHost.h"
0015 #include "CUDADataFormats/TrackingRecHit/interface/TrackingRecHitSoADevice.h"
0016 #include "CUDADataFormats/TrackingRecHit/interface/TrackingRecHitSoAHost.h"
0017 #include "DataFormats/Common/interface/Handle.h"
0018 #include "FWCore/Framework/interface/ConsumesCollector.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "FWCore/Utilities/interface/EDMException.h"
0024 #include "FWCore/Utilities/interface/isFinite.h"
0025 #include "HeterogeneousCore/CUDAServices/interface/CUDAInterface.h"
0026 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0027
0028 #include "CAHitNtupletGeneratorOnGPU.h"
0029
0030 namespace {
0031
0032   using namespace caHitNtupletGenerator;
0033   using namespace gpuPixelDoublets;
0034   using namespace pixelTopology;
0035   using namespace pixelTrack;
0036
0037   template <typename T>
0038   T sqr(T x) {
0039     return x * x;
0040   }
0041
0042   //Common Params
0043   AlgoParams makeCommonParams(edm::ParameterSet const& cfg) {
0044     return AlgoParams({cfg.getParameter<bool>("onGPU"),
0045                        cfg.getParameter<unsigned int>("minHitsForSharingCut"),
0046                        cfg.getParameter<bool>("useRiemannFit"),
0047                        cfg.getParameter<bool>("fitNas4"),
0048                        cfg.getParameter<bool>("includeJumpingForwardDoublets"),
0049                        cfg.getParameter<bool>("earlyFishbone"),
0050                        cfg.getParameter<bool>("lateFishbone"),
0051                        cfg.getParameter<bool>("fillStatistics"),
0052                        cfg.getParameter<bool>("doSharedHitCut"),
0053                        cfg.getParameter<bool>("dupPassThrough"),
0054                        cfg.getParameter<bool>("useSimpleTripletCleaner")});
0055   }
0056
0057   //This is needed to have the partial specialization for  isPhase1Topology/isPhase2Topology
0058   template <typename TrackerTraits, typename Enable = void>
0059   struct topologyCuts {};
0060
0061   template <typename TrackerTraits>
0062   struct topologyCuts<TrackerTraits, isPhase1Topology<TrackerTraits>> {
0063     static constexpr CAParamsT<TrackerTraits> makeCACuts(edm::ParameterSet const& cfg) {
0064       return CAParamsT<TrackerTraits>{{
0065           cfg.getParameter<unsigned int>("maxNumberOfDoublets"),
0066           cfg.getParameter<unsigned int>("minHitsPerNtuplet"),
0067           static_cast<float>(cfg.getParameter<double>("ptmin")),
0068           static_cast<float>(cfg.getParameter<double>("CAThetaCutBarrel")),
0069           static_cast<float>(cfg.getParameter<double>("CAThetaCutForward")),
0070           static_cast<float>(cfg.getParameter<double>("hardCurvCut")),
0071           static_cast<float>(cfg.getParameter<double>("dcaCutInnerTriplet")),
0072           static_cast<float>(cfg.getParameter<double>("dcaCutOuterTriplet")),
0073       }};
0074     };
0075
0076     static constexpr pixelTrack::QualityCutsT<TrackerTraits> makeQualityCuts(edm::ParameterSet const& pset) {
0077       auto coeff = pset.getParameter<std::array<double, 2>>("chi2Coeff");
0078       auto ptMax = pset.getParameter<double>("chi2MaxPt");
0079
0080       coeff[1] = (coeff[1] - coeff[0]) / log2(ptMax);
0081       return pixelTrack::QualityCutsT<TrackerTraits>{// polynomial coefficients for the pT-dependent chi2 cut
0082                                                      {(float)coeff[0], (float)coeff[1], 0.f, 0.f},
0083                                                      // max pT used to determine the chi2 cut
0084                                                      (float)ptMax,
0085                                                      // chi2 scale factor: 8 for broken line fit, ?? for Riemann fit
0086                                                      (float)pset.getParameter<double>("chi2Scale"),
0087                                                      // regional cuts for triplets
0088                                                      {(float)pset.getParameter<double>("tripletMaxTip"),
0089                                                       (float)pset.getParameter<double>("tripletMinPt"),
0090                                                       (float)pset.getParameter<double>("tripletMaxZip")},
0091                                                      // regional cuts for quadruplets
0092                                                      {(float)pset.getParameter<double>("quadrupletMaxTip"),
0093                                                       (float)pset.getParameter<double>("quadrupletMinPt"),
0094                                                       (float)pset.getParameter<double>("quadrupletMaxZip")}};
0095     }
0096   };
0097
0098   template <typename TrackerTraits>
0099   struct topologyCuts<TrackerTraits, isPhase2Topology<TrackerTraits>> {
0100     static constexpr CAParamsT<TrackerTraits> makeCACuts(edm::ParameterSet const& cfg) {
0101       return CAParamsT<TrackerTraits>{{cfg.getParameter<unsigned int>("maxNumberOfDoublets"),
0102                                        cfg.getParameter<unsigned int>("minHitsPerNtuplet"),
0103                                        static_cast<float>(cfg.getParameter<double>("ptmin")),
0104                                        static_cast<float>(cfg.getParameter<double>("CAThetaCutBarrel")),
0105                                        static_cast<float>(cfg.getParameter<double>("CAThetaCutForward")),
0106                                        static_cast<float>(cfg.getParameter<double>("hardCurvCut")),
0107                                        static_cast<float>(cfg.getParameter<double>("dcaCutInnerTriplet")),
0108                                        static_cast<float>(cfg.getParameter<double>("dcaCutOuterTriplet"))},
0109                                       {(bool)cfg.getParameter<bool>("includeFarForwards")}};
0110     }
0111
0112     static constexpr pixelTrack::QualityCutsT<TrackerTraits> makeQualityCuts(edm::ParameterSet const& pset) {
0113       return pixelTrack::QualityCutsT<TrackerTraits>{
0114           static_cast<float>(pset.getParameter<double>("maxChi2")),
0115           static_cast<float>(pset.getParameter<double>("minPt")),
0116           static_cast<float>(pset.getParameter<double>("maxTip")),
0117           static_cast<float>(pset.getParameter<double>("maxZip")),
0118       };
0119     }
0120   };
0121
0122   //Cell Cuts, as they are the cuts have the same logic for Phase2 and Phase1
0123   //keeping them separate would allow further differentiation in the future
0124   //moving them to topologyCuts and using the same syntax
0125   template <typename TrackerTraits>
0126   CellCutsT<TrackerTraits> makeCellCuts(edm::ParameterSet const& cfg) {
0127     return CellCutsT<TrackerTraits>{cfg.getParameter<bool>("doClusterCut"),
0128                                     cfg.getParameter<bool>("doZ0Cut"),
0129                                     cfg.getParameter<bool>("doPtCut"),
0130                                     cfg.getParameter<bool>("idealConditions"),
0131                                     (float)cfg.getParameter<double>("z0Cut"),
0132                                     (float)cfg.getParameter<double>("ptCut"),
0133                                     cfg.getParameter<std::vector<int>>("phiCuts")};
0134   }
0135
0136 }  // namespace
0137
0138 using namespace std;
0139
0140 template <typename TrackerTraits>
0141 CAHitNtupletGeneratorOnGPU<TrackerTraits>::CAHitNtupletGeneratorOnGPU(const edm::ParameterSet& cfg,
0142                                                                       edm::ConsumesCollector& iC)
0143     : m_params(makeCommonParams(cfg),
0144                makeCellCuts<TrackerTraits>(cfg),
0145                topologyCuts<TrackerTraits>::makeQualityCuts(cfg.getParameterSet("trackQualityCuts")),
0146                topologyCuts<TrackerTraits>::makeCACuts(cfg)) {
0147 #ifdef DUMP_GPU_TK_TUPLES
0148   printf("TK: %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s\n",
0149          "tid",
0150          "qual",
0151          "nh",
0152          "nl",
0153          "charge",
0154          "pt",
0155          "eta",
0156          "phi",
0157          "tip",
0158          "zip",
0159          "chi2",
0160          "h1",
0161          "h2",
0162          "h3",
0163          "h4",
0164          "h5",
0165          "hn");
0166 #endif
0167 }
0168
0169 template <typename TrackerTraits>
0170 void CAHitNtupletGeneratorOnGPU<TrackerTraits>::fillDescriptions(edm::ParameterSetDescription& desc) {
0171   fillDescriptionsCommon(desc);
0172   edm::LogWarning("CAHitNtupletGeneratorOnGPU::fillDescriptions")
0173       << "Note: this fillDescriptions is a dummy one. Most probably you are missing some parameters. \n"
0174          "please implement your TrackerTraits descriptions in CAHitNtupletGeneratorOnGPU. \n";
0175 }
0176
0177 template <>
0178 void CAHitNtupletGeneratorOnGPU<pixelTopology::Phase1>::fillDescriptions(edm::ParameterSetDescription& desc) {
0179   fillDescriptionsCommon(desc);
0180
0181   desc.add<bool>("idealConditions", true);
0182   desc.add<bool>("includeJumpingForwardDoublets", false);
0183   desc.add<double>("z0Cut", 12.0);
0184   desc.add<double>("ptCut", 0.5);
0185
0186   edm::ParameterSetDescription trackQualityCuts;
0187   trackQualityCuts.add<double>("chi2MaxPt", 10.)->setComment("max pT used to determine the pT-dependent chi2 cut");
0188   trackQualityCuts.add<std::vector<double>>("chi2Coeff", {0.9, 1.8})->setComment("chi2 at 1GeV and at ptMax above");
0189   trackQualityCuts.add<double>("chi2Scale", 8.)
0190       ->setComment(
0191           "Factor to multiply the pT-dependent chi2 cut (currently: 8 for the broken line fit, ?? for the Riemann "
0192           "fit)");
0193   trackQualityCuts.add<double>("tripletMinPt", 0.5)->setComment("Min pT for triplets, in GeV");
0194   trackQualityCuts.add<double>("tripletMaxTip", 0.3)->setComment("Max |Tip| for triplets, in cm");
0195   trackQualityCuts.add<double>("tripletMaxZip", 12.)->setComment("Max |Zip| for triplets, in cm");
0196   trackQualityCuts.add<double>("quadrupletMinPt", 0.3)->setComment("Min pT for quadruplets, in GeV");
0197   trackQualityCuts.add<double>("quadrupletMaxTip", 0.5)->setComment("Max |Tip| for quadruplets, in cm");
0198   trackQualityCuts.add<double>("quadrupletMaxZip", 12.)->setComment("Max |Zip| for quadruplets, in cm");
0199
0200   desc.add<std::vector<int>>(
0201           "phiCuts", std::vector<int>(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts)))
0202       ->setComment("Cuts in phi for cells");
0203
0204   desc.add<edm::ParameterSetDescription>("trackQualityCuts", trackQualityCuts)
0205       ->setComment(
0206           "Quality cuts based on the results of the track fit:\n  - apply a pT-dependent chi2 cut;\n  - apply \"region "
0207           "cuts\" based on the fit results (pT, Tip, Zip).");
0208 }
0209
0210 template <>
0211 void CAHitNtupletGeneratorOnGPU<pixelTopology::HIonPhase1>::fillDescriptions(edm::ParameterSetDescription& desc) {
0212   fillDescriptionsCommon(desc);
0213
0214   desc.add<bool>("idealConditions", false);
0215   desc.add<bool>("includeJumpingForwardDoublets", false);
0216   desc.add<double>("z0Cut", 10.0);
0217   desc.add<double>("ptCut", 0.0);
0218
0219   edm::ParameterSetDescription trackQualityCuts;
0220   trackQualityCuts.add<double>("chi2MaxPt", 10.)->setComment("max pT used to determine the pT-dependent chi2 cut");
0221   trackQualityCuts.add<std::vector<double>>("chi2Coeff", {0.9, 1.8})->setComment("chi2 at 1GeV and at ptMax above");
0222   trackQualityCuts.add<double>("chi2Scale", 8.)
0223       ->setComment(
0224           "Factor to multiply the pT-dependent chi2 cut (currently: 8 for the broken line fit, ?? for the Riemann "
0225           "fit)");
0226   trackQualityCuts.add<double>("tripletMinPt", 0.0)->setComment("Min pT for triplets, in GeV");
0227   trackQualityCuts.add<double>("tripletMaxTip", 0.1)->setComment("Max |Tip| for triplets, in cm");
0228   trackQualityCuts.add<double>("tripletMaxZip", 6.)->setComment("Max |Zip| for triplets, in cm");
0229   trackQualityCuts.add<double>("quadrupletMinPt", 0.0)->setComment("Min pT for quadruplets, in GeV");
0230   trackQualityCuts.add<double>("quadrupletMaxTip", 0.5)->setComment("Max |Tip| for quadruplets, in cm");
0231   trackQualityCuts.add<double>("quadrupletMaxZip", 6.)->setComment("Max |Zip| for quadruplets, in cm");
0232
0233   desc.add<std::vector<int>>(
0234           "phiCuts", std::vector<int>(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts)))
0235       ->setComment("Cuts in phi for cells");
0236
0237   desc.add<edm::ParameterSetDescription>("trackQualityCuts", trackQualityCuts)
0238       ->setComment(
0239           "Quality cuts based on the results of the track fit:\n  - apply a pT-dependent chi2 cut;\n  - apply \"region "
0240           "cuts\" based on the fit results (pT, Tip, Zip).");
0241 }
0242
0243 template <>
0244 void CAHitNtupletGeneratorOnGPU<pixelTopology::Phase2>::fillDescriptions(edm::ParameterSetDescription& desc) {
0245   fillDescriptionsCommon(desc);
0246
0247   desc.add<bool>("idealConditions", false);
0248   desc.add<bool>("includeFarForwards", true);
0249   desc.add<bool>("includeJumpingForwardDoublets", true);
0250   desc.add<double>("z0Cut", 7.5);
0251   desc.add<double>("ptCut", 0.85);
0252
0253   edm::ParameterSetDescription trackQualityCuts;
0254   trackQualityCuts.add<double>("maxChi2", 5.)->setComment("Max normalized chi2");
0255   trackQualityCuts.add<double>("minPt", 0.5)->setComment("Min pT in GeV");
0256   trackQualityCuts.add<double>("maxTip", 0.3)->setComment("Max |Tip| in cm");
0257   trackQualityCuts.add<double>("maxZip", 12.)->setComment("Max |Zip|, in cm");
0258
0259   desc.add<std::vector<int>>(
0260           "phiCuts", std::vector<int>(std::begin(phase2PixelTopology::phicuts), std::end(phase2PixelTopology::phicuts)))
0261       ->setComment("Cuts in phi for cells");
0262
0263   desc.add<edm::ParameterSetDescription>("trackQualityCuts", trackQualityCuts)
0264       ->setComment(
0265           "Quality cuts based on the results of the track fit:\n  - apply cuts based on the fit results (pT, Tip, "
0266           "Zip).");
0267 }
0268
0269 template <typename TrackerTraits>
0270 void CAHitNtupletGeneratorOnGPU<TrackerTraits>::fillDescriptionsCommon(edm::ParameterSetDescription& desc) {
0271   // 87 cm/GeV = 1/(3.8T * 0.3)
0272   // take less than radius given by the hardPtCut and reject everything below
0273   // auto hardCurvCut = 1.f/(0.35 * 87.f);
0274   desc.add<double>("ptmin", 0.9)->setComment("Cut on minimum pt");
0275   desc.add<double>("CAThetaCutBarrel", 0.002)->setComment("Cut on RZ alignement for Barrel");
0276   desc.add<double>("CAThetaCutForward", 0.003)->setComment("Cut on RZ alignment for Forward");
0277   desc.add<double>("hardCurvCut", 1. / (0.35 * 87.))->setComment("Cut on minimum curvature");
0278   desc.add<double>("dcaCutInnerTriplet", 0.15)->setComment("Cut on origin radius when the inner hit is on BPix1");
0279   desc.add<double>("dcaCutOuterTriplet", 0.25)->setComment("Cut on origin radius when the outer hit is on BPix1");
0280   desc.add<bool>("earlyFishbone", true);
0281   desc.add<bool>("lateFishbone", false);
0282   desc.add<bool>("fillStatistics", false);
0283   desc.add<unsigned int>("minHitsPerNtuplet", 4);
0284   desc.add<unsigned int>("maxNumberOfDoublets", TrackerTraits::maxNumberOfDoublets);
0285   desc.add<unsigned int>("minHitsForSharingCut", 10)
0286       ->setComment("Maximum number of hits in a tuple to clean also if the shared hit is on bpx1");
0287
0288   desc.add<bool>("fitNas4", false)->setComment("fit only 4 hits out of N");
0289   desc.add<bool>("doClusterCut", true);
0290   desc.add<bool>("doZ0Cut", true);
0291   desc.add<bool>("doPtCut", true);
0292   desc.add<bool>("useRiemannFit", false)->setComment("true for Riemann, false for BrokenLine");
0293   desc.add<bool>("doSharedHitCut", true)->setComment("Sharing hit nTuples cleaning");
0294   desc.add<bool>("dupPassThrough", false)->setComment("Do not reject duplicate");
0295   desc.add<bool>("useSimpleTripletCleaner", true)->setComment("use alternate implementation");
0296 }
0297
0298 template <typename TrackerTraits>
0299 void CAHitNtupletGeneratorOnGPU<TrackerTraits>::beginJob() {
0300   if (m_params.onGPU_) {
0301     // allocate pinned host memory only if CUDA is available
0302     edm::Service<CUDAInterface> cuda;
0303     if (cuda and cuda->enabled()) {
0304       cudaCheck(cudaMalloc(&m_counters, sizeof(Counters)));
0305       cudaCheck(cudaMemset(m_counters, 0, sizeof(Counters)));
0306     }
0307   } else {
0308     m_counters = new Counters();
0309     memset(m_counters, 0, sizeof(Counters));
0310   }
0311 }
0312
0313 template <typename TrackerTraits>
0314 void CAHitNtupletGeneratorOnGPU<TrackerTraits>::endJob() {
0315   if (m_params.onGPU_) {
0316     // print the gpu statistics and free pinned host memory only if CUDA is available
0317     edm::Service<CUDAInterface> cuda;
0318     if (cuda and cuda->enabled()) {
0319       if (m_params.doStats_) {
0320         // crash on multi-gpu processes
0321         CAHitNtupletGeneratorKernelsGPU<TrackerTraits>::printCounters(m_counters);
0322       }
0323       cudaFree(m_counters);
0324     }
0325   } else {
0326     if (m_params.doStats_) {
0327       CAHitNtupletGeneratorKernelsCPU<TrackerTraits>::printCounters(m_counters);
0328     }
0329     delete m_counters;
0330   }
0331 }
0332
0333 template <typename TrackerTraits>
0334 TrackSoAHeterogeneousDevice<TrackerTraits> CAHitNtupletGeneratorOnGPU<TrackerTraits>::makeTuplesAsync(
0335     HitsOnDevice const& hits_d, float bfield, cudaStream_t stream) const {
0336   using HelixFitOnGPU = HelixFitOnGPU<TrackerTraits>;
0337   using TrackSoA = TrackSoAHeterogeneousDevice<TrackerTraits>;
0338   using GPUKernels = CAHitNtupletGeneratorKernelsGPU<TrackerTraits>;
0339
0340   TrackSoA tracks(stream);
0341
0342   GPUKernels kernels(m_params);
0343   kernels.setCounters(m_counters);
0344   kernels.allocateOnGPU(hits_d.nHits(), stream);
0345
0346   kernels.buildDoublets(hits_d.view(), hits_d.offsetBPIX2(), stream);
0347
0348   kernels.launchKernels(hits_d.view(), tracks.view(), stream);
0349
0350   HelixFitOnGPU fitter(bfield, m_params.fitNas4_);
0351   fitter.allocateOnGPU(kernels.tupleMultiplicity(), tracks.view());
0352   if (m_params.useRiemannFit_) {
0353     fitter.launchRiemannKernels(hits_d.view(), hits_d.nHits(), TrackerTraits::maxNumberOfQuadruplets, stream);
0354   } else {
0355     fitter.launchBrokenLineKernels(hits_d.view(), hits_d.nHits(), TrackerTraits::maxNumberOfQuadruplets, stream);
0356   }
0357   kernels.classifyTuples(hits_d.view(), tracks.view(), stream);
0358 #ifdef GPU_DEBUG
0359   cudaDeviceSynchronize();
0360   cudaCheck(cudaGetLastError());
0361   std::cout << "finished building pixel tracks on GPU" << std::endl;
0362 #endif
0363
0364   return tracks;
0365 }
0366
0367 template <typename TrackerTraits>
0368 TrackSoAHeterogeneousHost<TrackerTraits> CAHitNtupletGeneratorOnGPU<TrackerTraits>::makeTuples(HitsOnHost const& hits_h,
0369                                                                                                float bfield) const {
0370   using HelixFitOnGPU = HelixFitOnGPU<TrackerTraits>;
0371   using TrackSoA = TrackSoAHeterogeneousHost<TrackerTraits>;
0372   using CPUKernels = CAHitNtupletGeneratorKernelsCPU<TrackerTraits>;
0373
0374   TrackSoA tracks;
0375
0376   CPUKernels kernels(m_params);
0377   kernels.setCounters(m_counters);
0378   kernels.allocateOnGPU(hits_h.nHits(), nullptr);
0379
0380   kernels.buildDoublets(hits_h.view(), hits_h.offsetBPIX2(), nullptr);
0381   kernels.launchKernels(hits_h.view(), tracks.view(), nullptr);
0382
0383   if (0 == hits_h.nHits())
0384     return tracks;
0385
0386   // now fit
0387   HelixFitOnGPU fitter(bfield, m_params.fitNas4_);
0388   fitter.allocateOnGPU(kernels.tupleMultiplicity(), tracks.view());
0389
0390   if (m_params.useRiemannFit_) {
0391     fitter.launchRiemannKernelsOnCPU(hits_h.view(), hits_h.nHits(), TrackerTraits::maxNumberOfQuadruplets);
0392   } else {
0393     fitter.launchBrokenLineKernelsOnCPU(hits_h.view(), hits_h.nHits(), TrackerTraits::maxNumberOfQuadruplets);
0394   }
0395
0396   kernels.classifyTuples(hits_h.view(), tracks.view(), nullptr);
0397
0398 #ifdef GPU_DEBUG
0399   std::cout << "finished building pixel tracks on CPU" << std::endl;
0400 #endif
0401
0402   // check that the fixed-size SoA does not overflow
0403   auto maxTracks = tracks.view().metadata().size();
0404   auto nTracks = tracks.view().nTracks();
0405   assert(nTracks < maxTracks);
0406   if (nTracks == maxTracks - 1) {
0407     edm::LogWarning("PixelTracks") << "Unsorted reconstructed pixel tracks truncated to " << maxTracks - 1
0408                                    << " candidates";
0409   }
0410
0411   return tracks;
0412 }
0413
0414 template class CAHitNtupletGeneratorOnGPU<pixelTopology::Phase1>;
0415 template class CAHitNtupletGeneratorOnGPU<pixelTopology::Phase2>;
0416 template class CAHitNtupletGeneratorOnGPU<pixelTopology::HIonPhase1>;