[HGCal] TICL v3 major upgrade

DataFormats/HGCalReco/interface/TICLLayerTile.h

@@ -33,6 +33,14 @@ class TICLLayerTileT {
     return phiBin;
   }
 
+  std::array<int, 4> searchBoxEtaPhi(float etaMin, float etaMax, float phiMin, float phiMax) const {
+    int etaBinMin = etaBin(etaMin);
+    int etaBinMax = etaBin(etaMax);
+    int phiBinMin = phiBin(phiMin);
+    int phiBinMax = phiBin(phiMax);
+    return std::array<int, 4>({{etaBinMin, etaBinMax, phiBinMin, phiBinMax}});
+  }
+
   int globalBin(int etaBin, int phiBin) const { return phiBin + etaBin * T::nPhiBins; }
 
   int globalBin(double eta, double phi) const { return phiBin(phi) + etaBin(eta) * T::nPhiBins; }

RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py

@@ -7,6 +7,7 @@
       'keep *_ticlMultiClustersFromTrackstersEM_*_*',
       'keep *_ticlMultiClustersFromTrackstersHAD_*_*',
       'keep *_ticlMultiClustersFromTrackstersTrk_*_*',
+      'keep *_ticlMultiClustersFromTrackstersTrkEM_*_*',
       'keep *_ticlMultiClustersFromTrackstersMIP_*_*',
       'keep *_ticlMultiClustersFromTrackstersMerge_*_*',
       )
@@ -15,16 +16,15 @@
 #RECO content
 TICL_RECO = cms.PSet(
     outputCommands = cms.untracked.vstring(
+      'keep *_ticlTrackstersTrkEM_*_*',
       'keep *_ticlTrackstersEM_*_*',
       'keep *_ticlTrackstersHAD_*_*',
       'keep *_ticlTrackstersTrk_*_*',
       'keep *_ticlTrackstersMIP_*_*',
       'keep *_ticlTrackstersMerge_*_*',
-      'keep *_ticlCandidateFromTracksters_*_*',
       'keep *_ticlTrackstersHFNoseEM_*_*',
       'keep *_ticlTrackstersHFNoseMIP_*_*',
       'keep *_ticlTrackstersHFNoseMerge_*_*',
-      'keep *_ticlCandidateFromTrackstersHFNose_*_*',
       'keep *_pfTICL_*_*'
       )
     )

RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h

@@ -0,0 +1,55 @@
+// Authors: Marco Rovere - [email protected], Felice Pantaleo - [email protected]
+// Date: 09/2020
+
+#ifndef RecoHGCal_TICL_ClusterFilterByAlgoAndSizeAndLayerRange_H__
+#define RecoHGCal_TICL_ClusterFilterByAlgoAndSizeAndLayerRange_H__
+
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "ClusterFilterBase.h"
+
+#include <memory>
+#include <utility>
+
+// Filter clusters that belong to a specific algorithm
+namespace ticl {
+  class ClusterFilterByAlgoAndSizeAndLayerRange final : public ClusterFilterBase {
+  public:
+    ClusterFilterByAlgoAndSizeAndLayerRange(const edm::ParameterSet& ps)
+        : ClusterFilterBase(ps),
+          algo_number_(ps.getParameter<int>("algo_number")),
+          min_cluster_size_(ps.getParameter<int>("min_cluster_size")),
+          max_cluster_size_(ps.getParameter<int>("max_cluster_size")),
+          min_layerId_(ps.getParameter<int>("min_layerId")),
+          max_layerId_(ps.getParameter<int>("max_layerId")) {}
+    ~ClusterFilterByAlgoAndSizeAndLayerRange() override{};
+
+    void filter(const std::vector<reco::CaloCluster>& layerClusters,
+                const HgcalClusterFilterMask& availableLayerClusters,
+                std::vector<float>& layerClustersMask,
+                hgcal::RecHitTools& rhtools) const override {
+      auto filteredLayerClusters = std::make_unique<HgcalClusterFilterMask>();
+      for (auto const& cl : availableLayerClusters) {
+        auto const& layerCluster = layerClusters[cl.first];
+        auto const& haf = layerCluster.hitsAndFractions();
+        auto layerId = rhtools.getLayerWithOffset(haf[0].first);
+
+        if (layerCluster.algo() == algo_number_ and layerId <= max_layerId_ and layerId >= min_layerId_ and
+            haf.size() <= max_cluster_size_ and
+            (haf.size() >= min_cluster_size_ or !(rhtools.isSilicon(haf[0].first)))) {
+          filteredLayerClusters->emplace_back(cl);
+        } else {
+          layerClustersMask[cl.first] = 0.;
+        }
+      }
+    }
+
+  private:
+    int algo_number_;
+    unsigned int min_cluster_size_;
+    unsigned int max_cluster_size_;
+    unsigned int min_layerId_;
+    unsigned int max_layerId_;
+  };
+}  // namespace ticl
+
+#endif

RecoHGCal/TICL/plugins/FilteredLayerClustersProducer.cc

@@ -65,6 +65,8 @@ void FilteredLayerClustersProducer::fillDescriptions(edm::ConfigurationDescripti
   desc.add<int>("algo_number", 9);
   desc.add<int>("min_cluster_size", 0);
   desc.add<int>("max_cluster_size", 9999);
+  desc.add<int>("min_layerId", 0);
+  desc.add<int>("max_layerId", 9999);
   descriptions.add("filteredLayerClustersProducer", desc);
 }
 

RecoHGCal/TICL/plugins/HGCGraph.cc

@@ -19,14 +19,18 @@ void HGCGraphT<TILES>::makeAndConnectDoublets(const TILES &histo,
                                               int deltaIPhi,
                                               float minCosTheta,
                                               float minCosPointing,
+                                              float root_doublet_max_distance_from_seed_squared,
                                               float etaLimitIncreaseWindow,
-                                              int missing_layers,
+                                              int skip_layers,
                                               int maxNumberOfLayers,
                                               float maxDeltaTime) {
   isOuterClusterOfDoublets_.clear();
   isOuterClusterOfDoublets_.resize(layerClusters.size());
   allDoublets_.clear();
   theRootDoublets_.clear();
+  bool checkDistanceRootDoubletVsSeed = root_doublet_max_distance_from_seed_squared < 9999;
+  float origin_eta;
+  float origin_phi;
   for (const auto &r : regions) {
     bool isGlobal = (r.index == -1);
     auto zSide = r.zSide;
@@ -37,19 +41,22 @@ void HGCGraphT<TILES>::makeAndConnectDoublets(const TILES &histo,
       startPhiBin = 0;
       endEtaBin = nEtaBins;
       endPhiBin = nPhiBins;
+      origin_eta = 0;
+      origin_phi = 0;
     } else {
       auto firstLayerOnZSide = maxNumberOfLayers * zSide;
       const auto &firstLayerHisto = histo[firstLayerOnZSide];
-
-      int entryEtaBin = firstLayerHisto.etaBin(r.origin.eta());
-      int entryPhiBin = firstLayerHisto.phiBin(r.origin.phi());
+      origin_eta = r.origin.eta();
+      origin_phi = r.origin.phi();
+      int entryEtaBin = firstLayerHisto.etaBin(origin_eta);
+      int entryPhiBin = firstLayerHisto.phiBin(origin_phi);
       // For track-seeded iterations, if the impact point is below a certain
       // eta-threshold, i.e., it has higher eta, make the initial search
       // window bigger in both eta and phi by one bin, to contain better low
       // energy showers.
       auto etaWindow = deltaIEta;
       auto phiWindow = deltaIPhi;
-      if (std::abs(r.origin.eta()) > etaLimitIncreaseWindow) {
+      if (std::abs(origin_eta) > etaLimitIncreaseWindow) {
         etaWindow++;
         phiWindow++;
         LogDebug("HGCGraph") << "Limit of Eta for increase: " << etaLimitIncreaseWindow
@@ -60,9 +67,9 @@ void HGCGraphT<TILES>::makeAndConnectDoublets(const TILES &histo,
       startPhiBin = entryPhiBin - phiWindow;
       endPhiBin = entryPhiBin + phiWindow + 1;
       if (verbosity_ > Guru) {
-        LogDebug("HGCGraph") << " Entrance eta, phi: " << r.origin.eta() << ", " << r.origin.phi()
+        LogDebug("HGCGraph") << " Entrance eta, phi: " << origin_eta << ", " << origin_phi
                              << " entryEtaBin: " << entryEtaBin << " entryPhiBin: " << entryPhiBin
-                             << " globalBin: " << firstLayerHisto.globalBin(r.origin.eta(), r.origin.phi())
+                             << " globalBin: " << firstLayerHisto.globalBin(origin_eta, origin_phi)
                              << " on layer: " << firstLayerOnZSide << " startEtaBin: " << startEtaBin
                              << " endEtaBin: " << endEtaBin << " startPhiBin: " << startPhiBin
                              << " endPhiBin: " << endPhiBin << " phiBin(0): " << firstLayerHisto.phiBin(0.)
@@ -75,7 +82,7 @@ void HGCGraphT<TILES>::makeAndConnectDoublets(const TILES &histo,
     }
 
     for (int il = 0; il < maxNumberOfLayers - 1; ++il) {
-      for (int outer_layer = 0; outer_layer < std::min(1 + missing_layers, maxNumberOfLayers - 1 - il); ++outer_layer) {
+      for (int outer_layer = 0; outer_layer < std::min(1 + skip_layers, maxNumberOfLayers - 1 - il); ++outer_layer) {
         int currentInnerLayerId = il + maxNumberOfLayers * zSide;
         int currentOuterLayerId = currentInnerLayerId + 1 + outer_layer;
         auto const &outerLayerHisto = histo[currentOuterLayerId];
@@ -171,8 +178,17 @@ void HGCGraphT<TILES>::makeAndConnectDoublets(const TILES &histo,
                                                                               minCosTheta,
                                                                               minCosPointing,
                                                                               verbosity_ > Advanced);
-                    if (isRootDoublet)
+                    if (isRootDoublet and checkDistanceRootDoubletVsSeed) {
+                      auto dEtaSquared = (layerClusters[innerClusterId].eta() - origin_eta);
+                      dEtaSquared *= dEtaSquared;
+                      auto dPhiSquared = (layerClusters[innerClusterId].phi() - origin_phi);
+                      dPhiSquared *= dPhiSquared;
+                      if (dEtaSquared + dPhiSquared > root_doublet_max_distance_from_seed_squared)
+                        isRootDoublet = false;
+                    }
+                    if (isRootDoublet) {
                       theRootDoublets_.push_back(doubletId);
+                    }
                   }
                 }
               }

RecoHGCal/TICL/plugins/HGCGraph.h

@@ -25,8 +25,9 @@ class HGCGraphT {
                               int deltaIPhi,
                               float minCosThetai,
                               float maxCosPointing,
+                              float root_doublet_max_distance_from_seed_squared,
                               float etaLimitIncreaseWindow,
-                              int missing_layers,
+                              int skip_layers,
                               int maxNumberOfLayers,
                               float maxDeltaTime);
 

RecoHGCal/TICL/plugins/MultiClustersFromTrackstersProducer.cc

@@ -63,13 +63,14 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev
   }
 
   std::for_each(std::begin(tracksters), std::end(tracksters), [&](auto const& trackster) {
-    // Do not create a multicluster if the trackster has no layer clusters.
+    // Create an empty multicluster if the trackster has no layer clusters.
     // This could happen when a seed leads to no trackster and a dummy one is produced.
+
+    std::array<double, 3> baricenter{{0., 0., 0.}};
+    double total_weight = 0.;
+    reco::HGCalMultiCluster temp;
+    int counter = 0;
     if (!trackster.vertices().empty()) {
-      std::array<double, 3> baricenter{{0., 0., 0.}};
-      double total_weight = 0.;
-      reco::HGCalMultiCluster temp;
-      int counter = 0;
       std::for_each(std::begin(trackster.vertices()), std::end(trackster.vertices()), [&](unsigned int idx) {
         temp.push_back(clusterPtrs[idx]);
         auto fraction = 1.f / trackster.vertex_multiplicity(counter++);
@@ -86,13 +87,13 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev
           std::begin(baricenter), std::end(baricenter), std::begin(baricenter), [&total_weight](double val) -> double {
             return val / total_weight;
           });
-      temp.setEnergy(total_weight);
-      temp.setCorrectedEnergy(total_weight);
-      temp.setPosition(math::XYZPoint(baricenter[0], baricenter[1], baricenter[2]));
-      temp.setAlgoId(reco::CaloCluster::hgcal_em);
-      temp.setTime(trackster.time(), trackster.timeError());
-      multiclusters->push_back(temp);
     }
+    temp.setEnergy(total_weight);
+    temp.setCorrectedEnergy(total_weight);
+    temp.setPosition(math::XYZPoint(baricenter[0], baricenter[1], baricenter[2]));
+    temp.setAlgoId(reco::CaloCluster::hgcal_em);
+    temp.setTime(trackster.time(), trackster.timeError());
+    multiclusters->push_back(temp);
   });
 
   evt.put(std::move(multiclusters));

RecoHGCal/TICL/plugins/PFTICLProducer.cc

@@ -36,7 +36,7 @@ PFTICLProducer::PFTICLProducer(const edm::ParameterSet& conf)
 
 void PFTICLProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
-  desc.add<edm::InputTag>("ticlCandidateSrc", edm::InputTag("ticlCandidateFromTracksters"));
+  desc.add<edm::InputTag>("ticlCandidateSrc", edm::InputTag("ticlTrackstersMerge"));
   descriptions.add("pfTICLProducer", desc);
 }
 
@@ -52,6 +52,15 @@ void PFTICLProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSet
     const auto abs_pdg_id = std::abs(ticl_cand.pdgId());
     const auto charge = ticl_cand.charge();
     const auto& four_mom = ticl_cand.p4();
+    double ecal_energy = 0.;
+
+    for (const auto& t : ticl_cand.tracksters()) {
+      double ecal_energy_fraction = t->raw_em_pt() / t->raw_pt();
+      ecal_energy += t->raw_energy() * ecal_energy_fraction;
+    }
+    double hcal_energy = ticl_cand.rawEnergy() - ecal_energy;
+    // fix for floating point rounding could go slightly below 0
+    hcal_energy = hcal_energy < 0 ? 0 : hcal_energy;
 
     reco::PFCandidate::ParticleType part_type;
     switch (abs_pdg_id) {
@@ -78,6 +87,8 @@ void PFTICLProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSet
     candidates->emplace_back(charge, four_mom, part_type);
 
     auto& candidate = candidates->back();
+    candidate.setEcalEnergy(ecal_energy, ecal_energy);
+    candidate.setHcalEnergy(hcal_energy, hcal_energy);
     if (candidate.charge()) {  // otherwise PFCandidate throws
       // Construct edm::Ref from edm::Ptr. As of now, assumes type to be reco::Track. To be extended (either via
       // dynamic type checking or configuration) if additional track types are needed.