SiStripMonitorSummary/plugins/SiStripCorrelateBadStripAndNoise.cc

0001 #include "SiStripCorrelateBadStripAndNoise.h"
0002 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0003 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0004 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0005 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0006 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0007
0008 SiStripCorrelateBadStripAndNoise::SiStripCorrelateBadStripAndNoise(const edm::ParameterSet &iConfig)
0009     : qualityToken_{esConsumes<edm::Transition::BeginRun>()},
0010       noiseToken_{esConsumes<edm::Transition::BeginRun>()},
0011       tTopoToken_{esConsumes<edm::Transition::BeginRun>()},
0012       tkGeomToken_{esConsumes<edm::Transition::BeginRun>()} {
0013   // now do what ever initialization is needed
0014   file = new TFile("correlTest.root", "RECREATE");
0015   tkmap = new TrackerMap();
0016 }
0017
0018 void SiStripCorrelateBadStripAndNoise::beginRun(const edm::Run &run, const edm::EventSetup &es) {
0019   auto newNoise = noiseWatcher_.check(es);
0020   auto newQuality = qualityWatcher_.check(es);
0021   if (newNoise || newQuality) {
0022     edm::LogInfo("") << "[SiStripCorrelateBadStripAndNoise::beginRun]" << std::endl;
0023
0024     quality_ = &es.getData(qualityToken_);
0025     noises_ = &es.getData(noiseToken_);
0026
0027     DoAnalysis(es);
0028   }
0029 }
0030
0031 void SiStripCorrelateBadStripAndNoise::DoAnalysis(const edm::EventSetup &es) {
0032   // Loop on quality bad stirps
0033   // for each strip, look at the noise
0034   // evalaute the mean apv noise and the ratio among strip noise and
0035   // meanApvNoise put the value in the histo in terms of ratio Vs percentage of
0036   // badStrips per APV
0037
0038   // Fill an histo per subdet and layer (and plus && minus for TEC/TID)
0039   edm::LogInfo("") << "[Doanalysis]";
0040   iterateOnDets(es.getData(tTopoToken_), es.getData(tkGeomToken_));
0041 }
0042
0043 void SiStripCorrelateBadStripAndNoise::iterateOnDets(const TrackerTopology &tTopo, const TrackerGeometry &tGeom) {
0044   const auto rbegin = quality_->getRegistryVectorBegin();
0045   const auto rend = quality_->getRegistryVectorEnd();
0046   for (auto rp = rbegin; rp != rend; ++rp) {
0047     const uint32_t detid = rp->detid;
0048
0049     auto sqrange =
0050         SiStripQuality::Range(quality_->getDataVectorBegin() + rp->ibegin, quality_->getDataVectorBegin() + rp->iend);
0051     iterateOnBadStrips(detid, tTopo, tGeom, sqrange);
0052   }
0053 }
0054
0055 void SiStripCorrelateBadStripAndNoise::iterateOnBadStrips(const uint32_t &detid,
0056                                                           const TrackerTopology &tTopo,
0057                                                           const TrackerGeometry &tGeom,
0058                                                           SiStripQuality::Range &sqrange) {
0059   float percentage = 0;
0060   for (int it = 0; it < sqrange.second - sqrange.first; it++) {
0061     unsigned int firstStrip = quality_->decode(*(sqrange.first + it)).firstStrip;
0062     unsigned int range = quality_->decode(*(sqrange.first + it)).range;
0063
0064     correlateWithNoise(detid, tTopo, firstStrip, range);
0065
0066     edm::LogInfo("range") << range;
0067     percentage += range;
0068   }
0069   if (percentage != 0)
0070     percentage /= dynamic_cast<const StripGeomDetUnit *>(tGeom.idToDet(detid))->specificTopology().nstrips();
0071   if (percentage > 1)
0072     edm::LogError("SiStripQualityStatistics") << "PROBLEM detid " << detid << " value " << percentage << std::endl;
0073
0074   //------- Global Statistics on percentage of bad components along the IOVs
0075   //------//
0076   if (percentage != 0)
0077     edm::LogInfo("") << "percentage " << detid << " " << percentage;
0078 }
0079
0080 void SiStripCorrelateBadStripAndNoise::correlateWithNoise(const uint32_t &detid,
0081                                                           const TrackerTopology &tTopo,
0082                                                           const uint32_t &firstStrip,
0083                                                           const uint32_t &range) {
0084   std::vector<TH2F *> histos;
0085
0086   SiStripNoises::Range noiseRange = noises_->getRange(detid);
0087   edm::LogInfo("Domenico") << "detid " << detid << " first " << firstStrip << " range " << range;
0088   float meanAPVNoise = getMeanNoise(noiseRange, firstStrip / 128, 128);
0089
0090   // float meanNoiseHotStrips=getMeanNoise(noiseRange,firstStrip,range);
0091   for (size_t theStrip = firstStrip; theStrip < firstStrip + range; theStrip++) {
0092     float meanNoiseHotStrips = getMeanNoise(noiseRange, theStrip, 1);
0093
0094     // Get the histogram for this detid
0095     getHistos(detid, tTopo, histos);
0096     float yvalue = range < 21 ? 1. * range : 21;
0097
0098     for (size_t i = 0; i < histos.size(); ++i)
0099       histos[i]->Fill(meanNoiseHotStrips / meanAPVNoise - 1., yvalue);
0100
0101     if (meanNoiseHotStrips / meanAPVNoise - 1. < -0.3)
0102       tkmap->fillc(detid, 0xFF0000);
0103     else
0104       tkmap->fillc(detid, 0x0000FF);
0105   }
0106 }
0107
0108 float SiStripCorrelateBadStripAndNoise::getMeanNoise(const SiStripNoises::Range &noiseRange,
0109                                                      const uint32_t &firstStrip,
0110                                                      const uint32_t &range) {
0111   float mean = 0;
0112   for (size_t istrip = firstStrip; istrip < firstStrip + range; istrip++) {
0113     mean += noises_->getNoise(istrip, noiseRange);
0114   }
0115   return mean / (1. * range);
0116 }
0117
0118 void SiStripCorrelateBadStripAndNoise::getHistos(const uint32_t &detid,
0119                                                  const TrackerTopology &tTopo,
0120                                                  std::vector<TH2F *> &histos) {
0121   histos.clear();
0122
0123   int subdet = -999;
0124   int component = -999;
0125   SiStripDetId a(detid);
0126   if (a.subdetId() == 3) {
0127     subdet = 0;
0128     component = tTopo.tibLayer(detid);
0129   } else if (a.subdetId() == 4) {
0130     subdet = 1;
0131     component = tTopo.tidSide(detid) == 2 ? tTopo.tidWheel(detid) : tTopo.tidWheel(detid) + 3;
0132   } else if (a.subdetId() == 5) {
0133     subdet = 2;
0134     component = tTopo.tobLayer(detid);
0135   } else if (a.subdetId() == 6) {
0136     subdet = 3;
0137     component = tTopo.tecSide(detid) == 2 ? tTopo.tecWheel(detid) : tTopo.tecWheel(detid) + 9;
0138   }
0139
0140   int index = 100 + subdet * 100 + component;
0141
0142   histos.push_back(getHisto(subdet));
0143   histos.push_back(getHisto(index));
0144 }
0145
0146 TH2F *SiStripCorrelateBadStripAndNoise::getHisto(const long unsigned int &index) {
0147   if (vTH2.size() < index + 1)
0148     vTH2.resize(index + 1, nullptr);
0149
0150   if (vTH2[index] == nullptr) {
0151     char name[128];
0152     sprintf(name, "%lu", index);
0153     edm::LogInfo("") << "[getHisto] creating index " << index << std::endl;
0154     vTH2[index] = new TH2F(name, name, 50, -2., 2., 21, 0.5, 21.5);
0155   }
0156
0157   return vTH2[index];
0158 }
0159
0160 void SiStripCorrelateBadStripAndNoise::endJob() {
0161   for (size_t i = 0; i < vTH2.size(); i++)
0162     if (vTH2[i] != nullptr)
0163       vTH2[i]->Write();
0164
0165   file->Write();
0166   file->Close();
0167
0168   tkmap->save(true, 0, 0, "testTkMap.png");
0169 }