Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-06-24 22:51:40 

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/SiStripHitEfficiencyHelpers.h"
0006 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0007 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0008 #include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" /* for STRIPS_PER_APV*/
0009 #include "DQM/SiStripCommon/interface/TkHistoMap.h"
0010 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0011 #include "FWCore/Framework/interface/ESWatcher.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "FWCore/ParameterSet/interface/ParameterDescription.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0019 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0020 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0021 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0022 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0023 
0024 //system includes
0025 #include <sstream>
0026 #include <numeric>  // for std::accumulate
0027 
0028 // ROOT includes
0029 #include "TEfficiency.h"
0030 
0031 // custom made printout
0032 #define LOGPRINT edm::LogPrint("SiStripHitEfficiencyHarvester")
0033 
0034 class SiStripHitEfficiencyHarvester : public DQMEDHarvester {
0035 public:
0036   explicit SiStripHitEfficiencyHarvester(const edm::ParameterSet&);
0037   ~SiStripHitEfficiencyHarvester() override = default;
0038   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0039 
0040   void endRun(edm::Run const&, edm::EventSetup const&) override;
0041   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0042 
0043 private:
0044   const std::string inputFolder_;
0045   const bool isAtPCL_;
0046   const bool showRings_, autoIneffModTagging_, doStoreOnDB_;
0047   const unsigned int nTEClayers_;
0048   const double threshold_;
0049   const int nModsMin_;
0050   const double tkMapMin_;
0051   const std::string title_, record_;
0052 
0053   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0054   const edm::ESGetToken<TkDetMap, TrackerTopologyRcd> tkDetMapToken_;
0055   const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> stripQualityToken_;
0056   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0057 
0058   std::unique_ptr<TrackerTopology> tTopo_;
0059   std::unique_ptr<TkDetMap> tkDetMap_;
0060   std::unique_ptr<SiStripQuality> stripQuality_;
0061   std::vector<DetId> stripDetIds_;
0062 
0063   void writeBadStripPayload(const SiStripQuality& quality) const;
0064   void printTotalStatistics(const std::array<long, bounds::k_END_OF_LAYERS>& layerFound,
0065                             const std::array<long, bounds::k_END_OF_LAYERS>& layerTotal) const;
0066   void printAndWriteBadModules(const SiStripQuality& quality, const SiStripDetInfo& detInfo) const;
0067   bool checkMapsValidity(const std::vector<MonitorElement*>& maps, const std::string& type) const;
0068   unsigned int countTotalHits(const std::vector<MonitorElement*>& maps); /* to check if TK was ON */
0069   void makeSummary(DQMStore::IGetter& getter, TFileService& fs) const;
0070   void makeSummaryVsBX(DQMStore::IGetter& getter, TFileService& fs) const;
0071   void makeSummaryVsLumi(DQMStore::IGetter& getter) const;
0072   void makeSummaryVsCM(DQMStore::IGetter& getter, TFileService& fs) const;
0073 };
0074 
0075 SiStripHitEfficiencyHarvester::SiStripHitEfficiencyHarvester(const edm::ParameterSet& conf)
0076     : inputFolder_(conf.getParameter<std::string>("inputFolder")),
0077       isAtPCL_(conf.getParameter<bool>("isAtPCL")),
0078       showRings_(conf.getUntrackedParameter<bool>("ShowRings", false)),
0079       autoIneffModTagging_(conf.getUntrackedParameter<bool>("AutoIneffModTagging", false)),
0080       doStoreOnDB_(conf.getParameter<bool>("doStoreOnDB")),
0081       nTEClayers_(showRings_ ? 7 : 9),  // number of rings or wheels
0082       threshold_(conf.getParameter<double>("Threshold")),
0083       nModsMin_(conf.getParameter<int>("nModsMin")),
0084       tkMapMin_(conf.getUntrackedParameter<double>("TkMapMin", 0.9)),
0085       title_(conf.getParameter<std::string>("Title")),
0086       record_(conf.getParameter<std::string>("Record")),
0087       tTopoToken_(esConsumes<edm::Transition::EndRun>()),
0088       tkDetMapToken_(esConsumes<edm::Transition::EndRun>()),
0089       stripQualityToken_(esConsumes<edm::Transition::EndRun>()),
0090       tkGeomToken_(esConsumes<edm::Transition::EndRun>()) {}
0091 
0092 void SiStripHitEfficiencyHarvester::endRun(edm::Run const&, edm::EventSetup const& iSetup) {
0093   if (!tTopo_) {
0094     tTopo_ = std::make_unique<TrackerTopology>(iSetup.getData(tTopoToken_));
0095   }
0096   if (!tkDetMap_) {
0097     tkDetMap_ = std::make_unique<TkDetMap>(iSetup.getData(tkDetMapToken_));
0098   }
0099   if (!stripQuality_) {
0100     stripQuality_ = std::make_unique<SiStripQuality>(iSetup.getData(stripQualityToken_));
0101   }
0102   if (stripDetIds_.empty()) {
0103     const auto& tkGeom = iSetup.getData(tkGeomToken_);
0104     for (const auto& det : tkGeom.detUnits()) {
0105       if (dynamic_cast<const StripGeomDetUnit*>(det)) {
0106         stripDetIds_.push_back(det->geographicalId());
0107       }
0108     }
0109   }
0110 }
0111 
0112 bool SiStripHitEfficiencyHarvester::checkMapsValidity(const std::vector<MonitorElement*>& maps,
0113                                                       const std::string& type) const {
0114   std::vector<bool> isAvailable;
0115   isAvailable.reserve(maps.size());
0116   std::transform(
0117       maps.begin() + 1, maps.end(), std::back_inserter(isAvailable), [](auto& x) { return !(x == nullptr); });
0118 
0119   int count{0};
0120   for (const auto& it : isAvailable) {
0121     count++;
0122     LogDebug("SiStripHitEfficiencyHarvester") << " layer: " << count << " " << it << std::endl;
0123     if (it)
0124       LogDebug("SiStripHitEfficiencyHarvester") << "resolving to " << maps[count]->getName() << std::endl;
0125   }
0126 
0127   // check on the input TkHistoMap
0128   bool areMapsAvailable{true};
0129   int layerCount{0};
0130   for (const auto& it : isAvailable) {
0131     layerCount++;
0132     if (!it) {
0133       edm::LogError("SiStripHitEfficiencyHarvester")
0134           << type << " TkHistoMap for layer " << layerCount << " was not found.\n -> Aborting!";
0135       areMapsAvailable = false;
0136       break;
0137     }
0138   }
0139   return areMapsAvailable;
0140 }
0141 
0142 unsigned int SiStripHitEfficiencyHarvester::countTotalHits(const std::vector<MonitorElement*>& maps) {
0143   return std::accumulate(maps.begin() + 1, maps.end(), 0, [](unsigned int total, MonitorElement* item) {
0144     return total + item->getEntries();
0145   });
0146 }
0147 
0148 void SiStripHitEfficiencyHarvester::dqmEndJob(DQMStore::IBooker& booker, DQMStore::IGetter& getter) {
0149   if (!autoIneffModTagging_)
0150     LOGPRINT << "A module is bad if efficiency < " << threshold_ << " and has at least " << nModsMin_ << " nModsMin.";
0151   else
0152     LOGPRINT << "A module is bad if the upper limit on the efficiency is < to the avg in the layer - " << threshold_
0153              << " and has at least " << nModsMin_ << " nModsMin.";
0154 
0155   auto h_module_total = std::make_unique<TkHistoMap>(tkDetMap_.get());
0156   h_module_total->loadTkHistoMap(fmt::format("{}/TkDetMaps", inputFolder_), "perModule_total");
0157   auto h_module_found = std::make_unique<TkHistoMap>(tkDetMap_.get());
0158   h_module_found->loadTkHistoMap(fmt::format("{}/TkDetMaps", inputFolder_), "perModule_found");
0159 
0160   // collect how many layers are missing
0161   const auto& totalMaps = h_module_total->getAllMaps();
0162   const auto& foundMaps = h_module_found->getAllMaps();
0163 
0164   LogDebug("SiStripHitEfficiencyHarvester")
0165       << "totalMaps.size(): " << totalMaps.size() << " foundMaps.size() " << foundMaps.size() << std::endl;
0166 
0167   // check on the input TkHistoMaps
0168   bool isTotalMapAvailable = this->checkMapsValidity(totalMaps, std::string("Total"));
0169   bool isFoundMapAvailable = this->checkMapsValidity(foundMaps, std::string("Found"));
0170 
0171   LogDebug("SiStripHitEfficiencyHarvester")
0172       << "isTotalMapAvailable: " << isTotalMapAvailable << " isFoundMapAvailable " << isFoundMapAvailable << std::endl;
0173 
0174   // no input TkHistoMaps -> early return
0175   if (!isTotalMapAvailable or !isFoundMapAvailable)
0176     return;
0177 
0178   LogDebug("SiStripHitEfficiencyHarvester")
0179       << "Entries in total TkHistoMap for layer 3: " << h_module_total->getMap(3)->getEntries() << ", found "
0180       << h_module_found->getMap(3)->getEntries();
0181 
0182   // count how many hits in the denominator we have
0183   const unsigned int totalHits = this->countTotalHits(totalMaps);
0184 
0185   // come back to the main folder
0186   booker.setCurrentFolder(inputFolder_);
0187 
0188   std::vector<MonitorElement*> hEffInLayer(std::size_t(1), nullptr);
0189   hEffInLayer.reserve(bounds::k_END_OF_LAYERS);
0190   for (std::size_t i = 1; i != bounds::k_END_OF_LAYERS; ++i) {
0191     const auto lyrName = ::layerName(i, showRings_, nTEClayers_);
0192     hEffInLayer.push_back(booker.book1D(
0193         Form("eff_layer%i", int(i)), Form("Module efficiency in layer %s", lyrName.c_str()), 201, 0, 1.005));
0194   }
0195   std::array<long, bounds::k_END_OF_LAYERS> layerTotal{};
0196   std::array<long, bounds::k_END_OF_LAYERS> layerFound{};
0197   layerTotal.fill(0);
0198   layerFound.fill(0);
0199 
0200   ////////////////////////////////////////////////////////////////
0201   // Effiency calculation, bad module tagging, and tracker maps //
0202   ////////////////////////////////////////////////////////////////
0203 
0204   TrackerMap tkMap{"  Detector Inefficiency  "};
0205   TrackerMap tkMapBad{"  Inefficient Modules  "};
0206   TrackerMap tkMapEff{title_};
0207   TrackerMap tkMapNum{" Detector numerator   "};
0208   TrackerMap tkMapDen{" Detector denominator "};
0209   std::map<unsigned int, double> badModules;
0210 
0211   for (auto det : stripDetIds_) {
0212     auto layer = ::checkLayer(det, tTopo_.get());
0213     const auto num = h_module_found->getValue(det);
0214     const auto denom = h_module_total->getValue(det);
0215     if (denom) {
0216       const auto eff = num / denom;
0217       hEffInLayer[layer]->Fill(eff);
0218       if (!autoIneffModTagging_) {
0219         if ((denom >= nModsMin_) && (eff < threshold_)) {
0220           // We have a bad module, put it in the list!
0221           badModules[det] = eff;
0222           tkMapBad.fillc(det, 255, 0, 0);
0223           LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ")  module "
0224                    << det.rawId() << " efficiency: " << eff << " , " << num << "/" << denom;
0225         } else {
0226           //Fill the bad list with empty results for every module
0227           tkMapBad.fillc(det, 255, 255, 255);
0228         }
0229         if (eff < threshold_)
0230           LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ")  module "
0231                    << det.rawId() << " efficiency: " << eff << " , " << num << "/" << denom;
0232 
0233         if (denom < nModsMin_) {
0234           LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ")  module "
0235                    << det.rawId() << " is under occupancy at " << denom;
0236         }
0237       }
0238       //Put any module into the TKMap
0239       tkMap.fill(det, 1. - eff);
0240       tkMapEff.fill(det, eff);
0241       tkMapNum.fill(det, num);
0242       tkMapDen.fill(det, denom);
0243 
0244       layerTotal[layer] += denom;
0245       layerFound[layer] += num;
0246     }
0247   }
0248 
0249   if (autoIneffModTagging_) {
0250     for (Long_t i = 1; i <= k_LayersAtTECEnd; i++) {
0251       //Compute threshold to use for each layer
0252       hEffInLayer[i]->getTH1()->GetXaxis()->SetRange(
0253           3, hEffInLayer[i]->getNbinsX() + 1);  // Remove from the avg modules below 1%
0254       const double layer_min_eff = hEffInLayer[i]->getMean() - std::max(2.5 * hEffInLayer[i]->getRMS(), threshold_);
0255       LOGPRINT << "Layer " << i << " threshold for bad modules: <" << layer_min_eff
0256                << "  (layer mean: " << hEffInLayer[i]->getMean() << " rms: " << hEffInLayer[i]->getRMS() << ")";
0257 
0258       hEffInLayer[i]->getTH1()->GetXaxis()->SetRange(1, hEffInLayer[i]->getNbinsX() + 1);
0259 
0260       for (auto det : stripDetIds_) {
0261         const auto layer = ::checkLayer(det, tTopo_.get());
0262         if (layer == i) {
0263           const auto num = h_module_found->getValue(det);
0264           const auto denom = h_module_total->getValue(det);
0265           if (denom) {
0266             const auto eff = num / denom;
0267             const auto eff_up = TEfficiency::Bayesian(denom, num, .99, 1, 1, true);
0268 
0269             if ((denom >= nModsMin_) && (eff_up < layer_min_eff)) {
0270               //We have a bad module, put it in the list!
0271               badModules[det] = eff;
0272               tkMapBad.fillc(det, 255, 0, 0);
0273             } else {
0274               //Fill the bad list with empty results for every module
0275               tkMapBad.fillc(det, 255, 255, 255);
0276             }
0277             if (eff_up < layer_min_eff + 0.08)  // printing message also for modules sligthly above (8%) the limit
0278 
0279               LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ")  module "
0280                        << det.rawId() << " efficiency: " << eff << " , " << num << "/" << denom
0281                        << " , upper limit: " << eff_up;
0282             if (denom < nModsMin_) {
0283               LOGPRINT << "Layer " << layer << " (" << ::layerName(layer, showRings_, nTEClayers_) << ")  module "
0284                        << det.rawId() << " layer " << layer << " is under occupancy at " << denom;
0285             }
0286           }
0287         }
0288       }
0289     }
0290   }
0291 
0292   tkMap.save(true, 0, 0, "SiStripHitEffTKMap_NEW.png");
0293   tkMapBad.save(true, 0, 0, "SiStripHitEffTKMapBad_NEW.png");
0294   tkMapEff.save(true, tkMapMin_, 1., "SiStripHitEffTKMapEff_NEW.png");
0295   tkMapNum.save(true, 0, 0, "SiStripHitEffTKMapNum_NEW.png");
0296   tkMapDen.save(true, 0, 0, "SiStripHitEffTKMapDen_NEW.png");
0297 
0298   const auto detInfo =
0299       SiStripDetInfoFileReader::read(edm::FileInPath{SiStripDetInfoFileReader::kDefaultFile}.fullPath());
0300   SiStripQuality pQuality{detInfo};
0301   //This is the list of the bad strips, use to mask out entire APVs
0302   //Now simply go through the bad hit list and mask out things that
0303   //are bad!
0304   for (const auto it : badModules) {
0305     const auto det = it.first;
0306     std::vector<unsigned int> badStripList;
0307     //We need to figure out how many strips are in this particular module
0308     //To Mask correctly!
0309     const auto nStrips = detInfo.getNumberOfApvsAndStripLength(det).first * sistrip::STRIPS_PER_APV;
0310     LOGPRINT << "Number of strips module " << det << " is " << nStrips;
0311     badStripList.push_back(pQuality.encode(0, nStrips, 0));
0312     //Now compact into a single bad module
0313     LOGPRINT << "ID1 shoudl match list of modules above " << det;
0314     pQuality.compact(det, badStripList);
0315     pQuality.put(det, SiStripQuality::Range(badStripList.begin(), badStripList.end()));
0316   }
0317   pQuality.fillBadComponents();
0318   if (doStoreOnDB_) {
0319     if (totalHits > 0u) {
0320       writeBadStripPayload(pQuality);
0321     } else {
0322       edm::LogPrint("SiStripHitEfficiencyHarvester")
0323           << __PRETTY_FUNCTION__ << " There are no SiStrip hits for a valid measurement, skipping!";
0324     }
0325   } else {
0326     edm::LogInfo("SiStripHitEfficiencyHarvester") << "Will not produce payload!";
0327   }
0328 
0329   printTotalStatistics(layerFound, layerTotal);  // statistics by layer and subdetector
0330   //LOGPRINT << "\n-----------------\nNew IOV starting from run " << e.id().run() << " event " << e.id().event()
0331   //     << " lumiBlock " << e.luminosityBlock() << " time " << e.time().value() << "\n-----------------\n";
0332   printAndWriteBadModules(pQuality, detInfo);  // TODO
0333 
0334   if (!isAtPCL_) {
0335     edm::Service<TFileService> fs;
0336     makeSummary(getter, *fs);      // TODO
0337     makeSummaryVsBX(getter, *fs);  // TODO
0338     makeSummaryVsCM(getter, *fs);  // TODO
0339   }
0340 
0341   makeSummaryVsLumi(getter);  // TODO
0342 }
0343 
0344 void SiStripHitEfficiencyHarvester::printTotalStatistics(
0345     const std::array<long, bounds::k_END_OF_LAYERS>& layerFound,
0346     const std::array<long, bounds::k_END_OF_LAYERS>& layerTotal) const {
0347   //Calculate the statistics by layer
0348   int totalfound = 0;
0349   int totaltotal = 0;
0350   double layereff;
0351   int subdetfound[5];
0352   int subdettotal[5];
0353 
0354   for (Long_t i = 1; i < 5; i++) {
0355     subdetfound[i] = 0;
0356     subdettotal[i] = 0;
0357   }
0358 
0359   for (Long_t i = 1; i <= bounds::k_LayersAtTECEnd; i++) {
0360     layereff = double(layerFound[i]) / double(layerTotal[i]);
0361     LOGPRINT << "Layer " << i << " (" << ::layerName(i, showRings_, nTEClayers_) << ") has total efficiency "
0362              << layereff << " " << layerFound[i] << "/" << layerTotal[i];
0363     totalfound += layerFound[i];
0364     totaltotal += layerTotal[i];
0365     if (i <= bounds::k_LayersAtTIBEnd) {
0366       subdetfound[1] += layerFound[i];
0367       subdettotal[1] += layerTotal[i];
0368     }
0369     if (i > bounds::k_LayersAtTIBEnd && i <= bounds::k_LayersAtTOBEnd) {
0370       subdetfound[2] += layerFound[i];
0371       subdettotal[2] += layerTotal[i];
0372     }
0373     if (i > bounds::k_LayersAtTOBEnd && i <= bounds::k_LayersAtTIDEnd) {
0374       subdetfound[3] += layerFound[i];
0375       subdettotal[3] += layerTotal[i];
0376     }
0377     if (i > bounds::k_LayersAtTIDEnd) {
0378       subdetfound[4] += layerFound[i];
0379       subdettotal[4] += layerTotal[i];
0380     }
0381   }
0382 
0383   LOGPRINT << "The total efficiency is " << double(totalfound) / double(totaltotal);
0384   LOGPRINT << "      TIB: " << double(subdetfound[1]) / subdettotal[1] << " " << subdetfound[1] << "/"
0385            << subdettotal[1];
0386   LOGPRINT << "      TOB: " << double(subdetfound[2]) / subdettotal[2] << " " << subdetfound[2] << "/"
0387            << subdettotal[2];
0388   LOGPRINT << "      TID: " << double(subdetfound[3]) / subdettotal[3] << " " << subdetfound[3] << "/"
0389            << subdettotal[3];
0390   LOGPRINT << "      TEC: " << double(subdetfound[4]) / subdettotal[4] << " " << subdetfound[4] << "/"
0391            << subdettotal[4];
0392 }
0393 
0394 void SiStripHitEfficiencyHarvester::writeBadStripPayload(const SiStripQuality& quality) const {
0395   SiStripBadStrip pBadStrip{};
0396   const auto pQdvBegin = quality.getDataVectorBegin();
0397   for (auto rIt = quality.getRegistryVectorBegin(); rIt != quality.getRegistryVectorEnd(); ++rIt) {
0398     const auto range = SiStripBadStrip::Range(pQdvBegin + rIt->ibegin, pQdvBegin + rIt->iend);
0399     if (!pBadStrip.put(rIt->detid, range))
0400       edm::LogError("SiStripHitEfficiencyHarvester") << "detid already exists in SiStripBadStrip";
0401   }
0402   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0403   if (poolDbService.isAvailable()) {
0404     poolDbService->writeOneIOV(pBadStrip, poolDbService->currentTime(), record_);
0405   } else {
0406     throw cms::Exception("PoolDBService required");
0407   }
0408 }
0409 
0410 void SiStripHitEfficiencyHarvester::makeSummary(DQMStore::IGetter& getter, TFileService& fs) const {
0411   // use goodlayer_total/found and alllayer_total/found, collapse side and/or ring if needed
0412 }
0413 
0414 void SiStripHitEfficiencyHarvester::makeSummaryVsBX(DQMStore::IGetter& getter, TFileService& fs) const {
0415   // use found/totalVsBx_layer%i [0,23)
0416 }
0417 
0418 void SiStripHitEfficiencyHarvester::makeSummaryVsLumi(DQMStore::IGetter& getter) const {
0419   for (unsigned int iLayer = 1; iLayer != (showRings_ ? 20 : 22); ++iLayer) {
0420     auto hfound = getter.get(fmt::format("{}/VsLumi/layerfound_vsLumi_layer_{}", inputFolder_, iLayer))->getTH1F();
0421     auto htotal = getter.get(fmt::format("{}/VsLumi/layertotal_vsLumi_layer_{}", inputFolder_, iLayer))->getTH1F();
0422 
0423     if (hfound == nullptr or htotal == nullptr) {
0424       if (hfound == nullptr)
0425         edm::LogError("SiStripHitEfficiencyHarvester")
0426             << fmt::format("{}/VsLumi/layerfound_vsLumi_layer_{}", inputFolder_, iLayer) << " was not found!";
0427       if (htotal == nullptr)
0428         edm::LogError("SiStripHitEfficiencyHarvester")
0429             << fmt::format("{}/VsLumi/layertotal_vsLumi_layer_{}", inputFolder_, iLayer) << " was not found!";
0430       // no input histograms -> continue in the loop
0431       continue;
0432     }
0433 
0434     if (!hfound->GetSumw2())
0435       hfound->Sumw2();
0436     if (!htotal->GetSumw2())
0437       htotal->Sumw2();
0438     for (Long_t i = 0; i != hfound->GetNbinsX() + 1; ++i) {
0439       if (hfound->GetBinContent(i) == 0)
0440         hfound->SetBinContent(i, 1e-6);
0441       if (htotal->GetBinContent(i) == 0)
0442         htotal->SetBinContent(i, 1);
0443     }
0444     LogDebug("SiStripHitEfficiencyHarvester")
0445         << "Total hits for layer " << iLayer << " (vs lumi): " << htotal->GetEntries() << ", found "
0446         << hfound->GetEntries();
0447   }
0448   // continue
0449 }
0450 
0451 void SiStripHitEfficiencyHarvester::makeSummaryVsCM(DQMStore::IGetter& getter, TFileService& fs) const {}
0452 
0453 namespace {
0454   void setBadComponents(int i,
0455                         int comp,
0456                         const SiStripQuality::BadComponent& bc,
0457                         std::stringstream ssV[4][19],
0458                         int nBad[4][19][4],
0459                         int nAPV) {
0460     ssV[i][comp] << "\n\t\t " << bc.detid << " \t " << bc.BadModule << " \t " << ((bc.BadFibers) & 0x1) << " ";
0461     if (nAPV == 4)
0462       ssV[i][comp] << "x " << ((bc.BadFibers >> 1) & 0x1);
0463 
0464     if (nAPV == 6)
0465       ssV[i][comp] << ((bc.BadFibers >> 1) & 0x1) << " " << ((bc.BadFibers >> 2) & 0x1);
0466     ssV[i][comp] << " \t " << ((bc.BadApvs) & 0x1) << " " << ((bc.BadApvs >> 1) & 0x1) << " ";
0467     if (nAPV == 4)
0468       ssV[i][comp] << "x x " << ((bc.BadApvs >> 2) & 0x1) << " " << ((bc.BadApvs >> 3) & 0x1);
0469     if (nAPV == 6)
0470       ssV[i][comp] << ((bc.BadApvs >> 2) & 0x1) << " " << ((bc.BadApvs >> 3) & 0x1) << " " << ((bc.BadApvs >> 4) & 0x1)
0471                    << " " << ((bc.BadApvs >> 5) & 0x1) << " ";
0472 
0473     if (bc.BadApvs) {
0474       nBad[i][0][2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
0475                        ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
0476       nBad[i][comp][2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
0477                           ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
0478     }
0479     if (bc.BadFibers) {
0480       nBad[i][0][1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
0481       nBad[i][comp][1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
0482     }
0483     if (bc.BadModule) {
0484       nBad[i][0][0]++;
0485       nBad[i][comp][0]++;
0486     }
0487   }
0488 }  // namespace
0489 
0490 void SiStripHitEfficiencyHarvester::printAndWriteBadModules(const SiStripQuality& quality,
0491                                                             const SiStripDetInfo& detInfo) const {
0492   ////////////////////////////////////////////////////////////////////////
0493   //try to write out what's in the quality record
0494   /////////////////////////////////////////////////////////////////////////////
0495   int nTkBadComp[4];  //k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
0496   int nBadComp[4][19][4];
0497   //legend: nBadComp[i][j][k]= SubSystem i, layer/disk/wheel j, BadModule/Fiber/Apv k
0498   //     i: 0=TIB, 1=TID, 2=TOB, 3=TEC
0499   //     k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips
0500   std::stringstream ssV[4][19];
0501 
0502   for (int i = 0; i < 4; ++i) {
0503     nTkBadComp[i] = 0;
0504     for (int j = 0; j < 19; ++j) {
0505       ssV[i][j].str("");
0506       for (int k = 0; k < 4; ++k)
0507         nBadComp[i][j][k] = 0;
0508     }
0509   }
0510 
0511   for (const auto& bc : quality.getBadComponentList()) {
0512     // Full Tk
0513     if (bc.BadModule)
0514       nTkBadComp[0]++;
0515     if (bc.BadFibers)
0516       nTkBadComp[1] += ((bc.BadFibers >> 2) & 0x1) + ((bc.BadFibers >> 1) & 0x1) + ((bc.BadFibers) & 0x1);
0517     if (bc.BadApvs)
0518       nTkBadComp[2] += ((bc.BadApvs >> 5) & 0x1) + ((bc.BadApvs >> 4) & 0x1) + ((bc.BadApvs >> 3) & 0x1) +
0519                        ((bc.BadApvs >> 2) & 0x1) + ((bc.BadApvs >> 1) & 0x1) + ((bc.BadApvs) & 0x1);
0520     // single subsystem
0521     DetId det(bc.detid);
0522     if ((det.subdetId() >= SiStripSubdetector::TIB) && (det.subdetId() <= SiStripSubdetector::TEC)) {
0523       const auto nAPV = detInfo.getNumberOfApvsAndStripLength(det).first;
0524       switch (det.subdetId()) {
0525         case SiStripSubdetector::TIB:
0526           setBadComponents(0, tTopo_->tibLayer(det), bc, ssV, nBadComp, nAPV);
0527           break;
0528         case SiStripSubdetector::TID:
0529           setBadComponents(1,
0530                            (tTopo_->tidSide(det) == 2 ? tTopo_->tidWheel(det) : tTopo_->tidWheel(det) + 3),
0531                            bc,
0532                            ssV,
0533                            nBadComp,
0534                            nAPV);
0535           break;
0536         case SiStripSubdetector::TOB:
0537           setBadComponents(2, tTopo_->tobLayer(det), bc, ssV, nBadComp, nAPV);
0538           break;
0539         case SiStripSubdetector::TEC:
0540           setBadComponents(3,
0541                            (tTopo_->tecSide(det) == 2 ? tTopo_->tecWheel(det) : tTopo_->tecWheel(det) + 9),
0542                            bc,
0543                            ssV,
0544                            nBadComp,
0545                            nAPV);
0546           break;
0547         default:
0548           break;
0549       }
0550     }
0551   }
0552   // single strip info
0553   for (auto rp = quality.getRegistryVectorBegin(); rp != quality.getRegistryVectorEnd(); ++rp) {
0554     DetId det{rp->detid};
0555     int subdet = -999;
0556     int component = -999;
0557     switch (det.subdetId()) {
0558       case SiStripSubdetector::TIB:
0559         subdet = 0;
0560         component = tTopo_->tibLayer(det);
0561         break;
0562       case SiStripSubdetector::TID:
0563         subdet = 1;
0564         component = tTopo_->tidSide(det) == 2 ? tTopo_->tidWheel(det) : tTopo_->tidWheel(det) + 3;
0565         break;
0566       case SiStripSubdetector::TOB:
0567         subdet = 2;
0568         component = tTopo_->tobLayer(det);
0569         break;
0570       case SiStripSubdetector::TEC:
0571         subdet = 3;
0572         component = tTopo_->tecSide(det) == 2 ? tTopo_->tecWheel(det) : tTopo_->tecWheel(det) + 9;
0573         break;
0574       default:
0575         break;
0576     }
0577 
0578     const auto pQdvBegin = quality.getDataVectorBegin();
0579     const auto sqrange = SiStripQuality::Range(pQdvBegin + rp->ibegin, pQdvBegin + rp->iend);
0580     float percentage = 0;
0581     for (int it = 0; it < sqrange.second - sqrange.first; it++) {
0582       unsigned int range = quality.decode(*(sqrange.first + it)).range;
0583       nTkBadComp[3] += range;
0584       nBadComp[subdet][0][3] += range;
0585       nBadComp[subdet][component][3] += range;
0586       percentage += range;
0587     }
0588     if (percentage != 0)
0589       percentage /= (sistrip::STRIPS_PER_APV * detInfo.getNumberOfApvsAndStripLength(det).first);
0590     if (percentage > 1)
0591       edm::LogError("SiStripHitEfficiencyHarvester") << "PROBLEM detid " << det.rawId() << " value " << percentage;
0592   }
0593 
0594   // printout
0595   std::ostringstream ss;
0596   ss << "\n-----------------\nGlobal Info\n-----------------";
0597   ss << "\nBadComp \t   Modules \tFibers "
0598         "\tApvs\tStrips\n----------------------------------------------------------------";
0599   ss << "\nTracker:\t\t" << nTkBadComp[0] << "\t" << nTkBadComp[1] << "\t" << nTkBadComp[2] << "\t" << nTkBadComp[3];
0600   ss << "\nTIB:\t\t\t" << nBadComp[0][0][0] << "\t" << nBadComp[0][0][1] << "\t" << nBadComp[0][0][2] << "\t"
0601      << nBadComp[0][0][3];
0602   ss << "\nTID:\t\t\t" << nBadComp[1][0][0] << "\t" << nBadComp[1][0][1] << "\t" << nBadComp[1][0][2] << "\t"
0603      << nBadComp[1][0][3];
0604   ss << "\nTOB:\t\t\t" << nBadComp[2][0][0] << "\t" << nBadComp[2][0][1] << "\t" << nBadComp[2][0][2] << "\t"
0605      << nBadComp[2][0][3];
0606   ss << "\nTEC:\t\t\t" << nBadComp[3][0][0] << "\t" << nBadComp[3][0][1] << "\t" << nBadComp[3][0][2] << "\t"
0607      << nBadComp[3][0][3];
0608   ss << "\n";
0609 
0610   for (int i = 1; i < 5; ++i)
0611     ss << "\nTIB Layer " << i << " :\t\t" << nBadComp[0][i][0] << "\t" << nBadComp[0][i][1] << "\t" << nBadComp[0][i][2]
0612        << "\t" << nBadComp[0][i][3];
0613   ss << "\n";
0614   for (int i = 1; i < 4; ++i)
0615     ss << "\nTID+ Disk " << i << " :\t\t" << nBadComp[1][i][0] << "\t" << nBadComp[1][i][1] << "\t" << nBadComp[1][i][2]
0616        << "\t" << nBadComp[1][i][3];
0617   for (int i = 4; i < 7; ++i)
0618     ss << "\nTID- Disk " << i - 3 << " :\t\t" << nBadComp[1][i][0] << "\t" << nBadComp[1][i][1] << "\t"
0619        << nBadComp[1][i][2] << "\t" << nBadComp[1][i][3];
0620   ss << "\n";
0621   for (int i = 1; i < 7; ++i)
0622     ss << "\nTOB Layer " << i << " :\t\t" << nBadComp[2][i][0] << "\t" << nBadComp[2][i][1] << "\t" << nBadComp[2][i][2]
0623        << "\t" << nBadComp[2][i][3];
0624   ss << "\n";
0625   for (int i = 1; i < 10; ++i)
0626     ss << "\nTEC+ Disk " << i << " :\t\t" << nBadComp[3][i][0] << "\t" << nBadComp[3][i][1] << "\t" << nBadComp[3][i][2]
0627        << "\t" << nBadComp[3][i][3];
0628   for (int i = 10; i < 19; ++i)
0629     ss << "\nTEC- Disk " << i - 9 << " :\t\t" << nBadComp[3][i][0] << "\t" << nBadComp[3][i][1] << "\t"
0630        << nBadComp[3][i][2] << "\t" << nBadComp[3][i][3];
0631   ss << "\n";
0632 
0633   ss << "\n----------------------------------------------------------------\n\t\t   Detid  \tModules Fibers "
0634         "Apvs\n----------------------------------------------------------------";
0635   for (int i = 1; i < 5; ++i)
0636     ss << "\nTIB Layer " << i << " :" << ssV[0][i].str();
0637   ss << "\n";
0638   for (int i = 1; i < 4; ++i)
0639     ss << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
0640   for (int i = 4; i < 7; ++i)
0641     ss << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
0642   ss << "\n";
0643   for (int i = 1; i < 7; ++i)
0644     ss << "\nTOB Layer " << i << " :" << ssV[2][i].str();
0645   ss << "\n";
0646   for (int i = 1; i < 10; ++i)
0647     ss << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
0648   for (int i = 10; i < 19; ++i)
0649     ss << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
0650 
0651   LOGPRINT << ss.str();
0652 
0653   // store also bad modules in log file
0654   std::ofstream badModules;
0655   badModules.open("BadModules_NEW.log");
0656   badModules << "\n----------------------------------------------------------------\n\t\t   Detid  \tModules Fibers "
0657                 "Apvs\n----------------------------------------------------------------";
0658   for (int i = 1; i < 5; ++i)
0659     badModules << "\nTIB Layer " << i << " :" << ssV[0][i].str();
0660   badModules << "\n";
0661   for (int i = 1; i < 4; ++i)
0662     badModules << "\nTID+ Disk " << i << " :" << ssV[1][i].str();
0663   for (int i = 4; i < 7; ++i)
0664     badModules << "\nTID- Disk " << i - 3 << " :" << ssV[1][i].str();
0665   badModules << "\n";
0666   for (int i = 1; i < 7; ++i)
0667     badModules << "\nTOB Layer " << i << " :" << ssV[2][i].str();
0668   badModules << "\n";
0669   for (int i = 1; i < 10; ++i)
0670     badModules << "\nTEC+ Disk " << i << " :" << ssV[3][i].str();
0671   for (int i = 10; i < 19; ++i)
0672     badModules << "\nTEC- Disk " << i - 9 << " :" << ssV[3][i].str();
0673   badModules.close();
0674 }
0675 
0676 void SiStripHitEfficiencyHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0677   edm::ParameterSetDescription desc;
0678   desc.add<std::string>("inputFolder", "AlCaReco/SiStripHitEfficiency");
0679   desc.add<bool>("isAtPCL", false);
0680   desc.add<bool>("doStoreOnDB", false);
0681   desc.add<std::string>("Record", "SiStripBadStrip");
0682   desc.add<double>("Threshold", 0.1);
0683   desc.add<std::string>("Title", "Hit Efficiency");
0684   desc.add<int>("nModsMin", 5);
0685   desc.addUntracked<bool>("AutoIneffModTagging", false);
0686   desc.addUntracked<double>("TkMapMin", 0.9);
0687   desc.addUntracked<bool>("ShowRings", false);
0688   descriptions.addWithDefaultLabel(desc);
0689 }
0690 
0691 #include "FWCore/Framework/interface/MakerMacros.h"
0692 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!