SiStripCorrelateNoise.cc

Go to the documentation of this file.

#include "DQM/SiStripMonitorSummary/plugins/SiStripCorrelateNoise.h"

#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "FWCore/Framework/interface/Run.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "TCanvas.h"

SiStripCorrelateNoise::SiStripCorrelateNoise(const edm::ParameterSet &iConfig)
    : refNoise(nullptr), oldGain(nullptr), newGain(nullptr), cacheID_noise(0xFFFFFFFF), cacheID_gain(0xFFFFFFFF) {
  // now do what ever initialization is needed
  if (!edm::Service<SiStripDetInfoFileReader>().isAvailable()) {
    edm::LogError("TkLayerMap") << "\n------------------------------------------"
                                   "\nUnAvailable Service SiStripDetInfoFileReader: please insert in "
                                   "the configuration file an instance like"
                                   "\n\tprocess.SiStripDetInfoFileReader = "
                                   "cms.Service(\"SiStripDetInfoFileReader\")"
                                   "\n------------------------------------------";
  }

  fr = edm::Service<SiStripDetInfoFileReader>().operator->();
  file = new TFile("correlTest.root", "RECREATE");

  file->cd();
  tkmap = new TrackerMap();
}

SiStripCorrelateNoise::~SiStripCorrelateNoise() {}

//

void SiStripCorrelateNoise::beginRun(const edm::Run &run, const edm::EventSetup &es) {
  if (getNoiseCache(es) == cacheID_noise)
    return;

  edm::LogInfo("") << "[SiStripCorrelateNoise::beginRun]  cacheID_noise " << cacheID_noise << std::endl;

  es.get<SiStripNoisesRcd>().get(noiseHandle_);
  SiStripNoises *aNoise = new SiStripNoises(*noiseHandle_.product());

  // Check if gain is the same from one noise iov to the other, otherwise cache
  // the new gain (and the old one) to rescale

  checkGainCache(es);

  if (cacheID_noise != 0xFFFFFFFF) {
    char dir[128];
    theRun = run.run();
    sprintf(dir, "Run_%d", theRun);
    file->cd("");
    file->mkdir(dir);
    file->cd(dir);
    DoAnalysis(es, *noiseHandle_.product(), *refNoise);
    DoPlots();
  }

  cacheID_noise = getNoiseCache(es);
  if (refNoise != nullptr)
    delete refNoise;
  refNoise = aNoise;
}

void SiStripCorrelateNoise::checkGainCache(const edm::EventSetup &es) {
  equalGain = true;
  if (getGainCache(es) != cacheID_gain) {
    es.get<SiStripApvGainRcd>().get(gainHandle_);
    if (oldGain != nullptr)
      delete oldGain;

    oldGain = newGain;
    newGain = new SiStripApvGain(*gainHandle_.product());

    if (cacheID_gain != 0xFFFFFFFF)
      equalGain = false;
    cacheID_gain = getGainCache(es);
    edm::LogInfo("") << "[SiStripCorrelateNoise::checkGainCache]  cacheID_gain " << cacheID_gain << std::endl;
  }
}

void SiStripCorrelateNoise::DoPlots() {
  TCanvas *C = new TCanvas();
  C->Divide(2, 2);

  char outName[128];
  sprintf(outName, "Run_%d.png", theRun);
  for (size_t i = 0; i < vTH1.size(); i++)
    if (vTH1[i] != nullptr) {
      if (i % 100 == 0) {
        C->cd(i / 100);
        vTH1[i]->SetLineColor(i / 100);
        vTH1[i]->Draw();
        C->cd(i / 100)->SetLogy();
      }
      vTH1[i]->Write();
    }

  C->Print(outName);
  delete C;

  vTH1.clear();
  file->cd("");

  char dir[128];
  sprintf(dir, "Run_%d_TkMap.png", theRun);
  tkmap->save(false, 0, 5, dir);
  delete tkmap;
  tkmap = new TrackerMap();
}

void SiStripCorrelateNoise::DoAnalysis(const edm::EventSetup &es,
                                       const SiStripNoises &_Noise,
                                       SiStripNoises &refNoise) {
  SiStripNoises Noise = _Noise;
  typedef std::vector<SiStripNoises::ratioData> collection;
  collection divNoise = Noise / refNoise;

  edm::LogInfo("") << "[Doanalysis]";

  // Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHandle;
  es.get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology *const tTopo = tTopoHandle.product();

  std::vector<TH1F *> histos;

  collection::const_iterator iter = divNoise.begin();
  collection::const_iterator iterE = divNoise.end();

  float value;
  float gainRatio = 1.;
  // Divide result by d
  for (; iter != iterE; ++iter) {
    getHistos(iter->detid, tTopo, histos);

    size_t strip = 0, stripE = iter->values.size();
    size_t apvNb = 7;

    for (; strip < stripE; ++strip) {
      if (!equalGain && strip / 128 != apvNb) {
        apvNb = strip / 128;
        if (apvNb < 6)
          gainRatio = getGainRatio(iter->detid, apvNb);
        else
          edm::LogInfo("") << "[Doanalysis] detid " << iter->detid << " strip " << strip << " apvNb " << apvNb;
      }
      // edm::LogInfo("") << "[Doanalysis] detid " << iter->detid << " strip "
      // << strip << " value " << iter->values[strip];
      value = iter->values[strip] * gainRatio;
      tkmap->fill(iter->detid, value);
      for (size_t i = 0; i < histos.size(); ++i)
        histos[i]->Fill(value);
    }
  }
}

float SiStripCorrelateNoise::getGainRatio(const uint32_t &detid, const uint16_t &apv) {
  SiStripApvGain::Range oldRange = oldGain->getRange(detid);
  SiStripApvGain::Range newRange = newGain->getRange(detid);

  if (oldRange.first == oldRange.second || newRange.first == newRange.second)
    return 1.;

  return oldGain->getApvGain(apv, oldRange) / newGain->getApvGain(apv, newRange);
}

float SiStripCorrelateNoise::getMeanNoise(const SiStripNoises::Range &noiseRange,
                                          const uint32_t &firstStrip,
                                          const uint32_t &range) {
  float mean = 0;
  for (size_t istrip = firstStrip; istrip < firstStrip + range; istrip++) {
    mean += noiseHandle_->getNoise(istrip, noiseRange);
  }
  return mean / (1. * range);
}

void SiStripCorrelateNoise::getHistos(const uint32_t &detid,
                                      const TrackerTopology *tTopo,
                                      std::vector<TH1F *> &histos) {
  histos.clear();

  int subdet = -999;
  int component = -999;
  SiStripDetId a(detid);
  if (a.subdetId() == 3) {
    subdet = 0;
    component = tTopo->tibLayer(detid);
  } else if (a.subdetId() == 4) {
    subdet = 1;
    component = tTopo->tidSide(detid) == 2 ? tTopo->tidWheel(detid) : tTopo->tidWheel(detid) + 3;
  } else if (a.subdetId() == 5) {
    subdet = 2;
    component = tTopo->tobLayer(detid);
  } else if (a.subdetId() == 6) {
    subdet = 3;
    component = tTopo->tecSide(detid) == 2 ? tTopo->tecWheel(detid) : tTopo->tecWheel(detid) + 9;
  }

  int index = 100 + subdet * 100 + component;

  histos.push_back(getHisto(100 + 100 * subdet));
  histos.push_back(getHisto(index));
}

TH1F *SiStripCorrelateNoise::getHisto(const long unsigned int &index) {
  if (vTH1.size() < index + 1)
    vTH1.resize(index + 1, nullptr);

  if (vTH1[index] == nullptr) {
    char name[128];
    std::string SubD;
    if (index < 200)
      SubD = "TIB";
    else if (index < 300)
      SubD = "TID";
    else if (index < 400)
      SubD = "TOB";
    else
      SubD = "TEC";
    sprintf(name, "%d_%lu__%s", theRun, index, SubD.c_str());
    edm::LogInfo("") << "[getHisto] creating index " << index << std::endl;
    vTH1[index] = new TH1F(name, name, 200, -0.5, 10.5);
  }

  return vTH1[index];
}

void SiStripCorrelateNoise::endJob() {
  file->Write();
  file->Close();
}