Add vector of tracks to TICL tracksters and candidates and add other validation plots

DataFormats/HGCalReco/interface/TICLCandidate.h

@@ -2,6 +2,7 @@
 #define DataFormats_HGCalReco_TICLCandidate_h
 
 #include "DataFormats/Candidate/interface/LeafCandidate.h"
+#include "DataFormats/Common/interface/Ptr.h"
 #include "DataFormats/Common/interface/Ref.h"
 #include "DataFormats/HGCalReco/interface/Trackster.h"
 #include "DataFormats/Math/interface/Point3D.h"
@@ -33,18 +34,19 @@ class TICLCandidate : public reco::LeafCandidate {
         rawEnergy_(0.f) {}
 
   TICLCandidate(const edm::Ptr<reco::Track> trackPtr, const edm::Ptr<ticl::Trackster>& tracksterPtr)
-      : LeafCandidate(), tracksters_{}, trackPtr_(trackPtr), time_(0.f), timeError_(-1.f) {
-    if (trackPtr_.isNull() and tracksterPtr.isNull())
+      : LeafCandidate(), tracksters_{}, trackPtrs_{}, time_(0.f), timeError_(-1.f) {
+    if (trackPtr.isNull() and tracksterPtr.isNull())
       throw cms::Exception("NullPointerError")
           << "TICLCandidate constructor: at least one between track and trackster must be valid";
 
     if (tracksterPtr.isNonnull()) {
       tracksters_.push_back(tracksterPtr);
       auto const& trackster = tracksters_[0].get();
       idProbabilities_ = trackster->id_probabilities();
-      if (trackPtr_.isNonnull()) {
+      if (trackPtr.isNonnull()) {
+        trackPtrs_.push_back(trackPtr);
         auto pdgId = trackster->isHadronic() ? 211 : 11;
-        auto const& tk = trackPtr_.get();
+        auto const& tk = trackPtr.get();
         setPdgId(pdgId * tk->charge());
         setCharge(tk->charge());
         rawEnergy_ = trackster->raw_energy();
@@ -69,7 +71,8 @@ class TICLCandidate : public reco::LeafCandidate {
       }
     } else {
       //candidate from track only
-      auto const& tk = trackPtr_.get();
+      trackPtrs_.push_back(trackPtr);
+      auto const& tk = trackPtr.get();
       setPdgId(211 * tk->charge());
       setCharge(tk->charge());
       const float energy = std::sqrt(tk->p() * tk->p() + ticl::mpion2);
@@ -95,15 +98,18 @@ class TICLCandidate : public reco::LeafCandidate {
     MTDtimeError_ = timeError;
   };
 
-  inline const edm::Ptr<reco::Track> trackPtr() const { return trackPtr_; }
-  void setTrackPtr(const edm::Ptr<reco::Track>& trackPtr) { trackPtr_ = trackPtr; }
+  inline const edm::Ptr<reco::Track> trackPtr(int index = 0) const {
+    return trackPtrs_.empty() ? edm::Ptr<reco::Track>() : trackPtrs_[index];
+  }
+  inline const std::vector<edm::Ptr<reco::Track>> trackPtrs() const { return trackPtrs_; }
+  void addTrackPtr(const edm::Ptr<reco::Track>& trackPtr) { trackPtrs_.push_back(trackPtr); }
 
   inline float rawEnergy() const { return rawEnergy_; }
   void setRawEnergy(float rawEnergy) { rawEnergy_ = rawEnergy; }
 
-  inline const std::vector<edm::Ptr<ticl::Trackster> > tracksters() const { return tracksters_; };
+  inline const std::vector<edm::Ptr<ticl::Trackster>> tracksters() const { return tracksters_; };
 
-  void setTracksters(const std::vector<edm::Ptr<ticl::Trackster> >& tracksters) { tracksters_ = tracksters; }
+  void setTracksters(const std::vector<edm::Ptr<ticl::Trackster>>& tracksters) { tracksters_ = tracksters; }
   void addTrackster(const edm::Ptr<ticl::Trackster>& trackster) {
     tracksters_.push_back(trackster);
     time_ = trackster->time();
@@ -129,8 +135,8 @@ class TICLCandidate : public reco::LeafCandidate {
 private:
   // vector of Ptr so Tracksters can come from different collections
   // and there can be derived classes
-  std::vector<edm::Ptr<ticl::Trackster> > tracksters_;
-  edm::Ptr<reco::Track> trackPtr_;
+  std::vector<edm::Ptr<ticl::Trackster>> tracksters_;
+  std::vector<edm::Ptr<reco::Track>> trackPtrs_;
   // Since it contains multiple tracksters, duplicate the probability interface
   std::array<float, 8> idProbabilities_;
 

DataFormats/HGCalReco/interface/Trackster.h

@@ -77,8 +77,8 @@ namespace ticl {
     inline void setRawEmPt(float value) { raw_em_pt_ = value; }
     inline void calculateRawEmPt() { raw_em_pt_ = raw_em_energy_ / std::cosh(barycenter_.eta()); }
     inline void setBarycenter(Vector value) { barycenter_ = value; }
-    inline void setTrackIdx(int index) { track_idx_ = index; }
-    int trackIdx() const { return track_idx_; }
+    inline void addTrackIdx(int index) { track_idxs_.push_back(index); }
+    int trackIdx(int index = 0) const { return track_idxs_.empty() ? -1 : track_idxs_[index]; }
     inline bool isHadronic(float th = 0.5f) const {
       return id_probability(Trackster::ParticleType::photon) + id_probability(Trackster::ParticleType::electron) < th;
     }
@@ -198,6 +198,7 @@ namespace ticl {
     inline const std::array<float, 3> &sigmasPCA() const { return sigmasPCA_; }
     inline const std::array<float, 8> &id_probabilities() const { return id_probabilities_; }
     inline const float id_probabilities(int index) const { return id_probabilities_[index]; }
+    inline const std::vector<int> &trackIdxs() const { return track_idxs_; }
 
     // convenience method to return the ID probability for a certain particle type
     inline float id_probability(ParticleType type) const {
@@ -236,7 +237,7 @@ namespace ticl {
     // created the trackster. For track-based seeding the pointer to the track
     // can be cooked using the previous ProductID and this index.
     int seedIndex_;
-    int track_idx_ = -1;
+    std::vector<int> track_idxs_;
 
     std::array<Vector, 3> eigenvectors_;
     std::array<float, 3> eigenvalues_;

DataFormats/HGCalReco/src/classes_def.xml

@@ -1,6 +1,7 @@
 
 <lcgdict>
-  <class name="ticl::Trackster" ClassVersion="11">
+  <class name="ticl::Trackster" ClassVersion="12">
+    <version ClassVersion="12" checksum="912917580"/>
     <version ClassVersion="11" checksum="2111726404"/>
     <version ClassVersion="10" checksum="556627704"/>
     <version ClassVersion="9" checksum="1001808235"/>
@@ -66,7 +67,8 @@
   <class name="edm::Wrapper<TICLSeedingRegion>" />
   <class name="edm::Wrapper<std::vector<TICLSeedingRegion> >" />
 
-  <class name="TICLCandidate" ClassVersion="4">
+  <class name="TICLCandidate" ClassVersion="5">
+   <version ClassVersion="5" checksum="809321696"/>
    <version ClassVersion="4" checksum="2260471800"/>
    <version ClassVersion="3" checksum="450468662"/>
   </class>

RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc

@@ -396,7 +396,7 @@ void LinkingAlgoByDirectionGeometric::linkTracksters(const edm::Handle<std::vect
   for (unsigned &i : candidateTrackIds) {
     if (tsNearTk[i].empty() && tsNearTkAtInt[i].empty()) {  // nothing linked to track, make charged hadrons
       TICLCandidate chargedHad;
-      chargedHad.setTrackPtr(edm::Ptr<reco::Track>(tkH, i));
+      chargedHad.addTrackPtr(edm::Ptr<reco::Track>(tkH, i));
       chargedHadronsFromTk.push_back(chargedHad);
       continue;
     }
@@ -506,11 +506,11 @@ void LinkingAlgoByDirectionGeometric::linkTracksters(const edm::Handle<std::vect
     // do not create a candidate if no tracksters were added to candidate
     // can happen if all the tracksters linked to that track were already masked
     if (!chargedCandidate.tracksters().empty()) {
-      chargedCandidate.setTrackPtr(edm::Ptr<reco::Track>(tkH, i));
+      chargedCandidate.addTrackPtr(edm::Ptr<reco::Track>(tkH, i));
       chargedCandidates.push_back(chargedCandidate);
     } else {  // create charged hadron
       TICLCandidate chargedHad;
-      chargedHad.setTrackPtr(edm::Ptr<reco::Track>(tkH, i));
+      chargedHad.addTrackPtr(edm::Ptr<reco::Track>(tkH, i));
       chargedHadronsFromTk.push_back(chargedHad);
     }
   }

RecoHGCal/TICL/plugins/SimTrackstersProducer.cc

@@ -435,9 +435,8 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
 
   makePUTrackster(inputClusterMask, *output_mask, *resultPU, caloParticles_h.id(), 0);
 
-  auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::pair<int, float> {
-    int trackIdx = -1;
-    float quality = 0.f;
+  auto simTrackToRecoTrack = [&](UniqueSimTrackId simTkId) -> std::vector<int> {
+    std::vector<int> trackIdx;
     auto ipos = simTrackToTPMap.mapping.find(simTkId);
     if (ipos != simTrackToTPMap.mapping.end()) {
       auto jpos = TPtoRecoTrackMap.find((ipos->second));
@@ -446,41 +445,55 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
         if (!associatedRecoTracks.empty()) {
           // associated reco tracks are sorted by decreasing quality
           if (associatedRecoTracks[0].second > qualityCutTrack_) {
-            trackIdx = &(*associatedRecoTracks[0].first) - &recoTracks[0];
-            quality = associatedRecoTracks[0].second;
+            trackIdx.push_back(&(*associatedRecoTracks[0].first) - &recoTracks[0]);
+          }
+        }
+      }
+      const auto& tp = (*ipos->second);
+      if (!tp.decayVertices().empty()) {
+        const auto& iTV = tp.decayVertices()[0];
+        for (auto iTP = iTV->daughterTracks_begin(); iTP != iTV->daughterTracks_end(); ++iTP) {
+          auto kpos = TPtoRecoTrackMap.find((*iTP));
+          if (kpos != TPtoRecoTrackMap.end()) {
+            auto& associatedRecoTracks = kpos->val;
+            if (!associatedRecoTracks.empty()) {
+              // associated reco tracks are sorted by decreasing quality
+              if (associatedRecoTracks[0].second > qualityCutTrack_) {
+                trackIdx.push_back(&(*associatedRecoTracks[0].first) - &recoTracks[0]);
+              }
+            }
           }
         }
       }
     }
-    return {trackIdx, quality};
+    return trackIdx;
   };
 
-  // Creating the map from TrackingParticle to SimTrackstersFromCP
+  // Set the reco track id to SimTrackstersFromCP
   auto& simTrackstersFromCP = *result_fromCP;
   for (unsigned int i = 0; i < simTrackstersFromCP.size(); ++i) {
     if (simTrackstersFromCP[i].vertices().empty())
       continue;
     const auto& simTrack = caloparticles[simTrackstersFromCP[i].seedIndex()].g4Tracks()[0];
     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
-    auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
-    if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
-      auto trackIndex = bestAssociatedRecoTrack.first;
-      simTrackstersFromCP[i].setTrackIdx(trackIndex);
+    auto bestAssociatedRecoTracks = simTrackToRecoTrack(simTkIds);
+    if (not bestAssociatedRecoTracks.empty()) {
+      for (auto const trackIndex : bestAssociatedRecoTracks)
+        simTrackstersFromCP[i].addTrackIdx(trackIndex);
     }
   }
 
   auto& simTracksters = *result;
-  // Creating the map from TrackingParticle to SimTrackster
-  std::unordered_map<unsigned int, std::vector<unsigned int>> TPtoSimTracksterMap;
+  // Set the reco track id to simTracksters
   for (unsigned int i = 0; i < simTracksters.size(); ++i) {
     const auto& simTrack = (simTracksters[i].seedID() == caloParticles_h.id())
                                ? caloparticles[simTracksters[i].seedIndex()].g4Tracks()[0]
                                : simclusters[simTracksters[i].seedIndex()].g4Tracks()[0];
     UniqueSimTrackId simTkIds(simTrack.trackId(), simTrack.eventId());
-    auto bestAssociatedRecoTrack = simTrackToRecoTrack(simTkIds);
-    if (bestAssociatedRecoTrack.first != -1 and bestAssociatedRecoTrack.second > qualityCutTrack_) {
-      auto trackIndex = bestAssociatedRecoTrack.first;
-      simTracksters[i].setTrackIdx(trackIndex);
+    auto bestAssociatedRecoTracks = simTrackToRecoTrack(simTkIds);
+    if (not bestAssociatedRecoTracks.empty()) {
+      for (auto const trackIndex : bestAssociatedRecoTracks)
+        simTracksters[i].addTrackIdx(trackIndex);
     }
   }
 
@@ -502,12 +515,15 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
                        : SimClusterToCaloParticleMap[simTrackster.seedIndex()];
     auto const& tCP = (*result_fromCP)[cp_index];
     if (!tCP.vertices().empty()) {
-      auto trackIndex = tCP.trackIdx();
+      auto trackIndices = tCP.trackIdxs();
 
       auto& cand = (*result_ticlCandidates)[cp_index];
       cand.addTrackster(edm::Ptr<Trackster>(simTracksters_h, i));
-      if (trackIndex != -1 && caloparticles[cp_index].charge() != 0)
-        cand.setTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
+      if (cand.trackPtrs().empty() and not trackIndices.empty() and caloparticles[cp_index].charge() != 0) {
+        for (const auto trackIndex : trackIndices) {
+          cand.addTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
+        }
+      }
       toKeep.push_back(cp_index);
     }
   }

RecoHGCal/TICL/plugins/TICLCandidateProducer.cc

@@ -434,15 +434,16 @@ void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &r
     float invTimeErr = 0.f;
     float timeErr = -1.f;
 
+    const int trackIndex =
+        cand.trackPtr().isNonnull() ? (cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get()) : -1;
     for (const auto &tr : cand.tracksters()) {
       if (tr->timeError() > 0) {
         const auto invTimeESq = pow(tr->timeError(), -2);
         const auto x = tr->barycenter().X();
         const auto y = tr->barycenter().Y();
         const auto z = tr->barycenter().Z();
         auto path = std::sqrt(x * x + y * y + z * z);
-        if (cand.trackPtr().get() != nullptr) {
-          const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
+        if (trackIndex != -1) {
           if (useMTDTiming_ and inputTimingView.timeErr()[trackIndex] > 0) {
             const auto xMtd = inputTimingView.posInMTD_x()[trackIndex];
             const auto yMtd = inputTimingView.posInMTD_y()[trackIndex];
@@ -473,7 +474,6 @@ void TICLCandidateProducer::assignTimeToCandidates(std::vector<TICLCandidate> &r
 
     if (useMTDTiming_ and cand.charge()) {
       // Check MTD timing availability
-      const auto &trackIndex = cand.trackPtr().get() - (edm::Ptr<reco::Track>(track_h, 0)).get();
       const bool assocQuality = inputTimingView.MVAquality()[trackIndex] > timingQualityThreshold_;
       if (assocQuality) {
         const auto timeHGC = cand.time();