Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-09 23:46:17 

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