More refined error treatment in trend of efficiency vs lumi, PU and BX in SiStripHitEfficiencyHarvester

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 {
@@ -215,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

CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEffFromCalibTree.cc

@@ -120,7 +120,7 @@ class SiStripHitEffFromCalibTree : public ConditionDBWriter<SiStripBadStrip> {
   TString getLayerSideName(Long_t k);
 
   // to be used everywhere
-  static constexpr int SiStripLayers_ = 22;
+  static constexpr int siStripLayers_ = 22;
   static constexpr double nBxInAnOrbit_ = 3565;
 
   edm::Service<TFileService> fs;
@@ -460,7 +460,7 @@ void SiStripHitEffFromCalibTree::algoAnalyze(const edm::Event& e, const edm::Eve
         badquality = true;
 
       // don't compute efficiencies in modules from TOB6 and TEC9
-      if (!showTOB6TEC9_ && (layer_wheel == 10 || layer_wheel == SiStripLayers_))
+      if (!showTOB6TEC9_ && (layer_wheel == 10 || layer_wheel == siStripLayers_))
         continue;
 
       // don't use bad modules given in the bad module list
@@ -616,7 +616,7 @@ void SiStripHitEffFromCalibTree::algoAnalyze(const edm::Event& e, const edm::Eve
               goodlayerfound[layer + 3]++;
             goodlayertotal[layer + 3]++;
           }
-        } else if (layer > 13 && layer <= SiStripLayers_) {
+        } else if (layer > 13 && layer <= siStripLayers_) {
           if (((id >> 18) & 0x3) == 1) {
             if (!badflag)
               goodlayerfound[layer + 3]++;
@@ -643,7 +643,7 @@ void SiStripHitEffFromCalibTree::algoAnalyze(const edm::Event& e, const edm::Eve
             alllayerfound[layer + 3]++;
           alllayertotal[layer + 3]++;
         }
-      } else if (layer > 13 && layer <= SiStripLayers_) {
+      } else if (layer > 13 && layer <= siStripLayers_) {
         if (((id >> 18) & 0x3) == 1) {
           if (!badflag)
             alllayerfound[layer + 3]++;
@@ -885,7 +885,7 @@ void SiStripHitEffFromCalibTree::makeHotColdMaps() {
   //Already have access to the data as a private variable
   //Create all of the histograms in the TFileService
   TH2F* temph2;
-  for (Long_t maplayer = 1; maplayer <= SiStripLayers_; maplayer++) {
+  for (Long_t maplayer = 1; maplayer <= siStripLayers_; maplayer++) {
     //Initialize all of the histograms
     if (maplayer > 0 && maplayer <= 4) {
       //We are in the TIB
@@ -959,7 +959,7 @@ void SiStripHitEffFromCalibTree::makeHotColdMaps() {
       HotColdMaps.push_back(temph2);
     }
   }
-  for (Long_t mylayer = 1; mylayer <= SiStripLayers_; mylayer++) {
+  for (Long_t mylayer = 1; mylayer <= siStripLayers_; mylayer++) {
     //Determine what kind of plot we want to write out
     //Loop through the entirety of each layer
     //Create an array of the histograms
@@ -1013,7 +1013,7 @@ void SiStripHitEffFromCalibTree::makeTKMap(bool autoTagging = false) {
   double myeff, mynum, myden, myeff_up;
   double layer_min_eff = 0;
 
-  for (Long_t i = 1; i <= SiStripLayers_; i++) {
+  for (Long_t i = 1; i <= siStripLayers_; i++) {
     //Loop over every layer, extracting the information from
     //the map of the efficiencies
     layertotal[i] = 0;
@@ -1152,7 +1152,7 @@ void SiStripHitEffFromCalibTree::totalStatistics() {
     subdettotal[i] = 0;
   }
 
-  for (Long_t i = 1; i <= SiStripLayers_; i++) {
+  for (Long_t i = 1; i <= siStripLayers_; i++) {
     layereff = double(layerfound[i]) / double(layertotal[i]);
     LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers) << ") has total efficiency " << layereff
              << " " << layerfound[i] << "/" << layertotal[i];
@@ -1195,7 +1195,7 @@ void SiStripHitEffFromCalibTree::makeSummary() {
     nLayers = 30;
   if (!showEndcapSides_) {
     if (!showRings_)
-      nLayers = SiStripLayers_;
+      nLayers = siStripLayers_;
     else
       nLayers = 20;
   }
@@ -1374,7 +1374,7 @@ void SiStripHitEffFromCalibTree::makeSummary() {
 void SiStripHitEffFromCalibTree::makeSummaryVsBx() {
   LOGPRINT << "Computing efficiency vs bx";
 
-  unsigned int nLayers = SiStripLayers_;
+  unsigned int nLayers = siStripLayers_;
   if (showRings_)
     nLayers = 20;
 
@@ -1453,7 +1453,7 @@ void SiStripHitEffFromCalibTree::makeSummaryVsBx() {
 }
 
 void SiStripHitEffFromCalibTree::computeEff(vector<TH1F*>& vhfound, vector<TH1F*>& vhtotal, string name) {
-  unsigned int nLayers = SiStripLayers_;
+  unsigned int nLayers = siStripLayers_;
   if (showRings_)
     nLayers = 20;
 
@@ -1499,7 +1499,7 @@ void SiStripHitEffFromCalibTree::makeSummaryVsLumi() {
 
   else {  // from infos per hit
 
-    unsigned int nLayers = SiStripLayers_;
+    unsigned int nLayers = siStripLayers_;
     if (showRings_)
       nLayers = 20;
     unsigned int nLayersForAvg = 0;

CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc

@@ -90,7 +90,7 @@ class SiStripHitEfficiencyHarvester : public DQMEDHarvester {
   void printAndWriteBadModules(const SiStripQuality& quality, const SiStripDetInfo& detInfo) const;
   bool checkMapsValidity(const std::vector<MonitorElement*>& maps, const std::string& type) const;
   unsigned int countTotalHits(const std::vector<MonitorElement*>& maps); /* to check if TK was ON */
-  void makeSummary(DQMStore::IGetter& getter, DQMStore::IBooker& booker) const;
+  void makeSummary(DQMStore::IGetter& getter, DQMStore::IBooker& booker, bool doProfiles = false) const;
   template <typename T>
   void setEffBinLabels(const T gr, const T gr2, const unsigned int nLayers) const;
   void makeSummaryVsVariable(DQMStore::IGetter& getter, DQMStore::IBooker& booker, ::projections theProj) const;
@@ -500,13 +500,8 @@ void SiStripHitEfficiencyHarvester::printTotalStatistics(
   int totalfound = 0;
   int totaltotal = 0;
   double layereff;
-  int subdetfound[5];
-  int subdettotal[5];
-
-  for (unsigned int i = 1; i < 5; i++) {
-    subdetfound[i] = 0;
-    subdettotal[i] = 0;
-  }
+  int subdetfound[5] = {0, 0, 0, 0, 0};
+  int subdettotal[5] = {0, 0, 0, 0, 0};
 
   for (unsigned int i = 1; i <= bounds::k_LayersAtTECEnd; i++) {
     layereff = double(layerFound[i]) / double(layerTotal[i]);
@@ -559,7 +554,9 @@ void SiStripHitEfficiencyHarvester::writeBadStripPayload(const SiStripQuality& q
   }
 }
 
-void SiStripHitEfficiencyHarvester::makeSummary(DQMStore::IGetter& getter, DQMStore::IBooker& booker) const {
+void SiStripHitEfficiencyHarvester::makeSummary(DQMStore::IGetter& getter,
+                                                DQMStore::IBooker& booker,
+                                                bool doProfiles) const {
   // use goodlayer_total/found and alllayer_total/found, collapse side and/or ring if needed
   unsigned int nLayers{34};  // default
   if (showRings_)
@@ -657,6 +654,23 @@ void SiStripHitEfficiencyHarvester::makeSummary(DQMStore::IGetter& getter, DQMSt
   MonitorElement* h_eff_good =
       booker.book1D("eff_good", "Strip hit efficiency for good modules", nLayers + 1, 0, nLayers + 1);
 
+  if (doProfiles) {
+    // now do the profile
+    TProfile* profile_all = ::computeEff(found2->getTH1F(), all2->getTH1F(), "all");
+    profile_all->SetMinimum(tkMapMin_);
+    profile_all->SetTitle("Strip hit efficiency for all modules");
+    booker.bookProfile(profile_all->GetName(), profile_all);
+
+    TProfile* profile_good = ::computeEff(found->getTH1F(), all->getTH1F(), "good");
+    profile_good->SetMinimum(tkMapMin_);
+    profile_good->SetTitle("Strip hit efficiency for good modules");
+    booker.bookProfile(profile_good->GetName(), profile_good);
+
+    // clean the house
+    delete profile_all;
+    delete profile_good;
+  }
+
   for (int i = 1; i < found->getNbinsX(); i++) {
     const auto& den_all = all2->getBinContent(i);
     const auto& num_all = found2->getBinContent(i);
@@ -665,18 +679,26 @@ void SiStripHitEfficiencyHarvester::makeSummary(DQMStore::IGetter& getter, DQMSt
 
     // fill all modules efficiency
     if (den_all > 0.) {
-      float eff_all = num_all / den_all;
-      float err_eff_all = (eff_all * (1 - eff_all)) / den_all;
-      h_eff_all->setBinContent(i, eff_all);
-      h_eff_all->setBinError(i, err_eff_all);
+      // naive binomial errors
+      //float eff_all = num_all / den_all;
+      //float err_eff_all = (eff_all * (1 - eff_all)) / den_all;
+
+      // use Clopper-Pearson errors
+      const auto& effPair_all = ::computeCPEfficiency(num_all, den_all);
+      h_eff_all->setBinContent(i, effPair_all.value());
+      h_eff_all->setBinError(i, effPair_all.error());
     }
 
     // fill good modules efficiency
     if (den_good > 0.) {
-      float eff_good = num_good / den_good;
-      float err_eff_good = (eff_good * (1 - eff_good)) / den_good;
-      h_eff_good->setBinContent(i, eff_good);
-      h_eff_good->setBinError(i, err_eff_good);
+      // naive binomial errors
+      //float eff_good = num_good / den_good;
+      //float err_eff_good = (eff_good * (1 - eff_good)) / den_good;
+
+      // use Clopper-Pearson errors
+      const auto& effPair_good = ::computeCPEfficiency(num_good, den_good);
+      h_eff_good->setBinContent(i, effPair_good.value());
+      h_eff_good->setBinError(i, effPair_good.error());
     }
   }
 
@@ -876,18 +898,29 @@ void SiStripHitEfficiencyHarvester::makeSummaryVsVariable(DQMStore::IGetter& get
       const auto& num = hfound->getBinContent(i);
 
       // fill all modules efficiency
-      // error on efficiency computed with binomial approximation
       if (den > 0.) {
-        float eff = num / den;
-        float err_eff = (eff * (1 - eff)) / den;
-        effVsVariable[iLayer]->setBinContent(i, eff);
-        effVsVariable[iLayer]->setBinError(i, err_eff);
+        const auto& effPair = ::computeCPEfficiency(num, den);
+        effVsVariable[iLayer]->setBinContent(i, effPair.value());
+        effVsVariable[iLayer]->setBinError(i, effPair.error());
+
+        LogDebug("SiStripHitEfficiencyHarvester")
+            << __PRETTY_FUNCTION__ << " " << lyrName << " bin:" << i << " err:" << effPair.error() << std::endl;
       }
     }
 
     // graphics adjustment
     effVsVariable[iLayer]->getTH1F()->SetMinimum(tkMapMin_);
 
+    // now do the profile
+    TProfile* profile = ::computeEff(hfound->getTH1F(), htotal->getTH1F(), lyrName);
+    TString title =
+        fmt::sprintf("Efficiency vs %s for layer %s;%s;SiStrip Hit efficiency", titleString, lyrName, titleXString);
+    profile->SetMinimum(tkMapMin_);
+
+    profile->SetTitle(title.Data());
+    booker.bookProfile(profile->GetName(), profile);
+
+    delete profile;
   }  // loop on layers
 }