[12.5.X] SiStripHitEfficiency PCL workflow, take into account modules with FEDErrors

CalibTracker/SiStripHitEfficiency/interface/SiStripHitEffData.h

@@ -8,17 +8,63 @@
 
 struct SiStripHitEffData {
 public:
-  SiStripHitEffData() : EventStats(), FEDErrorOccupancy(nullptr) {}
+  SiStripHitEffData() : EventStats(), FEDErrorOccupancy(nullptr), m_isLoaded(false) {}
 
   void fillTkMapFromMap() {
     for (const auto& [id, count] : fedErrorCounts) {
-      FEDErrorOccupancy.fill(id, count);
+      FEDErrorOccupancy->fill(id, count);
+    }
+  }
+
+  void fillMapFromTkMap(const int nevents, const float threshold, const std::vector<DetId>& stripDetIds) {
+    const auto& Maps = FEDErrorOccupancy->getAllMaps();
+    std::vector<bool> isThere;
+    isThere.reserve(Maps.size());
+    std::transform(Maps.begin() + 1, Maps.end(), std::back_inserter(isThere), [](auto& x) { return !(x == nullptr); });
+
+    int count{0};
+    for (const auto& it : isThere) {
+      count++;
+      LogTrace("SiStripHitEffData") << " layer: " << count << " " << it << std::endl;
+      if (it)
+        LogTrace("SiStripHitEffData") << "resolving to " << Maps[count]->getName()
+                                      << " with entries: " << Maps[count]->getEntries() << std::endl;
+      // color the map
+      Maps[count]->setOption("colz");
+    }
+
+    for (const auto& det : stripDetIds) {
+      const auto& counts = FEDErrorOccupancy->getValue(det);
+
+      if (counts > 0) {
+        float fraction = counts / nevents;
+
+        LogTrace("SiStripHitEffData") << det.rawId() << " has " << counts << " counts, " << fraction * 100
+                                      << "% fraction of the " << nevents << " events processed" << std::endl;
+
+        if (fraction > threshold) {
+          fedErrorCounts.insert(std::make_pair(det, 1));
+        }
+      }  // do not check functioning modules
+    }
+    // the map has been loaded
+    m_isLoaded = true;
+  }
+
+  const bool checkFedError(const DetId det) {
+    if (m_isLoaded) {
+      return (fedErrorCounts.find(det) == fedErrorCounts.end());
+    } else {
+      throw cms::Exception("LogicError") << __PRETTY_FUNCTION__ << "cannot check DetId when map not filled";
     }
   }
 
   dqm::reco::MonitorElement* EventStats;
   std::unordered_map<uint32_t, int> fedErrorCounts;
-  TkHistoMap FEDErrorOccupancy;
+  std::unique_ptr<TkHistoMap> FEDErrorOccupancy;
+
+private:
+  bool m_isLoaded;
 };
 
 #endif

CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc

@@ -146,13 +146,12 @@ void SiStripHitEfficiencyHarvester::endRun(edm::Run const&, edm::EventSetup cons
 
 bool SiStripHitEfficiencyHarvester::checkMapsValidity(const std::vector<MonitorElement*>& maps,
                                                       const std::string& type) const {
-  std::vector<bool> isAvailable;
-  isAvailable.reserve(maps.size());
-  std::transform(
-      maps.begin() + 1, maps.end(), std::back_inserter(isAvailable), [](auto& x) { return !(x == nullptr); });
+  std::vector<bool> isThere;
+  isThere.reserve(maps.size());
+  std::transform(maps.begin() + 1, maps.end(), std::back_inserter(isThere), [](auto& x) { return !(x == nullptr); });
 
   int count{0};
-  for (const auto& it : isAvailable) {
+  for (const auto& it : isThere) {
     count++;
     LogDebug("SiStripHitEfficiencyHarvester") << " layer: " << count << " " << it << std::endl;
     if (it)
@@ -162,7 +161,7 @@ bool SiStripHitEfficiencyHarvester::checkMapsValidity(const std::vector<MonitorE
   // check on the input TkHistoMap
   bool areMapsAvailable{true};
   int layerCount{0};
-  for (const auto& it : isAvailable) {
+  for (const auto& it : isThere) {
     layerCount++;
     if (!it) {
       edm::LogError("SiStripHitEfficiencyHarvester")
@@ -264,13 +263,28 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
   TrackerMap tkMapDen{" Detector denominator "};
   std::map<unsigned int, double> badModules;
 
+  // load the FEDError map
+  const auto& EventStats = getter.get(fmt::format("{}/EventInfo/EventStats", inputFolder_));
+  const int totalEvents = EventStats->getBinContent(1., 1.);  // first bin contains info on number of events run
+  calibData_.FEDErrorOccupancy = std::make_unique<TkHistoMap>(tkDetMap_.get());
+  calibData_.FEDErrorOccupancy->loadTkHistoMap(fmt::format("{}/FEDErrorTkDetMaps", inputFolder_),
+                                               "perModule_FEDErrors");
+
+  // tag as bad from FEDErrors the modules that have an error on 75% of the events
+  calibData_.fillMapFromTkMap(totalEvents, 0.75, stripDetIds_);
+
+  for (const auto& [badId, fraction] : calibData_.fedErrorCounts) {
+    LogDebug("SiStripHitEfficiencyHarvester")
+        << __PRETTY_FUNCTION__ << " bad module from FEDError " << badId << "," << fraction << std::endl;
+  }
+
   for (auto det : stripDetIds_) {
     auto layer = ::checkLayer(det, tTopo_.get());
     const auto num = h_module_found->getValue(det);
     const auto denom = h_module_total->getValue(det);
     if (denom) {
       // use only the "good" modules
-      if (stripQuality_->getBadApvs(det) == 0) {
+      if (stripQuality_->getBadApvs(det) == 0 && calibData_.checkFedError(det)) {
         const auto eff = num / denom;
         hEffInLayer[layer]->Fill(eff);
         if (!autoIneffModTagging_) {
@@ -314,7 +328,7 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
         layerFound[layer] += num;
 
         // for the summary
-        //Have to do the decoding for which side to go on (ugh)
+        // Have to do the decoding for which side to go on (ugh)
         if (layer <= bounds::k_LayersAtTOBEnd) {
           goodlayerfound[layer] += num;
           goodlayertotal[layer] += denom;
@@ -374,7 +388,8 @@ void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStor
       hEffInLayer[i]->getTH1()->GetXaxis()->SetRange(1, hEffInLayer[i]->getNbinsX() + 1);
 
       for (auto det : stripDetIds_) {
-        if (stripQuality_->getBadApvs(det) == 0) {
+        // use only the "good" modules
+        if (stripQuality_->getBadApvs(det) == 0 && calibData_.checkFedError(det)) {
           const auto layer = ::checkLayer(det, tTopo_.get());
           if (layer == i) {
             const auto num = h_module_found->getValue(det);

CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyWorker.cc

@@ -66,8 +66,6 @@ class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
 private:
-  void beginJob();  // TODO remove
-  void endJob();    // TODO remove
   void bookHistograms(DQMStore::IBooker& booker, const edm::Run& run, const edm::EventSetup& setup) override;
   void analyze(const edm::Event& e, const edm::EventSetup& c) override;
   void fillForTraj(const TrajectoryAtInvalidHit& tm,
@@ -132,9 +130,6 @@ class SiStripHitEfficiencyWorker : public DQMEDAnalyzer {
   // output file
   std::set<uint32_t> badModules_;
 
-  // counters
-  int events, EventTrackCKF;
-
   struct EffME1 {
     EffME1() : hTotal(nullptr), hFound(nullptr) {}
     EffME1(MonitorElement* total, MonitorElement* found) : hTotal(total), hFound(found) {}
@@ -248,12 +243,6 @@ SiStripHitEfficiencyWorker::SiStripHitEfficiencyWorker(const edm::ParameterSet&
   }
 }
 
-void SiStripHitEfficiencyWorker::beginJob() {
-  // TODO convert to counters, or simply remove?
-  events = 0;
-  EventTrackCKF = 0;
-}
-
 void SiStripHitEfficiencyWorker::bookHistograms(DQMStore::IBooker& booker,
                                                 const edm::Run& run,
                                                 const edm::EventSetup& setup) {
@@ -382,7 +371,8 @@ void SiStripHitEfficiencyWorker::bookHistograms(DQMStore::IBooker& booker,
   // fill the FED Errors
   booker.setCurrentFolder(dqmDir_);
   const auto FEDErrorMapFolder = fmt::format("{}/FEDErrorTkDetMaps", dqmDir_);
-  calibData_.FEDErrorOccupancy = TkHistoMap(tkDetMap, booker, FEDErrorMapFolder, "perModule_FEDErrors", 0, false, true);
+  calibData_.FEDErrorOccupancy =
+      std::make_unique<TkHistoMap>(tkDetMap, booker, FEDErrorMapFolder, "perModule_FEDErrors", 0, false, true);
 }
 
 void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSetup& es) {
@@ -439,6 +429,10 @@ void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSe
 
   // fill the calibData with the FEDErrors
   for (const auto& fedErr : *fedErrorIds) {
+    // fill the TkHistoMap occupancy map
+    calibData_.FEDErrorOccupancy->fill(fedErr.rawId(), 1.);
+
+    // fill the unordered map
     if (calibData_.fedErrorCounts.find(fedErr.rawId()) != calibData_.fedErrorCounts.end()) {
       calibData_.fedErrorCounts[fedErr.rawId()] += 1;
     } else {
@@ -457,8 +451,6 @@ void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSe
   const auto& chi2Estimator = es.getData(chi2EstimatorToken_);
   const auto& prop = es.getData(propagatorToken_);
 
-  ++events;
-
   // Tracking
   LogDebug("SiStripHitEfficiencyWorker") << "number ckf tracks found = " << tracksCKF->size();
 
@@ -477,8 +469,6 @@ void SiStripHitEfficiencyWorker::analyze(const edm::Event& e, const edm::EventSe
       LogDebug("SiStripHitEfficiencyWorker")
           << "starting checking good event with < " << trackMultiplicityCut_ << " tracks";
 
-    ++EventTrackCKF;
-
     // actually should do a loop over all the tracks in the event here
 
     // Looping over traj-track associations to be able to get traj & track informations
@@ -1028,14 +1018,6 @@ void SiStripHitEfficiencyWorker::fillForTraj(const TrajectoryAtInvalidHit& tm,
                                          << ", layers=" << layers_;
 }
 
-void SiStripHitEfficiencyWorker::endJob() {
-  LogDebug("SiStripHitEfficiencyWorker") << " Events Analysed             " << events;
-  LogDebug("SiStripHitEfficiencyWorker") << " Number Of Tracked events    " << EventTrackCKF;
-
-  // fill the TkMap Data
-  calibData_.fillTkMapFromMap();
-}
-
 void SiStripHitEfficiencyWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
   desc.add<std::string>("dqmDir", "AlCaReco/SiStripHitEfficiency");

CalibTracker/SiStripHitEfficiency/src/TrajectoryAtInvalidHit.cc

@@ -1,18 +1,19 @@
 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
-#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
-#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
-#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
-#include "Geometry/CommonDetUnit/interface/GeomDet.h"
-#include "Geometry/CommonTopologies/interface/StripTopology.h"
-#include "Geometry/CommonTopologies/interface/PixelTopology.h"
 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "Geometry/CommonDetUnit/interface/GeomDet.h"
 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
-// #include "RecoTracker/MeasurementDet/interface/RecHitPropagator.h"
-#include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
+#include "Geometry/CommonTopologies/interface/PixelTopology.h"
+#include "Geometry/CommonTopologies/interface/StripTopology.h"
 #include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
+#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
+#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
+#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
+#include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
+// #include "RecoTracker/MeasurementDet/interface/RecHitPropagator.h"
 
 using namespace std;
 TrajectoryAtInvalidHit::TrajectoryAtInvalidHit(const TrajectoryMeasurement& tm,
@@ -67,7 +68,7 @@ TrajectoryAtInvalidHit::TrajectoryAtInvalidHit(const TrajectoryMeasurement& tm,
     theCombinedPredictedState = propagator.propagate(theCombinedPredictedState, surface);
 
     if (!theCombinedPredictedState.isValid()) {
-      cout << "found invalid combinedpredictedstate after propagation" << endl;
+      edm::LogWarning("TrajectoryAtInvalidHit") << "found invalid combinedpredictedstate after propagation" << endl;
       return;
     }