diff --git a/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h b/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h index f74717c41e4d7..1c076e8705aa0 100644 --- a/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h +++ b/CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h @@ -1,14 +1,23 @@ #ifndef CUDADataFormats_Track_TrackHeterogeneousT_H #define CUDADataFormats_Track_TrackHeterogeneousT_H +#include +#include + #include "CUDADataFormats/Track/interface/TrajectoryStateSoAT.h" #include "HeterogeneousCore/CUDAUtilities/interface/HistoContainer.h" #include "CUDADataFormats/Common/interface/HeterogeneousSoA.h" namespace pixelTrack { - enum class Quality : uint8_t { bad = 0, dup, loose, strict, tight, highPurity }; -} + enum class Quality : uint8_t { bad = 0, edup, dup, loose, strict, tight, highPurity, notQuality }; + constexpr uint32_t qualitySize{uint8_t(Quality::notQuality)}; + const std::string qualityName[qualitySize]{"bad", "edup", "dup", "loose", "strict", "tight", "highPurity"}; + inline Quality qualityByName(std::string const &name) { + auto qp = std::find(qualityName, qualityName + qualitySize, name) - qualityName; + return static_cast(qp); + } +} // namespace pixelTrack template class TrackSoAHeterogeneousT { diff --git a/CUDADataFormats/Track/test/BuildFile.xml b/CUDADataFormats/Track/test/BuildFile.xml index 598b345d4709d..a16afd104442a 100644 --- a/CUDADataFormats/Track/test/BuildFile.xml +++ b/CUDADataFormats/Track/test/BuildFile.xml @@ -1,5 +1,10 @@ + + + + + diff --git a/CUDADataFormats/Track/test/TrackSoAHeterogeneous_t.cpp b/CUDADataFormats/Track/test/TrackSoAHeterogeneous_t.cpp new file mode 100644 index 0000000000000..9708b689dd05b --- /dev/null +++ b/CUDADataFormats/Track/test/TrackSoAHeterogeneous_t.cpp @@ -0,0 +1,21 @@ +#include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h" + +#include +#include + +int main() { + // test quality + + auto q = pixelTrack::qualityByName("tight"); + assert(pixelTrack::Quality::tight == q); + q = pixelTrack::qualityByName("toght"); + assert(pixelTrack::Quality::notQuality == q); + + for (uint32_t i = 0; i < pixelTrack::qualitySize; ++i) { + auto const qt = static_cast(i); + auto q = pixelTrack::qualityByName(pixelTrack::qualityName[i]); + assert(qt == q); + } + + return 0; +} diff --git a/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py b/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py index d5deba78b46c8..7c3519ed8786b 100644 --- a/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py +++ b/DQM/TrackingMonitorSource/python/pixelTracksMonitoring_cff.py @@ -8,7 +8,7 @@ pixelTracksMonitor.beamSpot = 'offlineBeamSpot' pixelTracksMonitor.primaryVertex = 'pixelVertices' pixelTracksMonitor.pvNDOF = 1 -pixelTracksMonitor.doAllPlots = True +pixelTracksMonitor.doAllPlots = False pixelTracksMonitor.doLumiAnalysis = True pixelTracksMonitor.doProfilesVsLS = True pixelTracksMonitor.doDCAPlots = True @@ -26,27 +26,54 @@ cut = cms.string("") ) -pixelTracksPt0to1 = _trackSelector.clone(cut = "pt >= 0 & pt < 1 ") -pixelTracksPt1 = _trackSelector.clone(cut = "pt >= 1 ") -from DQM.TrackingMonitorSource.TrackCollections2monitor_cff import highPurityPV0p1 as _highPurityPV0p1 -pixelTracksPV0p1 = _highPurityPV0p1.clone( - src = "pixelTracks", - quality = "", - vertexTag = "goodPixelVertices" -) +quality = { + "L" : "loose", + "T" : "tight", + "HP" : "highPurity", +} -pixelTracksMonitorPt0to1 = pixelTracksMonitor.clone( - TrackProducer = "pixelTracksPt0to1", - FolderName = "Tracking/PixelTrackParameters/pt_0to1" -) -pixelTracksMonitorPt1 = pixelTracksMonitor.clone( - TrackProducer = "pixelTracksPt1", - FolderName = "Tracking/PixelTrackParameters/pt_1" -) -pixelTracksMonitorPV0p1 = pixelTracksMonitor.clone( - TrackProducer = "pixelTracksPV0p1", - FolderName = "Tracking/PixelTrackParameters/dzPV0p1" -) +for key,value in quality.items(): + label = "pixelTrks"+key +# print label + cutstring = "quality('" + value + "')" +# print cutstring + if label not in globals(): + locals()[label] = _trackSelector.clone( cut = cutstring ) + locals()[label].setLabel(label) + else : + print(label,"already configured") + +for key,value in quality.items(): + label = "pixelTrksMonitor"+key + locals()[label] = pixelTracksMonitor.clone( + TrackProducer = "pixelTrks"+key, + FolderName = "Tracking/PixelTrackParameters/"+value + ) + locals()[label].setLabel(label) + +ntuplet = { + '3' : "3Hits", # ==3 + '4' : "4Hits" # >=4 +} +for kN,vN in ntuplet.items(): + for key,value in quality.items(): + label = "pixelTrks" + vN + key +# print label + + cutstring = "numberOfValidHits == " + kN + " & quality('" + value + "')" +# print cutstring + locals()[label] = _trackSelector.clone( cut = cutstring ) + locals()[label].setLabel(label) + +for kN,vN in ntuplet.items(): + for key,value in quality.items(): + label = "pixelTrks" + vN + "Monitor" + key +# print label + locals()[label] = pixelTracksMonitor.clone( + TrackProducer = "pixelTrks" + vN + key, + FolderName = "Tracking/PixelTrackParameters/" + vN + "/" + value + ) + locals()[label].setLabel(label) from CommonTools.ParticleFlow.goodOfflinePrimaryVertices_cfi import goodOfflinePrimaryVertices as _goodOfflinePrimaryVertices @@ -62,16 +89,27 @@ pixelTracksMonitoringTask = cms.Task( goodPixelVertices, - pixelTracksPt0to1, - pixelTracksPt1, - pixelTracksPV0p1, ) +for category in ["pixelTrks", "pixelTrks3Hits", "pixelTrks4Hits"]: + for key in quality: + label = category+key +# print label + pixelTracksMonitoringTask.add(locals()[label]) + +allPixelTracksMonitoring = cms.Sequence() +for category in ["pixelTrksMonitor", "pixelTrks3HitsMonitor", "pixelTrks4HitsMonitor" ]: + for key in quality: + label = category+key +# print label + allPixelTracksMonitoring += locals()[label] + pixelTracksMonitoring = cms.Sequence( - pixelTracksMonitor + - pixelTracksMonitorPt0to1 + - pixelTracksMonitorPt1 + - pixelTracksMonitorPV0p1 + + allPixelTracksMonitoring + pixelVertexResolution, pixelTracksMonitoringTask ) + + + + diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc index 91c3a44cc8643..eb25abc2d32ab 100644 --- a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc +++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.cc @@ -15,38 +15,60 @@ #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" -#include "PixelTrackProducer.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackReconstruction.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "storeTracks.h" +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +namespace edm { + class Event; + class EventSetup; + class ParameterSet; + class ConfigurationDescriptions; +} // namespace edm +class TrackerTopology; + using namespace pixeltrackfitting; using edm::ParameterSet; -PixelTrackProducer::PixelTrackProducer(const ParameterSet& cfg) : theReconstruction(cfg, consumesCollector()) { - edm::LogInfo("PixelTrackProducer") << " construction..."; - produces(); - produces(); - produces(); -} +class PixelTrackProducer : public edm::stream::EDProducer<> { +public: + explicit PixelTrackProducer(const edm::ParameterSet& cfg) + : theReconstruction(cfg, consumesCollector()), htTopoToken_(esConsumes()) { + edm::LogInfo("PixelTrackProducer") << " construction..."; + produces(); + produces(); + produces(); + } + + ~PixelTrackProducer() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; -PixelTrackProducer::~PixelTrackProducer() {} + desc.add("passLabel", "pixelTracks"); // What is this? It is not used anywhere in this code. + PixelTrackReconstruction::fillDescriptions(desc); -void PixelTrackProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - edm::ParameterSetDescription desc; + descriptions.add("pixelTracks", desc); + } - desc.add("passLabel", "pixelTracks"); // What is this? It is not used anywhere in this code. - PixelTrackReconstruction::fillDescriptions(desc); + void produce(edm::Event& ev, const edm::EventSetup& es) override { + LogDebug("PixelTrackProducer, produce") << "event# :" << ev.id(); - descriptions.add("pixelTracks", desc); -} + TracksWithTTRHs tracks; + theReconstruction.run(tracks, ev, es); + auto htTopo = es.getData(htTopoToken_); -void PixelTrackProducer::produce(edm::Event& ev, const edm::EventSetup& es) { - LogDebug("PixelTrackProducer, produce") << "event# :" << ev.id(); + // store tracks + storeTracks(ev, tracks, htTopo); + } - TracksWithTTRHs tracks; - theReconstruction.run(tracks, ev, es); - edm::ESHandle httopo; - es.get().get(httopo); +private: + PixelTrackReconstruction theReconstruction; + const edm::ESGetToken htTopoToken_; +}; - // store tracks - storeTracks(ev, tracks, *httopo); -} +DEFINE_FWK_MODULE(PixelTrackProducer); diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h deleted file mode 100644 index c38fd44c0d7f5..0000000000000 --- a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducer.h +++ /dev/null @@ -1,29 +0,0 @@ -#ifndef RecoPixelVertexing_PixelTrackFitting_plugins_PixelTrackProducer_h -#define RecoPixelVertexing_PixelTrackFitting_plugins_PixelTrackProducer_h - -#include "FWCore/Framework/interface/stream/EDProducer.h" -#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackReconstruction.h" - -namespace edm { - class Event; - class EventSetup; - class ParameterSet; - class ConfigurationDescriptions; -} // namespace edm -class TrackerTopology; - -class PixelTrackProducer : public edm::stream::EDProducer<> { -public: - explicit PixelTrackProducer(const edm::ParameterSet& conf); - - ~PixelTrackProducer() override; - - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - - void produce(edm::Event& ev, const edm::EventSetup& es) override; - -private: - PixelTrackReconstruction theReconstruction; -}; - -#endif // RecoPixelVertexing_PixelTrackFitting_plugins_PixelTrackProducer_h diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc index 60225eceebc00..cabe23c818470 100644 --- a/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc +++ b/RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc @@ -63,6 +63,7 @@ class PixelTrackProducerFromSoA : public edm::global::EDProducer<> { const edm::ESGetToken ttTopoToken_; int32_t const minNumberOfHits_; + pixelTrack::Quality const minQuality_; }; PixelTrackProducerFromSoA::PixelTrackProducerFromSoA(const edm::ParameterSet &iConfig) @@ -72,7 +73,16 @@ PixelTrackProducerFromSoA::PixelTrackProducerFromSoA(const edm::ParameterSet &iC hmsToken_(consumes(iConfig.getParameter("pixelRecHitLegacySrc"))), idealMagneticFieldToken_(esConsumes()), ttTopoToken_(esConsumes()), - minNumberOfHits_(iConfig.getParameter("minNumberOfHits")) { + minNumberOfHits_(iConfig.getParameter("minNumberOfHits")), + minQuality_(pixelTrack::qualityByName(iConfig.getParameter("minQuality"))) { + if (minQuality_ == pixelTrack::Quality::notQuality) { + throw cms::Exception("PixelTrackConfiguration") + << iConfig.getParameter("minQuality") + " is not a pixelTrack::Quality"; + } + if (minQuality_ < pixelTrack::Quality::dup) { + throw cms::Exception("PixelTrackConfiguration") + << iConfig.getParameter("minQuality") + " not supported"; + } produces(); produces(); produces(); @@ -85,13 +95,23 @@ void PixelTrackProducerFromSoA::fillDescriptions(edm::ConfigurationDescriptions desc.add("trackSrc", edm::InputTag("pixelTracksSoA")); desc.add("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplittingLegacy")); desc.add("minNumberOfHits", 0); - + desc.add("minQuality", "loose"); descriptions.addWithDefaultLabel(desc); } void PixelTrackProducerFromSoA::produce(edm::StreamID streamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const { + // enum class Quality : uint8_t { bad = 0, edup, dup, loose, strict, tight, highPurity }; + reco::TrackBase::TrackQuality recoQuality[] = {reco::TrackBase::undefQuality, + reco::TrackBase::undefQuality, + reco::TrackBase::discarded, + reco::TrackBase::loose, + reco::TrackBase::tight, + reco::TrackBase::tight, + reco::TrackBase::highPurity}; + assert(reco::TrackBase::highPurity == recoQuality[int(pixelTrack::Quality::highPurity)]); + // std::cout << "Converting gpu helix in reco tracks" << std::endl; auto indToEdmP = std::make_unique(); @@ -137,6 +157,8 @@ void PixelTrackProducerFromSoA::produce(edm::StreamID streamID, auto const &hitIndices = tsoa.hitIndices; auto maxTracks = tsoa.stride(); + tracks.reserve(maxTracks); + int32_t nt = 0; for (int32_t it = 0; it < maxTracks; ++it) { @@ -145,7 +167,7 @@ void PixelTrackProducerFromSoA::produce(edm::StreamID streamID, break; // this is a guard: maybe we need to move to nTracks... indToEdm.push_back(-1); auto q = quality[it]; - if (q != pixelTrack::Quality::loose) + if (q < minQuality_) continue; if (nHits < minNumberOfHits_) continue; @@ -192,6 +214,18 @@ void PixelTrackProducerFromSoA::produce(edm::StreamID streamID, math::XYZVector mom(pp.x(), pp.y(), pp.z()); auto track = std::make_unique(chi2, ndof, pos, mom, gp.charge(), CurvilinearTrajectoryError(mo)); + + // bad and edup not supported as fit not present or not reliable + auto tkq = recoQuality[int(q)]; + track->setQuality(tkq); + // loose,tight and HP are inclusive + if (reco::TrackBase::highPurity == tkq) { + track->setQuality(reco::TrackBase::tight); + track->setQuality(reco::TrackBase::loose); + } else if (reco::TrackBase::tight == tkq) { + track->setQuality(reco::TrackBase::loose); + } + track->setQuality(tkq); // filter??? tracks.emplace_back(track.release(), hits); } diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/SealModule.cc b/RecoPixelVertexing/PixelTrackFitting/plugins/SealModule.cc deleted file mode 100644 index 11a4c3e12c308..0000000000000 --- a/RecoPixelVertexing/PixelTrackFitting/plugins/SealModule.cc +++ /dev/null @@ -1,5 +0,0 @@ -#include "FWCore/PluginManager/interface/ModuleDef.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "PixelTrackProducer.h" -DEFINE_FWK_MODULE(PixelTrackProducer); diff --git a/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h b/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h index 59101b6ba5214..c16923a20c9f8 100644 --- a/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h +++ b/RecoPixelVertexing/PixelTrackFitting/plugins/storeTracks.h @@ -24,13 +24,16 @@ void storeTracks(Ev& ev, const TWH& tracksWithHits, const TrackerTopology& ttopo int cc = 0, nTracks = tracksWithHits.size(); + trackExtras->resize(nTracks); + tracks->reserve(nTracks); + recHits->reserve(4 * nTracks); + for (int i = 0; i < nTracks; i++) { reco::Track* track = tracksWithHits[i].first; const auto& hits = tracksWithHits[i].second; for (unsigned int k = 0; k < hits.size(); k++) { - auto* hit = hits[k]->clone(); - + auto* hit = hits[k]->clone(); // need to clone (at least if from SoA) track->appendHitPattern(*hit, ttopo); recHits->push_back(hit); } @@ -44,17 +47,16 @@ void storeTracks(Ev& ev, const TWH& tracksWithHits, const TrackerTopology& ttopo edm::RefProd hitCollProd(ohRH); for (int k = 0; k < nTracks; k++) { - reco::TrackExtra theTrackExtra{}; + auto& aTrackExtra = (*trackExtras)[k]; //fill the TrackExtra with TrackingRecHitRef - unsigned int nHits = tracks->at(k).numberOfValidHits(); - theTrackExtra.setHits(hitCollProd, cc, nHits); + unsigned int nHits = (*tracks)[k].numberOfValidHits(); + aTrackExtra.setHits(hitCollProd, cc, nHits); cc += nHits; AlgebraicVector5 v = AlgebraicVector5(0, 0, 0, 0, 0); reco::TrackExtra::TrajParams trajParams(nHits, LocalTrajectoryParameters(v, 1.)); reco::TrackExtra::Chi2sFive chi2s(nHits, 0); - theTrackExtra.setTrajParams(std::move(trajParams), std::move(chi2s)); - trackExtras->push_back(theTrackExtra); + aTrackExtra.setTrajParams(std::move(trajParams), std::move(chi2s)); } LogDebug("TrackProducer") << "put the collection of TrackExtra in the event" @@ -63,7 +65,7 @@ void storeTracks(Ev& ev, const TWH& tracksWithHits, const TrackerTopology& ttopo for (int k = 0; k < nTracks; k++) { const reco::TrackExtraRef theTrackExtraRef(ohTE, k); - (tracks->at(k)).setExtra(theTrackExtraRef); + (*tracks)[k].setExtra(theTrackExtraRef); } ev.put(std::move(tracks)); diff --git a/RecoPixelVertexing/PixelTrackFitting/src/PixelTrackReconstruction.cc b/RecoPixelVertexing/PixelTrackFitting/src/PixelTrackReconstruction.cc index d55d10fba4265..065fdfe2b25a0 100644 --- a/RecoPixelVertexing/PixelTrackFitting/src/PixelTrackReconstruction.cc +++ b/RecoPixelVertexing/PixelTrackFitting/src/PixelTrackReconstruction.cc @@ -77,7 +77,13 @@ void PixelTrackReconstruction::run(TracksWithTTRHs& tracks, edm::Event& ev, cons continue; } } - + /* roll back to verify if HLT tau is affected + // all legacy tracks are "highPurity" + // setting all others bits as well (as in ckf) + track->setQuality(reco::TrackBase::loose); + track->setQuality(reco::TrackBase::tight); + track->setQuality(reco::TrackBase::highPurity); + */ // add tracks tracks.emplace_back(track.release(), tuplet); } diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cc b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cc index 5a2157bf87b7c..36c8c1986d112 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cc +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cc @@ -125,7 +125,7 @@ void CAHitNtupletGeneratorKernelsCPU::launchKernels(HitsOnCPU const &hh, TkSoA * cms::cuda::finalizeBulk(device_hitTuple_apc_, tuples_d); // remove duplicates (tracks that share a doublet) - kernel_earlyDuplicateRemover(device_theCells_.get(), device_nCells_, tuples_d, quality_d); + kernel_earlyDuplicateRemover(device_theCells_.get(), device_nCells_, tuples_d, quality_d, params_.dupPassThrough_); kernel_countMultiplicity(tuples_d, quality_d, device_tupleMultiplicity_.get()); cms::cuda::launchFinalize(device_tupleMultiplicity_.get(), cudaStream); @@ -166,7 +166,7 @@ void CAHitNtupletGeneratorKernelsCPU::classifyTuples(HitsOnCPU const &hh, TkSoA } // remove duplicates (tracks that share a doublet) - kernel_fastDuplicateRemover(device_theCells_.get(), device_nCells_, tuples_d, tracks_d); + kernel_fastDuplicateRemover(device_theCells_.get(), device_nCells_, tuples_d, tracks_d, params_.dupPassThrough_); // fill hit->track "map" if (params_.doSharedHitCut_ || params_.doStats_) { @@ -175,12 +175,41 @@ void CAHitNtupletGeneratorKernelsCPU::classifyTuples(HitsOnCPU const &hh, TkSoA kernel_fillHitInTracks(tuples_d, quality_d, device_hitToTuple_.get()); } - // remove duplicates (tracks that share a hit) + // remove duplicates (tracks that share at least one hit) if (params_.doSharedHitCut_) { - kernel_sharedHitCleaner( - hh.view(), tuples_d, tracks_d, quality_d, params_.minHitsForSharingCut_, device_hitToTuple_.get()); + kernel_rejectDuplicate(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + + kernel_sharedHitCleaner(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + if (params_.useSimpleTripletCleaner_) { + kernel_simpleTripletCleaner(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + } else { + kernel_tripletCleaner(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + } } - if (params_.doStats_) { // counters (add flag???) kernel_doStatsForHitInTracks(device_hitToTuple_.get(), counters_); diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cu b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cu index ca334223f8fc7..ec11631876f46 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cu +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.cu @@ -103,7 +103,7 @@ void CAHitNtupletGeneratorKernelsGPU::launchKernels(HitsOnCPU const &hh, TkSoA * // remove duplicates (tracks that share a doublet) numberOfBlocks = nDoubletBlocks(blockSize); kernel_earlyDuplicateRemover<<>>( - device_theCells_.get(), device_nCells_, tuples_d, quality_d); + device_theCells_.get(), device_nCells_, tuples_d, quality_d, params_.dupPassThrough_); cudaCheck(cudaGetLastError()); blockSize = 128; @@ -244,16 +244,16 @@ void CAHitNtupletGeneratorKernelsGPU::classifyTuples(HitsOnCPU const &hh, TkSoA cudaCheck(cudaGetLastError()); } - // remove duplicates (tracks that share a doublet) + // mark duplicates (tracks that share a doublet) numberOfBlocks = nDoubletBlocks(blockSize); kernel_fastDuplicateRemover<<>>( - device_theCells_.get(), device_nCells_, tuples_d, tracks_d); + device_theCells_.get(), device_nCells_, tuples_d, tracks_d, params_.dupPassThrough_); cudaCheck(cudaGetLastError()); #ifdef GPU_DEBUG cudaCheck(cudaDeviceSynchronize()); #endif - if (params_.minHitsPerNtuplet_ < 4 || params_.doStats_) { + if (params_.doSharedHitCut_ || params_.doStats_) { // fill hit->track "map" assert(hitToTupleView_.offSize > nhits); numberOfBlocks = nQuadrupletBlocks(blockSize); @@ -272,11 +272,46 @@ void CAHitNtupletGeneratorKernelsGPU::classifyTuples(HitsOnCPU const &hh, TkSoA } if (params_.doSharedHitCut_) { - // remove duplicates (tracks that share a hit) + // mark duplicates (tracks that share at least one hit) numberOfBlocks = (hitToTupleView_.offSize + blockSize - 1) / blockSize; - kernel_sharedHitCleaner<<>>( - hh.view(), tuples_d, tracks_d, quality_d, params_.minHitsForSharingCut_, device_hitToTuple_.get()); + + kernel_rejectDuplicate<<>>(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + + kernel_sharedHitCleaner<<>>(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + + if (params_.useSimpleTripletCleaner_) { + kernel_simpleTripletCleaner<<>>(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + } else { + kernel_tripletCleaner<<>>(hh.view(), + tuples_d, + tracks_d, + quality_d, + params_.minHitsForSharingCut_, + params_.dupPassThrough_, + device_hitToTuple_.get()); + } cudaCheck(cudaGetLastError()); +#ifdef GPU_DEBUG + cudaCheck(cudaDeviceSynchronize()); +#endif } if (params_.doStats_) { diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.h b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.h index f5c2755ea8ac6..2e6a7b6293873 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernels.h @@ -17,6 +17,7 @@ namespace cAHitNtupletGenerator { unsigned long long nCells; unsigned long long nTuples; unsigned long long nFitTracks; + unsigned long long nLooseTracks; unsigned long long nGoodTracks; unsigned long long nUsedHits; unsigned long long nDupHits; @@ -51,7 +52,7 @@ namespace cAHitNtupletGenerator { Region quadruplet; }; - // params + // params (FIXME: thi si a POD: so no constructor no traling _ and no const as params_ is already const) struct Params { Params(bool onGPU, uint32_t minHitsPerNtuplet, @@ -68,6 +69,8 @@ namespace cAHitNtupletGenerator { bool doZ0Cut, bool doPtCut, bool doSharedHitCut, + bool dupPassThrough, + bool useSimpleTripletCleaner, float ptmin, float CAThetaCutBarrel, float CAThetaCutForward, @@ -91,6 +94,8 @@ namespace cAHitNtupletGenerator { doZ0Cut_(doZ0Cut), doPtCut_(doPtCut), doSharedHitCut_(doSharedHitCut), + dupPassThrough_(dupPassThrough), + useSimpleTripletCleaner_(useSimpleTripletCleaner), ptmin_(ptmin), CAThetaCutBarrel_(CAThetaCutBarrel), CAThetaCutForward_(CAThetaCutForward), @@ -114,6 +119,8 @@ namespace cAHitNtupletGenerator { const bool doZ0Cut_; const bool doPtCut_; const bool doSharedHitCut_; + const bool dupPassThrough_; + const bool useSimpleTripletCleaner_; const float ptmin_; const float CAThetaCutBarrel_; const float CAThetaCutForward_; diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsImpl.h b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsImpl.h index 550541fdca6fb..852eb6a5c1757 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsImpl.h +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorKernelsImpl.h @@ -7,6 +7,7 @@ #include #include +#include #include @@ -30,6 +31,14 @@ using Quality = pixelTrack::Quality; using TkSoA = pixelTrack::TrackSoA; using HitContainer = pixelTrack::HitContainer; +namespace { + + constexpr uint16_t tkNotFound = std::numeric_limits::max(); + constexpr float maxScore = std::numeric_limits::max(); + constexpr float nSigma2 = 25.f; + +} // namespace + __global__ void kernel_checkOverflows(HitContainer const *foundNtuplets, caConstants::TupleMultiplicity const *tupleMultiplicity, CAHitNtupletGeneratorKernelsGPU::HitToTuple const *hitToTuple, @@ -62,7 +71,7 @@ __global__ void kernel_checkOverflows(HitContainer const *foundNtuplets, apc->get().n, nHits, hitToTuple->totOnes()); - if (apc->get().m < caConstants::maxNumberOfQuadruplets()) { + if (apc->get().m < caConstants::maxNumberOfQuadruplets) { assert(foundNtuplets->size(apc->get().m) == 0); assert(foundNtuplets->size() == apc->get().n); } @@ -111,7 +120,7 @@ __global__ void kernel_checkOverflows(HitContainer const *foundNtuplets, } __global__ void kernel_fishboneCleaner(GPUCACell const *cells, uint32_t const *__restrict__ nCells, Quality *quality) { - constexpr auto bad = pixelTrack::Quality::bad; + constexpr auto reject = pixelTrack::Quality::dup; auto first = threadIdx.x + blockIdx.x * blockDim.x; for (int idx = first, nt = (*nCells); idx < nt; idx += gridDim.x * blockDim.x) { @@ -120,17 +129,19 @@ __global__ void kernel_fishboneCleaner(GPUCACell const *cells, uint32_t const *_ continue; for (auto it : thisCell.tracks()) - quality[it] = bad; + quality[it] = reject; } } +// remove shorter tracks if sharing a cell +// It does not seem to affect efficiency in any way! __global__ void kernel_earlyDuplicateRemover(GPUCACell const *cells, uint32_t const *__restrict__ nCells, HitContainer *foundNtuplets, - Quality *quality) { - // constexpr auto bad = trackQuality::bad; - constexpr auto dup = pixelTrack::Quality::dup; - // constexpr auto loose = trackQuality::loose; + Quality *quality, + bool dupPassThrough) { + // quality to mark rejected + constexpr auto reject = pixelTrack::Quality::edup; /// cannot be loose assert(nCells); auto first = threadIdx.x + blockIdx.x * blockDim.x; @@ -150,19 +161,24 @@ __global__ void kernel_earlyDuplicateRemover(GPUCACell const *cells, maxNh = std::max(nh, maxNh); } + // quad pass through (leave it her for tests) + // maxNh = std::min(4U, maxNh); + for (auto it : thisCell.tracks()) { - if (foundNtuplets->size(it) != maxNh) - quality[it] = dup; //no race: simple assignment of the same constant + if (foundNtuplets->size(it) < maxNh) + quality[it] = reject; //no race: simple assignment of the same constant } } } +// assume the above (so, short tracks already removed) __global__ void kernel_fastDuplicateRemover(GPUCACell const *__restrict__ cells, uint32_t const *__restrict__ nCells, HitContainer const *__restrict__ foundNtuplets, - TkSoA *__restrict__ tracks) { - constexpr auto bad = pixelTrack::Quality::bad; - constexpr auto dup = pixelTrack::Quality::dup; + TkSoA *__restrict__ tracks, + bool dupPassThrough) { + // quality to mark rejected + auto const reject = dupPassThrough ? pixelTrack::Quality::loose : pixelTrack::Quality::dup; constexpr auto loose = pixelTrack::Quality::loose; assert(nCells); @@ -174,25 +190,81 @@ __global__ void kernel_fastDuplicateRemover(GPUCACell const *__restrict__ cells, continue; // if (thisCell.theDoubletId < 0) continue; - float mc = 10000.f; - uint16_t im = 60000; + float mc = maxScore; + uint16_t im = tkNotFound; + /* chi2 penalize higher-pt tracks (try rescale it?) auto score = [&](auto it) { - return std::abs(tracks->tip(it)); // tip - // return tracks->chi2(it); //chi2 + return foundNtuplets->size(it) < 4 ? + std::abs(tracks->tip(it)) : // tip for triplets + tracks->chi2(it); //chi2 for quads }; + */ + + auto score = [&](auto it) { return std::abs(tracks->tip(it)); }; + + // full crazy combinatorics + int ntr = thisCell.tracks().size(); + for (int i = 0; i < ntr; ++i) { + auto it = thisCell.tracks()[i]; + auto qi = tracks->quality(it); + if (qi <= reject) + continue; + auto opi = tracks->stateAtBS.state(it)(2); + auto e2opi = tracks->stateAtBS.covariance(it)(9); + auto cti = tracks->stateAtBS.state(it)(3); + auto e2cti = tracks->stateAtBS.covariance(it)(12); + for (auto j = i + 1; j < ntr; ++j) { + auto jt = thisCell.tracks()[j]; + auto qj = tracks->quality(jt); + if (qj <= reject) + continue; +#ifdef GPU_DEBUG + if (foundNtuplets->size(it) != foundNtuplets->size(jt)) + printf(" a mess\n"); +#endif + auto opj = tracks->stateAtBS.state(jt)(2); + auto ctj = tracks->stateAtBS.state(jt)(3); + auto dct = nSigma2 * (tracks->stateAtBS.covariance(jt)(12) + e2cti); + if ((cti - ctj) * (cti - ctj) > dct) + continue; + auto dop = nSigma2 * (tracks->stateAtBS.covariance(jt)(9) + e2opi); + if ((opi - opj) * (opi - opj) > dop) + continue; + if ((qj < qi) || (qj == qi && score(it) < score(jt))) + tracks->quality(jt) = reject; + else { + tracks->quality(it) = reject; + break; + } + } + } + + // find maxQual + auto maxQual = reject; // no duplicate! + for (auto it : thisCell.tracks()) { + if (tracks->quality(it) > maxQual) + maxQual = tracks->quality(it); + } + + if (maxQual <= loose) + continue; // find min score for (auto it : thisCell.tracks()) { - if (tracks->quality(it) == loose && score(it) < mc) { + if (tracks->quality(it) == maxQual && score(it) < mc) { mc = score(it); im = it; } } - // mark all other duplicates + + if (tkNotFound == im) + continue; + + // mark all other duplicates (not yet, keep it loose) for (auto it : thisCell.tracks()) { - if (tracks->quality(it) != bad && it != im) - tracks->quality(it) = dup; //no race: simple assignment of the same constant + if (tracks->quality(it) > loose && it != im) + tracks->quality(it) = loose; //no race: simple assignment of the same constant } } } @@ -310,7 +382,7 @@ __global__ void kernel_countMultiplicity(HitContainer const *__restrict__ foundN auto nhits = foundNtuplets->size(it); if (nhits < 3) continue; - if (quality[it] == pixelTrack::Quality::dup) + if (quality[it] == pixelTrack::Quality::edup) continue; assert(quality[it] == pixelTrack::Quality::bad); if (nhits > 5) @@ -328,7 +400,7 @@ __global__ void kernel_fillMultiplicity(HitContainer const *__restrict__ foundNt auto nhits = foundNtuplets->size(it); if (nhits < 3) continue; - if (quality[it] == pixelTrack::Quality::dup) + if (quality[it] == pixelTrack::Quality::edup) continue; assert(quality[it] == pixelTrack::Quality::bad); if (nhits > 5) @@ -349,7 +421,7 @@ __global__ void kernel_classifyTracks(HitContainer const *__restrict__ tuples, break; // guard // if duplicate: not even fit - if (quality[it] == pixelTrack::Quality::dup) + if (quality[it] == pixelTrack::Quality::edup) continue; assert(quality[it] == pixelTrack::Quality::bad); @@ -370,6 +442,8 @@ __global__ void kernel_classifyTracks(HitContainer const *__restrict__ tuples, continue; } + quality[it] = pixelTrack::Quality::strict; + // compute a pT-dependent chi2 cut // default parameters: // - chi2MaxPt = 10 GeV @@ -392,6 +466,8 @@ __global__ void kernel_classifyTracks(HitContainer const *__restrict__ tuples, continue; } + quality[it] = pixelTrack::Quality::tight; + // impose "region cuts" based on the fit results (phi, Tip, pt, cotan(theta)), Zip) // default cuts: // - for triplets: |Tip| < 0.3 cm, pT > 0.5 GeV, |Zip| < 12.0 cm @@ -402,7 +478,7 @@ __global__ void kernel_classifyTracks(HitContainer const *__restrict__ tuples, (std::abs(tracks->zip(it)) < region.maxZip); if (isOk) - quality[it] = pixelTrack::Quality::loose; + quality[it] = pixelTrack::Quality::highPurity; } } @@ -413,7 +489,10 @@ __global__ void kernel_doStatsForTracks(HitContainer const *__restrict__ tuples, for (int idx = first, ntot = tuples->nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { if (tuples->size(idx) == 0) break; //guard - if (quality[idx] != pixelTrack::Quality::loose) + if (quality[idx] < pixelTrack::Quality::loose) + continue; + atomicAdd(&(counters->nLooseTracks), 1); + if (quality[idx] < pixelTrack::Quality::strict) continue; atomicAdd(&(counters->nGoodTracks), 1); } @@ -426,8 +505,6 @@ __global__ void kernel_countHitInTracks(HitContainer const *__restrict__ tuples, for (int idx = first, ntot = tuples->nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { if (tuples->size(idx) == 0) break; // guard - if (quality[idx] != pixelTrack::Quality::loose) - continue; for (auto h = tuples->begin(idx); h != tuples->end(idx); ++h) hitToTuple->count(*h); } @@ -440,8 +517,6 @@ __global__ void kernel_fillHitInTracks(HitContainer const *__restrict__ tuples, for (int idx = first, ntot = tuples->nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { if (tuples->size(idx) == 0) break; // guard - if (quality[idx] != pixelTrack::Quality::loose) - continue; for (auto h = tuples->begin(idx); h != tuples->end(idx); ++h) hitToTuple->fill(*h, idx); } @@ -477,19 +552,144 @@ __global__ void kernel_doStatsForHitInTracks(CAHitNtupletGeneratorKernelsGPU::Hi } } +__global__ void kernel_countSharedHit(int *__restrict__ nshared, + HitContainer const *__restrict__ ptuples, + Quality const *__restrict__ quality, + CAHitNtupletGeneratorKernelsGPU::HitToTuple const *__restrict__ phitToTuple) { + constexpr auto loose = pixelTrack::Quality::loose; + + auto &hitToTuple = *phitToTuple; + auto const &foundNtuplets = *ptuples; + + int first = blockDim.x * blockIdx.x + threadIdx.x; + for (int idx = first, ntot = hitToTuple.nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { + if (hitToTuple.size(idx) < 2) + continue; + + int nt = 0; + + // count "good" tracks + for (auto it = hitToTuple.begin(idx); it != hitToTuple.end(idx); ++it) { + if (quality[*it] < loose) + continue; + ++nt; + } + + if (nt < 2) + continue; + + // now mark each track triplet as sharing a hit + for (auto it = hitToTuple.begin(idx); it != hitToTuple.end(idx); ++it) { + if (foundNtuplets.size(*it) > 3) + continue; + atomicAdd(&nshared[*it], 1); + } + + } // hit loop +} + +__global__ void kernel_markSharedHit(int const *__restrict__ nshared, + HitContainer const *__restrict__ tuples, + Quality *__restrict__ quality, + bool dupPassThrough) { + // constexpr auto bad = pixelTrack::Quality::bad; + constexpr auto dup = pixelTrack::Quality::dup; + constexpr auto loose = pixelTrack::Quality::loose; + // constexpr auto strict = pixelTrack::Quality::strict; + + // quality to mark rejected + auto const reject = dupPassThrough ? loose : dup; + + int first = blockDim.x * blockIdx.x + threadIdx.x; + for (int idx = first, ntot = tuples->nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { + if (tuples->size(idx) == 0) + break; //guard + if (quality[idx] <= reject) + continue; + if (nshared[idx] > 2) + quality[idx] = reject; + } +} + +// mostly for very forward triplets..... +__global__ void kernel_rejectDuplicate(TrackingRecHit2DSOAView const *__restrict__ hhp, + HitContainer const *__restrict__ ptuples, + TkSoA const *__restrict__ ptracks, + Quality *__restrict__ quality, + uint16_t nmin, + bool dupPassThrough, + CAHitNtupletGeneratorKernelsGPU::HitToTuple const *__restrict__ phitToTuple) { + // quality to mark rejected + auto const reject = dupPassThrough ? pixelTrack::Quality::loose : pixelTrack::Quality::dup; + + auto &hitToTuple = *phitToTuple; + auto const &foundNtuplets = *ptuples; + auto const &tracks = *ptracks; + + int first = blockDim.x * blockIdx.x + threadIdx.x; + for (int idx = first, ntot = hitToTuple.nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { + if (hitToTuple.size(idx) < 2) + continue; + + /* chi2 is bad for large pt + auto score = [&](auto it, auto nh) { + return nh < 4 ? std::abs(tracks.tip(it)) : // tip for triplets + tracks.chi2(it); //chi2 + }; + */ + auto score = [&](auto it, auto nh) { return std::abs(tracks.tip(it)); }; + + // full combinatorics + for (auto ip = hitToTuple.begin(idx); ip != hitToTuple.end(idx); ++ip) { + auto const it = *ip; + auto qi = quality[it]; + if (qi <= reject) + continue; + auto opi = tracks.stateAtBS.state(it)(2); + auto e2opi = tracks.stateAtBS.covariance(it)(9); + auto cti = tracks.stateAtBS.state(it)(3); + auto e2cti = tracks.stateAtBS.covariance(it)(12); + auto nhi = foundNtuplets.size(it); + for (auto jp = ip + 1; jp != hitToTuple.end(idx); ++jp) { + auto const jt = *jp; + auto qj = quality[jt]; + if (qj <= reject) + continue; + auto opj = tracks.stateAtBS.state(jt)(2); + auto ctj = tracks.stateAtBS.state(jt)(3); + auto dct = nSigma2 * (tracks.stateAtBS.covariance(jt)(12) + e2cti); + if ((cti - ctj) * (cti - ctj) > dct) + continue; + auto dop = nSigma2 * (tracks.stateAtBS.covariance(jt)(9) + e2opi); + if ((opi - opj) * (opi - opj) > dop) + continue; + auto nhj = foundNtuplets.size(jt); + if (nhj < nhi || (nhj == nhi && (qj < qi || (qj == qi && score(it, nhi) < score(jt, nhj))))) + quality[jt] = reject; + else { + quality[it] = reject; + break; + } + } + } + } +} + __global__ void kernel_sharedHitCleaner(TrackingRecHit2DSOAView const *__restrict__ hhp, HitContainer const *__restrict__ ptuples, TkSoA const *__restrict__ ptracks, Quality *__restrict__ quality, uint16_t nmin, + bool dupPassThrough, CAHitNtupletGeneratorKernelsGPU::HitToTuple const *__restrict__ phitToTuple) { - constexpr auto bad = pixelTrack::Quality::bad; - constexpr auto dup = pixelTrack::Quality::dup; - // constexpr auto loose = trackQuality::loose; + // quality to mark rejected + auto const reject = dupPassThrough ? pixelTrack::Quality::loose : pixelTrack::Quality::dup; + // quality of longest track + auto const longTqual = pixelTrack::Quality::highPurity; auto &hitToTuple = *phitToTuple; auto const &foundNtuplets = *ptuples; - auto const &tracks = *ptracks; + // auto const &tracks = *ptracks; auto const &hh = *hhp; int l1end = hh.hitsLayerStart()[1]; @@ -499,16 +699,23 @@ __global__ void kernel_sharedHitCleaner(TrackingRecHit2DSOAView const *__restric if (hitToTuple.size(idx) < 2) continue; - float mc = 10000.f; - uint16_t im = 60000; uint32_t maxNh = 0; // find maxNh for (auto it = hitToTuple.begin(idx); it != hitToTuple.end(idx); ++it) { + if (quality[*it] < longTqual) + continue; uint32_t nh = foundNtuplets.size(*it); maxNh = std::max(nh, maxNh); } - // kill all tracks shorter than maxHn (only triplets???) + + if (maxNh < 4) + continue; + + // quad pass through (leave for tests) + // maxNh = std::min(4U, maxNh); + + // kill all tracks shorter than maxHn (only triplets??? for (auto it = hitToTuple.begin(idx); it != hitToTuple.end(idx); ++it) { uint32_t nh = foundNtuplets.size(*it); @@ -516,26 +723,115 @@ __global__ void kernel_sharedHitCleaner(TrackingRecHit2DSOAView const *__restric if (idx < l1end and nh > nmin) continue; - if (maxNh != nh) - quality[*it] = dup; + if (nh < maxNh && quality[*it] > reject) + quality[*it] = reject; + } + } +} + +__global__ void kernel_tripletCleaner(TrackingRecHit2DSOAView const *__restrict__ hhp, + HitContainer const *__restrict__ ptuples, + TkSoA const *__restrict__ ptracks, + Quality *__restrict__ quality, + uint16_t nmin, + bool dupPassThrough, + CAHitNtupletGeneratorKernelsGPU::HitToTuple const *__restrict__ phitToTuple) { + // quality to mark rejected + auto const reject = pixelTrack::Quality::loose; + /// min quality of good + auto const good = pixelTrack::Quality::strict; + + auto &hitToTuple = *phitToTuple; + auto const &foundNtuplets = *ptuples; + auto const &tracks = *ptracks; + + int first = blockDim.x * blockIdx.x + threadIdx.x; + for (int idx = first, ntot = hitToTuple.nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { + if (hitToTuple.size(idx) < 2) + continue; + + float mc = maxScore; + uint16_t im = tkNotFound; + uint32_t maxNh = 0; + + // find maxNh + for (auto it = hitToTuple.begin(idx); it != hitToTuple.end(idx); ++it) { + if (quality[*it] <= good) + continue; + uint32_t nh = foundNtuplets.size(*it); + maxNh = std::max(nh, maxNh); + } + + // only triplets + if (maxNh != 3) + continue; + + // for triplets choose best tip! (should we first find best quality???) + for (auto ip = hitToTuple.begin(idx); ip != hitToTuple.end(idx); ++ip) { + auto const it = *ip; + if (quality[it] >= good && std::abs(tracks.tip(it)) < mc) { + mc = std::abs(tracks.tip(it)); + im = it; + } } - if (maxNh > 3) + if (tkNotFound == im) continue; - // for triplets choose best tip! + + // mark worse ambiguities for (auto ip = hitToTuple.begin(idx); ip != hitToTuple.end(idx); ++ip) { auto const it = *ip; - if (quality[it] != bad && std::abs(tracks.tip(it)) < mc) { + if (quality[it] > reject && it != im) + quality[it] = reject; //no race: simple assignment of the same constant + } + + } // loop over hits +} + +__global__ void kernel_simpleTripletCleaner( + TrackingRecHit2DSOAView const *__restrict__ hhp, + HitContainer const *__restrict__ ptuples, + TkSoA const *__restrict__ ptracks, + Quality *__restrict__ quality, + uint16_t nmin, + bool dupPassThrough, + CAHitNtupletGeneratorKernelsGPU::HitToTuple const *__restrict__ phitToTuple) { + // quality to mark rejected + auto const reject = pixelTrack::Quality::loose; + /// min quality of good + auto const good = pixelTrack::Quality::loose; + + auto &hitToTuple = *phitToTuple; + auto const &foundNtuplets = *ptuples; + auto const &tracks = *ptracks; + + int first = blockDim.x * blockIdx.x + threadIdx.x; + for (int idx = first, ntot = hitToTuple.nOnes(); idx < ntot; idx += gridDim.x * blockDim.x) { + if (hitToTuple.size(idx) < 2) + continue; + + float mc = maxScore; + uint16_t im = tkNotFound; + + // choose best tip! (should we first find best quality???) + for (auto ip = hitToTuple.begin(idx); ip != hitToTuple.end(idx); ++ip) { + auto const it = *ip; + if (quality[it] >= good && std::abs(tracks.tip(it)) < mc) { mc = std::abs(tracks.tip(it)); im = it; } } - // mark duplicates + + if (tkNotFound == im) + continue; + + // mark worse ambiguities for (auto ip = hitToTuple.begin(idx); ip != hitToTuple.end(idx); ++ip) { auto const it = *ip; - if (quality[it] != bad && it != im) - quality[it] = dup; //no race: simple assignment of the same constant + if (quality[it] > reject && foundNtuplets.size(it) == 3 && it != im) + quality[it] = reject; //no race: simple assignment of the same constant } + } // loop over hits } @@ -576,14 +872,16 @@ __global__ void kernel_print_found_ntuplets(TrackingRecHit2DSOAView const *__res __global__ void kernel_printCounters(cAHitNtupletGenerator::Counters const *counters) { auto const &c = *counters; printf( - "||Counters | nEvents | nHits | nCells | nTuples | nFitTacks | nGoodTracks | nUsedHits | nDupHits | " + "||Counters | nEvents | nHits | nCells | nTuples | nFitTacks | nLooseTracks | nGoodTracks | nUsedHits | " + "nDupHits | " "nKilledCells | " "nEmptyCells | nZeroTrackCells ||\n"); - printf("Counters Raw %lld %lld %lld %lld %lld %lld %lld %lld %lld %lld %lld\n", + printf("Counters Raw %lld %lld %lld %lld %lld %lld %lld %lld %lld %lld %lld %lld\n", c.nEvents, c.nHits, c.nCells, c.nTuples, + c.nLooseTracks, c.nGoodTracks, c.nFitTracks, c.nUsedHits, @@ -591,12 +889,13 @@ __global__ void kernel_printCounters(cAHitNtupletGenerator::Counters const *coun c.nKilledCells, c.nEmptyCells, c.nZeroTrackCells); - printf("Counters Norm %lld || %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.3f| %.3f||\n", + printf("Counters Norm %lld || %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.1f| %.3f| %.3f||\n", c.nEvents, c.nHits / double(c.nEvents), c.nCells / double(c.nEvents), c.nTuples / double(c.nEvents), c.nFitTracks / double(c.nEvents), + c.nLooseTracks / double(c.nEvents), c.nGoodTracks / double(c.nEvents), c.nUsedHits / double(c.nEvents), c.nDupHits / double(c.nEvents), diff --git a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorOnGPU.cc b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorOnGPU.cc index 3f5ba1d04e7db..9d542c40f2c80 100644 --- a/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorOnGPU.cc +++ b/RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletGeneratorOnGPU.cc @@ -71,6 +71,8 @@ CAHitNtupletGeneratorOnGPU::CAHitNtupletGeneratorOnGPU(const edm::ParameterSet& cfg.getParameter("doZ0Cut"), cfg.getParameter("doPtCut"), cfg.getParameter("doSharedHitCut"), + cfg.getParameter("dupPassThrough"), + cfg.getParameter("useSimpleTripletCleaner"), cfg.getParameter("ptmin"), cfg.getParameter("CAThetaCutBarrel"), cfg.getParameter("CAThetaCutForward"), @@ -154,6 +156,8 @@ void CAHitNtupletGeneratorOnGPU::fillDescriptions(edm::ParameterSetDescription& desc.add("doPtCut", true); desc.add("useRiemannFit", false)->setComment("true for Riemann, false for BrokenLine"); desc.add("doSharedHitCut", true)->setComment("Sharing hit nTuples cleaning"); + desc.add("dupPassThrough", false)->setComment("Do not reject duplicate"); + desc.add("useSimpleTripletCleaner", false)->setComment("use alternate implementation"); edm::ParameterSetDescription trackQualityCuts; trackQualityCuts.add("chi2MaxPt", 10.)->setComment("max pT used to determine the pT-dependent chi2 cut"); diff --git a/RecoPixelVertexing/PixelVertexFinding/plugins/gpuVertexFinder.cc b/RecoPixelVertexing/PixelVertexFinding/plugins/gpuVertexFinder.cc index d685ced488233..8f12a7a596378 100644 --- a/RecoPixelVertexing/PixelVertexFinding/plugins/gpuVertexFinder.cc +++ b/RecoPixelVertexing/PixelVertexFinding/plugins/gpuVertexFinder.cc @@ -1,6 +1,3 @@ -#ifndef RecoPixelVertexing_PixelVertexFinding_plugins_gpuVertexFinderImpl_h -#define RecoPixelVertexing_PixelVertexFinding_plugins_gpuVertexFinderImpl_h - #include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h" #include "gpuClusterTracksByDensity.h" @@ -39,7 +36,7 @@ namespace gpuVertexFinder { if (nHits < 4) continue; // no triplets - if (quality[idx] != pixelTrack::Quality::loose) + if (quality[idx] < pixelTrack::Quality::highPurity) continue; auto pt = tracks.pt(idx); @@ -188,5 +185,3 @@ namespace gpuVertexFinder { } } // namespace gpuVertexFinder - -#endif // RecoPixelVertexing_PixelVertexFinding_plugins_gpuVertexFinderImpl_h diff --git a/Validation/RecoTrack/plugins/MultiTrackValidator.cc b/Validation/RecoTrack/plugins/MultiTrackValidator.cc index 41d2a2d578072..d7347a52ba992 100644 --- a/Validation/RecoTrack/plugins/MultiTrackValidator.cc +++ b/Validation/RecoTrack/plugins/MultiTrackValidator.cc @@ -1033,13 +1033,16 @@ void MultiTrackValidator::dqmAnalyze(const edm::Event& event, size_t n_selTrack_dr = 0; //calculate dR for tracks + declareDynArray(float, trackCollection.size(), dR_trk); + declareDynArray(float, trackCollection.size(), dR_trk_jet); +#ifndef NO_TRACK_DR + // this accounts for most of the time spent in MTV and it is used to fill just one histo that is of doubtful usefulness but (maybe) for the whole collection const edm::View* trackCollectionDr = &trackCollection; if (calculateDrSingleCollection_) { trackCollectionDr = trackCollectionForDrCalculation.product(); } - declareDynArray(float, trackCollection.size(), dR_trk); - declareDynArray(float, trackCollection.size(), dR_trk_jet); trackDR(trackCollection, *trackCollectionDr, dR_trk, dR_trk_jet, coresVector); +#endif for (View::size_type i = 0; i < trackCollection.size(); ++i) { auto track = trackCollection.refAt(i); diff --git a/Validation/RecoTrack/python/TrackValidation_cff.py b/Validation/RecoTrack/python/TrackValidation_cff.py index 74c91bd011922..b0a933365082e 100644 --- a/Validation/RecoTrack/python/TrackValidation_cff.py +++ b/Validation/RecoTrack/python/TrackValidation_cff.py @@ -978,7 +978,7 @@ def _uniqueFirstLayers(layerList): tracksValidationSeedSelectorsTrackingOnly ) - +#################################################################################################### ### Pixel tracking only mode (placeholder for now) trackingParticlePixelTrackAsssociation = trackingParticleRecoTrackAsssociation.clone( label_tr = "pixelTracks", @@ -992,13 +992,44 @@ def _uniqueFirstLayers(layerList): src = "pixelTracks", vertexTag = "pixelVertices", ) -pixelTracksPt09 = generalTracksPt09.clone(quality = ["undefQuality"], **_pixelTracksCustom) -pixelTracksFromPV = generalTracksFromPV.clone(quality = "undefQuality", **_pixelTracksCustom) -pixelTracksFromPVPt09 = pixelTracksPt09.clone(src = "pixelTracksFromPV") + +trackRefSelector = cms.EDFilter('TrackRefSelector', + src = cms.InputTag('pixelTracks'), + cut = cms.string("") +) + +trackSelector = cms.EDFilter('TrackSelector', + src = cms.InputTag('pixelTracks'), + cut = cms.string("") +) + +cutstring = "numberOfValidHits == 3" +pixelTracks3hits = trackRefSelector.clone( cut = cutstring ) + +cutstring = "numberOfValidHits >= 4" +pixelTracks4hits = trackRefSelector.clone( cut = cutstring ) + +cutstring = "pt > 0.9" +pixelTracksPt09 = trackRefSelector.clone( cut = cutstring ) +#pixelTracksPt09 = generalTracksPt09.clone(quality = ["undefQuality"], **_pixelTracksCustom) + +pixelTracksFromPV = generalTracksFromPV.clone(quality = "highPurity", ptMin = 0.0, ptMax = 99999., ptErrorCut = 99999., copyExtras = True, **_pixelTracksCustom) +#pixelTracksFromPVPt09 = generalTracksPt09.clone(quality = ["loose","tight","highPurity"], vertexTag = "pixelVertices", src = "pixelTracksFromPV") +pixelTracksFromPVPt09 = pixelTracksFromPV.clone(ptMin = 0.9) + +cutstring = "numberOfValidHits >= 4" +#pixelTracksFromPV4hits = trackRefSelector.clone( cut = cutstring, src = "pixelTracksFromPV" ) +pixelTracksFromPV4hits = pixelTracksFromPV.clone( numberOfValidPixelHits = 4 ) + trackValidatorPixelTrackingOnly = trackValidator.clone( dirName = "Tracking/PixelTrack/", - label = ["pixelTracks", "pixelTracksPt09"], + label = [ + "pixelTracks", "pixelTracksPt09", "pixelTracks3hits", "pixelTracks4hits", + "pixelTracksL", "pixelTracksPt09L", "pixelTracks3hitsL", "pixelTracks4hitsL", + "pixelTracksT", "pixelTracksPt09T", "pixelTracks3hitsT", "pixelTracks4hitsT", + "pixelTracksHP", "pixelTracksPt09HP", "pixelTracks3hitsHP", "pixelTracks4hitsHP", + ], doResolutionPlotsForLabels = [], trackCollectionForDrCalculation = "pixelTracks", associators = ["trackingParticlePixelTrackAsssociation"], @@ -1009,7 +1040,12 @@ def _uniqueFirstLayers(layerList): ) trackValidatorFromPVPixelTrackingOnly = trackValidatorPixelTrackingOnly.clone( dirName = "Tracking/PixelTrackFromPV/", - label = ["pixelTracksFromPV", "pixelTracksFromPVPt09"], + label = [ + "pixelTracksFromPV", "pixelTracksFromPVPt09", "pixelTracksFromPV4hits", + "pixelTracksFromPVL", "pixelTracksFromPVT", "pixelTracksFromPVHP", + "pixelTracksFromPVPt09L", "pixelTracksFromPVPt09T", "pixelTracksFromPVPt09HP", + "pixelTracksFromPV4hitsL", "pixelTracksFromPV4hitsT", "pixelTracksFromPV4hitsHP", + ], label_tp_effic = "trackingParticlesSignal", label_tp_fake = "trackingParticlesSignal", label_tp_effic_refvector = True, @@ -1030,6 +1066,12 @@ def _uniqueFirstLayers(layerList): ) trackValidatorBHadronPixelTrackingOnly = trackValidatorPixelTrackingOnly.clone( dirName = "Tracking/PixelTrackBHadron/", + label = [ + "pixelTracks", "pixelTracksPt09", + "pixelTracksL", "pixelTracks3hitsL", "pixelTracks4hitsL", + "pixelTracksT", "pixelTracks3hitsT", "pixelTracks4hitsT", + "pixelTracksHP", "pixelTracks3hitsHP", "pixelTracks4hitsHP", + ], label_tp_effic = "trackingParticlesBHadron", label_tp_effic_refvector = True, doSimPlots = True, @@ -1041,14 +1083,64 @@ def _uniqueFirstLayers(layerList): tracksValidationTruthPixelTrackingOnly.replace(trackingParticleRecoTrackAsssociation, trackingParticlePixelTrackAsssociation) tracksValidationTruthPixelTrackingOnly.replace(VertexAssociatorByPositionAndTracks, PixelVertexAssociatorByPositionAndTracks) tracksValidationTruthPixelTrackingOnly.add(trackingParticlesBHadron) +tracksValidationTruthPixelTrackingOnly.add( pixelTracks3hits ) +tracksValidationTruthPixelTrackingOnly.add( pixelTracks4hits ) +tracksValidationTruthPixelTrackingOnly.add( pixelTracksPt09 ) +tracksValidationTruthPixelTrackingOnly.add( pixelTracksFromPV ) +tracksValidationTruthPixelTrackingOnly.add( pixelTracksFromPVPt09 ) +tracksValidationTruthPixelTrackingOnly.add( pixelTracksFromPV4hits ) tracksPreValidationPixelTrackingOnly = cms.Task( tracksValidationTruthPixelTrackingOnly, - trackingParticlesSignal, - pixelTracksPt09, - pixelTracksFromPV, - pixelTracksFromPVPt09, -) + trackingParticlesSignal) + +##https://cmssdt.cern.ch/lxr/source/DataFormats/TrackReco/interface/TrackBase.h#0150 +quality = { + "L" : (1,"loose", ["loose","tight","highPurity"]), + "T" : (2,"tight", ["tight","highPurity"]), + "HP" : (4,"highPurity",["highPurity"]), +} + +for key,value in quality.items(): + qualityName = value[1] + qualityBit = value[0] + qualityList = value[2] + + label = "pixelTracks"+str(key) + cutstring = "qualityMask <= 7 & qualityMask >= " + str(qualityBit) + locals()[label] = trackRefSelector.clone( cut = cutstring ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) + + label = "pixelTracksFromPV"+key +# locals()[label] = generalTracksPt09.clone( ptMin = 0.0, vertexTag = "pixelVertices", src = "pixelTracksFromPV", quality = qualityList ) + locals()[label] = pixelTracksFromPV.clone( quality = qualityName ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) +#----------- + cutstring = "pt > 0.9 & qualityMask <= 7 & qualityMask >= " + str(qualityBit) + label = "pixelTracksPt09"+key + locals()[label] = trackRefSelector.clone( cut = cutstring ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) + + label = "pixelTracksFromPVPt09"+key + # locals()[label] = generalTracksPt09.clone( ptMin = 0.9, vertexTag = "pixelVertices", src = "pixelTracksFromPV", quality = qualityList ) + locals()[label] = pixelTracksFromPVPt09.clone( quality = qualityName ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) +#----------- + label = "pixelTracks4hits"+key + cutstring = "numberOfValidHits == 4 & qualityMask <= 7 & qualityMask >= " + str(qualityBit) + locals()[label] = trackRefSelector.clone( cut = cutstring ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) + + label = "pixelTracksFromPV4hits"+key +# locals()[label] = generalTracksPt09.clone( ptMin = 0.0, minHit = 4, vertexTag = "pixelVertices", src = "pixelTracksFromPV", quality = qualityList ) + locals()[label] = pixelTracksFromPV4hits.clone( quality = qualityName ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) +#-------- + label = "pixelTracks3hits"+key + cutstring = "numberOfValidHits == 3 & qualityMask <= 7 & qualityMask >= " + str(qualityBit) + locals()[label] = trackRefSelector.clone( cut = cutstring ) + tracksPreValidationPixelTrackingOnly.add(locals()[label]) + tracksValidationPixelTrackingOnly = cms.Sequence( trackValidatorPixelTrackingOnly + trackValidatorFromPVPixelTrackingOnly + @@ -1056,7 +1148,7 @@ def _uniqueFirstLayers(layerList): trackValidatorBHadronPixelTrackingOnly, tracksPreValidationPixelTrackingOnly ) - +#################################################################################################### ### Lite mode (only generalTracks and HP) trackValidatorLite = trackValidator.clone(