Merge pull request #46678 from AuroraPerego/timeHGCalCP

Add time to HGCAL CaloParticles

RecoHGCal/TICL/plugins/SimTrackstersProducer.cc

@@ -32,15 +32,15 @@
 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTracksterFwd.h"
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 
-#include "SimDataFormats/Vertex/interface/SimVertex.h"
-
 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
 #include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h"
 
 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
 
 #include "DataFormats/HGCalReco/interface/Common.h"
 
+#include <CLHEP/Units/SystemOfUnits.h>
+
 #include "TrackstersPCA.h"
 #include <vector>
 #include <map>
@@ -96,7 +96,6 @@ class SimTrackstersProducer : public edm::stream::EDProducer<> {
   const float fractionCut_;
   const float qualityCutTrack_;
   const edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticleToken_;
-  const edm::EDGetTokenT<std::vector<SimVertex>> simVerticesToken_;
 
   const edm::EDGetTokenT<std::vector<reco::Track>> recoTracksToken_;
   const StringCutObjectSelector<reco::Track> cutTk_;
@@ -127,7 +126,6 @@ SimTrackstersProducer::SimTrackstersProducer(const edm::ParameterSet& ps)
       qualityCutTrack_(ps.getParameter<double>("qualityCutTrack")),
       trackingParticleToken_(
           consumes<std::vector<TrackingParticle>>(ps.getParameter<edm::InputTag>("trackingParticles"))),
-      simVerticesToken_(consumes<std::vector<SimVertex>>(ps.getParameter<edm::InputTag>("simVertices"))),
       recoTracksToken_(consumes<std::vector<reco::Track>>(ps.getParameter<edm::InputTag>("recoTracks"))),
       cutTk_(ps.getParameter<std::string>("cutTk")),
       associatormapStRsToken_(consumes(ps.getParameter<edm::InputTag>("tpToTrack"))),
@@ -162,7 +160,6 @@ void SimTrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& des
   desc.add<edm::InputTag>("tpToTrack", edm::InputTag("trackingParticleRecoTrackAsssociation"));
 
   desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
-  desc.add<edm::InputTag>("simVertices", edm::InputTag("g4SimHits"));
 
   desc.add<edm::InputTag>("simTrackToTPMap", edm::InputTag("simHitTPAssocProducer", "simTrackToTP"));
   desc.add<double>("fractionCut", 0.);
@@ -225,7 +222,7 @@ void SimTrackstersProducer::addTrackster(
   tmpTrackster.setRegressedEnergy(energy);
   tmpTrackster.setIteration(iter);
   tmpTrackster.setSeed(seed, index);
-  tmpTrackster.setBoundaryTime(time * 1e9);
+  tmpTrackster.setBoundaryTime(time);
   if (add) {
     result[index] = tmpTrackster;
     loop_index += 1;
@@ -259,7 +256,6 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
 
   const auto& simClustersToRecoColl = evt.get(associatorMapSimClusterToReco_token_);
   const auto& caloParticlesToRecoColl = evt.get(associatorMapCaloParticleToReco_token_);
-  const auto& simVertices = evt.get(simVerticesToken_);
 
   edm::Handle<std::vector<TrackingParticle>> trackingParticles_h;
   evt.getByToken(trackingParticleToken_, trackingParticles_h);
@@ -293,7 +289,8 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
     // Create a Trackster from the object entering HGCal
     if (cp.g4Tracks()[0].crossedBoundary()) {
       regr_energy = cp.g4Tracks()[0].getMomentumAtBoundary().energy();
-      float time = cp.g4Tracks()[0].getPositionAtBoundary().t();
+      float time = cp.g4Tracks()[0].getPositionAtBoundary().t() *
+                   CLHEP::s;  // Geant4 time is in seconds, convert to ns (CLHEP::s = 1e9)
       addTrackster(cpIndex,
                    lcVec,
                    inputClusterMask,
@@ -323,7 +320,8 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
                      sc.g4Tracks()[0].getMomentumAtBoundary().energy(),
                      sc.pdgId(),
                      sc.charge(),
-                     sc.g4Tracks()[0].getPositionAtBoundary().t(),
+                     sc.g4Tracks()[0].getPositionAtBoundary().t() *
+                         CLHEP::s,  // Geant4 time is in seconds, convert to ns (CLHEP::s = 1e9)
                      scRef.id(),
                      ticl::Trackster::SIM,
                      *output_mask,
@@ -339,7 +337,7 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
       }
       scSimTracksterIdx.shrink_to_fit();
     }
-    float time = simVertices[cp.g4Tracks()[0].vertIndex()].position().t();
+    float time = cp.simTime();
     // Create a Trackster from any CP
     addTrackster(cpIndex,
                  lcVec,
@@ -482,7 +480,7 @@ void SimTrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es)
       cand.setMTDTime(MTDst->time(), 0);
     }
 
-    cand.setTime(simVertices[cp.g4Tracks()[0].vertIndex()].position().t() * pow(10, 9), 0);
+    cand.setTime(cp.simTime(), 0);
 
     for (const auto& trackster : cand.tracksters()) {
       rawEnergy += trackster->raw_energy();

SimDataFormats/CaloAnalysis/interface/CaloParticle.h

@@ -165,6 +165,9 @@ class CaloParticle {
   /** @brief Gives the total number of SimHits, in the cluster */
   int numberOfRecHits() const { return hits_.size(); }
 
+  /** @brief returns the time in ns of the caloparticle */
+  float simTime() const { return time_; }
+
   /** @brief add rechit with fraction */
   void addRecHitAndFraction(uint32_t hit, float fraction) {
     hits_.emplace_back(hit);
@@ -189,12 +192,16 @@ class CaloParticle {
     ++nsimhits_;
   }
 
+  /** @brief add vertex time to the caloparticle */
+  void setSimTime(const float time) { time_ = time; }
+
 protected:
   uint64_t nsimhits_{0};
   EncodedEventId event_;
 
   uint32_t particleId_{0};
   float simhit_energy_{0.f};
+  float time_{std::numeric_limits<float>::lowest()};
   std::vector<uint32_t> hits_;
   std::vector<float> fractions_;
 

SimDataFormats/CaloAnalysis/interface/MtdCaloParticle.h

@@ -33,19 +33,13 @@ class MtdCaloParticle : public CaloParticle {
   const MtdSimClusterRefVector &simClusters() const { return mtdsimClusters_; }
   void clearSimClusters() { mtdsimClusters_.clear(); }
 
-  /** @brief returns the time of the caloparticle */
-  float simTime() const { return simhit_time_; }
-
-  void addSimTime(const float time) { simhit_time_ = time; }
-
   /** @brief add simhit's energy to cluster */
   void addSimHit(PSimHit &hit) {
     simhit_energy_ += hit.energyLoss();
     ++nsimhits_;
   }
 
 private:
-  float simhit_time_{-99.f};
   MtdSimClusterRefVector mtdsimClusters_;
 };
 

SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc

@@ -47,6 +47,8 @@
 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
 #include "Geometry/Records/interface/CaloGeometryRecord.h"
 
+#include <CLHEP/Units/SystemOfUnits.h>
+
 namespace {
   using Index_t = unsigned;
   using Barcode_t = int;
@@ -151,11 +153,13 @@ namespace {
                              CaloTruthAccumulator::calo_particles &caloParticles,
                              std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex,
                              std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
+                             std::unordered_map<uint32_t, float> &vertex_time_map,
                              Selector selector)
         : output_(output),
           caloParticles_(caloParticles),
           simHitBarcodeToIndex_(simHitBarcodeToIndex),
           simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap),
+          vertex_time_map_(vertex_time_map),
           selector_(selector) {}
     template <typename Vertex, typename Graph>
     void discover_vertex(Vertex u, const Graph &g) {
@@ -191,6 +195,7 @@ namespace {
         if (selector_(edge_property)) {
           IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
           output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
+          output_.pCaloParticles->back().setSimTime(vertex_time_map_[(edge_property.simTrack)->vertIndex()]);
           caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
         }
       }
@@ -213,6 +218,7 @@ namespace {
     CaloTruthAccumulator::calo_particles &caloParticles_;
     std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex_;
     std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap_;
+    std::unordered_map<uint32_t, float> &vertex_time_map_;
     Selector selector_;
   };
 }  // namespace
@@ -425,6 +431,12 @@ void CaloTruthAccumulator::accumulateEvent(const T &event,
     idx++;
   }
 
+  std::unordered_map<uint32_t, float> vertex_time_map;
+  for (uint32_t i = 0; i < vertices.size(); i++) {
+    // Geant4 time is in seconds, convert to ns (CLHEP::s = 1e9)
+    vertex_time_map[i] = vertices[i].position().t() * CLHEP::s;
+  }
+
   /**
   Build the main decay graph and assign the SimTrack to each edge. The graph
   built here will only contain the particles that have a decay vertex
@@ -512,6 +524,7 @@ void CaloTruthAccumulator::accumulateEvent(const T &event,
       m_caloParticles,
       m_simHitBarcodeToIndex,
       simTrackDetIdEnergyMap,
+      vertex_time_map,
       [&](EdgeProperty &edge_property) -> bool {
         // Apply selection on SimTracks in order to promote them to be
         // CaloParticles. The function returns TRUE if the particle satisfies

SimGeneral/CaloAnalysis/plugins/MtdTruthAccumulator.cc

@@ -58,6 +58,8 @@
 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
 
+#include <CLHEP/Units/SystemOfUnits.h>
+
 namespace {
   using Index_t = unsigned;
   using Barcode_t = int;
@@ -106,8 +108,6 @@ class MtdTruthAccumulator : public DigiAccumulatorMixMod {
   */
   const unsigned int maximumSubsequentBunchCrossing_;
 
-  const unsigned int bunchSpacing_;
-
   const edm::InputTag simTrackLabel_;
   const edm::InputTag simVertexLabel_;
   edm::Handle<std::vector<SimTrack>> hSimTracks;
@@ -223,7 +223,7 @@ namespace {
         if (selector_(edge_property)) {
           IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
           output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
-          output_.pCaloParticles->back().addSimTime(vertex_time_map_[(edge_property.simTrack)->vertIndex()]);
+          output_.pCaloParticles->back().setSimTime(vertex_time_map_[(edge_property.simTrack)->vertIndex()]);
           caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
         }
       }
@@ -257,7 +257,6 @@ MtdTruthAccumulator::MtdTruthAccumulator(const edm::ParameterSet &config,
     : messageCategory_("MtdTruthAccumulator"),
       maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
       maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
-      bunchSpacing_(config.getParameter<unsigned int>("bunchspace")),
       simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
       simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
       collectionTags_(),
@@ -619,7 +618,8 @@ void MtdTruthAccumulator::accumulateEvent(const T &event,
   DecayChain decay;
 
   for (uint32_t i = 0; i < vertices.size(); i++) {
-    vertex_time_map[i] = vertices[i].position().t() * 1e9 + event.bunchCrossing() * static_cast<int>(bunchSpacing_);
+    // Geant4 time is in seconds, convert to ns (CLHEP::s = 1e9)
+    vertex_time_map[i] = vertices[i].position().t() * CLHEP::s;
   }
 
   IfLogDebug(DEBUG, messageCategory_) << " TRACKS" << std::endl;

SimGeneral/MixingModule/python/mtdTruthProducer_cfi.py

@@ -7,7 +7,6 @@
     premixStage1 = cms.bool(False),
     maximumPreviousBunchCrossing = cms.uint32(0),
     maximumSubsequentBunchCrossing = cms.uint32(0),
-    bunchspace = cms.uint32(25), #ns
 
     simHitCollections = cms.PSet(
         mtdCollections = cms.VInputTag(