File indexing completed on 2024-04-06 11:59:55
0001
0002 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0003 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0004 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0005 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEffData.h"
0006 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h"
0007 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0008 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0009 #include "DQM/SiStripCommon/interface/TkHistoMap.h"
0010 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0011 #include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" /* for STRIPS_PER_APV*/
0012 #include "FWCore/Framework/interface/ESWatcher.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterDescription.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0022 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0023 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0024
0025
0026 #include <boost/type_index.hpp>
0027 #include <fmt/printf.h>
0028 #include <numeric> // for std::accumulate
0029 #include <sstream>
0030
0031
0032 #include "TCanvas.h"
0033 #include "TEfficiency.h"
0034 #include "TGraphAsymmErrors.h"
0035 #include "TLegend.h"
0036 #include "TStyle.h"
0037 #include "TTree.h"
0038
0039
0040 #define LOGPRINT edm::LogPrint("SiStripHitEfficiencyHarvester")
0041
0042 class SiStripHitEfficiencyHarvester : public DQMEDHarvester {
0043 public:
0044 explicit SiStripHitEfficiencyHarvester(const edm::ParameterSet&);
0045 ~SiStripHitEfficiencyHarvester() override = default;
0046 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048 void endRun(edm::Run const&, edm::EventSetup const&) override;
0049 void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0050
0051 private:
0052 SiStripHitEffData calibData_;
0053 const std::string inputFolder_;
0054 const bool isAtPCL_;
0055 const bool autoIneffModTagging_, doStoreOnDB_;
0056 const bool doStoreOnTree_;
0057 const bool showRings_, showEndcapSides_, showTOB6TEC9_, showOnlyGoodModules_;
0058 const unsigned int nTEClayers_;
0059 const double threshold_;
0060 const int nModsMin_;
0061 const float effPlotMin_;
0062 const double tkMapMin_;
0063 const std::string title_, record_;
0064
0065 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0066 const edm::ESGetToken<TkDetMap, TrackerTopologyRcd> tkDetMapToken_;
0067 const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> stripQualityToken_;
0068 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0069
0070 std::unique_ptr<TrackerTopology> tTopo_;
0071 std::unique_ptr<TkDetMap> tkDetMap_;
0072 std::unique_ptr<SiStripQuality> stripQuality_;
0073 std::vector<DetId> stripDetIds_;
0074
0075 int goodlayertotal[bounds::k_END_OF_LAYS_AND_RINGS];
0076 int goodlayerfound[bounds::k_END_OF_LAYS_AND_RINGS];
0077 int alllayertotal[bounds::k_END_OF_LAYS_AND_RINGS];
0078 int alllayerfound[bounds::k_END_OF_LAYS_AND_RINGS];
0079
0080
0081 TTree* tree;
0082 unsigned int t_DetId, t_found, t_total;
0083 unsigned char t_layer;
0084 bool t_isTaggedIneff;
0085 float t_threshold;
0086
0087 void writeBadStripPayload(const SiStripQuality& quality) const;
0088 void printTotalStatistics(const std::array<long, bounds::k_END_OF_LAYERS>& layerFound,
0089 const std::array<long, bounds::k_END_OF_LAYERS>& layerTotal) const;
0090 void printAndWriteBadModules(const SiStripQuality& quality, const SiStripDetInfo& detInfo) const;
0091 bool checkMapsValidity(const std::vector<MonitorElement*>& maps, const std::string& type) const;
0092 unsigned int countTotalHits(const std::vector<MonitorElement*>& maps);
0093 void makeSummary(DQMStore::IGetter& getter, DQMStore::IBooker& booker, bool doProfiles = false) const;
0094 void makeSummaryVsVariable(DQMStore::IGetter& getter,
0095 DQMStore::IBooker& booker,
0096 ::projections theProj,
0097 bool doProfiles = true) const;
0098 template <typename T>
0099 void setEffBinLabels(const T gr, const T gr2, const unsigned int nLayers) const;
0100 };
0101
0102 SiStripHitEfficiencyHarvester::SiStripHitEfficiencyHarvester(const edm::ParameterSet& conf)
0103 : inputFolder_(conf.getParameter<std::string>("inputFolder")),
0104 isAtPCL_(conf.getParameter<bool>("isAtPCL")),
0105 autoIneffModTagging_(conf.getUntrackedParameter<bool>("AutoIneffModTagging", false)),
0106 doStoreOnDB_(conf.getParameter<bool>("doStoreOnDB")),
0107 doStoreOnTree_(conf.getUntrackedParameter<bool>("doStoreOnTree")),
0108 showRings_(conf.getUntrackedParameter<bool>("ShowRings", false)),
0109 showEndcapSides_(conf.getUntrackedParameter<bool>("ShowEndcapSides", true)),
0110 showTOB6TEC9_(conf.getUntrackedParameter<bool>("ShowTOB6TEC9", false)),
0111 showOnlyGoodModules_(conf.getUntrackedParameter<bool>("ShowOnlyGoodModules", false)),
0112 nTEClayers_(showRings_ ? 7 : 9),
0113 threshold_(conf.getParameter<double>("Threshold")),
0114 nModsMin_(conf.getParameter<int>("nModsMin")),
0115 effPlotMin_(conf.getUntrackedParameter<double>("EffPlotMin", 0.9)),
0116 tkMapMin_(conf.getUntrackedParameter<double>("TkMapMin", 0.9)),
0117 title_(conf.getParameter<std::string>("Title")),
0118 record_(conf.getParameter<std::string>("Record")),
0119 tTopoToken_(esConsumes<edm::Transition::EndRun>()),
0120 tkDetMapToken_(esConsumes<edm::Transition::EndRun>()),
0121 stripQualityToken_(esConsumes<edm::Transition::EndRun>()),
0122 tkGeomToken_(esConsumes<edm::Transition::EndRun>()) {
0123
0124 for (int l = 0; l < bounds::k_END_OF_LAYS_AND_RINGS; l++) {
0125 goodlayertotal[l] = 0;
0126 goodlayerfound[l] = 0;
0127 alllayertotal[l] = 0;
0128 alllayerfound[l] = 0;
0129 }
0130 }
0131
0132 void SiStripHitEfficiencyHarvester::endRun(edm::Run const&, edm::EventSetup const& iSetup) {
0133 if (!tTopo_) {
0134 tTopo_ = std::make_unique<TrackerTopology>(iSetup.getData(tTopoToken_));
0135 }
0136 if (!tkDetMap_) {
0137 tkDetMap_ = std::make_unique<TkDetMap>(iSetup.getData(tkDetMapToken_));
0138 }
0139 if (!stripQuality_) {
0140 stripQuality_ = std::make_unique<SiStripQuality>(iSetup.getData(stripQualityToken_));
0141 }
0142 if (stripDetIds_.empty()) {
0143 const auto& tkGeom = iSetup.getData(tkGeomToken_);
0144 for (const auto& det : tkGeom.detUnits()) {
0145 if (dynamic_cast<const StripGeomDetUnit*>(det)) {
0146 stripDetIds_.push_back(det->geographicalId());
0147 }
0148 }
0149 }
0150 }
0151
0152 bool SiStripHitEfficiencyHarvester::checkMapsValidity(const std::vector<MonitorElement*>& maps,
0153 const std::string& type) const {
0154 std::vector<bool> isThere;
0155 isThere.reserve(maps.size());
0156 std::transform(maps.begin() + 1, maps.end(), std::back_inserter(isThere), [](auto& x) { return !(x == nullptr); });
0157
0158 int count{0};
0159 for (const auto& it : isThere) {
0160 count++;
0161 LogDebug("SiStripHitEfficiencyHarvester") << " layer: " << count << " " << it << std::endl;
0162 if (it)
0163 LogDebug("SiStripHitEfficiencyHarvester") << "resolving to " << maps[count]->getName() << std::endl;
0164 }
0165
0166
0167 bool areMapsAvailable{true};
0168 int layerCount{0};
0169 for (const auto& it : isThere) {
0170 layerCount++;
0171 if (!it) {
0172 edm::LogError("SiStripHitEfficiencyHarvester")
0173 << type << " TkHistoMap for layer " << layerCount << " was not found.\n -> Aborting!";
0174 areMapsAvailable = false;
0175 break;
0176 }
0177 }
0178 return areMapsAvailable;
0179 }
0180
0181 unsigned int SiStripHitEfficiencyHarvester::countTotalHits(const std::vector<MonitorElement*>& maps) {
0182 return std::accumulate(maps.begin() + 1, maps.end(), 0, [](unsigned int total, MonitorElement* item) {
0183 return total + item->getEntries();
0184 });
0185 }
0186
0187 void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStore::IGetter& getter) {
0188 if (!isAtPCL_) {
0189 edm::Service<TFileService> fs;
0190 if (!fs.isAvailable()) {
0191 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0192 << "please add it to config file";
0193 }
0194
0195 if (doStoreOnTree_) {
0196
0197 tree = fs->make<TTree>("ModEff", "ModEff");
0198 tree->Branch("DetId", &t_DetId, "DetId/i");
0199 tree->Branch("Layer", &t_layer, "Layer/b");
0200 tree->Branch("FoundHits", &t_found, "FoundHits/i");
0201 tree->Branch("AllHits", &t_total, "AllHits/i");
0202 tree->Branch("IsTaggedIneff", &t_isTaggedIneff, "IsTaggedIneff/O");
0203 tree->Branch("TagThreshold", &t_threshold, "TagThreshold/F");
0204 }
0205 }
0206
0207 if (!autoIneffModTagging_)
0208 LOGPRINT << "A module is bad if efficiency < " << threshold_ << " and has at least " << nModsMin_ << " nModsMin.";
0209 else
0210 LOGPRINT << "A module is bad if the upper limit on the efficiency is < to the avg in the layer - " << threshold_
0211 << " and has at least " << nModsMin_ << " nModsMin.";
0212
0213 auto h_module_total = std::make_unique<TkHistoMap>(tkDetMap_.get());
0214 h_module_total->loadTkHistoMap(fmt::format("{}/TkDetMaps", inputFolder_), "perModule_total");
0215 auto h_module_found = std::make_unique<TkHistoMap>(tkDetMap_.get());
0216 h_module_found->loadTkHistoMap(fmt::format("{}/TkDetMaps", inputFolder_), "perModule_found");
0217
0218
0219 const auto& totalMaps = h_module_total->getAllMaps();
0220 const auto& foundMaps = h_module_found->getAllMaps();
0221
0222 LogDebug("SiStripHitEfficiencyHarvester")
0223 << "totalMaps.size(): " << totalMaps.size() << " foundMaps.size() " << foundMaps.size() << std::endl;
0224
0225
0226 bool isTotalMapAvailable = this->checkMapsValidity(totalMaps, std::string("Total"));
0227 bool isFoundMapAvailable = this->checkMapsValidity(foundMaps, std::string("Found"));
0228
0229 LogDebug("SiStripHitEfficiencyHarvester")
0230 << "isTotalMapAvailable: " << isTotalMapAvailable << " isFoundMapAvailable " << isFoundMapAvailable << std::endl;
0231
0232
0233 if (!isTotalMapAvailable or !isFoundMapAvailable)
0234 return;
0235
0236 LogDebug("SiStripHitEfficiencyHarvester")
0237 << "Entries in total TkHistoMap for layer 3: " << h_module_total->getMap(3)->getEntries() << ", found "
0238 << h_module_found->getMap(3)->getEntries();
0239
0240
0241 const unsigned int totalHits = this->countTotalHits(totalMaps);
0242
0243
0244 for (size_t i = 1; i < totalMaps.size(); i++) {
0245 h_module_total->getMap(i)->setOption("colz");
0246 h_module_found->getMap(i)->setOption("colz");
0247 }
0248
0249
0250 booker.setCurrentFolder(inputFolder_);
0251
0252 std::vector<MonitorElement*> hEffInLayer(std::size_t(1), nullptr);
0253 hEffInLayer.reserve(bounds::k_END_OF_LAYERS);
0254 for (std::size_t i = 1; i != bounds::k_END_OF_LAYERS; ++i) {
0255 const auto lyrName = ::layerName(i, showRings_, nTEClayers_);
0256 hEffInLayer.push_back(booker.book1D(
0257 Form("eff_layer%i", int(i)), Form("Module efficiency in layer %s", lyrName.c_str()), 201, 0, 1.005));
0258 }
0259 std::array<long, bounds::k_END_OF_LAYERS> layerTotal{};
0260 std::array<long, bounds::k_END_OF_LAYERS> layerFound{};
0261 layerTotal.fill(0);
0262 layerFound.fill(0);
0263
0264
0265
0266
0267
0268 TrackerMap tkMap{" Detector Inefficiency "};
0269 TrackerMap tkMapBad{" Inefficient Modules "};
0270 TrackerMap tkMapEff{title_};
0271 TrackerMap tkMapNum{" Detector numerator "};
0272 TrackerMap tkMapDen{" Detector denominator "};
0273 std::map<unsigned int, double> badModules;
0274
0275
0276 const auto& EventStats = getter.get(fmt::format("{}/EventInfo/EventStats", inputFolder_));
0277 const int totalEvents = EventStats->getBinContent(1., 1.);
0278 calibData_.FEDErrorOccupancy = std::make_unique<TkHistoMap>(tkDetMap_.get());
0279 calibData_.FEDErrorOccupancy->loadTkHistoMap(fmt::format("{}/FEDErrorTkDetMaps", inputFolder_),
0280 "perModule_FEDErrors");
0281
0282
0283 calibData_.fillMapFromTkMap(totalEvents, 0.75, stripDetIds_);
0284
0285 for (const auto& [badId, fraction] : calibData_.fedErrorCounts) {
0286 LogDebug("SiStripHitEfficiencyHarvester")
0287 << __PRETTY_FUNCTION__ << " bad module from FEDError " << badId << "," << fraction << std::endl;
0288 }
0289
0290 for (auto det : stripDetIds_) {
0291 auto layer = ::checkLayer(det, tTopo_.get());
0292 const auto num = h_module_found->getValue(det);
0293 const auto denom = h_module_total->getValue(det);
0294 if (denom) {
0295
0296 if (stripQuality_->getBadApvs(det) == 0 && calibData_.checkFedError(det)) {
0297 const auto eff = num / denom;
0298 hEffInLayer[layer]->Fill(eff);
0299 if (!autoIneffModTagging_) {
0300 if ((denom >= nModsMin_) && (eff < threshold_)) {
0301
0302 badModules[det] = eff;
0303 tkMapBad.fillc(det, 255, 0, 0);
0304 LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ") module "
0305 << det.rawId() << " efficiency: " << eff << " , " << num << "/" << denom;
0306 } else {
0307
0308 tkMapBad.fillc(det, 255, 255, 255);
0309 }
0310 if (eff < threshold_)
0311 LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ") module "
0312 << det.rawId() << " efficiency: " << eff << " , " << num << "/" << denom;
0313
0314 if (denom < nModsMin_) {
0315 LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ") module "
0316 << det.rawId() << " is under occupancy at " << denom;
0317 }
0318
0319 if (doStoreOnTree_ && !isAtPCL_) {
0320 t_DetId = det.rawId();
0321 t_layer = layer;
0322 t_found = num;
0323 t_total = denom;
0324 t_isTaggedIneff = false;
0325 t_threshold = 0;
0326 tree->Fill();
0327 }
0328 }
0329 }
0330 }
0331 }
0332
0333 if (autoIneffModTagging_) {
0334 for (unsigned int i = 1; i <= k_LayersAtTECEnd; i++) {
0335
0336 hEffInLayer[i]->getTH1()->GetXaxis()->SetRange(
0337 3, hEffInLayer[i]->getNbinsX() + 1);
0338 const double layer_min_eff = hEffInLayer[i]->getMean() - std::max(2.5 * hEffInLayer[i]->getRMS(), threshold_);
0339 LOGPRINT << "Layer " << i << " threshold for bad modules: <" << layer_min_eff
0340 << " (layer mean: " << hEffInLayer[i]->getMean() << " rms: " << hEffInLayer[i]->getRMS() << ")";
0341
0342 hEffInLayer[i]->getTH1()->GetXaxis()->SetRange(1, hEffInLayer[i]->getNbinsX() + 1);
0343
0344 for (auto det : stripDetIds_) {
0345 const auto layer = ::checkLayer(det, tTopo_.get());
0346 if (layer == i) {
0347 const auto num = h_module_found->getValue(det);
0348 const auto denom = h_module_total->getValue(det);
0349 if (denom) {
0350 assert(num <= denom);
0351
0352 if (stripQuality_->getBadApvs(det) == 0 && calibData_.checkFedError(det)) {
0353 const auto eff = num / denom;
0354 const auto eff_up = TEfficiency::Bayesian(denom, num, .99, 1, 1, true);
0355
0356 if ((denom >= nModsMin_) && (eff_up < layer_min_eff)) {
0357
0358 badModules[det] = eff;
0359 tkMapBad.fillc(det, 255, 0, 0);
0360 if (!isAtPCL_ && doStoreOnTree_) {
0361 t_isTaggedIneff = true;
0362 }
0363 } else {
0364
0365 tkMapBad.fillc(det, 255, 255, 255);
0366 if (!isAtPCL_ && doStoreOnTree_) {
0367 t_isTaggedIneff = false;
0368 }
0369 }
0370 if (eff_up < layer_min_eff + 0.08) {
0371
0372 LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ") module "
0373 << det.rawId() << " efficiency: " << eff << " , " << num << "/" << denom
0374 << " , upper limit: " << eff_up;
0375 }
0376 if (denom < nModsMin_) {
0377 LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ") module "
0378 << det.rawId() << " layer " << layer << " is under occupancy at " << denom;
0379 }
0380
0381 if (!isAtPCL_ && doStoreOnTree_) {
0382 t_DetId = det.rawId();
0383 t_layer = layer;
0384 t_found = num;
0385 t_total = denom;
0386 t_threshold = layer_min_eff;
0387 tree->Fill();
0388 }
0389
0390
0391 tkMap.fill(det, 1. - eff);
0392 tkMapEff.fill(det, eff);
0393 tkMapNum.fill(det, num);
0394 tkMapDen.fill(det, denom);
0395
0396 layerTotal[layer] += denom;
0397 layerFound[layer] += num;
0398
0399
0400
0401 if (layer <= bounds::k_LayersAtTOBEnd) {
0402 goodlayerfound[layer] += num;
0403 goodlayertotal[layer] += denom;
0404 } else if (layer <= bounds::k_LayersAtTIDEnd) {
0405 if (tTopo_->tidSide(det) == 1) {
0406 goodlayerfound[layer] += num;
0407 goodlayertotal[layer] += denom;
0408 } else if (tTopo_->tidSide(det) == 2) {
0409 goodlayerfound[layer + 3] += num;
0410 goodlayertotal[layer + 3] += denom;
0411 }
0412 } else if (layer <= bounds::k_LayersAtTECEnd) {
0413 if (tTopo_->tecSide(det) == 1) {
0414 goodlayerfound[layer + 3] += num;
0415 goodlayertotal[layer + 3] += denom;
0416 } else if (tTopo_->tecSide(det) == 2) {
0417 goodlayerfound[layer + 3 + nTEClayers_] += num;
0418 goodlayertotal[layer + 3 + nTEClayers_] += denom;
0419 }
0420 }
0421 }
0422
0423
0424 if (layer <= bounds::k_LayersAtTOBEnd) {
0425 alllayerfound[layer] += num;
0426 alllayertotal[layer] += denom;
0427 } else if (layer <= bounds::k_LayersAtTIDEnd) {
0428 if (tTopo_->tidSide(det) == 1) {
0429 alllayerfound[layer] += num;
0430 alllayertotal[layer] += denom;
0431 } else if (tTopo_->tidSide(det) == 2) {
0432 alllayerfound[layer + 3] += num;
0433 alllayertotal[layer + 3] += denom;
0434 }
0435 } else if (layer <= bounds::k_LayersAtTECEnd) {
0436 if (tTopo_->tecSide(det) == 1) {
0437 alllayerfound[layer + 3] += num;
0438 alllayertotal[layer + 3] += denom;
0439 } else if (tTopo_->tecSide(det) == 2) {
0440 alllayerfound[layer + 3 + nTEClayers_] += num;
0441 alllayertotal[layer + 3 + nTEClayers_] += denom;
0442 }
0443 }
0444 }
0445 }
0446 }
0447 }
0448 }
0449
0450 tkMap.save(true, 0, 0, "SiStripHitEffTKMap_NEW.png");
0451 tkMapBad.save(true, 0, 0, "SiStripHitEffTKMapBad_NEW.png");
0452 tkMapEff.save(true, tkMapMin_, 1., "SiStripHitEffTKMapEff_NEW.png");
0453 tkMapNum.save(true, 0, 0, "SiStripHitEffTKMapNum_NEW.png");
0454 tkMapDen.save(true, 0, 0, "SiStripHitEffTKMapDen_NEW.png");
0455
0456 const auto detInfo =
0457 SiStripDetInfoFileReader::read(edm::FileInPath{SiStripDetInfoFileReader::kDefaultFile}.fullPath());
0458 SiStripQuality pQuality{detInfo};
0459
0460
0461
0462 for (const auto it : badModules) {
0463 const auto det = it.first;
0464 std::vector<unsigned int> badStripList;
0465
0466
0467 const auto nStrips = detInfo.getNumberOfApvsAndStripLength(det).first * sistrip::STRIPS_PER_APV;
0468 LOGPRINT << "Number of strips module " << det << " is " << nStrips;
0469 badStripList.push_back(pQuality.encode(0, nStrips, 0));
0470
0471 LOGPRINT << "ID1 should match list of modules above " << det;
0472 pQuality.compact(det, badStripList);
0473 pQuality.put(det, SiStripQuality::Range(badStripList.begin(), badStripList.end()));
0474 }
0475 pQuality.fillBadComponents();
0476 if (doStoreOnDB_) {
0477 if (totalHits > 0u) {
0478 writeBadStripPayload(pQuality);
0479 } else {
0480 edm::LogPrint("SiStripHitEfficiencyHarvester")
0481 << __PRETTY_FUNCTION__ << " There are no SiStrip hits for a valid measurement, skipping!";
0482 }
0483 } else {
0484 edm::LogInfo("SiStripHitEfficiencyHarvester") << "Will not produce payload!";
0485 }
0486
0487 printTotalStatistics(layerFound, layerTotal);
0488
0489
0490 printAndWriteBadModules(pQuality, detInfo);
0491
0492
0493 makeSummary(getter, booker);
0494 makeSummaryVsVariable(getter, booker, projections::k_vs_LUMI);
0495 makeSummaryVsVariable(getter, booker, projections::k_vs_PU);
0496 makeSummaryVsVariable(getter, booker, projections::k_vs_BX);
0497 }
0498
0499 void SiStripHitEfficiencyHarvester::printTotalStatistics(
0500 const std::array<long, bounds::k_END_OF_LAYERS>& layerFound,
0501 const std::array<long, bounds::k_END_OF_LAYERS>& layerTotal) const {
0502
0503 int totalfound = 0;
0504 int totaltotal = 0;
0505 double layereff;
0506 int subdetfound[5] = {0, 0, 0, 0, 0};
0507 int subdettotal[5] = {0, 0, 0, 0, 0};
0508
0509 for (unsigned int i = 1; i <= bounds::k_LayersAtTECEnd; i++) {
0510 layereff = double(layerFound[i]) / double(layerTotal[i]);
0511 LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers_) << ") has total efficiency "
0512 << layereff << " " << layerFound[i] << "/" << layerTotal[i];
0513 totalfound += layerFound[i];
0514 totaltotal += layerTotal[i];
0515 if (i <= bounds::k_LayersAtTIBEnd) {
0516 subdetfound[1] += layerFound[i];
0517 subdettotal[1] += layerTotal[i];
0518 }
0519 if (i > bounds::k_LayersAtTIBEnd && i <= bounds::k_LayersAtTOBEnd) {
0520 subdetfound[2] += layerFound[i];
0521 subdettotal[2] += layerTotal[i];
0522 }
0523 if (i > bounds::k_LayersAtTOBEnd && i <= bounds::k_LayersAtTIDEnd) {
0524 subdetfound[3] += layerFound[i];
0525 subdettotal[3] += layerTotal[i];
0526 }
0527 if (i > bounds::k_LayersAtTIDEnd) {
0528 subdetfound[4] += layerFound[i];
0529 subdettotal[4] += layerTotal[i];
0530 }
0531 }
0532
0533 LOGPRINT << "The total efficiency is " << double(totalfound) / double(totaltotal);
0534 LOGPRINT << " TIB: " << double(subdetfound[1]) / subdettotal[1] << " " << subdetfound[1] << "/"
0535 << subdettotal[1];
0536 LOGPRINT << " TOB: " << double(subdetfound[2]) / subdettotal[2] << " " << subdetfound[2] << "/"
0537 << subdettotal[2];
0538 LOGPRINT << " TID: " << double(subdetfound[3]) / subdettotal[3] << " " << subdetfound[3] << "/"
0539 << subdettotal[3];
0540 LOGPRINT << " TEC: " << double(subdetfound[4]) / subdettotal[4] << " " << subdetfound[4] << "/"
0541 << subdettotal[4];
0542 }
0543
0544 void SiStripHitEfficiencyHarvester::writeBadStripPayload(const SiStripQuality& quality) const {
0545 SiStripBadStrip pBadStrip{};
0546 const auto pQdvBegin = quality.getDataVectorBegin();
0547 for (auto rIt = quality.getRegistryVectorBegin(); rIt != quality.getRegistryVectorEnd(); ++rIt) {
0548 const auto range = SiStripBadStrip::Range(pQdvBegin + rIt->ibegin, pQdvBegin + rIt->iend);
0549 if (!pBadStrip.put(rIt->detid, range))
0550 edm::LogError("SiStripHitEfficiencyHarvester") << "detid already exists in SiStripBadStrip";
0551 }
0552 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0553 if (poolDbService.isAvailable()) {
0554 poolDbService->writeOneIOV(pBadStrip, poolDbService->currentTime(), record_);
0555 } else {
0556 throw cms::Exception("PoolDBService required");
0557 }
0558 }
0559
0560 void SiStripHitEfficiencyHarvester::makeSummary(DQMStore::IGetter& getter,
0561 DQMStore::IBooker& booker,
0562 bool doProfiles) const {
0563
0564 unsigned int nLayers{34};
0565 if (showRings_)
0566 nLayers = 30;
0567 if (!showEndcapSides_) {
0568 if (!showRings_)
0569 nLayers = 22;
0570 else
0571 nLayers = 20;
0572 }
0573
0574
0575 booker.setCurrentFolder(fmt::format("{}/EfficiencySummary", inputFolder_));
0576 MonitorElement* found_good =
0577 booker.book1D("found_good", "found hits per layer (good modules only)", nLayers + 1, 0, nLayers + 1);
0578 MonitorElement* all_good =
0579 booker.book1D("all_good", "all hits per layer (good modules only)", nLayers + 1, 0, nLayers + 1);
0580 MonitorElement* found_all =
0581 booker.book1D("found_all", "found hit per layer (all modules)", nLayers + 1, 0, nLayers + 1);
0582 MonitorElement* all_all = booker.book1D("all_all", "all hits per layer (all modules)", nLayers + 1, 0, nLayers + 1);
0583
0584
0585 found_good->setBinContent(0, -1);
0586 all_good->setBinContent(0, 1);
0587
0588
0589 for (unsigned int i = 1; i < nLayers + 2; ++i) {
0590 found_good->setBinContent(i, 1e-6);
0591 all_good->setBinContent(i, 1);
0592 found_all->setBinContent(i, 1e-6);
0593 all_all->setBinContent(i, 1);
0594 }
0595
0596 TCanvas* c7 = new TCanvas("c7", " test ", 10, 10, 800, 600);
0597 c7->SetFillColor(0);
0598 c7->SetGrid();
0599
0600 unsigned int nLayers_max = nLayers + 1;
0601 if (!showEndcapSides_)
0602 nLayers_max = 11;
0603 for (unsigned int i = 1; i < nLayers_max; ++i) {
0604 LOGPRINT << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i]
0605 << " B = " << goodlayertotal[i];
0606 if (goodlayertotal[i] > 5) {
0607 found_good->setBinContent(i, goodlayerfound[i]);
0608 all_good->setBinContent(i, goodlayertotal[i]);
0609 }
0610
0611 LOGPRINT << "Filling all modules layer " << i << ": S = " << alllayerfound[i] << " B = " << alllayertotal[i];
0612 if (alllayertotal[i] > 5) {
0613 found_all->setBinContent(i, alllayerfound[i]);
0614 all_all->setBinContent(i, alllayertotal[i]);
0615 }
0616 }
0617
0618
0619 if (!showEndcapSides_) {
0620 for (unsigned int i = 11; i < 14; ++i) {
0621 LOGPRINT << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] + goodlayerfound[i + 3]
0622 << " B = " << goodlayertotal[i] + goodlayertotal[i + 3];
0623 if (goodlayertotal[i] + goodlayertotal[i + 3] > 5) {
0624 found_good->setBinContent(i, goodlayerfound[i] + goodlayerfound[i + 3]);
0625 all_good->setBinContent(i, goodlayertotal[i] + goodlayertotal[i + 3]);
0626 }
0627 LOGPRINT << "Filling all modules layer " << i << ": S = " << alllayerfound[i] + alllayerfound[i + 3]
0628 << " B = " << alllayertotal[i] + alllayertotal[i + 3];
0629 if (alllayertotal[i] + alllayertotal[i + 3] > 5) {
0630 found_all->setBinContent(i, alllayerfound[i] + alllayerfound[i + 3]);
0631 all_all->setBinContent(i, alllayertotal[i] + alllayertotal[i + 3]);
0632 }
0633 }
0634 for (unsigned int i = 17; i < 17 + nTEClayers_; ++i) {
0635 LOGPRINT << "Fill only good modules layer " << i - 3
0636 << ": S = " << goodlayerfound[i] + goodlayerfound[i + nTEClayers_]
0637 << " B = " << goodlayertotal[i] + goodlayertotal[i + nTEClayers_];
0638 if (goodlayertotal[i] + goodlayertotal[i + nTEClayers_] > 5) {
0639 found_good->setBinContent(i - 3, goodlayerfound[i] + goodlayerfound[i + nTEClayers_]);
0640 all_good->setBinContent(i - 3, goodlayertotal[i] + goodlayertotal[i + nTEClayers_]);
0641 }
0642 LOGPRINT << "Filling all modules layer " << i - 3
0643 << ": S = " << alllayerfound[i] + alllayerfound[i + nTEClayers_]
0644 << " B = " << alllayertotal[i] + alllayertotal[i + nTEClayers_];
0645 if (alllayertotal[i] + alllayertotal[i + nTEClayers_] > 5) {
0646 found_all->setBinContent(i - 3, alllayerfound[i] + alllayerfound[i + nTEClayers_]);
0647 all_all->setBinContent(i - 3, alllayertotal[i] + alllayertotal[i + nTEClayers_]);
0648 }
0649 }
0650 }
0651
0652 found_good->getTH1F()->Sumw2();
0653 all_good->getTH1F()->Sumw2();
0654
0655 found_all->getTH1F()->Sumw2();
0656 all_all->getTH1F()->Sumw2();
0657
0658 MonitorElement* h_eff_all =
0659 booker.book1D("eff_all", "Strip hit efficiency for all modules", nLayers + 1, 0, nLayers + 1);
0660 MonitorElement* h_eff_good =
0661 booker.book1D("eff_good", "Strip hit efficiency for good modules", nLayers + 1, 0, nLayers + 1);
0662
0663 if (doProfiles) {
0664
0665 TProfile* profile_all = ::computeEff(found_all->getTH1F(), all_all->getTH1F(), "all");
0666 profile_all->SetMinimum(tkMapMin_);
0667 profile_all->SetTitle("Strip hit efficiency for all modules");
0668 booker.bookProfile(profile_all->GetName(), profile_all);
0669
0670 TProfile* profile_good = ::computeEff(found_good->getTH1F(), all_good->getTH1F(), "good");
0671 profile_good->SetMinimum(tkMapMin_);
0672 profile_good->SetTitle("Strip hit efficiency for good modules");
0673 booker.bookProfile(profile_good->GetName(), profile_good);
0674
0675
0676 delete profile_all;
0677 delete profile_good;
0678 }
0679
0680 for (int i = 1; i < found_good->getNbinsX(); i++) {
0681 const auto& den_all = all_all->getBinContent(i);
0682 const auto& num_all = found_all->getBinContent(i);
0683 const auto& den_good = all_good->getBinContent(i);
0684 const auto& num_good = found_good->getBinContent(i);
0685
0686
0687 if (den_all > 0.) {
0688
0689
0690
0691
0692
0693 const auto& effPair_all = ::computeCPEfficiency(num_all, den_all);
0694 h_eff_all->setBinContent(i, effPair_all.value());
0695 h_eff_all->setBinError(i, effPair_all.error());
0696 }
0697
0698
0699 if (den_good > 0.) {
0700
0701
0702
0703
0704
0705 const auto& effPair_good = ::computeCPEfficiency(num_good, den_good);
0706 h_eff_good->setBinContent(i, effPair_good.value());
0707 h_eff_good->setBinError(i, effPair_good.error());
0708 }
0709 }
0710
0711 h_eff_all->getTH1F()->SetMinimum(effPlotMin_);
0712 h_eff_good->getTH1F()->SetMinimum(effPlotMin_);
0713
0714
0715 this->setEffBinLabels(h_eff_all->getTH1F(), h_eff_good->getTH1F(), nLayers);
0716
0717 if (!isAtPCL_) {
0718
0719 edm::Service<TFileService> fs;
0720 if (!fs.isAvailable()) {
0721 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0722 << "please add it to config file";
0723 }
0724
0725 TGraphAsymmErrors* gr = (*fs).make<TGraphAsymmErrors>(nLayers + 1);
0726 gr->SetName("eff_good");
0727 gr->BayesDivide(found_good->getTH1F(), all_good->getTH1F());
0728
0729 TGraphAsymmErrors* gr2 = (*fs).make<TGraphAsymmErrors>(nLayers + 1);
0730 gr2->SetName("eff_all");
0731 gr2->BayesDivide(found_all->getTH1F(), all_all->getTH1F());
0732
0733 for (unsigned int j = 0; j < nLayers + 1; j++) {
0734 gr->SetPointError(j, 0., 0., gr->GetErrorYlow(j), gr->GetErrorYhigh(j));
0735 gr2->SetPointError(j, 0., 0., gr2->GetErrorYlow(j), gr2->GetErrorYhigh(j));
0736 }
0737
0738 this->setEffBinLabels(gr, gr2, nLayers);
0739
0740 gr->GetXaxis()->SetLimits(0, nLayers);
0741 gr->SetMarkerColor(2);
0742 gr->SetMarkerSize(1.2);
0743 gr->SetLineColor(2);
0744 gr->SetLineWidth(4);
0745 gr->SetMarkerStyle(20);
0746 gr->SetMinimum(effPlotMin_);
0747 gr->SetMaximum(1.001);
0748 gr->GetYaxis()->SetTitle("Efficiency");
0749 gStyle->SetTitleFillColor(0);
0750 gStyle->SetTitleBorderSize(0);
0751 gr->SetTitle("SiStripHitEfficiency by Layer");
0752
0753 gr2->GetXaxis()->SetLimits(0, nLayers);
0754 gr2->SetMarkerColor(1);
0755 gr2->SetMarkerSize(1.2);
0756 gr2->SetLineColor(1);
0757 gr2->SetLineWidth(4);
0758 gr2->SetMarkerStyle(21);
0759 gr2->SetMinimum(effPlotMin_);
0760 gr2->SetMaximum(1.001);
0761 gr2->GetYaxis()->SetTitle("Efficiency");
0762 gr2->SetTitle("SiStripHitEfficiency by Layer");
0763
0764 gr->Draw("AP");
0765 gr->GetXaxis()->SetNdivisions(36);
0766
0767 c7->cd();
0768 TPad* overlay = new TPad("overlay", "", 0, 0, 1, 1);
0769 overlay->SetFillStyle(4000);
0770 overlay->SetFillColor(0);
0771 overlay->SetFrameFillStyle(4000);
0772 overlay->Draw("same");
0773 overlay->cd();
0774 if (!showOnlyGoodModules_)
0775 gr2->Draw("AP");
0776
0777 TLegend* leg = new TLegend(0.70, 0.27, 0.88, 0.40);
0778 leg->AddEntry(gr, "Good Modules", "p");
0779 if (!showOnlyGoodModules_)
0780 leg->AddEntry(gr2, "All Modules", "p");
0781 leg->SetTextSize(0.020);
0782 leg->SetFillColor(0);
0783 leg->Draw("same");
0784
0785 c7->SaveAs("Summary.png");
0786 c7->SaveAs("Summary.root");
0787 }
0788 }
0789
0790 template <typename T>
0791 void SiStripHitEfficiencyHarvester::setEffBinLabels(const T gr, const T gr2, const unsigned int nLayers) const {
0792 LogDebug("SiStripHitEfficiencyHarvester")
0793 << "nLayers = " << nLayers << " number of bins, gr1: " << gr->GetXaxis()->GetNbins()
0794 << " number of bins, gr2: " << gr2->GetXaxis()->GetNbins() << " showRings: " << showRings_
0795 << " showEndcapSides: " << showEndcapSides_ << " type of object is "
0796 << boost::typeindex::type_id<T>().pretty_name();
0797
0798 for (unsigned int k = 1; k < nLayers + 1; k++) {
0799 std::string label{};
0800 if (showEndcapSides_)
0801 label = ::layerSideName(k, showRings_, nTEClayers_);
0802 else
0803 label = ::layerName(k, showRings_, nTEClayers_);
0804 if (!showTOB6TEC9_) {
0805 if (k == 10)
0806 label = "";
0807 if (!showRings_ && k == nLayers)
0808 label = "";
0809 if (!showRings_ && showEndcapSides_ && k == 25)
0810 label = "";
0811 }
0812
0813 int bin{-1};
0814 if constexpr (std::is_same_v<T, TGraphAsymmErrors*>) {
0815 edm::LogInfo("SiStripHitEfficiencyHarvester")
0816 << "class name: " << gr->ClassName() << " expected TGraphAsymErrors" << std::endl;
0817 if (!showRings_) {
0818 if (showEndcapSides_) {
0819 bin = (((k + 1) * 100 + 2) / (nLayers)-4);
0820 } else {
0821 bin = ((k + 1) * 100 / (nLayers)-6);
0822 }
0823 } else {
0824 if (showEndcapSides_) {
0825 bin = ((k + 1) * 100 / (nLayers)-4);
0826 } else {
0827 bin = ((k + 1) * 100 / (nLayers)-7);
0828 }
0829 }
0830 } else {
0831 edm::LogInfo("SiStripHitEfficiencyHarvester")
0832 << "class name: " << gr->ClassName() << " expected TH1F" << std::endl;
0833 bin = k;
0834 }
0835 gr->GetXaxis()->SetBinLabel(bin, label.data());
0836 gr2->GetXaxis()->SetBinLabel(bin, label.data());
0837 }
0838 }
0839
0840 void SiStripHitEfficiencyHarvester::makeSummaryVsVariable(DQMStore::IGetter& getter,
0841 DQMStore::IBooker& booker,
0842 ::projections theProj,
0843 bool doProfiles) const {
0844 std::vector<MonitorElement*> effVsVariable;
0845 effVsVariable.reserve(showRings_ ? 20 : 22);
0846
0847 const auto& folderString = ::projFolder[theProj];
0848 const auto& foundHistoString = ::projFoundHisto[theProj];
0849 const auto& totalHistoString = ::projTotalHisto[theProj];
0850 const auto& titleString = ::projTitle[theProj];
0851 const auto& titleXString = ::projXtitle[theProj];
0852
0853 LogDebug("SiStripHitEfficiencyHarvester")
0854 << " inside" << __PRETTY_FUNCTION__ << " from " << ::projFolder[theProj] << " " << __LINE__ << std::endl;
0855
0856 for (unsigned int iLayer = 1; iLayer != (showRings_ ? 20 : 22); ++iLayer) {
0857 LogDebug("SiStripHitEfficiencyHarvester")
0858 << "iLayer " << iLayer << " " << fmt::format("{}/{}/{}{}", inputFolder_, folderString, foundHistoString, iLayer)
0859 << std::endl;
0860
0861 const auto lyrName = ::layerName(iLayer, showRings_, nTEClayers_);
0862 auto hfound = getter.get(fmt::format("{}/{}/{}{}", inputFolder_, folderString, foundHistoString, iLayer));
0863 auto htotal = getter.get(fmt::format("{}/{}/{}{}", inputFolder_, folderString, totalHistoString, iLayer));
0864
0865 if (hfound == nullptr or htotal == nullptr) {
0866 if (hfound == nullptr)
0867 edm::LogError("SiStripHitEfficiencyHarvester")
0868 << fmt::format("{}/{}/{}{}", inputFolder_, folderString, foundHistoString, iLayer) << " was not found!";
0869 if (htotal == nullptr)
0870 edm::LogError("SiStripHitEfficiencyHarvester")
0871 << fmt::format("{}/{}/{}{}", inputFolder_, folderString, totalHistoString, iLayer) << " was not found!";
0872
0873 continue;
0874 }
0875
0876
0877 if (!hfound->getTH1F()->GetSumw2())
0878 hfound->getTH1F()->Sumw2();
0879 if (!htotal->getTH1F()->GetSumw2())
0880 htotal->getTH1F()->Sumw2();
0881
0882
0883 for (int i = 0; i != hfound->getNbinsX() + 1; ++i) {
0884 if (hfound->getBinContent(i) == 0)
0885 hfound->setBinContent(i, 1e-6);
0886 if (htotal->getBinContent(i) == 0)
0887 htotal->setBinContent(i, 1);
0888 }
0889 LogDebug("SiStripHitEfficiencyHarvester") << "Total hits for layer " << iLayer << " (" << folderString
0890 << "): " << htotal->getEntries() << ", found " << hfound->getEntries();
0891
0892 booker.setCurrentFolder(fmt::format("{}/EfficiencySummary{}", inputFolder_, folderString));
0893 effVsVariable[iLayer] = booker.book1D(
0894 fmt::sprintf("eff%sLayer%s", folderString, lyrName),
0895 fmt::sprintf("Efficiency vs %s for layer %s;%s;SiStrip Hit efficiency", titleString, lyrName, titleXString),
0896 hfound->getNbinsX(),
0897 hfound->getAxisMin(),
0898 hfound->getAxisMax());
0899
0900 effVsVariable[iLayer]->setOption("e");
0901
0902 LogDebug("SiStripHitEfficiencyHarvester")
0903 << " bin 0 " << hfound->getAxisMin() << " bin last: " << hfound->getAxisMax() << std::endl;
0904
0905 for (int i = 0; i != hfound->getNbinsX() + 1; ++i) {
0906 const auto& den = htotal->getBinContent(i);
0907 const auto& num = hfound->getBinContent(i);
0908
0909
0910 if (den > 0.) {
0911 const auto& effPair = ::computeCPEfficiency(num, den);
0912 effVsVariable[iLayer]->setBinContent(i, effPair.value());
0913 effVsVariable[iLayer]->setBinError(i, effPair.error());
0914
0915 LogDebug("SiStripHitEfficiencyHarvester")
0916 << __PRETTY_FUNCTION__ << " " << lyrName << " bin:" << i << " err:" << effPair.error() << std::endl;
0917 }
0918 }
0919
0920
0921 effVsVariable[iLayer]->getTH1F()->SetMinimum(tkMapMin_);
0922
0923 if (doProfiles) {
0924
0925 TProfile* profile = ::computeEff(hfound->getTH1F(), htotal->getTH1F(), lyrName);
0926 profile->SetOption("s");
0927 TString title =
0928 fmt::sprintf("Efficiency vs %s for layer %s;%s;SiStrip Hit efficiency", titleString, lyrName, titleXString);
0929 profile->SetMinimum(tkMapMin_);
0930
0931 profile->SetTitle(title.Data());
0932 booker.bookProfile(profile->GetName(), profile);
0933
0934 delete profile;
0935 }
0936 }
0937 }
0938
0939 namespace {
0940 void setBadComponents(int i,
0941 int comp,
0942 const SiStripQuality::BadComponent& bc,
0943 std::stringstream ssV[4][19],
0944 int nBad[4][19][4],
0945 int nAPV) {
0946 ssV[i][comp] << "\n\t\t " << bc.detid << " \t " << bc.BadModule << " \t " << ((bc.BadFibers) & 0x1) << " ";
0947 if (nAPV == 4)
0948 ssV[i][comp] << "x " << ((bc.BadFibers >> 1) & 0x1);
0949
0950 if (nAPV == 6)
0951 ssV[i][comp] << ((bc.BadFibers >> 1) & 0x1) << " " << ((bc.BadFibers >> 2) & 0x1);
0952 ssV[i][comp] << " \t " << ((bc.BadApvs) & 0x1) << " " << ((bc.BadApvs >> 1) & 0x1) << " ";
0953 if (nAPV == 4)
0954 ssV[i][comp] << "x x " << ((bc.BadApvs >> 2) & 0x1) << " " << ((bc.BadApvs >> 3) & 0x1);
0955 if (nAPV == 6)
0956 ssV[i][comp] << ((bc.BadApvs >> 2) & 0x1) << " " << ((bc.BadApvs >> 3) & 0x1) << " " << ((bc.BadApvs >> 4) & 0x1)
0957 << " " << ((bc.BadApvs >> 5) & 0x1) << " ";
0958
0959 if (bc.BadApvs) {
0960 nBad[i][0][2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
0961 ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
0962 nBad[i][comp][2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
0963 ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
0964 }
0965 if (bc.BadFibers) {
0966 nBad[i][0][1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
0967 nBad[i][comp][1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
0968 }
0969 if (bc.BadModule) {
0970 nBad[i][0][0]++;
0971 nBad[i][comp][0]++;
0972 }
0973 }
0974 }
0975
0976 void SiStripHitEfficiencyHarvester::printAndWriteBadModules(const SiStripQuality& quality,
0977 const SiStripDetInfo& detInfo) const {
0978
0979
0980
0981 int nTkBadComp[4];
0982 int nBadComp[4][19][4];
0983
0984
0985
0986 std::stringstream ssV[4][19];
0987
0988 for (int i = 0; i < 4; ++i) {
0989 nTkBadComp[i] = 0;
0990 for (int j = 0; j < 19; ++j) {
0991 ssV[i][j].str("");
0992 for (int k = 0; k < 4; ++k)
0993 nBadComp[i][j][k] = 0;
0994 }
0995 }
0996
0997 for (const auto& bc : quality.getBadComponentList()) {
0998
0999 if (bc.BadModule)
1000 nTkBadComp[0]++;
1001 if (bc.BadFibers)
1002 nTkBadComp[1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
1003 if (bc.BadApvs)
1004 nTkBadComp[2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
1005 ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
1006
1007 DetId det(bc.detid);
1008 if ((det.subdetId() >= SiStripSubdetector::TIB) && (det.subdetId() <= SiStripSubdetector::TEC)) {
1009 const auto nAPV = detInfo.getNumberOfApvsAndStripLength(det).first;
1010 switch (det.subdetId()) {
1011 case SiStripSubdetector::TIB:
1012 setBadComponents(0, tTopo_->tibLayer(det), bc, ssV, nBadComp, nAPV);
1013 break;
1014 case SiStripSubdetector::TID:
1015 setBadComponents(1,
1016 (tTopo_->tidSide(det) == 2 ? tTopo_->tidWheel(det) : tTopo_->tidWheel(det) + 3),
1017 bc,
1018 ssV,
1019 nBadComp,
1020 nAPV);
1021 break;
1022 case SiStripSubdetector::TOB:
1023 setBadComponents(2, tTopo_->tobLayer(det), bc, ssV, nBadComp, nAPV);
1024 break;
1025 case SiStripSubdetector::TEC:
1026 setBadComponents(3,
1027 (tTopo_->tecSide(det) == 2 ? tTopo_->tecWheel(det) : tTopo_->tecWheel(det) + 9),
1028 bc,
1029 ssV,
1030 nBadComp,
1031 nAPV);
1032 break;
1033 default:
1034 break;
1035 }
1036 }
1037 }
1038
1039 for (auto rp = quality.getRegistryVectorBegin(); rp != quality.getRegistryVectorEnd(); ++rp) {
1040 DetId det{rp->detid};
1041 int subdet = -999;
1042 int component = -999;
1043 switch (det.subdetId()) {
1044 case SiStripSubdetector::TIB:
1045 subdet = 0;
1046 component = tTopo_->tibLayer(det);
1047 break;
1048 case SiStripSubdetector::TID:
1049 subdet = 1;
1050 component = tTopo_->tidSide(det) == 2 ? tTopo_->tidWheel(det) : tTopo_->tidWheel(det) + 3;
1051 break;
1052 case SiStripSubdetector::TOB:
1053 subdet = 2;
1054 component = tTopo_->tobLayer(det);
1055 break;
1056 case SiStripSubdetector::TEC:
1057 subdet = 3;
1058 component = tTopo_->tecSide(det) == 2 ? tTopo_->tecWheel(det) : tTopo_->tecWheel(det) + 9;
1059 break;
1060 default:
1061 break;
1062 }
1063
1064 const auto pQdvBegin = quality.getDataVectorBegin();
1065 const auto sqrange = SiStripQuality::Range(pQdvBegin + rp->ibegin, pQdvBegin + rp->iend);
1066 float percentage = 0;
1067 for (int it = 0; it < sqrange.second - sqrange.first; it++) {
1068 unsigned int range = quality.decode(*(sqrange.first + it)).range;
1069 nTkBadComp[3] += range;
1070 nBadComp[subdet][0][3] += range;
1071 nBadComp[subdet][component][3] += range;
1072 percentage += range;
1073 }
1074 if (percentage != 0)
1075 percentage /= (sistrip::STRIPS_PER_APV * detInfo.getNumberOfApvsAndStripLength(det).first);
1076 if (percentage > 1)
1077 edm::LogError("SiStripHitEfficiencyHarvester") << "PROBLEM detid " << det.rawId() << " value " << percentage;
1078 }
1079
1080
1081 std::ostringstream ss;
1082 ss << "\n-----------------\nGlobal Info\n-----------------";
1083 ss << "\nBadComp \t Modules \tFibers "
1084 "\tApvs\tStrips\n----------------------------------------------------------------";
1085 ss << "\nTracker:\t\t" << nTkBadComp[0] << "\t" << nTkBadComp[1] << "\t" << nTkBadComp[2] << "\t" << nTkBadComp[3];
1086 ss << "\nTIB:\t\t\t" << nBadComp[0][0][0] << "\t" << nBadComp[0][0][1] << "\t" << nBadComp[0][0][2] << "\t"
1087 << nBadComp[0][0][3];
1088 ss << "\nTID:\t\t\t" << nBadComp[1][0][0] << "\t" << nBadComp[1][0][1] << "\t" << nBadComp[1][0][2] << "\t"
1089 << nBadComp[1][0][3];
1090 ss << "\nTOB:\t\t\t" << nBadComp[2][0][0] << "\t" << nBadComp[2][0][1] << "\t" << nBadComp[2][0][2] << "\t"
1091 << nBadComp[2][0][3];
1092 ss << "\nTEC:\t\t\t" << nBadComp[3][0][0] << "\t" << nBadComp[3][0][1] << "\t" << nBadComp[3][0][2] << "\t"
1093 << nBadComp[3][0][3];
1094 ss << "\n";
1095
1096 for (int i = 1; i < 5; ++i)
1097 ss << "\nTIB Layer " << i << " :\t\t" << nBadComp[0][i][0] << "\t" << nBadComp[0][i][1] << "\t" << nBadComp[0][i][2]
1098 << "\t" << nBadComp[0][i][3];
1099 ss << "\n";
1100 for (int i = 1; i < 4; ++i)
1101 ss << "\nTID+ Disk " << i << " :\t\t" << nBadComp[1][i][0] << "\t" << nBadComp[1][i][1] << "\t" << nBadComp[1][i][2]
1102 << "\t" << nBadComp[1][i][3];
1103 for (int i = 4; i < 7; ++i)
1104 ss << "\nTID- Disk " << i - 3 << " :\t\t" << nBadComp[1][i][0] << "\t" << nBadComp[1][i][1] << "\t"
1105 << nBadComp[1][i][2] << "\t" << nBadComp[1][i][3];
1106 ss << "\n";
1107 for (int i = 1; i < 7; ++i)
1108 ss << "\nTOB Layer " << i << " :\t\t" << nBadComp[2][i][0] << "\t" << nBadComp[2][i][1] << "\t" << nBadComp[2][i][2]
1109 << "\t" << nBadComp[2][i][3];
1110 ss << "\n";
1111 for (int i = 1; i < 10; ++i)
1112 ss << "\nTEC+ Disk " << i << " :\t\t" << nBadComp[3][i][0] << "\t" << nBadComp[3][i][1] << "\t" << nBadComp[3][i][2]
1113 << "\t" << nBadComp[3][i][3];
1114 for (int i = 10; i < 19; ++i)
1115 ss << "\nTEC- Disk " << i - 9 << " :\t\t" << nBadComp[3][i][0] << "\t" << nBadComp[3][i][1] << "\t"
1116 << nBadComp[3][i][2] << "\t" << nBadComp[3][i][3];
1117 ss << "\n";
1118
1119 ss << "\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
1120 "Apvs\n----------------------------------------------------------------";
1121 for (int i = 1; i < 5; ++i)
1122 ss << "\nTIB Layer " << i << " :" << ssV[0][i].str();
1123 ss << "\n";
1124 for (int i = 1; i < 4; ++i)
1125 ss << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
1126 for (int i = 4; i < 7; ++i)
1127 ss << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
1128 ss << "\n";
1129 for (int i = 1; i < 7; ++i)
1130 ss << "\nTOB Layer " << i << " :" << ssV[2][i].str();
1131 ss << "\n";
1132 for (int i = 1; i < 10; ++i)
1133 ss << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
1134 for (int i = 10; i < 19; ++i)
1135 ss << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
1136
1137 LOGPRINT << ss.str();
1138
1139
1140 std::ofstream badModules;
1141 badModules.open("BadModules_NEW.log");
1142 badModules << "\n----------------------------------------------------------------\n\t\t Detid \tModules Fibers "
1143 "Apvs\n----------------------------------------------------------------";
1144 for (int i = 1; i < 5; ++i)
1145 badModules << "\nTIB Layer " << i << " :" << ssV[0][i].str();
1146 badModules << "\n";
1147 for (int i = 1; i < 4; ++i)
1148 badModules << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
1149 for (int i = 4; i < 7; ++i)
1150 badModules << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
1151 badModules << "\n";
1152 for (int i = 1; i < 7; ++i)
1153 badModules << "\nTOB Layer " << i << " :" << ssV[2][i].str();
1154 badModules << "\n";
1155 for (int i = 1; i < 10; ++i)
1156 badModules << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
1157 for (int i = 10; i < 19; ++i)
1158 badModules << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
1159 badModules.close();
1160 }
1161
1162 void SiStripHitEfficiencyHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1163 edm::ParameterSetDescription desc;
1164 desc.add<std::string>("inputFolder", "AlCaReco/SiStripHitEfficiency");
1165 desc.add<bool>("isAtPCL", false);
1166 desc.add<bool>("doStoreOnDB", false);
1167 desc.add<std::string>("Record", "SiStripBadStrip");
1168 desc.add<double>("Threshold", 0.1);
1169 desc.add<std::string>("Title", "Hit Efficiency");
1170 desc.add<int>("nModsMin", 5);
1171 desc.addUntracked<bool>("doStoreOnTree", false);
1172 desc.addUntracked<bool>("AutoIneffModTagging", false);
1173 desc.addUntracked<double>("TkMapMin", 0.9);
1174 desc.addUntracked<double>("EffPlotMin", 0.9);
1175 desc.addUntracked<bool>("ShowRings", false);
1176 desc.addUntracked<bool>("ShowEndcapSides", true);
1177 desc.addUntracked<bool>("ShowTOB6TEC9", false);
1178 desc.addUntracked<bool>("ShowOnlyGoodModules", false);
1179 descriptions.addWithDefaultLabel(desc);
1180 }
1181
1182 #include "FWCore/Framework/interface/MakerMacros.h"
1183 DEFINE_FWK_MODULE(SiStripHitEfficiencyHarvester);