GeneralInterpretation (cms-sw#19)

* ticlv5: optimize and add new constructor and method to TICLCandidate * ticlv5: add isHadronic method to Trackster dataformat * ticlv5: add Input struct for timing information * ticlv5: implement general interpretation algo for linking tracks to tracksters --------- Co-authored-by: AuroraPerego <[email protected]>
waredjeb · May 30, 2024 · 0070428 · 0070428
1 parent 0d80d5e
commit 0070428
Show file tree

Hide file tree

Showing 28 changed files with 1,074 additions and 213 deletions.
diff --git a/CommonTools/RecoAlgos/interface/MultiVectorManager.h b/CommonTools/RecoAlgos/interface/MultiVectorManager.h
@@ -10,84 +10,70 @@
 template <typename T>
 class MultiVectorManager {
 public:
-    void addVector(const std::vector<T>& vec) {
-        vectors.emplace_back(vec.begin(), vec.end());
-        offsets.push_back(totalSize);
-        totalSize += vec.size();
-    }
+  void addVector(const std::vector<T>& vec) {
+    vectors.emplace_back(vec.begin(), vec.end());
+    offsets.push_back(totalSize);
+    totalSize += vec.size();
+  }
 
-    T& operator[](size_t globalIndex) {
-        return const_cast<T&>(static_cast<const MultiVectorManager*>(this)->operator[](globalIndex));
-    }
+  T& operator[](size_t globalIndex) {
+    return const_cast<T&>(static_cast<const MultiVectorManager*>(this)->operator[](globalIndex));
+  }
 
-    const T& operator[](size_t globalIndex) const {
-        assert(globalIndex < totalSize && "Global index out of range");
+  const T& operator[](size_t globalIndex) const {
+    assert(globalIndex < totalSize && "Global index out of range");
 
-        auto it = std::upper_bound(offsets.begin(), offsets.end(), globalIndex);
-        size_t vectorIndex = std::distance(offsets.begin(), it) - 1;
-        size_t localIndex = globalIndex - offsets[vectorIndex];
+    auto it = std::upper_bound(offsets.begin(), offsets.end(), globalIndex);
+    size_t vectorIndex = std::distance(offsets.begin(), it) - 1;
+    size_t localIndex = globalIndex - offsets[vectorIndex];
 
-        return vectors[vectorIndex][localIndex];
-    }
+    return vectors[vectorIndex][localIndex];
+  }
 
-    size_t getGlobalIndex(size_t vectorIndex, size_t localIndex) const {
-        assert(vectorIndex < vectors.size() && "Vector index out of range");
+  size_t getGlobalIndex(size_t vectorIndex, size_t localIndex) const {
+    assert(vectorIndex < vectors.size() && "Vector index out of range");
 
-        const auto& vec = vectors[vectorIndex];
-        assert(localIndex < vec.size() && "Local index out of range");
+    const auto& vec = vectors[vectorIndex];
+    assert(localIndex < vec.size() && "Local index out of range");
 
-        return offsets[vectorIndex] + localIndex;
-    }
+    return offsets[vectorIndex] + localIndex;
+  }
 
-    size_t size() const {
-        return totalSize;
-    }
+  size_t size() const { return totalSize; }
 
-    std::pair<size_t, size_t> getVectorAndLocalIndex(size_t globalIndex) const {
-        assert(globalIndex < totalSize && "Global index out of range");
+  std::pair<size_t, size_t> getVectorAndLocalIndex(size_t globalIndex) const {
+    assert(globalIndex < totalSize && "Global index out of range");
 
-        auto it = std::upper_bound(offsets.begin(), offsets.end(), globalIndex);
-        size_t vectorIndex = std::distance(offsets.begin(), it) - 1;
-        size_t localIndex = globalIndex - offsets[vectorIndex];
+    auto it = std::upper_bound(offsets.begin(), offsets.end(), globalIndex);
+    size_t vectorIndex = std::distance(offsets.begin(), it) - 1;
+    size_t localIndex = globalIndex - offsets[vectorIndex];
 
-        return {vectorIndex, localIndex};
-    }
+    return {vectorIndex, localIndex};
+  }
 
-    class Iterator {
-    public:
-        Iterator(const MultiVectorManager& manager, size_t index)
-            : manager(manager), currentIndex(index) {}
+  class Iterator {
+  public:
+    Iterator(const MultiVectorManager& manager, size_t index) : manager(manager), currentIndex(index) {}
 
-        bool operator!=(const Iterator& other) const {
-            return currentIndex != other.currentIndex;
-        }
+    bool operator!=(const Iterator& other) const { return currentIndex != other.currentIndex; }
 
-        T& operator*() const {
-            return const_cast<T&>(manager[currentIndex]);
-        }
+    T& operator*() const { return const_cast<T&>(manager[currentIndex]); }
 
-        void operator++() {
-            ++currentIndex;
-        }
+    void operator++() { ++currentIndex; }
 
-    private:
-        const MultiVectorManager& manager;
-        size_t currentIndex;
-    };
+  private:
+    const MultiVectorManager& manager;
+    size_t currentIndex;
+  };
 
-    Iterator begin() const {
-        return Iterator(*this, 0);
-    }
+  Iterator begin() const { return Iterator(*this, 0); }
 
-    Iterator end() const {
-        return Iterator(*this, totalSize);
-    }
+  Iterator end() const { return Iterator(*this, totalSize); }
 
 private:
-    std::vector<std::vector<T>> vectors;
-    std::vector<size_t> offsets;
-    size_t totalSize = 0;
+  std::vector<std::vector<T>> vectors;
+  std::vector<size_t> offsets;
+  size_t totalSize = 0;
 };
 
-
-#endif
+#endif
diff --git a/DataFormats/HGCalReco/interface/TICLCandidate.h b/DataFormats/HGCalReco/interface/TICLCandidate.h
@@ -7,6 +7,7 @@
 #include "DataFormats/Math/interface/Point3D.h"
 #include "DataFormats/TrackReco/interface/Track.h"
 #include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "RecoHGCal/TICL/interface/commons.h"
 
 // A TICLCandidate is a lightweight physics object made from one or multiple Tracksters.
 
@@ -15,17 +16,68 @@ class TICLCandidate : public reco::LeafCandidate {
   typedef ticl::Trackster::ParticleType ParticleType;
 
   TICLCandidate(Charge q, const LorentzVector& p4)
-      : LeafCandidate(q, p4), time_(0.f), timeError_(-1.f), rawEnergy_(0.f), idProbabilities_{} {}
+      : LeafCandidate(q, p4), idProbabilities_{}, time_(0.f), timeError_(-1.f), rawEnergy_(0.f) {}
 
-  TICLCandidate() : LeafCandidate(), time_(0.f), timeError_(-1.f), rawEnergy_(0.f), idProbabilities_{} {}
+  TICLCandidate() : LeafCandidate(), idProbabilities_{}, time_(0.f), timeError_(-1.f), rawEnergy_(0.f) {}
 
   TICLCandidate(const edm::Ptr<ticl::Trackster>& trackster)
       : LeafCandidate(),
+        idProbabilities_{},
+        tracksters_({trackster}),
         time_(trackster->time()),
         timeError_(trackster->timeError()),
-        rawEnergy_(0.f),
-        tracksters_({trackster}),
-        idProbabilities_{} {}
+        rawEnergy_(0.f) {}
+
+  TICLCandidate(const edm::Ptr<reco::Track> trackPtr, edm::Ptr<ticl::Trackster>& tracksterPtr)
+      : LeafCandidate(),
+        tracksters_({std::move(tracksterPtr)}),
+        trackPtr_(std::move(trackPtr)),
+        time_(0.f),
+        timeError_(-1.f) {
+    //TODO: Raise Error
+    assert(trackPtr_.isNonnull() or tracksters_[0].isNonnull());
+
+    if (tracksters_[0].isNonnull()) {
+      auto const& trackster = tracksters_[0].get();
+      idProbabilities_ = trackster->id_probabilities();
+      if (trackPtr_.isNonnull()) {
+        auto pdgId = trackster->isHadronic() ? 211 : 11;
+        auto const& tk = trackPtr_.get();
+        setPdgId(pdgId * tk->charge());
+        setCharge(tk->charge());
+        rawEnergy_ = trackster->raw_energy();
+        auto const& regrE = trackster->regressed_energy();
+        math::XYZTLorentzVector p4(regrE * tk->momentum().unit().x(),
+                                   regrE * tk->momentum().unit().y(),
+                                   regrE * tk->momentum().unit().z(),
+                                   regrE);
+        setP4(p4);
+
+      } else {
+        auto pdgId = trackster->isHadronic() ? 130 : 22;
+        setPdgId(pdgId);
+        setCharge(0);
+        rawEnergy_ = trackster->raw_energy();
+        const float& regrE = trackster->regressed_energy();
+        math::XYZTLorentzVector p4(regrE * trackster->barycenter().unit().x(),
+                                   regrE * trackster->barycenter().unit().y(),
+                                   regrE * trackster->barycenter().unit().z(),
+                                   regrE);
+        setP4(p4);
+      }
+    }
+
+    else {
+      //candidate from track only
+      auto const& tk = trackPtr_.get();
+      setPdgId(211 * tk->charge());
+      setCharge(tk->charge());
+      const float energy = std::sqrt(tk->p() * tk->p() + ticl::mpion2);
+      setRawEnergy(energy);
+      math::PtEtaPhiMLorentzVector p4Polar(tk->pt(), tk->eta(), tk->phi(), ticl::mpion);
+      setP4(p4Polar);
+    }
+  }
 
   inline float time() const { return time_; }
   inline float timeError() const { return timeError_; }
@@ -65,17 +117,17 @@ class TICLCandidate : public reco::LeafCandidate {
   inline void setIdProbability(ParticleType type, float value) { idProbabilities_[int(type)] = value; }
 
 private:
-  float time_;
-  float timeError_;
+  std::array<float, 8> idProbabilities_;
+  std::vector<edm::Ptr<ticl::Trackster> > tracksters_;
   edm::Ptr<reco::Track> trackPtr_;
 
+  float time_;
+  float timeError_;
   float rawEnergy_;
 
   // vector of Ptr so Tracksters can come from different collections
   // and there can be derived classes
-  std::vector<edm::Ptr<ticl::Trackster> > tracksters_;
 
   // Since it contains multiple tracksters, duplicate the probability interface
-  std::array<float, 8> idProbabilities_;
 };
 #endif
diff --git a/DataFormats/HGCalReco/interface/TICLLayerTile.h b/DataFormats/HGCalReco/interface/TICLLayerTile.h
@@ -12,9 +12,7 @@ class TICLLayerTileT {
 public:
   typedef T type;
 
-  void fill(float eta, float phi, unsigned int layerClusterId) {
-    tile_[globalBin(eta, phi)].push_back(layerClusterId);
-  }
+  void fill(float eta, float phi, unsigned int layerClusterId) { tile_[globalBin(eta, phi)].push_back(layerClusterId); }
 
   int etaBin(float eta) const {
     constexpr float etaRange = T::maxEta - T::minEta;

diff --git a/DataFormats/HGCalReco/interface/Trackster.h b/DataFormats/HGCalReco/interface/Trackster.h
@@ -74,20 +74,23 @@ namespace ticl {
     inline void setBarycenter(Vector value) { barycenter_ = value; }
     inline void setTrackIdx(int index) { track_idx_ = index; }
     int trackIdx() const { return track_idx_; }
-    inline void mergeTracksters(const Trackster& other){
+    inline bool isHadronic(float th = 0.5f) const {
+      return id_probability(Trackster::ParticleType::photon) + id_probability(Trackster::ParticleType::electron) < th;
+    }
+    inline void mergeTracksters(const Trackster &other) {
       *this += other;
 
-      //remove duplicates 
+      //remove duplicates
       removeDuplicates();
       zeroProbabilities();
     }
 
-    inline void mergeTracksters(const std::vector<Trackster>& others){
-      for(auto& other: others){
+    inline void mergeTracksters(const std::vector<Trackster> &others) {
+      for (auto &other : others) {
         *this += other;
       }
 
-      //remove duplicates 
+      //remove duplicates
       removeDuplicates();
       zeroProbabilities();
     }
@@ -222,14 +225,12 @@ namespace ticl {
           auto iDup = std::find(std::next(firstEl), vertices_.end(), lcIdx);
           while (iDup != vertices_.end()) {
             vertices_.erase(iDup);
-            vertex_multiplicity_.erase(vertex_multiplicity_.begin() +
-                                                     std::distance(std::begin(vertices_), iDup));
+            vertex_multiplicity_.erase(vertex_multiplicity_.begin() + std::distance(std::begin(vertices_), iDup));
             vertex_multiplicity_[firstPos] -= 1;
             iDup = std::find(std::next(firstEl), vertices_.end(), lcIdx);
           };
         }
-    }
-
+      }
     }
     inline void operator+=(const Trackster &other) {
       // use getters on other
@@ -239,9 +240,9 @@ namespace ticl {
       raw_em_pt_ += other.raw_em_pt();
       // add vertices and multiplicities
       std::copy(std::begin(other.vertices()), std::end(other.vertices()), std::back_inserter(vertices_));
-      std::copy(std::begin(other.vertex_multiplicity()), std::end(other.vertex_multiplicity()),
+      std::copy(std::begin(other.vertex_multiplicity()),
+                std::end(other.vertex_multiplicity()),
                 std::back_inserter(vertex_multiplicity_));
-
     }
   };
 

diff --git a/DataFormats/HGCalReco/src/classes_def.xml b/DataFormats/HGCalReco/src/classes_def.xml
@@ -51,7 +51,8 @@
   <class name="edm::Wrapper<TICLSeedingRegion>" />
   <class name="edm::Wrapper<std::vector<TICLSeedingRegion> >" />
 
-  <class name="TICLCandidate" ClassVersion="3">
+  <class name="TICLCandidate" ClassVersion="4">
+   <version ClassVersion="4" checksum="3155265136"/>
    <version ClassVersion="3" checksum="450468662"/>
   </class>
   <class name="std::vector<TICLCandidate>" />

diff --git a/DataFormats/TrackReco/src/classes_def.xml b/DataFormats/TrackReco/src/classes_def.xml
@@ -513,4 +513,7 @@
    <class name="std::vector<SeedStopInfo>" persistent="false"/>
    <class name="edm::Wrapper<std::vector<SeedStopInfo> >" persistent="false"/>
 
+   <class name="edm::ValueMap<Point3DBase<float,GlobalTag> >"/>
+   <class name="edm::Wrapper<edm::ValueMap<Point3DBase<float,GlobalTag> > >"/>
+
   </lcgdict>
diff --git a/Fireworks/Calo/plugins/FWCaloClusterProxyBuilder.cc b/Fireworks/Calo/plugins/FWCaloClusterProxyBuilder.cc
@@ -114,8 +114,8 @@ void FWCaloClusterProxyBuilder::build(const reco::CaloCluster &iData,
       continue;
 
     // HGCal
-    if (iData.algo() == reco::CaloCluster::hgcal_em || iData.algo() == reco::CaloCluster::hgcal_had || iData.algo() == reco::CaloCluster::hgcal_scintillator ||
-        (type >= 8 && type <= 10)) {
+    if (iData.algo() == reco::CaloCluster::hgcal_em || iData.algo() == reco::CaloCluster::hgcal_had ||
+        iData.algo() == reco::CaloCluster::hgcal_scintillator || (type >= 8 && type <= 10)) {
       if (heatmap && hitmap->find(it->first) == hitmap->end())
         continue;