Improve MC truth information in TICL + Associators + TICLDumper

- Add best associated RecoTrack to SimTracksters, by using SimTrack-TrackingParticle associator and TrackingParticle-RecoTrack associator. Add SimTrackster from PU by considering as PU all the remaining contribution from the LCs after the SimTracksterCP/SC building - Creates SimTICLCandidates, object represeting the best PF reconstruction possible with the available recoObjects. It contains all the SimTracksters beloning to the same CP and the recoTrack associated to the CP - Associators: compute score between SimTracksters and RecoTrackster - TICLDumper: EDAnalyzer to save TICL information in root file Co-authored-by: Wahid Redjeb <[email protected]> Co-authored-by: Aurora Perego <[email protected]>
cms-sw · Aug 9, 2023 · d20c3e1 · d20c3e1
1 parent 63ed0a8
commit d20c3e1
Show file tree

Hide file tree

Showing 34 changed files with 3,297 additions and 119 deletions.
diff --git a/DataFormats/HGCalReco/interface/TICLCandidate.h b/DataFormats/HGCalReco/interface/TICLCandidate.h
@@ -42,18 +42,23 @@ class TICLCandidate : public reco::LeafCandidate {
   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 addTrackster(const edm::Ptr<ticl::Trackster>& trackster) {
-    tracksters_.push_back(trackster);
-    time_ = trackster->time();
-    timeError_ = trackster->timeError();
-  }
+  void addTrackster(const edm::Ptr<ticl::Trackster>& trackster) { tracksters_.push_back(trackster); }
   // convenience method to return the ID probability for a certain particle type
   inline float id_probability(ParticleType type) const {
     // probabilities are stored in the same order as defined in the ParticleType enum
     return idProbabilities_[(int)type];
   }
 
+  inline const std::array<float, 8>& idProbabilities() const { return idProbabilities_; }
+
+  void zeroProbabilities() {
+    for (auto& p : idProbabilities_) {
+      p = 0.f;
+    }
+  }
+
   void setIdProbabilities(const std::array<float, 8>& idProbs) { idProbabilities_ = idProbs; }
+  inline void setIdProbability(ParticleType type, float value) { idProbabilities_[int(type)] = value; }
 
 private:
   float time_;

diff --git a/DataFormats/HGCalReco/interface/Trackster.h b/DataFormats/HGCalReco/interface/Trackster.h
@@ -18,7 +18,7 @@
 namespace ticl {
   class Trackster {
   public:
-    typedef math::XYZVector Vector;
+    typedef math::XYZVectorF Vector;
 
     enum IterationIndex { TRKEM = 0, EM, TRKHAD, HAD, MIP, SIM, SIM_CP };
 
@@ -37,21 +37,20 @@ namespace ticl {
     enum class PCAOrdering { ascending = 0, descending };
 
     Trackster()
-        : iterationIndex_(0),
-          seedIndex_(-1),
-          time_(0.f),
-          timeError_(-1.f),
+        : barycenter_({0.f, 0.f, 0.f}),
           regressed_energy_(0.f),
           raw_energy_(0.f),
+          time_(0.f),
+          timeError_(-1.f),
           raw_em_energy_(0.f),
+          id_probabilities_{},
           raw_pt_(0.f),
           raw_em_pt_(0.f),
-          barycenter_({0., 0., 0.}),
-          eigenvalues_{{0.f, 0.f, 0.f}},
-          sigmas_{{0.f, 0.f, 0.f}},
-          sigmasPCA_{{0.f, 0.f, 0.f}} {
-      zeroProbabilities();
-    }
+          seedIndex_(-1),
+          eigenvalues_{},
+          sigmas_{},
+          sigmasPCA_{},
+          iterationIndex_(0) {}
 
     inline void setIteration(const Trackster::IterationIndex i) { iterationIndex_ = i; }
     std::vector<unsigned int> &vertices() { return vertices_; }
@@ -73,6 +72,9 @@ namespace ticl {
     inline void setRawPt(float value) { raw_pt_ = value; }
     inline void setRawEmPt(float value) { raw_em_pt_ = value; }
     inline void setBarycenter(Vector value) { barycenter_ = value; }
+    inline void setTrackIdx(int index) { track_idx_ = index; }
+    int trackIdx() const { return track_idx_; }
+
     inline void fillPCAVariables(Eigen::Vector3d &eigenvalues,
                                  Eigen::Matrix3d &eigenvectors,
                                  Eigen::Vector3d &sigmas,
@@ -147,20 +149,23 @@ namespace ticl {
     }
 
   private:
-    // TICL iteration producing the trackster
-    uint8_t iterationIndex_;
+    Vector barycenter_;
+    float regressed_energy_;
+    float raw_energy_;
+    // -99, -1 if not available. ns units otherwise
+    float time_;
+    float timeError_;
+    float raw_em_energy_;
+
+    // trackster ID probabilities
+    std::array<float, 8> id_probabilities_;
 
     // The vertices of the DAG are the indices of the
     // 2d objects in the global collection
     std::vector<unsigned int> vertices_;
     std::vector<float> vertex_multiplicity_;
-
-    // The edges connect two vertices together in a directed doublet
-    // ATTENTION: order matters!
-    // A doublet generator should create edges in which:
-    // the first element is on the inner layer and
-    // the outer element is on the outer layer.
-    std::vector<std::array<unsigned int, 2> > edges_;
+    float raw_pt_;
+    float raw_em_pt_;
 
     // Product ID of the seeding collection used to create the Trackster.
     // For GlobalSeeding the ProductID is set to 0. For track-based seeding
@@ -173,31 +178,22 @@ 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;
 
-    // We also need the pointer to the original seeding region ??
-    // something like:
-    // int seedingRegionIdx;
-
-    // -99, -1 if not available. ns units otherwise
-    float time_;
-    float timeError_;
-
-    // regressed energy
-    float regressed_energy_;
-    float raw_energy_;
-    float raw_em_energy_;
-    float raw_pt_;
-    float raw_em_pt_;
-
-    // PCA Variables
-    Vector barycenter_;
-    std::array<float, 3> eigenvalues_;
     std::array<Vector, 3> eigenvectors_;
+    std::array<float, 3> eigenvalues_;
     std::array<float, 3> sigmas_;
     std::array<float, 3> sigmasPCA_;
 
-    // trackster ID probabilities
-    std::array<float, 8> id_probabilities_;
+    // The edges connect two vertices together in a directed doublet
+    // ATTENTION: order matters!
+    // A doublet generator should create edges in which:
+    // the first element is on the inner layer and
+    // the outer element is on the outer layer.
+    std::vector<std::array<unsigned int, 2> > edges_;
+
+    // TICL iteration producing the trackster
+    uint8_t iterationIndex_;
   };
 
   typedef std::vector<Trackster> TracksterCollection;

diff --git a/DataFormats/HGCalReco/src/classes_def.xml b/DataFormats/HGCalReco/src/classes_def.xml
@@ -1,6 +1,7 @@
 
 <lcgdict>
-  <class name="ticl::Trackster" ClassVersion="9">
+  <class name="ticl::Trackster" ClassVersion="10">
+    <version ClassVersion="10" checksum="556627704"/>
     <version ClassVersion="9" checksum="1001808235"/>
     <version ClassVersion="8" checksum="749412802"/>
     <version ClassVersion="7" checksum="1072014319"/>

diff --git a/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py b/RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
@@ -19,7 +19,11 @@
        'keep *_ticlTrackstersHFNoseMIP_*_*',
        'keep *_ticlTrackstersHFNoseHAD_*_*',
        'keep *_ticlTrackstersHFNoseMerge_*_*',] +
-      ['keep *_pfTICL_*_*']
+      ['keep *_pfTICL_*_*'] +
+      ['keep *_layerClusterSimClusterAssociationProducer_*_*','keep *_layerClusterCaloParticleAssociationProducer_*_*', 'keep *_layerClusterSimTracksterAssociationProducer_*_*'] + 
+      ['keep *_tracksterSimTracksterAssociationLinking_*_*' ,'keep *_tracksterSimTracksterAssociationPR_*_*'] +
+      ['keep *_tracksterSimTracksterAssociationLinkingPU_*_*' ,'keep *_tracksterSimTracksterAssociationPRPU_*_*'] +
+      ['keep *_tracksterSimTracksterAssociationLinkingbyCLUE3D_*_*', 'keep *_tracksterSimTracksterAssociationPRbyCLUE3D_*_*'] 
       )
     )
 TICL_RECO.outputCommands.extend(TICL_AOD.outputCommands)
@@ -28,6 +32,7 @@
 TICL_FEVT = cms.PSet(
     outputCommands = cms.untracked.vstring(
       'keep *_ticlSimTracksters_*_*',
+      'keep *_ticlSimTICLCandidates_*_*',
       'keep *_ticlSimTrackstersFromCP_*_*',
       )
     )
@@ -48,6 +53,13 @@ def cleanOutputAndSet(outputModule, ticl_outputCommads):
                                             'keep *_layerClusterSimClusterAssociationProducer_*_*',
                                             'keep *_layerClusterCaloParticleAssociationProducer_*_*',
                                             'keep *_randomEngineStateProducer_*_*',
+                                            'keep *_layerClusterSimTracksterAssociationProducer_*_*',
+                                            'keep *_tracksterSimTracksterAssociationLinking_*_*',
+                                            'keep *_tracksterSimTracksterAssociationPR_*_*', 
+                                            'keep *_tracksterSimTracksterAssociationLinkingPU_*_*',
+                                            'keep *_tracksterSimTracksterAssociationPRPU_*_*', 
+                                            'keep *_tracksterSimTracksterAssociationLinkingbyCLUE3D_*_*',
+                                            'keep *_tracksterSimTracksterAssociationPRbyCLUE3D_*_*', 
                                             ])
 
     if hasattr(process, 'FEVTDEBUGEventContent'):

diff --git a/RecoHGCal/TICL/interface/commons.h b/RecoHGCal/TICL/interface/commons.h
@@ -4,12 +4,28 @@
 #include "DataFormats/Common/interface/Ref.h"
 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
 #include "DataFormats/HGCalReco/interface/Trackster.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 
 namespace ticl {
 
   //constants
   constexpr double mpion = 0.13957;
   constexpr float mpion2 = mpion * mpion;
+  typedef math::XYZVectorF Vector;
+  enum LayerType {
+
+    CE_E_120 = 0,
+    CE_E_200 = 1,
+    CE_E_300 = 2,
+    CE_H_120_F = 3,
+    CE_H_200_F = 4,
+    CE_H_300_F = 5,
+    CE_H_200_C = 6,
+    CE_H_300_C = 7,
+    CE_H_SCINT_C = 8,
+    EnumSize = 9
+
+  };
 
   inline Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge) {
     if (pdgId == 111) {
@@ -40,6 +56,44 @@ namespace ticl {
   // verbosity levels for ticl algorithms
   enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru };
 
+  inline int returnClusterType(DetId& lc_seed, const hgcal::RecHitTools& rhtools_) {
+    auto layer_number = rhtools_.getLayerWithOffset(lc_seed);
+    auto thickness = rhtools_.getSiThickIndex(lc_seed);
+    auto isEELayer = (layer_number <= rhtools_.lastLayerEE(false));
+    auto isScintillator = rhtools_.isScintillator(lc_seed);
+    auto isFine = (layer_number <= rhtools_.lastLayerEE(false) + 7);
+
+    if (isScintillator) {
+      return CE_H_SCINT_C;
+    }
+    if (isEELayer) {
+      if (thickness == 0) {
+        return CE_E_120;
+      } else if (thickness == 1) {
+        return CE_E_200;
+      } else if (thickness == 2) {
+        return CE_E_300;
+      }
+    } else {
+      if (isFine) {
+        if (thickness == 0) {
+          return CE_H_120_F;
+        } else if (thickness == 1) {
+          return CE_H_200_F;
+        } else if (thickness == 2) {
+          return CE_H_300_F;
+        }
+      } else {
+        if (thickness == 1) {
+          return CE_H_200_C;
+        } else if (thickness == 2) {
+          return CE_H_300_C;
+        }
+      }
+    }
+    return -1;
+  };
+
 }  // namespace ticl
 
 #endif
diff --git a/RecoHGCal/TICL/plugins/BuildFile.xml b/RecoHGCal/TICL/plugins/BuildFile.xml
@@ -29,6 +29,8 @@
 <use name="TrackingTools/GeomPropagators"/>
 <use name="TrackingTools/Records"/>
 <use name="TrackingTools/TrajectoryState"/>
+<use name="PhysicsTools/UtilAlgos"/>
+<use name="FWCore/ServiceRegistry"/>
 <use name="fastjet"/>
 <library file="*.cc" name="RecoHGCalTICLPlugins">
   <flags EDM_PLUGIN="1"/>

diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc
@@ -10,6 +10,7 @@
 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
 
 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
+#include "DataFormats/HGCalReco/interface/Trackster.h"
 
 using namespace ticl;
 
@@ -36,10 +37,10 @@ void LinkingAlgoByDirectionGeometric::initialize(const HGCalDDDConstants *hgcons
   propagator_ = propH;
 }
 
-math::XYZVector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t,
-                                                                    const unsigned idx,
-                                                                    float zVal,
-                                                                    std::array<TICLLayerTile, 2> &tracksterTiles) {
+Vector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t,
+                                                           const unsigned idx,
+                                                           float zVal,
+                                                           std::array<TICLLayerTile, 2> &tracksterTiles) {
   // needs only the positive Z co-ordinate of the surface to propagate to
   // the correct sign is calculated inside according to the barycenter of trackster
   Vector const &baryc = t.barycenter();

diff --git a/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h b/RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h
@@ -20,6 +20,7 @@
 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
 
 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+#include "DataFormats/HGCalReco/interface/Trackster.h"
 
 namespace ticl {
   class LinkingAlgoByDirectionGeometric final : public LinkingAlgoBase {
@@ -44,14 +45,14 @@ namespace ticl {
     static void fillPSetDescription(edm::ParameterSetDescription &desc);
 
   private:
-    typedef math::XYZVector Vector;
+    using Vector = ticl::Trackster::Vector;
 
     void buildLayers();
 
-    math::XYZVector propagateTrackster(const Trackster &t,
-                                       const unsigned idx,
-                                       float zVal,
-                                       std::array<TICLLayerTile, 2> &tracksterTiles);
+    Vector propagateTrackster(const Trackster &t,
+                              const unsigned idx,
+                              float zVal,
+                              std::array<TICLLayerTile, 2> &tracksterTiles);
 
     void findTrackstersInWindow(const std::vector<std::pair<Vector, unsigned>> &seedingCollection,
                                 const std::array<TICLLayerTile, 2> &tracksterTiles,