[13.1.X] Introduce missing hits recovery in the SiStripHitEfficiency PCL and CalibTree workflows

CalibTracker/SiStripHitEfficiency/plugins/HitEff.cc

@@ -104,6 +104,10 @@ HitEff::HitEff(const edm::ParameterSet& conf)
   useLastMeas_ = conf_.getUntrackedParameter<bool>("useLastMeas", false);
   useAllHitsFromTracksWithMissingHits_ =
       conf_.getUntrackedParameter<bool>("useAllHitsFromTracksWithMissingHits", false);
+  doMissingHitsRecovery_ = conf_.getUntrackedParameter<bool>("doMissingHitsRecovery", false);
+
+  hitRecoveryCounters.resize(k_END_OF_LAYERS, 0);
+  hitTotalCounters.resize(k_END_OF_LAYERS, 0);
 }
 
 void HitEff::beginJob() {
@@ -163,6 +167,9 @@ void HitEff::beginJob() {
 
   events = 0;
   EventTrackCKF = 0;
+
+  totalNbHits = 0;
+  missHitPerLayer.resize(k_END_OF_LAYERS, 0);
 }
 
 void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
@@ -358,6 +365,7 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
       highPurity = itrack->quality(reco::TrackBase::TrackQuality::highPurity);
 
       std::vector<TrajectoryMeasurement> TMeas = itraj->measurements();
+      totalNbHits += int(TMeas.size());
       vector<TrajectoryMeasurement>::iterator itm;
       double xloc = 0.;
       double yloc = 0.;
@@ -369,15 +377,54 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
 
       // Check whether the trajectory has some missing hits
       bool hasMissingHits = false;
-      for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {
-        auto theHit = (*itm).recHit();
+      int previous_layer = 999;
+      vector<unsigned int> missedLayers;
+      for (const auto& itm : TMeas) {
+        auto theHit = itm.recHit();
+        unsigned int iidd = theHit->geographicalId().rawId();
+        int layer = ::checkLayer(iidd, tTopo);
+        int missedLayer = (layer + 1);
+        int diffPreviousLayer = (layer - previous_layer);
+        if (doMissingHitsRecovery_) {
+          //Layers from TIB + TOB
+          if (diffPreviousLayer == -2 && missedLayer > k_LayersStart && missedLayer < k_LayersAtTOBEnd) {
+            missHitPerLayer[missedLayer] += 1;
+            hasMissingHits = true;
+          }
+          //Layers from TID
+          else if (diffPreviousLayer == -2 && (missedLayer > k_LayersAtTOBEnd + 1 && missedLayer <= k_LayersAtTIDEnd)) {
+            missHitPerLayer[missedLayer] += 1;
+            hasMissingHits = true;
+          }
+          //Layers from TEC
+          else if (diffPreviousLayer == -2 && missedLayer > k_LayersAtTIDEnd && missedLayer <= k_LayersAtTECEnd) {
+            missHitPerLayer[missedLayer] += 1;
+            hasMissingHits = true;
+          }
+
+          //##### TID Layer 11 in transition TID -> TIB : layer is in TIB, previous layer  = 12
+          if ((layer > k_LayersStart && layer <= k_LayersAtTIBEnd) && (previous_layer == 12)) {
+            missHitPerLayer[11] += 1;
+            hasMissingHits = true;
+          }
+
+          //##### TEC Layer 14 in transition TEC -> TOB : layer is in TOB, previous layer =  15
+          if ((layer > k_LayersAtTIBEnd && layer <= k_LayersAtTOBEnd) && (previous_layer == 15)) {
+            missHitPerLayer[14] += 1;
+            hasMissingHits = true;
+          }
+        }
         if (theHit->getType() == TrackingRecHit::Type::missing)
           hasMissingHits = true;
+
+        if (hasMissingHits)
+          missedLayers.push_back(layer);
+        previous_layer = layer;
       }
 
       // Loop on each measurement and take it into consideration
       //--------------------------------------------------------
-
+      unsigned int prev_TKlayers = 0;
       for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {
         auto theInHit = (*itm).recHit();
 
@@ -430,12 +477,141 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
           // the matched layer should be added to the study as well
           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 1));
           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator, 2));
-          LogDebug("SiStripHitEfficiency:HitEff") << " found a hit with a missing partner" << endl;
+          LogDebug("SiStripHitEfficiency:HitEff") << " found a hit with a missing partner";
         } else {
           //only add one TM for the single surface and the other will be added in the next iteration
           TMs.push_back(TrajectoryAtInvalidHit(*itm, tTopo, tkgeom, propagator));
         }
+        bool missingHitAdded = false;
+
+        vector<TrajectoryMeasurement> tmpTmeas;
+        unsigned int misLayer = TKlayers + 1;
+        //Use bool doMissingHitsRecovery to add possible missing hits based on actual/previous hit
+        if (doMissingHitsRecovery_) {
+          if (int(TKlayers) - int(prev_TKlayers) == -2) {
+            const DetLayer* detlayer = itm->layer();
+            const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
+            const TrajectoryStateOnSurface tsos = itm->updatedState();
+            std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
+
+            if (misLayer > k_LayersAtTIDEnd && misLayer < k_LayersAtTECEnd) {  //TEC
+              std::vector<ForwardDetLayer const*> negTECLayers =
+                  measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
+              std::vector<ForwardDetLayer const*> posTECLayers =
+                  measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
+              const DetLayer* tecLayerneg = negTECLayers[misLayer - k_LayersAtTIDEnd - 1];
+              const DetLayer* tecLayerpos = posTECLayers[misLayer - k_LayersAtTIDEnd - 1];
+              if (tTopo->tecSide(iidd) == 1) {
+                tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
+              } else if (tTopo->tecSide(iidd) == 2) {
+                tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
+              }
+            }
+
+            else if (misLayer == (k_LayersAtTIDEnd - 1) ||
+                     misLayer == k_LayersAtTIDEnd) {  // This is for TID layers 12 and 13
+
+              std::vector<ForwardDetLayer const*> negTIDLayers =
+                  measurementTrackerHandle->geometricSearchTracker()->negTidLayers();
+              std::vector<ForwardDetLayer const*> posTIDLayers =
+                  measurementTrackerHandle->geometricSearchTracker()->posTidLayers();
+              const DetLayer* tidLayerneg = negTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
+              const DetLayer* tidLayerpos = posTIDLayers[misLayer - k_LayersAtTOBEnd - 1];
+
+              if (tTopo->tidSide(iidd) == 1) {
+                tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, *thePropagator, *estimator);
+              } else if (tTopo->tidSide(iidd) == 2) {
+                tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, *thePropagator, *estimator);
+              }
+            }
 
+            if (misLayer > k_LayersStart && misLayer < k_LayersAtTOBEnd) {  // Barrel
+
+              std::vector<BarrelDetLayer const*> barrelTIBLayers =
+                  measurementTrackerHandle->geometricSearchTracker()->tibLayers();
+              std::vector<BarrelDetLayer const*> barrelTOBLayers =
+                  measurementTrackerHandle->geometricSearchTracker()->tobLayers();
+
+              if (misLayer > k_LayersStart && misLayer <= k_LayersAtTIBEnd) {
+                const DetLayer* tibLayer = barrelTIBLayers[misLayer - k_LayersStart - 1];
+                tmpTmeas = layerMeasurements.measurements(*tibLayer, tsos, *thePropagator, *estimator);
+              } else if (misLayer > k_LayersAtTIBEnd && misLayer < k_LayersAtTOBEnd) {
+                const DetLayer* tobLayer = barrelTOBLayers[misLayer - k_LayersAtTIBEnd - 1];
+                tmpTmeas = layerMeasurements.measurements(*tobLayer, tsos, *thePropagator, *estimator);
+              }
+            }
+          }
+          if ((int(TKlayers) > k_LayersStart && int(TKlayers) <= k_LayersAtTIBEnd) && int(prev_TKlayers) == 12) {
+            const DetLayer* detlayer = itm->layer();
+            const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
+            const TrajectoryStateOnSurface tsos = itm->updatedState();
+            std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
+            std::vector<ForwardDetLayer const*> negTIDLayers =
+                measurementTrackerHandle->geometricSearchTracker()->negTidLayers();
+            std::vector<ForwardDetLayer const*> posTIDLayers =
+                measurementTrackerHandle->geometricSearchTracker()->posTidLayers();
+
+            const DetLayer* tidLayerneg = negTIDLayers[k_LayersStart];
+            const DetLayer* tidLayerpos = posTIDLayers[k_LayersStart];
+            if (tTopo->tidSide(iidd) == 1) {
+              tmpTmeas = layerMeasurements.measurements(*tidLayerneg, tsos, *thePropagator, *estimator);
+            } else if (tTopo->tidSide(iidd) == 2) {
+              tmpTmeas = layerMeasurements.measurements(*tidLayerpos, tsos, *thePropagator, *estimator);
+            }
+          }
+
+          if ((int(TKlayers) > k_LayersAtTIBEnd && int(TKlayers) <= k_LayersAtTOBEnd) && int(prev_TKlayers) == 15) {
+            const DetLayer* detlayer = itm->layer();
+            const LayerMeasurements layerMeasurements{*measurementTrackerHandle, *measurementTrackerEvent};
+            const TrajectoryStateOnSurface tsos = itm->updatedState();
+            std::vector<DetLayer::DetWithState> compatDets = detlayer->compatibleDets(tsos, *thePropagator, *estimator);
+
+            std::vector<ForwardDetLayer const*> negTECLayers =
+                measurementTrackerHandle->geometricSearchTracker()->negTecLayers();
+            std::vector<ForwardDetLayer const*> posTECLayers =
+                measurementTrackerHandle->geometricSearchTracker()->posTecLayers();
+
+            const DetLayer* tecLayerneg = negTECLayers[k_LayersStart];
+            const DetLayer* tecLayerpos = posTECLayers[k_LayersStart];
+            if (tTopo->tecSide(iidd) == 1) {
+              tmpTmeas = layerMeasurements.measurements(*tecLayerneg, tsos, *thePropagator, *estimator);
+            } else if (tTopo->tecSide(iidd) == 2) {
+              tmpTmeas = layerMeasurements.measurements(*tecLayerpos, tsos, *thePropagator, *estimator);
+            }
+          }
+
+          if (!tmpTmeas.empty()) {
+            TrajectoryMeasurement TM_tmp(tmpTmeas.back());
+            unsigned int iidd_tmp = TM_tmp.recHit()->geographicalId().rawId();
+            if (iidd_tmp != 0) {
+              LogDebug("SiStripHitEfficiency:HitEff") << " hit actually being added to TM vector";
+              if ((!useAllHitsFromTracksWithMissingHits_ || (!useFirstMeas_ && isFirstMeas)))
+                TMs.clear();
+              if (::isDoubleSided(iidd_tmp, tTopo)) {
+                TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 1));
+                TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator, 2));
+              } else
+                TMs.push_back(TrajectoryAtInvalidHit(TM_tmp, tTopo, tkgeom, propagator));
+              missingHitAdded = true;
+              hitRecoveryCounters[misLayer] += 1;
+            }
+          }
+        }
+
+        prev_TKlayers = TKlayers;
+        if (!useFirstMeas_ && isFirstMeas && !missingHitAdded)
+          continue;
+        if (!useLastMeas_ && isLastMeas)
+          continue;
+        bool hitsWithBias = false;
+        for (auto ilayer : missedLayers) {
+          if (ilayer < TKlayers)
+            hitsWithBias = true;
+        }
+        if (hasMissingHits && theInHit->getType() != TrackingRecHit::Type::missing && !missingHitAdded &&
+            hitsWithBias && !useAllHitsFromTracksWithMissingHits_) {
+          continue;
+        }
         //////////////////////////////////////////////
         //Now check for tracks at TOB6 and TEC9
 
@@ -537,6 +713,7 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
             }
           }  //else LOGPRINT << "tec9 tmp empty" << endl;
         }
+        hitTotalCounters[TKlayers] += 1;
 
         ////////////////////////////////////////////////////////
 
@@ -872,6 +1049,93 @@ void HitEff::endJob() {
 
   LogDebug("SiStripHitEfficiency:HitEff") << " Events Analysed             " << events << endl;
   LogDebug("SiStripHitEfficiency:HitEff") << " Number Of Tracked events    " << EventTrackCKF << endl;
+
+  if (doMissingHitsRecovery_) {
+    float totTIB = 0.0;
+    float totTOB = 0.0;
+    float totTID = 0.0;
+    float totTEC = 0.0;
+
+    float totTIBrepro = 0.0;
+    float totTOBrepro = 0.0;
+    float totTIDrepro = 0.0;
+    float totTECrepro = 0.0;
+
+    edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TIB :";
+    for (int i = 0; i <= k_LayersAtTIBEnd; i++) {
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
+          << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
+      totTIB += missHitPerLayer[i];
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
+          << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
+          << " % of missing hit";
+      totTIBrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
+    }
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "TOTAL % of missing hits within TIB :" << (totTIB * 1.0 / totalNbHits) * 100 << "%";
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "AFTER repropagation :" << (totTIBrepro * 1.0 / totalNbHits) * 100 << "%";
+
+    edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TOB :";
+    for (int i = k_LayersAtTIBEnd + 1; i <= k_LayersAtTOBEnd; i++) {
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
+          << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
+      totTOB += missHitPerLayer[i];
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
+          << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
+          << " % of missing hit";
+      totTOBrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
+    }
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "TOTAL % of missing hits within TOB :" << (totTOB * 1.0 / totalNbHits) * 100 << "%";
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "AFTER repropagation :" << (totTOBrepro * 1.0 / totalNbHits) * 100 << "%";
+
+    edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TID :";
+    for (int i = k_LayersAtTOBEnd + 1; i <= k_LayersAtTIDEnd; i++) {
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
+          << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
+      totTID += missHitPerLayer[i];
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
+          << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
+          << " % of missing hit";
+      totTIDrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
+    }
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "TOTAL % of missing hits within TID :" << (totTID * 1.0 / totalNbHits) * 100 << "%";
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "AFTER repropagation :" << (totTIDrepro * 1.0 / totalNbHits) * 100 << "%";
+
+    edm::LogInfo("SiStripHitEfficiency:HitEff") << "Within TEC :";
+    for (int i = k_LayersAtTIDEnd + 1; i < k_END_OF_LAYERS; i++) {
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Layer " << i << " has : " << missHitPerLayer[i] << "/" << totalNbHits << " = "
+          << (missHitPerLayer[i] * 1.0 / totalNbHits) * 100 << " % of missing hit";
+      totTEC += missHitPerLayer[i];
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << "Removing recovered hits : layer " << i << " has : " << missHitPerLayer[i] - hitRecoveryCounters[i] << "/"
+          << totalNbHits << " = " << ((missHitPerLayer[i] - hitRecoveryCounters[i]) * 1.0 / totalNbHits) * 100
+          << " % of missing hit";
+      totTECrepro += (missHitPerLayer[i] - hitRecoveryCounters[i]);
+    }
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "TOTAL % of missing hits within TEC :" << (totTEC * 1.0 / totalNbHits) * 100 << "%";
+    edm::LogInfo("SiStripHitEfficiency:HitEff")
+        << "AFTER repropagation :" << (totTECrepro * 1.0 / totalNbHits) * 100 << "%";
+
+    edm::LogInfo("SiStripHitEfficiency:HitEff") << " Hit recovery summary:";
+
+    for (int ilayer = 0; ilayer < k_END_OF_LAYERS; ilayer++) {
+      edm::LogInfo("SiStripHitEfficiency:HitEff")
+          << " layer " << ilayer << ": " << hitRecoveryCounters[ilayer] << " / " << hitTotalCounters[ilayer];
+    }
+  }
 }
 
 //define this as a plug-in

CalibTracker/SiStripHitEfficiency/plugins/HitEff.h

@@ -74,6 +74,7 @@ class HitEff : public edm::one::EDAnalyzer<edm::one::SharedResources> {
   bool useFirstMeas_;
   bool useLastMeas_;
   bool useAllHitsFromTracksWithMissingHits_;
+  bool doMissingHitsRecovery_;
 
   const edm::EDGetTokenT<reco::TrackCollection> combinatorialTracks_token_;
   const edm::EDGetTokenT<std::vector<Trajectory> > trajectories_token_;
@@ -103,6 +104,8 @@ class HitEff : public edm::one::EDAnalyzer<edm::one::SharedResources> {
   bool DEBUG;
   unsigned int whatlayer;
 
+  std::vector<unsigned int> hitRecoveryCounters;
+  std::vector<unsigned int> hitTotalCounters;
 // Tree declarations
 // Trajectory positions for modules included in the study
 #ifdef ExtendedCALIBTree
@@ -113,6 +116,8 @@ class HitEff : public edm::one::EDAnalyzer<edm::one::SharedResources> {
   int nLostHits;
   float p, chi2;
 #endif
+  int totalNbHits;
+  std::vector<int> missHitPerLayer;
   float TrajGlbX, TrajGlbY, TrajGlbZ;
   float TrajLocX, TrajLocY, TrajLocAngleX, TrajLocAngleY;
   float TrajLocErrX, TrajLocErrY;

CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc

@@ -401,6 +401,7 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
             const auto num = h_module_found->getValue(det);
             const auto denom = h_module_total->getValue(det);
             if (denom) {
+              assert(num <= denom);  // can't have this happen
               const auto eff = num / denom;
               const auto eff_up = TEfficiency::Bayesian(denom, num, .99, 1, 1, true);
 
@@ -465,7 +466,7 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
     LOGPRINT << "Number of strips module " << det << " is " << nStrips;
     badStripList.push_back(pQuality.encode(0, nStrips, 0));
     //Now compact into a single bad module
-    LOGPRINT << "ID1 shoudl match list of modules above " << det;
+    LOGPRINT << "ID1 should match list of modules above " << det;
     pQuality.compact(det, badStripList);
     pQuality.put(det, SiStripQuality::Range(badStripList.begin(), badStripList.end()));
   }