Merge pull request #41126 from mmusich/mm_devel_130X_modernizeCalibTr…

…ackerSiStripHitEfficiency [13_0_X] Combined backport of developments in `CalibTracker/SiStripHitEfficiency`
cms-sw · Mar 22, 2023 · 23b415a · 23b415a
2 parents 619a0dc + 8307e38
commit 23b415a
Show file tree

Hide file tree

Showing 6 changed files with 423 additions and 283 deletions.
diff --git a/CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h b/CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h
@@ -4,13 +4,22 @@
 // A bunch of helper functions to deal with menial tasks in the
 // hit efficiency computation for the PCL workflow
 
-#include "TString.h"
-#include <string>
+// system includes
 #include <fmt/printf.h>
+#include <string>
+
+// user includes
 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
+#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
 
+// ROOT includes
+#include "TEfficiency.h"
+#include "TProfile.h"
+#include "TString.h"
+
 namespace {
 
   enum bounds {
@@ -23,6 +32,20 @@ namespace {
     k_END_OF_LAYS_AND_RINGS = 35
   };
 
+  /*
+   * for the trend plots of efficiency vs some variable
+   */
+  enum projections { k_vs_LUMI = 0, k_vs_PU = 1, k_vs_BX = 2, k_SIZE = 3 };
+
+  const std::array<std::string, projections::k_SIZE> projFolder = {{"VsLumi", "VsPu", "VsBx"}};
+  const std::array<std::string, projections::k_SIZE> projFoundHisto = {
+      {"layerfound_vsLumi_layer_", "layerfound_vsPU_layer_", "foundVsBx_layer"}};
+  const std::array<std::string, projections::k_SIZE> projTotalHisto = {
+      {"layertotal_vsLumi_layer_", "layertotal_vsPU_layer_", "totalVsBx_layer"}};
+  const std::array<std::string, projections::k_SIZE> projTitle = {{"Inst Lumi", "Pile-Up", "Bunch Crossing"}};
+  const std::array<std::string, projections::k_SIZE> projXtitle = {
+      {"instantaneous luminosity [Hz/cm^{2}]", "Pile-Up events", "Bunch Crossing number"}};
+
   inline void replaceInString(std::string& str, const std::string& from, const std::string& to) {
     if (from.empty())
       return;
@@ -201,5 +224,50 @@ namespace {
     return phi;
   }
 
+  inline TProfile* computeEff(const TH1F* num, const TH1F* denum, const std::string nameHist) {
+    std::string name = "eff_" + nameHist;
+    std::string title = "SiStrip Hit Efficiency" + std::string(num->GetTitle());
+    TProfile* efficHist = new TProfile(name.c_str(),
+                                       title.c_str(),
+                                       denum->GetXaxis()->GetNbins(),
+                                       denum->GetXaxis()->GetXmin(),
+                                       denum->GetXaxis()->GetXmax());
+
+    for (int i = 1; i <= denum->GetNbinsX(); i++) {
+      double nNum = num->GetBinContent(i);
+      double nDenum = denum->GetBinContent(i);
+      if (nDenum == 0 || nNum == 0) {
+        continue;
+      }
+      if (nNum > nDenum) {
+        edm::LogWarning("SiStripHitEfficiencyHelpers")
+            << "Alert! specific bin's num is bigger than denum " << i << " " << nNum << " " << nDenum;
+        nNum = nDenum;  // set the efficiency to 1
+      }
+      const double effVal = nNum / nDenum;
+      efficHist->SetBinContent(i, effVal);
+      efficHist->SetBinEntries(i, 1);
+      const double errLo = TEfficiency::ClopperPearson((int)nDenum, (int)nNum, 0.683, false);
+      const double errUp = TEfficiency::ClopperPearson((int)nDenum, (int)nNum, 0.683, true);
+      const double errVal = (effVal - errLo > errUp - effVal) ? effVal - errLo : errUp - effVal;
+      efficHist->SetBinError(i, sqrt(effVal * effVal + errVal * errVal));
+
+      LogDebug("SiStripHitEfficiencyHelpers") << __PRETTY_FUNCTION__ << " " << nameHist << " bin:" << i
+                                              << " err:" << sqrt(effVal * effVal + errVal * errVal);
+    }
+    return efficHist;
+  }
+
+  inline Measurement1D computeCPEfficiency(const double num, const double den) {
+    if (den > 0) {
+      const double effVal = num / den;
+      const double errLo = TEfficiency::ClopperPearson((int)den, (int)num, 0.683, false);
+      const double errUp = TEfficiency::ClopperPearson((int)den, (int)num, 0.683, true);
+      const double errVal = (effVal - errLo > errUp - effVal) ? effVal - errLo : errUp - effVal;
+      return Measurement1D(effVal, errVal);
+    } else {
+      return Measurement1D(0., 0.);
+    }
+  }
 }  // namespace
 #endif
diff --git a/CalibTracker/SiStripHitEfficiency/plugins/HitEff.cc b/CalibTracker/SiStripHitEfficiency/plugins/HitEff.cc
@@ -92,6 +92,7 @@ HitEff::HitEff(const edm::ParameterSet& conf)
       chi2MeasurementEstimatorToken_(esConsumes(edm::ESInputTag("", "Chi2"))),
       propagatorToken_(esConsumes(edm::ESInputTag("", "PropagatorWithMaterial"))),
       conf_(conf) {
+  usesResource(TFileService::kSharedResource);
   compSettings = conf_.getUntrackedParameter<int>("CompressionSettings", -1);
   layers = conf_.getParameter<int>("Layer");
   DEBUG = conf_.getParameter<bool>("Debug");
@@ -243,7 +244,7 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
   //e.getByLabel("siStripDigis", fedErrorIds );
   e.getByToken(digis_token_, fedErrorIds);
 
-  ESHandle<MeasurementTracker> measurementTrackerHandle = es.getHandle(measurementTkToken_);
+  edm::ESHandle<MeasurementTracker> measurementTrackerHandle = es.getHandle(measurementTkToken_);
 
   edm::Handle<MeasurementTrackerEvent> measurementTrackerEvent;
   //e.getByLabel("MeasurementTrackerEvent", measurementTrackerEvent);
@@ -298,7 +299,7 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
 
 #ifdef ExtendedCALIBTree
     //get dEdx info if available
-    Handle<ValueMap<DeDxData> > dEdxUncalibHandle;
+    edm::Handle<ValueMap<DeDxData> > dEdxUncalibHandle;
     if (e.getByLabel("dedxMedianCTF", dEdxUncalibHandle)) {
       const ValueMap<DeDxData> dEdxTrackUncalib = *dEdxUncalibHandle.product();
 
@@ -311,7 +312,7 @@ void HitEff::analyze(const edm::Event& e, const edm::EventSetup& es) {
     }
 
     //get muon and ecal timing info if available
-    Handle<MuonCollection> muH;
+    edm::Handle<MuonCollection> muH;
     if (e.getByLabel("muonsWitht0Correction", muH)) {
       const MuonCollection& muonsT0 = *muH.product();
       if (!muonsT0.empty()) {

diff --git a/CalibTracker/SiStripHitEfficiency/plugins/HitEff.h b/CalibTracker/SiStripHitEfficiency/plugins/HitEff.h
@@ -49,7 +49,7 @@
 
 class TrackerTopology;
 
-class HitEff : public edm::one::EDAnalyzer<> {
+class HitEff : public edm::one::EDAnalyzer<edm::one::SharedResources> {
 public:
   explicit HitEff(const edm::ParameterSet& conf);
   ~HitEff() override = default;