Fix bugs in HYBRID_NEWKF hit pattern

L1Trigger/TrackFindingTracklet/interface/HitPatternHelper.h

@@ -16,6 +16,9 @@
 // considers truncaton effect and is bit/clock cycle accurate.
 // With this version of the KF, HitPatternHelper makes predictions based on the spatial coordinates of tracks and detector modules.
 //
+// WARNING: This code needs major revision. There is several bugs, and half the functions here only work
+// correctly for HYBRID_NEWKF, and the other half only work correctly for HYBRID.
+//
 //  Created by J.Li on 1/23/21.
 //
 
@@ -66,9 +69,15 @@ namespace hph {
     int etaRegion(double z0, double cot, bool useNewKF) const;
     const std::vector<int>& layerEncoding(double zT) const { return layerEncoding_->layerEncoding(zT); }
     // LayerEncoding call filling numPS, num2S, numMissingPS and numMissingPS for given hitPattern and trajectory
-    void analyze(
-        int hitpattern, double cot, double z0, int& numPS, int& num2S, int& numMissingPS, int& numMissing2S) const {
-      layerEncoding_->analyze(hitpattern, cot, z0, numPS, num2S, numMissingPS, numMissing2S);
+    void analyze(int hitpattern,
+                 double cot,
+                 double z0,
+                 int& rzSect,
+                 int& numPS,
+                 int& num2S,
+                 int& numMissingPS,
+                 int& numMissing2S) const {
+      layerEncoding_->analyze(hitpattern, cot, z0, rzSect, numPS, num2S, numMissingPS, numMissing2S);
     }
 
   private:
@@ -117,6 +126,7 @@ namespace hph {
     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
+    bool newKF() const { return useNewKF_; }
 
   private:
     const Setup* setup_;

L1Trigger/TrackFindingTracklet/src/HitPatternHelper.cc

@@ -27,7 +27,7 @@ namespace hph {
         nKalmanLayers_(tmtt::KFbase::nKFlayer_) {
     if (useNewKF_) {
       chosenRofZ_ = chosenRofZNewKF_;
-      etaRegions_ = etaRegionsNewKF_;
+      etaRegions_ = etaRegionsNewKF_;  // TO FIX: Not defined
     } else {
       chosenRofZ_ = iConfig.chosenRofZ_;
       etaRegions_ = iConfig.etaRegions_;
@@ -60,7 +60,7 @@ namespace hph {
         layerEncoding_(setup->layerEncoding(zT_)),
         numExpLayer_(layerEncoding_.size()),
         hitpattern_(hitpattern),
-        etaSector_(setup_->etaRegion(z0, cot, useNewKF_)),
+        etaSector_(setup_->etaRegion(z0, cot, useNewKF_)),  // Only works for old KF
         numMissingLayer_(0),
         numMissingPS_(0),
         numMissing2S_(0),
@@ -70,14 +70,17 @@ namespace hph {
         numMissingInterior2_(0),
         binary_(11, 0),  //there are 11 unique layer IDs, as defined in variable "layerIds"
         bonusFeatures_() {
-    setup->analyze(hitpattern, cot, z0, numPS_, num2S_, numMissingPS_, numMissing2S_);
+    // FIX: Function analyze is only designed to work correctly for NEWKF.
+    int rzSect = 0;
+    setup->analyze(hitpattern, cot, z0, rzSect, numPS_, num2S_, numMissingPS_, numMissing2S_);
+    if (useNewKF_)
+      etaSector_ = rzSect;
     int kf_eta_reg = etaSector_;
     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;
     }
-
     int nbits = floor(log2(hitpattern_)) + 1;
     int lay_i = 0;
     bool seq = false;
@@ -105,36 +108,40 @@ namespace hph {
 
     if (hphDebug_) {
       if (useNewKF_) {
-        edm::LogVerbatim("TrackTriggerHPH") << "Running with New KF";
+        edm::LogVerbatim("Tracklet") << "Running with New KF";
       } else {
-        edm::LogVerbatim("TrackTriggerHPH") << "Running with Old KF";
+        edm::LogVerbatim("Tracklet") << "Running with Old KF";
       }
-      edm::LogVerbatim("TrackTriggerHPH") << "======================================================";
-      edm::LogVerbatim("TrackTriggerHPH")
-          << "Looking at hitpattern " << std::bitset<7>(hitpattern_) << "; Looping over KF layers:";
+      edm::LogVerbatim("Tracklet") << "======================================================";
+      edm::LogVerbatim("Tracklet") << "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;
+          edm::LogVerbatim("Tracklet") << "--------------------------";
+          edm::LogVerbatim("Tracklet") << "Looking at KF layer " << i;
           if (layerEncoding_[i] < 10) {
-            edm::LogVerbatim("TrackTriggerHPH") << "KF expects L" << layerEncoding_[i];
+            edm::LogVerbatim("Tracklet") << "KF expects L" << layerEncoding_[i];
           } else {
-            edm::LogVerbatim("TrackTriggerHPH") << "KF expects D" << layerEncoding_[i] - 10;
+            edm::LogVerbatim("Tracklet") << "KF expects D" << layerEncoding_[i] - 10;
           }
         }
 
         if (((1 << i) & hitpattern_) >> i) {
           if (hphDebug_) {
-            edm::LogVerbatim("TrackTriggerHPH") << "Layer found in hitpattern";
+            edm::LogVerbatim("Tracklet") << "Layer found in hitpattern";
+          }
+          if (layerEncoding_[i] < 0) {
+            edm::LogError("Tracklet") << "Track's hit-pattern shows stub in non-traversed layer";
+            continue;
           }
           binary_[reducedId(layerEncoding_[i])] = 1;
         } else {
           if (hphDebug_) {
-            edm::LogVerbatim("TrackTriggerHPH") << "Layer missing in hitpattern";
+            edm::LogVerbatim("Tracklet") << "Layer missing in hitpattern";
           }
         }
       }
@@ -144,13 +151,13 @@ namespace hph {
       for (int i = 0; i < nKalmanLayers_; i++) {  //Loop over each digit of hitpattern
 
         if (hphDebug_) {
-          edm::LogVerbatim("TrackTriggerHPH") << "--------------------------";
-          edm::LogVerbatim("TrackTriggerHPH") << "Looking at KF layer " << i;
+          edm::LogVerbatim("Tracklet") << "--------------------------";
+          edm::LogVerbatim("Tracklet") << "Looking at KF layer " << i;
         }
 
         if (layermap_[kf_eta_reg][i].empty()) {
           if (hphDebug_) {
-            edm::LogVerbatim("TrackTriggerHPH") << "KF does not expect this layer";
+            edm::LogVerbatim("Tracklet") << "KF does not expect this layer";
           }
 
           continue;
@@ -161,45 +168,45 @@ namespace hph {
 
           if (hphDebug_) {
             if (j < 10) {
-              edm::LogVerbatim("TrackTriggerHPH") << "KF expects L" << j;
+              edm::LogVerbatim("Tracklet") << "KF expects L" << j;
             } else {
-              edm::LogVerbatim("TrackTriggerHPH") << "KF expects D" << j - 10;
+              edm::LogVerbatim("Tracklet") << "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";
+              edm::LogVerbatim("Tracklet") << "Rejected by sensor modules";
             }
 
             continue;
           }
 
           if (hphDebug_) {
-            edm::LogVerbatim("TrackTriggerHPH") << "Confirmed by sensor modules";
+            edm::LogVerbatim("Tracklet") << "Confirmed by sensor modules";
           }
           //prediction is accepted
           if (((1 << i) & hitpattern_) >> i) {
             if (hphDebug_) {
-              edm::LogVerbatim("TrackTriggerHPH") << "Layer found in hitpattern";
+              edm::LogVerbatim("Tracklet") << "Layer found in hitpattern";
             }
             binary_[reducedId(j)] = 1;
           } else {
             if (hphDebug_) {
-              edm::LogVerbatim("TrackTriggerHPH") << "Layer missing in hitpattern";
+              edm::LogVerbatim("Tracklet") << "Layer missing in hitpattern";
             }
           }
         }
       }
     }
 
     if (hphDebug_) {
-      edm::LogVerbatim("TrackTriggerHPH") << "------------------------------";
-      edm::LogVerbatim("TrackTriggerHPH") << "numPS = " << numPS_ << ", num2S = " << num2S_
-                                          << ", missingPS = " << numMissingPS_ << ", missing2S = " << numMissing2S_;
-      edm::LogVerbatim("TrackTriggerHPH") << "======================================================";
+      edm::LogVerbatim("Tracklet") << "------------------------------";
+      edm::LogVerbatim("Tracklet") << "numPS = " << numPS_ << ", num2S = " << num2S_
+                                   << ", missingPS = " << numMissingPS_ << ", missing2S = " << numMissing2S_;
+      edm::LogVerbatim("Tracklet") << "======================================================";
     }
   }
 
@@ -228,7 +235,7 @@ namespace hph {
 
   int HitPatternHelper::reducedId(int layerId) {
     if (hphDebug_ && (layerId > 15 || layerId < 1)) {
-      edm::LogVerbatim("TrackTriggerHPH") << "Warning: invalid layer id !";
+      edm::LogVerbatim("Tracklet") << "Warning: invalid layer id !";
     }
     if (layerId <= 6) {
       layerId = layerId - 1;

L1Trigger/TrackFindingTracklet/src/TrackFindingProcessor.cc

@@ -81,7 +81,9 @@ namespace trklet {
       return;
     // create bit vectors
     std::string s = trackTQ.hitPattern().str();
-    s.resize(TTTrack_TrackWord::TrackBitWidths::kHitPatternSize);
+    std::reverse(s.begin(), s.end());
+    // Drop outermost (8th) track layer, as data format foresees only 7 bits.
+    s.erase(0, 1);
     hitPattern_ = TTBV(s);
     const TTBV other = TTBV(0, 2 * TTTrack_TrackWord::TrackBitWidths::kMVAQualitySize);
     const TTBV chi2bend = TTBV(0, TTTrack_TrackWord::TrackBitWidths::kBendChi2Size);

L1Trigger/TrackFindingTracklet/src/TrackQuality.cc

@@ -89,7 +89,7 @@ namespace trklet {
       // transform double to AP_FIXED_BDT
       static const double d = std::pow(2., 10);
       const std::vector<AP_FIXED_BDT> features({nStubs, zT / d, cot / d, chi20 / d, chi21 / d, nGaps});
-      // BDT Inference
+      // Run the Track Quality BDT calculation
       const AP_FIXED_BDT mvaFixed = bdt_->decision_function(features).at(0);
       const AP_INT_BDT mvaInt = mvaFixed.range(mvaFixed.width - 1, 0);
       // bin mva

L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker.cc

@@ -177,6 +177,7 @@ class L1TrackNtupleMaker : public one::EDAnalyzer<one::WatchRuns, one::SharedRes
   std::vector<int>* m_trk_lhits;
   std::vector<int>* m_trk_dhits;
   std::vector<int>* m_trk_seed;
+  // WARNING - info unpacked from the hit pattern is unrelable due to bugs in HitPatternHelper
   std::vector<int>* m_trk_hitpattern;
   std::vector<int>* m_trk_lhits_hitpattern;  // 6-digit hit mask (barrel layer only) dervied from hitpattern
   std::vector<int>* m_trk_dhits_hitpattern;  // disk only
@@ -1101,6 +1102,7 @@ void L1TrackNtupleMaker::analyze(const edm::Event& iEvent, const edm::EventSetup
 
       int tmp_trk_hitpattern = 0;
       tmp_trk_hitpattern = (int)iterL1Track->hitPattern();
+      // TO FIX -- The HitPatternHelper code contains several bugs & design flaws.
       hph::HitPatternHelper hph(hphSetup, tmp_trk_hitpattern, tmp_trk_tanL, tmp_trk_z0);
       std::vector<int> hitpattern_expanded_binary = hph.binary();
       int tmp_trk_lhits_hitpattern = 0;
@@ -1165,6 +1167,8 @@ void L1TrackNtupleMaker::analyze(const edm::Event& iEvent, const edm::EventSetup
       //float tmp_trk_bend_chi2 = 0;
       int tmp_trk_dhits = 0;
       int tmp_trk_lhits = 0;
+      set<int> hitLayPS;
+      set<int> hitLay2S;
 
       if (true) {
         // loop over stubs
@@ -1181,22 +1185,53 @@ void L1TrackNtupleMaker::analyze(const edm::Event& iEvent, const edm::EventSetup
           double z = posStub.z();
 
           int layer = -999999;
-          if (detIdStub.subdetId() == StripSubdetector::TOB) {
+          bool barrel = (detIdStub.subdetId() == StripSubdetector::TOB);
+          if (barrel) {
             layer = static_cast<int>(tTopo->layer(detIdStub));
             if (DebugMode)
               edm::LogVerbatim("Tracklet")
                   << "   stub in layer " << layer << " at position x y z = " << x << " " << y << " " << z;
             tmp_trk_lhits += pow(10, layer - 1);
-          } else if (detIdStub.subdetId() == StripSubdetector::TID) {
+          } else {
             layer = static_cast<int>(tTopo->layer(detIdStub));
             if (DebugMode)
               edm::LogVerbatim("Tracklet")
                   << "   stub in disk " << layer << " at position x y z = " << x << " " << y << " " << z;
             tmp_trk_dhits += pow(10, layer - 1);
           }
 
+          bool psMod = (theTrackerGeom->getDetectorType(detIdStub) == TrackerGeometry::ModuleType::Ph2PSP);
+          int layerdisk = barrel ? layer : 10 + layer;
+          if (psMod) {
+            hitLayPS.insert(layerdisk);
+          } else {
+            hitLay2S.insert(layerdisk);
+          }
         }  //end loop over stubs
       }
+
+      // Check accuracy of hit pattern info.
+      if (DebugMode && hph.newKF()) {
+        // Only bother for New KF, as HitPatternHelper buggy for Old KF.
+        if (hph.numPS() != int(hitLayPS.size()) || hph.num2S() != int(hitLay2S.size())) {
+          // Some inaccuracy expected, as estimating number of PS stubs from hit pattern is approximate,
+          // due to r-boundary of PS-2S not being well defined, and due to truncation of hit pattern
+          // from 8 to 7 bits.
+          std::stringstream ss;
+          ss << "Number of layers with stubs estimated from hit pattern is inaccurate: (PS,2S) = (" << hph.numPS()
+             << "," << hph.num2S() << ") vs (" << hitLayPS.size() << "," << hitLay2S.size()
+             << "), eta sect = " << hph.etaSector() << ", hitpattern = " << std::bitset<8>(tmp_trk_hitpattern)
+             << ", hit layers = PS(";
+          for (const auto& lay : hitLayPS)
+            ss << " " << lay;
+          ss << ") + 2S(";
+          for (const auto& lay : hitLay2S)
+            ss << " " << lay;
+          ss << ")";
+          edm::LogWarning("Tracklet") << ss.str();
+        }
+      }
+
       // ----------------------------------------------------------------------------------------------
 
       int tmp_trk_genuine = 0;

L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py

@@ -40,8 +40,8 @@
 process.load('FWCore.MessageService.MessageLogger_cfi')
 process.MessageLogger.L1track = dict(limit = -1)
 process.MessageLogger.Tracklet = dict(limit = -1)
-process.MessageLogger.TrackTriggerHPH = dict(limit = -1)
-
+process.MessageLogger.cout.enableStatistics = True
+process.MessageLogger.cerr.enableStatistics = True
 
 print("using geometry " + GEOMETRY + " (tilted)")
 process.load('Configuration.Geometry.GeometryExtendedRun4' + GEOMETRY + 'Reco_cff')

L1Trigger/TrackerTFP/interface/LayerEncoding.h

@@ -32,9 +32,15 @@ namespace trackerTFP {
     int maybePS(int zT) const;
     // encoded layer id which may be PS or 2S for given zT in cm
     int maybePS(double zT) const;
-    // fills numPS, num2S, numMissingPS and numMissingPS for given hitPattern and trajectory
-    void analyze(
-        int hitpattern, double cot, double z0, int& numPS, int& num2S, int& numMissingPS, int& numMissing2S) const;
+    // fills binZT (unsigned), numPS, num2S, numMissingPS and numMissingPS for given hitPattern and trajectory
+    void analyze(int hitpattern,
+                 double cot,
+                 double z0,
+                 int& binZT,
+                 int& numPS,
+                 int& num2S,
+                 int& numMissingPS,
+                 int& numMissing2S) const;
 
   private:
     // helper class providing run-time constants

L1Trigger/TrackerTFP/src/LayerEncoding.cc

@@ -159,11 +159,18 @@ namespace trackerTFP {
     return maybePS(binZT);
   }
 
-  // fills numPS, num2S, numMissingPS and numMissingPS for given hitPattern and trajectory
-  void LayerEncoding::analyze(
-      int hitpattern, double cot, double z0, int& numPS, int& num2S, int& numMissingPS, int& numMissing2S) const {
+  // fills binZT (unsigned), numPS, num2S, numMissingPS and numMissingPS for given hitPattern and trajectory
+  void LayerEncoding::analyze(int hitpattern,
+                              double cot,
+                              double z0,
+                              int& binZT,
+                              int& numPS,
+                              int& num2S,
+                              int& numMissingPS,
+                              int& numMissing2S) const {
     // look up layer encoding nad maybe pattern
     const double zT = z0 + setup_->chosenRofZ() * cot;
+    binZT = zT_->toUnsigned(zT_->integer(zT));
     const std::vector<int>& le = this->layerEncoding(zT);
     const TTBV& mp = this->maybePattern(zT);
     const TTBV hp(hitpattern, setup_->numLayers());
@@ -179,17 +186,17 @@ namespace trackerTFP {
         const int diskId = layerId - setup_->offsetLayerDisks() - setup_->offsetLayerId();
         // avergae disk z position
         const double z = setup_->hybridDiskZ(diskId) * (cot < 0. ? -1. : 1.);
-        // innermost edge of 2S modules
-        const double rLimit = setup_->disk2SR(diskId, 0) - setup_->pitchCol2S();
-        // trajectory radius at avergae disk z position
+        // smallest stub radii from 2S disks
+        double rLimit = setup_->disk2SR(diskId, 0) - .5 * setup_->pitchCol2S();
+        // trajectory radius at average disk z position
         const double r = (z - z0) / cot;
         // compare with innermost edge of 2S modules to identify PS
         if (r < rLimit)
           ps = true;
       }
       if (hp.test(layerIdKF))  // layer is hit
         ps ? numPS++ : num2S++;
-      else if (!mp.test(layerIdKF))  // layer is not hit but should have been hitted (roughly by) trajectory
+      else if (!mp.test(layerIdKF))  // layer is not hit but should have been hit (roughly) by trajectory
         ps ? numMissingPS++ : numMissing2S++;
     }
   }