[NGT] Add MC truth labels support to NanoAOD reco::Tracks and introduce PixelTracks+RecHits ntuples

HLTrigger/NGTScouting/plugins/HLTTracksRecHitsTableProducer.cc

@@ -0,0 +1,124 @@
+/**  \class HLTTracksRecHitsTableProducer
+ *
+ *   \brief Produces a nanoAOD flat table with recHits information for HLT tracks
+ *
+ *   This producer creates a nanoAOD flat table containing the recHits global positions and errors
+ *   starting from a collection of reco::Track.
+ *   This data can be added as an extension to the HLTTracks table.
+ *   The maximum number of recHits per track is fixed to a configurable value;
+ *   if a track has more recHits, a warning is issued and the extra recHits are ignored.
+ *
+ *   \author Luca Ferragina (INFN BO), 2025
+ */
+
+#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
+#include "DataFormats/BeamSpot/interface/BeamSpot.h"
+#include "DataFormats/Common/interface/ValueMap.h"
+#include "DataFormats/NanoAOD/interface/FlatTable.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+class HLTTracksRecHitsTableProducer : public edm::stream::EDProducer<> {
+public:
+  explicit HLTTracksRecHitsTableProducer(const edm::ParameterSet&);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  const edm::EDGetTokenT<std::vector<reco::Track>> tracks_;
+
+  const std::string tableName_;
+  const unsigned int maxRecHits_;
+  const unsigned int precision_;
+  const bool skipNonExistingSrc_;
+};
+
+HLTTracksRecHitsTableProducer::HLTTracksRecHitsTableProducer(const edm::ParameterSet& params)
+    : tracks_(consumes<std::vector<reco::Track>>(params.getParameter<edm::InputTag>("tracksSrc"))),
+      tableName_(params.getParameter<std::string>("tableName")),
+      maxRecHits_(params.getParameter<uint>("maxRecHits")),
+      precision_(params.getParameter<int>("precision")),
+      skipNonExistingSrc_(params.getParameter<bool>("skipNonExistingSrc")) {
+  produces<nanoaod::FlatTable>(tableName_);
+}
+
+void HLTTracksRecHitsTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  auto tracksIn = iEvent.getHandle(tracks_);
+  const size_t nTracks = tracksIn.isValid() ? (*tracksIn).size() : 0;
+
+  static constexpr float default_value = std::numeric_limits<float>::quiet_NaN();
+
+  std::vector<float> globalX(maxRecHits_ * nTracks, default_value);
+  std::vector<float> globalY(maxRecHits_ * nTracks, default_value);
+  std::vector<float> globalZ(maxRecHits_ * nTracks, default_value);
+  std::vector<float> globalErrX(maxRecHits_ * nTracks, default_value);
+  std::vector<float> globalErrY(maxRecHits_ * nTracks, default_value);
+  std::vector<float> globalErrZ(maxRecHits_ * nTracks, default_value);
+
+  if (tracksIn.isValid() || !(this->skipNonExistingSrc_)) {
+    const auto& tracks = *tracksIn;
+    for (size_t tkIndex = 0; tkIndex < nTracks; ++tkIndex) {
+      const auto& track = tracks[tkIndex];
+      for (auto it = track.recHitsBegin(); it != track.recHitsEnd(); ++it) {
+        auto hit = *it;
+        auto globalPoint = hit->globalPosition();
+        auto globalError = hit->globalPositionError();
+        auto hitIndex = std::distance(track.recHitsBegin(), it);
+        if (hitIndex >= maxRecHits_) {
+          edm::LogWarning("HLTTracksRecHitsTableProducer")
+              << " Track " << tkIndex << " has more (" << track.recHitsSize() << ") than " << maxRecHits_
+              << " recHits, skipping the rest.";
+          break;
+        }
+        globalX[tkIndex * maxRecHits_ + hitIndex] = globalPoint.x();
+        globalY[tkIndex * maxRecHits_ + hitIndex] = globalPoint.y();
+        globalZ[tkIndex * maxRecHits_ + hitIndex] = globalPoint.z();
+        globalErrX[tkIndex * maxRecHits_ + hitIndex] = globalError.cxx();
+        globalErrY[tkIndex * maxRecHits_ + hitIndex] = globalError.cyy();
+        globalErrZ[tkIndex * maxRecHits_ + hitIndex] = globalError.czz();
+      }
+    }
+  } else {
+    edm::LogWarning("HLTTracksRecHitsTableProducer")
+        << " Invalid handle for " << tableName_ << " in tracks input collection";
+  }
+
+  assert(globalX.size() == globalY.size() && globalX.size() == globalZ.size() && globalX.size() == globalErrX.size() &&
+         globalX.size() == globalErrY.size() && globalX.size() == globalErrZ.size());
+
+  // Table for tracks recHits
+  auto tracksTable =
+      std::make_unique<nanoaod::FlatTable>(nTracks * maxRecHits_, tableName_, /*singleton*/ false, /*extension*/ false);
+  tracksTable->addColumn<float>("globalX", globalX, "RecHits global x coordinate", precision_);
+  tracksTable->addColumn<float>("globalY", globalY, "RecHits global y coordinate", precision_);
+  tracksTable->addColumn<float>("globalZ", globalZ, "RecHits global z coordinate", precision_);
+  tracksTable->addColumn<float>("globalErrX", globalErrX, "RecHits global x error", precision_);
+  tracksTable->addColumn<float>("globalErrY", globalErrY, "RecHits global y error", precision_);
+  tracksTable->addColumn<float>("globalErrZ", globalErrZ, "RecHits global z error", precision_);
+
+  iEvent.put(std::move(tracksTable), tableName_);
+}
+
+void HLTTracksRecHitsTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<uint>("maxRecHits", 16)->setComment("maximum number of recHits per track to be stored in the table");
+  desc.add<bool>("skipNonExistingSrc", false)
+      ->setComment("whether or not to skip producing the table on absent input product");
+  desc.add<std::string>("tableName")->setComment("name of the flat table ouput");
+  desc.add<edm::InputTag>("tracksSrc")->setComment("std::vector<reco::Track> input collection");
+  desc.add<int>("precision", 7);
+  descriptions.addWithDefaultLabel(desc);
+}
+
+// Define this as a plug-in
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(HLTTracksRecHitsTableProducer);

HLTrigger/NGTScouting/python/HLTNanoProducer_cff.py

@@ -77,8 +77,10 @@
 
 # Store PixelTracks objects
 NanoPixelTables = cms.Sequence(
-    hltPixelTrackTable
+    pixelTrackAssoc
+    + hltPixelTrackTable
     + hltPixelTrackExtTable
+    + hltPixelTrackRecHitsTable
 )
 
 # Store variables and associators for validation purposes
@@ -124,6 +126,12 @@
     + NanoValTables
 )
 
+# Phase-2 HLT Nano flavour for pixel tracking validation / optimization
+hltPixelOnlyNanoFlavour = cms.Sequence(
+    NanoGenTables
+    + NanoPixelTables
+)
+
 ######################################
 # Customization
 ######################################

HLTrigger/NGTScouting/python/hltTracks_cfi.py

@@ -1,14 +1,44 @@
 import FWCore.ParameterSet.Config as cms
 from PhysicsTools.NanoAOD.common_cff import *
 
+from PhysicsTools.NanoAOD.trackingAssocValueMapsProducer_cfi import trackingAssocValueMapsProducer
+
+tpSelectorPixelTracks = cms.PSet(
+    lip = cms.double(30.0),
+    chargedOnly = cms.bool(True),
+    stableOnly = cms.bool(False),
+    pdgId = cms.vint32(),
+    signalOnly = cms.bool(False),
+    intimeOnly = cms.bool(True),
+    minRapidity = cms.double(-4.5),
+    minHit = cms.int32(0),
+    ptMin = cms.double(0.9),
+    ptMax = cms.double(1e100),
+    maxRapidity = cms.double(4.5),
+    tip = cms.double(2.5),
+    invertRapidityCut = cms.bool(False),
+    minPhi = cms.double(-3.2),
+    maxPhi = cms.double(3.2),
+)
+
+pixelTrackAssoc = trackingAssocValueMapsProducer.clone(
+    trackCollection =  cms.InputTag("hltPhase2PixelTracks"),
+    associator = cms.InputTag("hltTrackAssociatorByHits"),
+    trackingParticles = cms.InputTag("mix", "MergedTrackTruth"),
+    tpSelectorPSet = tpSelectorPixelTracks,
+    storeTPKinematics = cms.bool(True),
+)
+
+from Configuration.ProcessModifiers.phase2CAExtension_cff import phase2CAExtension
+phase2CAExtension.toModify(pixelTrackAssoc, trackCollection = "hltPhase2PixelTracksCAExtension")
+
 hltPixelTrackTable = cms.EDProducer(
     "SimpleTriggerTrackFlatTableProducer",
     skipNonExistingSrc = cms.bool(True),
     src = cms.InputTag("hltPhase2PixelTracks"),
     cut = cms.string(""),
     name = cms.string("hltPixelTrack"),
     doc = cms.string("HLT Pixel Track information"),
-    extension = cms.bool(False),
     variables = cms.PSet(
         pt = Var("pt()", "float", doc = "p_T (GeV)"),
         eta = Var("eta()", "float", doc = "#eta"),
@@ -29,6 +59,37 @@
         #isLoose = Var("quality('loose')", "bool", doc = "Loose track flag"),
         isTight = Var("quality('tight')", "bool", doc = "Tight track flag"),
         isHighPurity = Var("quality('highPurity')", "bool", doc = "High-purity track flag"),
+        qoverp = Var("qoverp()", "float", doc = "q/p"),
+        dsz = Var("dsz()", "float", doc = "dsz (cm)"),
+        qoverpErr = Var("qoverpError()", "float", doc = ""),
+        ptErr = Var("ptError()", "float", doc = ""),
+        lambdaErr = Var("lambdaError()", "float", doc = ""),
+        dszErr = Var("dszError()", "float", doc = ""),
+        etaErr = Var("etaError()", "float", doc = ""),
+        phiErr = Var("phiError()", "float", doc = ""),
+    ),
+        externalVariables = cms.PSet(
+        matched   = cms.PSet(src = cms.InputTag("pixelTrackAssoc","matched"),
+                             doc = cms.string("1 if matched to a TrackingParticle"),
+                             type = cms.string("uint8")),
+        duplicate = cms.PSet(src = cms.InputTag("pixelTrackAssoc","duplicate"),
+                             doc = cms.string("1 if multiple reco tracks map to same TP"),
+                             type = cms.string("uint8")),
+        tpPdgId  = cms.PSet(src = cms.InputTag("pixelTrackAssoc","tpPdgId"),
+                             doc = cms.string("pdgId of matched TrackingParticle"),
+                             type = cms.string("int16")),
+        tpCharge = cms.PSet(src = cms.InputTag("pixelTrackAssoc","tpCharge"),
+                             doc = cms.string("charge of matched TrackingParticle"),
+                             type = cms.string("int16")),
+        tpPt     = cms.PSet(src = cms.InputTag("pixelTrackAssoc","tpPt"),
+                             doc = cms.string("pt of matched TrackingParticle"),
+                             type = cms.string("float")),
+        tpEta    = cms.PSet(src = cms.InputTag("pixelTrackAssoc","tpEta"),
+                             doc = cms.string("eta of matched TrackingParticle"),
+                             type = cms.string("float")),
+        tpPhi    = cms.PSet(src = cms.InputTag("pixelTrackAssoc","tpPhi"),
+                             doc = cms.string("phi of matched TrackingParticle"),
+                             type = cms.string("float"))
     )
 )
 
@@ -39,6 +100,19 @@
                                        beamSpot = cms.InputTag("hltOnlineBeamSpot"),
                                        precision = cms.int32(7))
 
+hltPixelTrackRecHitsTable = cms.EDProducer("HLTTracksRecHitsTableProducer",
+                                            tableName = cms.string("hltPixelTrackRecHits"),
+                                            skipNonExistingSrc = cms.bool(True),
+                                            tracksSrc = cms.InputTag("hltPhase2PixelTracks"),
+                                            maxRecHits = cms.uint32(16),
+                                            precision = cms.int32(7)
+)
+
+
+phase2CAExtension.toModify(hltPixelTrackTable, src = "hltPhase2PixelTracksCAExtension")
+phase2CAExtension.toModify(hltPixelTrackExtTable, tracksSrc = "hltPhase2PixelTracksCAExtension")
+phase2CAExtension.toModify(hltPixelTrackRecHitsTable, tracksSrc = "hltPhase2PixelTracksCAExtension")
+
 hltGeneralTrackTable = cms.EDProducer(
     "SimpleTriggerTrackFlatTableProducer",
     skipNonExistingSrc = cms.bool(True),

PhysicsTools/NanoAOD/plugins/BuildFile.xml

@@ -18,6 +18,7 @@
 <use name="DataFormats/OnlineMetaData"/>
 <use name="DataFormats/PatCandidates"/>
 <use name="DataFormats/ProtonReco"/>
+<use name="DataFormats/TrackReco"/>
 <use name="FWCore/Framework"/>
 <use name="FWCore/ParameterSet"/>
 <use name="FWCore/ServiceRegistry"/>
@@ -36,6 +37,9 @@
   <use name="RecoJets/JetAlgorithms"/>
   <use name="TrackingTools/Records"/>
   <use name="roottmva"/>
+  <use name="SimDataFormats/Associations"/>
+  <use name="SimDataFormats/TrackingAnalysis"/>
+  <use name="SimTracker/Common"/>
   <flags EDM_PLUGIN="1"/>
 </library>
 <!-- Disable RNuple plugin: https://github.com/cms-sw/cmsdist/pull/9451#issuecomment-2393861542