Merge pull request #48935 from cms-ngt-hlt/ngtScoutingTrackExtras

[NGT] Add extra NGT scouting track table producer, when using `NANO:@NGTScoutingVal`

Author: cmsbuild

HLTrigger/NGTScouting/plugins/HLTTracksExtraTableProducer.cc

@@ -0,0 +1,103 @@
+#include <memory>
+
+// user include files
+#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/MakerMacros.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+//
+// class declaration
+//
+
+class HLTTracksExtraTableProducer : public edm::stream::EDProducer<> {
+public:
+  explicit HLTTracksExtraTableProducer(const edm::ParameterSet&);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  // ----------member data ---------------------------
+  const bool skipNonExistingSrc_;
+  const std::string tableName_;
+  const unsigned int precision_;
+  const edm::EDGetTokenT<std::vector<reco::Track>> tracks_;
+  const edm::EDGetTokenT<reco::BeamSpot> beamSpot_;
+};
+
+//
+// constructors
+//
+
+HLTTracksExtraTableProducer::HLTTracksExtraTableProducer(const edm::ParameterSet& params)
+    : skipNonExistingSrc_(params.getParameter<bool>("skipNonExistingSrc")),
+      tableName_(params.getParameter<std::string>("tableName")),
+      precision_(params.getParameter<int>("precision")),
+      tracks_(consumes<std::vector<reco::Track>>(params.getParameter<edm::InputTag>("tracksSrc"))),
+      beamSpot_(consumes<reco::BeamSpot>(params.getParameter<edm::InputTag>("beamSpot"))) {
+  produces<nanoaod::FlatTable>(tableName_);
+}
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data  ------------
+void HLTTracksExtraTableProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  using namespace edm;
+
+  //vertex collection
+  auto tracksIn = iEvent.getHandle(tracks_);
+  auto beamSpotIn = iEvent.getHandle(beamSpot_);
+  const size_t nTracks = tracksIn.isValid() ? (*tracksIn).size() : 0;
+
+  static constexpr float default_value = std::numeric_limits<float>::quiet_NaN();
+
+  // initialize to quiet Nans
+  std::vector<float> v_dxyBS(nTracks, default_value);
+  std::vector<float> v_dzBS(nTracks, default_value);
+
+  if ((tracksIn.isValid() && beamSpotIn.isValid()) || !(this->skipNonExistingSrc_)) {
+    const auto& tracks = *tracksIn;
+    const auto& beamSpot = *beamSpotIn;
+    math::XYZPoint point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
+    for (size_t tk_index = 0; tk_index < nTracks; ++tk_index) {
+      v_dxyBS[tk_index] = tracks[tk_index].dxy(point);
+      v_dzBS[tk_index] = tracks[tk_index].dz(point);
+    }
+  } else {
+    edm::LogWarning("HLTTracksExtraTableProducer")
+        << " Invalid handle for " << tableName_ << " in tracks input collection";
+  }
+
+  //table for all primary vertices
+  auto tracksTable = std::make_unique<nanoaod::FlatTable>(nTracks, tableName_, /*singleton*/ false, /*extension*/ true);
+  tracksTable->addColumn<float>("dzBS", v_dzBS, "tracks dz() w.r.t Beam Spot", precision_);
+  tracksTable->addColumn<float>("dxyBS", v_dxyBS, "tracks dxy() w.r.t Beam Spot", precision_);
+  iEvent.put(std::move(tracksTable), tableName_);
+}
+
+// ------------ fill 'descriptions' with the allowed parameters for the module ------------
+void HLTTracksExtraTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+
+  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 collections");
+  desc.add<edm::InputTag>("beamSpot")->setComment("input BeamSpot collection");
+  desc.add<int>("precision", 7);
+  descriptions.addWithDefaultLabel(desc);
+}
+
+// ------------ define this as a plug-in ------------
+DEFINE_FWK_MODULE(HLTTracksExtraTableProducer);

HLTrigger/NGTScouting/plugins/HLTVertexTableProducer.cc

@@ -2,7 +2,6 @@
 
 // user include files
 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
-#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
 #include "DataFormats/Common/interface/ValueMap.h"
 #include "DataFormats/NanoAOD/interface/FlatTable.h"
 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
@@ -37,6 +36,7 @@ class HLTVertexTableProducer : public edm::stream::EDProducer<> {
   const edm::EDGetTokenT<edm::ValueMap<float>> pvsScore_;
   const StringCutObjectSelector<reco::Vertex> goodPvCut_;
   const std::string goodPvCutString_;
+  const bool usePF_;
   const std::string pvName_;
   const double dlenMin_, dlenSigMin_;
 };
@@ -52,11 +52,11 @@ HLTVertexTableProducer::HLTVertexTableProducer(const edm::ParameterSet& params)
       pvsScore_(consumes<edm::ValueMap<float>>(params.getParameter<edm::InputTag>("pvSrc"))),
       goodPvCut_(params.getParameter<std::string>("goodPvCut"), true),
       goodPvCutString_(params.getParameter<std::string>("goodPvCut")),
+      usePF_(params.getParameter<bool>("usePF")),
       pvName_(params.getParameter<std::string>("pvName")),
       dlenMin_(params.getParameter<double>("dlenMin")),
       dlenSigMin_(params.getParameter<double>("dlenSigMin")) {
   produces<nanoaod::FlatTable>("PV");
-  produces<edm::PtrVector<reco::VertexCompositePtrCandidate>>();
 }
 
 //
@@ -69,6 +69,7 @@ void HLTVertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup&
 
   //vertex collection
   auto pvsIn = iEvent.getHandle(pvs_);
+  auto pvScoreIn = iEvent.getHandle(pvsScore_);
   const size_t nPVs = pvsIn.isValid() ? (*pvsIn).size() : 0;
 
   static constexpr float default_value = std::numeric_limits<float>::quiet_NaN();
@@ -90,7 +91,6 @@ void HLTVertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup&
 
   if (pvsIn.isValid() || !(this->skipNonExistingSrc_)) {
     const auto& pvs = *pvsIn;
-    const auto& pvsScoreProd = iEvent.get(pvsScore_);
 
     auto pfcIn = iEvent.getHandle(pfc_);
     const bool isPfcValid = pfcIn.isValid();
@@ -109,40 +109,55 @@ void HLTVertexTableProducer::produce(edm::Event& iEvent, const edm::EventSetup&
       v_zError[i] = pv.zError();
       v_nTracks[i] = pv.nTracks();
       v_is_good[i] = goodPvCut_(pv);
-      v_pv_score[i] = pvsScoreProd.get(pvsIn.id(), i);
 
-      float sumpt2 = 0.f, sumpx = 0.f, sumpy = 0.f;
+      if (pvScoreIn.isValid() || !(this->skipNonExistingSrc_)) {
+        const auto& pvsScoreProd = *pvScoreIn;
+        v_pv_score[i] = pvsScoreProd.get(pvsIn.id(), i);
+      }
 
-      if (isPfcValid || !(this->skipNonExistingSrc_)) {
-        for (const auto& obj : *pfcIn) {
-          if (obj.charge() == 0 || !obj.trackRef().isNonnull())
-            continue;
+      float sumpt2 = 0.f, sumpx = 0.f, sumpy = 0.f;
 
-          const auto dz = std::abs(obj.trackRef()->dz(pos));
-          if (dz >= 0.2)
-            continue;
+      if (usePF_) {
+        if (isPfcValid || !(this->skipNonExistingSrc_)) {
+          for (const auto& obj : *pfcIn) {
+            if (obj.charge() == 0 || !obj.trackRef().isNonnull())
+              continue;
 
-          bool isClosest = true;
-          for (size_t j = 0; j < nPVs; ++j) {
-            if (j == i)
+            const auto dz = std::abs(obj.trackRef()->dz(pos));
+            if (dz >= 0.2)
               continue;
-            const auto dz_j = std::abs(obj.trackRef()->dz(pvs[j].position()));
-            if (dz_j < dz) {
-              isClosest = false;
-              break;
+
+            bool isClosest = true;
+            for (size_t j = 0; j < nPVs; ++j) {
+              if (j == i)
+                continue;
+              const auto dz_j = std::abs(obj.trackRef()->dz(pvs[j].position()));
+              if (dz_j < dz) {
+                isClosest = false;
+                break;
+              }
             }
-          }
 
-          if (isClosest) {
-            const float pt = obj.pt();
-            sumpt2 += pt * pt;
-            sumpx += obj.px();
-            sumpy += obj.py();
+            if (isClosest) {
+              const float pt = obj.pt();
+              sumpt2 += pt * pt;
+              sumpx += obj.px();
+              sumpy += obj.py();
+            }
           }
+        } else {
+          edm::LogWarning("HLTVertexTableProducer")
+              << " Invalid handle for " << pvName_ << " in PF candidate input collection";
         }
       } else {
-        edm::LogWarning("HLTVertexTableProducer")
-            << " Invalid handle for " << pvName_ << " in PF candidate input collection";
+        // Loop over tracks used in PV fit
+        for (auto t = pv.tracks_begin(); t != pv.tracks_end(); ++t) {
+          const auto& trk = **t;  // trk is a reco::TrackBase
+          const float pt = trk.pt();
+          sumpt2 += pt * pt;
+          sumpx += trk.px();
+          sumpy += trk.py();
+        }
       }
       v_pv_sumpt2[i] = sumpt2;
       v_pv_sumpx[i] = sumpx;
@@ -186,6 +201,8 @@ void HLTVertexTableProducer::fillDescriptions(edm::ConfigurationDescriptions& de
   desc.add<std::string>("pvName")->setComment("name of the flat table ouput");
   desc.add<edm::InputTag>("pvSrc")->setComment(
       "std::vector<reco::Vertex> and ValueMap<float> primary vertex input collections");
+  desc.add<bool>("usePF", true)
+      ->setComment("if true, use PF candidate-based association; if false, use only tracks used in PV fit");
   desc.add<edm::InputTag>("pfSrc")->setComment("reco::PFCandidateCollection PF candidates input collections");
   desc.add<std::string>("goodPvCut")->setComment("selection on the primary vertex");
   desc.add<double>("dlenMin")->setComment("minimum value of dl to select secondary vertex");

HLTrigger/NGTScouting/python/HLTNanoProducer_cff.py

@@ -47,6 +47,7 @@
     #+ hltTriggerAcceptFilter
     + hltVertexTable
     + hltPixelTrackTable
+    + hltPixelVertexTable
     + hltGeneralTrackTable
     + hltEgammaPacker
     + hltPhotonTable
@@ -70,6 +71,7 @@
     + dstTriggerAcceptFilter
     + hltVertexTable
     + hltPixelTrackTable
+    + hltPixelVertexTable
     + hltGeneralTrackTable
     + hltEgammaPacker
     + hltPhotonTable
@@ -88,6 +90,11 @@
     + HTTable
 )
 
+trackingExtraNanoProducer = cms.Sequence(
+    hltPixelTrackExtTable+
+    hltGeneralTrackExtTable
+)
+
 def hltNanoCustomize(process):
 
     if hasattr(process, "NANOAODSIMoutput"):
@@ -108,6 +115,10 @@ def hltNanoCustomize(process):
 def hltNanoValCustomize(process):
     if hasattr(process, "dstNanoProducer"):
 
-        process.dstNanoProducer += (process.hltTiclAssociationsTableSequence + process.hltSimTracksterSequence + process.hltSimTiclCandidateTable + process.hltSimTiclCandidateExtraTable )
+        process.dstNanoProducer += (process.hltTiclAssociationsTableSequence +
+                                    process.hltSimTracksterSequence +
+                                    process.hltSimTiclCandidateTable +
+                                    process.hltSimTiclCandidateExtraTable +
+                                    process.trackingExtraNanoProducer)
 
     return process

HLTrigger/NGTScouting/python/hltTracks_cfi.py

@@ -15,6 +15,8 @@
         phi = Var("phi()", "float", doc = "#phi (rad)"),
         dXY = Var("dxy()", "float", doc = "dXY (cm)"),
         dZ = Var("dz()", "float", doc = "dZ (cm)"),
+        dxyError = Var("dxyError()", "float", doc = "dxyError (cm)"),
+        dZError = Var("dzError()", "float", doc = "dzError (cm)"),
         t0 = Var("t0()", "float", doc = "t0 (ns)"),
         vx = Var("vx()", "float", doc = "vx (cm)"),
         vy = Var("vy()", "float", doc = "vy (cm)"),
@@ -30,6 +32,13 @@
     )
 )
 
+hltPixelTrackExtTable = cms.EDProducer("HLTTracksExtraTableProducer",
+                                       tableName = cms.string("hltPixelTrack"),                                    
+                                       skipNonExistingSrc = cms.bool(True),
+                                       tracksSrc = cms.InputTag("hltPhase2PixelTracks"),
+                                       beamSpot = cms.InputTag("hltOnlineBeamSpot"),
+                                       precision = cms.int32(7))
+
 hltGeneralTrackTable = cms.EDProducer(
     "SimpleTriggerTrackFlatTableProducer",
     skipNonExistingSrc = cms.bool(True),
@@ -55,3 +64,10 @@
         ndof = Var("ndof()", "float", doc = "Number of degrees of freedom"),
     )
 )
+
+hltGeneralTrackExtTable = cms.EDProducer("HLTTracksExtraTableProducer",
+                                         tableName = cms.string("hltGeneralTrack"),
+                                         skipNonExistingSrc = cms.bool(True),
+                                         tracksSrc = cms.InputTag("hltGeneralTracks"),
+                                         beamSpot = cms.InputTag("hltOnlineBeamSpot"),
+                                         precision = cms.int32(7))

HLTrigger/NGTScouting/python/hltVertices_cfi.py

@@ -8,5 +8,14 @@
                                 pfSrc = cms.InputTag("hltParticleFlowTmp"),
                                 dlenMin = cms.double(0),
                                 dlenSigMin = cms.double(3),
-                                pvName = cms.string("hltPrimaryVertex"),
-                                )
+                                pvName = cms.string("hltPrimaryVertex"))
+
+hltPixelVertexTable = cms.EDProducer("HLTVertexTableProducer",
+                                     skipNonExistingSrc = cms.bool(True),
+                                     pvSrc = cms.InputTag("hltPhase2PixelVertices"),
+                                     goodPvCut = cms.string(""),
+                                     usePF = cms.bool(False), # use directly the tracks from PV fit 
+                                     pfSrc = cms.InputTag(""),
+                                     dlenMin = cms.double(0),
+                                     dlenSigMin = cms.double(3),
+                                     pvName = cms.string("hltPixelVertex"))