Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-05 22:39:24

0001 //#define GPU_DEBUG
0002 //#define DUMP_GPU_TK_TUPLES
0003 
0004 #include <array>
0005 #include <cassert>
0006 #include <functional>
0007 #include <vector>
0008 
0009 #include <alpaka/alpaka.hpp>
0010 
0011 #include "DataFormats/TrackSoA/interface/TracksDevice.h"
0012 #include "DataFormats/TrackSoA/interface/TracksHost.h"
0013 #include "DataFormats/TrackSoA/interface/alpaka/TracksSoACollection.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "FWCore/Utilities/interface/Exception.h"
0017 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0018 
0019 #include "CAHitNtupletGenerator.h"
0020 #include "CAHitNtupletGeneratorKernels.h"
0021 #include "CAPixelDoublets.h"
0022 #include "CAPixelDoubletsAlgos.h"
0023 
0024 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0025   namespace {
0026 
0027     using namespace caHitNtupletGenerator;
0028     using namespace caPixelDoublets;
0029     using namespace pixelTopology;
0030     using namespace pixelTrack;
0031 
0032     template <typename T>
0033     T sqr(T x) {
0034       return x * x;
0035     }
0036 
0037     //Common Params
0038     void fillDescriptionsCommon(edm::ParameterSetDescription& desc) {
0039       // 87 cm/GeV = 1/(3.8T * 0.3)
0040       // take less than radius given by the hardPtCut and reject everything below
0041       // auto hardCurvCut = 1.f/(0.35 * 87.f);
0042       desc.add<double>("ptmin", 0.9f)->setComment("Cut on minimum pt");
0043       desc.add<double>("CAThetaCutBarrel", 0.002f)->setComment("Cut on RZ alignement for Barrel");
0044       desc.add<double>("CAThetaCutForward", 0.003f)->setComment("Cut on RZ alignment for Forward");
0045       desc.add<double>("hardCurvCut", 1.f / (0.35 * 87.f))
0046           ->setComment("Cut on minimum curvature, used in DCA ntuplet selection");
0047       desc.add<double>("dcaCutInnerTriplet", 0.15f)->setComment("Cut on origin radius when the inner hit is on BPix1");
0048       desc.add<double>("dcaCutOuterTriplet", 0.25f)->setComment("Cut on origin radius when the outer hit is on BPix1");
0049       desc.add<bool>("earlyFishbone", true);
0050       desc.add<bool>("lateFishbone", false);
0051       desc.add<bool>("fillStatistics", false);
0052       desc.add<unsigned int>("minHitsPerNtuplet", 4);
0053       desc.add<unsigned int>("minHitsForSharingCut", 10)
0054           ->setComment("Maximum number of hits in a tuple to clean also if the shared hit is on bpx1");
0055 
0056       desc.add<bool>("fitNas4", false)->setComment("fit only 4 hits out of N");
0057       desc.add<bool>("doClusterCut", true);
0058       desc.add<bool>("doZ0Cut", true);
0059       desc.add<bool>("doPtCut", true);
0060       desc.add<bool>("useRiemannFit", false)->setComment("true for Riemann, false for BrokenLine");
0061       desc.add<bool>("doSharedHitCut", true)->setComment("Sharing hit nTuples cleaning");
0062       desc.add<bool>("dupPassThrough", false)->setComment("Do not reject duplicate");
0063       desc.add<bool>("useSimpleTripletCleaner", true)->setComment("use alternate implementation");
0064     }
0065 
0066     AlgoParams makeCommonParams(edm::ParameterSet const& cfg) {
0067       return AlgoParams({cfg.getParameter<unsigned int>("minHitsForSharingCut"),
0068                          cfg.getParameter<bool>("useRiemannFit"),
0069                          cfg.getParameter<bool>("fitNas4"),
0070                          cfg.getParameter<bool>("includeJumpingForwardDoublets"),
0071                          cfg.getParameter<bool>("earlyFishbone"),
0072                          cfg.getParameter<bool>("lateFishbone"),
0073                          cfg.getParameter<bool>("fillStatistics"),
0074                          cfg.getParameter<bool>("doSharedHitCut"),
0075                          cfg.getParameter<bool>("dupPassThrough"),
0076                          cfg.getParameter<bool>("useSimpleTripletCleaner")});
0077     }
0078 
0079     //This is needed to have the partial specialization for isPhase1Topology/isPhase2Topology
0080     template <typename TrackerTraits, typename Enable = void>
0081     struct TopologyCuts {};
0082 
0083     template <typename TrackerTraits>
0084     struct TopologyCuts<TrackerTraits, isPhase1Topology<TrackerTraits>> {
0085       static constexpr CAParamsT<TrackerTraits> makeCACuts(edm::ParameterSet const& cfg) {
0086         return CAParamsT<TrackerTraits>{{cfg.getParameter<unsigned int>("maxNumberOfDoublets"),
0087                                          cfg.getParameter<unsigned int>("minHitsPerNtuplet"),
0088                                          (float)cfg.getParameter<double>("ptmin"),
0089                                          (float)cfg.getParameter<double>("CAThetaCutBarrel"),
0090                                          (float)cfg.getParameter<double>("CAThetaCutForward"),
0091                                          (float)cfg.getParameter<double>("hardCurvCut"),
0092                                          (float)cfg.getParameter<double>("dcaCutInnerTriplet"),
0093                                          (float)cfg.getParameter<double>("dcaCutOuterTriplet")}};
0094       };
0095 
0096       static constexpr ::pixelTrack::QualityCutsT<TrackerTraits> makeQualityCuts(edm::ParameterSet const& pset) {
0097         auto coeff = pset.getParameter<std::array<double, 2>>("chi2Coeff");
0098         auto ptMax = pset.getParameter<double>("chi2MaxPt");
0099 
0100         coeff[1] = (coeff[1] - coeff[0]) / log2(ptMax);
0101         return ::pixelTrack::QualityCutsT<TrackerTraits>{// polynomial coefficients for the pT-dependent chi2 cut
0102                                                          {(float)coeff[0], (float)coeff[1], 0.f, 0.f},
0103                                                          // max pT used to determine the chi2 cut
0104                                                          (float)ptMax,
0105                                                          // chi2 scale factor: 8 for broken line fit, ?? for Riemann fit
0106                                                          (float)pset.getParameter<double>("chi2Scale"),
0107                                                          // regional cuts for triplets
0108                                                          {(float)pset.getParameter<double>("tripletMaxTip"),
0109                                                           (float)pset.getParameter<double>("tripletMinPt"),
0110                                                           (float)pset.getParameter<double>("tripletMaxZip")},
0111                                                          // regional cuts for quadruplets
0112                                                          {(float)pset.getParameter<double>("quadrupletMaxTip"),
0113                                                           (float)pset.getParameter<double>("quadrupletMinPt"),
0114                                                           (float)pset.getParameter<double>("quadrupletMaxZip")}};
0115       }
0116     };
0117 
0118     template <typename TrackerTraits>
0119     struct TopologyCuts<TrackerTraits, isPhase2Topology<TrackerTraits>> {
0120       static constexpr CAParamsT<TrackerTraits> makeCACuts(edm::ParameterSet const& cfg) {
0121         return CAParamsT<TrackerTraits>{{cfg.getParameter<unsigned int>("maxNumberOfDoublets"),
0122                                          cfg.getParameter<unsigned int>("minHitsPerNtuplet"),
0123                                          (float)cfg.getParameter<double>("ptmin"),
0124                                          (float)cfg.getParameter<double>("CAThetaCutBarrel"),
0125                                          (float)cfg.getParameter<double>("CAThetaCutForward"),
0126                                          (float)cfg.getParameter<double>("hardCurvCut"),
0127                                          (float)cfg.getParameter<double>("dcaCutInnerTriplet"),
0128                                          (float)cfg.getParameter<double>("dcaCutOuterTriplet")},
0129                                         {(bool)cfg.getParameter<bool>("includeFarForwards")}};
0130       }
0131 
0132       static constexpr ::pixelTrack::QualityCutsT<TrackerTraits> makeQualityCuts(edm::ParameterSet const& pset) {
0133         return ::pixelTrack::QualityCutsT<TrackerTraits>{
0134             static_cast<float>(pset.getParameter<double>("maxChi2")),
0135             static_cast<float>(pset.getParameter<double>("minPt")),
0136             static_cast<float>(pset.getParameter<double>("maxTip")),
0137             static_cast<float>(pset.getParameter<double>("maxZip")),
0138         };
0139       }
0140     };
0141 
0142     //Cell Cuts, as they are the cuts have the same logic for Phase2 and Phase1
0143     //keeping them separate would allow further differentiation in the future
0144     //moving them to TopologyCuts and using the same syntax
0145     template <typename TrackerTraits>
0146     CellCutsT<TrackerTraits> makeCellCuts(edm::ParameterSet const& cfg) {
0147       return CellCutsT<TrackerTraits>{cfg.getParameter<bool>("doClusterCut"),
0148                                       cfg.getParameter<bool>("doZ0Cut"),
0149                                       cfg.getParameter<bool>("doPtCut"),
0150                                       cfg.getParameter<bool>("idealConditions"),
0151                                       (float)cfg.getParameter<double>("cellZ0Cut"),
0152                                       (float)cfg.getParameter<double>("cellPtCut"),
0153                                       cfg.getParameter<std::vector<int>>("phiCuts")};
0154     }
0155 
0156   }  // namespace
0157 
0158   using namespace std;
0159 
0160   template <typename TrackerTraits>
0161   CAHitNtupletGenerator<TrackerTraits>::CAHitNtupletGenerator(const edm::ParameterSet& cfg)
0162       : m_params(makeCommonParams(cfg),
0163                  makeCellCuts<TrackerTraits>(cfg),
0164                  TopologyCuts<TrackerTraits>::makeQualityCuts(cfg.getParameterSet("trackQualityCuts")),
0165                  TopologyCuts<TrackerTraits>::makeCACuts(cfg)) {
0166 #ifdef DUMP_GPU_TK_TUPLES
0167     printf("TK: %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s\n",
0168            "tid",
0169            "qual",
0170            "nh",
0171            "nl",
0172            "charge",
0173            "pt",
0174            "eta",
0175            "phi",
0176            "tip",
0177            "zip",
0178            "chi2",
0179            "h1",
0180            "h2",
0181            "h3",
0182            "h4",
0183            "h5",
0184            "hn");
0185 #endif
0186   }
0187 
0188   template <typename TrackerTraits>
0189   void CAHitNtupletGenerator<TrackerTraits>::fillPSetDescription(edm::ParameterSetDescription& desc) {
0190     static_assert(sizeof(TrackerTraits) == 0,
0191                   "Note: this fillPSetDescription is a dummy one. Please specialise it for the correct version of "
0192                   "CAHitNtupletGenerator<TrackerTraits>.");
0193   }
0194 
0195   template <>
0196   void CAHitNtupletGenerator<pixelTopology::Phase1>::fillPSetDescription(edm::ParameterSetDescription& desc) {
0197     fillDescriptionsCommon(desc);
0198 
0199     desc.add<unsigned int>("maxNumberOfDoublets", pixelTopology::Phase1::maxNumberOfDoublets);
0200     desc.add<bool>("idealConditions", true);
0201     desc.add<bool>("includeJumpingForwardDoublets", false);
0202     desc.add<double>("cellZ0Cut", 12.0);
0203     desc.add<double>("cellPtCut", 0.5);
0204 
0205     edm::ParameterSetDescription trackQualityCuts;
0206     trackQualityCuts.add<double>("chi2MaxPt", 10.)->setComment("max pT used to determine the pT-dependent chi2 cut");
0207     trackQualityCuts.add<std::vector<double>>("chi2Coeff", {0.9, 1.8})->setComment("chi2 at 1GeV and at ptMax above");
0208     trackQualityCuts.add<double>("chi2Scale", 8.)
0209         ->setComment(
0210             "Factor to multiply the pT-dependent chi2 cut (currently: 8 for the broken line fit, ?? for the Riemann "
0211             "fit)");
0212     trackQualityCuts.add<double>("tripletMinPt", 0.5)->setComment("Min pT for triplets, in GeV");
0213     trackQualityCuts.add<double>("tripletMaxTip", 0.3)->setComment("Max |Tip| for triplets, in cm");
0214     trackQualityCuts.add<double>("tripletMaxZip", 12.)->setComment("Max |Zip| for triplets, in cm");
0215     trackQualityCuts.add<double>("quadrupletMinPt", 0.3)->setComment("Min pT for quadruplets, in GeV");
0216     trackQualityCuts.add<double>("quadrupletMaxTip", 0.5)->setComment("Max |Tip| for quadruplets, in cm");
0217     trackQualityCuts.add<double>("quadrupletMaxZip", 12.)->setComment("Max |Zip| for quadruplets, in cm");
0218     desc.add<edm::ParameterSetDescription>("trackQualityCuts", trackQualityCuts)
0219         ->setComment(
0220             "Quality cuts based on the results of the track fit:\n  - apply a pT-dependent chi2 cut;\n  - apply "
0221             "\"region "
0222             "cuts\" based on the fit results (pT, Tip, Zip).");
0223 
0224     desc.add<std::vector<int>>(
0225             "phiCuts",
0226             std::vector<int>(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts)))
0227         ->setComment("Cuts in phi for cells");
0228   }
0229 
0230   template <>
0231   void CAHitNtupletGenerator<pixelTopology::HIonPhase1>::fillPSetDescription(edm::ParameterSetDescription& desc) {
0232     fillDescriptionsCommon(desc);
0233 
0234     desc.add<unsigned int>("maxNumberOfDoublets", pixelTopology::HIonPhase1::maxNumberOfDoublets);
0235     desc.add<bool>("idealConditions", false);
0236     desc.add<bool>("includeJumpingForwardDoublets", false);
0237     desc.add<double>("cellZ0Cut", 10.0);
0238     desc.add<double>("cellPtCut", 0.0);
0239 
0240     edm::ParameterSetDescription trackQualityCuts;
0241     trackQualityCuts.add<double>("chi2MaxPt", 10.)->setComment("max pT used to determine the pT-dependent chi2 cut");
0242     trackQualityCuts.add<std::vector<double>>("chi2Coeff", {0.9, 1.8})->setComment("chi2 at 1GeV and at ptMax above");
0243     trackQualityCuts.add<double>("chi2Scale", 8.)
0244         ->setComment(
0245             "Factor to multiply the pT-dependent chi2 cut (currently: 8 for the broken line fit, ?? for the Riemann "
0246             "fit)");
0247     trackQualityCuts.add<double>("tripletMinPt", 0.0)->setComment("Min pT for triplets, in GeV");
0248     trackQualityCuts.add<double>("tripletMaxTip", 0.1)->setComment("Max |Tip| for triplets, in cm");
0249     trackQualityCuts.add<double>("tripletMaxZip", 6.)->setComment("Max |Zip| for triplets, in cm");
0250     trackQualityCuts.add<double>("quadrupletMinPt", 0.0)->setComment("Min pT for quadruplets, in GeV");
0251     trackQualityCuts.add<double>("quadrupletMaxTip", 0.5)->setComment("Max |Tip| for quadruplets, in cm");
0252     trackQualityCuts.add<double>("quadrupletMaxZip", 6.)->setComment("Max |Zip| for quadruplets, in cm");
0253 
0254     desc.add<edm::ParameterSetDescription>("trackQualityCuts", trackQualityCuts)
0255         ->setComment(
0256             "Quality cuts based on the results of the track fit:\n  - apply a pT-dependent chi2 cut;\n  - apply "
0257             "\"region "
0258             "cuts\" based on the fit results (pT, Tip, Zip).");
0259 
0260     desc.add<std::vector<int>>(
0261             "phiCuts",
0262             std::vector<int>(std::begin(phase1PixelTopology::phicuts), std::end(phase1PixelTopology::phicuts)))
0263         ->setComment("Cuts in phi for cells");
0264   }
0265 
0266   template <>
0267   void CAHitNtupletGenerator<pixelTopology::Phase2>::fillPSetDescription(edm::ParameterSetDescription& desc) {
0268     fillDescriptionsCommon(desc);
0269 
0270     desc.add<unsigned int>("maxNumberOfDoublets", pixelTopology::Phase2::maxNumberOfDoublets);
0271     desc.add<bool>("idealConditions", false);
0272     desc.add<bool>("includeFarForwards", true);
0273     desc.add<bool>("includeJumpingForwardDoublets", true);
0274     desc.add<double>("cellZ0Cut", 7.5);
0275     desc.add<double>("cellPtCut", 0.85);
0276 
0277     edm::ParameterSetDescription trackQualityCuts;
0278     trackQualityCuts.add<double>("maxChi2", 5.)->setComment("Max normalized chi2");
0279     trackQualityCuts.add<double>("minPt", 0.5)->setComment("Min pT in GeV");
0280     trackQualityCuts.add<double>("maxTip", 0.3)->setComment("Max |Tip| in cm");
0281     trackQualityCuts.add<double>("maxZip", 12.)->setComment("Max |Zip|, in cm");
0282     desc.add<edm::ParameterSetDescription>("trackQualityCuts", trackQualityCuts)
0283         ->setComment(
0284             "Quality cuts based on the results of the track fit:\n  - apply cuts based on the fit results (pT, Tip, "
0285             "Zip).");
0286 
0287     desc.add<std::vector<int>>(
0288             "phiCuts",
0289             std::vector<int>(std::begin(phase2PixelTopology::phicuts), std::end(phase2PixelTopology::phicuts)))
0290         ->setComment("Cuts in phi for cells");
0291   }
0292 
0293   template <typename TrackerTraits>
0294   TracksSoACollection<TrackerTraits> CAHitNtupletGenerator<TrackerTraits>::makeTuplesAsync(
0295       HitsOnDevice const& hits_d, ParamsOnDevice const* cpeParams, float bfield, Queue& queue) const {
0296     using HelixFit = HelixFit<TrackerTraits>;
0297     using TrackSoA = TracksSoACollection<TrackerTraits>;
0298     using GPUKernels = CAHitNtupletGeneratorKernels<TrackerTraits>;
0299 
0300     TrackSoA tracks(queue);
0301 
0302     // Don't bother if less than 2 this
0303     if (hits_d.view().metadata().size() < 2) {
0304       const auto device = alpaka::getDev(queue);
0305       auto ntracks_d = cms::alpakatools::make_device_view(device, tracks.view().nTracks());
0306       alpaka::memset(queue, ntracks_d, 0);
0307       return tracks;
0308     }
0309     GPUKernels kernels(m_params, hits_d.view().metadata().size(), hits_d.offsetBPIX2(), queue);
0310 
0311     kernels.buildDoublets(hits_d.view(), hits_d.offsetBPIX2(), queue);
0312     kernels.launchKernels(hits_d.view(), hits_d.offsetBPIX2(), tracks.view(), queue);
0313 
0314     HelixFit fitter(bfield, m_params.fitNas4_);
0315     fitter.allocate(kernels.tupleMultiplicity(), tracks.view());
0316     if (m_params.useRiemannFit_) {
0317       fitter.launchRiemannKernels(
0318           hits_d.view(), cpeParams, hits_d.view().metadata().size(), TrackerTraits::maxNumberOfQuadruplets, queue);
0319     } else {
0320       fitter.launchBrokenLineKernels(
0321           hits_d.view(), cpeParams, hits_d.view().metadata().size(), TrackerTraits::maxNumberOfQuadruplets, queue);
0322     }
0323     kernels.classifyTuples(hits_d.view(), tracks.view(), queue);
0324 #ifdef GPU_DEBUG
0325     alpaka::wait(queue);
0326     std::cout << "finished building pixel tracks on GPU" << std::endl;
0327 #endif
0328 
0329     return tracks;
0330   }
0331 
0332   template class CAHitNtupletGenerator<pixelTopology::Phase1>;
0333   template class CAHitNtupletGenerator<pixelTopology::Phase2>;
0334   template class CAHitNtupletGenerator<pixelTopology::HIonPhase1>;
0335 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE