SiStripMonitorSummary/plugins/SiStripCorrelateNoise.cc

0001 #include "DQM/SiStripMonitorSummary/plugins/SiStripCorrelateNoise.h"
0002 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0003 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0004 #include "FWCore/Framework/interface/Run.h"
0005 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0006 #include "TCanvas.h"
0007
0008 SiStripCorrelateNoise::SiStripCorrelateNoise(const edm::ParameterSet &iConfig)
0009     : noiseToken_{esConsumes<edm::Transition::BeginRun>()},
0010       gainToken_{esConsumes<edm::Transition::BeginRun>()},
0011       tTopoToken_{esConsumes<edm::Transition::BeginRun>()},
0012       refNoise(nullptr),
0013       oldGain(nullptr),
0014       newGain(nullptr) {
0015   // now do what ever initialization is needed
0016   file = new TFile("correlTest.root", "RECREATE");
0017
0018   file->cd();
0019   tkmap = new TrackerMap();
0020 }
0021
0022 void SiStripCorrelateNoise::beginRun(const edm::Run &run, const edm::EventSetup &es) {
0023   if (noiseWatcher_.check(es)) {
0024     auto aNoise = std::make_unique<SiStripNoises>(es.getData(noiseToken_));
0025     // Check if gain is the same from one noise iov to the other, otherwise cache
0026     // the new gain (and the old one) to rescale
0027     checkGainCache(es);
0028     if (refNoise) {
0029       char dir[128];
0030       theRun = run.run();
0031       sprintf(dir, "Run_%d", theRun);
0032       file->cd("");
0033       file->mkdir(dir);
0034       file->cd(dir);
0035       DoAnalysis(es, *aNoise, *refNoise);
0036       DoPlots();
0037     }
0038     refNoise = std::move(aNoise);
0039   }
0040 }
0041
0042 void SiStripCorrelateNoise::checkGainCache(const edm::EventSetup &es) {
0043   equalGain = true;
0044   if (gainWatcher_.check(es)) {
0045     if (oldGain) {
0046       equalGain = false;
0047     }
0048     oldGain = std::move(newGain);
0049     newGain = std::make_unique<SiStripApvGain>(es.getData(gainToken_));
0050   }
0051 }
0052
0053 void SiStripCorrelateNoise::DoPlots() {
0054   TCanvas *C = new TCanvas();
0055   C->Divide(2, 2);
0056
0057   char outName[128];
0058   sprintf(outName, "Run_%d.png", theRun);
0059   for (size_t i = 0; i < vTH1.size(); i++)
0060     if (vTH1[i] != nullptr) {
0061       if (i % 100 == 0) {
0062         C->cd(i / 100);
0063         vTH1[i]->SetLineColor(i / 100);
0064         vTH1[i]->Draw();
0065         C->cd(i / 100)->SetLogy();
0066       }
0067       vTH1[i]->Write();
0068     }
0069
0070   C->Print(outName);
0071   delete C;
0072
0073   vTH1.clear();
0074   file->cd("");
0075
0076   char dir[128];
0077   sprintf(dir, "Run_%d_TkMap.png", theRun);
0078   tkmap->save(false, 0, 5, dir);
0079   delete tkmap;
0080   tkmap = new TrackerMap();
0081 }
0082
0083 void SiStripCorrelateNoise::DoAnalysis(const edm::EventSetup &es,
0084                                        const SiStripNoises &_Noise,
0085                                        SiStripNoises &refNoise) {
0086   SiStripNoises Noise = _Noise;
0087   typedef std::vector<SiStripNoises::ratioData> collection;
0088   collection divNoise = Noise / refNoise;
0089
0090   edm::LogInfo("") << "[Doanalysis]";
0091
0092   const auto tTopo = &es.getData(tTopoToken_);
0093
0094   std::vector<TH1F *> histos;
0095
0096   collection::const_iterator iter = divNoise.begin();
0097   collection::const_iterator iterE = divNoise.end();
0098
0099   float value;
0100   float gainRatio = 1.;
0101   // Divide result by d
0102   for (; iter != iterE; ++iter) {
0103     getHistos(iter->detid, tTopo, histos);
0104
0105     size_t strip = 0, stripE = iter->values.size();
0106     size_t apvNb = 7;
0107
0108     for (; strip < stripE; ++strip) {
0109       if (!equalGain && strip / 128 != apvNb) {
0110         apvNb = strip / 128;
0111         if (apvNb < 6)
0112           gainRatio = getGainRatio(iter->detid, apvNb);
0113         else
0114           edm::LogInfo("") << "[Doanalysis] detid " << iter->detid << " strip " << strip << " apvNb " << apvNb;
0115       }
0116       // edm::LogInfo("") << "[Doanalysis] detid " << iter->detid << " strip "
0117       // << strip << " value " << iter->values[strip];
0118       value = iter->values[strip] * gainRatio;
0119       tkmap->fill(iter->detid, value);
0120       for (size_t i = 0; i < histos.size(); ++i)
0121         histos[i]->Fill(value);
0122     }
0123   }
0124 }
0125
0126 float SiStripCorrelateNoise::getGainRatio(const uint32_t &detid, const uint16_t &apv) {
0127   SiStripApvGain::Range oldRange = oldGain->getRange(detid);
0128   SiStripApvGain::Range newRange = newGain->getRange(detid);
0129
0130   if (oldRange.first == oldRange.second || newRange.first == newRange.second)
0131     return 1.;
0132
0133   return oldGain->getApvGain(apv, oldRange) / newGain->getApvGain(apv, newRange);
0134 }
0135
0136 void SiStripCorrelateNoise::getHistos(const uint32_t &detid,
0137                                       const TrackerTopology *tTopo,
0138                                       std::vector<TH1F *> &histos) {
0139   histos.clear();
0140
0141   int subdet = -999;
0142   int component = -999;
0143   SiStripDetId a(detid);
0144   if (a.subdetId() == 3) {
0145     subdet = 0;
0146     component = tTopo->tibLayer(detid);
0147   } else if (a.subdetId() == 4) {
0148     subdet = 1;
0149     component = tTopo->tidSide(detid) == 2 ? tTopo->tidWheel(detid) : tTopo->tidWheel(detid) + 3;
0150   } else if (a.subdetId() == 5) {
0151     subdet = 2;
0152     component = tTopo->tobLayer(detid);
0153   } else if (a.subdetId() == 6) {
0154     subdet = 3;
0155     component = tTopo->tecSide(detid) == 2 ? tTopo->tecWheel(detid) : tTopo->tecWheel(detid) + 9;
0156   }
0157
0158   int index = 100 + subdet * 100 + component;
0159
0160   histos.push_back(getHisto(100 + 100 * subdet));
0161   histos.push_back(getHisto(index));
0162 }
0163
0164 TH1F *SiStripCorrelateNoise::getHisto(const long unsigned int &index) {
0165   if (vTH1.size() < index + 1)
0166     vTH1.resize(index + 1, nullptr);
0167
0168   if (vTH1[index] == nullptr) {
0169     char name[128];
0170     std::string SubD;
0171     if (index < 200)
0172       SubD = "TIB";
0173     else if (index < 300)
0174       SubD = "TID";
0175     else if (index < 400)
0176       SubD = "TOB";
0177     else
0178       SubD = "TEC";
0179     sprintf(name, "%d_%lu__%s", theRun, index, SubD.c_str());
0180     edm::LogInfo("") << "[getHisto] creating index " << index << std::endl;
0181     vTH1[index] = new TH1F(name, name, 200, -0.5, 10.5);
0182   }
0183
0184   return vTH1[index];
0185 }
0186
0187 void SiStripCorrelateNoise::endJob() {
0188   file->Write();
0189   file->Close();
0190 }