Merge pull request #2 from felicepantaleo/simTICLCandidates_13_0_0_pre2

SimTICLCandidate

DataFormats/HGCalReco/interface/TICLCandidate.h

@@ -42,18 +42,21 @@ 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];
   }
 
+  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_;

RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py

@@ -28,6 +28,7 @@
 TICL_FEVT = cms.PSet(
     outputCommands = cms.untracked.vstring(
       'keep *_ticlSimTracksters_*_*',
+      'keep *_ticlSimTICLCandidates_*_*',
       'keep *_ticlSimTrackstersFromCP_*_*',
       )
     )

RecoHGCal/TICL/plugins/SimTICLCandidatesProducer.cc

@@ -0,0 +1,195 @@
+// Author: Felice Pantaleo, Wahid Redjeb - felice.pantaleo@cern.ch, wahid.redjeb@cern.ch
+// Date: 01/2023
+
+// user include files
+
+#include "FWCore/Framework/interface/ESHandle.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/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "DataFormats/HGCalReco/interface/Trackster.h"
+#include "DataFormats/HGCalReco/interface/TICLCandidate.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
+
+#include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "RecoHGCal/TICL/interface/commons.h"
+#include "SimDataFormats/Track/interface/UniqueSimTrackId.h"
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
+
+// #include "TrackstersPCA.h"
+#include <vector>
+#include <map>
+#include <iterator>
+#include <algorithm>
+
+using namespace ticl;
+
+class SimTICLCandidatesProducer : public edm::stream::EDProducer<> {
+public:
+  explicit SimTICLCandidatesProducer(const edm::ParameterSet&);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+private:
+  std::string detector_;
+  const bool doNose_ = false;
+
+  const edm::EDGetTokenT<std::vector<SimCluster>> simclusters_token_;
+  const edm::EDGetTokenT<std::vector<CaloParticle>> caloparticles_token_;
+  const edm::EDGetTokenT<std::vector<Trackster>> simTrackstersToken_;
+  const edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticleToken_;
+  const edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken_;
+  const StringCutObjectSelector<reco::Track> cutTk_;
+
+  const edm::EDGetTokenT<reco::SimToRecoCollection> associatormapStRsToken_;
+  const edm::EDGetTokenT<reco::RecoToSimCollection> associatormapRtSsToken_;
+  const edm::EDGetTokenT<SimTrackToTPMap> associationSimTrackToTPToken_;
+};
+DEFINE_FWK_MODULE(SimTICLCandidatesProducer);
+
+SimTICLCandidatesProducer::SimTICLCandidatesProducer(const edm::ParameterSet& ps)
+    : detector_(ps.getParameter<std::string>("detector")),
+      doNose_(detector_ == "HFNose"),
+      simclusters_token_(consumes(ps.getParameter<edm::InputTag>("simclusters"))),
+      caloparticles_token_(consumes(ps.getParameter<edm::InputTag>("caloparticles"))),
+      simTrackstersToken_(consumes(ps.getParameter<edm::InputTag>("simTracksters"))),
+      trackingParticleToken_(
+          consumes<std::vector<TrackingParticle>>(ps.getParameter<edm::InputTag>("trackingParticles"))),
+      recoTracksToken_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("recoTracks"))),
+      cutTk_(ps.getParameter<std::string>("cutTk")),
+      associatormapStRsToken_(consumes(ps.getParameter<edm::InputTag>("tpToTrack"))),
+      associatormapRtSsToken_(consumes(ps.getParameter<edm::InputTag>("trackToTp"))),
+      associationSimTrackToTPToken_(consumes(ps.getParameter<edm::InputTag>("simTrackToTPMap"))) {
+  produces<std::vector<TICLCandidate>>();
+}
+
+void SimTICLCandidatesProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("detector", "HGCAL");
+  desc.add<edm::InputTag>("recoTracks", edm::InputTag("generalTracks"));
+  desc.add<std::string>("cutTk",
+                        "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
+                        "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
+  desc.add<edm::InputTag>("tpToTrack", edm::InputTag("trackingParticleRecoTrackAsssociation"));
+  desc.add<edm::InputTag>("trackToTp", edm::InputTag("trackingParticleRecoTrackAsssociation"));
+  desc.add<edm::InputTag>("simclusters", edm::InputTag("mix", "MergedCaloTruth"));
+  desc.add<edm::InputTag>("caloparticles", edm::InputTag("mix", "MergedCaloTruth"));
+  desc.add<edm::InputTag>("simTracksters", edm::InputTag("ticlSimTracksters"));
+  desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
+  desc.add<edm::InputTag>("simTrackToTPMap", edm::InputTag("simHitTPAssocProducer", "simTrackToTP"));
+  descriptions.addWithDefaultLabel(desc);
+}
+
+void SimTICLCandidatesProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
+  auto result = std::make_unique<std::vector<TICLCandidate>>();
+
+  edm::Handle<std::vector<CaloParticle>> caloParticles_h;
+  evt.getByToken(caloparticles_token_, caloParticles_h);
+  edm::Handle<std::vector<TrackingParticle>> trackingParticles_h;
+  evt.getByToken(trackingParticleToken_, trackingParticles_h);
+  edm::Handle<std::vector<reco::Track>> recoTracks_h;
+  evt.getByToken(recoTracksToken_, recoTracks_h);
+  edm::Handle<std::vector<Trackster>> simTracksters_h;
+  evt.getByToken(simTrackstersToken_, simTracksters_h);
+  const auto& simClusters = evt.get(simclusters_token_);
+  const auto& caloParticles = *caloParticles_h;
+  const auto& simTracksters = *simTracksters_h;
+  const auto& trackingParticles = *trackingParticles_h;
+  const auto& TPtoRecoTrackMap = evt.get(associatormapStRsToken_);
+  const auto& simTrackToTPMap = evt.get(associationSimTrackToTPToken_);
+  const auto& recoTracks = *recoTracks_h;
+  result->reserve(simTracksters.size());
+
+  std::vector<int> usedSimTrackster(simTracksters.size(), 0);
+
+  // Creating the map from TrackingParticle to SimTrackster
+  std::unordered_map<unsigned int, std::vector<unsigned int>> TPtoSimTracksterMap;
+  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 ipos = simTrackToTPMap.mapping.find(simTkIds);
+    if (ipos != simTrackToTPMap.mapping.end()) {
+      auto tpIdx = (ipos->second).get() - (edm::Ref<std::vector<TrackingParticle>>(trackingParticles_h, 0)).get();
+      TPtoSimTracksterMap[tpIdx].push_back(i);
+    }
+  }
+
+  for (const auto& [tpKey, associatedRecoTracks] : TPtoRecoTrackMap) {
+    auto const& tp = *(tpKey);
+    auto tpIndex = &tp - &trackingParticles[0];
+
+    if (!associatedRecoTracks.empty()) {
+      if (associatedRecoTracks[0].second > 0.75f) {
+        auto const& track = associatedRecoTracks[0].first;
+        auto trackIndex = &(*track) - &recoTracks[0];
+
+        if (cutTk_(*track)) {
+          if (!TPtoSimTracksterMap[tpIndex].empty()) {
+            TICLCandidate tmpCand;
+            tmpCand.zeroProbabilities();
+            auto const particleType = tracksterParticleTypeFromPdgId(tp.pdgId(), 1);
+            tmpCand.setIdProbability(particleType, 1.f);
+            tmpCand.setTrackPtr(edm::Ptr<reco::Track>(recoTracks_h, trackIndex));
+            tmpCand.setPdgId(tp.pdgId());
+            tmpCand.setCharge(tp.charge());
+            float rawEnergy = 0.f;
+            float regressedEnergy = 0.f;
+            for (auto simTracksterIndex : TPtoSimTracksterMap[tpIndex]) {
+              if (usedSimTrackster[simTracksterIndex] == 0) {
+                usedSimTrackster[simTracksterIndex] = 1;
+                rawEnergy += simTracksters[simTracksterIndex].raw_energy();
+                regressedEnergy += simTracksters[simTracksterIndex].regressed_energy();
+                tmpCand.addTrackster(edm::Ptr<Trackster>(simTracksters_h, simTracksterIndex));
+              }
+            }
+            math::XYZTLorentzVector p4(regressedEnergy * track->momentum().unit().x(),
+                                       regressedEnergy * track->momentum().unit().y(),
+                                       regressedEnergy * track->momentum().unit().z(),
+                                       regressedEnergy);
+            tmpCand.setP4(p4);
+            tmpCand.setRawEnergy(rawEnergy);
+            result->push_back(tmpCand);
+          }
+        }
+      }
+    }
+  }
+
+  for (unsigned int i = 0; i < simTracksters.size(); ++i) {
+    if (usedSimTrackster[i] == 0) {
+      usedSimTrackster[i] = 1;
+      TICLCandidate tmpCand;
+      const auto& simTrackster = simTracksters[i];
+      tmpCand.setIdProbabilities(simTrackster.id_probabilities());
+      const auto pdgId = (simTracksters[i].seedID() == caloParticles_h.id())
+                             ? caloParticles[simTrackster.seedIndex()].pdgId()
+                             : simClusters[simTrackster.seedIndex()].pdgId();
+      tmpCand.setPdgId(pdgId);
+      tmpCand.setCharge(0);
+      tmpCand.setRawEnergy(simTrackster.raw_energy());
+      float regressedEnergy = simTrackster.regressed_energy();
+      math::XYZTLorentzVector p4(regressedEnergy * simTrackster.barycenter().unit().x(),
+                                 regressedEnergy * simTrackster.barycenter().unit().y(),
+                                 regressedEnergy * simTrackster.barycenter().unit().z(),
+                                 regressedEnergy);
+      tmpCand.setP4(p4);
+      result->push_back(tmpCand);
+    }
+  }
+
+  result->shrink_to_fit();
+
+  evt.put(std::move(result));
+}

RecoHGCal/TICL/plugins/SimTrackstersProducer.cc

@@ -246,7 +246,8 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
       (*cpToSc_SimTrackstersMap)[index] = scSimTracksterIdx;
     }
   }
-
+  // TODO: remove time computation from PCA calculation and
+  //       store time from boundary position in simTracksters
   ticl::assignPCAtoTracksters(
       *result, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z());
   result->shrink_to_fit();

RecoHGCal/TICL/python/SimTICLCandidates_cff.py

@@ -0,0 +1,14 @@
+import FWCore.ParameterSet.Config as cms
+
+from RecoHGCal.TICL.simTICLCandidatesProducer_cfi import simTICLCandidatesProducer as _simTICLCandidatesProducer
+
+ticlSimTICLCandidates = _simTICLCandidatesProducer.clone(
+)
+
+from Configuration.ProcessModifiers.premix_stage2_cff import premix_stage2
+premix_stage2.toModify(ticlSimTICLCandidates,
+    simclusters = "mixData:MergedCaloTruth",
+    caloparticles = "mixData:MergedCaloTruth",
+)
+
+ticlSimTICLCandidatesTask = cms.Task(ticlSimTICLCandidates)

SimDataFormats/Track/BuildFile.xml

@@ -1,3 +1,4 @@
+<use name="hepmc"/>
 <use name="DataFormats/Common"/>
 <use name="DataFormats/Math"/>
 <use name="SimDataFormats/EncodedEventId"/>

SimDataFormats/Track/interface/UniqueSimTrackId.h

@@ -2,7 +2,8 @@
 #define SimDataFormatsTrackUniqueSimTrackId_H
 
 #include "FWCore/Utilities/interface/hash_combine.h"
-
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
 #include <tuple>
 
@@ -13,4 +14,7 @@ struct UniqueSimTrackIdHash {
   }
 };
 
+struct SimTrackToTPMap {
+  std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash> mapping;
+};
 #endif

SimDataFormats/Track/src/classes.h

@@ -1,6 +1,7 @@
 #include "SimDataFormats/Track/interface/CoreSimTrack.h"
 #include "SimDataFormats/Track/interface/SimTrack.h"
 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
+#include "SimDataFormats/Track/interface/UniqueSimTrackId.h"
 #include "DataFormats/Common/interface/Wrapper.h"
 
 #include <vector>

SimDataFormats/Track/src/classes_def.xml

@@ -13,4 +13,8 @@
   <class name="edm::RefProd<std::vector<SimTrack> >"/>
   <class name="edm::Ref<std::vector<SimTrack>,SimTrack,edm::refhelper::FindUsingAdvance<std::vector<SimTrack>,SimTrack> >"/>
   <class name="edm::RefVector<std::vector<SimTrack>,SimTrack,edm::refhelper::FindUsingAdvance<std::vector<SimTrack>,SimTrack> >"/>
+  <class name="std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash>"/>
+  <class name="edm::Wrapper<std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash>>"/>
+  <class name="SimTrackToTPMap"/>
+  <class name="edm::Wrapper<SimTrackToTPMap>"/>
 </lcgdict>

SimGeneral/TrackingAnalysis/plugins/SimHitTPAssociationProducer.cc

@@ -27,6 +27,7 @@ SimHitTPAssociationProducer::SimHitTPAssociationProducer(const edm::ParameterSet
       _trackingParticleSrc(
           consumes<TrackingParticleCollection>(cfg.getParameter<edm::InputTag>("trackingParticleSrc"))) {
   produces<SimHitTPAssociationList>();
+  produces<SimTrackToTPMap>("simTrackToTP");
   std::vector<edm::InputTag> tags = cfg.getParameter<std::vector<edm::InputTag>>("simHitSrc");
   _simHitSrc.reserve(tags.size());
   for (auto const &tag : tags) {
@@ -44,7 +45,7 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con
   iEvent.getByToken(_trackingParticleSrc, TPCollectionH);
 
   // prepare temporary map between SimTrackId and TrackingParticle index
-  std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash> mapping;
+  auto simTrackToTPMap = std::make_unique<SimTrackToTPMap>();
   auto const &tpColl = *TPCollectionH.product();
   for (TrackingParticleCollection::size_type itp = 0, size = tpColl.size(); itp < size; ++itp) {
     auto const &trackingParticle = tpColl[itp];
@@ -53,7 +54,7 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con
     EncodedEventId eid(trackingParticle.eventId());
     for (auto const &trk : trackingParticle.g4Tracks()) {
       UniqueSimTrackId trkid(trk.trackId(), eid);
-      mapping.insert(std::make_pair(trkid, trackingParticleRef));
+      simTrackToTPMap->mapping.insert(std::make_pair(trkid, trackingParticleRef));
     }
   }
 
@@ -66,15 +67,16 @@ void SimHitTPAssociationProducer::produce(edm::StreamID, edm::Event &iEvent, con
       TrackPSimHitRef pSimHitRef(PSimHitCollectionH, psimHitI);
       auto const &pSimHit = pSimHitCollection[psimHitI];
       UniqueSimTrackId simTkIds(pSimHit.trackId(), pSimHit.eventId());
-      auto ipos = mapping.find(simTkIds);
-      if (ipos != mapping.end()) {
+      auto ipos = simTrackToTPMap->mapping.find(simTkIds);
+      if (ipos != simTrackToTPMap->mapping.end()) {
         simHitTPList->emplace_back(ipos->second, pSimHitRef);
       }
     }
   }
 
   std::sort(simHitTPList->begin(), simHitTPList->end(), simHitTPAssociationListGreater);
   iEvent.put(std::move(simHitTPList));
+  iEvent.put(std::move(simTrackToTPMap), "simTrackToTP");
 }
 
 #include "FWCore/Framework/interface/MakerMacros.h"

SimTracker/Configuration/python/SimTracker_EventContent_cff.py

@@ -15,7 +15,8 @@
         'keep *_assocOutInConversionTracks_*_*', 
         'keep *_assocInOutConversionTracks_*_*',
         'keep *_TTClusterAssociatorFromPixelDigis_*_*',
-        'keep *_TTStubAssociatorFromPixelDigis_*_*')
+        'keep *_TTStubAssociatorFromPixelDigis_*_*',
+        'keep *_simHitTPAssocProducer_*_*')
 
 )
 # For phase2 premixing switch the sim digi collections to the ones including pileup

Validation/Configuration/python/globalValidation_cff.py

@@ -200,7 +200,7 @@
 )
 
 globalValidationHGCal = cms.Sequence(hgcalValidation)
-globalPrevalidationHGCal = cms.Sequence(hgcalAssociators, ticlSimTrackstersTask)
+globalPrevalidationHGCal = cms.Sequence(hgcalAssociators, ticlSimTrackstersTask, ticlSimTICLCandidatesTask)
 
 globalValidationMTD = cms.Sequence()
 

Validation/Configuration/python/hgcalSimValid_cff.py

@@ -11,6 +11,8 @@
 from Validation.HGCalValidation.rechitValidation_cff    import *
 from Validation.HGCalValidation.hgcalHitValidation_cff  import *
 from RecoHGCal.TICL.SimTracksters_cff import *
+from RecoHGCal.TICL.SimTICLCandidates_cff import *
+
 
 from Validation.HGCalValidation.HGCalValidator_cfi import hgcalValidator
 from Validation.RecoParticleFlow.PFJetValidation_cff import pfJetValidation1 as _hgcalPFJetValidation