EcalMonitorClient/src/OccupancyClient.cc

0001 #include "DQM/EcalMonitorClient/interface/OccupancyClient.h"
0002
0003 #include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
0004
0005 #include "CondFormats/EcalObjects/interface/EcalDQMStatusHelper.h"
0006
0007 #include "DQM/EcalCommon/interface/EcalDQMCommonUtils.h"
0008
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010
0011 namespace ecaldqm {
0012   OccupancyClient::OccupancyClient() : DQWorkerClient(), minHits_(0), deviationThreshold_(0.) {
0013     qualitySummaries_.insert("QualitySummary");
0014   }
0015
0016   void OccupancyClient::setParams(edm::ParameterSet const& _params) {
0017     minHits_ = _params.getUntrackedParameter<int>("minHits");
0018     deviationThreshold_ = _params.getUntrackedParameter<double>("deviationThreshold");
0019   }
0020
0021   void OccupancyClient::producePlots(ProcessType) {
0022     using namespace std;
0023
0024     // number of allowed ieta indices
0025     // EE-: -28 to -1 with -27, -25 empty
0026     // EE+: 1 to 28 with 26, 28 empty
0027     unsigned const nPhiRings(56);
0028
0029     MESet& meQualitySummary(MEs_.at("QualitySummary"));
0030     //     MESet& meHotDigi(MEs_.at("HotDigi"));
0031     //     MESet& meHotRecHitThr(MEs_.at("HotRecHitThr"));
0032     //     MESet& meHotTPDigiThr(MEs_.at("HotTPDigiThr"));
0033
0034     MESet const& sDigi(sources_.at("DigiAll"));
0035     MESet const& sRecHitThr(sources_.at("RecHitThrAll"));
0036     MESet const& sTPDigiThr(sources_.at("TPDigiThrAll"));
0037
0038     uint32_t mask(1 << EcalDQMStatusHelper::PEDESTAL_ONLINE_HIGH_GAIN_RMS_ERROR |
0039                   1 << EcalDQMStatusHelper::PHYSICS_BAD_CHANNEL_WARNING |
0040                   1 << EcalDQMStatusHelper::PHYSICS_BAD_CHANNEL_ERROR);
0041
0042     double digiPhiRingMean[nPhiRings];
0043     std::fill_n(digiPhiRingMean, nPhiRings, 0.);
0044     double rechitPhiRingMean[nPhiRings];
0045     std::fill_n(rechitPhiRingMean, nPhiRings, 0.);
0046     int numCrystals[nPhiRings];  // this is static, but is easier to count now
0047     std::fill_n(numCrystals, nPhiRings, 0);
0048
0049     MESet::const_iterator dEnd(sDigi.end(GetElectronicsMap()));
0050     MESet::const_iterator rItr(GetElectronicsMap(), sRecHitThr);
0051     for (MESet::const_iterator dItr(sDigi.beginChannel(GetElectronicsMap())); dItr != dEnd;
0052          dItr.toNextChannel(GetElectronicsMap())) {
0053       rItr = dItr;
0054
0055       float entries(dItr->getBinContent());
0056       float rhentries(rItr->getBinContent());
0057
0058       DetId id(dItr->getId());
0059       int ieta(0);
0060       if (id.subdetId() == EcalTriggerTower)  // barrel
0061         ieta = EcalTrigTowerDetId(id).ieta();
0062       else {
0063         std::vector<DetId> ids(scConstituents(EcalScDetId(id)));
0064         if (ids.empty())
0065           continue;
0066         ieta = GetTrigTowerMap()->towerOf(ids[0]).ieta();
0067       }
0068
0069       unsigned index(ieta < 0 ? ieta + 28 : ieta + 27);
0070
0071       digiPhiRingMean[index] += entries;
0072       rechitPhiRingMean[index] += rhentries;
0073       numCrystals[index] += 1;
0074     }
0075
0076     for (unsigned ie(0); ie < nPhiRings; ie++) {
0077       digiPhiRingMean[ie] /= numCrystals[ie];
0078       rechitPhiRingMean[ie] /= numCrystals[ie];
0079     }
0080
0081     // Store # of entries for Occupancy analysis
0082     std::vector<float> Nentries(nDCC, 0.);    // digis
0083     std::vector<float> Nrhentries(nDCC, 0.);  // (filtered) rechits
0084
0085     // second round to find hot towers
0086     for (MESet::const_iterator dItr(sDigi.beginChannel(GetElectronicsMap())); dItr != dEnd;
0087          dItr.toNextChannel(GetElectronicsMap())) {
0088       DetId id(dItr->getId());
0089
0090       bool doMask(meQualitySummary.maskMatches(id, mask, statusManager_, GetTrigTowerMap()));
0091
0092       rItr = dItr;
0093
0094       float entries(dItr->getBinContent());
0095       float rhentries(rItr->getBinContent());
0096
0097       int ieta(0);
0098       if (id.subdetId() == EcalTriggerTower)  // barrel
0099         ieta = EcalTrigTowerDetId(id).ieta();
0100       else {
0101         std::vector<DetId> ids(scConstituents(EcalScDetId(id)));
0102         if (ids.empty())
0103           continue;
0104         ieta = GetTrigTowerMap()->towerOf(ids[0]).ieta();
0105       }
0106
0107       unsigned index(ieta < 0 ? ieta + 28 : ieta + 27);
0108
0109       int quality(doMask ? kMGood : kGood);
0110
0111       if (entries > minHits_ && entries > digiPhiRingMean[index] * deviationThreshold_) {
0112         //        meHotDigi->fill(id);
0113         quality = doMask ? kMBad : kBad;
0114       }
0115       if (rhentries > minHits_ && rhentries > rechitPhiRingMean[index] * deviationThreshold_) {
0116         //        meHotRecHitThr->fill(id);
0117         quality = doMask ? kMBad : kBad;
0118       }
0119
0120       meQualitySummary.setBinContent(getEcalDQMSetupObjects(), id, double(quality));
0121
0122       // Keep count of digis & rechits for Occupancy analysis
0123       unsigned iDCC(dccId(id, GetElectronicsMap()) - 1);
0124       if (entries > minHits_)
0125         Nentries[iDCC] += entries;
0126       if (rhentries > minHits_)
0127         Nrhentries[iDCC] += rhentries;
0128     }
0129
0130     double tpdigiPhiRingMean[nPhiRings];
0131     std::fill_n(tpdigiPhiRingMean, nPhiRings, 0.);
0132
0133     for (unsigned iTT(0); iTT < EcalTrigTowerDetId::kSizeForDenseIndexing; ++iTT) {
0134       EcalTrigTowerDetId ttid(EcalTrigTowerDetId::detIdFromDenseIndex(iTT));
0135       float entries(sTPDigiThr.getBinContent(getEcalDQMSetupObjects(), ttid));
0136
0137       unsigned index(ttid.ieta() < 0 ? ttid.ieta() + 28 : ttid.ieta() + 27);
0138
0139       tpdigiPhiRingMean[index] += entries;
0140     }
0141
0142     for (int ie(0); ie < 28; ie++) {
0143       float denom(-1.);
0144       if (ie < 27)
0145         denom = 72.;
0146       else
0147         denom = 36.;
0148       tpdigiPhiRingMean[ie] /= denom;
0149       tpdigiPhiRingMean[55 - ie] /= denom;
0150     }
0151
0152     for (unsigned iTT(0); iTT < EcalTrigTowerDetId::kSizeForDenseIndexing; ++iTT) {
0153       EcalTrigTowerDetId ttid(EcalTrigTowerDetId::detIdFromDenseIndex(iTT));
0154
0155       float entries(sTPDigiThr.getBinContent(getEcalDQMSetupObjects(), ttid));
0156
0157       unsigned index(ttid.ieta() < 0 ? ttid.ieta() + 28 : ttid.ieta() + 27);
0158
0159       int quality(kGood);
0160
0161       if (entries > minHits_ && entries > tpdigiPhiRingMean[index] * deviationThreshold_) {
0162         //        meHotTPDigiThr.fill(ttid);
0163         quality = kBad;
0164       }
0165
0166       if (quality != kBad)
0167         continue;
0168
0169       std::vector<DetId> ids(GetTrigTowerMap()->constituentsOf(ttid));
0170       for (unsigned iD(0); iD < ids.size(); ++iD) {
0171         DetId& id(ids[iD]);
0172
0173         int quality(meQualitySummary.getBinContent(getEcalDQMSetupObjects(), id));
0174         if (quality == kMBad || quality == kBad)
0175           continue;
0176
0177         meQualitySummary.setBinContent(
0178             getEcalDQMSetupObjects(),
0179             id,
0180             meQualitySummary.maskMatches(id, mask, statusManager_, GetTrigTowerMap()) ? kMBad : kBad);
0181       }
0182     }
0183
0184     // Quality check: set entire FED to BAD if its occupancy begins to vanish
0185     // Fill FED statistics from (filtered) RecHit Occupancy
0186     float meanFEDEB(0), meanFEDEE(0), rmsFEDEB(0), rmsFEDEE(0);
0187     unsigned int nFEDEB(0), nFEDEE(0);
0188     for (unsigned iDCC(0); iDCC < nDCC; iDCC++) {
0189       if (iDCC >= kEBmLow && iDCC <= kEBpHigh) {
0190         meanFEDEB += Nrhentries[iDCC];
0191         rmsFEDEB += Nrhentries[iDCC] * Nrhentries[iDCC];
0192         nFEDEB++;
0193       } else {
0194         meanFEDEE += Nrhentries[iDCC];
0195         rmsFEDEE += Nrhentries[iDCC] * Nrhentries[iDCC];
0196         nFEDEE++;
0197       }
0198     }
0199     meanFEDEB /= float(nFEDEB);
0200     rmsFEDEB /= float(nFEDEB);
0201     meanFEDEE /= float(nFEDEE);
0202     rmsFEDEE /= float(nFEDEE);
0203     rmsFEDEB = sqrt(abs(rmsFEDEB - meanFEDEB * meanFEDEB));
0204     rmsFEDEE = sqrt(abs(rmsFEDEE - meanFEDEE * meanFEDEE));
0205     // Analyze FED statistics
0206     float meanFED(0.), rmsFED(0.), nRMS(5.);
0207     for (MESet::iterator qsItr(meQualitySummary.beginChannel(GetElectronicsMap()));
0208          qsItr != meQualitySummary.end(GetElectronicsMap());
0209          qsItr.toNextChannel(GetElectronicsMap())) {
0210       DetId id(qsItr->getId());
0211       unsigned iDCC(dccId(id, GetElectronicsMap()) - 1);
0212       if (iDCC >= kEBmLow && iDCC <= kEBpHigh) {
0213         meanFED = meanFEDEB;
0214         rmsFED = rmsFEDEB;
0215       } else {
0216         meanFED = meanFEDEE;
0217         rmsFED = rmsFEDEE;
0218       }
0219       float threshold(meanFED < nRMS * rmsFED ? minHits_ : meanFED - nRMS * rmsFED);
0220       if (meanFED > 1000. && Nrhentries[iDCC] < threshold)
0221         meQualitySummary.setBinContent(
0222             getEcalDQMSetupObjects(),
0223             id,
0224             meQualitySummary.maskMatches(id, mask, statusManager_, GetTrigTowerMap()) ? kMBad : kBad);
0225     }
0226
0227   }  // producePlots()
0228
0229   DEFINE_ECALDQM_WORKER(OccupancyClient);
0230 }  // namespace ecaldqm