CastorMonitor/src/CastorDigiMonitor.cc

0001 //****************************************************//
0002 //********** CastorDigiMonitor: ******************//
0003 //********** Author: Dmytro Volyanskyy   *************//
0004 //********** Date  : 29.08.2008 (first version) ******//
0005 ////---- digi values in Castor r/o channels
0006 //// last revision: 31.05.2011 (Panos Katsas) to remove selecting N events for
0007 /// filling the histograms
0008 //****************************************************//
0009 //---- critical revision 26.06.2014 (Vladimir Popov)
0010 //     add rms check, DB   15.04.2015 (Vladimir Popov)
0011 //==================================================================//
0012
0013 #include "DQM/CastorMonitor/interface/CastorDigiMonitor.h"
0014 #include "DQM/CastorMonitor/interface/CastorLEDMonitor.h"
0015 #include "DQMServices/Core/interface/DQMStore.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/LuminosityBlock.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/Utilities/interface/Transition.h"
0021
0022 #include <string>
0023
0024 using namespace std;
0025 using namespace edm;
0026
0027 namespace {
0028   vector<std::string> HltPaths_;
0029   int StatusBadChannel = 1;
0030   int ChannelStatus[14][16]{};
0031   int N_GoodChannels = 224;
0032   int EtowerLastModule = 5;
0033   int TrigIndexMax = 0;
0034 }  // namespace
0035
0036 CastorDigiMonitor::CastorDigiMonitor(const edm::ParameterSet &ps, edm::ConsumesCollector &&iC)
0037     : castorChannelQualityToken_{
0038           iC.esConsumes<CastorChannelQuality, CastorChannelQualityRcd, edm::Transition::BeginRun>()} {
0039   fVerbosity = ps.getUntrackedParameter<int>("debug", 0);
0040   subsystemname_ = ps.getUntrackedParameter<std::string>("subSystemFolder", "Castor");
0041   EtowerLastModule = ps.getUntrackedParameter<int>("towerLastModule", 6);
0042   RatioThresh1 = ps.getUntrackedParameter<double>("ratioThreshold", 0.9);
0043   Qrms_DEAD = ps.getUntrackedParameter<double>("QrmsDead", 0.01);  // fC
0044   HltPaths_ = ps.getParameter<vector<string>>("HltPaths");
0045
0046   Qrms_DEAD = Qrms_DEAD * Qrms_DEAD;
0047   TS_MAX = ps.getUntrackedParameter<double>("qieTSmax", 6);
0048   StatusBadChannel = CastorChannelStatus::StatusBit::BAD;
0049   if (fVerbosity > 0)
0050     LogPrint("CastorDigi") << "enum CastorChannelStatus::StatusBit::BAD=" << StatusBadChannel
0051                            << "EtowerLastModule = " << EtowerLastModule << endl;
0052 }
0053
0054 CastorDigiMonitor::~CastorDigiMonitor() {}
0055
0056 void CastorDigiMonitor::bookHistograms(DQMStore::IBooker &ibooker,
0057                                        const edm::Run &iRun,
0058                                        const edm::EventSetup &iSetup) {
0059   char s[60];
0060   string st;
0061   if (fVerbosity > 0)
0062     LogPrint("CastorMonitorModule") << "Digi bookHist(start)";
0063
0064   getDbData(iSetup);
0065
0066   char sTileIndex[50];
0067   sprintf(sTileIndex, "Cell(=moduleZ*16+sector#phi)");
0068
0069   ievt_ = 0;
0070
0071   ibooker.setCurrentFolder(subsystemname_);
0072   hBX = ibooker.bookProfile(
0073       "average E(digi) in BX", "Castor average E (digi);Event.BX;fC", 3601, -0.5, 3600.5, 0., 1.e10, "");
0074   hBX->getTProfile()->SetOption("hist");
0075
0076   string trname = HltPaths_[0];
0077   hpBXtrig = ibooker.bookProfile("average E(digi) in BXtrig",
0078                                  "Castor average E (digi) trigger:'" + trname + "';Event.BX;fC",
0079                                  3601,
0080                                  -0.5,
0081                                  3600.5,
0082                                  0.,
0083                                  1.e10,
0084                                  "");
0085   hpBXtrig->getTProfile()->SetOption("hist");
0086
0087   hpTrigRes = ibooker.bookProfile(
0088       "E(digi)vsTriggerIndex", "Castor average E(digi) by triggerIndex;triggerIndex;fC", 512, 0., 512, 0., 1.e10, "");
0089   hpTrigRes->getTProfile()->SetOption("hist");
0090
0091   ibooker.setCurrentFolder(subsystemname_ + "/CastorDigiMonitor");
0092
0093   std::string s2 = "CASTOR QIE_capID+er+dv";
0094   h2digierr = ibooker.bookProfile2D(s2, s2, 14, 0., 14., 16, 0., 16., 100, 0, 1.e10, "");
0095   h2digierr->getTProfile2D()->GetXaxis()->SetTitle("Module Z");
0096   h2digierr->getTProfile2D()->GetYaxis()->SetTitle("Sector #phi");
0097   h2digierr->getTProfile2D()->SetMaximum(1.);
0098   h2digierr->getTProfile2D()->SetMinimum(QIEerrThreshold);
0099   h2digierr->getTProfile2D()->SetOption("colz");
0100
0101   sprintf(s, "CASTORreportSummaryMap");
0102   h2repsum = ibooker.bookProfile2D(s, s, 14, 0., 14., 16, 0., 16., 100, 0, 1.e10, "");
0103   h2repsum->getTProfile2D()->GetXaxis()->SetTitle("Module Z");
0104   h2repsum->getTProfile2D()->GetYaxis()->SetTitle("Sector #phi");
0105   h2repsum->getTProfile2D()->SetMaximum(1.);
0106   h2repsum->getTProfile2D()->SetMinimum(QIEerrThreshold);
0107   h2repsum->getTProfile2D()->SetOption("colz");
0108
0109   sprintf(s, "CASTOR BadChannelsMap");
0110   h2status = ibooker.book2D(s, s, 14, 0., 14., 16, 0., 16.);
0111   h2status->setAxisTitle("Module Z");
0112   h2status->setAxisTitle("Sector #phi", /* axis */ 2);
0113   h2status->setOption("colz");
0114
0115   sprintf(s, "CASTOR TSmax Significance Map");
0116   h2TSratio = ibooker.book2D(s, s, 14, 0., 14., 16, 0., 16.);
0117   h2TSratio->setAxisTitle("Module Z");
0118   h2TSratio->setAxisTitle("Sector #phi", /* axis */ 2);
0119   h2TSratio->setOption("colz");
0120
0121   sprintf(s, "CASTOR TSmax Significance All chan");
0122   hTSratio = ibooker.book1D(s, s, 105, 0., 1.05);
0123
0124   sprintf(s, "DigiSize");
0125   hdigisize = ibooker.book1DD(s, s, 20, 0., 20.);
0126   sprintf(s, "ModuleZ(fC)_allTS");
0127   hModule = ibooker.book1D(s, s, 14, 0., 14.);
0128   hModule->setAxisTitle("ModuleZ");
0129   hModule->setAxisTitle("QIE(fC)", /* axis */ 2);
0130   sprintf(s, "Sector #phi(fC)_allTS");
0131   hSector = ibooker.book1D(s, s, 16, 0., 16.);
0132   hSector->setAxisTitle("Sector #phi");
0133   hSector->setAxisTitle("QIE(fC)", /* axis */ 2);
0134
0135   st = "Castor cells avr digi(fC) per event Map TS vs Channel";
0136   h2QmeantsvsCh =
0137       ibooker.bookProfile2D(st, st + ";" + string(sTileIndex) + ";TS", 224, 0., 224., 10, 0., 10., 0., 1.e10, "");
0138   h2QmeantsvsCh->getTProfile2D()->SetOption("colz");
0139
0140   st = "Castor cells avr digiRMS(fC) per event Map TS vs Channel";
0141   h2QrmsTSvsCh = ibooker.book2D(st, st + ";" + string(sTileIndex) + ";TS", 224, 0., 224., 10, 0., 10.);
0142   h2QrmsTSvsCh->setOption("colz");
0143
0144   sprintf(s, "CASTOR data quality");
0145   h2qualityMap = ibooker.book2D(s, s, 14, 0, 14, 16, 0, 16);
0146   h2qualityMap->setAxisTitle("module Z");
0147   h2qualityMap->setAxisTitle("Sector #phi", /* axis */ 2);
0148   h2qualityMap->setOption("colz");
0149
0150   hReport = ibooker.bookFloat("CASTOR reportSummary");
0151
0152   sprintf(s, "QmeanfC_map(allTS)");
0153   h2QmeanMap = ibooker.book2D(s, s, 14, 0., 14., 16, 0., 16.);
0154   h2QmeanMap->setAxisTitle("Module Z");
0155   h2QmeanMap->setAxisTitle("Sector #phi", /* axis */ 2);
0156   h2QmeanMap->setOption("textcolz");
0157
0158   const int NEtow = 20;
0159   float EhadTow[NEtow + 1];
0160   float EMTow[NEtow + 1];
0161   float ETower[NEtow + 2];
0162   double E0tow = 500. / 1024.;
0163   EMTow[0] = 0.;
0164   EMTow[1] = E0tow;
0165   EhadTow[0] = 0.;
0166   EhadTow[1] = E0tow;
0167   ETower[0] = 0.;
0168   ETower[1] = E0tow;
0169   double lnBtow = log(1.8);  // 2.
0170   for (int j = 1; j < NEtow; j++)
0171     EMTow[j + 1] = E0tow * exp(j * lnBtow);
0172   for (int j = 1; j < NEtow; j++)
0173     EhadTow[j + 1] = E0tow * exp(j * lnBtow);
0174   for (int j = 1; j <= NEtow; j++)
0175     ETower[j + 1] = E0tow * exp(j * lnBtow);
0176
0177   sprintf(s, "CASTOR_Tower_EMvsEhad(fC)");
0178   h2towEMvsHAD = ibooker.book2D(s, s, NEtow, EhadTow, NEtow, EMTow);
0179   h2towEMvsHAD->setAxisTitle("Ehad [fC]");
0180   h2towEMvsHAD->setAxisTitle("EM [fC]", /* axis */ 2);
0181   h2towEMvsHAD->setOption("colz");
0182
0183   sprintf(s, "CASTOR_TowerTotalEnergy(fC)");
0184   htowE = ibooker.book1D(s, s, NEtow + 1, ETower);
0185   htowE->setAxisTitle("fC");
0186
0187   for (int ts = 0; ts <= 1; ts++) {
0188     sprintf(s, "QIErms_TS=%d", ts);
0189     hQIErms[ts] = ibooker.book1D(s, s, 1000, 0., 100.);
0190     hQIErms[ts]->setAxisTitle("QIErms(fC)");
0191   }
0192
0193   for (int ind = 0; ind < 224; ind++)
0194     for (int ts = 0; ts < 10; ts++)
0195       QrmsTS[ind][ts] = QmeanTS[ind][ts] = 0.;
0196
0197   return;
0198 }
0199
0200 void CastorDigiMonitor::processEvent(edm::Event const &event,
0201                                      const CastorDigiCollection &castorDigis,
0202                                      const edm::TriggerResults &TrigResults,
0203                                      const CastorDbService &cond) {
0204   if (fVerbosity > 1)
0205     LogPrint("CastorDigiMonitor") << "processEvent(begin)";
0206
0207   if (castorDigis.empty()) {
0208     for (int mod = 0; mod < 14; mod++)
0209       for (int sec = 0; sec < 16; sec++)
0210         h2repsum->Fill(mod, sec, 0.);
0211     hBX->Fill(event.bunchCrossing(), 0.);
0212     fillTrigRes(event, TrigResults, 0.);
0213     return;
0214   }
0215
0216   float Ecell[14][16]{};
0217   for (CastorDigiCollection::const_iterator j = castorDigis.begin(); j != castorDigis.end(); j++) {
0218     const CastorDataFrame digi = (const CastorDataFrame)(*j);
0219
0220     int module = digi.id().module() - 1;
0221     int sector = digi.id().sector() - 1;
0222     if (ChannelStatus[module][sector] == StatusBadChannel)
0223       continue;
0224
0225     int capid1 = digi.sample(0).capid();
0226     hdigisize->Fill(digi.size());
0227     int err = 0, err2 = 0;
0228     for (int i = 0; i < digi.size(); i++) {
0229       int capid = digi.sample(i).capid();
0230       int dv = digi.sample(i).dv();
0231       int er = digi.sample(i).er();
0232       int rawd = digi.sample(i).adc();
0233       rawd = rawd & 0x7F;
0234       err |= (capid != capid1) | er << 1 | (!dv) << 2;  // =0
0235       err2 += (capid != capid1) | er | (!dv);           // =0
0236       //     if(err !=0) continue;
0237       int ind = ModSecToIndex(module, sector);
0238       h2QmeantsvsCh->Fill(ind, i, LedMonAdc2fc[rawd]);
0239       float q = LedMonAdc2fc[rawd];
0240       Ecell[module][sector] = q;
0241       QrmsTS[ind][i] += (q * q);
0242       QmeanTS[ind][i] += q;
0243       if (err != 0 && fVerbosity > 0)
0244         LogPrint("CastorDigiMonitor") << "event/idigi=" << ievt_ << "/" << i << " cap=cap1_dv_er_err: " << capid << "="
0245                                       << capid1 << " " << dv << " " << er << " " << err;
0246       if (capid1 < 3)
0247         capid1 = capid + 1;
0248       else
0249         capid1 = 0;
0250     }
0251     h2digierr->Fill(module, sector, err);
0252     h2repsum->Fill(module, sector, 1. - err2 / digi.size());
0253   }  // end for(CastorDigiCollection::const_iterator ...
0254
0255   ievt_++;
0256
0257   double Etotal = 0.;
0258   for (int sec = 0; sec < 16; sec++)
0259     for (int mod = 0; mod < 14; mod++)
0260       Etotal = Ecell[mod][sec];
0261   hBX->Fill(event.bunchCrossing(), Etotal);
0262   fillTrigRes(event, TrigResults, Etotal);
0263
0264   for (int sec = 0; sec < 16; sec++) {
0265     float em = Ecell[0][sec] + Ecell[1][sec];
0266     double ehad = 0.;
0267     for (int mod = 2; mod < EtowerLastModule; mod++)
0268       ehad += Ecell[mod][sec];
0269     h2towEMvsHAD->Fill(em, ehad);
0270     htowE->Fill(em + ehad);
0271   }
0272
0273   const float repChanBAD = 0.9;
0274   const float repChanWarning = 0.95;
0275   if (ievt_ % 100 != 0)
0276     return;
0277
0278   float ModuleSum[14], SectorSum[16];
0279   for (int m = 0; m < 14; m++)
0280     ModuleSum[m] = 0.;
0281   for (int s = 0; s < 16; s++)
0282     SectorSum[s] = 0.;
0283   for (int mod = 0; mod < 14; mod++)
0284     for (int sec = 0; sec < 16; sec++) {
0285       for (int ts = 0; ts <= 1; ts++) {
0286         int ind = ModSecToIndex(mod, sec);
0287         double Qmean = QmeanTS[ind][ts] / ievt_;
0288         double Qrms = sqrt(QrmsTS[ind][ts] / ievt_ - Qmean * Qmean);
0289         hQIErms[ts]->Fill(Qrms);
0290       }
0291
0292       double sum = 0.;
0293       for (int ts = 1; ts <= TS_MAX; ts++) {
0294         int ind = ModSecToIndex(mod, sec) + 1;
0295         double a =  //(1) h2QtsvsCh->getTH2D()->GetBinContent(ind,ts);
0296             h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind, ts);
0297         sum += a;
0298         double Qmean = QmeanTS[ind - 1][ts - 1] / ievt_;
0299         double Qrms = QrmsTS[ind - 1][ts - 1] / ievt_ - Qmean * Qmean;
0300         h2QrmsTSvsCh->setBinContent(ind, ts, sqrt(Qrms));
0301       }
0302       ModuleSum[mod] += sum;
0303       SectorSum[sec] += sum;
0304       float isum = float(int(sum * 10. + 0.5)) / 10.;
0305       if (ChannelStatus[mod][sec] != StatusBadChannel)
0306         h2QmeanMap->setBinContent(mod + 1, sec + 1, isum);
0307     }  // end for(int mod=0; mod<14; mod++) for(int sec=0;...
0308
0309   for (int mod = 0; mod < 14; mod++)
0310     hModule->setBinContent(mod + 1, ModuleSum[mod]);
0311   for (int sec = 0; sec < 16; sec++)
0312     hSector->setBinContent(sec + 1, SectorSum[sec]);
0313
0314   int nGoodCh = 0;
0315   hTSratio->Reset();
0316   for (int mod = 0; mod < 14; mod++)
0317     for (int sec = 0; sec < 16; sec++) {
0318       if (ChannelStatus[mod][sec] == StatusBadChannel)
0319         continue;
0320       int ind = ModSecToIndex(mod, sec);
0321       double Qmean = QmeanTS[ind][TSped] / ievt_;
0322       double Qrms = QrmsTS[ind][TSped] / ievt_ - Qmean * Qmean;
0323       float ChanStatus = 0.;
0324       if (Qrms < Qrms_DEAD)
0325         ChanStatus = 1.;
0326       h2status->setBinContent(mod + 1, sec + 1, ChanStatus);
0327
0328       float am = 0.;
0329       for (int ts = 0; ts < TS_MAX - 1; ts++) {
0330         float a = h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 1) +
0331                   h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 2);
0332         if (am < a)
0333           am = a;
0334       }
0335
0336       double sum = 0.;
0337       for (int ts = 0; ts < TS_MAX; ts++)
0338         sum += h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 1);
0339
0340       float r = 0.;  // worth case - no peak
0341       if (am > 0.)
0342         r = 1. - (sum - am) / (TS_MAX - 2) / am * 2.;
0343       // if(r<0.|| r>1.) cout<<"ievt="<<ievt<<" r="<<r<<" amax= "<<am<<"
0344       // sum="<<sum<<endl;
0345       h2TSratio->setBinContent(mod + 1, sec + 1, r);
0346       hTSratio->Fill(r);
0347
0348       float statusTS = 1.0;
0349       if (r > RatioThresh1)
0350         statusTS = repChanWarning;
0351       else if (r > 0.99)
0352         statusTS = repChanBAD;
0353       float gChanStatus = statusTS;
0354       if (ChanStatus > 0.)
0355         gChanStatus = repChanBAD;  // RMS
0356       h2qualityMap->setBinContent(mod + 1, sec + 1, gChanStatus);
0357       if (gChanStatus > repChanBAD)
0358         ++nGoodCh;
0359     }
0360   hReport->Fill(float(nGoodCh) / N_GoodChannels);
0361   return;
0362 }
0363
0364 void CastorDigiMonitor::endRun() {
0365   if (fVerbosity > 0)
0366     LogPrint("CastorDigiMonitor") << "DigiMonitor::endRun: trigger max index = " << TrigIndexMax
0367                                   << "  TriggerIndexies(N):" << endl;
0368   for (int i = 1; i < hpTrigRes->getTProfile()->GetNbinsX(); i++)
0369     if (hpTrigRes->getTProfile()->GetBinContent(i) > 0)
0370       LogPrint("CastorDigiMonitor") << i - 1 << "(" << hpTrigRes->getTProfile()->GetBinContent(i) << ") ";
0371 }
0372
0373 void CastorDigiMonitor::fillTrigRes(edm::Event const &event, const edm::TriggerResults &TrigResults, double Etotal) {
0374   int nTriggers = TrigResults.size();
0375   const edm::TriggerNames &trigName = event.triggerNames(TrigResults);
0376   bool event_triggered = false;
0377   if (nTriggers > 0)
0378     for (int iTrig = 0; iTrig < nTriggers; ++iTrig) {
0379       if (TrigResults.accept(iTrig)) {
0380         int index = trigName.triggerIndex(trigName.triggerName(iTrig));
0381         if (TrigIndexMax < index)
0382           TrigIndexMax = index;
0383         if (fVerbosity > 0)
0384           LogPrint("CastorDigi") << "trigger[" << iTrig << "] name:" << trigName.triggerName(iTrig)
0385                                  << " index= " << index << endl;
0386         hpTrigRes->Fill(index, Etotal);
0387         for (int n = 0; n < int(HltPaths_.size()); n++) {
0388           if (trigName.triggerName(iTrig).find(HltPaths_[n]) != std::string::npos)
0389             event_triggered = true;
0390         }
0391       }  // end if(TrigResults.accept(iTrig)
0392     }
0393
0394   if (event_triggered)
0395     hpBXtrig->Fill(event.bunchCrossing(), Etotal);
0396   return;
0397 }
0398
0399 void CastorDigiMonitor::getDbData(const edm::EventSetup &iSetup) {
0400   edm::ESHandle<CastorChannelQuality> dbChQuality = iSetup.getHandle(castorChannelQualityToken_);
0401   if (fVerbosity > 0) {
0402     LogPrint("CastorDigiMonitor") << " CastorChQuality in CondDB=" << dbChQuality.isValid() << endl;
0403   }
0404
0405   int chInd = 0;
0406   for (int mod = 0; mod < 14; mod++)
0407     for (int sec = 0; sec < 16; sec++)
0408       ChannelStatus[mod][sec] = 0;
0409   std::vector<DetId> channels = dbChQuality->getAllChannels();
0410   N_GoodChannels = 224 - channels.size();
0411   if (fVerbosity > 0)
0412     LogPrint("CastorDigiMonitor") << "CastorDigiMonitor::getDBData: QualityRcdSize=" << channels.size();
0413   for (std::vector<DetId>::iterator ch = channels.begin(); ch != channels.end(); ch++) {
0414     const CastorChannelStatus *quality = dbChQuality->getValues(*ch);
0415     int value = quality->getValue();
0416     int rawId = quality->rawId();
0417     chInd++;
0418     int mod = HcalCastorDetId(*ch).module() - 1;
0419     int sec = HcalCastorDetId(*ch).sector() - 1;
0420     if (mod > 0 && mod < 16 && sec > 0 && sec < 16)
0421       ChannelStatus[mod][sec] = value;
0422     if (fVerbosity > 0)
0423       LogPrint("CastorDigiMonitor") << chInd << " module=" << mod << " sec=" << sec << " rawId=" << rawId
0424                                     << " value=" << value << endl;
0425   }  // end for(std::vector<DetId>::it...
0426   return;
0427 }
0428
0429 int CastorDigiMonitor::ModSecToIndex(int module, int sector) {
0430   int ind = sector + module * 16;
0431   if (ind > 223)
0432     ind = 223;
0433   return (ind);
0434 }