CastorDigiMonitor.cc

Go to the documentation of this file.

//****************************************************//
//********** CastorDigiMonitor: ******************//
//********** Author: Dmytro Volyanskyy   *************//
//********** Date  : 29.08.2008 (first version) ******//
//****************************************************//
//---- critical revision 26.06.2014 (Vladimir Popov)
//     add rms check, DB   15.04.2015 (Vladimir Popov)
//==================================================================//

#include "CondFormats/CastorObjects/interface/CastorChannelQuality.h"
#include "DQM/CastorMonitor/interface/CastorDigiMonitor.h"
#include "DQM/CastorMonitor/interface/CastorLEDMonitor.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/LuminosityBlock.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <string>

using namespace std;
using namespace edm;

namespace {
  vector<std::string> HltPaths_;
  int StatusBadChannel = 1;
  int ChannelStatus[14][16]{};
  int N_GoodChannels = 224;
  int EtowerLastModule = 5;
  int TrigIndexMax = 0;
}  // namespace

CastorDigiMonitor::CastorDigiMonitor(const edm::ParameterSet &ps) {
  fVerbosity = ps.getUntrackedParameter<int>("debug", 0);
  subsystemname_ = ps.getUntrackedParameter<std::string>("subSystemFolder", "Castor");
  EtowerLastModule = ps.getUntrackedParameter<int>("towerLastModule", 6);
  RatioThresh1 = ps.getUntrackedParameter<double>("ratioThreshold", 0.9);
  Qrms_DEAD = ps.getUntrackedParameter<double>("QrmsDead", 0.01);  // fC
  HltPaths_ = ps.getParameter<vector<string>>("HltPaths");

  Qrms_DEAD = Qrms_DEAD * Qrms_DEAD;
  TS_MAX = ps.getUntrackedParameter<double>("qieTSmax", 6);
  StatusBadChannel = CastorChannelStatus::StatusBit::BAD;
  if (fVerbosity > 0)
    LogPrint("CastorDigi") << "enum CastorChannelStatus::StatusBit::BAD=" << StatusBadChannel
                           << "EtowerLastModule = " << EtowerLastModule << endl;
}

CastorDigiMonitor::~CastorDigiMonitor() {}

void CastorDigiMonitor::bookHistograms(DQMStore::IBooker &ibooker,
                                       const edm::Run &iRun,
                                       const edm::EventSetup &iSetup) {
  char s[60];
  string st;
  if (fVerbosity > 0)
    LogPrint("CastorMonitorModule") << "Digi bookHist(start)";

  getDbData(iSetup);

  char sTileIndex[50];
  sprintf(sTileIndex, "Cell(=moduleZ*16+sector#phi)");

  ievt_ = 0;

  ibooker.setCurrentFolder(subsystemname_);
  hBX = ibooker.bookProfile(
      "average E(digi) in BX", "Castor average E (digi);Event.BX;fC", 3601, -0.5, 3600.5, 0., 1.e10, "");
  hBX->getTProfile()->SetOption("hist");

  string trname = HltPaths_[0];
  hpBXtrig = ibooker.bookProfile("average E(digi) in BXtrig",
                                 "Castor average E (digi) trigger:'" + trname + "';Event.BX;fC",
                                 3601,
                                 -0.5,
                                 3600.5,
                                 0.,
                                 1.e10,
                                 "");
  hpBXtrig->getTProfile()->SetOption("hist");

  hpTrigRes = ibooker.bookProfile(
      "E(digi)vsTriggerIndex", "Castor average E(digi) by triggerIndex;triggerIndex;fC", 512, 0., 512, 0., 1.e10, "");
  hpTrigRes->getTProfile()->SetOption("hist");

  ibooker.setCurrentFolder(subsystemname_ + "/CastorDigiMonitor");

  std::string s2 = "CASTOR QIE_capID+er+dv";
  h2digierr = ibooker.bookProfile2D(s2, s2, 14, 0., 14., 16, 0., 16., 100, 0, 1.e10, "");
  h2digierr->getTProfile2D()->GetXaxis()->SetTitle("Module Z");
  h2digierr->getTProfile2D()->GetYaxis()->SetTitle("Sector #phi");
  h2digierr->getTProfile2D()->SetMaximum(1.);
  h2digierr->getTProfile2D()->SetMinimum(QIEerrThreshold);
  h2digierr->getTProfile2D()->SetOption("colz");

  sprintf(s, "CASTORreportSummaryMap");
  h2repsum = ibooker.bookProfile2D(s, s, 14, 0., 14., 16, 0., 16., 100, 0, 1.e10, "");
  h2repsum->getTProfile2D()->GetXaxis()->SetTitle("Module Z");
  h2repsum->getTProfile2D()->GetYaxis()->SetTitle("Sector #phi");
  h2repsum->getTProfile2D()->SetMaximum(1.);
  h2repsum->getTProfile2D()->SetMinimum(QIEerrThreshold);
  h2repsum->getTProfile2D()->SetOption("colz");

  sprintf(s, "CASTOR BadChannelsMap");
  h2status = ibooker.book2D(s, s, 14, 0., 14., 16, 0., 16.);
  h2status->setAxisTitle("Module Z");
  h2status->setAxisTitle("Sector #phi", /* axis */ 2);
  h2status->setOption("colz");

  sprintf(s, "CASTOR TSmax Significance Map");
  h2TSratio = ibooker.book2D(s, s, 14, 0., 14., 16, 0., 16.);
  h2TSratio->setAxisTitle("Module Z");
  h2TSratio->setAxisTitle("Sector #phi", /* axis */ 2);
  h2TSratio->setOption("colz");

  sprintf(s, "CASTOR TSmax Significance All chan");
  hTSratio = ibooker.book1D(s, s, 105, 0., 1.05);

  sprintf(s, "DigiSize");
  hdigisize = ibooker.book1DD(s, s, 20, 0., 20.);
  sprintf(s, "ModuleZ(fC)_allTS");
  hModule = ibooker.book1D(s, s, 14, 0., 14.);
  hModule->setAxisTitle("ModuleZ");
  hModule->setAxisTitle("QIE(fC)", /* axis */ 2);
  sprintf(s, "Sector #phi(fC)_allTS");
  hSector = ibooker.book1D(s, s, 16, 0., 16.);
  hSector->setAxisTitle("Sector #phi");
  hSector->setAxisTitle("QIE(fC)", /* axis */ 2);

  st = "Castor cells avr digi(fC) per event Map TS vs Channel";
  h2QmeantsvsCh =
      ibooker.bookProfile2D(st, st + ";" + string(sTileIndex) + ";TS", 224, 0., 224., 10, 0., 10., 0., 1.e10, "");
  h2QmeantsvsCh->getTProfile2D()->SetOption("colz");

  st = "Castor cells avr digiRMS(fC) per event Map TS vs Channel";
  h2QrmsTSvsCh = ibooker.book2D(st, st + ";" + string(sTileIndex) + ";TS", 224, 0., 224., 10, 0., 10.);
  h2QrmsTSvsCh->setOption("colz");

  sprintf(s, "CASTOR data quality");
  h2qualityMap = ibooker.book2D(s, s, 14, 0, 14, 16, 0, 16);
  h2qualityMap->setAxisTitle("module Z");
  h2qualityMap->setAxisTitle("Sector #phi", /* axis */ 2);
  h2qualityMap->setOption("colz");

  hReport = ibooker.bookFloat("CASTOR reportSummary");

  sprintf(s, "QmeanfC_map(allTS)");
  h2QmeanMap = ibooker.book2D(s, s, 14, 0., 14., 16, 0., 16.);
  h2QmeanMap->setAxisTitle("Module Z");
  h2QmeanMap->setAxisTitle("Sector #phi", /* axis */ 2);
  h2QmeanMap->setOption("textcolz");

  const int NEtow = 20;
  float EhadTow[NEtow + 1];
  float EMTow[NEtow + 1];
  float ETower[NEtow + 2];
  double E0tow = 500. / 1024.;
  EMTow[0] = 0.;
  EMTow[1] = E0tow;
  EhadTow[0] = 0.;
  EhadTow[1] = E0tow;
  ETower[0] = 0.;
  ETower[1] = E0tow;
  double lnBtow = log(1.8);  // 2.
  for (int j = 1; j < NEtow; j++)
    EMTow[j + 1] = E0tow * exp(j * lnBtow);
  for (int j = 1; j < NEtow; j++)
    EhadTow[j + 1] = E0tow * exp(j * lnBtow);
  for (int j = 1; j <= NEtow; j++)
    ETower[j + 1] = E0tow * exp(j * lnBtow);

  sprintf(s, "CASTOR_Tower_EMvsEhad(fC)");
  h2towEMvsHAD = ibooker.book2D(s, s, NEtow, EhadTow, NEtow, EMTow);
  h2towEMvsHAD->setAxisTitle("Ehad [fC]");
  h2towEMvsHAD->setAxisTitle("EM [fC]", /* axis */ 2);
  h2towEMvsHAD->setOption("colz");

  sprintf(s, "CASTOR_TowerTotalEnergy(fC)");
  htowE = ibooker.book1D(s, s, NEtow + 1, ETower);
  htowE->setAxisTitle("fC");

  for (int ts = 0; ts <= 1; ts++) {
    sprintf(s, "QIErms_TS=%d", ts);
    hQIErms[ts] = ibooker.book1D(s, s, 1000, 0., 100.);
    hQIErms[ts]->setAxisTitle("QIErms(fC)");
  }

  for (int ind = 0; ind < 224; ind++)
    for (int ts = 0; ts < 10; ts++)
      QrmsTS[ind][ts] = QmeanTS[ind][ts] = 0.;

  return;
}

void CastorDigiMonitor::processEvent(edm::Event const &event,
                                     const CastorDigiCollection &castorDigis,
                                     const edm::TriggerResults &TrigResults,
                                     const CastorDbService &cond) {
  if (fVerbosity > 1)
    LogPrint("CastorDigiMonitor") << "processEvent(begin)";

  if (castorDigis.empty()) {
    for (int mod = 0; mod < 14; mod++)
      for (int sec = 0; sec < 16; sec++)
        h2repsum->Fill(mod, sec, 0.);
    hBX->Fill(event.bunchCrossing(), 0.);
    fillTrigRes(event, TrigResults, 0.);
    return;
  }

  float Ecell[14][16]{};
  for (CastorDigiCollection::const_iterator j = castorDigis.begin(); j != castorDigis.end(); j++) {
    const CastorDataFrame digi = (const CastorDataFrame)(*j);

    int module = digi.id().module() - 1;
    int sector = digi.id().sector() - 1;
    if (ChannelStatus[module][sector] == StatusBadChannel)
      continue;

    int capid1 = digi.sample(0).capid();
    hdigisize->Fill(digi.size());
    double sum = 0.;
    int err = 0, err2 = 0;
    for (int i = 0; i < digi.size(); i++) {
      int capid = digi.sample(i).capid();
      int dv = digi.sample(i).dv();
      int er = digi.sample(i).er();
      int rawd = digi.sample(i).adc();
      rawd = rawd & 0x7F;
      err |= (capid != capid1) | er << 1 | (!dv) << 2;  // =0
      err2 += (capid != capid1) | er | (!dv);           // =0
      //     if(err !=0) continue;
      int ind = ModSecToIndex(module, sector);
      h2QmeantsvsCh->Fill(ind, i, LedMonAdc2fc[rawd]);
      float q = LedMonAdc2fc[rawd];
      Ecell[module][sector] = q;
      sum += q;  //     sum += LedMonAdc2fc[rawd];
      QrmsTS[ind][i] += (q * q);
      QmeanTS[ind][i] += q;
      if (err != 0 && fVerbosity > 0)
        LogPrint("CastorDigiMonitor") << "event/idigi=" << ievt_ << "/" << i << " cap=cap1_dv_er_err: " << capid << "="
                                      << capid1 << " " << dv << " " << er << " " << err;
      if (capid1 < 3)
        capid1 = capid + 1;
      else
        capid1 = 0;
    }
    h2digierr->Fill(module, sector, err);
    h2repsum->Fill(module, sector, 1. - err2 / digi.size());
  }  // end for(CastorDigiCollection::const_iterator ...

  ievt_++;

  double Etotal = 0.;
  for (int sec = 0; sec < 16; sec++)
    for (int mod = 0; mod < 14; mod++)
      Etotal = Ecell[mod][sec];
  hBX->Fill(event.bunchCrossing(), Etotal);
  fillTrigRes(event, TrigResults, Etotal);

  for (int sec = 0; sec < 16; sec++) {
    float em = Ecell[0][sec] + Ecell[1][sec];
    double ehad = 0.;
    for (int mod = 2; mod < EtowerLastModule; mod++)
      ehad += Ecell[mod][sec];
    h2towEMvsHAD->Fill(em, ehad);
    htowE->Fill(em + ehad);
  }

  const float repChanBAD = 0.9;
  const float repChanWarning = 0.95;
  if (ievt_ % 100 != 0)
    return;

  float ModuleSum[14], SectorSum[16];
  for (int m = 0; m < 14; m++)
    ModuleSum[m] = 0.;
  for (int s = 0; s < 16; s++)
    SectorSum[s] = 0.;
  for (int mod = 0; mod < 14; mod++)
    for (int sec = 0; sec < 16; sec++) {
      for (int ts = 0; ts <= 1; ts++) {
        int ind = ModSecToIndex(mod, sec);
        double Qmean = QmeanTS[ind][ts] / ievt_;
        double Qrms = sqrt(QrmsTS[ind][ts] / ievt_ - Qmean * Qmean);
        hQIErms[ts]->Fill(Qrms);
      }

      double sum = 0.;
      for (int ts = 1; ts <= TS_MAX; ts++) {
        int ind = ModSecToIndex(mod, sec) + 1;
        double a =  //(1) h2QtsvsCh->getTH2D()->GetBinContent(ind,ts);
            h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind, ts);
        sum += a;
        double Qmean = QmeanTS[ind - 1][ts - 1] / ievt_;
        double Qrms = QrmsTS[ind - 1][ts - 1] / ievt_ - Qmean * Qmean;
        h2QrmsTSvsCh->setBinContent(ind, ts, sqrt(Qrms));
      }
      ModuleSum[mod] += sum;
      SectorSum[sec] += sum;
      float isum = float(int(sum * 10. + 0.5)) / 10.;
      if (ChannelStatus[mod][sec] != StatusBadChannel)
        h2QmeanMap->setBinContent(mod + 1, sec + 1, isum);
    }  // end for(int mod=0; mod<14; mod++) for(int sec=0;...

  for (int mod = 0; mod < 14; mod++)
    hModule->setBinContent(mod + 1, ModuleSum[mod]);
  for (int sec = 0; sec < 16; sec++)
    hSector->setBinContent(sec + 1, SectorSum[sec]);

  int nGoodCh = 0;
  hTSratio->Reset();
  for (int mod = 0; mod < 14; mod++)
    for (int sec = 0; sec < 16; sec++) {
      if (ChannelStatus[mod][sec] == StatusBadChannel)
        continue;
      int ind = ModSecToIndex(mod, sec);
      double Qmean = QmeanTS[ind][TSped] / ievt_;
      double Qrms = QrmsTS[ind][TSped] / ievt_ - Qmean * Qmean;
      float ChanStatus = 0.;
      if (Qrms < Qrms_DEAD)
        ChanStatus = 1.;
      h2status->setBinContent(mod + 1, sec + 1, ChanStatus);

      float am = 0.;
      for (int ts = 0; ts < TS_MAX - 1; ts++) {
        float a = h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 1) +
                  h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 2);
        if (am < a)
          am = a;
      }

      double sum = 0.;
      for (int ts = 0; ts < TS_MAX; ts++)
        sum += h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 1);

      float r = 0.;  // worth case - no peak
      if (am > 0.)
        r = 1. - (sum - am) / (TS_MAX - 2) / am * 2.;
      // if(r<0.|| r>1.) cout<<"ievt="<<ievt<<" r="<<r<<" amax= "<<am<<"
      // sum="<<sum<<endl;
      h2TSratio->setBinContent(mod + 1, sec + 1, r);
      hTSratio->Fill(r);

      float statusTS = 1.0;
      if (r > RatioThresh1)
        statusTS = repChanWarning;
      else if (r > 0.99)
        statusTS = repChanBAD;
      float gChanStatus = statusTS;
      if (ChanStatus > 0.)
        gChanStatus = repChanBAD;  // RMS
      h2qualityMap->setBinContent(mod + 1, sec + 1, gChanStatus);
      if (gChanStatus > repChanBAD)
        ++nGoodCh;
    }
  hReport->Fill(float(nGoodCh) / N_GoodChannels);
  return;
}

void CastorDigiMonitor::endRun() {
  if (fVerbosity > 0)
    LogPrint("CastorDigiMonitor") << "DigiMonitor::endRun: trigger max index = " << TrigIndexMax
                                  << "  TriggerIndexies(N):" << endl;
  for (int i = 1; i < hpTrigRes->getTProfile()->GetNbinsX(); i++)
    if (hpTrigRes->getTProfile()->GetBinContent(i) > 0)
      LogPrint("CastorDigiMonitor") << i - 1 << "(" << hpTrigRes->getTProfile()->GetBinContent(i) << ") ";
}

void CastorDigiMonitor::fillTrigRes(edm::Event const &event, const edm::TriggerResults &TrigResults, double Etotal) {
  int nTriggers = TrigResults.size();
  const edm::TriggerNames &trigName = event.triggerNames(TrigResults);
  bool event_triggered = false;
  if (nTriggers > 0)
    for (int iTrig = 0; iTrig < nTriggers; ++iTrig) {
      if (TrigResults.accept(iTrig)) {
        int index = trigName.triggerIndex(trigName.triggerName(iTrig));
        if (TrigIndexMax < index)
          TrigIndexMax = index;
        if (fVerbosity > 0)
          LogPrint("CastorDigi") << "trigger[" << iTrig << "] name:" << trigName.triggerName(iTrig)
                                 << " index= " << index << endl;
        hpTrigRes->Fill(index, Etotal);
        for (int n = 0; n < int(HltPaths_.size()); n++) {
          if (trigName.triggerName(iTrig).find(HltPaths_[n]) != std::string::npos)
            event_triggered = true;
        }
      }  // end if(TrigResults.accept(iTrig)
    }

  if (event_triggered)
    hpBXtrig->Fill(event.bunchCrossing(), Etotal);
  return;
}

void CastorDigiMonitor::getDbData(const edm::EventSetup &iSetup) {
  edm::ESHandle<CastorChannelQuality> dbChQuality;
  iSetup.get<CastorChannelQualityRcd>().get(dbChQuality);
  if (fVerbosity > 0) {
    LogPrint("CastorDigiMonitor") << " CastorChQuality in CondDB=" << dbChQuality.isValid() << endl;
  }

  int chInd = 0;
  for (int mod = 0; mod < 14; mod++)
    for (int sec = 0; sec < 16; sec++)
      ChannelStatus[mod][sec] = 0;
  std::vector<DetId> channels = dbChQuality->getAllChannels();
  N_GoodChannels = 224 - channels.size();
  if (fVerbosity > 0)
    LogPrint("CastorDigiMonitor") << "CastorDigiMonitor::getDBData: QualityRcdSize=" << channels.size();
  for (std::vector<DetId>::iterator ch = channels.begin(); ch != channels.end(); ch++) {
    const CastorChannelStatus *quality = dbChQuality->getValues(*ch);
    int value = quality->getValue();
    int rawId = quality->rawId();
    chInd++;
    int mod = HcalCastorDetId(*ch).module() - 1;
    int sec = HcalCastorDetId(*ch).sector() - 1;
    if (mod > 0 && mod < 16 && sec > 0 && sec < 16)
      ChannelStatus[mod][sec] = value;
    if (fVerbosity > 0)
      LogPrint("CastorDigiMonitor") << chInd << " module=" << mod << " sec=" << sec << " rawId=" << rawId
                                    << " value=" << value << endl;
  }  // end for(std::vector<DetId>::it...
  return;
}

int CastorDigiMonitor::ModSecToIndex(int module, int sector) {
  int ind = sector + module * 16;
  if (ind > 223)
    ind = 223;
  return (ind);
}