Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-25 23:40:02

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