File indexing completed on 2021-10-13 23:03:34
0001
0002
0003
0004
0005
0006
0007 #include "EcalSelectiveReadoutValidation.h"
0008
0009 #include "ecalDccMap.h"
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/ConsumesCollector.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017
0018 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
0019 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0023 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025
0026 #include "DataFormats/Common/interface/Handle.h"
0027
0028 using namespace cms;
0029 using namespace edm;
0030 using namespace std;
0031
0032 const double EcalSelectiveReadoutValidation::rad2deg = 45. / atan(1.);
0033
0034 const int EcalSelectiveReadoutValidation::nDccRus_[nDccs_] = {
0035
0036
0037 34,
0038 32,
0039 33,
0040 33,
0041 32,
0042 34,
0043 33,
0044 34,
0045 33,
0046
0047
0048 68,
0049 68,
0050 68,
0051 68,
0052 68,
0053 68,
0054 68,
0055 68,
0056 68,
0057 68,
0058 68,
0059 68,
0060 68,
0061 68,
0062 68,
0063 68,
0064 68,
0065 68,
0066
0067
0068 68,
0069 68,
0070 68,
0071 68,
0072 68,
0073 68,
0074 68,
0075 68,
0076 68,
0077 68,
0078 68,
0079 68,
0080 68,
0081 68,
0082 68,
0083 68,
0084 68,
0085 68,
0086
0087
0088 32,
0089 33,
0090 33,
0091 32,
0092 34,
0093 33,
0094 34,
0095 33,
0096 34};
0097
0098 EcalSelectiveReadoutValidation::EcalSelectiveReadoutValidation(const ParameterSet& ps)
0099 : geoToken(esConsumes()),
0100 ecalmapping(esConsumes<edm::Transition::BeginRun>()),
0101 hTriggerTowerMap(esConsumes<edm::Transition::BeginRun>()),
0102 physHandle(esConsumes()),
0103 lutGrpHandle(esConsumes()),
0104 lutMapHandle(esConsumes()),
0105 collNotFoundWarn_(ps.getUntrackedParameter<bool>("warnIfCollectionNotFound", true)),
0106 ebDigis_(ps.getParameter<edm::InputTag>("EbDigiCollection"), false, collNotFoundWarn_),
0107 eeDigis_(ps.getParameter<edm::InputTag>("EeDigiCollection"), false, collNotFoundWarn_),
0108 ebNoZsDigis_(ps.getParameter<edm::InputTag>("EbUnsuppressedDigiCollection"), false, false ),
0109 eeNoZsDigis_(ps.getParameter<edm::InputTag>("EeUnsuppressedDigiCollection"), false, false ),
0110 ebSrFlags_(ps.getParameter<edm::InputTag>("EbSrFlagCollection"), false, collNotFoundWarn_),
0111 eeSrFlags_(ps.getParameter<edm::InputTag>("EeSrFlagCollection"), false, collNotFoundWarn_),
0112 ebComputedSrFlags_(
0113 ps.getParameter<edm::InputTag>("EbSrFlagFromTTCollection"), false, false ),
0114 eeComputedSrFlags_(
0115 ps.getParameter<edm::InputTag>("EeSrFlagFromTTCollection"), false, false ),
0116 ebSimHits_(ps.getParameter<edm::InputTag>("EbSimHitCollection"), false, false ),
0117 eeSimHits_(ps.getParameter<edm::InputTag>("EeSimHitCollection"), false, false ),
0118 tps_(ps.getParameter<edm::InputTag>("TrigPrimCollection"), false, collNotFoundWarn_),
0119 ebRecHits_(ps.getParameter<edm::InputTag>("EbRecHitCollection"), false, false ),
0120 eeRecHits_(ps.getParameter<edm::InputTag>("EeRecHitCollection"), false, false ),
0121 fedRaw_(ps.getParameter<edm::InputTag>("FEDRawCollection"), false, false ),
0122 tmax(0),
0123 tmin(numeric_limits<int64_t>::max()),
0124 l1aOfTmin(0),
0125 l1aOfTmax(0),
0126 triggerTowerMap_(nullptr),
0127 localReco_(ps.getParameter<bool>("LocalReco")),
0128 weights_(ps.getParameter<vector<double>>("weights")),
0129 tpInGeV_(ps.getParameter<bool>("tpInGeV")),
0130 firstFIRSample_(ps.getParameter<int>("ecalDccZs1stSample")),
0131 useEventRate_(ps.getParameter<bool>("useEventRate")),
0132 logErrForDccs_(nDccs_, false),
0133 ievt_(0),
0134 allHists_(false),
0135 histDir_(ps.getParameter<string>("histDir")),
0136 withEeSimHit_(false),
0137 withEbSimHit_(false) {
0138 edm::ConsumesCollector collector(consumesCollector());
0139 ebDigis_.setToken(collector);
0140 eeDigis_.setToken(collector);
0141 ebNoZsDigis_.setToken(collector);
0142 eeNoZsDigis_.setToken(collector);
0143 ebSrFlags_.setToken(collector);
0144 eeSrFlags_.setToken(collector);
0145 ebComputedSrFlags_.setToken(collector);
0146 eeComputedSrFlags_.setToken(collector);
0147 ebSimHits_.setToken(collector);
0148 eeSimHits_.setToken(collector);
0149 tps_.setToken(collector);
0150 ebRecHits_.setToken(collector);
0151 eeRecHits_.setToken(collector);
0152 fedRaw_.setToken(collector);
0153
0154 double ebZsThr = ps.getParameter<double>("ebZsThrADCCount");
0155 double eeZsThr = ps.getParameter<double>("eeZsThrADCCount");
0156
0157 ebZsThr_ = lround(ebZsThr * 4);
0158 eeZsThr_ = lround(eeZsThr * 4);
0159
0160
0161 srpAlgoErrorLogFileName_ = ps.getUntrackedParameter<string>("srpAlgoErrorLogFile", "");
0162 logSrpAlgoErrors_ = (!srpAlgoErrorLogFileName_.empty());
0163
0164
0165 srApplicationErrorLogFileName_ = ps.getUntrackedParameter<string>("srApplicationErrorLogFile", "");
0166 logSrApplicationErrors_ = (!srApplicationErrorLogFileName_.empty());
0167
0168
0169 configFirWeights(ps.getParameter<vector<double>>("dccWeights"));
0170
0171
0172 outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
0173
0174 if (!outputFile_.empty()) {
0175 LogInfo("OutputInfo") << " Ecal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
0176 } else {
0177 LogInfo("OutputInfo") << " Ecal Digi Task histograms will NOT be saved";
0178 }
0179
0180
0181 verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0182
0183
0184
0185 vector<string> hists(ps.getUntrackedParameter<vector<string>>("histograms", vector<string>(1, "all")));
0186
0187 for (vector<string>::iterator it = hists.begin(); it != hists.end(); ++it)
0188 histList_.insert(*it);
0189 if (histList_.find("all") != histList_.end())
0190 allHists_ = true;
0191 }
0192
0193 void EcalSelectiveReadoutValidation::updateL1aRate(const edm::Event& event) {
0194 const int32_t bx = event.bunchCrossing();
0195 if (bx < 1 || bx > 3564)
0196 return;
0197
0198 int64_t t = event.bunchCrossing() + (event.orbitNumber() - 1) * 3564;
0199
0200 if (t < tmin) {
0201 tmin = t;
0202 l1aOfTmin = event.id().event();
0203 }
0204
0205 if (t > tmax) {
0206 tmax = t;
0207 l1aOfTmax = event.id().event();
0208 }
0209 }
0210
0211 double EcalSelectiveReadoutValidation::getL1aRate() const {
0212 LogDebug("EcalSrValid") << __FILE__ << ":" << __LINE__ << ": "
0213 << "Tmax = " << tmax << " x 25ns; Tmin = " << tmin << " x 25ns; L1A(Tmax) = " << l1aOfTmax
0214 << "; L1A(Tmin) = " << l1aOfTmin << "\n";
0215 return (double)(l1aOfTmax - l1aOfTmin) / ((tmax - tmin) * 25e-9);
0216 }
0217
0218 void EcalSelectiveReadoutValidation::analyze(Event const& event, EventSetup const& es) {
0219 updateL1aRate(event);
0220
0221
0222 readAllCollections(event);
0223
0224 withEeSimHit_ = (!eeSimHits_->empty());
0225 withEbSimHit_ = (!ebSimHits_->empty());
0226
0227 if (ievt_ < 10) {
0228 edm::LogInfo("EcalSrValid") << "Size of TP collection: " << tps_->size() << std::endl
0229 << "Size of EB SRF collection read from data: " << ebSrFlags_->size() << std::endl
0230 << "Size of EB SRF collection computed from data TTFs: " << ebComputedSrFlags_->size()
0231 << std::endl
0232 << "Size of EE SRF collection read from data: " << eeSrFlags_->size() << std::endl
0233 << "Size of EE SRF collection computed from data TTFs: " << eeComputedSrFlags_->size()
0234 << std::endl;
0235 }
0236
0237 if (ievt_ == 0) {
0238 selectFedsForLog();
0239 }
0240
0241
0242 setTtEtSums(es, *ebNoZsDigis_, *eeNoZsDigis_);
0243
0244
0245 analyzeDataVolume(event, es);
0246
0247
0248
0249
0250 analyzeEB(event, es);
0251
0252
0253
0254
0255 analyzeEE(event, es);
0256
0257 fill(meFullRoCnt_, nEeFROCnt_ + nEbFROCnt_);
0258 fill(meEbFullRoCnt_, nEbFROCnt_);
0259 fill(meEeFullRoCnt_, nEeFROCnt_);
0260
0261 fill(meEbZsErrCnt_, nEbZsErrors_);
0262 fill(meEeZsErrCnt_, nEeZsErrors_);
0263 fill(meZsErrCnt_, nEbZsErrors_ + nEeZsErrors_);
0264
0265 fill(meEbZsErrType1Cnt_, nEbZsErrorsType1_);
0266 fill(meEeZsErrType1Cnt_, nEeZsErrorsType1_);
0267 fill(meZsErrType1Cnt_, nEbZsErrorsType1_ + nEeZsErrorsType1_);
0268
0269
0270 analyzeTP(event, es);
0271
0272 if (!ebComputedSrFlags_->empty()) {
0273 compareSrfColl(event, *ebSrFlags_, *ebComputedSrFlags_);
0274 }
0275 if (!eeComputedSrFlags_->empty()) {
0276 compareSrfColl(event, *eeSrFlags_, *eeComputedSrFlags_);
0277 }
0278 nDroppedFRO_ = 0;
0279 nIncompleteFRO_ = 0;
0280 nCompleteZS_ = 0;
0281 checkSrApplication(event, *ebSrFlags_);
0282 checkSrApplication(event, *eeSrFlags_);
0283 fill(meDroppedFROCnt_, nDroppedFRO_);
0284 fill(meIncompleteFROCnt_, nIncompleteFRO_);
0285 fill(meCompleteZSCnt_, nCompleteZS_);
0286 ++ievt_;
0287 }
0288
0289 void EcalSelectiveReadoutValidation::analyzeEE(const edm::Event& event, const edm::EventSetup& es) {
0290 bool eventError = false;
0291 nEeZsErrors_ = 0;
0292 nEeZsErrorsType1_ = 0;
0293
0294
0295
0296
0297 std::unique_ptr<std::array<std::array<std::array<energiesEe_t, nEeY>, nEeX>, nEndcaps>> eeEnergies =
0298 std::make_unique<std::array<std::array<std::array<energiesEe_t, nEeY>, nEeX>, nEndcaps>>();
0299
0300 for (int iZ0 = 0; iZ0 < nEndcaps; ++iZ0) {
0301 for (int iX0 = 0; iX0 < nEeX; ++iX0) {
0302 for (int iY0 = 0; iY0 < nEeY; ++iY0) {
0303 (*eeEnergies)[iZ0][iX0][iY0].noZsRecE = -numeric_limits<double>::max();
0304 (*eeEnergies)[iZ0][iX0][iY0].recE = -numeric_limits<double>::max();
0305 (*eeEnergies)[iZ0][iX0][iY0].simE = 0;
0306 (*eeEnergies)[iZ0][iX0][iY0].simHit = 0;
0307 (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
0308 }
0309 }
0310 }
0311
0312
0313 auto geoHandle = es.getHandle(geoToken);
0314 const CaloSubdetectorGeometry* geometry_p = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0315
0316
0317
0318 for (unsigned int digis = 0; digis < eeNoZsDigis_->size(); ++digis) {
0319 EEDataFrame frame = (*eeNoZsDigis_)[digis];
0320 int iX0 = iXY2cIndex(frame.id().ix());
0321 int iY0 = iXY2cIndex(frame.id().iy());
0322 int iZ0 = frame.id().zside() > 0 ? 1 : 0;
0323
0324 if (iX0 < 0 || iX0 >= nEeX) {
0325 edm::LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
0326 << "[0," << nEeX - 1 << "]\n";
0327 }
0328 if (iY0 < 0 || iY0 >= nEeY) {
0329 edm::LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
0330 << "[0," << nEeY - 1 << "]\n";
0331 }
0332
0333 (*eeEnergies)[iZ0][iX0][iY0].noZsRecE = frame2Energy(frame);
0334
0335 (*eeEnergies)[iZ0][iX0][iY0].gain12 = true;
0336 for (int i = 0; i < frame.size(); ++i) {
0337 const int gain12Code = 0x1;
0338 if (frame[i].gainId() != gain12Code)
0339 (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
0340 }
0341
0342 const GlobalPoint xtalPos = geometry_p->getGeometry(frame.id())->getPosition();
0343
0344 (*eeEnergies)[iZ0][iX0][iY0].phi = rad2deg * ((double)xtalPos.phi());
0345 (*eeEnergies)[iZ0][iX0][iY0].eta = xtalPos.eta();
0346 }
0347
0348
0349 if (!localReco_) {
0350 for (RecHitCollection::const_iterator it = eeRecHits_->begin(); it != eeRecHits_->end(); ++it) {
0351 const RecHit& hit = *it;
0352 int iX0 = iXY2cIndex(static_cast<const EEDetId&>(hit.id()).ix());
0353 int iY0 = iXY2cIndex(static_cast<const EEDetId&>(hit.id()).iy());
0354 int iZ0 = static_cast<const EEDetId&>(hit.id()).zside() > 0 ? 1 : 0;
0355
0356 if (iX0 < 0 || iX0 >= nEeX) {
0357 LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
0358 << "[0," << nEeX - 1 << "]\n";
0359 }
0360 if (iY0 < 0 || iY0 >= nEeY) {
0361 LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
0362 << "[0," << nEeY - 1 << "]\n";
0363 }
0364
0365 (*eeEnergies)[iZ0][iX0][iY0].recE = hit.energy();
0366 }
0367 }
0368
0369
0370 for (vector<PCaloHit>::const_iterator it = eeSimHits_->begin(); it != eeSimHits_->end(); ++it) {
0371 const PCaloHit& simHit = *it;
0372 EEDetId detId(simHit.id());
0373 int iX = detId.ix();
0374 int iX0 = iXY2cIndex(iX);
0375 int iY = detId.iy();
0376 int iY0 = iXY2cIndex(iY);
0377 int iZ0 = detId.zside() > 0 ? 1 : 0;
0378 (*eeEnergies)[iZ0][iX0][iY0].simE += simHit.energy();
0379 ++(*eeEnergies)[iZ0][iX0][iY0].simHit;
0380 }
0381
0382 bool EEcrystalShot[nEeX][nEeY][2];
0383 pair<int, int> EExtalCoor[nEeX][nEeY][2];
0384
0385 for (int iEeZ = 0; iEeZ < 2; ++iEeZ) {
0386 for (int iEeX = 0; iEeX < nEeX; ++iEeX) {
0387 for (int iEeY = 0; iEeY < nEeY; ++iEeY) {
0388 EEcrystalShot[iEeX][iEeY][iEeZ] = false;
0389 EExtalCoor[iEeX][iEeY][iEeZ] = make_pair(0, 0);
0390 }
0391 }
0392 }
0393
0394
0395 for (EEDigiCollection::const_iterator it = eeDigis_->begin(); it != eeDigis_->end(); ++it) {
0396 const EEDataFrame& frame = *it;
0397 int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
0398 int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
0399 int iZ0 = static_cast<const EEDetId&>(frame.id()).zside() > 0 ? 1 : 0;
0400 if (iX0 < 0 || iX0 >= nEeX) {
0401 LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
0402 << "[0," << nEeX - 1 << "]\n";
0403 }
0404 if (iY0 < 0 || iY0 >= nEeY) {
0405 LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
0406 << "[0," << nEeY - 1 << "]\n";
0407 }
0408
0409 if (!EEcrystalShot[iX0][iY0][iZ0]) {
0410 EEcrystalShot[iX0][iY0][iZ0] = true;
0411 EExtalCoor[iX0][iY0][iZ0] = make_pair(xtalGraphX(frame.id()), xtalGraphY(frame.id()));
0412 } else {
0413 cout << "Error: several digi for same crystal!";
0414 abort();
0415 }
0416
0417 if (localReco_) {
0418 (*eeEnergies)[iZ0][iX0][iY0].recE = frame2Energy(frame);
0419 }
0420
0421 (*eeEnergies)[iZ0][iX0][iY0].gain12 = true;
0422 for (int i = 0; i < frame.size(); ++i) {
0423 const int gain12Code = 0x1;
0424 if (frame[i].gainId() != gain12Code) {
0425 (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
0426 }
0427 }
0428
0429 EESrFlagCollection::const_iterator srf = eeSrFlags_->find(readOutUnitOf(frame.id()));
0430
0431 bool highInterest = false;
0432
0433 if (srf == eeSrFlags_->end())
0434 continue;
0435
0436 if (srf != eeSrFlags_->end()) {
0437 highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK) == EcalSrFlag::SRF_FULL);
0438 }
0439
0440 if (highInterest) {
0441 fill(meEeHiZsFir_, dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr));
0442 } else {
0443 int v = dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr);
0444 fill(meEeLiZsFir_, v);
0445 if (v < eeZsThr_) {
0446 eventError = true;
0447 ++nEeZsErrors_;
0448 pair<int, int> ru = dccCh(frame.id());
0449 if (isRuComplete_[ru.first - 1][ru.second - 1])
0450 ++nEeZsErrorsType1_;
0451 if (nEeZsErrors_ < 3) {
0452 srApplicationErrorLog_ << event.id() << ", "
0453 << "RU " << frame.id() << ", "
0454 << "DCC " << ru.first << " Ch : " << ru.second << ": "
0455 << "LI channel under ZS threshold.\n";
0456 }
0457 if (nEeZsErrors_ == 3) {
0458 srApplicationErrorLog_ << event.id() << ": "
0459 << "more ZS errors for this event...\n";
0460 }
0461 }
0462 }
0463 }
0464
0465 for (int iEeZ = 0; iEeZ < 2; ++iEeZ) {
0466 for (int iEeX = 0; iEeX < nEeX; ++iEeX) {
0467 for (int iEeY = 0; iEeY < nEeY; ++iEeY) {
0468 fill(meChOcc_,
0469 EExtalCoor[iEeX][iEeY][iEeZ].first,
0470 EExtalCoor[iEeX][iEeY][iEeZ].second,
0471 EEcrystalShot[iEeX][iEeY][iEeZ] ? 1 : 0);
0472 }
0473 }
0474 }
0475
0476 for (int iZ0 = 0; iZ0 < nEndcaps; ++iZ0) {
0477 for (int iX0 = 0; iX0 < nEeX; ++iX0) {
0478 for (int iY0 = 0; iY0 < nEeY; ++iY0) {
0479 double recE = (*eeEnergies)[iZ0][iX0][iY0].recE;
0480 if (recE == -numeric_limits<double>::max())
0481 continue;
0482 fill(meEeRecE_, (*eeEnergies)[iZ0][iX0][iY0].recE);
0483
0484 fill(meEeEMean_, ievt_ + 1, (*eeEnergies)[iZ0][iX0][iY0].recE);
0485
0486 if (withEeSimHit_) {
0487 if (!(*eeEnergies)[iZ0][iX0][iY0].simHit) {
0488 fill(meEeNoise_, (*eeEnergies)[iZ0][iX0][iY0].noZsRecE);
0489 } else {
0490 fill(meEeSimE_, (*eeEnergies)[iZ0][iX0][iY0].simE);
0491 fill(meEeRecEHitXtal_, (*eeEnergies)[iZ0][iX0][iY0].recE);
0492 }
0493 fill(meEeRecVsSimE_, (*eeEnergies)[iZ0][iX0][iY0].simE, (*eeEnergies)[iZ0][iX0][iY0].recE);
0494 fill(meEeNoZsRecVsSimE_, (*eeEnergies)[iZ0][iX0][iY0].simE, (*eeEnergies)[iZ0][iX0][iY0].noZsRecE);
0495 }
0496 }
0497 }
0498 }
0499
0500 int EEZs1RuCount[2][20][20];
0501 int EEFullRuCount[2][20][20];
0502 int EEForcedRuCount[2][20][20];
0503 for (int iZ(0); iZ < 2; iZ++) {
0504 for (int iX(0); iX < 20; iX++) {
0505 for (int iY(0); iY < 20; iY++) {
0506 EEZs1RuCount[iZ][iX][iY] = 0;
0507 EEFullRuCount[iZ][iX][iY] = 0;
0508 EEForcedRuCount[iZ][iX][iY] = 0;
0509 }
0510 }
0511 }
0512
0513 nEeFROCnt_ = 0;
0514 char eeSrfMark[2][20][20];
0515 bzero(eeSrfMark, sizeof(eeSrfMark));
0516
0517 for (EESrFlagCollection::const_iterator it = eeSrFlags_->begin(); it != eeSrFlags_->end(); ++it) {
0518 const EESrFlag& srf = *it;
0519
0520 int iX = srf.id().ix();
0521 int iY = srf.id().iy();
0522 int zside = srf.id().zside();
0523 if (iX < 1 || iY > 100)
0524 throw cms::Exception("EcalSelectiveReadoutValidation")
0525 << "Found an endcap SRF with an invalid det ID: " << srf.id() << ".\n";
0526 ++eeSrfMark[zside > 0 ? 1 : 0][iX - 1][iY - 1];
0527 if (eeSrfMark[zside > 0 ? 1 : 0][iX - 1][iY - 1] > 1)
0528 throw cms::Exception("EcalSelectiveReadoutValidation") << "Duplicate SRF for supercrystal " << srf.id() << ".\n";
0529 int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
0530 if (flag == EcalSrFlag::SRF_ZS1) {
0531 EEZs1RuCount[zside > 0 ? 1 : 0][iX - 1][iY - 1] += 1;
0532 }
0533
0534 if (flag == EcalSrFlag::SRF_FULL) {
0535 EEFullRuCount[zside > 0 ? 1 : 0][iX - 1][iY - 1] += 1;
0536 ++nEeFROCnt_;
0537 }
0538
0539 if (srf.value() & EcalSrFlag::SRF_FORCED_MASK) {
0540 EEForcedRuCount[zside > 0 ? 1 : 0][iX - 1][iY - 1] += 1;
0541 }
0542 }
0543 for (int iZ(0); iZ < 2; iZ++) {
0544 int xOffset(iZ == 0 ? -40 : 20);
0545 for (int iX(0); iX < 20; iX++) {
0546 for (int iY(0); iY < 20; iY++) {
0547 int GraphX = (iX + 1) + xOffset;
0548 int GraphY = (iY + 1);
0549 fill(meZs1Ru_, GraphX, GraphY, EEZs1RuCount[iZ][iX][iY]);
0550 fill(meFullRoRu_, GraphX, GraphY, EEFullRuCount[iZ][iX][iY]);
0551 fill(meForcedRu_, GraphX, GraphY, EEForcedRuCount[iZ][iX][iY]);
0552 }
0553 }
0554 }
0555
0556 if (eventError)
0557 srApplicationErrorLog_ << event.id() << ": " << nEeZsErrors_
0558 << " ZS-flagged EE channels under "
0559 "the ZS threshold, whose "
0560 << nEeZsErrorsType1_ << " in a complete RU.\n";
0561 }
0562
0563 void EcalSelectiveReadoutValidation::analyzeEB(const edm::Event& event, const edm::EventSetup& es) {
0564 bool eventError = false;
0565 nEbZsErrors_ = 0;
0566 nEbZsErrorsType1_ = 0;
0567 vector<pair<int, int>> xtalEtaPhi;
0568
0569 xtalEtaPhi.reserve(nEbPhi * nEbEta);
0570
0571
0572
0573
0574 std::unique_ptr<std::array<std::array<energiesEb_t, nEbPhi>, nEbEta>> ebEnergies =
0575 std::make_unique<std::array<std::array<energiesEb_t, nEbPhi>, nEbEta>>();
0576
0577 for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
0578 for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0) {
0579 (*ebEnergies)[iEta0][iPhi0].noZsRecE = -numeric_limits<double>::max();
0580 (*ebEnergies)[iEta0][iPhi0].recE = -numeric_limits<double>::max();
0581 (*ebEnergies)[iEta0][iPhi0].simE = 0;
0582 (*ebEnergies)[iEta0][iPhi0].simHit = 0;
0583 (*ebEnergies)[iEta0][iPhi0].gain12 = false;
0584 xtalEtaPhi.push_back(pair<int, int>(iEta0, iPhi0));
0585 }
0586 }
0587
0588
0589 auto geoHandle = es.getHandle(geoToken);
0590 const CaloSubdetectorGeometry* geometry_p = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0591
0592
0593
0594 for (EBDigiCollection::const_iterator it = ebNoZsDigis_->begin(); it != ebNoZsDigis_->end(); ++it) {
0595 const EBDataFrame& frame = *it;
0596 int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(frame.id()).ieta());
0597 int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(frame.id()).iphi());
0598 if (iEta0 < 0 || iEta0 >= nEbEta) {
0599 stringstream s;
0600 s << "EcalSelectiveReadoutValidation: "
0601 << "iEta0 (= " << iEta0 << ") is out of range ("
0602 << "[0," << nEbEta - 1 << "]\n";
0603 throw cms::Exception(s.str());
0604 }
0605 if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
0606 stringstream s;
0607 s << "EcalSelectiveReadoutValidation: "
0608 << "iPhi0 (= " << iPhi0 << ") is out of range ("
0609 << "[0," << nEbPhi - 1 << "]\n";
0610 throw cms::Exception(s.str());
0611 }
0612
0613 (*ebEnergies)[iEta0][iPhi0].noZsRecE = frame2Energy(frame);
0614 (*ebEnergies)[iEta0][iPhi0].gain12 = true;
0615 for (int i = 0; i < frame.size(); ++i) {
0616 const int gain12Code = 0x1;
0617 if (frame[i].gainId() != gain12Code)
0618 (*ebEnergies)[iEta0][iPhi0].gain12 = false;
0619 }
0620
0621 const GlobalPoint xtalPos = geometry_p->getGeometry(frame.id())->getPosition();
0622
0623 (*ebEnergies)[iEta0][iPhi0].phi = rad2deg * ((double)xtalPos.phi());
0624 (*ebEnergies)[iEta0][iPhi0].eta = xtalPos.eta();
0625 }
0626
0627
0628 for (vector<PCaloHit>::const_iterator it = ebSimHits_->begin(); it != ebSimHits_->end(); ++it) {
0629 const PCaloHit& simHit = *it;
0630 EBDetId detId(simHit.id());
0631 int iEta = detId.ieta();
0632 int iEta0 = iEta2cIndex(iEta);
0633 int iPhi = detId.iphi();
0634 int iPhi0 = iPhi2cIndex(iPhi);
0635 (*ebEnergies)[iEta0][iPhi0].simE += simHit.energy();
0636 ++(*ebEnergies)[iEta0][iPhi0].simHit;
0637 }
0638
0639 bool crystalShot[nEbEta][nEbPhi];
0640 pair<int, int> EBxtalCoor[nEbEta][nEbPhi];
0641
0642 for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
0643 for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0) {
0644 crystalShot[iEta0][iPhi0] = false;
0645 EBxtalCoor[iEta0][iPhi0] = make_pair(0, 0);
0646 }
0647 }
0648
0649 int nEbDigi = 0;
0650
0651 for (EBDigiCollection::const_iterator it = ebDigis_->begin(); it != ebDigis_->end(); ++it) {
0652 ++nEbDigi;
0653 const EBDataFrame& frame = *it;
0654 int iEta = static_cast<const EBDetId&>(frame.id()).ieta();
0655 int iPhi = static_cast<const EBDetId&>(frame.id()).iphi();
0656 int iEta0 = iEta2cIndex(iEta);
0657 int iPhi0 = iPhi2cIndex(iPhi);
0658 if (iEta0 < 0 || iEta0 >= nEbEta) {
0659 throw(cms::Exception("EcalSelectiveReadoutValidation") << "iEta0 (= " << iEta0 << ") is out of range ("
0660 << "[0," << nEbEta - 1 << "]");
0661 }
0662 if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
0663 throw(cms::Exception("EcalSelectiveReadoutValidation") << "iPhi0 (= " << iPhi0 << ") is out of range ("
0664 << "[0," << nEbPhi - 1 << "]");
0665 }
0666 assert(iEta0 >= 0 && iEta0 < nEbEta);
0667 assert(iPhi0 >= 0 && iPhi0 < nEbPhi);
0668 if (!crystalShot[iEta0][iPhi0]) {
0669 crystalShot[iEta0][iPhi0] = true;
0670 EBxtalCoor[iEta0][iPhi0] = make_pair(xtalGraphX(frame.id()), xtalGraphY(frame.id()));
0671 } else {
0672 cout << "Error: several digi for same crystal!";
0673 abort();
0674 }
0675 if (localReco_) {
0676 (*ebEnergies)[iEta0][iPhi0].recE = frame2Energy(frame);
0677 }
0678
0679 (*ebEnergies)[iEta0][iPhi0].gain12 = true;
0680 for (int i = 0; i < frame.size(); ++i) {
0681 const int gain12Code = 0x1;
0682 if (frame[i].gainId() != gain12Code) {
0683 (*ebEnergies)[iEta0][iPhi0].gain12 = false;
0684 }
0685 }
0686
0687 EBSrFlagCollection::const_iterator srf = ebSrFlags_->find(readOutUnitOf(frame.id()));
0688
0689 bool highInterest = false;
0690
0691
0692
0693
0694
0695
0696 if (srf != ebSrFlags_->end()) {
0697 highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK) == EcalSrFlag::SRF_FULL);
0698 }
0699
0700 if (highInterest) {
0701 fill(meEbHiZsFir_, dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr));
0702 } else {
0703 int v = dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr);
0704 fill(meEbLiZsFir_, v);
0705 if (v < ebZsThr_) {
0706 eventError = true;
0707 ++nEbZsErrors_;
0708 pair<int, int> ru = dccCh(frame.id());
0709 if (isRuComplete_[ru.first - 1][ru.second - 1])
0710 ++nEbZsErrorsType1_;
0711 if (nEbZsErrors_ < 3) {
0712 srApplicationErrorLog_ << event.id() << ", "
0713 << "RU " << frame.id() << ", "
0714 << "DCC " << ru.first << " Ch : " << ru.second << ": "
0715 << "LI channel under ZS threshold.\n";
0716 }
0717 if (nEbZsErrors_ == 3) {
0718 srApplicationErrorLog_ << event.id() << ": "
0719 << "more ZS errors for this event...\n";
0720 }
0721 }
0722 }
0723 }
0724
0725 for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
0726 for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0)
0727 fill(meChOcc_,
0728 EBxtalCoor[iEta0][iPhi0].first,
0729 EBxtalCoor[iEta0][iPhi0].second,
0730 crystalShot[iEta0][iPhi0] ? 1. : 0.);
0731 }
0732
0733 if (!localReco_) {
0734 for (RecHitCollection::const_iterator it = ebRecHits_->begin(); it != ebRecHits_->end(); ++it) {
0735 ++nEbDigi;
0736 const RecHit& hit = *it;
0737 int iEta = static_cast<const EBDetId&>(hit.id()).ieta();
0738 int iPhi = static_cast<const EBDetId&>(hit.id()).iphi();
0739 int iEta0 = iEta2cIndex(iEta);
0740 int iPhi0 = iPhi2cIndex(iPhi);
0741 if (iEta0 < 0 || iEta0 >= nEbEta) {
0742 LogError("EcalSrValid") << "iEta0 (= " << iEta0 << ") is out of range ("
0743 << "[0," << nEbEta - 1 << "]\n";
0744 }
0745 if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
0746 LogError("EcalSrValid") << "iPhi0 (= " << iPhi0 << ") is out of range ("
0747 << "[0," << nEbPhi - 1 << "]\n";
0748 }
0749 (*ebEnergies)[iEta0][iPhi0].recE = hit.energy();
0750 }
0751 }
0752
0753 for (unsigned int i = 0; i < xtalEtaPhi.size(); ++i) {
0754 int iEta0 = xtalEtaPhi[i].first;
0755 int iPhi0 = xtalEtaPhi[i].second;
0756 energiesEb_t& energies = (*ebEnergies)[iEta0][iPhi0];
0757
0758 double recE = energies.recE;
0759 if (recE != -numeric_limits<double>::max()) {
0760 fill(meEbRecE_, (*ebEnergies)[iEta0][iPhi0].recE);
0761 fill(meEbEMean_, ievt_ + 1, recE);
0762 }
0763
0764 if (withEbSimHit_) {
0765 if (!energies.simHit) {
0766 fill(meEbNoise_, energies.noZsRecE);
0767 } else {
0768 fill(meEbSimE_, energies.simE);
0769 fill(meEbRecEHitXtal_, energies.recE);
0770 }
0771 fill(meEbRecVsSimE_, energies.simE, energies.recE);
0772 fill(meEbNoZsRecVsSimE_, energies.simE, energies.noZsRecE);
0773 }
0774 }
0775
0776 int EBZs1RuCount[2][17][72];
0777 int EBFullRuCount[2][17][72];
0778 int EBForcedRuCount[2][17][72];
0779 std::pair<int, int> EBtowerCoor[2][17][72];
0780 for (int iZ(0); iZ < 2; iZ++) {
0781 for (int iEta(0); iEta < 17; iEta++) {
0782 for (int iPhi(0); iPhi < 72; iPhi++) {
0783 EBZs1RuCount[iZ][iEta][iPhi] = 0;
0784 EBFullRuCount[iZ][iEta][iPhi] = 0;
0785 EBForcedRuCount[iZ][iEta][iPhi] = 0;
0786 }
0787 }
0788 }
0789
0790
0791 nEbFROCnt_ = 0;
0792 char ebSrfMark[2][17][72];
0793 bzero(ebSrfMark, sizeof(ebSrfMark));
0794
0795 for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
0796 const EBSrFlag& srf = *it;
0797 int iEtaAbs = srf.id().ietaAbs();
0798 int iPhi = srf.id().iphi();
0799 int iZ = srf.id().zside();
0800
0801
0802
0803
0804 if (iEtaAbs < 1 || iEtaAbs > 17 || iPhi < 1 || iPhi > 72)
0805 throw cms::Exception("EcalSelectiveReadoutValidation")
0806 << "Found a barrel SRF with an invalid det ID: " << srf.id() << ".\n";
0807 ++ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1];
0808 if (ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] > 1)
0809 throw cms::Exception("EcalSelectiveReadoutValidation") << "Duplicate SRF for RU " << srf.id() << ".\n";
0810
0811 EBtowerCoor[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] = std::pair<int, int>(srf.id().ieta(), srf.id().iphi());
0812
0813 int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
0814 if (flag == EcalSrFlag::SRF_ZS1) {
0815 EBZs1RuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
0816 }
0817 if (flag == EcalSrFlag::SRF_FULL) {
0818 EBFullRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
0819 ++nEbFROCnt_;
0820 }
0821 if (srf.value() & EcalSrFlag::SRF_FORCED_MASK) {
0822 EBForcedRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
0823 }
0824 }
0825 for (int iZ(0); iZ < 2; iZ++) {
0826 for (int iEta(0); iEta < 17; iEta++) {
0827 for (int iPhi(0); iPhi < 72; iPhi++) {
0828 float x(EBtowerCoor[iZ][iEta][iPhi].first);
0829 float y(EBtowerCoor[iZ][iEta][iPhi].second);
0830 fill(meZs1Ru_, x, y, EBZs1RuCount[iZ][iEta][iPhi]);
0831 fill(meFullRoRu_, x, y, EBFullRuCount[iZ][iEta][iPhi]);
0832 fill(meForcedRu_, x, y, EBForcedRuCount[iZ][iEta][iPhi]);
0833 }
0834 }
0835 }
0836
0837 if (eventError)
0838 srApplicationErrorLog_ << event.id() << ": " << nEbZsErrors_
0839 << " ZS-flagged EB channels under "
0840 "the ZS threshold, whose "
0841 << nEbZsErrorsType1_ << " in a complete RU.\n";
0842 }
0843
0844 EcalSelectiveReadoutValidation::~EcalSelectiveReadoutValidation() {}
0845
0846 void EcalSelectiveReadoutValidation::dqmBeginRun(edm::Run const& r, edm::EventSetup const& es) {
0847
0848 triggerTowerMap_ = &es.getData(hTriggerTowerMap);
0849
0850
0851 elecMap_ = &es.getData(ecalmapping);
0852
0853 initAsciiFile();
0854 }
0855
0856 void EcalSelectiveReadoutValidation::dqmEndRun(const edm::Run& r, const edm::EventSetup& es) {
0857 meL1aRate_->Fill(getL1aRate());
0858 }
0859
0860 void EcalSelectiveReadoutValidation::bookHistograms(DQMStore::IBooker& ibooker,
0861 edm::Run const&,
0862 edm::EventSetup const&) {
0863 ibooker.setCurrentFolder("EcalDigisV/SelectiveReadout");
0864
0865 {
0866 auto scope = DQMStore::IBooker::UseRunScope(ibooker);
0867 meL1aRate_ = bookFloat(ibooker, "l1aRate_");
0868 }
0869
0870 meDccVol_ = bookProfile(ibooker,
0871 "hDccVol",
0872 "ECAL DCC event fragment size;Dcc id; "
0873 "<Event size> (kB)",
0874 nDccs_,
0875 .5,
0876 .5 + nDccs_);
0877
0878 meDccLiVol_ = bookProfile(ibooker,
0879 "hDccLiVol",
0880 "LI channel payload per DCC;Dcc id; "
0881 "<Event size> (kB)",
0882 nDccs_,
0883 .5,
0884 .5 + nDccs_);
0885
0886 meDccHiVol_ = bookProfile(ibooker,
0887 "hDccHiVol",
0888 "HI channel payload per DCC;Dcc id; "
0889 "<Event size> (kB)",
0890 nDccs_,
0891 .5,
0892 .5 + nDccs_);
0893
0894 meDccVolFromData_ = bookProfile(ibooker,
0895 "hDccVolFromData",
0896 "ECAL DCC event fragment size;Dcc id; "
0897 "<Event size> (kB)",
0898 nDccs_,
0899 .5,
0900 .5 + nDccs_);
0901
0902 meVolBLI_ = book1D(ibooker,
0903 "hVolBLI",
0904 "ECAL Barrel low interest crystal data payload;"
0905 "Event size (kB);Nevts",
0906 100,
0907 0.,
0908 200.);
0909
0910 meVolELI_ = book1D(ibooker,
0911 "hVolELI",
0912 "Endcap low interest crystal data payload;"
0913 "Event size (kB);Nevts",
0914 100,
0915 0.,
0916 200.);
0917
0918 meVolLI_ = book1D(ibooker,
0919 "hVolLI",
0920 "ECAL low interest crystal data payload;"
0921 "Event size (kB);Nevts",
0922 100,
0923 0.,
0924 200.);
0925
0926 meVolBHI_ = book1D(ibooker,
0927 "hVolBHI",
0928 "Barrel high interest crystal data payload;"
0929 "Event size (kB);Nevts",
0930 100,
0931 0.,
0932 200.);
0933
0934 meVolEHI_ = book1D(ibooker,
0935 "hVolEHI",
0936 "Endcap high interest crystal data payload;"
0937 "Event size (kB);Nevts",
0938 100,
0939 0.,
0940 200.);
0941
0942 meVolHI_ = book1D(ibooker,
0943 "hVolHI",
0944 "ECAL high interest crystal data payload;"
0945 "Event size (kB);Nevts",
0946 100,
0947 0.,
0948 200.);
0949
0950 meVolB_ = book1D(ibooker,
0951 "hVolB",
0952 "Barrel data volume;Event size (kB);Nevts",
0953 100,
0954 0.,
0955 200.);
0956
0957 meVolE_ = book1D(ibooker,
0958 "hVolE",
0959 "Endcap data volume;Event size (kB);Nevts",
0960 100,
0961 0.,
0962 200.);
0963
0964 meVol_ = book1D(ibooker,
0965 "hVol",
0966 "ECAL data volume;Event size (kB);Nevts",
0967 100,
0968 0.,
0969 200.);
0970
0971 meChOcc_ = bookProfile2D(ibooker,
0972 "h2ChOcc",
0973 "ECAL crystal channel occupancy after zero suppression;"
0974 "iX -200 / iEta / iX + 100;"
0975 "iY / iPhi (starting from -10^{o}!);"
0976 "Event count rate",
0977 401,
0978 -200.5,
0979 200.5,
0980 360,
0981 .5,
0982 360.5);
0983
0984
0985 string tpUnit;
0986 if (tpInGeV_)
0987 tpUnit = string("GeV");
0988 else
0989 tpUnit = string("TP hw unit");
0990 string title;
0991 title = string("Trigger primitive TT E_{T};E_{T} ") + tpUnit + string(";Event Count");
0992 meTp_ = bookProfile(ibooker,
0993 "hTp",
0994 title,
0995 (tpInGeV_ ? 100 : 40),
0996 0.,
0997 (tpInGeV_ ? 10. : 40.));
0998
0999 meTtf_ = bookProfile(ibooker,
1000 "hTtf",
1001 "Trigger primitive TT flag;Flag number;Event count",
1002 8,
1003 -.5,
1004 7.5);
1005
1006 title = string("Trigger tower flag vs TP;E_{T}(TT) (") + tpUnit + string(");Flag number");
1007 meTtfVsTp_ = book2D(ibooker, "h2TtfVsTp", title, 100, 0., (tpInGeV_ ? 10. : 40.), 8, -.5, 7.5);
1008
1009 meTtfVsEtSum_ = book2D(ibooker,
1010 "h2TtfVsEtSum",
1011 "Trigger tower flag vs #sumE_{T};"
1012 "E_{T}(TT) (GeV);"
1013 "TTF",
1014 100,
1015 0.,
1016 10.,
1017 8,
1018 -.5,
1019 7.5);
1020 title = string(
1021 "Trigger primitive Et (TP) vs #sumE_{T};"
1022 "E_{T} (sum) (GeV);"
1023 "E_{T} (TP) (") +
1024 tpUnit + string(")");
1025
1026 meTpVsEtSum_ = book2D(ibooker, "h2TpVsEtSum", title, 100, 0., 10., 100, 0., (tpInGeV_ ? 10. : 40.));
1027
1028 title = string(
1029 "Trigger primitive E_{T};"
1030 "iEta;"
1031 "iPhi;"
1032 "E_{T} (TP) (") +
1033 tpUnit + string(")");
1034 meTpMap_ = bookProfile2D(ibooker, "h2Tp", title, 57, -28.5, 28.5, 72, .5, 72.5);
1035
1036
1037 meFullRoRu_ = book2D(ibooker,
1038 "h2FRORu",
1039 "Full Read-out readout unit;"
1040 "iX - 40 / iEta / iX + 20;"
1041 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1042 "Event count",
1043 80,
1044 -39.5,
1045 40.5,
1046 72,
1047 .5,
1048 72.5);
1049
1050 meFullRoCnt_ = book1D(ibooker,
1051 "hFROCnt",
1052 "Number of Full-readout-flagged readout units;"
1053 "FRO RU count;Event count",
1054 300,
1055 -.5,
1056 299.5);
1057
1058 meEbFullRoCnt_ = book1D(ibooker,
1059 "hEbFROCnt",
1060 "Number of EB Full-readout-flagged readout units;"
1061 "FRO RU count;Event count",
1062 200,
1063 -.5,
1064 199.5);
1065
1066 meEeFullRoCnt_ = book1D(ibooker,
1067 "hEeFROCnt",
1068 "Number of EE Full-readout-flagged readout units;"
1069 "FRO RU count;Event count",
1070 200,
1071 -.5,
1072 199.5);
1073
1074 meZs1Ru_ = book2D(ibooker,
1075 "h2Zs1Ru",
1076 "Readout unit with ZS-thr-1 flag;"
1077 "iX - 40 / iEta / iX + 20;"
1078 "iY0 / iPhi0 (iPhi = 1 at phi = 0 rad);"
1079 "Event count",
1080 80,
1081 -39.5,
1082 40.5,
1083 72,
1084 .5,
1085 72.5);
1086
1087 meForcedRu_ = book2D(ibooker,
1088 "h2ForcedRu",
1089 "ECAL readout unit with forced bit of SR flag on;"
1090 "iX - 40 / iEta / iX + 20;"
1091 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1092 "Event count",
1093 80,
1094 -39.5,
1095 40.5,
1096 72,
1097 .5,
1098 72.5);
1099
1100 meLiTtf_ = bookProfile2D(ibooker,
1101 "h2LiTtf",
1102 "Low interest trigger tower flags;"
1103 "iEta;"
1104 "iPhi;"
1105 "Event count",
1106 57,
1107 -28.5,
1108 28.5,
1109 72,
1110 .5,
1111 72.5);
1112
1113 meMiTtf_ = bookProfile2D(ibooker,
1114 "h2MiTtf",
1115 "Mid interest trigger tower flags;"
1116 "iEta;"
1117 "iPhi;"
1118 "Event count",
1119 57,
1120 -28.5,
1121 28.5,
1122 72,
1123 .5,
1124 72.5);
1125
1126 meHiTtf_ = bookProfile2D(ibooker,
1127 "h2HiTtf",
1128 "High interest trigger tower flags;"
1129 "iEta;"
1130 "iPhi;"
1131 "Event count",
1132 57,
1133 -28.5,
1134 28.5,
1135 72,
1136 .5,
1137 72.5);
1138
1139 meForcedTtf_ = book2D(ibooker,
1140 "h2ForcedTtf",
1141 "Trigger tower flags with forced bit set;"
1142 "iEta;"
1143 "iPhi;"
1144 "Event count",
1145 57,
1146 -28.5,
1147 28.5,
1148 72,
1149 .5,
1150 72.5);
1151
1152 const float ebMinNoise = -1.;
1153 const float ebMaxNoise = 1.;
1154
1155 const float eeMinNoise = -1.;
1156 const float eeMaxNoise = 1.;
1157
1158 const float ebMinE = ebMinNoise;
1159 const float ebMaxE = ebMaxNoise;
1160
1161 const float eeMinE = eeMinNoise;
1162 const float eeMaxE = eeMaxNoise;
1163
1164 const int evtMax = 500;
1165
1166 meEbRecE_ = book1D(ibooker, "hEbRecE", "Crystal reconstructed energy;E (GeV);Event count", 100, ebMinE, ebMaxE);
1167
1168 meEbEMean_ = bookProfile(ibooker, "hEbEMean", "EE <E_hit>;event #;<E_hit> (GeV)", evtMax, .5, evtMax + .5);
1169
1170 meEbNoise_ = book1D(ibooker,
1171 "hEbNoise",
1172 "Crystal noise "
1173 "(rec E of crystal without deposited energy)"
1174 ";Rec E (GeV);Event count",
1175 100,
1176 ebMinNoise,
1177 ebMaxNoise);
1178
1179 meEbLiZsFir_ = book1D(ibooker,
1180 "zsEbLiFIRemu",
1181 "Emulated ouput of ZS FIR filter for EB "
1182 "low interest crystals;"
1183 "ADC count*4;"
1184 "Event count",
1185 60,
1186 -30,
1187 30);
1188
1189 meEbHiZsFir_ = book1D(ibooker,
1190 "zsEbHiFIRemu",
1191 "Emulated ouput of ZS FIR filter for EB "
1192 "high interest crystals;"
1193 "ADC count*4;"
1194 "Event count",
1195 60,
1196 -30,
1197 30);
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207 meEbSimE_ = book1D(ibooker, "hEbSimE", "EB hit crystal simulated energy", 100, ebMinE, ebMaxE);
1208
1209 meEbRecEHitXtal_ = book1D(ibooker, "hEbRecEHitXtal", "EB rec energy of hit crystals", 100, ebMinE, ebMaxE);
1210
1211 meEbRecVsSimE_ = book2D(ibooker,
1212 "hEbRecVsSimE",
1213 "Crystal simulated vs reconstructed energy;"
1214 "Esim (GeV);Erec GeV);Event count",
1215 100,
1216 ebMinE,
1217 ebMaxE,
1218 100,
1219 ebMinE,
1220 ebMaxE);
1221
1222 meEbNoZsRecVsSimE_ = book2D(ibooker,
1223 "hEbNoZsRecVsSimE",
1224 "Crystal no-zs simulated vs reconstructed "
1225 "energy;"
1226 "Esim (GeV);Erec GeV);Event count",
1227 100,
1228 ebMinE,
1229 ebMaxE,
1230 100,
1231 ebMinE,
1232 ebMaxE);
1233
1234 meEeRecE_ = book1D(ibooker,
1235 "hEeRecE",
1236 "EE crystal reconstructed energy;E (GeV);"
1237 "Event count",
1238 100,
1239 eeMinE,
1240 eeMaxE);
1241
1242 meEeEMean_ = bookProfile(ibooker, "hEeEMean", "<E_{EE hit}>;event;<E_{hit}> (GeV)", evtMax, .5, evtMax + .5);
1243
1244 meEeNoise_ = book1D(ibooker,
1245 "hEeNoise",
1246 "EE crystal noise "
1247 "(rec E of crystal without deposited energy);"
1248 "E (GeV);Event count",
1249 200,
1250 eeMinNoise,
1251 eeMaxNoise);
1252
1253 meEeLiZsFir_ = book1D(ibooker,
1254 "zsEeLiFIRemu",
1255 "Emulated ouput of ZS FIR filter for EE "
1256 "low interest crystals;"
1257 "ADC count*4;"
1258 "Event count",
1259 60,
1260 -30,
1261 30);
1262
1263 meEeHiZsFir_ = book1D(ibooker,
1264 "zsEeHiFIRemu",
1265 "Emulated ouput of ZS FIR filter for EE "
1266 "high interest crystals;"
1267 "ADC count*4;"
1268 "Event count",
1269 60,
1270 -30,
1271 30);
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281 meEeSimE_ = book1D(ibooker, "hEeSimE", "EE hit crystal simulated energy", 100, eeMinE, eeMaxE);
1282
1283 meEeRecEHitXtal_ = book1D(ibooker, "hEeRecEHitXtal", "EE rec energy of hit crystals", 100, eeMinE, eeMaxE);
1284
1285 meEeRecVsSimE_ = book2D(ibooker,
1286 "hEeRecVsSimE",
1287 "EE crystal simulated vs reconstructed energy;"
1288 "Esim (GeV);Erec GeV);Event count",
1289 100,
1290 eeMinE,
1291 eeMaxE,
1292 100,
1293 eeMinE,
1294 eeMaxE);
1295
1296 meEeNoZsRecVsSimE_ = book2D(ibooker,
1297 "hEeNoZsRecVsSimE",
1298 "EE crystal no-zs simulated vs "
1299 "reconstructed "
1300 "energy;Esim (GeV);Erec GeV);Event count",
1301 100,
1302 eeMinE,
1303 eeMaxE,
1304 100,
1305 eeMinE,
1306 eeMaxE);
1307
1308 meSRFlagsConsistency_ = book2D(ibooker,
1309 "hSRAlgoErrorMap",
1310 "TTFlags and SR Flags mismatch;"
1311 "iX - 40 / iEta / iX + 20;"
1312 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1313 "Event count",
1314 80,
1315 -39.5,
1316 40.5,
1317 72,
1318 .5,
1319 72.5);
1320
1321
1322 meIncompleteFROMap_ = book2D(ibooker,
1323 "hIncompleteFROMap",
1324 "Incomplete full-readout-flagged readout units;"
1325 "iX - 40 / iEta / iX + 20;"
1326 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1327 "Event count",
1328 80,
1329 -39.5,
1330 40.5,
1331 72,
1332 .5,
1333 72.5);
1334
1335 meIncompleteFROCnt_ = book1D(ibooker,
1336 "hIncompleteFROCnt",
1337 "Number of incomplete full-readout-flagged "
1338 "readout units;"
1339 "Number of RUs;Event count;",
1340 200,
1341 -.5,
1342 199.5);
1343
1344 meIncompleteFRORateMap_ = bookProfile2D(ibooker,
1345 "hIncompleteFRORateMap",
1346 "Incomplete full-readout-flagged readout units;"
1347 "iX - 40 / iEta / iX + 20;"
1348 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1349 "Incomplete error rate",
1350 80,
1351 -39.5,
1352 40.5,
1353 72,
1354 .5,
1355 72.5);
1356
1357 meDroppedFROMap_ = book2D(ibooker,
1358 "hDroppedFROMap",
1359 "Dropped full-readout-flagged readout units;"
1360 "iX - 40 / iEta / iX + 20;"
1361 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1362 "Event count",
1363 80,
1364 -39.5,
1365 40.5,
1366 72,
1367 .5,
1368 72.5);
1369
1370 meDroppedFROCnt_ = book1D(ibooker,
1371 "hDroppedFROCnt",
1372 "Number of dropped full-readout-flagged "
1373 "RU count;RU count;Event count",
1374 200,
1375 -.5,
1376 199.5);
1377
1378 meCompleteZSCnt_ = book1D(ibooker,
1379 "hCompleteZsCnt",
1380 "Number of zero-suppressed-flagged RU fully "
1381 "readout;"
1382 "RU count;Event count",
1383 200,
1384 -.5,
1385 199.5);
1386
1387 stringstream buf;
1388 buf << "Number of LI EB channels below the " << ebZsThr_ / 4.
1389 << " ADC count ZS threshold;"
1390 "Channel count;Event count",
1391 meEbZsErrCnt_ = book1D(ibooker, "hEbZsErrCnt", buf.str(), 200, -.5, 199.5);
1392
1393 buf.str("");
1394 buf << "Number of LI EE channels below the " << eeZsThr_ / 4.
1395 << " ADC count ZS theshold;"
1396 "Channel count;Event count",
1397 meEeZsErrCnt_ = book1D(ibooker, "hEeZsErrCnt", buf.str(), 200, -.5, 199.5);
1398
1399 meZsErrCnt_ = book1D(ibooker,
1400 "hZsErrCnt",
1401 "Number of LI channels below the ZS threshold;"
1402 "Channel count;Event count",
1403 200,
1404 -.5,
1405 199.5);
1406
1407 meEbZsErrType1Cnt_ = book1D(ibooker,
1408 "hEbZsErrType1Cnt",
1409 "Number of EB channels below the ZS "
1410 "threshold in a LI but fully readout RU;"
1411 "Channel count;Event count;",
1412 200,
1413 -.5,
1414 199.5);
1415
1416 meEeZsErrType1Cnt_ = book1D(ibooker,
1417 "hEeZsErrType1Cnt",
1418 "Number EE channels below the ZS threshold"
1419 " in a LI but fully readout RU;"
1420 "Channel count;Event count",
1421 200,
1422 -.5,
1423 199.5);
1424
1425 meZsErrType1Cnt_ = book1D(ibooker,
1426 "hZsErrType1Cnt",
1427 "Number of LI channels below the ZS threshold "
1428 "in a LI but fully readout RU;"
1429 "Channel count;Event count",
1430 200,
1431 -.5,
1432 199.5);
1433
1434 meDroppedFRORateMap_ = bookProfile2D(ibooker,
1435 "hDroppedFRORateMap",
1436 "Dropped full-readout-flagged readout units"
1437 "iX - 40 / iEta / iX + 20;"
1438 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1439 "Dropping rate",
1440 80,
1441 -39.5,
1442 40.5,
1443 72,
1444 .5,
1445 72.5);
1446
1447 meCompleteZSMap_ = book2D(ibooker,
1448 "hCompleteZSMap",
1449 "Complete zero-suppressed-flagged readout units;"
1450 "iX - 40 / iEta / iX + 20;"
1451 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1452 "Event count",
1453 80,
1454 -39.5,
1455 40.5,
1456 72,
1457 .5,
1458 72.5);
1459
1460 meCompleteZSRateMap_ = bookProfile2D(ibooker,
1461 "hCompleteZSRate",
1462 "Complete zero-suppressed-flagged readout units;"
1463 "iX - 40 / iEta / iX + 20;"
1464 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1465 "Completeness rate",
1466 80,
1467 -39.5,
1468 40.5,
1469 72,
1470 .5,
1471 72.5);
1472
1473
1474
1475 printAvailableHists();
1476
1477
1478 stringstream s;
1479 for (set<string>::iterator it = histList_.begin(); it != histList_.end(); ++it) {
1480 if (*it != string("all") && availableHistList_.find(*it) == availableHistList_.end()) {
1481 s << (s.str().empty() ? "" : ", ") << *it;
1482 }
1483 }
1484 if (!s.str().empty()) {
1485 LogWarning("Configuration") << "Parameter 'histList' contains some unknown histogram(s). "
1486 "Check spelling. Following name were not found: "
1487 << s.str();
1488 }
1489 }
1490
1491 void EcalSelectiveReadoutValidation::analyzeTP(edm::Event const& event, edm::EventSetup const& es) {
1492 int TTFlagCount[8];
1493 int LiTTFlagCount[nTtEta][nTtPhi];
1494 int MiTTFlagCount[nTtEta][nTtPhi];
1495 int HiTTFlagCount[nTtEta][nTtPhi];
1496 for (int iTTFlag(0); iTTFlag < 8; iTTFlag++) {
1497 TTFlagCount[iTTFlag] = 0;
1498 }
1499 for (int iTtEta(0); iTtEta < nTtEta; iTtEta++) {
1500 for (int iTtPhi(0); iTtPhi < nTtPhi; iTtPhi++) {
1501 LiTTFlagCount[iTtEta][iTtPhi] = 0;
1502 MiTTFlagCount[iTtEta][iTtPhi] = 0;
1503 HiTTFlagCount[iTtEta][iTtPhi] = 0;
1504 }
1505 }
1506 int tpEtCount[100];
1507 for (int iEt(0); iEt < 100; iEt++) {
1508 tpEtCount[iEt] = 0;
1509 }
1510
1511 const EcalTPGPhysicsConstMap& physMap = es.getData(physHandle).getMap();
1512
1513 const EcalTPGGroups::EcalTPGGroupsMap& lutGrpMap = es.getData(lutGrpHandle).getMap();
1514
1515 const EcalTPGLutIdMap::EcalTPGLutMap& lutMap = es.getData(lutMapHandle).getMap();
1516
1517 EcalTPGPhysicsConstMapIterator ebItr(physMap.find(DetId(DetId::Ecal, EcalBarrel).rawId()));
1518 double lsb10bitsEB(ebItr == physMap.end() ? 0. : ebItr->second.EtSat / 1024.);
1519 EcalTPGPhysicsConstMapIterator eeItr(physMap.find(DetId(DetId::Ecal, EcalEndcap).rawId()));
1520 double lsb10bitsEE(eeItr == physMap.end() ? 0. : eeItr->second.EtSat / 1024.);
1521
1522 for (EcalTrigPrimDigiCollection::const_iterator it = tps_->begin(); it != tps_->end(); ++it) {
1523 double tpEt;
1524 if (tpInGeV_) {
1525 EcalTrigTowerDetId const& towerId(it->id());
1526 unsigned int ADC = it->compressedEt();
1527
1528 double lsb10bits(0.);
1529 if (towerId.subDet() == EcalBarrel)
1530 lsb10bits = lsb10bitsEB;
1531 else if (towerId.subDet() == EcalEndcap)
1532 lsb10bits = lsb10bitsEE;
1533
1534 int tpg10bits = 0;
1535 EcalTPGGroups::EcalTPGGroupsMapItr itgrp = lutGrpMap.find(towerId.rawId());
1536 uint32_t lutGrp = 999;
1537 if (itgrp != lutGrpMap.end())
1538 lutGrp = itgrp->second;
1539
1540 EcalTPGLutIdMap::EcalTPGLutMapItr itLut = lutMap.find(lutGrp);
1541 if (itLut != lutMap.end()) {
1542 const unsigned int* lut = (itLut->second).getLut();
1543 for (unsigned int i = 0; i < 1024; i++)
1544 if (ADC == (0xff & lut[i])) {
1545 tpg10bits = i;
1546 break;
1547 }
1548 }
1549
1550 tpEt = lsb10bits * tpg10bits;
1551 } else {
1552 tpEt = it->compressedEt();
1553 }
1554 int iEta = it->id().ieta();
1555 int iEta0 = iTtEta2cIndex(iEta);
1556 int iPhi = it->id().iphi();
1557 int iPhi0 = iTtPhi2cIndex(iPhi);
1558 double etSum = ttEtSums[iEta0][iPhi0];
1559
1560 int iE = meTp_->getTProfile()->FindFixBin(tpEt);
1561 if ((iE >= 0) && (iE < 100)) {
1562 ++tpEtCount[iE];
1563 } else {
1564
1565
1566
1567 }
1568
1569 fill(meTpVsEtSum_, etSum, tpEt);
1570 ++TTFlagCount[it->ttFlag()];
1571 if ((it->ttFlag() & 0x3) == 0) {
1572 LiTTFlagCount[iEta0][iPhi0] += 1;
1573 } else if ((it->ttFlag() & 0x3) == 1) {
1574 MiTTFlagCount[iEta0][iPhi0] += 1;
1575 } else if ((it->ttFlag() & 0x3) == 3) {
1576 HiTTFlagCount[iEta0][iPhi0] += 1;
1577 }
1578 if ((it->ttFlag() & 0x4)) {
1579 fill(meForcedTtf_, iEta, iPhi);
1580 }
1581
1582 fill(meTtfVsTp_, tpEt, it->ttFlag());
1583 fill(meTtfVsEtSum_, etSum, it->ttFlag());
1584 fill(meTpMap_, iEta, iPhi, tpEt, 1.);
1585 }
1586
1587 for (int ittflag(0); ittflag < 8; ittflag++) {
1588 fill(meTtf_, ittflag, TTFlagCount[ittflag]);
1589 }
1590 for (int iTtEta(0); iTtEta < nTtEta; iTtEta++) {
1591 for (int iTtPhi(0); iTtPhi < nTtPhi; iTtPhi++) {
1592 fill(meLiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), LiTTFlagCount[iTtEta][iTtPhi]);
1593 fill(meMiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), MiTTFlagCount[iTtEta][iTtPhi]);
1594 fill(meHiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), HiTTFlagCount[iTtEta][iTtPhi]);
1595 }
1596 }
1597 if (tpInGeV_) {
1598 for (int iE(0); iE < 100; iE++) {
1599 fill(meTp_, iE, tpEtCount[iE]);
1600 }
1601 } else {
1602 for (int iE(0); iE < 40; iE++) {
1603 fill(meTp_, iE, tpEtCount[iE]);
1604 }
1605 }
1606 }
1607
1608 void EcalSelectiveReadoutValidation::analyzeDataVolume(const Event& e, const EventSetup& es) {
1609 anaDigiInit();
1610
1611
1612 for (int iDcc = minDccId_; iDcc <= maxDccId_; ++iDcc) {
1613 for (int iCh = 1; iCh < nDccRus_[iDcc - minDccId_]; ++iCh) {
1614 isRuComplete_[iDcc - minDccId_][iCh - 1] = (nPerRu_[iDcc - minDccId_][iCh - 1] == getCrystalCount(iDcc, iCh));
1615 }
1616 }
1617
1618
1619 for (unsigned int digis = 0; digis < ebDigis_->size(); ++digis) {
1620 EBDataFrame ebdf = (*ebDigis_)[digis];
1621 anaDigi(ebdf, *ebSrFlags_);
1622 }
1623
1624
1625 for (unsigned int digis = 0; digis < eeDigis_->size(); ++digis) {
1626 EEDataFrame eedf = (*eeDigis_)[digis];
1627 anaDigi(eedf, *eeSrFlags_);
1628 }
1629
1630
1631 for (unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0) {
1632 fill(meDccVol_, iDcc0 + 1, getDccEventSize(iDcc0, nPerDcc_[iDcc0]) / kByte_);
1633 fill(meDccLiVol_, iDcc0 + 1, getDccSrDependentPayload(iDcc0, nLiRuPerDcc_[iDcc0], nLiPerDcc_[iDcc0]) / kByte_);
1634 fill(meDccHiVol_, iDcc0 + 1, getDccSrDependentPayload(iDcc0, nHiRuPerDcc_[iDcc0], nHiPerDcc_[iDcc0]) / kByte_);
1635 const FEDRawDataCollection& raw = *fedRaw_;
1636 fill(meDccVolFromData_, iDcc0 + 1, ((double)raw.FEDData(601 + iDcc0).size()) / kByte_);
1637 }
1638
1639
1640 double a = nEbLI_ * getBytesPerCrystal() / kByte_;
1641 fill(meVolBLI_, a);
1642 double b = nEeLI_ * getBytesPerCrystal() / kByte_;
1643 fill(meVolELI_, b);
1644 fill(meVolLI_, a + b);
1645
1646
1647 a = nEbHI_ * getBytesPerCrystal() / kByte_;
1648 fill(meVolBHI_, a);
1649 b = nEeHI_ * getBytesPerCrystal() / kByte_;
1650 fill(meVolEHI_, b);
1651 fill(meVolHI_, a + b);
1652
1653
1654 a = getEbEventSize(nEb_) / kByte_;
1655 fill(meVolB_, a);
1656 b = getEeEventSize(nEe_) / kByte_;
1657 fill(meVolE_, b);
1658 fill(meVol_, a + b);
1659 }
1660
1661 template <class T, class U>
1662 void EcalSelectiveReadoutValidation::anaDigi(const T& frame, const U& srFlagColl) {
1663 const DetId& xtalId = frame.id();
1664 typedef typename U::key_type RuDetId;
1665 const RuDetId& ruId = readOutUnitOf(frame.id());
1666 typename U::const_iterator srf = srFlagColl.find(ruId);
1667
1668 bool highInterest = false;
1669 int flag = 0;
1670
1671 if (srf != srFlagColl.end()) {
1672 flag = srf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
1673
1674 highInterest = (flag == EcalSrFlag::SRF_FULL);
1675 }
1676
1677 bool barrel = (xtalId.subdetId() == EcalBarrel);
1678
1679 pair<int, int> ch = dccCh(xtalId);
1680
1681 if (barrel) {
1682 ++nEb_;
1683 if (highInterest) {
1684 ++nEbHI_;
1685 } else {
1686 ++nEbLI_;
1687 }
1688 int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(xtalId).ieta());
1689 int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(xtalId).iphi());
1690 if (!ebRuActive_[iEta0 / ebTtEdge][iPhi0 / ebTtEdge]) {
1691 ++nRuPerDcc_[ch.first - minDccId_];
1692 if (highInterest) {
1693 ++nHiRuPerDcc_[ch.first - minDccId_];
1694 } else {
1695 ++nLiRuPerDcc_[ch.first - minDccId_];
1696 }
1697
1698 ebRuActive_[iEta0 / ebTtEdge][iPhi0 / ebTtEdge] = true;
1699 }
1700 } else {
1701 ++nEe_;
1702 if (highInterest) {
1703 ++nEeHI_;
1704 } else {
1705 ++nEeLI_;
1706 }
1707 int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
1708 int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
1709 int iZ0 = static_cast<const EEDetId&>(frame.id()).zside() > 0 ? 1 : 0;
1710
1711 if (!eeRuActive_[iZ0][iX0 / scEdge][iY0 / scEdge]) {
1712 ++nRuPerDcc_[ch.first - minDccId_];
1713 if (highInterest) {
1714 ++nHiRuPerDcc_[ch.first - minDccId_];
1715 } else {
1716 ++nLiRuPerDcc_[ch.first - minDccId_];
1717 }
1718
1719 eeRuActive_[iZ0][iX0 / scEdge][iY0 / scEdge] = true;
1720 }
1721 }
1722
1723 if (ch.second < 1 || ch.second > 68) {
1724 throw cms::Exception("EcalSelectiveReadoutValidation")
1725 << "Error in DCC channel retrieval for crystal with detId " << xtalId.rawId()
1726 << "DCC channel out of allowed range [1..68]\n";
1727 }
1728 ++nPerDcc_[ch.first - minDccId_];
1729 ++nPerRu_[ch.first - minDccId_][ch.second - 1];
1730 if (highInterest) {
1731 ++nHiPerDcc_[ch.first - minDccId_];
1732 } else {
1733 ++nLiPerDcc_[ch.first - minDccId_];
1734 }
1735 }
1736
1737 void EcalSelectiveReadoutValidation::anaDigiInit() {
1738 nEb_ = 0;
1739 nEe_ = 0;
1740 nEeLI_ = 0;
1741 nEeHI_ = 0;
1742 nEbLI_ = 0;
1743 nEbHI_ = 0;
1744 bzero(nPerDcc_, sizeof(nPerDcc_));
1745 bzero(nLiPerDcc_, sizeof(nLiPerDcc_));
1746 bzero(nHiPerDcc_, sizeof(nHiPerDcc_));
1747 bzero(nRuPerDcc_, sizeof(nRuPerDcc_));
1748 bzero(ebRuActive_, sizeof(ebRuActive_));
1749 bzero(eeRuActive_, sizeof(eeRuActive_));
1750 bzero(nPerRu_, sizeof(nPerRu_));
1751 bzero(nLiRuPerDcc_, sizeof(nLiRuPerDcc_));
1752 bzero(nHiRuPerDcc_, sizeof(nHiRuPerDcc_));
1753 }
1754
1755 double EcalSelectiveReadoutValidation::frame2Energy(const EcalDataFrame& frame) const {
1756 static std::atomic<bool> firstCall{true};
1757 bool expected = true;
1758 if (firstCall.compare_exchange_strong(expected, false)) {
1759 stringstream buf;
1760 buf << "Weights:";
1761 for (unsigned i = 0; i < weights_.size(); ++i) {
1762 buf << "\t" << weights_[i];
1763 }
1764 edm::LogInfo("EcalSrValid") << buf.str() << "\n";
1765 firstCall = false;
1766 }
1767 double adc2GeV = 0.;
1768
1769 if (typeid(EBDataFrame) == typeid(frame)) {
1770 adc2GeV = .035;
1771 } else if (typeid(EEDataFrame) == typeid(frame)) {
1772 adc2GeV = 0.06;
1773 } else {
1774 assert(false);
1775 }
1776
1777 double acc = 0;
1778
1779 const int n = min(frame.size(), (int)weights_.size());
1780
1781 double gainInv[] = {12., 1., 6., 12.};
1782
1783 for (int i = 0; i < n; ++i) {
1784 acc += weights_[i] * frame[i].adc() * gainInv[frame[i].gainId()] * adc2GeV;
1785 }
1786 return acc;
1787 }
1788
1789 int EcalSelectiveReadoutValidation::getRuCount(int iDcc0) const { return nRuPerDcc_[iDcc0]; }
1790
1791 pair<int, int> EcalSelectiveReadoutValidation::dccCh(const DetId& detId) const {
1792 if (detId.det() != DetId::Ecal) {
1793 throw cms::Exception("InvalidParameter") << "Wrong type of DetId passed to the "
1794 "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
1795 "An ECAL DetId was expected.\n";
1796 }
1797
1798 DetId xtalId;
1799 switch (detId.subdetId()) {
1800 case EcalTriggerTower:
1801 {
1802 const EcalTrigTowerDetId tt = detId;
1803
1804
1805
1806 const int iTtPhi0 = iTtPhi2cIndex(tt.iphi());
1807 const int iTtEta0 = iTtEta2cIndex(tt.ieta());
1808 const int oneXtalPhi0 = iTtPhi0 * 5;
1809 const int oneXtalEta0 = (iTtEta0 - nOneEeTtEta) * 5;
1810
1811 xtalId = EBDetId(cIndex2iEta(oneXtalEta0), cIndex2iPhi(oneXtalPhi0));
1812 } break;
1813 case EcalEndcap:
1814 if (detId.rawId() & 0x8000) {
1815 return elecMap_->getDCCandSC(EcalScDetId(detId));
1816 } else {
1817 xtalId = detId;
1818 }
1819 break;
1820 case EcalBarrel:
1821 xtalId = detId;
1822 break;
1823 default:
1824 throw cms::Exception("InvalidParameter")
1825 << "Wrong type of DetId passed to the method "
1826 "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
1827 "A valid EcalTriggerTower, EcalBarrel or EcalEndcap DetId was expected. "
1828 "detid = "
1829 << xtalId.rawId() << ".\n";
1830 }
1831
1832 const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
1833
1834 pair<int, int> result;
1835 result.first = EcalElecId.dccId();
1836
1837 if (result.first < minDccId_ || result.second > maxDccId_) {
1838 throw cms::Exception("OutOfRange") << "Got an invalid DCC ID, DCCID = " << result.first << " for DetId 0x" << hex
1839 << detId.rawId() << " and 0x" << xtalId.rawId() << dec << "\n";
1840 }
1841
1842 result.second = EcalElecId.towerId();
1843
1844 if (result.second < 1 || result.second > 68) {
1845 throw cms::Exception("OutOfRange") << "Got an invalid DCC channel ID, DCC_CH = " << result.second << " for DetId 0x"
1846 << hex << detId.rawId() << " and 0x" << xtalId.rawId() << dec << "\n";
1847 }
1848
1849 return result;
1850 }
1851
1852 EcalTrigTowerDetId EcalSelectiveReadoutValidation::readOutUnitOf(const EBDetId& xtalId) const {
1853 return triggerTowerMap_->towerOf(xtalId);
1854 }
1855
1856 EcalScDetId EcalSelectiveReadoutValidation::readOutUnitOf(const EEDetId& xtalId) const {
1857 const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
1858 int iDCC = EcalElecId.dccId();
1859 int iDccChan = EcalElecId.towerId();
1860 const bool ignoreSingle = true;
1861 const vector<EcalScDetId> id = elecMap_->getEcalScDetId(iDCC, iDccChan, ignoreSingle);
1862 return !id.empty() ? id[0] : EcalScDetId();
1863 }
1864
1865 void EcalSelectiveReadoutValidation::setTtEtSums(const edm::EventSetup& es,
1866 const EBDigiCollection& ebDigis,
1867 const EEDigiCollection& eeDigis) {
1868
1869 const CaloSubdetectorGeometry* eeGeometry = nullptr;
1870 const CaloSubdetectorGeometry* ebGeometry = nullptr;
1871 if (eeGeometry == nullptr || ebGeometry == nullptr) {
1872 es.getData(geoToken);
1873 auto geoHandle = es.getHandle(geoToken);
1874 eeGeometry = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
1875 ebGeometry = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
1876 }
1877
1878
1879 for (int iEta0 = 0; iEta0 < nTtEta; ++iEta0) {
1880 for (int iPhi0 = 0; iPhi0 < nTtPhi; ++iPhi0) {
1881 ttEtSums[iEta0][iPhi0] = 0.;
1882 }
1883 }
1884
1885 for (EBDigiCollection::const_iterator it = ebDigis_->begin(); it != ebDigis_->end(); ++it) {
1886 const EBDataFrame& frame = *it;
1887 const EcalTrigTowerDetId& ttId = triggerTowerMap_->towerOf(frame.id());
1888
1889 const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
1890 const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
1891 double theta = ebGeometry->getGeometry(frame.id())->getPosition().theta();
1892 double e = frame2EnergyForTp(frame);
1893 if ((frame2EnergyForTp(frame, -1) < e) && (frame2EnergyForTp(frame, 1) < e)) {
1894 ttEtSums[iTtEta0][iTtPhi0] += e * sin(theta);
1895 }
1896 }
1897
1898 for (EEDigiCollection::const_iterator it = eeDigis.begin(); it != eeDigis.end(); ++it) {
1899 const EEDataFrame& frame = *it;
1900 const EcalTrigTowerDetId& ttId = triggerTowerMap_->towerOf(frame.id());
1901 const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
1902 const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
1903
1904 double theta = eeGeometry->getGeometry(frame.id())->getPosition().theta();
1905 double e = frame2EnergyForTp(frame);
1906 if ((frame2EnergyForTp(frame, -1) < e) && (frame2EnergyForTp(frame, 1) < e)) {
1907 ttEtSums[iTtEta0][iTtPhi0] += e * sin(theta);
1908 }
1909 }
1910
1911
1912 int innerTTEtas[] = {0, 1, 54, 55};
1913 for (unsigned iRing = 0; iRing < sizeof(innerTTEtas) / sizeof(innerTTEtas[0]); ++iRing) {
1914 int iTtEta0 = innerTTEtas[iRing];
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930 for (unsigned iTtPhi0 = 0; iTtPhi0 < nTtPhi - 1; iTtPhi0 += 2) {
1931 double et = .5 * (ttEtSums[iTtEta0][iTtPhi0] + ttEtSums[iTtEta0][iTtPhi0 + 1]);
1932
1933
1934
1935 ttEtSums[iTtEta0][iTtPhi0] = et;
1936 ttEtSums[iTtEta0][iTtPhi0 + 1] = et;
1937 }
1938 }
1939 }
1940
1941 template <class T>
1942 double EcalSelectiveReadoutValidation::frame2EnergyForTp(const T& frame, int offset) const {
1943
1944
1945 double weights[] = {0., -1 / 3., -1 / 3., -1 / 3., 0., 1.};
1946
1947 double adc2GeV = 0.;
1948 if (typeid(frame) == typeid(EBDataFrame)) {
1949 adc2GeV = 0.035;
1950 } else if (typeid(frame) == typeid(EEDataFrame)) {
1951 adc2GeV = 0.060;
1952 } else {
1953
1954 throw cms::Exception("Severe Error") << __FILE__ << ":" << __LINE__ << ": "
1955 << "this is a bug. Please report it.\n";
1956 }
1957
1958 double acc = 0;
1959
1960 const int n = min<int>(frame.size(), sizeof(weights) / sizeof(weights[0]));
1961
1962 double gainInv[] = {12., 1., 6., 12};
1963
1964 for (int i = offset; i < n; ++i) {
1965 int iframe = i + offset;
1966 if (iframe >= 0 && iframe < frame.size()) {
1967 acc += weights[i] * frame[iframe].adc() * gainInv[frame[iframe].gainId()] * adc2GeV;
1968 }
1969 }
1970
1971 return acc;
1972 }
1973
1974 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::bookFloat(DQMStore::IBooker& ibook,
1975 const std::string& name) {
1976 if (!registerHist(name, ""))
1977 return nullptr;
1978 MonitorElement* result = ibook.bookFloat(name);
1979 if (result == nullptr) {
1980 throw cms::Exception("DQM") << "Failed to book integer DQM monitor element" << name;
1981 }
1982 return result;
1983 }
1984
1985 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::book1D(
1986 DQMStore::IBooker& ibook, const std::string& name, const std::string& title, int nbins, double xmin, double xmax) {
1987 if (!registerHist(name, title))
1988 return nullptr;
1989 MonitorElement* result = ibook.book1D(name, title, nbins, xmin, xmax);
1990 if (result == nullptr) {
1991 throw cms::Exception("Histo") << "Failed to book histogram " << name;
1992 }
1993 return result;
1994 }
1995
1996 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::book2D(DQMStore::IBooker& ibook,
1997 const std::string& name,
1998 const std::string& title,
1999 int nxbins,
2000 double xmin,
2001 double xmax,
2002 int nybins,
2003 double ymin,
2004 double ymax) {
2005 if (!registerHist(name, title))
2006 return nullptr;
2007 MonitorElement* result = ibook.book2D(name, title, nxbins, xmin, xmax, nybins, ymin, ymax);
2008 if (result == nullptr) {
2009 throw cms::Exception("Histo") << "Failed to book histogram " << name;
2010 }
2011 return result;
2012 }
2013
2014 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::bookProfile(
2015 DQMStore::IBooker& ibook, const std::string& name, const std::string& title, int nbins, double xmin, double xmax) {
2016 if (!registerHist(name, title))
2017 return nullptr;
2018 MonitorElement* result = ibook.bookProfile(name, title, nbins, xmin, xmax, 0, 0, 0);
2019 if (result == nullptr) {
2020 throw cms::Exception("Histo") << "Failed to book histogram " << name;
2021 }
2022 return result;
2023 }
2024
2025 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::bookProfile2D(DQMStore::IBooker& ibook,
2026 const std::string& name,
2027 const std::string& title,
2028 int nbinx,
2029 double xmin,
2030 double xmax,
2031 int nbiny,
2032 double ymin,
2033 double ymax,
2034 const char* option) {
2035 if (!registerHist(name, title))
2036 return nullptr;
2037 MonitorElement* result = ibook.bookProfile2D(name, title, nbinx, xmin, xmax, nbiny, ymin, ymax, 0, 0, 0, option);
2038 if (result == nullptr) {
2039 throw cms::Exception("Histo") << "Failed to book histogram " << name;
2040 }
2041 return result;
2042 }
2043
2044 bool EcalSelectiveReadoutValidation::registerHist(const std::string& name, const std::string& title) {
2045 availableHistList_.insert(pair<string, string>(name, title));
2046 return allHists_ || histList_.find(name) != histList_.end();
2047 }
2048
2049 void EcalSelectiveReadoutValidation::readAllCollections(const edm::Event& event) {
2050 ebRecHits_.read(event);
2051 eeRecHits_.read(event);
2052 ebDigis_.read(event);
2053 eeDigis_.read(event);
2054 ebNoZsDigis_.read(event);
2055 eeNoZsDigis_.read(event);
2056 ebSrFlags_.read(event);
2057 eeSrFlags_.read(event);
2058 ebComputedSrFlags_.read(event);
2059 eeComputedSrFlags_.read(event);
2060 ebSimHits_.read(event);
2061 eeSimHits_.read(event);
2062 tps_.read(event);
2063 fedRaw_.read(event);
2064 }
2065
2066 void EcalSelectiveReadoutValidation::printAvailableHists() {
2067 LogInfo log("HistoList");
2068 log << "Avalailable histograms (DQM monitor elements): \n";
2069 for (map<string, string>::iterator it = availableHistList_.begin(); it != availableHistList_.end(); ++it) {
2070 log << it->first << ": " << it->second << "\n";
2071 }
2072 log << "\nTo include an histogram add its name in the vstring parameter "
2073 "'histograms' of the EcalSelectiveReadoutValidation module\n";
2074 }
2075
2076 double EcalSelectiveReadoutValidation::getEbEventSize(double nReadXtals) const {
2077 double ruHeaderPayload = 0.;
2078 const int firstEbDcc0 = nEeDccs / 2;
2079 for (int iDcc0 = firstEbDcc0; iDcc0 < firstEbDcc0 + nEbDccs; ++iDcc0) {
2080 ruHeaderPayload += getRuCount(iDcc0) * 8.;
2081 }
2082
2083 return getDccOverhead(EB) * nEbDccs + nReadXtals * getBytesPerCrystal() + ruHeaderPayload;
2084 }
2085
2086 double EcalSelectiveReadoutValidation::getEeEventSize(double nReadXtals) const {
2087 double ruHeaderPayload = 0.;
2088 const unsigned firstEbDcc0 = nEeDccs / 2;
2089 for (unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0) {
2090
2091 if (iDcc0 == firstEbDcc0)
2092 iDcc0 += nEbDccs;
2093 ruHeaderPayload += getRuCount(iDcc0) * 8.;
2094 }
2095 return getDccOverhead(EE) * nEeDccs + nReadXtals * getBytesPerCrystal() + ruHeaderPayload;
2096 }
2097
2098
2099
2100 int EcalSelectiveReadoutValidation::dccZsFIR(const EcalDataFrame& frame,
2101 const std::vector<int>& firWeights,
2102 int firstFIRSample,
2103 bool* saturated) {
2104 const int nFIRTaps = 6;
2105
2106 const vector<int>& w = firWeights;
2107
2108
2109 int acc = 0;
2110 bool gain12saturated = false;
2111 const int gain12 = 0x01;
2112 const int lastFIRSample = firstFIRSample + nFIRTaps - 1;
2113
2114 int iWeight = 0;
2115 for (int iSample = firstFIRSample - 1; iSample < lastFIRSample; ++iSample, ++iWeight) {
2116 if (iSample >= 0 && iSample < frame.size()) {
2117 EcalMGPASample sample(frame[iSample]);
2118 if (sample.gainId() != gain12)
2119 gain12saturated = true;
2120 LogTrace("DccFir") << (iSample >= firstFIRSample ? "+" : "") << sample.adc() << "*(" << w[iWeight] << ")";
2121 acc += sample.adc() * w[iWeight];
2122 } else {
2123 edm::LogWarning("DccFir") << __FILE__ << ":" << __LINE__
2124 << ": Not enough samples in data frame or 'ecalDccZs1stSample' module "
2125 "parameter is not valid...";
2126 }
2127 }
2128 LogTrace("DccFir") << "\n";
2129
2130
2131
2132 acc = (acc >= 0) ? (acc >> 8) : -(-acc >> 8);
2133
2134
2135
2136
2137 LogTrace("DccFir") << "acc: " << acc << "\n"
2138 << "saturated: " << (gain12saturated ? "yes" : "no") << "\n";
2139
2140 if (saturated) {
2141 *saturated = gain12saturated;
2142 }
2143
2144 return gain12saturated ? numeric_limits<int>::max() : acc;
2145 }
2146
2147 std::vector<int> EcalSelectiveReadoutValidation::getFIRWeights(const std::vector<double>& normalizedWeights) {
2148 const int nFIRTaps = 6;
2149 vector<int> firWeights(nFIRTaps, 0);
2150 const static int maxWeight = 0xEFF;
2151 for (unsigned i = 0; i < min((size_t)nFIRTaps, normalizedWeights.size()); ++i) {
2152 firWeights[i] = lround(normalizedWeights[i] * (1 << 10));
2153 if (abs(firWeights[i]) > maxWeight) {
2154 firWeights[i] = firWeights[i] < 0 ? -maxWeight : maxWeight;
2155 }
2156 }
2157 return firWeights;
2158 }
2159
2160 void EcalSelectiveReadoutValidation::configFirWeights(const vector<double>& weightsForZsFIR) {
2161 bool notNormalized = false;
2162 bool notInt = false;
2163 for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2164 if (weightsForZsFIR[i] > 1.)
2165 notNormalized = true;
2166 if ((int)weightsForZsFIR[i] != weightsForZsFIR[i])
2167 notInt = true;
2168 }
2169 if (notInt && notNormalized) {
2170 throw cms::Exception("InvalidParameter") << "weigtsForZsFIR paramater values are not valid: they "
2171 << "must either be integer and uses the hardware representation "
2172 << "of the weights or less or equal than 1 and used the normalized "
2173 << "representation.";
2174 }
2175 LogInfo log("DccFir");
2176 if (notNormalized) {
2177 firWeights_ = vector<int>(weightsForZsFIR.size());
2178 for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2179 firWeights_[i] = (int)weightsForZsFIR[i];
2180 }
2181 } else {
2182 firWeights_ = getFIRWeights(weightsForZsFIR);
2183 }
2184
2185 log << "Input weights for FIR: ";
2186 for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2187 log << weightsForZsFIR[i] << "\t";
2188 }
2189
2190 double s2 = 0.;
2191 log << "\nActual FIR weights: ";
2192 for (unsigned i = 0; i < firWeights_.size(); ++i) {
2193 log << firWeights_[i] << "\t";
2194 s2 += firWeights_[i] * firWeights_[i];
2195 }
2196
2197 s2 = sqrt(s2);
2198 log << "\nNormalized FIR weights after hw representation rounding: ";
2199 for (unsigned i = 0; i < firWeights_.size(); ++i) {
2200 log << firWeights_[i] / (double)(1 << 10) << "\t";
2201 }
2202
2203 log << "\nFirst FIR sample: " << firstFIRSample_;
2204 }
2205
2206 void EcalSelectiveReadoutValidation::initAsciiFile() {
2207 if (logSrpAlgoErrors_) {
2208 srpAlgoErrorLog_.open(srpAlgoErrorLogFileName_.c_str(), ios::out | ios::trunc);
2209 if (!srpAlgoErrorLog_.good()) {
2210 throw cms::Exception("Output") << "Failed to open the log file '" << srpAlgoErrorLogFileName_
2211 << "' for SRP algorithm result check.\n";
2212 }
2213 }
2214
2215 if (logSrApplicationErrors_) {
2216 srApplicationErrorLog_.open(srApplicationErrorLogFileName_.c_str(), ios::out | ios::trunc);
2217 if (!srApplicationErrorLog_.good()) {
2218 throw cms::Exception("Output") << "Failed to open the log file '" << srApplicationErrorLogFileName_
2219 << "' for Selective Readout decision application check.\n";
2220 }
2221 }
2222 }
2223
2224
2225
2226
2227 template <class T>
2228 void EcalSelectiveReadoutValidation::compareSrfColl(const edm::Event& event, T& srfFromData, T& computedSrf) {
2229 typedef typename T::const_iterator SrFlagCollectionConstIt;
2230 typedef typename T::key_type MyRuDetIdType;
2231 SrFlagCollectionConstIt itSrfFromData = srfFromData.begin();
2232 SrFlagCollectionConstIt itComputedSr = computedSrf.begin();
2233
2234 while (itSrfFromData != srfFromData.end() || itComputedSr != computedSrf.end()) {
2235 MyRuDetIdType inconsistentRu = 0;
2236 bool inconsistent = false;
2237 if (itComputedSr == computedSrf.end() ||
2238 (itSrfFromData != srfFromData.end() && itSrfFromData->id() < itComputedSr->id())) {
2239
2240 pair<int, int> ch = dccCh(itSrfFromData->id());
2241 srpAlgoErrorLog_ << event.id() << ": " << itSrfFromData->id() << ", DCC " << ch.first << " ch " << ch.second
2242 << " found in data (SRF:" << itSrfFromData->flagName()
2243 << ") but not in the set of SRFs computed from the data TTF.\n";
2244 inconsistentRu = itSrfFromData->id();
2245 inconsistent = true;
2246 ++itSrfFromData;
2247 } else if (itSrfFromData == srfFromData.end() ||
2248 (itComputedSr != computedSrf.end() && itComputedSr->id() < itSrfFromData->id())) {
2249
2250 pair<int, int> ch = dccCh(itComputedSr->id());
2251 if (logErrForDccs_[ch.first - minDccId_]) {
2252 srpAlgoErrorLog_ << event.id() << ": " << itComputedSr->id() << ", DCC " << ch.first << " ch " << ch.second
2253 << " not found in data. Computed SRF: " << itComputedSr->flagName() << ".\n";
2254 inconsistentRu = itComputedSr->id();
2255 inconsistent = true;
2256 }
2257 ++itComputedSr;
2258 } else {
2259
2260 if (itComputedSr->value() != itSrfFromData->value()) {
2261 pair<int, int> ch = dccCh(itSrfFromData->id());
2262 srpAlgoErrorLog_ << event.id() << ", " << itSrfFromData->id() << ", DCC " << ch.first << " ch " << ch.second
2263 << ", SRF inconsistency: "
2264 << "from data: " << itSrfFromData->flagName()
2265 << ", computed from TTF: " << itComputedSr->flagName() << "\n";
2266 inconsistentRu = itComputedSr->id();
2267 inconsistent = true;
2268 }
2269 if (itComputedSr != computedSrf.end())
2270 ++itComputedSr;
2271 if (itSrfFromData != srfFromData.end())
2272 ++itSrfFromData;
2273 }
2274
2275 if (inconsistent)
2276 fill(meSRFlagsConsistency_, ruGraphX(inconsistentRu), ruGraphY(inconsistentRu));
2277 }
2278 }
2279
2280 int EcalSelectiveReadoutValidation::dccId(const EcalScDetId& detId) const { return elecMap_->getDCCandSC(detId).first; }
2281
2282 int EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId& detId) const {
2283 if (detId.ietaAbs() > 17) {
2284 throw cms::Exception("InvalidArgument")
2285 << "Argument of EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId&) "
2286 << "must be a barrel trigger tower Id\n";
2287 }
2288 return dccCh(detId).first;
2289 }
2290
2291 void EcalSelectiveReadoutValidation::selectFedsForLog() {
2292 logErrForDccs_ = vector<bool>(nDccs_, false);
2293
2294 for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
2295 int iDcc = dccId(it->id()) - minDccId_;
2296
2297 logErrForDccs_.at(iDcc) = true;
2298 }
2299
2300 for (EESrFlagCollection::const_iterator it = eeSrFlags_->begin(); it != eeSrFlags_->end(); ++it) {
2301 int iDcc = dccId(it->id()) - minDccId_;
2302
2303 logErrForDccs_.at(iDcc) = true;
2304 }
2305
2306 stringstream buf;
2307 buf << "List of DCCs found in the first processed event: ";
2308 bool first = true;
2309 for (unsigned iDcc = 0; iDcc < nDccs_; ++iDcc) {
2310 if (logErrForDccs_[iDcc]) {
2311 buf << (first ? "" : ", ") << (iDcc + minDccId_);
2312 first = false;
2313 }
2314 }
2315 buf << "\nOnly DCCs from this list will be considered for error logging\n";
2316 srpAlgoErrorLog_ << buf.str();
2317 srApplicationErrorLog_ << buf.str();
2318 LogInfo("EcalSrValid") << buf.str();
2319 }
2320
2321 template <class T>
2322 void EcalSelectiveReadoutValidation::checkSrApplication(const edm::Event& event, T& srfs) {
2323 typedef typename T::const_iterator SrFlagCollectionConstIt;
2324 typedef typename T::key_type MyRuDetIdType;
2325
2326 for (SrFlagCollectionConstIt itSrf = srfs.begin(); itSrf != srfs.end(); ++itSrf) {
2327 int flag = itSrf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
2328 pair<int, int> ru = dccCh(itSrf->id());
2329
2330 if (flag == EcalSrFlag::SRF_FULL) {
2331 if (nPerRu_[ru.first - minDccId_][ru.second - 1] == getCrystalCount(ru.first, ru.second)) {
2332 fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2333 fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2334 } else if (nPerRu_[ru.first - minDccId_][ru.second - 1] == 0) {
2335 fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2336 fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2337 fill(meDroppedFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2338 ++nDroppedFRO_;
2339 srApplicationErrorLog_ << event.id() << ": Flag of RU " << itSrf->id() << " (DCC " << ru.first << " ch "
2340 << ru.second << ") is 'Full readout' "
2341 << "while none of its channel was read out\n";
2342 } else {
2343 fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2344 fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2345 fill(meIncompleteFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2346 ++nIncompleteFRO_;
2347 srApplicationErrorLog_ << event.id() << ": Flag of RU" << itSrf->id() << " (DCC " << ru.first << " ch "
2348 << ru.second << ") is 'Full readout' "
2349 << "while only " << nPerRu_[ru.first - minDccId_][ru.second - 1] << " / "
2350 << getCrystalCount(ru.first, ru.second) << " channels were read out.\n";
2351 }
2352 }
2353
2354 if (flag == EcalSrFlag::SRF_ZS1 || flag == EcalSrFlag::SRF_ZS2) {
2355 if (nPerRu_[ru.first - minDccId_][ru.second - 1] == getCrystalCount(ru.first, ru.second)) {
2356
2357
2358 fill(meCompleteZSMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()));
2359 fill(meCompleteZSRateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2360
2361 ++nCompleteZS_;
2362 } else {
2363 fill(meCompleteZSRateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2364 }
2365 }
2366 }
2367 }
2368
2369 int EcalSelectiveReadoutValidation::getCrystalCount(int iDcc, int iDccCh) {
2370 if (iDcc < minDccId_ || iDcc > maxDccId_) {
2371 return 0;
2372 } else if (10 <= iDcc && iDcc <= 45) {
2373 return 25;
2374 } else {
2375 int iDccPhi;
2376 if (iDcc < 10)
2377 iDccPhi = iDcc;
2378 else
2379 iDccPhi = iDcc - 45;
2380 switch (iDccPhi * 100 + iDccCh) {
2381 case 110:
2382 case 232:
2383 case 312:
2384 case 412:
2385 case 532:
2386 case 610:
2387 case 830:
2388 case 806:
2389
2390 return 20;
2391 case 134:
2392 case 634:
2393 case 827:
2394 case 803:
2395 return 10;
2396 case 330:
2397 case 430:
2398 return 20;
2399 case 203:
2400 case 503:
2401 case 721:
2402 case 921:
2403 return 21;
2404 default:
2405 return 25;
2406 }
2407 }
2408 }
