Reorganization of kdtrees for PFBlockAlgo and optimize track-hcal links

DataFormats/ParticleFlowReco/interface/PFBlockElement.h

@@ -111,12 +111,17 @@ namespace reco {
     virtual bool isLinkedToDisplacedVertex() const { return false; }
 
     // Glowinski & Gouzevitch
-    void setMultilinks(const PFMultiLinksTC& ml) { multilinks_ = ml; }
-    void setIsValidMultilinks(bool isVal) { multilinks_.isValid = isVal; }
-    void setMultilinksList(const PFMultilinksType& links) { multilinks_.linkedClusters = links; }
-
-    bool isMultilinksValide() const { return multilinks_.isValid; }
-    const PFMultilinksType& getMultilinks() const { return multilinks_.linkedClusters; }
+    void setMultilinks(const PFMultiLinksTC& ml, Type type) { multilinks_[type] = ml; }
+    void setIsValidMultilinks(bool isVal, Type type) { multilinks_[type].isValid = isVal; }
+
+    bool isMultilinksValide(Type type) const {
+      const auto& it = multilinks_.find(type);
+      if (it != multilinks_.end())
+        return it->second.isValid;
+      else
+        return false;  // no multilinks_ for the specified type
+    }
+    const PFMultilinksType& getMultilinks(Type type) const { return multilinks_.at(type).linkedClusters; }
     // ! Glowinski & Gouzevitch
 
     /// do we have a valid time information
@@ -145,7 +150,8 @@ namespace reco {
     unsigned index_;
 
     // Glowinski & Gouzevitch
-    PFMultiLinksTC multilinks_;
+    // PFMultiLinks for each different link target type
+    std::map<reco::PFBlockElement::Type, PFMultiLinksTC> multilinks_;
     // ! Glowinski & Gouzevitch
 
     /// timing information (valid if timeError_ >= 0)

DataFormats/ParticleFlowReco/src/classes_def_2.xml

@@ -102,8 +102,6 @@
    <version ClassVersion="12" checksum="2850643502"/>
    <version ClassVersion="11" checksum="2924298092"/>
    <version ClassVersion="10" checksum="2979491836"/>
-    <!-- <field name="doPropagation_" transient="true"/> -->
-    <!-- <field name="color_" transient="true"/> -->
   </class>
   <class name="std::vector<reco::PFRecTrack>"/>
   <class name="edm::Wrapper<std::vector<reco::PFRecTrack> >"/>
@@ -113,8 +111,6 @@
    <version ClassVersion="12" checksum="2203231294"/>
    <version ClassVersion="11" checksum="2592328732"/>
    <version ClassVersion="10" checksum="3411096264"/>
-    <!-- <field name="doPropagation_" transient="true"/> -->
-    <!-- <field name="color_" transient="true"/> -->
   </class>
   <class name="std::vector<reco::GsfPFRecTrack>"/>
   <class name="edm::Wrapper<std::vector<reco::GsfPFRecTrack> >"/>
@@ -136,7 +132,8 @@
   <class name="std::vector<reco::PFSimParticle>"/>
   <class name="edm::Wrapper<std::vector<reco::PFSimParticle> >"/>
 
-  <class name="reco::PFBlockElement" ClassVersion="13">
+  <class name="reco::PFBlockElement" ClassVersion="14">
+   <version ClassVersion="14" checksum="3982859549"/>
    <version ClassVersion="13" checksum="4257111441"/>
    <version ClassVersion="12" checksum="3101496393"/>
    <version ClassVersion="11" checksum="3944635097"/>
@@ -149,11 +146,13 @@
    <field name="nullMuon_" transient="true"/>
    <field name="nullVertex_" transient="true"/>
   </class>
+  <class name="std::map<reco::PFBlockElement::Type,reco::PFMultiLinksTC>" />
   <class name="std::vector<reco::PFBlockElement *>" />
   <class name="edm::OwnVector<reco::PFBlockElement, edm::ClonePolicy<reco::PFBlockElement> >"/>
   <class name="edm::Wrapper<edm::OwnVector<reco::PFBlockElement, edm::ClonePolicy<reco::PFBlockElement> > >" />
 
-  <class name="reco::PFBlockElementTrack" ClassVersion="14">
+  <class name="reco::PFBlockElementTrack" ClassVersion="15">
+   <version ClassVersion="15" checksum="1577934967"/>
    <version ClassVersion="14" checksum="3952424659"/>
    <version ClassVersion="13" checksum="1449646331"/>
    <version ClassVersion="12" checksum="417575083"/>
@@ -164,7 +163,8 @@
 -->
   </class>
 
-  <class name="reco::PFBlockElementGsfTrack" ClassVersion="13">
+  <class name="reco::PFBlockElementGsfTrack" ClassVersion="14">
+   <version ClassVersion="14" checksum="3921871079"/>
    <version ClassVersion="13" checksum="2453155267"/>
    <version ClassVersion="12" checksum="1521047275"/>
    <version ClassVersion="11" checksum="1367777435"/>
@@ -175,7 +175,8 @@
 -->
   </class>
 
-  <class name="reco::PFBlockElementBrem" ClassVersion="13">
+  <class name="reco::PFBlockElementBrem" ClassVersion="14">
+   <version ClassVersion="14" checksum="3702597792"/>
    <version ClassVersion="13" checksum="157699540"/>
    <version ClassVersion="12" checksum="112730380"/>
    <version ClassVersion="11" checksum="2362699420"/>
@@ -187,7 +188,8 @@
   </class>
 
 
-  <class name="reco::PFBlockElementCluster" ClassVersion="14">
+  <class name="reco::PFBlockElementCluster" ClassVersion="15">
+   <version ClassVersion="15" checksum="3252388"/>
    <version ClassVersion="14" checksum="407383072"/>
    <version ClassVersion="13" checksum="775826696"/>
    <version ClassVersion="12" checksum="1263691576"/>
@@ -265,7 +267,8 @@
   <class name="edm::Wrapper<edm::ValueMap<reco::PreIdRef> >"/>
   <class name="edm::Ref<std::vector<reco::PreId>,reco::PreId,edm::refhelper::FindUsingAdvance<std::vector<reco::PreId>,reco::PreId> >" />
 
-  <class name="reco::PFBlockElementSuperCluster" ClassVersion="14">
+  <class name="reco::PFBlockElementSuperCluster" ClassVersion="15">
+   <version ClassVersion="15" checksum="1074027768"/>
    <version ClassVersion="14" checksum="285242772"/>
    <version ClassVersion="13" checksum="3749373244"/>
    <version ClassVersion="12" checksum="1656578540"/>

HLTrigger/Configuration/python/customizeHLTforCMSSW.py

@@ -170,10 +170,27 @@ def customiseFor2017DtUnpacking(process):
 
     return process
 
+def customiseFor31295(process):
+    """Reorganization of kdtrees for PFBlockAlgo and optimize track-hcal links"""
+
+    # for PFBlockProducer
+    for producer in producers_by_type(process, "PFBlockProducer"):
+        if hasattr(producer,'linkDefinitions'):
+            for ps in producer.linkDefinitions.value():
+                if hasattr(ps,'linkerName') and (ps.linkerName == 'TrackAndHCALLinker'):
+                    if not hasattr(ps,'nMaxHcalLinksPerTrack'):
+                        ps.nMaxHcalLinksPerTrack = cms.int32(1)
+                if hasattr(ps,'linkerName') and (ps.linkerName == 'ECALAndHCALLinker'):
+                    if not hasattr(ps,'minAbsEtaEcal'):
+                        ps.minAbsEtaEcal = cms.double(2.5)
+
+    return process
+
 # CMSSW version specific customizations
 def customizeHLTforCMSSW(process, menuType="GRun"):
 
     # add call to action function in proper order: newest last!
     # process = customiseFor12718(process)
+    process = customiseFor31295(process)
 
     return process

RecoParticleFlow/PFClusterProducer/plugins/Basic2DClusterForEachSeed.cc

@@ -31,7 +31,6 @@ void Basic2DClusterForEachSeed::buildClusters(const edm::Handle<reco::PFRecHitCo
     cluster.setPosition(math::XYZPoint(refhit->position().x(), refhit->position().y(), refhit->position().z()));
     cluster.calculatePositionREP();
     cluster.setDepth(refhit->depth());
-    cluster.calculatePositionREP();
 
     output.push_back(cluster);
 

RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc

@@ -8,6 +8,8 @@ using Vector2D = Basic2DVector<double>::MathVector;
 namespace {
   const Vector2D one2D = BVector2D(1.0, 1.0).v;
   const Vector2D fivepercent2D = BVector2D(0.05, 0.05).v;
+  // rechit with fraction below this value will be ignored in LinkByRecHit
+  const float minPFRecHitFrac = 1E-4;
 }  // namespace
 
 // to enable debugs
@@ -216,7 +218,7 @@ double LinkByRecHit::testTrackAndClusterByRecHit(const reco::PFRecTrack& track,
   for (unsigned int rhit = 0; rhit < fracs.size(); ++rhit) {
     const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
     double fraction = fracs[rhit].fraction();
-    if (fraction < 1E-4)
+    if (fraction < minPFRecHitFrac)
       continue;
     if (rh.isNull())
       continue;
@@ -415,7 +417,7 @@ double LinkByRecHit::testECALAndPSByRecHit(const reco::PFCluster& clusterECAL,
   for (unsigned int rhit = 0; rhit < fracs.size(); ++rhit) {
     const auto& rh = fracs[rhit].recHitRef();
     double fraction = fracs[rhit].fraction();
-    if (fraction < 1E-4)
+    if (fraction < minPFRecHitFrac)
       continue;
     if (rh.isNull())
       continue;
@@ -496,13 +498,13 @@ double LinkByRecHit::testHFEMAndHFHADByRecHit(const reco::PFCluster& clusterHFEM
   const auto& posxyzEM = clusterHFEM.position();
   const auto& posxyzHAD = clusterHFHAD.position();
 
-  double dX = posxyzEM.X() - posxyzHAD.X();
-  double dY = posxyzEM.Y() - posxyzHAD.Y();
   double sameZ = posxyzEM.Z() * posxyzHAD.Z();
-
   if (sameZ < 0)
     return -1.;
 
+  double dX = posxyzEM.X() - posxyzHAD.X();
+  double dY = posxyzEM.Y() - posxyzHAD.Y();
+
   double dist2 = dX * dX + dY * dY;
 
   if (dist2 < 0.1) {

RecoParticleFlow/PFProducer/interface/KDTreeLinkerBase.h

@@ -88,6 +88,16 @@ class KDTreeLinkerBase {
 
   // Debug boolean. Not used until now.
   bool debug_ = false;
+
+  // Sorted indexes
+  template <typename T>
+  static std::vector<size_t> sort_indexes(const std::vector<T> &v) {
+    std::vector<size_t> idx(v.size());
+    for (size_t i = 0; i != idx.size(); ++i)
+      idx[i] = i;
+    std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return v[i1] < v[i2]; });
+    return idx;
+  }
 };
 
 #include "FWCore/PluginManager/interface/PluginFactory.h"

RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerPSEcal.cc

@@ -143,8 +143,6 @@ void KDTreeLinkerPSEcal::searchLinks() {
 
   // We iterate over the PS clusters.
   for (BlockEltSet::iterator it = targetSet_.begin(); it != targetSet_.end(); it++) {
-    (*it)->setIsValidMultilinks(true);
-
     reco::PFClusterRef clusterPSRef = (*it)->clusterRef();
     const reco::PFCluster &clusterPS = *clusterPSRef;
 
@@ -233,18 +231,25 @@ void KDTreeLinkerPSEcal::searchLinks() {
 void KDTreeLinkerPSEcal::updatePFBlockEltWithLinks() {
   //TODO YG : Check if cluster positionREP() is valid ?
 
-  // Here we save in each track the list of phi/eta values of linked clusters.
+  // Here we save in each PS the list of phi/eta values of linked ECAL clusters.
   for (BlockElt2BlockEltMap::iterator it = target2ClusterLinks_.begin(); it != target2ClusterLinks_.end(); ++it) {
+    const auto &psElt = it->first;
+    const auto &ecalEltSet = it->second;
     reco::PFMultiLinksTC multitracks(true);
 
-    for (BlockEltSet::iterator jt = it->second.begin(); jt != it->second.end(); ++jt) {
-      double clusterphi = (*jt)->clusterRef()->positionREP().phi();
-      double clustereta = (*jt)->clusterRef()->positionREP().eta();
+    for (const auto &ecalElt : ecalEltSet) {
+      double clusterphi = ecalElt->clusterRef()->positionREP().phi();
+      double clustereta = ecalElt->clusterRef()->positionREP().eta();
 
       multitracks.linkedClusters.push_back(std::make_pair(clusterphi, clustereta));
+
+      // We set the multilinks flag of the ECAL element (for links to PS) to true. It will allow us to
+      // use it in an optimized way in prefilter
+      ecalElt->setIsValidMultilinks(true, _targetType);
     }
 
-    it->first->setMultilinks(multitracks);
+    // We set multilinks of the PS element (for links to ECAL)
+    psElt->setMultilinks(multitracks, _fieldType);
   }
 }
 

RecoParticleFlow/PFProducer/plugins/kdtrees/KDTreeLinkerTrackEcal.cc

@@ -1,4 +1,5 @@
 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
+#include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerBase.h"
 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
 
@@ -12,7 +13,7 @@ class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
   KDTreeLinkerTrackEcal(const edm::ParameterSet &conf);
   ~KDTreeLinkerTrackEcal() override;
 
-  // With this method, we create the list of psCluster that we want to link.
+  // With this method, we create the list of track that we want to link.
   void insertTargetElt(reco::PFBlockElement *track) override;
 
   // Here, we create the list of ecalCluster that we want to link. From ecalCluster
@@ -29,7 +30,7 @@ class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
   // all closest ecalClusters.
   void searchLinks() override;
 
-  // Here, we will store all PS/ECAL founded links in the PFBlockElement class
+  // Here, we will store all Track/ECAL founded links in the PFBlockElement class
   // of each psCluster in the PFmultilinks field.
   void updatePFBlockEltWithLinks() override;
 
@@ -45,7 +46,7 @@ class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
   RecHitSet rechitsSet_;
 
   // Map of linked Track/ECAL clusters.
-  BlockElt2BlockEltMap target2ClusterLinks_;
+  BlockElt2BlockEltMap cluster2TargetLinks_;
 
   // Map of the ECAL clusters associated to a rechit.
   RecHit2BlockEltMap rechit2ClusterLinks_;
@@ -141,10 +142,6 @@ void KDTreeLinkerTrackEcal::searchLinks() {
   for (BlockEltSet::iterator it = targetSet_.begin(); it != targetSet_.end(); it++) {
     reco::PFRecTrackRef trackref = (*it)->trackRefPF();
 
-    // We set the multilinks flag of the track to true. It will allow us to
-    // use in an optimized way our algo results in the recursive linking algo.
-    (*it)->setIsValidMultilinks(true);
-
     const reco::PFTrajectoryPoint &atECAL = trackref->extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax);
 
     // The track didn't reach ecal
@@ -204,7 +201,7 @@ void KDTreeLinkerTrackEcal::searchLinks() {
 
           // Check if the track and the cluster are linked
           if (deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
-            target2ClusterLinks_[*it].insert(*clusterIt);
+            cluster2TargetLinks_[*clusterIt].insert(*it);
 
         } else {  // ENDCAP
 
@@ -231,7 +228,7 @@ void KDTreeLinkerTrackEcal::searchLinks() {
 
           // Check if the track and the cluster are linked
           if (isinside)
-            target2ClusterLinks_[*it].insert(*clusterIt);
+            cluster2TargetLinks_[*clusterIt].insert(*it);
         }
       }
     }
@@ -241,18 +238,27 @@ void KDTreeLinkerTrackEcal::searchLinks() {
 void KDTreeLinkerTrackEcal::updatePFBlockEltWithLinks() {
   //TODO YG : Check if cluster positionREP() is valid ?
 
-  // Here we save in each track the list of phi/eta values of linked clusters.
-  for (BlockElt2BlockEltMap::iterator it = target2ClusterLinks_.begin(); it != target2ClusterLinks_.end(); ++it) {
+  // Here we save in each ECAL cluster the list of phi/eta values of linked tracks.
+  for (BlockElt2BlockEltMap::iterator it = cluster2TargetLinks_.begin(); it != cluster2TargetLinks_.end(); ++it) {
+    const auto &ecalElt = it->first;
+    const auto &trackEltSet = it->second;
     reco::PFMultiLinksTC multitracks(true);
 
-    for (BlockEltSet::iterator jt = it->second.begin(); jt != it->second.end(); ++jt) {
-      double clusterphi = (*jt)->clusterRef()->positionREP().phi();
-      double clustereta = (*jt)->clusterRef()->positionREP().eta();
+    for (const auto &trackElt : trackEltSet) {
+      const reco::PFRecTrackRef &trackref = trackElt->trackRefPF();
+      const reco::PFTrajectoryPoint &atECAL = trackref->extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax);
+      double tracketa = atECAL.positionREP().eta();
+      double trackphi = atECAL.positionREP().phi();
+
+      multitracks.linkedClusters.push_back(std::make_pair(trackphi, tracketa));
 
-      multitracks.linkedClusters.push_back(std::make_pair(clusterphi, clustereta));
+      // We set the multilinks flag of the track (for links to ECAL) to true. It will allow us to
+      // use it in an optimized way in prefilter
+      trackElt->setIsValidMultilinks(true, _fieldType);
     }
 
-    it->first->setMultilinks(multitracks);
+    // We set multilinks of the ECAL element (for links to tracks)
+    ecalElt->setMultilinks(multitracks, _targetType);
   }
 }
 
@@ -263,7 +269,7 @@ void KDTreeLinkerTrackEcal::clear() {
   rechitsSet_.clear();
 
   rechit2ClusterLinks_.clear();
-  target2ClusterLinks_.clear();
+  cluster2TargetLinks_.clear();
 
   tree_.clear();
 }