From e0007acf4dd705015dac2b84eb9a661721376656 Mon Sep 17 00:00:00 2001
From: Jack Li <43122048+Jingyan95@users.noreply.github.com>
Date: Fri, 15 Oct 2021 18:27:12 +0200
Subject: [PATCH] L1 tk dev 12 0 0 pre4 hph (#94)

* First Commit

* Code Format

* Fix Typo

* Naming Rule

* Chris' 2nd Comment

Co-authored-by: Jack Li <jingyan.li@cern.ch>
---
 L1Trigger/TrackFindingTMTT/src/KFbase.cc      |   6 +-
 .../plugins/L1FPGATrackProducer.cc            |  10 +-
 .../python/L1HybridEmulationTracks_cff.py     |   2 +
 .../test/L1TrackNtupleMaker_cfg.py            |   1 +
 .../TrackTrigger/interface/HitPatternHelper.h | 169 ++++++++++
 .../interface/HitPatternHelperRcd.h           |  24 ++
 .../TrackTrigger/interface/TrackQuality.h     |   6 +
 L1Trigger/TrackTrigger/plugins/ProducerHPH.cc |  47 +++
 .../TrackTrigger/python/ProducerHPH_cff.py    |   5 +
 .../TrackTrigger/python/ProducerHPH_cfi.py    |  10 +
 .../TrackTrigger/src/ES_HitPatternHelper.cc   |   8 +
 .../TrackTrigger/src/HitPatternHelper.cc      | 294 ++++++++++++++++++
 .../TrackTrigger/src/HitPatternHelperRcd.cc   |   8 +
 L1Trigger/TrackTrigger/src/TrackQuality.cc    |  63 ++--
 14 files changed, 627 insertions(+), 26 deletions(-)
 create mode 100644 L1Trigger/TrackTrigger/interface/HitPatternHelper.h
 create mode 100644 L1Trigger/TrackTrigger/interface/HitPatternHelperRcd.h
 create mode 100644 L1Trigger/TrackTrigger/plugins/ProducerHPH.cc
 create mode 100644 L1Trigger/TrackTrigger/python/ProducerHPH_cff.py
 create mode 100644 L1Trigger/TrackTrigger/python/ProducerHPH_cfi.py
 create mode 100644 L1Trigger/TrackTrigger/src/ES_HitPatternHelper.cc
 create mode 100644 L1Trigger/TrackTrigger/src/HitPatternHelper.cc
 create mode 100644 L1Trigger/TrackTrigger/src/HitPatternHelperRcd.cc

diff --git a/L1Trigger/TrackFindingTMTT/src/KFbase.cc b/L1Trigger/TrackFindingTMTT/src/KFbase.cc
index 7096e8c088272..4de598fac34a6 100644
--- a/L1Trigger/TrackFindingTMTT/src/KFbase.cc
+++ b/L1Trigger/TrackFindingTMTT/src/KFbase.cc
@@ -729,9 +729,9 @@ namespace tmtt {
     // Fixes to layermap when "maybe layer" used
     if (settings_->kfUseMaybeLayers()) {
       switch (kfEtaReg) {
-        case 5:  //case 5: B1 B2 (B3+B4)* D1 D2 D3+D4 D5+D6  -- B3 is combined with B4 and is flagged as "maybe layer"
+        case 5:  //case 5: B1 B2 (B3+B4)* D1 D2 D3 D4+D5  -- B3 is combined with B4 and is flagged as "maybe layer"
           if (layerIDreduced == 6) {
-            kalmanLay = 5;
+            kalmanLay = 6;
           }
           break;
         case 6:  //case 6: B1* B2* D1 D2 D3 D4 D5 -- B1 and B2 are flagged as "maybe layer"
@@ -758,7 +758,7 @@ namespace tmtt {
           }
           break;
           //case 5:  // B1 B2 B3+B4 D1 D2 D3 D4/D5
-        case 5:  // B1 B2 B3 D1+B4 D2 D3 D4/D5
+        case 5:
           if (layerIDreduced == 5) {
             kalmanLay = 5;
           } else if (layerIDreduced == 7) {
diff --git a/L1Trigger/TrackFindingTracklet/plugins/L1FPGATrackProducer.cc b/L1Trigger/TrackFindingTracklet/plugins/L1FPGATrackProducer.cc
index 1c9c248b2ad12..82bf85c13bdad 100644
--- a/L1Trigger/TrackFindingTracklet/plugins/L1FPGATrackProducer.cc
+++ b/L1Trigger/TrackFindingTracklet/plugins/L1FPGATrackProducer.cc
@@ -94,6 +94,7 @@
 
 #include "L1Trigger/TrackTrigger/interface/StubPtConsistency.h"
 #include "L1Trigger/TrackTrigger/interface/TrackQuality.h"
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelper.h"
 
 //////////////
 // STD HEADERS
@@ -182,9 +183,11 @@ class L1FPGATrackProducer : public edm::one::EDProducer<edm::one::WatchRuns> {
 
   // helper class to store DTC configuration
   trackerDTC::Setup setup_;
+  const hph::Setup* setupHPH_;
 
   // Setup token
-  edm::ESGetToken<trackerDTC::Setup, trackerDTC::SetupRcd> esGetToken_;
+  const edm::ESGetToken<trackerDTC::Setup, trackerDTC::SetupRcd> esGetToken_;
+  const edm::ESGetToken<hph::Setup, hph::SetupRcd> esGetTokenHPH_;
   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
 
   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tGeomToken_;
@@ -237,6 +240,7 @@ L1FPGATrackProducer::L1FPGATrackProducer(edm::ParameterSet const& iConfig)
 
   // book ES product
   esGetToken_ = esConsumes<trackerDTC::Setup, trackerDTC::SetupRcd, edm::Transition::BeginRun>();
+  esGetTokenHPH_ = esConsumes<hph::Setup, hph::SetupRcd, edm::Transition::BeginRun>();
 
   // --------------------------------------------------------------------------------
   // set options in Settings based on inputs from configuration files
@@ -304,6 +308,10 @@ void L1FPGATrackProducer::beginRun(const edm::Run& run, const edm::EventSetup& i
   settings.setBfield(mMagneticFieldStrength);
 
   setup_ = iSetup.getData(esGetToken_);
+  setupHPH_ = &iSetup.getData(esGetTokenHPH_);
+  if (trackQuality_) {
+    trackQualityModel_->setHPHSetup(setupHPH_);
+  }
 
   // initialize the tracklet event processing (this sets all the processing & memory modules, wiring, etc)
   eventProcessor.init(settings);
diff --git a/L1Trigger/TrackFindingTracklet/python/L1HybridEmulationTracks_cff.py b/L1Trigger/TrackFindingTracklet/python/L1HybridEmulationTracks_cff.py
index e977bd04faa73..dcaea4d59c4ed 100644
--- a/L1Trigger/TrackFindingTracklet/python/L1HybridEmulationTracks_cff.py
+++ b/L1Trigger/TrackFindingTracklet/python/L1HybridEmulationTracks_cff.py
@@ -6,6 +6,8 @@
 
 from SimTracker.TrackTriggerAssociation.TrackTriggerAssociator_cff import *
 
+from L1Trigger.TrackTrigger.ProducerHPH_cff import *
+
 # prompt hybrid emulation
 TTTrackAssociatorFromPixelDigis.TTTracks = cms.VInputTag(cms.InputTag("TTTracksFromTrackletEmulation", "Level1TTTracks") )
 
diff --git a/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py b/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py
index 333cd14103252..a0b767205942d 100644
--- a/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py
+++ b/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py
@@ -30,6 +30,7 @@
 process.load('FWCore.MessageService.MessageLogger_cfi')
 process.MessageLogger.L1track = dict(limit = -1)
 process.MessageLogger.Tracklet = dict(limit = -1)
+process.MessageLogger.TrackTriggerHPH = dict(limit = -1)
 
 if GEOMETRY == "D49": 
     print("using geometry " + GEOMETRY + " (tilted)")
diff --git a/L1Trigger/TrackTrigger/interface/HitPatternHelper.h b/L1Trigger/TrackTrigger/interface/HitPatternHelper.h
new file mode 100644
index 0000000000000..b274166f37746
--- /dev/null
+++ b/L1Trigger/TrackTrigger/interface/HitPatternHelper.h
@@ -0,0 +1,169 @@
+// This is a helper function that can be used to decode hitpattern, which is a 7-bit integer produced by KF.
+//
+// There are three classes declared in HitPatternHelper (HPH) namesapce:
+// 1)SensorModule: This is used to store important information about the sensor modules. For example r,z coordinates.
+// 2)Setup: This is used to produce a collection of <SensorModule> needed by HPH.
+// 3)HitPatternHelper: This is used to decode hitpattern with the help of the information from sensor modules and layermap.
+//
+// Two predictions on which layers particles will hit are made using different information:
+// i)Loop over sensor modules and make predictions based on spatial coordinates of tracks. This prediction is considered more accurate.
+// ii)Make predictions based on a hard-coded layermap. This prediction is considered less accurate and is used by Old KF to encode hitpattern.
+//
+//
+//  Created by J.Li on 1/23/21.
+//
+
+#ifndef L1Trigger_TrackTrigger_interface_HitPatternHelper_h
+#define L1Trigger_TrackTrigger_interface_HitPatternHelper_h
+
+#include "FWCore/Framework/interface/data_default_record_trait.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
+#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
+#include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
+#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelperRcd.h"
+
+#include <bitset>
+#include <iostream>
+#include <vector>
+
+using namespace std;
+using namespace edm;
+
+namespace hph {
+
+  class SensorModule {
+  public:
+    SensorModule() {}
+    SensorModule(
+        bool isBarrel, bool isPS, int numColumns, int layerId, double r, double z, double pitchCol, double tilt);
+    ~SensorModule() {}
+
+    bool isBarrel() const { return isBarrel_; }
+    bool isPS() const { return isPS_; }
+    bool isMaybe() const { return isMaybe_; }
+    int numColumns() const { return numColumns_; }
+    int layerId() const { return layerId_; }
+    double r() const { return r_; }
+    double z() const { return z_; }
+    double pitchCol() const { return pitchCol_; }
+    double tilt() const { return tilt_; }
+    double sin() const { return sin_; }
+    double cos() const { return cos_; }
+
+    void setMaybe() { isMaybe_ = true; }
+
+  private:
+    bool isBarrel_;
+    bool isPS_;
+    bool isMaybe_;
+    int numColumns_;
+    int layerId_;
+    double r_;
+    double z_;
+    double pitchCol_;
+    double tilt_;
+    double sin_;
+    double cos_;
+  };
+
+  class Setup {
+  public:
+    Setup() {}
+    Setup(const edm::ParameterSet& iConfig,
+          const TrackerGeometry& trackerGeometry,
+          const TrackerTopology& trackerTopology);
+    ~Setup() {}
+
+    static auto smallerR(SensorModule lhs, SensorModule rhs) { return lhs.r() < rhs.r(); }
+    static auto smallerZ(SensorModule lhs, SensorModule rhs) { return lhs.z() < rhs.z(); }
+    static auto equalRZ(SensorModule lhs, SensorModule rhs) {
+      return abs(lhs.r() - rhs.r()) < delta_ && abs(lhs.z() - rhs.z()) < delta_;
+    }
+    std::vector<SensorModule> sensorModules() const { return sensorModules_; }
+
+    bool hphDebug() const { return iConfig_.getParameter<bool>("hphDebug"); }
+    bool useNewKF() const { return iConfig_.getParameter<bool>("useNewKF"); }
+    double chosenRofZ() const { return iConfig_.getParameter<double>("chosenRofZ"); }
+    double deltaTanL() const { return iConfig_.getParameter<double>("deltaTanL"); }
+
+  private:
+    edm::ParameterSet iConfig_;
+    const TrackerGeometry* trackerGeometry_;
+    const TrackerTopology* trackerTopology_;
+    static constexpr double delta_ = 1.e-3;
+    std::vector<SensorModule> sensorModules_;
+  };
+
+  class HitPatternHelper {
+  public:
+    HitPatternHelper() {}
+    HitPatternHelper(const Setup* setup, int hitpattern, double cot, double z0);
+    ~HitPatternHelper() {}
+
+    int etaSector() { return etaSector_; }        //Eta sectors defined in KF
+    int numExpLayer() { return numExpLayer_; }    //The number of layers KF expects
+    int numMissingPS() { return numMissingPS_; }  //The number of PS layers that are missing
+    int numMissing2S() { return numMissing2S_; }  //The number of 2S layers that are missing
+    int numPS() { return numPS_; }                //The number of PS layers are found in hitpattern
+    int num2S() { return num2S_; }                //The number of 2S layers are found in hitpattern
+    int numMissingInterior1() {
+      return numMissingInterior1_;
+    }  //The number of missing interior layers (using only hitpattern)
+    int numMissingInterior2() {
+      return numMissingInterior2_;
+    }  //The number of missing interior layers (using hitpattern and sensor modules)
+    std::vector<int> binary() { return binary_; }  //11-bit hitmask needed by TrackQuality.cc (0~5->L1~L6;6~10->D1~D5)
+    static auto smallerID(SensorModule lhs, SensorModule rhs) { return lhs.layerId() < rhs.layerId(); }
+    static auto equalID(SensorModule lhs, SensorModule rhs) { return lhs.layerId() == rhs.layerId(); }
+
+    int ReducedId(
+        int layerId);  //Converts layer id (1~6->L1~L6;11~15->D1~D5) to reduced layer id (0~5->L1~L6;6~10->D1~D5)
+    int findLayer(int layerId);  //Search for a layer id from sensor modules
+
+  private:
+    int etaSector_;
+    int hitpattern_;
+    int numExpLayer_;
+    int numMissingLayer_;
+    int numMissingPS_;
+    int numMissing2S_;
+    int numPS_;
+    int num2S_;
+    int numMissingInterior1_;
+    int numMissingInterior2_;
+    double cot_;
+    double z0_;
+    const Setup* setup_;
+    std::vector<SensorModule> layers_;  //Sensor modules that particles are expected to hit
+    std::vector<int> binary_;
+    bool hphDebug_;
+    bool useNewKF_;
+    float chosenRofZ_;
+    float deltaTanL_;
+    std::vector<float> etaRegions_ = {
+        -2.4, -2.08, -1.68, -1.26, -0.90, -0.62, -0.41, -0.20, 0.0, 0.20, 0.41, 0.62, 0.90, 1.26, 1.68, 2.08, 2.4};
+
+    //Layermap used in Old KF
+    //Ultimate config is assumed (with maybe layer)
+    //Index across is kalman layer
+    //Index down is eta sector
+    //Element is layer id where barrel layers=1,2,3,4,5,6 & endcap wheels=11,12,13,14,15; 0 is invalid.
+    std::vector<int> hitmap_[8][7] = {
+        {{1}, {2}, {3}, {4}, {5}, {6}, {0}},
+        {{1}, {2}, {3}, {4}, {5}, {6}, {0}},
+        {{1}, {2}, {3}, {4}, {5}, {6}, {0}},
+        {{1}, {2}, {3}, {4}, {5}, {6}, {0}},
+        {{1}, {2}, {3}, {4}, {5, 11}, {6, 12}, {13}},
+        {{1}, {2}, {3, 4}, {11}, {12}, {13}, {14, 15}},
+        {{1}, {2}, {11}, {12}, {13}, {14}, {15}},
+        {{1}, {11}, {12}, {13}, {14}, {15}, {0}},
+    };
+  };
+
+}  // namespace hph
+
+EVENTSETUP_DATA_DEFAULT_RECORD(hph::Setup, hph::SetupRcd);
+
+#endif
diff --git a/L1Trigger/TrackTrigger/interface/HitPatternHelperRcd.h b/L1Trigger/TrackTrigger/interface/HitPatternHelperRcd.h
new file mode 100644
index 0000000000000..76938ff2f39c8
--- /dev/null
+++ b/L1Trigger/TrackTrigger/interface/HitPatternHelperRcd.h
@@ -0,0 +1,24 @@
+//
+//  Created by J.Li on 1/23/21.
+//
+
+#ifndef L1Trigger_TrackTrigger_interface_HitPatternHelperRcd_h
+#define L1Trigger_TrackTrigger_interface_HitPatternHelperRcd_h
+
+#include "FWCore/Framework/interface/DependentRecordImplementation.h"
+
+#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
+#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+
+#include "FWCore/Utilities/interface/mplVector.h"
+
+namespace hph {
+
+  typedef edm::mpl::Vector<TrackerDigiGeometryRecord, TrackerTopologyRcd> Rcds;
+
+  // record of hph::SetupRcd
+  class SetupRcd : public edm::eventsetup::DependentRecordImplementation<SetupRcd, Rcds> {};
+
+}  // namespace hph
+
+#endif
diff --git a/L1Trigger/TrackTrigger/interface/TrackQuality.h b/L1Trigger/TrackTrigger/interface/TrackQuality.h
index 4e241c607acb5..a8558ce8d4984 100644
--- a/L1Trigger/TrackTrigger/interface/TrackQuality.h
+++ b/L1Trigger/TrackTrigger/interface/TrackQuality.h
@@ -26,6 +26,8 @@ C.Brown 28/07/20
 #include "DataFormats/L1TrackTrigger/interface/TTTrack.h"
 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
 
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelper.h"
+
 class TrackQuality {
 public:
   // Enum class used for determining prediction behaviour in setTrackQuality
@@ -60,6 +62,8 @@ class TrackQuality {
                     std::string const& ONNXInputName,
                     std::vector<std::string> const& featureNames);
 
+  void setHPHSetup(const hph::Setup* setup);
+
 private:
   // Private Member Data
   QualityAlgorithm qualityAlgorithm_ = QualityAlgorithm::None;
@@ -73,5 +77,7 @@ class TrackQuality {
   float minPt_;
   int nStubsmin_;
   float ONNXInvRScaling_;
+  const hph::Setup* setup_;
+  bool useHPH;
 };
 #endif
diff --git a/L1Trigger/TrackTrigger/plugins/ProducerHPH.cc b/L1Trigger/TrackTrigger/plugins/ProducerHPH.cc
new file mode 100644
index 0000000000000..56f8ecd4a8911
--- /dev/null
+++ b/L1Trigger/TrackTrigger/plugins/ProducerHPH.cc
@@ -0,0 +1,47 @@
+//
+//  Created by J.Li on 1/23/21.
+//
+
+#include "FWCore/Framework/interface/ModuleFactory.h"
+#include "FWCore/Framework/interface/ESProducer.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/ESGetToken.h"
+#include "FWCore/Utilities/interface/ESInputTag.h"
+#include "DataFormats/Provenance/interface/ParameterSetID.h"
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelper.h"
+
+#include <memory>
+
+using namespace std;
+using namespace edm;
+
+namespace hph {
+
+  class ProducerHPH : public ESProducer {
+  public:
+    ProducerHPH(const ParameterSet& iConfig);
+    ~ProducerHPH() override {}
+    unique_ptr<Setup> produce(const SetupRcd& Rcd);
+
+  private:
+    const ParameterSet iConfig_;
+    ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> getTokenTrackerGeometry_;
+    ESGetToken<TrackerTopology, TrackerTopologyRcd> getTokenTrackerTopology_;
+  };
+
+  ProducerHPH::ProducerHPH(const ParameterSet& iConfig) : iConfig_(iConfig) {
+    auto cc = setWhatProduced(this);
+    getTokenTrackerGeometry_ = cc.consumes();
+    getTokenTrackerTopology_ = cc.consumes();
+  }
+
+  unique_ptr<Setup> ProducerHPH::produce(const SetupRcd& Rcd) {
+    const TrackerGeometry& trackerGeometry = Rcd.get(getTokenTrackerGeometry_);
+    const TrackerTopology& trackerTopology = Rcd.get(getTokenTrackerTopology_);
+    return make_unique<Setup>(iConfig_, trackerGeometry, trackerTopology);
+  }
+
+}  // namespace hph
+
+DEFINE_FWK_EVENTSETUP_MODULE(hph::ProducerHPH);
diff --git a/L1Trigger/TrackTrigger/python/ProducerHPH_cff.py b/L1Trigger/TrackTrigger/python/ProducerHPH_cff.py
new file mode 100644
index 0000000000000..d788e09a4748b
--- /dev/null
+++ b/L1Trigger/TrackTrigger/python/ProducerHPH_cff.py
@@ -0,0 +1,5 @@
+import FWCore.ParameterSet.Config as cms
+
+from L1Trigger.TrackTrigger.ProducerHPH_cfi import HitPatternHelper_params
+
+HitPatternHelperSetup = cms.ESProducer("hph::ProducerHPH", HitPatternHelper_params)
diff --git a/L1Trigger/TrackTrigger/python/ProducerHPH_cfi.py b/L1Trigger/TrackTrigger/python/ProducerHPH_cfi.py
new file mode 100644
index 0000000000000..028fa18118e09
--- /dev/null
+++ b/L1Trigger/TrackTrigger/python/ProducerHPH_cfi.py
@@ -0,0 +1,10 @@
+import FWCore.ParameterSet.Config as cms
+
+HitPatternHelper_params = cms.PSet (
+
+  HPHdebug   = cms.bool(False),   
+  useNewKF   = cms.bool(False),
+  chosenRofZ = cms.double(50.0),
+  deltaTanL = cms.double(0.125)
+
+)
diff --git a/L1Trigger/TrackTrigger/src/ES_HitPatternHelper.cc b/L1Trigger/TrackTrigger/src/ES_HitPatternHelper.cc
new file mode 100644
index 0000000000000..24881410dc8f7
--- /dev/null
+++ b/L1Trigger/TrackTrigger/src/ES_HitPatternHelper.cc
@@ -0,0 +1,8 @@
+//
+//  Created by J.Li on 1/23/21.
+//
+
+#include "FWCore/Utilities/interface/typelookup.h"
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelper.h"
+
+TYPELOOKUP_DATA_REG(hph::Setup);
diff --git a/L1Trigger/TrackTrigger/src/HitPatternHelper.cc b/L1Trigger/TrackTrigger/src/HitPatternHelper.cc
new file mode 100644
index 0000000000000..b0991307a5ed0
--- /dev/null
+++ b/L1Trigger/TrackTrigger/src/HitPatternHelper.cc
@@ -0,0 +1,294 @@
+//
+//  Created by J.Li on 1/23/21.
+//
+
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelper.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include <algorithm>
+#include <cmath>
+
+namespace hph {
+
+  SensorModule::SensorModule(
+      bool isBarrel, bool isPS, int numColumns, int layerId, double r, double z, double pitchCol, double tilt)
+      : isBarrel_(isBarrel),
+        isPS_(isPS),
+        isMaybe_(false),
+        numColumns_(numColumns),
+        layerId_(layerId),
+        r_(r),
+        z_(z),
+        pitchCol_(pitchCol),
+        tilt_(tilt) {
+    sin_ = std::sin(tilt_);
+    cos_ = std::cos(tilt_);
+  }
+
+  Setup::Setup(const edm::ParameterSet& iConfig,
+               const TrackerGeometry& trackerGeometry,
+               const TrackerTopology& trackerTopology)
+      : trackerGeometry_(&trackerGeometry), trackerTopology_(&trackerTopology) {
+    iConfig_ = iConfig;
+
+    for (const auto& gd : trackerGeometry_->dets()) {
+      DetId detid = (*gd).geographicalId();
+      if (detid.subdetId() != StripSubdetector::TOB && detid.subdetId() != StripSubdetector::TID)
+        continue;  // only run on OT
+      if (!trackerTopology_->isLower(detid))
+        continue;  // loop on the stacks: choose the lower arbitrarily
+
+      // Get the DetSets of the Clusters
+      const GeomDetUnit* det0 = trackerGeometry_->idToDetUnit(detid);
+      const auto* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(det0);
+      const PixelTopology* topol = dynamic_cast<const PixelTopology*>(&(theGeomDet->specificTopology()));
+      const GlobalPoint pos0 = det0->position();
+      const GlobalPoint pos1 = trackerGeometry_->idToDetUnit(trackerTopology_->partnerDetId(detid))->position();
+      double r = pos0.perp();
+      double z = pos0.z();
+
+      bool flipped = pos0.mag() > pos1.mag();
+      bool isBarrel = detid.subdetId() == StripSubdetector::TOB;
+      bool isPS = trackerGeometry_->getDetectorType(detid) == TrackerGeometry::ModuleType::Ph2PSP;
+      double tilt = flipped ? atan2(pos1.z() - pos0.z(), pos0.perp() - pos1.perp())
+                            : atan2(pos0.z() - pos1.z(), pos1.perp() - pos0.perp());
+
+      int layerId = isBarrel ? trackerTopology_->layer(detid) : trackerTopology_->layer(detid) + 10;
+      int numColumns = topol->ncolumns();
+      double pitchCol = topol->pitch().second;
+
+      sensorModules_.emplace_back(isBarrel, isPS, numColumns, layerId, r, z, pitchCol, tilt);
+    }
+
+    sort(sensorModules_.begin(), sensorModules_.end(), smallerR);
+    sort(sensorModules_.begin(), sensorModules_.end(), smallerZ);
+    sensorModules_.erase(unique(sensorModules_.begin(), sensorModules_.end(), equalRZ), sensorModules_.end());
+  }
+
+  HitPatternHelper::HitPatternHelper(const Setup* setup, int hitpattern, double cot, double z0)
+      : hitpattern_(hitpattern),
+        numExpLayer_(0),
+        numMissingLayer_(0),
+        numMissingPS_(0),
+        numMissing2S_(0),
+        numPS_(0),
+        num2S_(0),
+        numMissingInterior1_(0),
+        numMissingInterior2_(0),
+        cot_(cot),
+        z0_(z0),
+        setup_(setup),
+        layers_(),
+        binary_(11, 0),
+        hphDebug_(setup_->hphDebug()),
+        useNewKF_(setup_->useNewKF()),
+        chosenRofZ_(setup_->chosenRofZ()),
+        deltaTanL_(setup_->deltaTanL()) {
+    //Calculating eta sector based on cot and z0
+    float kfzRef = z0_ + chosenRofZ_ * cot_;
+    int kf_eta_reg = 0;
+    for (int iEtaSec = 1; iEtaSec < ((int)etaRegions_.size() - 1); iEtaSec++) {  // Doesn't apply eta < 2.4 cut.
+      float etaMax = etaRegions_[iEtaSec];
+      float zRefMax = chosenRofZ_ / tan(2. * atan(exp(-etaMax)));
+      if (kfzRef > zRefMax) {
+        kf_eta_reg = iEtaSec;
+      }
+    }
+    etaSector_ = kf_eta_reg;
+    if (kf_eta_reg < ((int)etaRegions_.size() - 1) / 2) {
+      kf_eta_reg = ((int)etaRegions_.size() - 1) / 2 - 1 - kf_eta_reg;
+    } else {
+      kf_eta_reg = kf_eta_reg - (int)(etaRegions_.size() - 1) / 2;
+    }
+    //Looping over sensor modules to make predictions on which layers particles are expected to hit
+    for (SensorModule sm : setup_->sensorModules()) {
+      double d = (z0_ - sm.z() + sm.r() * cot_) / (sm.cos() - sm.sin() * cot_);
+      double d_p = (z0_ - sm.z() + sm.r() * (cot_ + deltaTanL_ / 2)) / (sm.cos() - sm.sin() * (cot_ + deltaTanL_ / 2));
+      double d_m = (z0_ - sm.z() + sm.r() * (cot_ - deltaTanL_ / 2)) / (sm.cos() - sm.sin() * (cot_ - deltaTanL_ / 2));
+      if (!(abs(d_p) < sm.numColumns() * sm.pitchCol() / 2. && abs(d_m) < sm.numColumns() * sm.pitchCol() / 2.))
+        sm.setMaybe();
+      if (useNewKF_ &&
+          (abs(d_p) < sm.numColumns() * sm.pitchCol() / 2. || abs(d_m) < sm.numColumns() * sm.pitchCol() / 2.)) {
+        layers_.push_back(sm);
+      }
+      if (!useNewKF_ && abs(d) < sm.numColumns() * sm.pitchCol() / 2.) {
+        layers_.push_back(sm);
+      }
+    }
+    //layers_ constains all the sensor modules that particles are expected to hit
+    sort(layers_.begin(), layers_.end(), smallerID);
+    layers_.erase(unique(layers_.begin(), layers_.end(), equalID), layers_.end());  //Keep only one sensor per layer
+
+    numExpLayer_ = layers_.size();
+
+    int nbits = floor(log2(hitpattern_)) + 1;
+    int lay_i = 0;
+    bool seq = false;
+    for (int i = 0; i < nbits; i++) {
+      lay_i = ((1 << i) & hitpattern_) >> i;  //0 or 1 in ith bit (right to left)
+
+      if (lay_i && !seq)
+        seq = true;  //sequence starts when first 1 found
+      if (!lay_i && seq) {
+        numMissingInterior1_++;  //This is the same as the "tmp_trk_nlaymiss_interior" calculated in Trackquality.cc
+      }
+      if (!lay_i) {
+        bool realhit = false;
+        for (int j : hitmap_[kf_eta_reg][i]) {
+          if (j < 1)
+            continue;
+          int k = findLayer(j);
+          if (k > 0)
+            realhit = true;
+        }
+        if (realhit)
+          numMissingInterior2_++;
+      }
+    }
+
+    if (hphDebug_) {
+      if (useNewKF_) {
+        edm::LogVerbatim("TrackTriggerHPH") << "Running with New KF";
+      } else {
+        edm::LogVerbatim("TrackTriggerHPH") << "Running with Old KF";
+      }
+      edm::LogVerbatim("TrackTriggerHPH") << "======================================================";
+      edm::LogVerbatim("TrackTriggerHPH")
+          << "Looking at hitpattern " << std::bitset<7>(hitpattern_) << "; Looping over KF layers:";
+    }
+
+    if (useNewKF_) {
+      //New KF uses sensor modules to determine the hitmask already
+      for (int i = 0; i < numExpLayer_; i++) {
+        if (hphDebug_) {
+          edm::LogVerbatim("TrackTriggerHPH") << "--------------------------";
+          edm::LogVerbatim("TrackTriggerHPH") << "Looking at KF layer " << i;
+          if (layers_[i].layerId() < 10) {
+            edm::LogVerbatim("TrackTriggerHPH") << "KF expects L" << layers_[i].layerId();
+          } else {
+            edm::LogVerbatim("TrackTriggerHPH") << "KF expects D" << layers_[i].layerId() - 10;
+          }
+        }
+
+        if (((1 << i) & hitpattern_) >> i) {
+          if (hphDebug_) {
+            edm::LogVerbatim("TrackTriggerHPH") << "Layer found in hitpattern";
+          }
+
+          binary_[ReducedId(layers_[i].layerId())] = 1;
+          if (layers_[i].isPS()) {
+            numPS_++;
+          } else {
+            num2S_++;
+          }
+        } else {
+          if (hphDebug_) {
+            edm::LogVerbatim("TrackTriggerHPH") << "Layer missing in hitpattern";
+          }
+
+          if (layers_[i].isPS()) {
+            numMissingPS_++;
+          } else {
+            numMissing2S_++;
+          }
+        }
+      }
+
+    } else {
+      //Old KF uses the hard coded layermap to determien hitmask
+      for (int i = 0; i < 7; i++) {  //Loop over each digit of hitpattern
+
+        if (hphDebug_) {
+          edm::LogVerbatim("TrackTriggerHPH") << "--------------------------";
+          edm::LogVerbatim("TrackTriggerHPH") << "Looking at KF layer " << i;
+        }
+
+        for (int j :
+             hitmap_[kf_eta_reg][i]) {  //Find out which layer the Old KF is dealing with when hitpattern is encoded
+          if (j < 1) {
+            if (hphDebug_) {
+              edm::LogVerbatim("TrackTriggerHPH") << "KF does not expect this layer";
+            }
+
+            continue;
+          }
+
+          if (hphDebug_) {
+            if (j < 10) {
+              edm::LogVerbatim("TrackTriggerHPH") << "KF expects L" << j;
+            } else {
+              edm::LogVerbatim("TrackTriggerHPH") << "KF expects D" << j - 10;
+            }
+          }
+
+          int k = findLayer(j);
+          if (k < 0) {
+            //k<0 means even though layer j is predicted by Old KF, this prediction is rejected because it contradicts
+            if (hphDebug_) {  //a more accurate prediction made with the help of information from sensor modules.
+              edm::LogVerbatim("TrackTriggerHPH") << "Rejected by sensor modules";
+            }
+
+            continue;
+          }
+
+          if (hphDebug_) {
+            edm::LogVerbatim("TrackTriggerHPH") << "Confirmed by sensor modules";
+          }
+          //prediction is accepted
+          if (((1 << i) & hitpattern_) >> i) {
+            if (hphDebug_) {
+              edm::LogVerbatim("TrackTriggerHPH") << "Layer found in hitpattern";
+            }
+
+            binary_[ReducedId(j)] = 1;
+            if (layers_[k].isPS()) {
+              numPS_++;
+            } else {
+              num2S_++;
+            }
+          } else {
+            if (hphDebug_) {
+              edm::LogVerbatim("TrackTriggerHPH") << "Layer missing in hitpattern";
+            }
+
+            if (layers_[k].isPS()) {
+              numMissingPS_++;
+            } else {
+              numMissing2S_++;
+            }
+          }
+        }
+      }
+    }
+
+    if (hphDebug_) {
+      edm::LogVerbatim("TrackTriggerHPH") << "------------------------------";
+      edm::LogVerbatim("TrackTriggerHPH") << "numPS = " << numPS_ << ", num2S = " << num2S_
+                                          << ", missingPS = " << numMissingPS_ << ", missing2S = " << numMissing2S_;
+      edm::LogVerbatim("TrackTriggerHPH") << "======================================================";
+    }
+  }
+
+  int HitPatternHelper::ReducedId(int layerId) {
+    if (hphDebug_ && (layerId > 15 || layerId < 1)) {
+      edm::LogVerbatim("TrackTriggerHPH") << "Warning: invalid layer id !";
+    }
+    if (layerId <= 6) {
+      layerId = layerId - 1;
+      return layerId;
+    } else {
+      layerId = layerId - 5;
+      return layerId;
+    }
+  };
+
+  int HitPatternHelper::findLayer(int layerId) {
+    for (int i = 0; i < (int)layers_.size(); i++) {
+      if (layerId == (int)layers_[i].layerId()) {
+        return i;
+      }
+    }
+    return -1;
+  }
+
+}  // namespace hph
diff --git a/L1Trigger/TrackTrigger/src/HitPatternHelperRcd.cc b/L1Trigger/TrackTrigger/src/HitPatternHelperRcd.cc
new file mode 100644
index 0000000000000..b502d531ccc05
--- /dev/null
+++ b/L1Trigger/TrackTrigger/src/HitPatternHelperRcd.cc
@@ -0,0 +1,8 @@
+//
+//  Created by J.Li on 1/23/21.
+//
+
+#include "L1Trigger/TrackTrigger/interface/HitPatternHelperRcd.h"
+#include "FWCore/Framework/interface/eventsetuprecord_registration_macro.h"
+
+EVENTSETUP_RECORD_REG(hph::SetupRcd);
diff --git a/L1Trigger/TrackTrigger/src/TrackQuality.cc b/L1Trigger/TrackTrigger/src/TrackQuality.cc
index 2f2615ce1608b..f32cc65df2ee6 100644
--- a/L1Trigger/TrackTrigger/src/TrackQuality.cc
+++ b/L1Trigger/TrackTrigger/src/TrackQuality.cc
@@ -11,7 +11,7 @@ C.Brown & C.Savard 07/2020
 
 TrackQuality::TrackQuality() {}
 
-TrackQuality::TrackQuality(const edm::ParameterSet& qualityParams) {
+TrackQuality::TrackQuality(const edm::ParameterSet& qualityParams) : setup_(), useHPH(false) {
   std::string AlgorithmString = qualityParams.getParameter<std::string>("qualityAlgorithm");
   // Unpacks EDM parameter set itself to save unecessary processing within TrackProducers
   if (AlgorithmString == "Cut") {
@@ -39,7 +39,8 @@ std::vector<float> TrackQuality::featureTransform(TTTrack<Ref_Phase2TrackerDigi_
   // {"log_chi2","log_chi2rphi","log_chi2rz","log_bendchi2","nstubs","lay1_hits","lay2_hits",
   // "lay3_hits","lay4_hits","lay5_hits","lay6_hits","disk1_hits","disk2_hits","disk3_hits",
   // "disk4_hits","disk5_hits","rinv","tanl","z0","dtot","ltot","chi2","chi2rz","chi2rphi",
-  // "bendchi2","pt","eta","nlaymiss_interior","phi","bendchi2_bin","chi2rz_bin","chi2rphi_bin"}
+  // "bendchi2","pt","eta","nlaymiss_interior","phi","bendchi2_bin","chi2rz_bin","chi2rphi_bin",
+  // "nlaymiss_PS","nlaymiss_2S"}
 
   std::vector<float> transformedFeatures;
 
@@ -66,16 +67,14 @@ std::vector<float> TrackQuality::featureTransform(TTTrack<Ref_Phase2TrackerDigi_
 
   // iterate through bits of the hitpattern and compare to 1 filling the hitpattern binary vector
   int tmp_trk_hitpattern = aTrack.hitPattern();
-  for (int i = 6; i >= 0; i--) {
-    int k = tmp_trk_hitpattern >> i;
-    if (k & 1)
-      hitpattern_binary[i] = 1;
-  }
-
   // calculate number of missed interior layers from hitpattern
   int nbits = floor(log2(tmp_trk_hitpattern)) + 1;
   int lay_i = 0;
   int tmp_trk_nlaymiss_interior = 0;
+  int tmp_trk_nlaymiss_PS = 0;
+  int tmp_trk_nlaymiss_2S = 0;
+  double tmp_trk_tanL = aTrack.tanL();
+  double tmp_trk_z0 = aTrack.z0();
   bool seq = false;
   for (int i = 0; i < nbits; i++) {
     lay_i = ((1 << i) & tmp_trk_hitpattern) >> i;  //0 or 1 in ith bit (right to left)
@@ -86,18 +85,29 @@ std::vector<float> TrackQuality::featureTransform(TTTrack<Ref_Phase2TrackerDigi_
       tmp_trk_nlaymiss_interior++;
   }
 
-  float eta = abs(aTrack.eta());
-  int eta_size = static_cast<int>(eta_bins.size());
-  // First iterate through eta bins
-
-  for (int j = 1; j < eta_size; j++) {
-    if (eta < eta_bins[j] && eta >= eta_bins[j - 1])  // if track in eta bin
-    {
-      // Iterate through hitpattern binary
-      for (int k = 0; k <= 6; k++)
-        // Fill expanded binary entries using the expected hitmap table positions
-        hitpattern_expanded_binary[hitmap[j - 1][k]] = hitpattern_binary[k];
-      break;
+  if (useHPH) {
+    hph::HitPatternHelper hph(setup_, tmp_trk_hitpattern, tmp_trk_tanL, tmp_trk_z0);
+    hitpattern_expanded_binary = hph.binary();
+    tmp_trk_nlaymiss_PS = hph.numMissingPS();
+    tmp_trk_nlaymiss_2S = hph.numMissing2S();
+  } else {
+    for (int i = 6; i >= 0; i--) {
+      int k = tmp_trk_hitpattern >> i;
+      if (k & 1)
+        hitpattern_binary[i] = 1;
+    }
+    float eta = abs(aTrack.eta());
+    int eta_size = static_cast<int>(eta_bins.size());
+    // First iterate through eta bins
+    for (int j = 1; j < eta_size; j++) {
+      if (eta < eta_bins[j] && eta >= eta_bins[j - 1])  // if track in eta bin
+      {
+        // Iterate through hitpattern binary
+        for (int k = 0; k <= 6; k++)
+          // Fill expanded binary entries using the expected hitmap table positions
+          hitpattern_expanded_binary[hitmap[j - 1][k]] = hitpattern_binary[k];
+        break;
+      }
     }
   }
 
@@ -162,8 +172,8 @@ std::vector<float> TrackQuality::featureTransform(TTTrack<Ref_Phase2TrackerDigi_
   // fill the feature map
   feature_map["nstub"] = stubRefs.size();
   feature_map["rinv"] = ONNXInvRScaling_ * abs(aTrack.rInv());
-  feature_map["tanl"] = abs(aTrack.tanL());
-  feature_map["z0"] = aTrack.z0();
+  feature_map["tanl"] = abs(tmp_trk_tanL);
+  feature_map["z0"] = tmp_trk_z0;
   feature_map["phi"] = aTrack.phi();
   feature_map["pt"] = aTrack.momentum().perp();
   feature_map["eta"] = aTrack.eta();
@@ -204,6 +214,10 @@ std::vector<float> TrackQuality::featureTransform(TTTrack<Ref_Phase2TrackerDigi_
   feature_map["chi2rphi_bin"] = tmp_trk_chi2rphi_bin;
   feature_map["chi2rz_bin"] = tmp_trk_chi2rz_bin;
 
+  //Bonus features from hitpattern
+  feature_map["nlaymiss_PS"] = float(tmp_trk_nlaymiss_PS);
+  feature_map["nlaymiss_2S"] = float(tmp_trk_nlaymiss_2S);
+
   // fill tensor with track params
   transformedFeatures.reserve(featureNames.size());
   for (const std::string& feature : featureNames)
@@ -304,3 +318,8 @@ void TrackQuality::setONNXModel(std::string const& AlgorithmString,
   ONNXInputName_ = ONNXInputName;
   featureNames_ = featureNames;
 }
+
+void TrackQuality::setHPHSetup(const hph::Setup* setup) {
+  setup_ = setup;
+  useHPH = true;
+}