Authors: Caleb Fangmeier
Date: August 16, 2018

[TOC]

# Fixing hit-skipping in ElectronNHitSeedProducer

When the `ElectronNHitSeedProducer` was written for HLT, it had a limitation wherein hit-skipping was not allowed. What does this mean? The pixel-matching procedure has two sets of inputs.

  1. A collection of ECAL superclusters
  2. A collection of pixel seeds which are pairs/triplets/quadruplets of individual
  pixel hits on different layers. The individual sets of hits may or may not be disjoint
  depending on the configuration of the seed creation steps. Nominally, however, they *are* disjoint.

A supercluster matches with a pixel seed if the pixel-match procedure matches the supercluster with at least N hits in the seed, where N is

  - 3 if the track passes through 4+ active layers of pixels
  - 2 otherwise.

  > Note: This isn't strictly true since the old pixel match routine only ever requires 2 hits.

The old (current Offline) matching scheme (using `ElectronSeedProducer`), allows for skipping hits in the pixel seed to get to this minimal number of matches. For example, if the seed contains four barrel hits denominated by B1, B2, and B3, the matching procedure could match with B1 and B3 (skipping B2) to get the required two matched hits.

The new pixel matching scheme lacked the ability to skip B2 so it would match only B1, falling short of the required two hits. The hack to avoid this issue is to create pixel seeds from all pairs and triplets of hits. So, in the above example, there would be, excepting a few edge cased, four seeds created:

  1. B1, B2
  2. B2, B3
  3. B1, B3
  4. B1, B2, B3

Remember, B2 is the hit that fails to match with the super-cluster, but the supercluster would still end up matching with seed 3 so now this track will not be missed.

There are a couple issues with this hack, however.

  - There are an awful lot more seeds that get created since all combinations are tried.
  - The hits in the seeds are not disjoint so cross-cleaning should be employed downstream to avoid double counting.

A fix to avoid these issues is to just enable hit-skipping in the new pixel matching code. Appendices [A](#a-hit-skipping-patch) and [B](#b-hit-skipping-code) has a diff and the straight code to implement this. There is also a complementary change to the configuration to get rid of the "make all pairs/triplets" step and instead hook up the `ElectronNHitSeedProducer` with the standard set of pixel seeds from general tracking.

# Performance Studies

I'm going to first compare the effect of the hit skipping on the new seeding. Then I'll compare performance with the old seeding.

<div class="row">
<div class="col-*-6 row_fig">
fig!::
</div>
<div class="col-*-6 row_fig">
</div>
</div>

fig!::tracking_roc_curve_dR|Tracking Efficiency with $\Delta R$ truth-matching


# Appendices

## A: Hit-Skipping Patch

Here is the patch for `/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc` to enable hit-skipping.

```patch
diff --git a/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc b/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc
index 506895a..75dad73 100644
--- a/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc
+++ b/RecoEgamma/EgammaElectronAlgos/src/TrajSeedMatcher.cc
@@ -26,6 +26,7 @@ TrajSeedMatcher::TrajSeedMatcher(const edm::ParameterSet& pset):
   minNrHitsValidLayerBins_(pset.getParameter<std::vector<int> >("minNrHitsValidLayerBins"))
 {
   useRecoVertex_ = pset.getParameter<bool>("useRecoVertex");
+  allowHitSkipping_ = pset.getParameter<bool>("allowHitSkipping");
   navSchoolLabel_ = pset.getParameter<std::string>("navSchool");
   detLayerGeomLabel_ = pset.getParameter<std::string>("detLayerGeom");
   const auto cutsPSets=pset.getParameter<std::vector<edm::ParameterSet> >("matchingCuts");
@@ -52,6 +53,7 @@ edm::ParameterSetDescription TrajSeedMatcher::makePSetDescription()
 {
   edm::ParameterSetDescription desc;
   desc.add<bool>("useRecoVertex",false);
+  desc.add<bool>("allowHitSkipping", false);
   desc.add<std::string>("navSchool","SimpleNavigationSchool");
   desc.add<std::string>("detLayerGeom","hltESPGlobalDetLayerGeometry");
   desc.add<std::vector<int> >("minNrHitsValidLayerBins",{4});
@@ -113,6 +115,7 @@ TrajSeedMatcher::compatibleSeeds(const TrajectorySeedCollection& seeds, const Gl
   for(const auto& seed : seeds) {
     std::vector<SCHitMatch> matchedHitsNeg = processSeed(seed,candPos,vprim,energy,-1);
     std::vector<SCHitMatch> matchedHitsPos = processSeed(seed,candPos,vprim,energy,+1);
+
     int nrValidLayersPos = 0;
     int nrValidLayersNeg = 0;
     if(matchedHitsNeg.size()>=2){
@@ -128,28 +131,15 @@ TrajSeedMatcher::compatibleSeeds(const TrajectorySeedCollection& seeds, const Gl
 
     int nrValidLayers = std::max(nrValidLayersNeg,nrValidLayersPos);
     size_t nrHitsRequired = getNrHitsRequired(nrValidLayers);
-    //so we require the number of hits to exactly match, this is because we currently do not
-    //do any duplicate cleaning for the input seeds
-    //this means is a hit pair with a 3rd hit will appear twice (or however many hits it has)
-    //so if you did >=nrHitsRequired, you would get the same seed multiple times
-    //ideally we should fix this and clean our input seed collection so each seed is only in once
-    //also it should be studied what impact having a 3rd hit has on a GsfTrack
-    //ie will we get a significantly different result seeding with a hit pair 
-    //and the same the hit pair with a 3rd hit added
-    //in that case, perhaps it should be >=
-    if(matchedHitsNeg.size()==nrHitsRequired ||
-       matchedHitsPos.size()==nrHitsRequired){
+
+    if(matchedHitsNeg.size()>=nrHitsRequired ||
+       matchedHitsPos.size()>=nrHitsRequired){
       matchedSeeds.push_back({seed,matchedHitsPos,matchedHitsNeg,nrValidLayers});
     }
-    
-
   }
   return matchedSeeds;
 }
 
-//checks if the hits of the seed match within requested selection
-//matched hits are required to be consecutive, as soon as hit isnt matched,
-//the function returns, it doesnt allow skipping hits
 std::vector<TrajSeedMatcher::SCHitMatch>
 TrajSeedMatcher::processSeed(const TrajectorySeed& seed, const GlobalPoint& candPos,
                              const GlobalPoint & vprim, const float energy, const int charge)
@@ -159,32 +149,36 @@ TrajSeedMatcher::processSeed(const TrajectorySeed& seed, const GlobalPoint& cand
 
   FreeTrajectoryState trajStateFromVtx = FTSFromVertexToPointFactory::get(*magField_, candPos, vprim, energy, charge);
   PerpendicularBoundPlaneBuilder bpb;
-  TrajectoryStateOnSurface initialTrajState(trajStateFromVtx,*bpb(trajStateFromVtx.position(), 
-								  trajStateFromVtx.momentum()));
- 
+  TrajectoryStateOnSurface initialTrajState(trajStateFromVtx, *bpb(trajStateFromVtx.position(), trajStateFromVtx.momentum()));
   std::vector<SCHitMatch> matchedHits;
+  FreeTrajectoryState firstHitFreeTraj;
+  GlobalPoint prevHitPos;
+  GlobalPoint vertex;
+  for (size_t hitNr=0;hitNr<matchingCuts_.size() && hitNr<seed.nHits();hitNr++) {
+    bool matched = false;
+    if (matchedHits.size() == 0) {  // First hit is special
       SCHitMatch firstHit = matchFirstHit(seed, initialTrajState, vprim, *backwardPropagator_);
       firstHit.setExtra(candEt, candEta, candPos.phi(), charge, 1);
       if (passesMatchSel(firstHit, 0)) {
+        matched = true;
         matchedHits.push_back(firstHit);
-
         //now we can figure out the z vertex
         double zVertex = useRecoVertex_ ? vprim.z() : getZVtxFromExtrapolation(vprim,firstHit.hitPos(),candPos);
-    GlobalPoint vertex(vprim.x(),vprim.y(),zVertex);
-    
-    FreeTrajectoryState firstHitFreeTraj = FTSFromVertexToPointFactory::get(*magField_, firstHit.hitPos(), 
-									    vertex, energy, charge) ;
- 
-    GlobalPoint prevHitPos = firstHit.hitPos();
-    for(size_t hitNr=1;hitNr<matchingCuts_.size() && hitNr<seed.nHits();hitNr++){
+        vertex = GlobalPoint(vprim.x(), vprim.y(), zVertex);
+        firstHitFreeTraj = FTSFromVertexToPointFactory::get(*magField_, firstHit.hitPos(), vertex, energy, charge);
+        prevHitPos = firstHit.hitPos();
+      }
+    } else {  // All subsequent hits are handled the same
       SCHitMatch hit = match2ndToNthHit(seed, firstHitFreeTraj, hitNr, prevHitPos, vertex, *forwardPropagator_);
       hit.setExtra(candEt,candEta,candPos.phi(),charge,1);
-      if(passesMatchSel(hit,hitNr)){
+      if (passesMatchSel(hit, matchedHits.size())) {
+        matched = true;
         matchedHits.push_back(hit);
         prevHitPos = hit.hitPos();
-      }else break;
       }
     }
+    if (!matched and !allowHitSkipping_) break;
+  }
   return matchedHits;
 }
 
```

## B: Hit-Skipping Code

And Here is the relevant code on its own.

```c++
std::vector<TrajSeedMatcher::SeedWithInfo>
TrajSeedMatcher::compatibleSeeds(const TrajectorySeedCollection& seeds, const GlobalPoint& candPos,
                                 const GlobalPoint & vprim, const float energy)
{
  if(!forwardPropagator_ || !backwardPropagator_ || !magField_.isValid()){
    throw cms::Exception("LogicError") <<__FUNCTION__<<" can not make pixel seeds as event setup has not properly been called";
  }

  clearCache();

  std::vector<SeedWithInfo> matchedSeeds;
  for(const auto& seed : seeds) {
    std::vector<SCHitMatch> matchedHitsNeg = processSeed(seed,candPos,vprim,energy,-1);
    std::vector<SCHitMatch> matchedHitsPos = processSeed(seed,candPos,vprim,energy,+1);

    int nrValidLayersPos = 0;
    int nrValidLayersNeg = 0;
    if(matchedHitsNeg.size()>=2){
      nrValidLayersNeg = getNrValidLayersAlongTraj(matchedHitsNeg[0],
                                                   matchedHitsNeg[1],
                                                   candPos,vprim,energy,-1);
    }
    if(matchedHitsPos.size()>=2){
      nrValidLayersPos = getNrValidLayersAlongTraj(matchedHitsPos[0],
                                                   matchedHitsPos[1],
                                                   candPos,vprim,energy,+1);
    }

    int nrValidLayers = std::max(nrValidLayersNeg,nrValidLayersPos);
    size_t nrHitsRequired = getNrHitsRequired(nrValidLayers);

    if(matchedHitsNeg.size()>=nrHitsRequired ||
       matchedHitsPos.size()>=nrHitsRequired){
      matchedSeeds.push_back({seed,matchedHitsPos,matchedHitsNeg,nrValidLayers});
    }
  }
  return matchedSeeds;
}

std::vector<TrajSeedMatcher::SCHitMatch>
TrajSeedMatcher::processSeed(const TrajectorySeed& seed, const GlobalPoint& candPos,
                             const GlobalPoint & vprim, const float energy, const int charge)
{
  const float candEta = candPos.eta();
  const float candEt = energy*std::sin(candPos.theta());

  FreeTrajectoryState trajStateFromVtx = FTSFromVertexToPointFactory::get(*magField_, candPos, vprim, energy, charge);
  PerpendicularBoundPlaneBuilder bpb;
  TrajectoryStateOnSurface initialTrajState(trajStateFromVtx, *bpb(trajStateFromVtx.position(), trajStateFromVtx.momentum()));
  std::vector<SCHitMatch> matchedHits;
  FreeTrajectoryState firstHitFreeTraj;
  GlobalPoint prevHitPos;
  GlobalPoint vertex;
  for (size_t hitNr=0;hitNr<matchingCuts_.size() && hitNr<seed.nHits();hitNr++) {
    bool matched = false;
    if (matchedHits.size() == 0) {  // First hit is special
      SCHitMatch firstHit = matchFirstHit(seed, initialTrajState, vprim, *backwardPropagator_);
      firstHit.setExtra(candEt, candEta, candPos.phi(), charge, 1);
      if (passesMatchSel(firstHit, 0)) {
        matched = true;
        matchedHits.push_back(firstHit);
        //now we can figure out the z vertex
        double zVertex = useRecoVertex_ ? vprim.z() : getZVtxFromExtrapolation(vprim,firstHit.hitPos(),candPos);
        vertex = GlobalPoint(vprim.x(), vprim.y(), zVertex);
        firstHitFreeTraj = FTSFromVertexToPointFactory::get(*magField_, firstHit.hitPos(), vertex, energy, charge);
        prevHitPos = firstHit.hitPos();
      }
    } else {  // All subsequent hits are handled the same
      SCHitMatch hit = match2ndToNthHit(seed, firstHitFreeTraj, hitNr, prevHitPos, vertex, *forwardPropagator_);
      hit.setExtra(candEt,candEta,candPos.phi(),charge,1);
      if (passesMatchSel(hit, matchedHits.size())) {
        matched = true;
        matchedHits.push_back(hit);
        prevHitPos = hit.hitPos();
      }
    }
    if (!matched and !allowHitSkipping_) break;
  }
  return matchedHits;
}
```