Merge pull request #33318 from felicepantaleo/TICL_TrackstersFromSimC…

…lustersProducer_11_3_X_2021-03-30-2300

[HGCAL] TICL SimTracksters

DataFormats/HGCalReco/interface/Trackster.h

@@ -114,7 +114,7 @@ namespace ticl {
         p = *(probs++);
       }
     }
-    inline void setIdProbability(ParticleType type, float value) { id_probabilities_[int(type)] = 1.f; }
+    inline void setIdProbability(ParticleType type, float value) { id_probabilities_[int(type)] = value; }
 
     inline const Trackster::IterationIndex ticlIteration() const { return (IterationIndex)iterationIndex_; }
     inline const std::vector<unsigned int> &vertices() const { return vertices_; }

RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py

@@ -3,13 +3,15 @@
 
 #AOD content
 TICL_AOD = cms.PSet(
+  # 13/04/2021 Felice: MultiClusters will be deprecated soon
     outputCommands = cms.untracked.vstring(
       'keep *_ticlMultiClustersFromTrackstersEM_*_*',
       'keep *_ticlMultiClustersFromTrackstersHAD_*_*',
       'keep *_ticlMultiClustersFromTrackstersTrk_*_*',
       'keep *_ticlMultiClustersFromTrackstersTrkEM_*_*',
       'keep *_ticlMultiClustersFromTrackstersMIP_*_*',
       'keep *_ticlMultiClustersFromTrackstersMerge_*_*',
+      'keep *_ticlMultiClustersFromSimTracksters_*_*',
       )
     )
 
@@ -35,6 +37,7 @@
 # FEVT Content
 TICL_FEVT = cms.PSet(
     outputCommands = cms.untracked.vstring(
+      'keep *_ticlSimTracksters_*_*',
       )
     )
 TICL_FEVT.outputCommands.extend(TICL_RECO.outputCommands)

RecoHGCal/TICL/plugins/BuildFile.xml

@@ -1,18 +1,36 @@
 <use name="CommonTools/Utils"/>
-<use name="RecoLocalCalo/HGCalRecProducers"/>
-<use name="RecoLocalCalo/HGCalRecAlgos"/>
-<use name="FWCore/Framework"/>
-<use name="FWCore/ParameterSet"/>
-<use name="FWCore/PluginManager"/>
 <use name="DataFormats/CaloRecHit"/>
-<use name="DataFormats/ParticleFlowCandidate"/>
+<use name="DataFormats/Candidate"/>
+<use name="DataFormats/Common"/>
+<use name="DataFormats/GeometrySurface"/>
 <use name="DataFormats/HGCalReco"/>
-<use name="FWCore/MessageLogger"/>
+<use name="DataFormats/L1TCorrelator"/>
+<use name="DataFormats/L1Trigger"/>
+<use name="DataFormats/Math"/>
+<use name="DataFormats/ParticleFlowCandidate"/>
+<use name="DataFormats/ParticleFlowReco"/>
 <use name="DataFormats/TrackReco"/>
-<use name="TrackingTools/Records"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/MessageLogger"/>
+<use name="FWCore/ParameterSet"/>
+<use name="FWCore/PluginManager"/>
+<use name="FWCore/Utilities"/>
+<use name="Geometry/CaloGeometry"/>
+<use name="Geometry/CommonDetUnit"/>
+<use name="Geometry/HGCalCommonData"/>
+<use name="Geometry/Records"/>
+<use name="MagneticField/Engine"/>
+<use name="MagneticField/Records"/>
 <use name="PhysicsTools/TensorFlow"/>
+<use name="RecoLocalCalo/HGCalRecAlgos"/>
+<use name="RecoLocalCalo/HGCalRecProducers"/>
 <use name="RecoHGCal/TICL"/>
 <use name="RecoParticleFlow/PFProducer"/>
+<use name="SimDataFormats/Associations"/>
+<use name="SimDataFormats/CaloAnalysis"/>
+<use name="TrackingTools/GeomPropagators"/>
+<use name="TrackingTools/Records"/>
+<use name="TrackingTools/TrajectoryState"/>
 <library file="*.cc" name="RecoHGCalTICLPlugins">
   <flags EDM_PLUGIN="1"/>
 </library>

RecoHGCal/TICL/plugins/TrackstersFromSimClustersProducer.cc

@@ -0,0 +1,151 @@
+// Author: Felice Pantaleo - [email protected]
+// Date: 02/2021
+
+// 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/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
+
+#include "DataFormats/HGCalReco/interface/Trackster.h"
+
+#include "DataFormats/Common/interface/ValueMap.h"
+#include "SimDataFormats/Associations/interface/LayerClusterToSimClusterAssociator.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+
+#include "TrackstersPCA.h"
+#include <vector>
+#include <iterator>
+#include <algorithm>
+using namespace ticl;
+
+namespace {
+  Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge) {
+    if (pdgId == 111) {
+      return Trackster::ParticleType::neutral_pion;
+    } else {
+      pdgId = std::abs(pdgId);
+      if (pdgId == 22) {
+        return Trackster::ParticleType::photon;
+      } else if (pdgId == 11) {
+        return Trackster::ParticleType::electron;
+      } else if (pdgId == 13) {
+        return Trackster::ParticleType::muon;
+      } else {
+        bool isHadron = (pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000);
+        if (isHadron) {
+          if (charge != 0) {
+            return Trackster::ParticleType::charged_hadron;
+          } else {
+            return Trackster::ParticleType::neutral_hadron;
+          }
+        } else {
+          return Trackster::ParticleType::unknown;
+        }
+      }
+    }
+  }
+}  // namespace
+
+class TrackstersFromSimClustersProducer : public edm::stream::EDProducer<> {
+public:
+  explicit TrackstersFromSimClustersProducer(const edm::ParameterSet&);
+  ~TrackstersFromSimClustersProducer() override {}
+  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<reco::CaloCluster>> clusters_token_;
+  const edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
+  const edm::EDGetTokenT<std::vector<float>> filtered_layerclusters_mask_token_;
+
+  edm::EDGetTokenT<std::vector<SimCluster>> simclusters_token_;
+
+  edm::InputTag associatorLayerClusterSimCluster_;
+  edm::EDGetTokenT<hgcal::SimToRecoCollectionWithSimClusters> associatorMapSimToReco_token_;
+  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geom_token_;
+  hgcal::RecHitTools rhtools_;
+};
+DEFINE_FWK_MODULE(TrackstersFromSimClustersProducer);
+
+TrackstersFromSimClustersProducer::TrackstersFromSimClustersProducer(const edm::ParameterSet& ps)
+    : detector_(ps.getParameter<std::string>("detector")),
+      doNose_(detector_ == "HFNose"),
+      clusters_token_(consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"))),
+      clustersTime_token_(
+          consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("time_layerclusters"))),
+      filtered_layerclusters_mask_token_(consumes<std::vector<float>>(ps.getParameter<edm::InputTag>("filtered_mask"))),
+      simclusters_token_(consumes<std::vector<SimCluster>>(ps.getParameter<edm::InputTag>("simclusters"))),
+      associatorLayerClusterSimCluster_(ps.getUntrackedParameter<edm::InputTag>("layerClusterSimClusterAssociator")),
+      associatorMapSimToReco_token_(
+          consumes<hgcal::SimToRecoCollectionWithSimClusters>(associatorLayerClusterSimCluster_)),
+      geom_token_(esConsumes()) {
+  produces<std::vector<Trackster>>();
+  produces<std::vector<float>>();
+}
+
+void TrackstersFromSimClustersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  // hgcalMultiClusters
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("detector", "HGCAL");
+  desc.add<edm::InputTag>("layer_clusters", edm::InputTag("hgcalLayerClusters"));
+  desc.add<edm::InputTag>("time_layerclusters", edm::InputTag("hgcalLayerClusters", "timeLayerCluster"));
+  desc.add<edm::InputTag>("filtered_mask", edm::InputTag("filteredLayerClustersSimTracksters", "ticlSimTracksters"));
+  desc.add<edm::InputTag>("simclusters", edm::InputTag("mix", "MergedCaloTruth"));
+  desc.addUntracked<edm::InputTag>("layerClusterSimClusterAssociator",
+                                   edm::InputTag("layerClusterSimClusterAssociationProducer"));
+  descriptions.add("trackstersFromSimClustersProducer", desc);
+}
+
+void TrackstersFromSimClustersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
+  auto result = std::make_unique<std::vector<Trackster>>();
+  auto output_mask = std::make_unique<std::vector<float>>();
+  const auto& layerClusters = evt.get(clusters_token_);
+  const auto& layerClustersTimes = evt.get(clustersTime_token_);
+  const auto& inputClusterMask = evt.get(filtered_layerclusters_mask_token_);
+  output_mask->resize(layerClusters.size(), 1.f);
+
+  const auto& simclusters = evt.get(simclusters_token_);
+  const auto& simToRecoColl = evt.get(associatorMapSimToReco_token_);
+
+  const auto& geom = es.getData(geom_token_);
+  rhtools_.setGeometry(geom);
+  auto num_simclusters = simclusters.size();
+  result->reserve(num_simclusters);
+  for (const auto& [key, values] : simToRecoColl) {
+    auto const& sc = *(key);
+    auto simClusterIndex = &sc - &simclusters[0];
+    Trackster tmpTrackster;
+    tmpTrackster.zeroProbabilities();
+    tmpTrackster.vertices().reserve(values.size());
+    tmpTrackster.vertex_multiplicity().reserve(values.size());
+
+    for (auto const& [lc, energyScorePair] : values) {
+      if (inputClusterMask[lc.index()] > 0) {
+        tmpTrackster.vertices().push_back(lc.index());
+        double fraction = energyScorePair.first / lc->energy();
+        (*output_mask)[lc.index()] -= fraction;
+        tmpTrackster.vertex_multiplicity().push_back(static_cast<uint8_t>(std::clamp(1. / fraction, 0., 255.)));
+      }
+    }
+    tmpTrackster.setIdProbability(tracksterParticleTypeFromPdgId(sc.pdgId(), sc.charge()), 1.f);
+    tmpTrackster.setSeed(key.id(), simClusterIndex);
+    result->emplace_back(tmpTrackster);
+  }
+  ticl::assignPCAtoTracksters(
+      *result, layerClusters, layerClustersTimes, rhtools_.getPositionLayer(rhtools_.lastLayerEE(doNose_)).z());
+  evt.put(std::move(result));
+  evt.put(std::move(output_mask));
+}

RecoHGCal/TICL/python/SimTracksters_cff.py

@@ -0,0 +1,31 @@
+import FWCore.ParameterSet.Config as cms
+
+from RecoHGCal.TICL.trackstersFromSimClustersProducer_cfi import trackstersFromSimClustersProducer as _trackstersFromSimClustersProducer
+from RecoHGCal.TICL.multiClustersFromTrackstersProducer_cfi import multiClustersFromTrackstersProducer as _multiClustersFromTrackstersProducer
+from RecoHGCal.TICL.filteredLayerClustersProducer_cfi import filteredLayerClustersProducer as _filteredLayerClustersProducer
+
+
+# CA - PATTERN RECOGNITION
+
+
+filteredLayerClustersSimTracksters = _filteredLayerClustersProducer.clone(
+    clusterFilter = "ClusterFilterByAlgoAndSize",
+    algo_number = 8,
+    min_cluster_size = 0, # inclusive
+    iteration_label = "ticlSimTracksters"
+)
+
+ticlSimTracksters = _trackstersFromSimClustersProducer.clone(
+)
+
+from Configuration.ProcessModifiers.premix_stage2_cff import premix_stage2
+premix_stage2.toModify(ticlSimTracksters,
+    simclusters = "mixData:MergedCaloTruth"
+)
+
+ticlMultiClustersFromSimTracksters = _multiClustersFromTrackstersProducer.clone(
+    Tracksters = "ticlSimTracksters"
+)
+
+ticlSimTrackstersTask = cms.Task(filteredLayerClustersSimTracksters, ticlSimTracksters, ticlMultiClustersFromSimTracksters)
+

RecoHGCal/TICL/python/iterativeTICL_cff.py

@@ -5,6 +5,7 @@
 from RecoHGCal.TICL.TrkStep_cff import *
 from RecoHGCal.TICL.EMStep_cff import *
 from RecoHGCal.TICL.HADStep_cff import *
+
 from RecoHGCal.TICL.ticlLayerTileProducer_cfi import ticlLayerTileProducer
 from RecoHGCal.TICL.pfTICLProducer_cfi import pfTICLProducer as _pfTICLProducer
 from RecoHGCal.TICL.trackstersMergeProducer_cfi import trackstersMergeProducer as _trackstersMergeProducer

SimDataFormats/Associations/interface/LayerClusterToSimClusterAssociator.h

@@ -19,8 +19,11 @@ namespace hgcal {
     LayerClusterToSimClusterAssociator() = default;
     LayerClusterToSimClusterAssociator(LayerClusterToSimClusterAssociator &&) = default;
     LayerClusterToSimClusterAssociator &operator=(LayerClusterToSimClusterAssociator &&) = default;
-    ~LayerClusterToSimClusterAssociator() = default;
+    LayerClusterToSimClusterAssociator(const LayerClusterToSimClusterAssociator &) = delete;  // stop default
 
+    ~LayerClusterToSimClusterAssociator() = default;
+    const LayerClusterToSimClusterAssociator &operator=(const LayerClusterToSimClusterAssociator &) =
+        delete;  // stop default
     // ---------- const member functions ---------------------
     /// Associate a LayerCluster to SimClusters
     hgcal::RecoToSimCollectionWithSimClusters associateRecoToSim(const edm::Handle<reco::CaloClusterCollection> &cCCH,
@@ -35,11 +38,6 @@ namespace hgcal {
     }
 
   private:
-    LayerClusterToSimClusterAssociator(const LayerClusterToSimClusterAssociator &) = delete;  // stop default
-
-    const LayerClusterToSimClusterAssociator &operator=(const LayerClusterToSimClusterAssociator &) =
-        delete;  // stop default
-
     // ---------- member data --------------------------------
     std::unique_ptr<LayerClusterToSimClusterAssociatorBaseImpl> m_impl;
   };

Validation/Configuration/python/globalValidation_cff.py

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

Validation/Configuration/python/hgcalSimValid_cff.py

@@ -10,6 +10,7 @@
 from Validation.HGCalValidation.digiValidation_cff      import *
 from Validation.HGCalValidation.rechitValidation_cff    import *
 from Validation.HGCalValidation.hgcalHitValidation_cfi  import *
+from RecoHGCal.TICL.SimTracksters_cff import *
 
 from Validation.HGCalValidation.HGCalValidator_cfi import hgcalValidator
 from Validation.RecoParticleFlow.PFJetValidation_cff import pfJetValidation1 as _hgcalPFJetValidation

Validation/HGCalValidation/python/HGCalValidator_cfi.py

@@ -11,8 +11,9 @@
 from RecoHGCal.TICL.iterativeTICL_cff import ticlIterLabels, ticlIterLabelsMerge
 
 labelMcl = [cms.InputTag("ticlMultiClustersFromTracksters"+iteration) for iteration in ticlIterLabelsMerge]
+labelMcl.extend(["ticlMultiClustersFromSimTracksters"])
 lcInputMask = [cms.InputTag("ticlTracksters"+iteration) for iteration in ticlIterLabels]
-
+lcInputMask.extend(["ticlSimTracksters"])
 hgcalValidator = DQMEDAnalyzer(
     "HGCalValidator",
 

Validation/HGCalValidation/python/PostProcessorHGCAL_cfi.py

@@ -35,6 +35,7 @@
 eff_simclusters.extend(["merge_phi_layer{:02d} 'LayerCluster Merge Rate vs #phi Layer{:02d} in z-' NumMerge_LayerCluster_in_SimCluster_Phi_perlayer{:02d} Denom_LayerCluster_in_SimCluster_Phi_perlayer{:02d}".format(i, i%maxlayerzm+1, i, i) if (i<maxlayerzm) else "merge_phi_layer{:02d} 'LayerCluster Merge Rate vs #phi Layer{:02d} in z+' NumMerge_LayerCluster_in_SimCluster_Phi_perlayer{:02d} Denom_LayerCluster_in_SimCluster_Phi_perlayer{:02d}".format(i, i%maxlayerzm+1, i, i) for i in range(maxlayerzp) ])
 
 subdirsSim = ['HGCAL/HGCalValidator/simClusters/ticlTracksters'+iteration+'/' for iteration in ticlIterLabelsMerge]
+subdirsSim.extend(['HGCAL/HGCalValidator/simClusters/ticlSimTracksters'])
 postProcessorHGCALsimclusters= DQMEDHarvester('DQMGenericClient',
     subDirs = cms.untracked.vstring(subdirsSim),
     efficiency = cms.vstring(eff_simclusters),
@@ -53,7 +54,7 @@
 eff_multiclusters.extend(["merge_eta 'MultiCluster Merge Rate vs #eta' NumMerge_MultiCluster_Eta Denom_MultiCluster_Eta"])
 eff_multiclusters.extend(["merge_phi 'MultiCluster Merge Rate vs #phi' NumMerge_MultiCluster_Phi Denom_MultiCluster_Phi"])
 
-subdirsMult = ['HGCAL/HGCalValidator/hgcalMultiClusters/']
+subdirsMult = ['HGCAL/HGCalValidator/hgcalMultiClusters/','HGCAL/HGCalValidator/ticlMultiClustersFromSimTracksters/']
 subdirsMult.extend('HGCAL/HGCalValidator/ticlMultiClustersFromTracksters'+iteration+'/' for iteration in ticlIterLabelsMerge)
 
 postProcessorHGCALmulticlusters = DQMEDHarvester('DQMGenericClient',