Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:35:10 

0001 // user includes
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 //system includes
0026 #include <boost/type_index.hpp>
0027 #include <fmt/printf.h>
0028 #include <numeric>  // for std::accumulate
0029 #include <sstream>
0030 
0031 // ROOT includes
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 // custom made printout
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   // information for the TTree
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); /* to check if TK was ON */
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),  // number of rings or wheels
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   // zero in all counts
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   // check on the input TkHistoMap
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       // store information per DetId in the output tree
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   // collect how many layers are missing
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   // check on the input TkHistoMaps
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   // no input TkHistoMaps -> early return
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   // count how many hits in the denominator we have
0241   const unsigned int totalHits = this->countTotalHits(totalMaps);
0242 
0243   // set colz
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   // come back to the main folder
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   // Effiency calculation, bad module tagging, and tracker maps //
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   // load the FEDError map
0276   const auto& EventStats = getter.get(fmt::format("{}/EventInfo/EventStats", inputFolder_));
0277   const int totalEvents = EventStats->getBinContent(1., 1.);  // first bin contains info on number of events run
0278   calibData_.FEDErrorOccupancy = std::make_unique<TkHistoMap>(tkDetMap_.get());
0279   calibData_.FEDErrorOccupancy->loadTkHistoMap(fmt::format("{}/FEDErrorTkDetMaps", inputFolder_),
0280                                                "perModule_FEDErrors");
0281 
0282   // tag as bad from FEDErrors the modules that have an error on 75% of the events
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       // use only the "good" modules
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             // We have a bad module, put it in the list!
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             //Fill the bad list with empty results for every module
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           }  // if storing tree
0328         }  // if not autoInefModTagging
0329       }  // if there are no bad APVs
0330     }  // if denom
0331   }  // loop on DetIds
0332 
0333   if (autoIneffModTagging_) {
0334     for (unsigned int i = 1; i <= k_LayersAtTECEnd; i++) {
0335       //Compute threshold to use for each layer
0336       hEffInLayer[i]->getTH1()->GetXaxis()->SetRange(
0337           3, hEffInLayer[i]->getNbinsX() + 1);  // Remove from the avg modules below 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);  // can't have this happen
0351             // use only the "good" modules
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                 //We have a bad module, put it in the list!
0358                 badModules[det] = eff;
0359                 tkMapBad.fillc(det, 255, 0, 0);
0360                 if (!isAtPCL_ && doStoreOnTree_) {
0361                   t_isTaggedIneff = true;
0362                 }
0363               } else {
0364                 //Fill the bad list with empty results for every module
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                 // printing message also for modules sligthly above (8%) the limit
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               }  // if storing tree
0389 
0390               //Put modules into the TKMap
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               // for the summary
0400               // Have to do the decoding for which side to go on (ugh)
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             }  // if the module is good!
0422 
0423             //Do the one where we don't exclude bad modules!
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           }  // if denom
0445         }  // layer = i
0446       }  // loop on detids
0447     }  // loop on layers
0448   }  // if auto tagging
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   //This is the list of the bad strips, use to mask out entire APVs
0460   //Now simply go through the bad hit list and mask out things that
0461   //are bad!
0462   for (const auto it : badModules) {
0463     const auto det = it.first;
0464     std::vector<unsigned int> badStripList;
0465     //We need to figure out how many strips are in this particular module
0466     //To Mask correctly!
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     //Now compact into a single bad module
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);  // statistics by layer and subdetector
0488   //LOGPRINT << "\n-----------------\nNew IOV starting from run " << e.id().run() << " event " << e.id().event()
0489   //     << " lumiBlock " << e.luminosityBlock() << " time " << e.time().value() << "\n-----------------\n";
0490   printAndWriteBadModules(pQuality, detInfo);  // TODO
0491 
0492   // make summary plots
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   //Calculate the statistics by layer
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   // use goodlayer_total/found and alllayer_total/found, collapse side and/or ring if needed
0564   unsigned int nLayers{34};  // default
0565   if (showRings_)
0566     nLayers = 30;
0567   if (!showEndcapSides_) {
0568     if (!showRings_)
0569       nLayers = 22;
0570     else
0571       nLayers = 20;
0572   }
0573 
0574   // come back to the main folder and create a final efficiency folder
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   // first bin only to keep real data off the y axis so set to -1
0585   found_good->setBinContent(0, -1);
0586   all_good->setBinContent(0, 1);
0587 
0588   // new ROOT version: TGraph::Divide don't handle null or negative values
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;  // barrel+endcap
0601   if (!showEndcapSides_)
0602     nLayers_max = 11;  // barrel
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   // endcap - merging sides
0619   if (!showEndcapSides_) {
0620     for (unsigned int i = 11; i < 14; ++i) {  // TID disks
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) {  // TEC disks
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     // now do the profile
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     // clean the house
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     // fill all modules efficiency
0687     if (den_all > 0.) {
0688       // naive binomial errors
0689       //float eff_all = num_all / den_all;
0690       //float err_eff_all = (eff_all * (1 - eff_all)) / den_all;
0691 
0692       // use Clopper-Pearson errors
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     // fill good modules efficiency
0699     if (den_good > 0.) {
0700       // naive binomial errors
0701       //float eff_good = num_good / den_good;
0702       //float err_eff_good = (eff_good * (1 - eff_good)) / den_good;
0703 
0704       // use Clopper-Pearson errors
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   // set the histogram bin labels
0715   this->setEffBinLabels(h_eff_all->getTH1F(), h_eff_good->getTH1F(), nLayers);
0716 
0717   if (!isAtPCL_) {
0718     // if TFileService is not avaible, just go on
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   }  // if it's not run at PCL
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       // no input histograms -> continue in the loop
0873       continue;
0874     }
0875 
0876     // in order to display correct errors when taking the ratio
0877     if (!hfound->getTH1F()->GetSumw2())
0878       hfound->getTH1F()->Sumw2();
0879     if (!htotal->getTH1F()->GetSumw2())
0880       htotal->getTH1F()->Sumw2();
0881 
0882     // prevent dividing by 0
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       // fill all modules efficiency
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     // graphics adjustment
0921     effVsVariable[iLayer]->getTH1F()->SetMinimum(tkMapMin_);
0922 
0923     if (doProfiles) {
0924       // now do the profile
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   }  // loop on layers
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 }  // namespace
0975 
0976 void SiStripHitEfficiencyHarvester::printAndWriteBadModules(const SiStripQuality& quality,
0977                                                             const SiStripDetInfo& detInfo) const {
0978   ////////////////////////////////////////////////////////////////////////
0979   //try to write out what's in the quality record
0980   /////////////////////////////////////////////////////////////////////////////
0981   int nTkBadComp[4];  //k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
0982   int nBadComp[4][19][4];
0983   //legend: nBadComp[i][j][k]= SubSystem i, layer/disk/wheel j, BadModule/Fiber/Apv k
0984   //     i: 0=TIB, 1=TID, 2=TOB, 3=TEC
0985   //     k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
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     // Full Tk
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     // single subsystem
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   // single strip info
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   // printout
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   // store also bad modules in log file
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!