28 #include "TEfficiency.h"
31 #define LOGPRINT edm::LogPrint("SiStripHitEfficiencyHarvester")
56 std::unique_ptr<TrackerTopology>
tTopo_;
62 void printTotalStatistics(
const std::array<long, 23>& layerFound,
const std::array<long, 23>& layerTotal)
const;
72 : isAtPCL_(conf.getParameter<bool>(
"isAtPCL")),
73 showRings_(conf.getUntrackedParameter<bool>(
"ShowRings",
false)),
74 autoIneffModTagging_(conf.getUntrackedParameter<bool>(
"AutoIneffModTagging",
false)),
75 doStoreOnDB_(conf.getParameter<bool>(
"doStoreOnDB")),
76 nTEClayers_(showRings_ ? 7 : 9),
77 threshold_(conf.getParameter<double>(
"Threshold")),
78 nModsMin_(conf.getParameter<int>(
"nModsMin")),
79 tkMapMin_(conf.getUntrackedParameter<double>(
"TkMapMin", 0.9)),
80 title_(conf.getParameter<std::
string>(
"Title")),
81 record_(conf.getParameter<std::
string>(
"Record")),
99 for (
const auto& det : tkGeom.detUnits()) {
100 if (dynamic_cast<const StripGeomDetUnit*>(det)) {
109 std::vector<bool> isAvailable;
110 isAvailable.reserve(maps.size());
112 maps.begin() + 1, maps.end(), std::back_inserter(isAvailable), [](
auto&
x) {
return !(
x ==
nullptr); });
115 for (
const auto& it : isAvailable) {
117 LogDebug(
"SiStripHitEfficiencyHarvester") <<
" layer: " <<
count <<
" " << it << std::endl;
119 LogDebug(
"SiStripHitEfficiencyHarvester") <<
"resolving to " << maps[
count]->getName() << std::endl;
123 bool areMapsAvailable{
true};
125 for (
const auto& it : isAvailable) {
129 << type <<
" TkHistoMap for layer " << layerCount <<
" was not found.\n -> Aborting!";
130 areMapsAvailable =
false;
134 return areMapsAvailable;
141 LOGPRINT <<
"A module is bad if the upper limit on the efficiency is < to the avg in the layer - " <<
threshold_
142 <<
" and has at least " << nModsMin_ <<
" nModsMin.";
144 auto h_module_total = std::make_unique<TkHistoMap>(
tkDetMap_.get());
145 h_module_total->loadTkHistoMap(
"AlCaReco/SiStripHitEfficiency",
"perModule_total");
146 auto h_module_found = std::make_unique<TkHistoMap>(
tkDetMap_.get());
147 h_module_found->loadTkHistoMap(
"AlCaReco/SiStripHitEfficiency",
"perModule_found");
150 const auto& totalMaps = h_module_total->getAllMaps();
151 const auto& foundMaps = h_module_found->getAllMaps();
153 LogDebug(
"SiStripHitEfficiencyHarvester")
154 <<
"totalMaps.size(): " << totalMaps.size() <<
" foundMaps.size() " << foundMaps.size() << std::endl;
160 LogDebug(
"SiStripHitEfficiencyHarvester")
161 <<
"isTotalMapAvailable: " << isTotalMapAvailable <<
" isFoundMapAvailable " << isFoundMapAvailable << std::endl;
164 if (!isTotalMapAvailable
or !isFoundMapAvailable)
167 LogDebug(
"SiStripHitEfficiencyHarvester")
168 <<
"Entries in total TkHistoMap for layer 3: " << h_module_total->getMap(3)->getEntries() <<
", found "
169 << h_module_found->getMap(3)->getEntries();
171 std::vector<MonitorElement*> hEffInLayer(std::size_t(1),
nullptr);
172 hEffInLayer.reserve(23);
173 for (std::size_t
i = 1;
i != 23; ++
i) {
174 hEffInLayer.push_back(
175 booker.
book1D(Form(
"eff_layer%i",
int(
i)), Form(
"Module efficiency in layer %i",
int(
i)), 201, 0, 1.005));
177 std::array<long, 23> layerTotal{};
178 std::array<long, 23> layerFound{};
190 TrackerMap tkMapDen{
" Detector denominator "};
191 std::map<unsigned int, double> badModules;
195 const auto num = h_module_found->getValue(det);
196 const auto denom = h_module_total->getValue(det);
198 const auto eff = num /
denom;
199 hEffInLayer[
layer]->Fill(eff);
203 badModules[det] = eff;
204 tkMapBad.fillc(det, 255, 0, 0);
206 << det.rawId() <<
" efficiency: " << eff <<
" , " << num <<
"/" <<
denom;
209 tkMapBad.fillc(det, 255, 255, 255);
213 << det.rawId() <<
" efficiency: " << eff <<
" , " << num <<
"/" <<
denom;
215 if (denom < nModsMin_) {
217 << det.rawId() <<
" is under occupancy at " <<
denom;
221 tkMap.fill(det, 1. - eff);
222 tkMapEff.fill(det, eff);
223 tkMapNum.fill(det, num);
224 tkMapDen.fill(det,
denom);
232 for (Long_t
i = 1;
i <= 22;
i++) {
234 hEffInLayer[
i]->getTH1()->GetXaxis()->SetRange(
235 3, hEffInLayer[
i]->getNbinsX() + 1);
236 const double layer_min_eff = hEffInLayer[
i]->getMean() -
std::max(2.5 * hEffInLayer[
i]->getRMS(),
threshold_);
237 LOGPRINT <<
"Layer " <<
i <<
" threshold for bad modules: <" << layer_min_eff
238 <<
" (layer mean: " << hEffInLayer[
i]->getMean() <<
" rms: " << hEffInLayer[
i]->getRMS() <<
")";
240 hEffInLayer[
i]->getTH1()->GetXaxis()->SetRange(1, hEffInLayer[
i]->getNbinsX() + 1);
242 for (
auto det : stripDetIds_) {
243 const auto layer = ::checkLayer(det,
tTopo_.get());
245 const auto num = h_module_found->getValue(det);
246 const auto denom = h_module_total->getValue(det);
249 const auto eff_up = TEfficiency::Bayesian(
denom,
num, .99, 1, 1,
true);
251 if ((
denom >= nModsMin_) && (eff_up < layer_min_eff)) {
253 badModules[det] = eff;
254 tkMapBad.fillc(det, 255, 0, 0);
257 tkMapBad.fillc(det, 255, 255, 255);
259 if (eff_up < layer_min_eff + 0.08)
262 << det.rawId() <<
" efficiency: " << eff <<
" , " <<
num <<
"/" <<
denom
263 <<
" , upper limit: " << eff_up;
264 if (
denom < nModsMin_) {
266 << det.rawId() <<
" layer " <<
layer <<
" is under occupancy at " <<
denom;
274 tkMap.save(
true, 0, 0,
"SiStripHitEffTKMap_NEW.png");
275 tkMapBad.save(
true, 0, 0,
"SiStripHitEffTKMapBad_NEW.png");
276 tkMapEff.save(
true,
tkMapMin_, 1.,
"SiStripHitEffTKMapEff_NEW.png");
277 tkMapNum.save(
true, 0, 0,
"SiStripHitEffTKMapNum_NEW.png");
278 tkMapDen.save(
true, 0, 0,
"SiStripHitEffTKMapDen_NEW.png");
286 for (
const auto it : badModules) {
287 const auto det = it.first;
288 std::vector<unsigned int> badStripList;
293 badStripList.push_back(pQuality.encode(0, nStrips, 0));
295 LOGPRINT <<
"ID1 shoudl match list of modules above " << det;
296 pQuality.compact(det, badStripList);
299 pQuality.fillBadComponents();
303 edm::LogInfo(
"SiStripHitEfficiencyHarvester") <<
"Will not produce payload!";
322 const std::array<long, 23>& layerTotal)
const {
330 for (Long_t
i = 1;
i < 5;
i++) {
335 for (Long_t
i = 1;
i <= 22;
i++) {
336 layereff = double(layerFound[
i]) / double(layerTotal[i]);
338 << layereff <<
" " << layerFound[
i] <<
"/" << layerTotal[
i];
339 totalfound += layerFound[
i];
340 totaltotal += layerTotal[
i];
342 subdetfound[1] += layerFound[
i];
343 subdettotal[1] += layerTotal[
i];
345 if (i >= 5 && i < 11) {
346 subdetfound[2] += layerFound[
i];
347 subdettotal[2] += layerTotal[
i];
349 if (i >= 11 && i < 14) {
350 subdetfound[3] += layerFound[
i];
351 subdettotal[3] += layerTotal[
i];
354 subdetfound[4] += layerFound[
i];
355 subdettotal[4] += layerTotal[
i];
359 LOGPRINT <<
"The total efficiency is " << double(totalfound) / double(totaltotal);
360 LOGPRINT <<
" TIB: " << double(subdetfound[1]) / subdettotal[1] <<
" " << subdetfound[1] <<
"/"
362 LOGPRINT <<
" TOB: " << double(subdetfound[2]) / subdettotal[2] <<
" " << subdetfound[2] <<
"/"
364 LOGPRINT <<
" TID: " << double(subdetfound[3]) / subdettotal[3] <<
" " << subdetfound[3] <<
"/"
366 LOGPRINT <<
" TEC: " << double(subdetfound[4]) / subdettotal[4] <<
" " << subdetfound[4] <<
"/"
375 if (!pBadStrip.put(rIt->detid,
range))
376 edm::LogError(
"SiStripHitEfficiencyHarvester") <<
"detid already exists in SiStripBadStrip";
395 for (
unsigned int iLayer = 1; iLayer != (
showRings_ ? 20 : 22); ++iLayer) {
397 getter.
get(
fmt::format(
"AlCaReco/SiStripHitEfficiency/layerfound_vsLumi_layer_{}", iLayer))->
getTH1F();
399 getter.
get(
fmt::format(
"AlCaReco/SiStripHitEfficiency/layertotal_vsLumi_layer_{}", iLayer))->
getTH1F();
401 if (hfound ==
nullptr or htotal ==
nullptr) {
402 if (hfound ==
nullptr)
404 <<
fmt::format(
"AlCaReco/SiStripHitEfficiency/layerfound_vsLumi_layer_{}", iLayer) <<
" was not found!";
405 if (htotal ==
nullptr)
407 <<
fmt::format(
"AlCaReco/SiStripHitEfficiency/layertotal_vsLumi_layer_{}", iLayer) <<
" was not found!";
412 if (!hfound->GetSumw2())
414 if (!htotal->GetSumw2())
416 for (Long_t
i = 0;
i != hfound->GetNbinsX() + 1; ++
i) {
417 if (hfound->GetBinContent(
i) == 0)
418 hfound->SetBinContent(
i, 1
e-6);
419 if (htotal->GetBinContent(
i) == 0)
420 htotal->SetBinContent(
i, 1);
422 LogDebug(
"SiStripHitEfficiencyHarvester")
423 <<
"Total hits for layer " << iLayer <<
" (vs lumi): " << htotal->GetEntries() <<
", found "
424 << hfound->GetEntries();
435 std::stringstream ssV[4][19],
444 ssV[
i][
comp] <<
" \t " << ((bc.
BadApvs) & 0x1) <<
" " << ((bc.BadApvs >> 1) & 0x1) <<
" ";
446 ssV[
i][
comp] <<
"x x " << ((bc.BadApvs >> 2) & 0x1) <<
" " << ((bc.BadApvs >> 3) & 0x1);
448 ssV[
i][
comp] << ((bc.BadApvs >> 2) & 0x1) <<
" " << ((bc.BadApvs >> 3) & 0x1) <<
" " << ((bc.BadApvs >> 4) & 0x1)
449 <<
" " << ((bc.BadApvs >> 5) & 0x1) <<
" ";
452 nBad[
i][0][2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
453 ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
454 nBad[
i][
comp][2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
455 ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
458 nBad[
i][0][1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
459 nBad[
i][
comp][1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
474 int nBadComp[4][19][4];
478 std::stringstream ssV[4][19];
480 for (
int i = 0;
i < 4; ++
i) {
482 for (
int j = 0;
j < 19; ++
j) {
484 for (
int k = 0;
k < 4; ++
k)
485 nBadComp[
i][
j][
k] = 0;
494 nTkBadComp[1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
496 nTkBadComp[2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
497 ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
532 DetId det{rp->detid};
534 int component = -999;
535 switch (det.subdetId()) {
538 component =
tTopo_->tibLayer(det);
542 component =
tTopo_->tidSide(det) == 2 ?
tTopo_->tidWheel(det) :
tTopo_->tidWheel(det) + 3;
546 component =
tTopo_->tobLayer(det);
550 component =
tTopo_->tecSide(det) == 2 ?
tTopo_->tecWheel(det) :
tTopo_->tecWheel(det) + 9;
558 float percentage = 0;
559 for (
int it = 0; it < sqrange.second - sqrange.first; it++) {
561 nTkBadComp[3] +=
range;
562 nBadComp[subdet][0][3] +=
range;
563 nBadComp[subdet][component][3] +=
range;
569 edm::LogError(
"SiStripHitEfficiencyHarvester") <<
"PROBLEM detid " << det.rawId() <<
" value " << percentage;
573 std::ostringstream
ss;
574 ss <<
"\n-----------------\nGlobal Info\n-----------------";
575 ss <<
"\nBadComp \t Modules \tFibers "
576 "\tApvs\tStrips\n----------------------------------------------------------------";
577 ss <<
"\nTracker:\t\t" << nTkBadComp[0] <<
"\t" << nTkBadComp[1] <<
"\t" << nTkBadComp[2] <<
"\t" << nTkBadComp[3];
578 ss <<
"\nTIB:\t\t\t" << nBadComp[0][0][0] <<
"\t" << nBadComp[0][0][1] <<
"\t" << nBadComp[0][0][2] <<
"\t"
579 << nBadComp[0][0][3];
580 ss <<
"\nTID:\t\t\t" << nBadComp[1][0][0] <<
"\t" << nBadComp[1][0][1] <<
"\t" << nBadComp[1][0][2] <<
"\t"
581 << nBadComp[1][0][3];
582 ss <<
"\nTOB:\t\t\t" << nBadComp[2][0][0] <<
"\t" << nBadComp[2][0][1] <<
"\t" << nBadComp[2][0][2] <<
"\t"
583 << nBadComp[2][0][3];
584 ss <<
"\nTEC:\t\t\t" << nBadComp[3][0][0] <<
"\t" << nBadComp[3][0][1] <<
"\t" << nBadComp[3][0][2] <<
"\t"
585 << nBadComp[3][0][3];
588 for (
int i = 1;
i < 5; ++
i)
589 ss <<
"\nTIB Layer " <<
i <<
" :\t\t" << nBadComp[0][
i][0] <<
"\t" << nBadComp[0][
i][1] <<
"\t" << nBadComp[0][
i][2]
590 <<
"\t" << nBadComp[0][
i][3];
592 for (
int i = 1;
i < 4; ++
i)
593 ss <<
"\nTID+ Disk " <<
i <<
" :\t\t" << nBadComp[1][
i][0] <<
"\t" << nBadComp[1][
i][1] <<
"\t" << nBadComp[1][
i][2]
594 <<
"\t" << nBadComp[1][
i][3];
595 for (
int i = 4;
i < 7; ++
i)
596 ss <<
"\nTID- Disk " <<
i - 3 <<
" :\t\t" << nBadComp[1][
i][0] <<
"\t" << nBadComp[1][
i][1] <<
"\t"
597 << nBadComp[1][
i][2] <<
"\t" << nBadComp[1][
i][3];
599 for (
int i = 1;
i < 7; ++
i)
600 ss <<
"\nTOB Layer " <<
i <<
" :\t\t" << nBadComp[2][
i][0] <<
"\t" << nBadComp[2][
i][1] <<
"\t" << nBadComp[2][
i][2]
601 <<
"\t" << nBadComp[2][
i][3];
603 for (
int i = 1;
i < 10; ++
i)
604 ss <<
"\nTEC+ Disk " <<
i <<
" :\t\t" << nBadComp[3][
i][0] <<
"\t" << nBadComp[3][
i][1] <<
"\t" << nBadComp[3][
i][2]
605 <<
"\t" << nBadComp[3][
i][3];
606 for (
int i = 10;
i < 19; ++
i)
607 ss <<
"\nTEC- Disk " <<
i - 9 <<
" :\t\t" << nBadComp[3][
i][0] <<
"\t" << nBadComp[3][
i][1] <<
"\t"
608 << nBadComp[3][
i][2] <<
"\t" << nBadComp[3][
i][3];
611 ss <<
"\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
612 "Apvs\n----------------------------------------------------------------";
613 for (
int i = 1;
i < 5; ++
i)
614 ss <<
"\nTIB Layer " <<
i <<
" :" << ssV[0][
i].
str();
616 for (
int i = 1;
i < 4; ++
i)
617 ss <<
"\nTID+ Disk " <<
i <<
" :" << ssV[1][
i].
str();
618 for (
int i = 4;
i < 7; ++
i)
619 ss <<
"\nTID- Disk " <<
i - 3 <<
" :" << ssV[1][
i].
str();
621 for (
int i = 1;
i < 7; ++
i)
622 ss <<
"\nTOB Layer " <<
i <<
" :" << ssV[2][
i].
str();
624 for (
int i = 1;
i < 10; ++
i)
625 ss <<
"\nTEC+ Disk " <<
i <<
" :" << ssV[3][
i].
str();
626 for (
int i = 10;
i < 19; ++
i)
627 ss <<
"\nTEC- Disk " <<
i - 9 <<
" :" << ssV[3][
i].
str();
632 std::ofstream badModules;
633 badModules.open(
"BadModules_NEW.log");
634 badModules <<
"\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
635 "Apvs\n----------------------------------------------------------------";
636 for (
int i = 1;
i < 5; ++
i)
637 badModules <<
"\nTIB Layer " <<
i <<
" :" << ssV[0][
i].
str();
639 for (
int i = 1;
i < 4; ++
i)
640 badModules <<
"\nTID+ Disk " <<
i <<
" :" << ssV[1][
i].
str();
641 for (
int i = 4;
i < 7; ++
i)
642 badModules <<
"\nTID- Disk " <<
i - 3 <<
" :" << ssV[1][
i].
str();
644 for (
int i = 1;
i < 7; ++
i)
645 badModules <<
"\nTOB Layer " <<
i <<
" :" << ssV[2][
i].
str();
647 for (
int i = 1;
i < 10; ++
i)
648 badModules <<
"\nTEC+ Disk " <<
i <<
" :" << ssV[3][
i].
str();
649 for (
int i = 10;
i < 19; ++
i)
650 badModules <<
"\nTEC- Disk " <<
i - 9 <<
" :" << ssV[3][
i].
str();
656 desc.
add<
bool>(
"isAtPCL",
false);
657 desc.
add<
bool>(
"doStoreOnDB",
false);
659 desc.
add<
double>(
"Threshold", 0.1);
661 desc.
add<
int>(
"nModsMin", 5);
bool checkMapsValidity(const std::vector< MonitorElement * > &maps, const std::string &type) const
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
const std::vector< BadComponent > & getBadComponentList() const
constexpr char const * layerName[numberOfLayers]
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
const edm::ESGetToken< TrackerTopology, TrackerTopologyRcd > tTopoToken_
std::unique_ptr< TrackerTopology > tTopo_
void endRun(edm::Run const &, edm::EventSetup const &) override
#define DEFINE_FWK_MODULE(type)
virtual TH1F * getTH1F() const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const bool autoIneffModTagging_
void writeBadStripPayload(const SiStripQuality &quality) const
Log< level::Error, false > LogError
tuple nStrips
1.2 is to make the matching window safely the two nearest strips 0.35 is the size of an ME0 chamber i...
void makeSummary(DQMStore::IGetter &getter, TFileService &fs) const
const edm::ESGetToken< SiStripQuality, SiStripQualityRcd > stripQualityToken_
constexpr std::array< uint8_t, layerIndexSize > layer
const uint16_t range(const Frame &aFrame)
bool getData(T &iHolder) const
const edm::ESGetToken< TrackerGeometry, TrackerDigiGeometryRecord > tkGeomToken_
RegistryIterator getRegistryVectorEnd() const
std::unique_ptr< TkDetMap > tkDetMap_
void printTotalStatistics(const std::array< long, 23 > &layerFound, const std::array< long, 23 > &layerTotal) const
std::vector< DetId > stripDetIds_
virtual MonitorElement * get(std::string const &fullpath) const
void makeSummaryVsBX(DQMStore::IGetter &getter, TFileService &fs) const
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
Hash writeOneIOV(const T &payload, Time_t time, const std::string &recordName)
const std::pair< unsigned short, double > getNumberOfApvsAndStripLength(uint32_t detId) const
SiStripDetInfo read(std::string filePath)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
SiStripHitEfficiencyHarvester(const edm::ParameterSet &)
ContainerIterator getDataVectorBegin() const
Log< level::Info, false > LogInfo
void setBadComponents(int i, int component, const SiStripQuality::BadComponent &BC, int NBadComponent[4][19][4])
void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override
Constants and enumerated types for FED/FEC systems.
~SiStripHitEfficiencyHarvester() override=default
const unsigned int nTEClayers_
RegistryIterator getRegistryVectorBegin() const
static const uint16_t STRIPS_PER_APV
void makeSummaryVsCM(DQMStore::IGetter &getter, TFileService &fs) const
cond::Time_t currentTime() const
std::pair< ContainerIterator, ContainerIterator > Range
std::unique_ptr< SiStripQuality > stripQuality_
static constexpr char const *const kDefaultFile
void makeSummaryVsLumi(DQMStore::IGetter &getter) const
void printAndWriteBadModules(const SiStripQuality &quality, const SiStripDetInfo &detInfo) const
const edm::ESGetToken< TkDetMap, TrackerTopologyRcd > tkDetMapToken_
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
const std::string record_
data decode(const unsigned int &value) const