Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:09:26

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     double sum = 0.;
0228     int err = 0, err2 = 0;
0229     for (int i = 0; i < digi.size(); i++) {
0230       int capid = digi.sample(i).capid();
0231       int dv = digi.sample(i).dv();
0232       int er = digi.sample(i).er();
0233       int rawd = digi.sample(i).adc();
0234       rawd = rawd & 0x7F;
0235       err |= (capid != capid1) | er << 1 | (!dv) << 2;  // =0
0236       err2 += (capid != capid1) | er | (!dv);           // =0
0237       //     if(err !=0) continue;
0238       int ind = ModSecToIndex(module, sector);
0239       h2QmeantsvsCh->Fill(ind, i, LedMonAdc2fc[rawd]);
0240       float q = LedMonAdc2fc[rawd];
0241       Ecell[module][sector] = q;
0242       sum += q;  //     sum += LedMonAdc2fc[rawd];
0243       QrmsTS[ind][i] += (q * q);
0244       QmeanTS[ind][i] += q;
0245       if (err != 0 && fVerbosity > 0)
0246         LogPrint("CastorDigiMonitor") << "event/idigi=" << ievt_ << "/" << i << " cap=cap1_dv_er_err: " << capid << "="
0247                                       << capid1 << " " << dv << " " << er << " " << err;
0248       if (capid1 < 3)
0249         capid1 = capid + 1;
0250       else
0251         capid1 = 0;
0252     }
0253     h2digierr->Fill(module, sector, err);
0254     h2repsum->Fill(module, sector, 1. - err2 / digi.size());
0255   }  // end for(CastorDigiCollection::const_iterator ...
0256 
0257   ievt_++;
0258 
0259   double Etotal = 0.;
0260   for (int sec = 0; sec < 16; sec++)
0261     for (int mod = 0; mod < 14; mod++)
0262       Etotal = Ecell[mod][sec];
0263   hBX->Fill(event.bunchCrossing(), Etotal);
0264   fillTrigRes(event, TrigResults, Etotal);
0265 
0266   for (int sec = 0; sec < 16; sec++) {
0267     float em = Ecell[0][sec] + Ecell[1][sec];
0268     double ehad = 0.;
0269     for (int mod = 2; mod < EtowerLastModule; mod++)
0270       ehad += Ecell[mod][sec];
0271     h2towEMvsHAD->Fill(em, ehad);
0272     htowE->Fill(em + ehad);
0273   }
0274 
0275   const float repChanBAD = 0.9;
0276   const float repChanWarning = 0.95;
0277   if (ievt_ % 100 != 0)
0278     return;
0279 
0280   float ModuleSum[14], SectorSum[16];
0281   for (int m = 0; m < 14; m++)
0282     ModuleSum[m] = 0.;
0283   for (int s = 0; s < 16; s++)
0284     SectorSum[s] = 0.;
0285   for (int mod = 0; mod < 14; mod++)
0286     for (int sec = 0; sec < 16; sec++) {
0287       for (int ts = 0; ts <= 1; ts++) {
0288         int ind = ModSecToIndex(mod, sec);
0289         double Qmean = QmeanTS[ind][ts] / ievt_;
0290         double Qrms = sqrt(QrmsTS[ind][ts] / ievt_ - Qmean * Qmean);
0291         hQIErms[ts]->Fill(Qrms);
0292       }
0293 
0294       double sum = 0.;
0295       for (int ts = 1; ts <= TS_MAX; ts++) {
0296         int ind = ModSecToIndex(mod, sec) + 1;
0297         double a =  //(1) h2QtsvsCh->getTH2D()->GetBinContent(ind,ts);
0298             h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind, ts);
0299         sum += a;
0300         double Qmean = QmeanTS[ind - 1][ts - 1] / ievt_;
0301         double Qrms = QrmsTS[ind - 1][ts - 1] / ievt_ - Qmean * Qmean;
0302         h2QrmsTSvsCh->setBinContent(ind, ts, sqrt(Qrms));
0303       }
0304       ModuleSum[mod] += sum;
0305       SectorSum[sec] += sum;
0306       float isum = float(int(sum * 10. + 0.5)) / 10.;
0307       if (ChannelStatus[mod][sec] != StatusBadChannel)
0308         h2QmeanMap->setBinContent(mod + 1, sec + 1, isum);
0309     }  // end for(int mod=0; mod<14; mod++) for(int sec=0;...
0310 
0311   for (int mod = 0; mod < 14; mod++)
0312     hModule->setBinContent(mod + 1, ModuleSum[mod]);
0313   for (int sec = 0; sec < 16; sec++)
0314     hSector->setBinContent(sec + 1, SectorSum[sec]);
0315 
0316   int nGoodCh = 0;
0317   hTSratio->Reset();
0318   for (int mod = 0; mod < 14; mod++)
0319     for (int sec = 0; sec < 16; sec++) {
0320       if (ChannelStatus[mod][sec] == StatusBadChannel)
0321         continue;
0322       int ind = ModSecToIndex(mod, sec);
0323       double Qmean = QmeanTS[ind][TSped] / ievt_;
0324       double Qrms = QrmsTS[ind][TSped] / ievt_ - Qmean * Qmean;
0325       float ChanStatus = 0.;
0326       if (Qrms < Qrms_DEAD)
0327         ChanStatus = 1.;
0328       h2status->setBinContent(mod + 1, sec + 1, ChanStatus);
0329 
0330       float am = 0.;
0331       for (int ts = 0; ts < TS_MAX - 1; ts++) {
0332         float a = h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 1) +
0333                   h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 2);
0334         if (am < a)
0335           am = a;
0336       }
0337 
0338       double sum = 0.;
0339       for (int ts = 0; ts < TS_MAX; ts++)
0340         sum += h2QmeantsvsCh->getTProfile2D()->GetBinContent(ind + 1, ts + 1);
0341 
0342       float r = 0.;  // worth case - no peak
0343       if (am > 0.)
0344         r = 1. - (sum - am) / (TS_MAX - 2) / am * 2.;
0345       // if(r<0.|| r>1.) cout<<"ievt="<<ievt<<" r="<<r<<" amax= "<<am<<"
0346       // sum="<<sum<<endl;
0347       h2TSratio->setBinContent(mod + 1, sec + 1, r);
0348       hTSratio->Fill(r);
0349 
0350       float statusTS = 1.0;
0351       if (r > RatioThresh1)
0352         statusTS = repChanWarning;
0353       else if (r > 0.99)
0354         statusTS = repChanBAD;
0355       float gChanStatus = statusTS;
0356       if (ChanStatus > 0.)
0357         gChanStatus = repChanBAD;  // RMS
0358       h2qualityMap->setBinContent(mod + 1, sec + 1, gChanStatus);
0359       if (gChanStatus > repChanBAD)
0360         ++nGoodCh;
0361     }
0362   hReport->Fill(float(nGoodCh) / N_GoodChannels);
0363   return;
0364 }
0365 
0366 void CastorDigiMonitor::endRun() {
0367   if (fVerbosity > 0)
0368     LogPrint("CastorDigiMonitor") << "DigiMonitor::endRun: trigger max index = " << TrigIndexMax
0369                                   << "  TriggerIndexies(N):" << endl;
0370   for (int i = 1; i < hpTrigRes->getTProfile()->GetNbinsX(); i++)
0371     if (hpTrigRes->getTProfile()->GetBinContent(i) > 0)
0372       LogPrint("CastorDigiMonitor") << i - 1 << "(" << hpTrigRes->getTProfile()->GetBinContent(i) << ") ";
0373 }
0374 
0375 void CastorDigiMonitor::fillTrigRes(edm::Event const &event, const edm::TriggerResults &TrigResults, double Etotal) {
0376   int nTriggers = TrigResults.size();
0377   const edm::TriggerNames &trigName = event.triggerNames(TrigResults);
0378   bool event_triggered = false;
0379   if (nTriggers > 0)
0380     for (int iTrig = 0; iTrig < nTriggers; ++iTrig) {
0381       if (TrigResults.accept(iTrig)) {
0382         int index = trigName.triggerIndex(trigName.triggerName(iTrig));
0383         if (TrigIndexMax < index)
0384           TrigIndexMax = index;
0385         if (fVerbosity > 0)
0386           LogPrint("CastorDigi") << "trigger[" << iTrig << "] name:" << trigName.triggerName(iTrig)
0387                                  << " index= " << index << endl;
0388         hpTrigRes->Fill(index, Etotal);
0389         for (int n = 0; n < int(HltPaths_.size()); n++) {
0390           if (trigName.triggerName(iTrig).find(HltPaths_[n]) != std::string::npos)
0391             event_triggered = true;
0392         }
0393       }  // end if(TrigResults.accept(iTrig)
0394     }
0395 
0396   if (event_triggered)
0397     hpBXtrig->Fill(event.bunchCrossing(), Etotal);
0398   return;
0399 }
0400 
0401 void CastorDigiMonitor::getDbData(const edm::EventSetup &iSetup) {
0402   edm::ESHandle<CastorChannelQuality> dbChQuality = iSetup.getHandle(castorChannelQualityToken_);
0403   if (fVerbosity > 0) {
0404     LogPrint("CastorDigiMonitor") << " CastorChQuality in CondDB=" << dbChQuality.isValid() << endl;
0405   }
0406 
0407   int chInd = 0;
0408   for (int mod = 0; mod < 14; mod++)
0409     for (int sec = 0; sec < 16; sec++)
0410       ChannelStatus[mod][sec] = 0;
0411   std::vector<DetId> channels = dbChQuality->getAllChannels();
0412   N_GoodChannels = 224 - channels.size();
0413   if (fVerbosity > 0)
0414     LogPrint("CastorDigiMonitor") << "CastorDigiMonitor::getDBData: QualityRcdSize=" << channels.size();
0415   for (std::vector<DetId>::iterator ch = channels.begin(); ch != channels.end(); ch++) {
0416     const CastorChannelStatus *quality = dbChQuality->getValues(*ch);
0417     int value = quality->getValue();
0418     int rawId = quality->rawId();
0419     chInd++;
0420     int mod = HcalCastorDetId(*ch).module() - 1;
0421     int sec = HcalCastorDetId(*ch).sector() - 1;
0422     if (mod > 0 && mod < 16 && sec > 0 && sec < 16)
0423       ChannelStatus[mod][sec] = value;
0424     if (fVerbosity > 0)
0425       LogPrint("CastorDigiMonitor") << chInd << " module=" << mod << " sec=" << sec << " rawId=" << rawId
0426                                     << " value=" << value << endl;
0427   }  // end for(std::vector<DetId>::it...
0428   return;
0429 }
0430 
0431 int CastorDigiMonitor::ModSecToIndex(int module, int sector) {
0432   int ind = sector + module * 16;
0433   if (ind > 223)
0434     ind = 223;
0435   return (ind);
0436 }