Skip to content

Commit

Permalink
add provision to create TProfiles of efficiency and use Clopper-Pears…
Browse files Browse the repository at this point in the history 
…on errors in trend plots
  • Loading branch information
@mmusich
mmusich committed Mar 22, 2023
1 parent b227118 commit 7046634
Show file tree
Hide file tree
Showing 2 changed files with 109 additions and 17 deletions.
58 changes: 56 additions & 2 deletions CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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
68 changes: 53 additions & 15 deletions CalibTracker/SiStripHitEfficiency/plugins/SiStripHitEfficiencyHarvester.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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_)
Expand Down Expand Up @@ -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);
Expand All @@ -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());
}
}

Expand Down Expand Up @@ -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
}

Expand Down

0 comments on commit 7046634

Please sign in to comment.