diff --git a/CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h b/CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h index 76bdb9bd1de95..b25afcd32803c 100644 --- 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 +// system includes #include +#include + +// 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 diff --git a/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc b/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc index 6f063234271c2..fd9ffc9eab6f3 100644 --- a/CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc +++ b/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& maps, const std::string& type) const; unsigned int countTotalHits(const std::vector& 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 void setEffBinLabels(const T gr, const T gr2, const unsigned int nLayers) const; void makeSummaryVsVariable(DQMStore::IGetter& getter, DQMStore::IBooker& booker, ::projections theProj) const; @@ -559,7 +559,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 +659,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 +684,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 +903,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 }