Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-25 05:06:24

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