Introduce Quality for PixelTracks. Improve Patatrack duplicate & ambiguity resolution

CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h

@@ -1,14 +1,23 @@
 #ifndef CUDADataFormats_Track_TrackHeterogeneousT_H
 #define CUDADataFormats_Track_TrackHeterogeneousT_H
 
+#include <string>
+#include <algorithm>
+
 #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 };
-}
+  constexpr uint32_t qualitySize{7};
+  enum class Quality : uint8_t { bad = 0, edup, dup, loose, strict, tight, highPurity, notQuality = qualitySize };
+  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<Quality>(qp);
+  }
+}  // namespace pixelTrack
 
 template <int32_t S>
 class TrackSoAHeterogeneousT {

CUDADataFormats/Track/test/BuildFile.xml

@@ -1,5 +1,10 @@
 <use name="HeterogeneousCore/CUDAUtilities"/>
 
+<bin file="TrackSoAHeterogeneous_t.cpp">
+  <use name="eigen"/>
+  <flags CXXFLAGS="-g -DGPU_DEBUG"/>
+</bin>
+
 <bin file="TrajectoryStateSOA_t.cpp" name="cpuTrajectoryStateSOA_t">
   <use name="eigen"/>
   <flags CXXFLAGS="-g -DGPU_DEBUG"/>

CUDADataFormats/Track/test/TrackSoAHeterogeneous_t.cpp

@@ -0,0 +1,21 @@
+#include "CUDADataFormats/Track/interface/TrackSoAHeterogeneousT.h"
+
+#include <iostream>
+#include <cassert>
+
+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<pixelTrack::Quality>(i);
+    auto q = pixelTrack::qualityByName(pixelTrack::qualityName[i]);
+    assert(qt == q);
+  }
+
+  return 0;
+}

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
 )
+
+
+
+

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<reco::TrackCollection>();
-  produces<TrackingRecHitCollection>();
-  produces<reco::TrackExtraCollection>();
-}
+class PixelTrackProducer : public edm::stream::EDProducer<> {
+public:
+  explicit PixelTrackProducer(const edm::ParameterSet& cfg)
+      : theReconstruction(cfg, consumesCollector()), htTopoToken_(esConsumes()) {
+    edm::LogInfo("PixelTrackProducer") << " construction...";
+    produces<reco::TrackCollection>();
+    produces<TrackingRecHitCollection>();
+    produces<reco::TrackExtraCollection>();
+  }
+
+  ~PixelTrackProducer() override {}
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+    edm::ParameterSetDescription desc;
 
-PixelTrackProducer::~PixelTrackProducer() {}
+    desc.add<std::string>("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<std::string>("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<TrackerTopology> httopo;
-  es.get<TrackerTopologyRcd>().get(httopo);
+private:
+  PixelTrackReconstruction theReconstruction;
+  const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> htTopoToken_;
+};
 
-  // store tracks
-  storeTracks(ev, tracks, *httopo);
-}
+DEFINE_FWK_MODULE(PixelTrackProducer);

RecoPixelVertexing/PixelTrackFitting/plugins/PixelTrackProducerFromSoA.cc

@@ -63,6 +63,7 @@ class PixelTrackProducerFromSoA : public edm::global::EDProducer<> {
   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> 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<HMSstorage>(iConfig.getParameter<edm::InputTag>("pixelRecHitLegacySrc"))),
       idealMagneticFieldToken_(esConsumes()),
       ttTopoToken_(esConsumes()),
-      minNumberOfHits_(iConfig.getParameter<int>("minNumberOfHits")) {
+      minNumberOfHits_(iConfig.getParameter<int>("minNumberOfHits")),
+      minQuality_(pixelTrack::qualityByName(iConfig.getParameter<std::string>("minQuality"))) {
+  if (minQuality_ == pixelTrack::Quality::notQuality) {
+    throw cms::Exception("PixelTrackConfiguration")
+        << iConfig.getParameter<std::string>("minQuality") + " is not a pixelTrack::Quality";
+  }
+  if (minQuality_ < pixelTrack::Quality::dup) {
+    throw cms::Exception("PixelTrackConfiguration")
+        << iConfig.getParameter<std::string>("minQuality") + " not supported";
+  }
   produces<reco::TrackCollection>();
   produces<TrackingRecHitCollection>();
   produces<reco::TrackExtraCollection>();
@@ -85,13 +95,23 @@ void PixelTrackProducerFromSoA::fillDescriptions(edm::ConfigurationDescriptions
   desc.add<edm::InputTag>("trackSrc", edm::InputTag("pixelTracksSoA"));
   desc.add<edm::InputTag>("pixelRecHitLegacySrc", edm::InputTag("siPixelRecHitsPreSplittingLegacy"));
   desc.add<int>("minNumberOfHits", 0);
-
+  desc.add<std::string>("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<IndToEdm>();
@@ -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<reco::Track>(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);
   }