TICL reconstruction update

DataFormats/HGCalReco/interface/Common.h

@@ -14,7 +14,20 @@ namespace ticl {
     static constexpr int nLayers = 104;
     static constexpr int iterations = 4;
     static constexpr int nBins = nEtaBins * nPhiBins;
+    static constexpr int type = 0;
   };
+
+  struct TileConstantsHFNose {
+    static constexpr float minEta = 3.0f;
+    static constexpr float maxEta = 4.2f;
+    static constexpr int nEtaBins = 24;
+    static constexpr int nPhiBins = 126;
+    static constexpr int nLayers = 16;  // 8x2
+    static constexpr int iterations = 4;
+    static constexpr int nBins = nEtaBins * nPhiBins;
+    static constexpr int type = 1;
+  };
+
 }  // namespace ticl
 
 namespace ticl {

DataFormats/HGCalReco/interface/TICLLayerTile.h

@@ -14,6 +14,8 @@ class TICLLayerTileT {
     tile_[globalBin(eta, phi)].push_back(layerClusterId);
   }
 
+  int typeT() const { return T::type; }
+
   int etaBin(float eta) const {
     constexpr float etaRange = T::maxEta - T::minEta;
     static_assert(etaRange >= 0.f);
@@ -31,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; }
@@ -51,6 +61,11 @@ namespace ticl {
   using TICLLayerTile = TICLLayerTileT<TileConstants>;
   using Tiles = std::array<TICLLayerTile, TileConstants::nLayers>;
   using TracksterTiles = std::array<TICLLayerTile, TileConstants::iterations>;
+
+  using TICLLayerTileHFNose = TICLLayerTileT<TileConstantsHFNose>;
+  using TilesHFNose = std::array<TICLLayerTileHFNose, TileConstantsHFNose::nLayers>;
+  using TracksterTilesHFNose = std::array<TICLLayerTileHFNose, TileConstantsHFNose::iterations>;
+
 }  // namespace ticl
 
 template <typename T>
@@ -68,5 +83,7 @@ class TICLGenericTile {
 
 using TICLLayerTiles = TICLGenericTile<ticl::Tiles>;
 using TICLTracksterTiles = TICLGenericTile<ticl::TracksterTiles>;
+using TICLLayerTilesHFNose = TICLGenericTile<ticl::TilesHFNose>;
+using TICLTracksterTilesHFNose = TICLGenericTile<ticl::TracksterTilesHFNose>;
 
 #endif

DataFormats/HGCalReco/src/classes_def.xml

@@ -25,6 +25,20 @@
   </class>
   <class name="edm::Wrapper<TICLGenericTile<ticl::TracksterTiles>>" persistent="false"/>
 
+  <class name="ticl::TICLLayerTileHFNose" persistent="false">
+    <field name="tile_" transient="true"/>
+  </class>
+
+  <class name="TICLGenericTile<ticl::TilesHFNose>" persistent="false">
+    <field name="tiles_" transient="true"/>
+  </class>
+  <class name="edm::Wrapper<TICLGenericTile<ticl::TilesHFNose>>" persistent="false"/>
+
+  <class name="TICLGenericTile<ticl::TracksterTilesHFNose>" persistent="false">
+    <field name="tiles_" transient="true"/>
+  </class>
+  <class name="edm::Wrapper<TICLGenericTile<ticl::TracksterTilesHFNose>>" persistent="false"/>
+
   <class name="TICLSeedingRegion" ClassVersion="3">
     <version ClassVersion="3" checksum="2409655320"/>
   </class>

RecoEgamma/EgammaTools/interface/HGCalEgammaIDHelper.h

@@ -14,6 +14,7 @@
 #include "FWCore/Framework/interface/ESHandle.h"
 #include "FWCore/Framework/interface/Event.h"
 #include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Utilities/interface/ESGetToken.h"
 #include "FWCore/Framework/interface/ConsumesCollector.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 
@@ -78,6 +79,7 @@ class HGCalEgammaIDHelper {
   edm::EDGetTokenT<HGCRecHitCollection> recHitsEE_;
   edm::EDGetTokenT<HGCRecHitCollection> recHitsFH_;
   edm::EDGetTokenT<HGCRecHitCollection> recHitsBH_;
+  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometry_;
   hgcal::RecHitTools recHitTools_;
   bool debug_;
 };

RecoEgamma/EgammaTools/src/HGCalEgammaIDHelper.cc

@@ -14,6 +14,7 @@ HGCalEgammaIDHelper::HGCalEgammaIDHelper(const edm::ParameterSet& iConfig, edm::
   recHitsEE_ = iC.consumes<HGCRecHitCollection>(eeRecHitInputTag_);
   recHitsFH_ = iC.consumes<HGCRecHitCollection>(fhRecHitInputTag_);
   recHitsBH_ = iC.consumes<HGCRecHitCollection>(bhRecHitInputTag_);
+  caloGeometry_ = iC.esConsumes<CaloGeometry, CaloGeometryRecord>();
   pcaHelper_.setdEdXWeights(dEdXWeights_);
   debug_ = iConfig.getUntrackedParameter<bool>("debug", false);
 }
@@ -26,7 +27,8 @@ void HGCalEgammaIDHelper::eventInit(const edm::Event& iEvent, const edm::EventSe
   edm::Handle<HGCRecHitCollection> recHitHandleBH;
   iEvent.getByToken(recHitsBH_, recHitHandleBH);
 
-  recHitTools_.getEventSetup(iSetup);
+  edm::ESHandle<CaloGeometry> geom = iSetup.getHandle(caloGeometry_);
+  recHitTools_.setGeometry(*geom);
   pcaHelper_.setRecHitTools(&recHitTools_);
   isoHelper_.setRecHitTools(&recHitTools_);
   pcaHelper_.fillHitMap(*recHitHandleEE, *recHitHandleFH, *recHitHandleBH);

RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py

@@ -7,6 +7,7 @@
       'keep *_ticlMultiClustersFromTrackstersEM_*_*',
       'keep *_ticlMultiClustersFromTrackstersHAD_*_*',
       'keep *_ticlMultiClustersFromTrackstersTrk_*_*',
+      'keep *_ticlMultiClustersFromTrackstersTrkEM_*_*',
       'keep *_ticlMultiClustersFromTrackstersMIP_*_*',
       'keep *_ticlMultiClustersFromTrackstersMerge_*_*',
       )
@@ -15,12 +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 *_pfTICL_*_*'
       )
     )

RecoHGCal/TICL/plugins/ClusterFilterByAlgoAndSizeAndLayerRange.h

@@ -0,0 +1,55 @@
+// Authors: Marco Rovere - marco.rovere@cern.ch, Felice Pantaleo - felice.pantaleo@cern.ch
+// 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

@@ -10,6 +10,7 @@
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Utilities/interface/ESGetToken.h"
 
 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
@@ -30,6 +31,7 @@ class FilteredLayerClustersProducer : public edm::stream::EDProducer<> {
 private:
   edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
   edm::EDGetTokenT<std::vector<float>> clustersMask_token_;
+  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometry_token_;
   std::string clusterFilter_;
   std::string iteration_label_;
   std::unique_ptr<const ticl::ClusterFilterBase> theFilter_;
@@ -39,26 +41,32 @@ class FilteredLayerClustersProducer : public edm::stream::EDProducer<> {
 DEFINE_FWK_MODULE(FilteredLayerClustersProducer);
 
 FilteredLayerClustersProducer::FilteredLayerClustersProducer(const edm::ParameterSet& ps) {
-  clusters_token_ = consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("HGCLayerClusters"));
+  clusters_token_ = consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("LayerClusters"));
   clustersMask_token_ = consumes<std::vector<float>>(ps.getParameter<edm::InputTag>("LayerClustersInputMask"));
+  caloGeometry_token_ = esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>();
   clusterFilter_ = ps.getParameter<std::string>("clusterFilter");
   theFilter_ = ClusterFilterFactory::get()->create(clusterFilter_, ps);
   iteration_label_ = ps.getParameter<std::string>("iteration_label");
   produces<std::vector<float>>(iteration_label_);
 }
 
-void FilteredLayerClustersProducer::beginRun(edm::Run const&, edm::EventSetup const& es) { rhtools_.getEventSetup(es); }
+void FilteredLayerClustersProducer::beginRun(edm::Run const&, edm::EventSetup const& es) {
+  edm::ESHandle<CaloGeometry> geom = es.getHandle(caloGeometry_token_);
+  rhtools_.setGeometry(*geom);
+}
 
 void FilteredLayerClustersProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   // hgcalMultiClusters
   edm::ParameterSetDescription desc;
-  desc.add<edm::InputTag>("HGCLayerClusters", edm::InputTag("hgcalLayerClusters"));
+  desc.add<edm::InputTag>("LayerClusters", edm::InputTag("hgcalLayerClusters"));
   desc.add<edm::InputTag>("LayerClustersInputMask", edm::InputTag("hgcalLayerClusters", "InitialLayerClustersMask"));
   desc.add<std::string>("iteration_label", "iterationLabelGoesHere");
   desc.add<std::string>("clusterFilter", "ClusterFilterByAlgoAndSize");
   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/HGCDoublet.cc

@@ -81,14 +81,17 @@ int HGCDoublet::areAligned(double xi,
 
   // inner product
   auto dot = dx1 * dx2 + dy1 * dy2 + dz1 * dz2;
+  auto dotsq = dot * dot;
   // magnitudes
-  auto mag1 = std::sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
-  auto mag2 = std::sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
-  // angle between the vectors
-  auto cosTheta = dot / (mag1 * mag2);
+  auto mag1sq = dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
+  auto mag2sq = dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
+
+  auto minCosTheta_sq = minCosTheta * minCosTheta;
+  bool isWithinLimits = (dotsq > minCosTheta_sq * mag1sq * mag2sq);
+
   if (debug) {
-    LogDebug("HGCDoublet") << "-- Are Aligned -- dot: " << dot << " mag1: " << mag1 << " mag2: " << mag2
-                           << " cosTheta: " << cosTheta << " isWithinLimits: " << (cosTheta > minCosTheta) << std::endl;
+    LogDebug("HGCDoublet") << "-- Are Aligned -- dotsq: " << dotsq << " mag1sq: " << mag1sq << " mag2sq: " << mag2sq
+                           << "minCosTheta_sq:" << minCosTheta_sq << " isWithinLimits: " << isWithinLimits << std::endl;
   }
 
   // Now check the compatibility with the pointing origin.
@@ -102,15 +105,18 @@ int HGCDoublet::areAligned(double xi,
   const GlobalVector pointingDir = (seedIndex_ == -1) ? GlobalVector(innerX(), innerY(), innerZ()) : refDir;
 
   auto dot_pointing = pointingDir.dot(firstDoublet);
-  auto mag_pointing = sqrt(pointingDir.mag2());
-  auto cosTheta_pointing = dot_pointing / (mag2 * mag_pointing);
+  auto dot_pointing_sq = dot_pointing * dot_pointing;
+  auto mag_pointing_sq = pointingDir.mag2();
+  auto minCosPointing_sq = minCosPointing * minCosPointing;
+  bool isWithinLimitsPointing = (dot_pointing_sq > minCosPointing_sq * mag_pointing_sq * mag2sq);
   if (debug) {
-    LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing: " << dot_pointing << " mag_pointing: " << mag_pointing
-                           << " mag2: " << mag2 << " cosTheta_pointing: " << cosTheta_pointing
-                           << " isWithinLimits: " << (cosTheta_pointing > minCosPointing) << std::endl;
+    LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing_sq: " << dot_pointing * dot_pointing
+                           << " mag_pointing_sq: " << mag_pointing_sq << " mag2sq: " << mag2sq
+                           << " isWithinLimits: " << isWithinLimitsPointing << std::endl;
   }
-
-  return (cosTheta > minCosTheta) && (cosTheta_pointing > minCosPointing);
+  // by squaring cosTheta and multiplying by the squares of the magnitudes
+  // an equivalent comparison is made without the division and square root which are costly FP ops.
+  return isWithinLimits && isWithinLimitsPointing;
 }
 
 void HGCDoublet::findNtuplets(std::vector<HGCDoublet> &allDoublets,