Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:06:09

0001 /*
0002  * \file EcalSelectiveReadoutValidation.cc
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     //EE- DCCs:
0036     //   1  2   3   4   5   6   7   8   9
0037     /**/ 34,
0038     32,
0039     33,
0040     33,
0041     32,
0042     34,
0043     33,
0044     34,
0045     33,
0046     //EB- DCCs:
0047     //  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27
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     //EB+ DCCs:
0067     //  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45
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     //EE+ DCCs:
0087     //  46  47  48  49  50  51  52  53  54
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 /*collNotFoundWarn_*/),
0109       eeNoZsDigis_(ps.getParameter<edm::InputTag>("EeUnsuppressedDigiCollection"), false, false /*collNotFoundWarn_*/),
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 /*collNotFoundWarn_*/),
0114       eeComputedSrFlags_(
0115           ps.getParameter<edm::InputTag>("EeSrFlagFromTTCollection"), false, false /*collNotFoundWarn_*/),
0116       ebSimHits_(ps.getParameter<edm::InputTag>("EbSimHitCollection"), false, false /*collNotFoundWarn_*/),
0117       eeSimHits_(ps.getParameter<edm::InputTag>("EeSimHitCollection"), false, false /*collNotFoundWarn_*/),
0118       tps_(ps.getParameter<edm::InputTag>("TrigPrimCollection"), false, collNotFoundWarn_),
0119       ebRecHits_(ps.getParameter<edm::InputTag>("EbRecHitCollection"), false, false /*collNotFoundWarn_*/),
0120       eeRecHits_(ps.getParameter<edm::InputTag>("EeRecHitCollection"), false, false /*collNotFoundWarn_*/),
0121       fedRaw_(ps.getParameter<edm::InputTag>("FEDRawCollection"), false, false /*collNotFoundWarn_*/),
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   //File to log SRP algorithem inconsistency
0161   srpAlgoErrorLogFileName_ = ps.getUntrackedParameter<string>("srpAlgoErrorLogFile", "");
0162   logSrpAlgoErrors_ = (!srpAlgoErrorLogFileName_.empty());
0163 
0164   //File to log SRP decision application inconsistency
0165   srApplicationErrorLogFileName_ = ps.getUntrackedParameter<string>("srApplicationErrorLogFile", "");
0166   logSrApplicationErrors_ = (!srApplicationErrorLogFileName_.empty());
0167 
0168   //FIR ZS weights
0169   configFirWeights(ps.getParameter<vector<double>>("dccWeights"));
0170 
0171   // DQM ROOT output
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   // verbosity switch
0181   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0182 
0183   // get hold of back-end interface
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   //retrieves event products:
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();  //note: must be called after readAllCollection
0239   }
0240 
0241   //computes Et sum trigger tower crystals:
0242   setTtEtSums(es, *ebNoZsDigis_, *eeNoZsDigis_);
0243 
0244   //Data Volume
0245   analyzeDataVolume(event, es);
0246 
0247   //EB digis
0248   //must be called after analyzeDataVolume because it uses
0249   //isRuComplete_ array that this method fills
0250   analyzeEB(event, es);
0251 
0252   //EE digis
0253   //must be called after analyzeDataVolume because it uses
0254   //isRuComplete_ array that this method fills
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   //TP
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   /** Energy deposited in ECAL endcap crystals. Endcap index is 0 for EE- and
0295    * 1 for EE+. X and Y index starts at x and y minimum in std CMS coordinate
0296    * system.*/
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;  //must be set to zero.
0306         (*eeEnergies)[iZ0][iX0][iY0].simHit = 0;
0307         (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
0308       }
0309     }
0310   }
0311 
0312   // gets the endcap geometry:
0313   auto geoHandle = es.getHandle(geoToken);
0314   const CaloSubdetectorGeometry* geometry_p = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0315   //CaloSubdetectorGeometry const& geometry = *geometry_p;
0316 
0317   //EE unsupressed digis:
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     //    cout << "EE no ZS energy computation..." ;
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   //EE rec hits:
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       //    cout << "EE no ZS energy computation..." ;
0365       (*eeEnergies)[iZ0][iX0][iY0].recE = hit.energy();
0366     }
0367   }
0368 
0369   //EE sim hits:
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   //EE suppressed digis
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   }  //next ZS digi.
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;  //not a crystal or ZS
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) {  //noise only crystal channel
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   //Filling RU histo
0517   for (EESrFlagCollection::const_iterator it = eeSrFlags_->begin(); it != eeSrFlags_->end(); ++it) {
0518     const EESrFlag& srf = *it;
0519     // srf.id() is EcalScDetId; 1 <= ix <= 20 1 <= iy <= 20
0520     int iX = srf.id().ix();
0521     int iY = srf.id().iy();
0522     int zside = srf.id().zside();  //-1 for EE-, +1 for EE+
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 }  //end of analyzeEE
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   /** Energy deposited in ECAL barrel crystals. Eta index starts from 0 at
0572    * eta minimum and phi index starts at phi=0+ in CMS std coordinate system.
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;  //must be zero.
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   // get the barrel geometry:
0589   auto geoHandle = es.getHandle(geoToken);
0590   const CaloSubdetectorGeometry* geometry_p = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0591   //CaloSubdetectorGeometry const& geometry = *geometry_p;
0592 
0593   //EB unsuppressed digis:
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   }  //next non-zs digi
0626 
0627   //EB sim hits
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   for (EBDigiCollection::const_iterator it = ebDigis_->begin(); it != ebDigis_->end(); ++it) {
0650     const EBDataFrame& frame = *it;
0651     int iEta = static_cast<const EBDetId&>(frame.id()).ieta();
0652     int iPhi = static_cast<const EBDetId&>(frame.id()).iphi();
0653     int iEta0 = iEta2cIndex(iEta);
0654     int iPhi0 = iPhi2cIndex(iPhi);
0655     if (iEta0 < 0 || iEta0 >= nEbEta) {
0656       throw(cms::Exception("EcalSelectiveReadoutValidation") << "iEta0 (= " << iEta0 << ") is out of range ("
0657                                                              << "[0," << nEbEta - 1 << "]");
0658     }
0659     if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
0660       throw(cms::Exception("EcalSelectiveReadoutValidation") << "iPhi0 (= " << iPhi0 << ") is out of range ("
0661                                                              << "[0," << nEbPhi - 1 << "]");
0662     }
0663     assert(iEta0 >= 0 && iEta0 < nEbEta);
0664     assert(iPhi0 >= 0 && iPhi0 < nEbPhi);
0665     if (!crystalShot[iEta0][iPhi0]) {
0666       crystalShot[iEta0][iPhi0] = true;
0667       EBxtalCoor[iEta0][iPhi0] = make_pair(xtalGraphX(frame.id()), xtalGraphY(frame.id()));
0668     } else {
0669       cout << "Error: several digi for same crystal!";
0670       abort();
0671     }
0672     if (localReco_) {
0673       (*ebEnergies)[iEta0][iPhi0].recE = frame2Energy(frame);
0674     }
0675 
0676     (*ebEnergies)[iEta0][iPhi0].gain12 = true;
0677     for (int i = 0; i < frame.size(); ++i) {
0678       const int gain12Code = 0x1;
0679       if (frame[i].gainId() != gain12Code) {
0680         (*ebEnergies)[iEta0][iPhi0].gain12 = false;
0681       }
0682     }
0683 
0684     EBSrFlagCollection::const_iterator srf = ebSrFlags_->find(readOutUnitOf(frame.id()));
0685 
0686     bool highInterest = false;
0687 
0688     // if(srf == ebSrFlags_->end()){
0689     //       throw cms::Exception("EcalSelectiveReadoutValidation")
0690     //  << __FILE__ << ":" << __LINE__ << ": SR flag not found";
0691     //}
0692 
0693     if (srf != ebSrFlags_->end()) {
0694       highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK) == EcalSrFlag::SRF_FULL);
0695     }
0696 
0697     if (highInterest) {
0698       fill(meEbHiZsFir_, dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr));
0699     } else {
0700       int v = dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr);
0701       fill(meEbLiZsFir_, v);
0702       if (v < ebZsThr_) {
0703         eventError = true;
0704         ++nEbZsErrors_;
0705         pair<int, int> ru = dccCh(frame.id());
0706         if (isRuComplete_[ru.first - 1][ru.second - 1])
0707           ++nEbZsErrorsType1_;
0708         if (nEbZsErrors_ < 3) {
0709           srApplicationErrorLog_ << event.id() << ", "
0710                                  << "RU " << frame.id() << ", "
0711                                  << "DCC " << ru.first << " Ch : " << ru.second << ": "
0712                                  << "LI channel under ZS threshold.\n";
0713         }
0714         if (nEbZsErrors_ == 3) {
0715           srApplicationErrorLog_ << event.id() << ": "
0716                                  << "more ZS errors for this event...\n";
0717         }
0718       }
0719     }
0720   }  //next EB digi
0721 
0722   for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
0723     for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0)
0724       fill(meChOcc_,
0725            EBxtalCoor[iEta0][iPhi0].first,
0726            EBxtalCoor[iEta0][iPhi0].second,
0727            crystalShot[iEta0][iPhi0] ? 1. : 0.);
0728   }
0729 
0730   if (!localReco_) {
0731     for (RecHitCollection::const_iterator it = ebRecHits_->begin(); it != ebRecHits_->end(); ++it) {
0732       const RecHit& hit = *it;
0733       int iEta = static_cast<const EBDetId&>(hit.id()).ieta();
0734       int iPhi = static_cast<const EBDetId&>(hit.id()).iphi();
0735       int iEta0 = iEta2cIndex(iEta);
0736       int iPhi0 = iPhi2cIndex(iPhi);
0737       if (iEta0 < 0 || iEta0 >= nEbEta) {
0738         LogError("EcalSrValid") << "iEta0 (= " << iEta0 << ") is out of range ("
0739                                 << "[0," << nEbEta - 1 << "]\n";
0740       }
0741       if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
0742         LogError("EcalSrValid") << "iPhi0 (= " << iPhi0 << ") is out of range ("
0743                                 << "[0," << nEbPhi - 1 << "]\n";
0744       }
0745       (*ebEnergies)[iEta0][iPhi0].recE = hit.energy();
0746     }
0747   }
0748 
0749   for (unsigned int i = 0; i < xtalEtaPhi.size(); ++i) {
0750     int iEta0 = xtalEtaPhi[i].first;
0751     int iPhi0 = xtalEtaPhi[i].second;
0752     energiesEb_t& energies = (*ebEnergies)[iEta0][iPhi0];
0753 
0754     double recE = energies.recE;
0755     if (recE != -numeric_limits<double>::max()) {  //not zero suppressed
0756       fill(meEbRecE_, (*ebEnergies)[iEta0][iPhi0].recE);
0757       fill(meEbEMean_, ievt_ + 1, recE);
0758     }  //not zero suppressed
0759 
0760     if (withEbSimHit_) {
0761       if (!energies.simHit) {  //noise only crystal channel
0762         fill(meEbNoise_, energies.noZsRecE);
0763       } else {
0764         fill(meEbSimE_, energies.simE);
0765         fill(meEbRecEHitXtal_, energies.recE);
0766       }
0767       fill(meEbRecVsSimE_, energies.simE, energies.recE);
0768       fill(meEbNoZsRecVsSimE_, energies.simE, energies.noZsRecE);
0769     }
0770   }
0771 
0772   int EBZs1RuCount[2][17][72];
0773   int EBFullRuCount[2][17][72];
0774   int EBForcedRuCount[2][17][72];
0775   std::pair<int, int> EBtowerCoor[2][17][72];
0776   for (int iZ(0); iZ < 2; iZ++) {
0777     for (int iEta(0); iEta < 17; iEta++) {
0778       for (int iPhi(0); iPhi < 72; iPhi++) {
0779         EBZs1RuCount[iZ][iEta][iPhi] = 0;
0780         EBFullRuCount[iZ][iEta][iPhi] = 0;
0781         EBForcedRuCount[iZ][iEta][iPhi] = 0;
0782       }
0783     }
0784   }
0785 
0786   //SRF
0787   nEbFROCnt_ = 0;
0788   char ebSrfMark[2][17][72];
0789   bzero(ebSrfMark, sizeof(ebSrfMark));
0790   //      int idbg = 0;
0791   for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
0792     const EBSrFlag& srf = *it;
0793     int iEtaAbs = srf.id().ietaAbs();
0794     int iPhi = srf.id().iphi();
0795     int iZ = srf.id().zside();
0796 
0797     //  cout << "--> " << ++idbg << iEtaAbs << " " << iPhi << " "  << iZ
0798     //       << " " << srf.id() << "\n";
0799 
0800     if (iEtaAbs < 1 || iEtaAbs > 17 || iPhi < 1 || iPhi > 72)
0801       throw cms::Exception("EcalSelectiveReadoutValidation")
0802           << "Found a barrel SRF with an invalid det ID: " << srf.id() << ".\n";
0803     ++ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1];
0804     if (ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] > 1)
0805       throw cms::Exception("EcalSelectiveReadoutValidation") << "Duplicate SRF for RU " << srf.id() << ".\n";
0806 
0807     EBtowerCoor[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] = std::pair<int, int>(srf.id().ieta(), srf.id().iphi());
0808 
0809     int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
0810     if (flag == EcalSrFlag::SRF_ZS1) {
0811       EBZs1RuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
0812     }
0813     if (flag == EcalSrFlag::SRF_FULL) {
0814       EBFullRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
0815       ++nEbFROCnt_;
0816     }
0817     if (srf.value() & EcalSrFlag::SRF_FORCED_MASK) {
0818       EBForcedRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
0819     }
0820   }
0821   for (int iZ(0); iZ < 2; iZ++) {
0822     for (int iEta(0); iEta < 17; iEta++) {
0823       for (int iPhi(0); iPhi < 72; iPhi++) {
0824         float x(EBtowerCoor[iZ][iEta][iPhi].first);
0825         float y(EBtowerCoor[iZ][iEta][iPhi].second);
0826         fill(meZs1Ru_, x, y, EBZs1RuCount[iZ][iEta][iPhi]);
0827         fill(meFullRoRu_, x, y, EBFullRuCount[iZ][iEta][iPhi]);
0828         fill(meForcedRu_, x, y, EBForcedRuCount[iZ][iEta][iPhi]);
0829       }
0830     }
0831   }
0832 
0833   if (eventError)
0834     srApplicationErrorLog_ << event.id() << ": " << nEbZsErrors_
0835                            << " ZS-flagged EB channels under "
0836                               "the ZS threshold, whose "
0837                            << nEbZsErrorsType1_ << " in a complete RU.\n";
0838 }
0839 
0840 EcalSelectiveReadoutValidation::~EcalSelectiveReadoutValidation() {}
0841 
0842 void EcalSelectiveReadoutValidation::dqmBeginRun(edm::Run const& r, edm::EventSetup const& es) {
0843   // endcap mapping
0844   triggerTowerMap_ = &es.getData(hTriggerTowerMap);
0845 
0846   //electronics map
0847   elecMap_ = &es.getData(ecalmapping);
0848 
0849   initAsciiFile();
0850 }
0851 
0852 void EcalSelectiveReadoutValidation::dqmEndRun(const edm::Run& r, const edm::EventSetup& es) {
0853   meL1aRate_->Fill(getL1aRate());
0854 }
0855 
0856 void EcalSelectiveReadoutValidation::bookHistograms(DQMStore::IBooker& ibooker,
0857                                                     edm::Run const&,
0858                                                     edm::EventSetup const&) {
0859   ibooker.setCurrentFolder("EcalDigisV/SelectiveReadout");
0860 
0861   {
0862     auto scope = DQMStore::IBooker::UseRunScope(ibooker);
0863     meL1aRate_ = bookFloat(ibooker, "l1aRate_");
0864   }
0865 
0866   meDccVol_ = bookProfile(ibooker,
0867                           "hDccVol",  //"EcalDccEventSizeComputed",
0868                           "ECAL DCC event fragment size;Dcc id; "
0869                           "<Event size> (kB)",
0870                           nDccs_,
0871                           .5,
0872                           .5 + nDccs_);
0873 
0874   meDccLiVol_ = bookProfile(ibooker,
0875                             "hDccLiVol",
0876                             "LI channel payload per DCC;Dcc id; "
0877                             "<Event size> (kB)",
0878                             nDccs_,
0879                             .5,
0880                             .5 + nDccs_);
0881 
0882   meDccHiVol_ = bookProfile(ibooker,
0883                             "hDccHiVol",
0884                             "HI channel payload per DCC;Dcc id; "
0885                             "<Event size> (kB)",
0886                             nDccs_,
0887                             .5,
0888                             .5 + nDccs_);
0889 
0890   meDccVolFromData_ = bookProfile(ibooker,
0891                                   "hDccVolFromData",  //"EcalDccEventSize",
0892                                   "ECAL DCC event fragment size;Dcc id; "
0893                                   "<Event size> (kB)",
0894                                   nDccs_,
0895                                   .5,
0896                                   .5 + nDccs_);
0897 
0898   meVolBLI_ = book1D(ibooker,
0899                      "hVolBLI",  // "EBLowInterestPayload",
0900                      "ECAL Barrel low interest crystal data payload;"
0901                      "Event size (kB);Nevts",
0902                      100,
0903                      0.,
0904                      200.);
0905 
0906   meVolELI_ = book1D(ibooker,
0907                      "hVolELI",  //"EELowInterestPayload",
0908                      "Endcap low interest crystal data payload;"
0909                      "Event size (kB);Nevts",
0910                      100,
0911                      0.,
0912                      200.);
0913 
0914   meVolLI_ = book1D(ibooker,
0915                     "hVolLI",  //"EcalLowInterestPayload",
0916                     "ECAL low interest crystal data payload;"
0917                     "Event size (kB);Nevts",
0918                     100,
0919                     0.,
0920                     200.);
0921 
0922   meVolBHI_ = book1D(ibooker,
0923                      "hVolBHI",  //"EBHighInterestPayload",
0924                      "Barrel high interest crystal data payload;"
0925                      "Event size (kB);Nevts",
0926                      100,
0927                      0.,
0928                      200.);
0929 
0930   meVolEHI_ = book1D(ibooker,
0931                      "hVolEHI",  //"EEHighInterestPayload",
0932                      "Endcap high interest crystal data payload;"
0933                      "Event size (kB);Nevts",
0934                      100,
0935                      0.,
0936                      200.);
0937 
0938   meVolHI_ = book1D(ibooker,
0939                     "hVolHI",  //"EcalHighInterestPayload",
0940                     "ECAL high interest crystal data payload;"
0941                     "Event size (kB);Nevts",
0942                     100,
0943                     0.,
0944                     200.);
0945 
0946   meVolB_ = book1D(ibooker,
0947                    "hVolB",  //"EBEventSize",
0948                    "Barrel data volume;Event size (kB);Nevts",
0949                    100,
0950                    0.,
0951                    200.);
0952 
0953   meVolE_ = book1D(ibooker,
0954                    "hVolE",  //"EEEventSize",
0955                    "Endcap data volume;Event size (kB);Nevts",
0956                    100,
0957                    0.,
0958                    200.);
0959 
0960   meVol_ = book1D(ibooker,
0961                   "hVol",  //"EcalEventSize",
0962                   "ECAL data volume;Event size (kB);Nevts",
0963                   100,
0964                   0.,
0965                   200.);
0966 
0967   meChOcc_ = bookProfile2D(ibooker,
0968                            "h2ChOcc",  //"EcalChannelOccupancy",
0969                            "ECAL crystal channel occupancy after zero suppression;"
0970                            "iX -200 / iEta / iX + 100;"
0971                            "iY / iPhi (starting from -10^{o}!);"
0972                            "Event count rate",
0973                            401,
0974                            -200.5,
0975                            200.5,
0976                            360,
0977                            .5,
0978                            360.5);
0979 
0980   //TP
0981   string tpUnit;
0982   if (tpInGeV_)
0983     tpUnit = string("GeV");
0984   else
0985     tpUnit = string("TP hw unit");
0986   string title;
0987   title = string("Trigger primitive TT E_{T};E_{T} ") + tpUnit + string(";Event Count");
0988   meTp_ = bookProfile(ibooker,
0989                       "hTp",  //"EcalTriggerPrimitiveEt",
0990                       title,
0991                       (tpInGeV_ ? 100 : 40),
0992                       0.,
0993                       (tpInGeV_ ? 10. : 40.));
0994 
0995   meTtf_ = bookProfile(ibooker,
0996                        "hTtf",  //"EcalTriggerTowerFlag",
0997                        "Trigger primitive TT flag;Flag number;Event count",
0998                        8,
0999                        -.5,
1000                        7.5);
1001 
1002   title = string("Trigger tower flag vs TP;E_{T}(TT) (") + tpUnit + string(");Flag number");
1003   meTtfVsTp_ = book2D(ibooker, "h2TtfVsTp", title, 100, 0., (tpInGeV_ ? 10. : 40.), 8, -.5, 7.5);
1004 
1005   meTtfVsEtSum_ = book2D(ibooker,
1006                          "h2TtfVsEtSum",
1007                          "Trigger tower flag vs #sumE_{T};"
1008                          "E_{T}(TT) (GeV);"
1009                          "TTF",
1010                          100,
1011                          0.,
1012                          10.,
1013                          8,
1014                          -.5,
1015                          7.5);
1016   title = string(
1017               "Trigger primitive Et (TP) vs #sumE_{T};"
1018               "E_{T} (sum) (GeV);"
1019               "E_{T} (TP) (") +
1020           tpUnit + string(")");
1021 
1022   meTpVsEtSum_ = book2D(ibooker, "h2TpVsEtSum", title, 100, 0., 10., 100, 0., (tpInGeV_ ? 10. : 40.));
1023 
1024   title = string(
1025               "Trigger primitive E_{T};"
1026               "iEta;"
1027               "iPhi;"
1028               "E_{T} (TP) (") +
1029           tpUnit + string(")");
1030   meTpMap_ = bookProfile2D(ibooker, "h2Tp", title, 57, -28.5, 28.5, 72, .5, 72.5);
1031 
1032   //SRF
1033   meFullRoRu_ = book2D(ibooker,
1034                        "h2FRORu",  //"EcalFullReadoutSRFlagMap",
1035                        "Full Read-out readout unit;"
1036                        "iX - 40 / iEta / iX + 20;"
1037                        "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1038                        "Event count",
1039                        80,
1040                        -39.5,
1041                        40.5,
1042                        72,
1043                        .5,
1044                        72.5);
1045 
1046   meFullRoCnt_ = book1D(ibooker,
1047                         "hFROCnt",
1048                         "Number of Full-readout-flagged readout units;"
1049                         "FRO RU count;Event count",
1050                         300,
1051                         -.5,
1052                         299.5);
1053 
1054   meEbFullRoCnt_ = book1D(ibooker,
1055                           "hEbFROCnt",
1056                           "Number of EB Full-readout-flagged readout units;"
1057                           "FRO RU count;Event count",
1058                           200,
1059                           -.5,
1060                           199.5);
1061 
1062   meEeFullRoCnt_ = book1D(ibooker,
1063                           "hEeFROCnt",
1064                           "Number of EE Full-readout-flagged readout units;"
1065                           "FRO RU count;Event count",
1066                           200,
1067                           -.5,
1068                           199.5);
1069 
1070   meZs1Ru_ = book2D(ibooker,
1071                     "h2Zs1Ru",  //"EbZeroSupp1SRFlagMap",
1072                     "Readout unit with ZS-thr-1 flag;"
1073                     "iX - 40 / iEta / iX + 20;"
1074                     "iY0 / iPhi0 (iPhi = 1 at phi = 0 rad);"
1075                     "Event count",
1076                     80,
1077                     -39.5,
1078                     40.5,
1079                     72,
1080                     .5,
1081                     72.5);
1082 
1083   meForcedRu_ = book2D(ibooker,
1084                        "h2ForcedRu",  //"EcalReadoutUnitForcedBitMap",
1085                        "ECAL readout unit with forced bit of SR flag on;"
1086                        "iX - 40 / iEta / iX + 20;"
1087                        "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1088                        "Event count",
1089                        80,
1090                        -39.5,
1091                        40.5,
1092                        72,
1093                        .5,
1094                        72.5);
1095 
1096   meLiTtf_ = bookProfile2D(ibooker,
1097                            "h2LiTtf",  //"EcalLowInterestTriggerTowerFlagMap",
1098                            "Low interest trigger tower flags;"
1099                            "iEta;"
1100                            "iPhi;"
1101                            "Event count",
1102                            57,
1103                            -28.5,
1104                            28.5,
1105                            72,
1106                            .5,
1107                            72.5);
1108 
1109   meMiTtf_ = bookProfile2D(ibooker,
1110                            "h2MiTtf",  //"EcalMidInterestTriggerTowerFlagMap",
1111                            "Mid interest trigger tower flags;"
1112                            "iEta;"
1113                            "iPhi;"
1114                            "Event count",
1115                            57,
1116                            -28.5,
1117                            28.5,
1118                            72,
1119                            .5,
1120                            72.5);
1121 
1122   meHiTtf_ = bookProfile2D(ibooker,
1123                            "h2HiTtf",  //"EcalHighInterestTriggerTowerFlagMap",
1124                            "High interest trigger tower flags;"
1125                            "iEta;"
1126                            "iPhi;"
1127                            "Event count",
1128                            57,
1129                            -28.5,
1130                            28.5,
1131                            72,
1132                            .5,
1133                            72.5);
1134 
1135   meForcedTtf_ = book2D(ibooker,
1136                         "h2ForcedTtf",  //"EcalTtfForcedBitMap",
1137                         "Trigger tower flags with forced bit set;"
1138                         "iEta;"
1139                         "iPhi;"
1140                         "Event count",
1141                         57,
1142                         -28.5,
1143                         28.5,
1144                         72,
1145                         .5,
1146                         72.5);
1147 
1148   const float ebMinNoise = -1.;
1149   const float ebMaxNoise = 1.;
1150 
1151   const float eeMinNoise = -1.;
1152   const float eeMaxNoise = 1.;
1153 
1154   const float ebMinE = ebMinNoise;
1155   const float ebMaxE = ebMaxNoise;
1156 
1157   const float eeMinE = eeMinNoise;
1158   const float eeMaxE = eeMaxNoise;
1159 
1160   const int evtMax = 500;
1161 
1162   meEbRecE_ = book1D(ibooker, "hEbRecE", "Crystal reconstructed energy;E (GeV);Event count", 100, ebMinE, ebMaxE);
1163 
1164   meEbEMean_ = bookProfile(ibooker, "hEbEMean", "EE <E_hit>;event #;<E_hit> (GeV)", evtMax, .5, evtMax + .5);
1165 
1166   meEbNoise_ = book1D(ibooker,
1167                       "hEbNoise",
1168                       "Crystal noise "
1169                       "(rec E of crystal without deposited energy)"
1170                       ";Rec E (GeV);Event count",
1171                       100,
1172                       ebMinNoise,
1173                       ebMaxNoise);
1174 
1175   meEbLiZsFir_ = book1D(ibooker,
1176                         "zsEbLiFIRemu",
1177                         "Emulated ouput of ZS FIR filter for EB "
1178                         "low interest crystals;"
1179                         "ADC count*4;"
1180                         "Event count",
1181                         60,
1182                         -30,
1183                         30);
1184 
1185   meEbHiZsFir_ = book1D(ibooker,
1186                         "zsEbHiFIRemu",
1187                         "Emulated ouput of ZS FIR filter for EB "
1188                         "high interest crystals;"
1189                         "ADC count*4;"
1190                         "Event count",
1191                         60,
1192                         -30,
1193                         30);
1194 
1195   //TODO: Fill this histogram...
1196   //   meEbIncompleteRUZsFir_ = book1D(ibooker, "zsEbIncompleteRUFIRemu",
1197   //                                   "Emulated ouput of ZS FIR filter for EB "
1198   //                                   "incomplete FRO-flagged RU;"
1199   //                                   "ADC count*4;"
1200   //                                   "Event count",
1201   //                                   60, -30, 30);
1202 
1203   meEbSimE_ = book1D(ibooker, "hEbSimE", "EB hit crystal simulated energy", 100, ebMinE, ebMaxE);
1204 
1205   meEbRecEHitXtal_ = book1D(ibooker, "hEbRecEHitXtal", "EB rec energy of hit crystals", 100, ebMinE, ebMaxE);
1206 
1207   meEbRecVsSimE_ = book2D(ibooker,
1208                           "hEbRecVsSimE",
1209                           "Crystal simulated vs reconstructed energy;"
1210                           "Esim (GeV);Erec GeV);Event count",
1211                           100,
1212                           ebMinE,
1213                           ebMaxE,
1214                           100,
1215                           ebMinE,
1216                           ebMaxE);
1217 
1218   meEbNoZsRecVsSimE_ = book2D(ibooker,
1219                               "hEbNoZsRecVsSimE",
1220                               "Crystal no-zs simulated vs reconstructed "
1221                               "energy;"
1222                               "Esim (GeV);Erec GeV);Event count",
1223                               100,
1224                               ebMinE,
1225                               ebMaxE,
1226                               100,
1227                               ebMinE,
1228                               ebMaxE);
1229 
1230   meEeRecE_ = book1D(ibooker,
1231                      "hEeRecE",
1232                      "EE crystal reconstructed energy;E (GeV);"
1233                      "Event count",
1234                      100,
1235                      eeMinE,
1236                      eeMaxE);
1237 
1238   meEeEMean_ = bookProfile(ibooker, "hEeEMean", "<E_{EE hit}>;event;<E_{hit}> (GeV)", evtMax, .5, evtMax + .5);
1239 
1240   meEeNoise_ = book1D(ibooker,
1241                       "hEeNoise",
1242                       "EE crystal noise "
1243                       "(rec E of crystal without deposited energy);"
1244                       "E (GeV);Event count",
1245                       200,
1246                       eeMinNoise,
1247                       eeMaxNoise);
1248 
1249   meEeLiZsFir_ = book1D(ibooker,
1250                         "zsEeLiFIRemu",
1251                         "Emulated ouput of ZS FIR filter for EE "
1252                         "low interest crystals;"
1253                         "ADC count*4;"
1254                         "Event count",
1255                         60,
1256                         -30,
1257                         30);
1258 
1259   meEeHiZsFir_ = book1D(ibooker,
1260                         "zsEeHiFIRemu",
1261                         "Emulated ouput of ZS FIR filter for EE "
1262                         "high interest crystals;"
1263                         "ADC count*4;"
1264                         "Event count",
1265                         60,
1266                         -30,
1267                         30);
1268 
1269   //TODO: Fill this histogram...
1270   //   meEeIncompleteRUZsFir_ = book1D(ibooker, "zsEeIncompleteRUFIRemu",
1271   //                                     "Emulated ouput of ZS FIR filter for EE "
1272   //                                     "incomplete FRO-flagged RU;"
1273   //                                     "ADC count*4;"
1274   //                                     "Event count",
1275   //                                     60, -30, 30);
1276 
1277   meEeSimE_ = book1D(ibooker, "hEeSimE", "EE hit crystal simulated energy", 100, eeMinE, eeMaxE);
1278 
1279   meEeRecEHitXtal_ = book1D(ibooker, "hEeRecEHitXtal", "EE rec energy of hit crystals", 100, eeMinE, eeMaxE);
1280 
1281   meEeRecVsSimE_ = book2D(ibooker,
1282                           "hEeRecVsSimE",
1283                           "EE crystal simulated vs reconstructed energy;"
1284                           "Esim (GeV);Erec GeV);Event count",
1285                           100,
1286                           eeMinE,
1287                           eeMaxE,
1288                           100,
1289                           eeMinE,
1290                           eeMaxE);
1291 
1292   meEeNoZsRecVsSimE_ = book2D(ibooker,
1293                               "hEeNoZsRecVsSimE",
1294                               "EE crystal no-zs simulated vs "
1295                               "reconstructed "
1296                               "energy;Esim (GeV);Erec GeV);Event count",
1297                               100,
1298                               eeMinE,
1299                               eeMaxE,
1300                               100,
1301                               eeMinE,
1302                               eeMaxE);
1303 
1304   meSRFlagsConsistency_ = book2D(ibooker,
1305                                  "hSRAlgoErrorMap",
1306                                  "TTFlags and SR Flags mismatch;"
1307                                  "iX - 40 / iEta / iX + 20;"
1308                                  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1309                                  "Event count",
1310                                  80,
1311                                  -39.5,
1312                                  40.5,
1313                                  72,
1314                                  .5,
1315                                  72.5);
1316 
1317   //Readout Units histos (interest/Ncrystals)
1318   meIncompleteFROMap_ = book2D(ibooker,
1319                                "hIncompleteFROMap",
1320                                "Incomplete full-readout-flagged readout units;"
1321                                "iX - 40 / iEta / iX + 20;"
1322                                "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1323                                "Event count",
1324                                80,
1325                                -39.5,
1326                                40.5,
1327                                72,
1328                                .5,
1329                                72.5);
1330 
1331   meIncompleteFROCnt_ = book1D(ibooker,
1332                                "hIncompleteFROCnt",
1333                                "Number of incomplete full-readout-flagged "
1334                                "readout units;"
1335                                "Number of RUs;Event count;",
1336                                200,
1337                                -.5,
1338                                199.5);
1339 
1340   meIncompleteFRORateMap_ = bookProfile2D(ibooker,
1341                                           "hIncompleteFRORateMap",
1342                                           "Incomplete full-readout-flagged readout units;"
1343                                           "iX - 40 / iEta / iX + 20;"
1344                                           "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1345                                           "Incomplete error rate",
1346                                           80,
1347                                           -39.5,
1348                                           40.5,
1349                                           72,
1350                                           .5,
1351                                           72.5);
1352 
1353   meDroppedFROMap_ = book2D(ibooker,
1354                             "hDroppedFROMap",
1355                             "Dropped full-readout-flagged readout units;"
1356                             "iX - 40 / iEta / iX + 20;"
1357                             "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1358                             "Event count",
1359                             80,
1360                             -39.5,
1361                             40.5,
1362                             72,
1363                             .5,
1364                             72.5);
1365 
1366   meDroppedFROCnt_ = book1D(ibooker,
1367                             "hDroppedFROCnt",
1368                             "Number of dropped full-readout-flagged "
1369                             "RU count;RU count;Event count",
1370                             200,
1371                             -.5,
1372                             199.5);
1373 
1374   meCompleteZSCnt_ = book1D(ibooker,
1375                             "hCompleteZsCnt",
1376                             "Number of zero-suppressed-flagged RU fully "
1377                             "readout;"
1378                             "RU count;Event count",
1379                             200,
1380                             -.5,
1381                             199.5);
1382 
1383   stringstream buf;
1384   buf << "Number of LI EB channels below the " << ebZsThr_ / 4.
1385       << " ADC count ZS threshold;"
1386          "Channel count;Event count",
1387       meEbZsErrCnt_ = book1D(ibooker, "hEbZsErrCnt", buf.str(), 200, -.5, 199.5);
1388 
1389   buf.str("");
1390   buf << "Number of LI EE channels below the " << eeZsThr_ / 4.
1391       << " ADC count ZS theshold;"
1392          "Channel count;Event count",
1393       meEeZsErrCnt_ = book1D(ibooker, "hEeZsErrCnt", buf.str(), 200, -.5, 199.5);
1394 
1395   meZsErrCnt_ = book1D(ibooker,
1396                        "hZsErrCnt",
1397                        "Number of LI channels below the ZS threshold;"
1398                        "Channel count;Event count",
1399                        200,
1400                        -.5,
1401                        199.5);
1402 
1403   meEbZsErrType1Cnt_ = book1D(ibooker,
1404                               "hEbZsErrType1Cnt",
1405                               "Number of EB channels below the ZS "
1406                               "threshold in a LI but fully readout RU;"
1407                               "Channel count;Event count;",
1408                               200,
1409                               -.5,
1410                               199.5);
1411 
1412   meEeZsErrType1Cnt_ = book1D(ibooker,
1413                               "hEeZsErrType1Cnt",
1414                               "Number EE channels below the ZS threshold"
1415                               " in a LI but fully readout RU;"
1416                               "Channel count;Event count",
1417                               200,
1418                               -.5,
1419                               199.5);
1420 
1421   meZsErrType1Cnt_ = book1D(ibooker,
1422                             "hZsErrType1Cnt",
1423                             "Number of LI channels below the ZS threshold "
1424                             "in a LI but fully readout RU;"
1425                             "Channel count;Event count",
1426                             200,
1427                             -.5,
1428                             199.5);
1429 
1430   meDroppedFRORateMap_ = bookProfile2D(ibooker,
1431                                        "hDroppedFRORateMap",
1432                                        "Dropped full-readout-flagged readout units"
1433                                        "iX - 40 / iEta / iX + 20;"
1434                                        "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1435                                        "Dropping rate",
1436                                        80,
1437                                        -39.5,
1438                                        40.5,
1439                                        72,
1440                                        .5,
1441                                        72.5);
1442 
1443   meCompleteZSMap_ = book2D(ibooker,
1444                             "hCompleteZSMap",
1445                             "Complete zero-suppressed-flagged readout units;"
1446                             "iX - 40 / iEta / iX + 20;"
1447                             "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1448                             "Event count",
1449                             80,
1450                             -39.5,
1451                             40.5,
1452                             72,
1453                             .5,
1454                             72.5);
1455 
1456   meCompleteZSRateMap_ = bookProfile2D(ibooker,
1457                                        "hCompleteZSRate",
1458                                        "Complete zero-suppressed-flagged readout units;"
1459                                        "iX - 40 / iEta / iX + 20;"
1460                                        "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1461                                        "Completeness rate",
1462                                        80,
1463                                        -39.5,
1464                                        40.5,
1465                                        72,
1466                                        .5,
1467                                        72.5);
1468 
1469   //print list of available histograms (must be called after
1470   //the bookXX methods):
1471   printAvailableHists();
1472 
1473   //check the histList parameter:
1474   stringstream s;
1475   for (set<string>::iterator it = histList_.begin(); it != histList_.end(); ++it) {
1476     if (*it != string("all") && availableHistList_.find(*it) == availableHistList_.end()) {
1477       s << (s.str().empty() ? "" : ", ") << *it;
1478     }
1479   }
1480   if (!s.str().empty()) {
1481     LogWarning("Configuration") << "Parameter 'histList' contains some unknown histogram(s). "
1482                                    "Check spelling. Following name were not found: "
1483                                 << s.str();
1484   }
1485 }
1486 
1487 void EcalSelectiveReadoutValidation::analyzeTP(edm::Event const& event, edm::EventSetup const& es) {
1488   int TTFlagCount[8];
1489   int LiTTFlagCount[nTtEta][nTtPhi];
1490   int MiTTFlagCount[nTtEta][nTtPhi];
1491   int HiTTFlagCount[nTtEta][nTtPhi];
1492   for (int iTTFlag(0); iTTFlag < 8; iTTFlag++) {
1493     TTFlagCount[iTTFlag] = 0;
1494   }
1495   for (int iTtEta(0); iTtEta < nTtEta; iTtEta++) {
1496     for (int iTtPhi(0); iTtPhi < nTtPhi; iTtPhi++) {
1497       LiTTFlagCount[iTtEta][iTtPhi] = 0;
1498       MiTTFlagCount[iTtEta][iTtPhi] = 0;
1499       HiTTFlagCount[iTtEta][iTtPhi] = 0;
1500     }
1501   }
1502   int tpEtCount[100];
1503   for (int iEt(0); iEt < 100; iEt++) {
1504     tpEtCount[iEt] = 0;
1505   }
1506 
1507   const EcalTPGPhysicsConstMap& physMap = es.getData(physHandle).getMap();
1508 
1509   const EcalTPGGroups::EcalTPGGroupsMap& lutGrpMap = es.getData(lutGrpHandle).getMap();
1510 
1511   const EcalTPGLutIdMap::EcalTPGLutMap& lutMap = es.getData(lutMapHandle).getMap();
1512 
1513   EcalTPGPhysicsConstMapIterator ebItr(physMap.find(DetId(DetId::Ecal, EcalBarrel).rawId()));
1514   double lsb10bitsEB(ebItr == physMap.end() ? 0. : ebItr->second.EtSat / 1024.);
1515   EcalTPGPhysicsConstMapIterator eeItr(physMap.find(DetId(DetId::Ecal, EcalEndcap).rawId()));
1516   double lsb10bitsEE(eeItr == physMap.end() ? 0. : eeItr->second.EtSat / 1024.);
1517 
1518   for (EcalTrigPrimDigiCollection::const_iterator it = tps_->begin(); it != tps_->end(); ++it) {
1519     double tpEt;
1520     if (tpInGeV_) {
1521       EcalTrigTowerDetId const& towerId(it->id());
1522       unsigned int ADC = it->compressedEt();
1523 
1524       double lsb10bits(0.);
1525       if (towerId.subDet() == EcalBarrel)
1526         lsb10bits = lsb10bitsEB;
1527       else if (towerId.subDet() == EcalEndcap)
1528         lsb10bits = lsb10bitsEE;
1529 
1530       int tpg10bits = 0;
1531       EcalTPGGroups::EcalTPGGroupsMapItr itgrp = lutGrpMap.find(towerId.rawId());
1532       uint32_t lutGrp = 999;
1533       if (itgrp != lutGrpMap.end())
1534         lutGrp = itgrp->second;
1535 
1536       EcalTPGLutIdMap::EcalTPGLutMapItr itLut = lutMap.find(lutGrp);
1537       if (itLut != lutMap.end()) {
1538         const unsigned int* lut = (itLut->second).getLut();
1539         for (unsigned int i = 0; i < 1024; i++)
1540           if (ADC == (0xff & lut[i])) {
1541             tpg10bits = i;
1542             break;
1543           }
1544       }
1545 
1546       tpEt = lsb10bits * tpg10bits;
1547     } else {
1548       tpEt = it->compressedEt();
1549     }
1550     int iEta = it->id().ieta();
1551     int iEta0 = iTtEta2cIndex(iEta);
1552     int iPhi = it->id().iphi();
1553     int iPhi0 = iTtPhi2cIndex(iPhi);
1554     double etSum = ttEtSums[iEta0][iPhi0];
1555 
1556     int iE = meTp_->getTProfile()->FindFixBin(tpEt);
1557     if ((iE >= 0) && (iE < 100)) {
1558       ++tpEtCount[iE];
1559     } else {
1560       // FindFixBin might return an overflow bin (outside tpEtCount range).
1561       // To prevent a memory overflow / segfault, these values are ignored.
1562       //std::cout << "EcalSelectiveReadoutValidation: Invalid iE value: " << iE << std::endl;
1563     }
1564 
1565     fill(meTpVsEtSum_, etSum, tpEt);
1566     ++TTFlagCount[it->ttFlag()];
1567     if ((it->ttFlag() & 0x3) == 0) {
1568       LiTTFlagCount[iEta0][iPhi0] += 1;
1569     } else if ((it->ttFlag() & 0x3) == 1) {
1570       MiTTFlagCount[iEta0][iPhi0] += 1;
1571     } else if ((it->ttFlag() & 0x3) == 3) {
1572       HiTTFlagCount[iEta0][iPhi0] += 1;
1573     }
1574     if ((it->ttFlag() & 0x4)) {
1575       fill(meForcedTtf_, iEta, iPhi);
1576     }
1577 
1578     fill(meTtfVsTp_, tpEt, it->ttFlag());
1579     fill(meTtfVsEtSum_, etSum, it->ttFlag());
1580     fill(meTpMap_, iEta, iPhi, tpEt, 1.);
1581   }
1582 
1583   for (int ittflag(0); ittflag < 8; ittflag++) {
1584     fill(meTtf_, ittflag, TTFlagCount[ittflag]);
1585   }
1586   for (int iTtEta(0); iTtEta < nTtEta; iTtEta++) {
1587     for (int iTtPhi(0); iTtPhi < nTtPhi; iTtPhi++) {
1588       fill(meLiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), LiTTFlagCount[iTtEta][iTtPhi]);
1589       fill(meMiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), MiTTFlagCount[iTtEta][iTtPhi]);
1590       fill(meHiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), HiTTFlagCount[iTtEta][iTtPhi]);
1591     }
1592   }
1593   if (tpInGeV_) {
1594     for (int iE(0); iE < 100; iE++) {
1595       fill(meTp_, iE, tpEtCount[iE]);
1596     }
1597   } else {
1598     for (int iE(0); iE < 40; iE++) {
1599       fill(meTp_, iE, tpEtCount[iE]);
1600     }
1601   }
1602 }
1603 
1604 void EcalSelectiveReadoutValidation::analyzeDataVolume(const Event& e, const EventSetup& es) {
1605   anaDigiInit();
1606 
1607   //Complete RU, i.e. RU actually fully readout
1608   for (int iDcc = minDccId_; iDcc <= maxDccId_; ++iDcc) {
1609     for (int iCh = 1; iCh < nDccRus_[iDcc - minDccId_]; ++iCh) {
1610       isRuComplete_[iDcc - minDccId_][iCh - 1] = (nPerRu_[iDcc - minDccId_][iCh - 1] == getCrystalCount(iDcc, iCh));
1611     }
1612   }
1613 
1614   //Barrel
1615   for (unsigned int digis = 0; digis < ebDigis_->size(); ++digis) {
1616     EBDataFrame ebdf = (*ebDigis_)[digis];
1617     anaDigi(ebdf, *ebSrFlags_);
1618   }
1619 
1620   // Endcap
1621   for (unsigned int digis = 0; digis < eeDigis_->size(); ++digis) {
1622     EEDataFrame eedf = (*eeDigis_)[digis];
1623     anaDigi(eedf, *eeSrFlags_);
1624   }
1625 
1626   //histos
1627   for (unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0) {
1628     fill(meDccVol_, iDcc0 + 1, getDccEventSize(iDcc0, nPerDcc_[iDcc0]) / kByte_);
1629     fill(meDccLiVol_, iDcc0 + 1, getDccSrDependentPayload(iDcc0, nLiRuPerDcc_[iDcc0], nLiPerDcc_[iDcc0]) / kByte_);
1630     fill(meDccHiVol_, iDcc0 + 1, getDccSrDependentPayload(iDcc0, nHiRuPerDcc_[iDcc0], nHiPerDcc_[iDcc0]) / kByte_);
1631     const FEDRawDataCollection& raw = *fedRaw_;
1632     fill(meDccVolFromData_, iDcc0 + 1, ((double)raw.FEDData(601 + iDcc0).size()) / kByte_);
1633   }
1634 
1635   //low interesest channels:
1636   double a = nEbLI_ * getBytesPerCrystal() / kByte_;  //getEbEventSize(nEbLI_)/kByte_;
1637   fill(meVolBLI_, a);
1638   double b = nEeLI_ * getBytesPerCrystal() / kByte_;  //getEeEventSize(nEeLI_)/kByte_;
1639   fill(meVolELI_, b);
1640   fill(meVolLI_, a + b);
1641 
1642   //high interest chanels:
1643   a = nEbHI_ * getBytesPerCrystal() / kByte_;  //getEbEventSize(nEbHI_)/kByte_;
1644   fill(meVolBHI_, a);
1645   b = nEeHI_ * getBytesPerCrystal() / kByte_;  //getEeEventSize(nEeHI_)/kByte_;
1646   fill(meVolEHI_, b);
1647   fill(meVolHI_, a + b);
1648 
1649   //any-interest channels:
1650   a = getEbEventSize(nEb_) / kByte_;
1651   fill(meVolB_, a);
1652   b = getEeEventSize(nEe_) / kByte_;
1653   fill(meVolE_, b);
1654   fill(meVol_, a + b);
1655 }
1656 
1657 template <class T, class U>
1658 void EcalSelectiveReadoutValidation::anaDigi(const T& frame, const U& srFlagColl) {
1659   const DetId& xtalId = frame.id();
1660   typedef typename U::key_type RuDetId;
1661   const RuDetId& ruId = readOutUnitOf(frame.id());
1662   typename U::const_iterator srf = srFlagColl.find(ruId);
1663 
1664   bool highInterest = false;
1665   int flag = 0;
1666 
1667   if (srf != srFlagColl.end()) {
1668     flag = srf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
1669 
1670     highInterest = (flag == EcalSrFlag::SRF_FULL);
1671   }
1672 
1673   bool barrel = (xtalId.subdetId() == EcalBarrel);
1674 
1675   pair<int, int> ch = dccCh(xtalId);
1676 
1677   if (barrel) {
1678     ++nEb_;
1679     if (highInterest) {
1680       ++nEbHI_;
1681     } else {  //low interest
1682       ++nEbLI_;
1683     }
1684     int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(xtalId).ieta());
1685     int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(xtalId).iphi());
1686     if (!ebRuActive_[iEta0 / ebTtEdge][iPhi0 / ebTtEdge]) {
1687       ++nRuPerDcc_[ch.first - minDccId_];
1688       if (highInterest) {
1689         ++nHiRuPerDcc_[ch.first - minDccId_];
1690       } else {
1691         ++nLiRuPerDcc_[ch.first - minDccId_];
1692       }
1693 
1694       ebRuActive_[iEta0 / ebTtEdge][iPhi0 / ebTtEdge] = true;
1695     }
1696   } else {  //endcap
1697     ++nEe_;
1698     if (highInterest) {
1699       ++nEeHI_;
1700     } else {  //low interest
1701       ++nEeLI_;
1702     }
1703     int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
1704     int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
1705     int iZ0 = static_cast<const EEDetId&>(frame.id()).zside() > 0 ? 1 : 0;
1706 
1707     if (!eeRuActive_[iZ0][iX0 / scEdge][iY0 / scEdge]) {
1708       ++nRuPerDcc_[ch.first - minDccId_];
1709       if (highInterest) {
1710         ++nHiRuPerDcc_[ch.first - minDccId_];
1711       } else {
1712         ++nLiRuPerDcc_[ch.first - minDccId_];
1713       }
1714 
1715       eeRuActive_[iZ0][iX0 / scEdge][iY0 / scEdge] = true;
1716     }
1717   }
1718 
1719   if (ch.second < 1 || ch.second > 68) {
1720     throw cms::Exception("EcalSelectiveReadoutValidation")
1721         << "Error in DCC channel retrieval for crystal with detId " << xtalId.rawId()
1722         << "DCC channel out of allowed range [1..68]\n";
1723   }
1724   ++nPerDcc_[ch.first - minDccId_];
1725   ++nPerRu_[ch.first - minDccId_][ch.second - 1];
1726   if (highInterest) {
1727     ++nHiPerDcc_[ch.first - minDccId_];
1728   } else {  //low interest channel
1729     ++nLiPerDcc_[ch.first - minDccId_];
1730   }
1731 }
1732 
1733 void EcalSelectiveReadoutValidation::anaDigiInit() {
1734   nEb_ = 0;
1735   nEe_ = 0;
1736   nEeLI_ = 0;
1737   nEeHI_ = 0;
1738   nEbLI_ = 0;
1739   nEbHI_ = 0;
1740   bzero(nPerDcc_, sizeof(nPerDcc_));
1741   bzero(nLiPerDcc_, sizeof(nLiPerDcc_));
1742   bzero(nHiPerDcc_, sizeof(nHiPerDcc_));
1743   bzero(nRuPerDcc_, sizeof(nRuPerDcc_));
1744   bzero(ebRuActive_, sizeof(ebRuActive_));
1745   bzero(eeRuActive_, sizeof(eeRuActive_));
1746   bzero(nPerRu_, sizeof(nPerRu_));
1747   bzero(nLiRuPerDcc_, sizeof(nLiRuPerDcc_));
1748   bzero(nHiRuPerDcc_, sizeof(nHiRuPerDcc_));
1749 }
1750 
1751 double EcalSelectiveReadoutValidation::frame2Energy(const EcalDataFrame& frame) const {
1752   static std::atomic<bool> firstCall{true};
1753   bool expected = true;
1754   if (firstCall.compare_exchange_strong(expected, false)) {
1755     stringstream buf;
1756     buf << "Weights:";
1757     for (unsigned i = 0; i < weights_.size(); ++i) {
1758       buf << "\t" << weights_[i];
1759     }
1760     edm::LogInfo("EcalSrValid") << buf.str() << "\n";
1761     firstCall = false;
1762   }
1763   double adc2GeV = 0.;
1764 
1765   if (typeid(EBDataFrame) == typeid(frame)) {  //barrel APD
1766     adc2GeV = .035;
1767   } else if (typeid(EEDataFrame) == typeid(frame)) {  //endcap VPT
1768     adc2GeV = 0.06;
1769   } else {
1770     assert(false);
1771   }
1772 
1773   double acc = 0;
1774 
1775   const int n = min(frame.size(), (int)weights_.size());
1776 
1777   double gainInv[] = {12., 1., 6., 12.};
1778 
1779   for (int i = 0; i < n; ++i) {
1780     acc += weights_[i] * frame[i].adc() * gainInv[frame[i].gainId()] * adc2GeV;
1781   }
1782   return acc;
1783 }
1784 
1785 int EcalSelectiveReadoutValidation::getRuCount(int iDcc0) const { return nRuPerDcc_[iDcc0]; }
1786 
1787 pair<int, int> EcalSelectiveReadoutValidation::dccCh(const DetId& detId) const {
1788   if (detId.det() != DetId::Ecal) {
1789     throw cms::Exception("InvalidParameter") << "Wrong type of DetId passed to the "
1790                                                 "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
1791                                                 "An ECAL DetId was expected.\n";
1792   }
1793 
1794   DetId xtalId;
1795   switch (detId.subdetId()) {
1796     case EcalTriggerTower:  //Trigger tower
1797     {
1798       const EcalTrigTowerDetId tt = detId;
1799       //pick up one crystal of the trigger tower: they are however all readout by
1800       //the same DCC channel in the barrel.
1801       //Arithmetic is easier on the "c" indices:
1802       const int iTtPhi0 = iTtPhi2cIndex(tt.iphi());
1803       const int iTtEta0 = iTtEta2cIndex(tt.ieta());
1804       const int oneXtalPhi0 = iTtPhi0 * 5;
1805       const int oneXtalEta0 = (iTtEta0 - nOneEeTtEta) * 5;
1806 
1807       xtalId = EBDetId(cIndex2iEta(oneXtalEta0), cIndex2iPhi(oneXtalPhi0));
1808     } break;
1809     case EcalEndcap:
1810       if (detId.rawId() & 0x8000) {  //Supercrystal
1811         return elecMap_->getDCCandSC(EcalScDetId(detId));
1812       } else {  //EE crystal
1813         xtalId = detId;
1814       }
1815       break;
1816     case EcalBarrel:  //EB crystal
1817       xtalId = detId;
1818       break;
1819     default:
1820       throw cms::Exception("InvalidParameter")
1821           << "Wrong type of DetId passed to the method "
1822              "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
1823              "A valid EcalTriggerTower, EcalBarrel or EcalEndcap DetId was expected. "
1824              "detid = "
1825           << xtalId.rawId() << ".\n";
1826   }
1827 
1828   const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
1829 
1830   pair<int, int> result;
1831   result.first = EcalElecId.dccId();
1832 
1833   if (result.first < minDccId_ || result.second > maxDccId_) {
1834     throw cms::Exception("OutOfRange") << "Got an invalid DCC ID, DCCID = " << result.first << " for DetId 0x" << hex
1835                                        << detId.rawId() << " and 0x" << xtalId.rawId() << dec << "\n";
1836   }
1837 
1838   result.second = EcalElecId.towerId();
1839 
1840   if (result.second < 1 || result.second > 68) {
1841     throw cms::Exception("OutOfRange") << "Got an invalid DCC channel ID, DCC_CH = " << result.second << " for DetId 0x"
1842                                        << hex << detId.rawId() << " and 0x" << xtalId.rawId() << dec << "\n";
1843   }
1844 
1845   return result;
1846 }
1847 
1848 EcalTrigTowerDetId EcalSelectiveReadoutValidation::readOutUnitOf(const EBDetId& xtalId) const {
1849   return triggerTowerMap_->towerOf(xtalId);
1850 }
1851 
1852 EcalScDetId EcalSelectiveReadoutValidation::readOutUnitOf(const EEDetId& xtalId) const {
1853   const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
1854   int iDCC = EcalElecId.dccId();
1855   int iDccChan = EcalElecId.towerId();
1856   const bool ignoreSingle = true;
1857   const vector<EcalScDetId> id = elecMap_->getEcalScDetId(iDCC, iDccChan, ignoreSingle);
1858   return !id.empty() ? id[0] : EcalScDetId();
1859 }
1860 
1861 void EcalSelectiveReadoutValidation::setTtEtSums(const edm::EventSetup& es,
1862                                                  const EBDigiCollection& ebDigis,
1863                                                  const EEDigiCollection& eeDigis) {
1864   //ecal geometry:
1865   const CaloSubdetectorGeometry* eeGeometry = nullptr;
1866   const CaloSubdetectorGeometry* ebGeometry = nullptr;
1867   if (eeGeometry == nullptr || ebGeometry == nullptr) {
1868     es.getData(geoToken);
1869     auto geoHandle = es.getHandle(geoToken);
1870     eeGeometry = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
1871     ebGeometry = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
1872   }
1873 
1874   //init etSum array:
1875   for (int iEta0 = 0; iEta0 < nTtEta; ++iEta0) {
1876     for (int iPhi0 = 0; iPhi0 < nTtPhi; ++iPhi0) {
1877       ttEtSums[iEta0][iPhi0] = 0.;
1878     }
1879   }
1880 
1881   for (EBDigiCollection::const_iterator it = ebDigis_->begin(); it != ebDigis_->end(); ++it) {
1882     const EBDataFrame& frame = *it;
1883     const EcalTrigTowerDetId& ttId = triggerTowerMap_->towerOf(frame.id());
1884 
1885     const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
1886     const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
1887     double theta = ebGeometry->getGeometry(frame.id())->getPosition().theta();
1888     double e = frame2EnergyForTp(frame);
1889     if ((frame2EnergyForTp(frame, -1) < e) && (frame2EnergyForTp(frame, 1) < e)) {
1890       ttEtSums[iTtEta0][iTtPhi0] += e * sin(theta);
1891     }
1892   }
1893 
1894   for (EEDigiCollection::const_iterator it = eeDigis.begin(); it != eeDigis.end(); ++it) {
1895     const EEDataFrame& frame = *it;
1896     const EcalTrigTowerDetId& ttId = triggerTowerMap_->towerOf(frame.id());
1897     const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
1898     const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
1899 
1900     double theta = eeGeometry->getGeometry(frame.id())->getPosition().theta();
1901     double e = frame2EnergyForTp(frame);
1902     if ((frame2EnergyForTp(frame, -1) < e) && (frame2EnergyForTp(frame, 1) < e)) {
1903       ttEtSums[iTtEta0][iTtPhi0] += e * sin(theta);
1904     }
1905   }
1906 
1907   //dealing with pseudo-TT in two inner EE eta-ring:
1908   int innerTTEtas[] = {0, 1, 54, 55};
1909   for (unsigned iRing = 0; iRing < sizeof(innerTTEtas) / sizeof(innerTTEtas[0]); ++iRing) {
1910     int iTtEta0 = innerTTEtas[iRing];
1911     //this detector eta-section is divided in only 36 phi bins
1912     //For this eta regions,
1913     //current tower eta numbering scheme is inconsistent. For geometry
1914     //version 133:
1915     //- TT are numbered from 0 to 72 for 36 bins
1916     //- some TT have an even index, some an odd index
1917     //For geometry version 125, there are 72 phi bins.
1918     //The code below should handle both geometry definition.
1919     //If there are 72 input trigger primitives for each inner eta-ring,
1920     //then the average of the trigger primitive of the two pseudo-TT of
1921     //a pair (nEta, nEta+1) is taken as Et of both pseudo TTs.
1922     //If there are only 36 input TTs for each inner eta ring, then half
1923     //of the present primitive of a pseudo TT pair is used as Et of both
1924     //pseudo TTs.
1925 
1926     for (unsigned iTtPhi0 = 0; iTtPhi0 < nTtPhi - 1; iTtPhi0 += 2) {
1927       double et = .5 * (ttEtSums[iTtEta0][iTtPhi0] + ttEtSums[iTtEta0][iTtPhi0 + 1]);
1928       //divides the TT into 2 phi bins in order to match with 72 phi-bins SRP
1929       //scheme or average the Et on the two pseudo TTs if the TT is already
1930       //divided into two trigger primitives.
1931       ttEtSums[iTtEta0][iTtPhi0] = et;
1932       ttEtSums[iTtEta0][iTtPhi0 + 1] = et;
1933     }
1934   }
1935 }
1936 
1937 template <class T>
1938 double EcalSelectiveReadoutValidation::frame2EnergyForTp(const T& frame, int offset) const {
1939   //we have to start by 0 in order to handle offset=-1
1940   //(however Fenix FIR has AFAK only 5 taps)
1941   double weights[] = {0., -1 / 3., -1 / 3., -1 / 3., 0., 1.};
1942 
1943   double adc2GeV = 0.;
1944   if (typeid(frame) == typeid(EBDataFrame)) {
1945     adc2GeV = 0.035;
1946   } else if (typeid(frame) == typeid(EEDataFrame)) {
1947     adc2GeV = 0.060;
1948   } else {  //T is an invalid type!
1949     //TODO: replace message by a cms exception
1950     throw cms::Exception("Severe Error") << __FILE__ << ":" << __LINE__ << ": "
1951                                          << "this is a bug. Please report it.\n";
1952   }
1953 
1954   double acc = 0;
1955 
1956   const int n = min<int>(frame.size(), sizeof(weights) / sizeof(weights[0]));
1957 
1958   double gainInv[] = {12., 1., 6., 12};
1959 
1960   for (int i = offset; i < n; ++i) {
1961     int iframe = i + offset;
1962     if (iframe >= 0 && iframe < frame.size()) {
1963       acc += weights[i] * frame[iframe].adc() * gainInv[frame[iframe].gainId()] * adc2GeV;
1964     }
1965   }
1966   //cout << "\n";
1967   return acc;
1968 }
1969 
1970 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::bookFloat(DQMStore::IBooker& ibook,
1971                                                                                           const std::string& name) {
1972   if (!registerHist(name, ""))
1973     return nullptr;  //this histo is disabled
1974   MonitorElement* result = ibook.bookFloat(name);
1975   if (result == nullptr) {
1976     throw cms::Exception("DQM") << "Failed to book integer DQM monitor element" << name;
1977   }
1978   return result;
1979 }
1980 
1981 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::book1D(
1982     DQMStore::IBooker& ibook, const std::string& name, const std::string& title, int nbins, double xmin, double xmax) {
1983   if (!registerHist(name, title))
1984     return nullptr;  //this histo is disabled
1985   MonitorElement* result = ibook.book1D(name, title, nbins, xmin, xmax);
1986   if (result == nullptr) {
1987     throw cms::Exception("Histo") << "Failed to book histogram " << name;
1988   }
1989   return result;
1990 }
1991 
1992 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::book2D(DQMStore::IBooker& ibook,
1993                                                                                        const std::string& name,
1994                                                                                        const std::string& title,
1995                                                                                        int nxbins,
1996                                                                                        double xmin,
1997                                                                                        double xmax,
1998                                                                                        int nybins,
1999                                                                                        double ymin,
2000                                                                                        double ymax) {
2001   if (!registerHist(name, title))
2002     return nullptr;  //this histo is disabled
2003   MonitorElement* result = ibook.book2D(name, title, nxbins, xmin, xmax, nybins, ymin, ymax);
2004   if (result == nullptr) {
2005     throw cms::Exception("Histo") << "Failed to book histogram " << name;
2006   }
2007   return result;
2008 }
2009 
2010 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::bookProfile(
2011     DQMStore::IBooker& ibook, const std::string& name, const std::string& title, int nbins, double xmin, double xmax) {
2012   if (!registerHist(name, title))
2013     return nullptr;  //this histo is disabled
2014   MonitorElement* result = ibook.bookProfile(name, title, nbins, xmin, xmax, 0, 0, 0);
2015   if (result == nullptr) {
2016     throw cms::Exception("Histo") << "Failed to book histogram " << name;
2017   }
2018   return result;
2019 }
2020 
2021 EcalSelectiveReadoutValidation::MonitorElement* EcalSelectiveReadoutValidation::bookProfile2D(DQMStore::IBooker& ibook,
2022                                                                                               const std::string& name,
2023                                                                                               const std::string& title,
2024                                                                                               int nbinx,
2025                                                                                               double xmin,
2026                                                                                               double xmax,
2027                                                                                               int nbiny,
2028                                                                                               double ymin,
2029                                                                                               double ymax,
2030                                                                                               const char* option) {
2031   if (!registerHist(name, title))
2032     return nullptr;  //this histo is disabled
2033   MonitorElement* result = ibook.bookProfile2D(name, title, nbinx, xmin, xmax, nbiny, ymin, ymax, 0, 0, 0, option);
2034   if (result == nullptr) {
2035     throw cms::Exception("Histo") << "Failed to book histogram " << name;
2036   }
2037   return result;
2038 }
2039 
2040 bool EcalSelectiveReadoutValidation::registerHist(const std::string& name, const std::string& title) {
2041   availableHistList_.insert(pair<string, string>(name, title));
2042   return allHists_ || histList_.find(name) != histList_.end();
2043 }
2044 
2045 void EcalSelectiveReadoutValidation::readAllCollections(const edm::Event& event) {
2046   ebRecHits_.read(event);
2047   eeRecHits_.read(event);
2048   ebDigis_.read(event);
2049   eeDigis_.read(event);
2050   ebNoZsDigis_.read(event);
2051   eeNoZsDigis_.read(event);
2052   ebSrFlags_.read(event);
2053   eeSrFlags_.read(event);
2054   ebComputedSrFlags_.read(event);
2055   eeComputedSrFlags_.read(event);
2056   ebSimHits_.read(event);
2057   eeSimHits_.read(event);
2058   tps_.read(event);
2059   fedRaw_.read(event);
2060 }
2061 
2062 void EcalSelectiveReadoutValidation::printAvailableHists() {
2063   LogInfo log("HistoList");
2064   log << "Avalailable histograms (DQM monitor elements): \n";
2065   for (map<string, string>::iterator it = availableHistList_.begin(); it != availableHistList_.end(); ++it) {
2066     log << it->first << ": " << it->second << "\n";
2067   }
2068   log << "\nTo include an histogram add its name in the vstring parameter "
2069          "'histograms' of the EcalSelectiveReadoutValidation module\n";
2070 }
2071 
2072 double EcalSelectiveReadoutValidation::getEbEventSize(double nReadXtals) const {
2073   double ruHeaderPayload = 0.;
2074   const int firstEbDcc0 = nEeDccs / 2;
2075   for (int iDcc0 = firstEbDcc0; iDcc0 < firstEbDcc0 + nEbDccs; ++iDcc0) {
2076     ruHeaderPayload += getRuCount(iDcc0) * 8.;
2077   }
2078 
2079   return getDccOverhead(EB) * nEbDccs + nReadXtals * getBytesPerCrystal() + ruHeaderPayload;
2080 }
2081 
2082 double EcalSelectiveReadoutValidation::getEeEventSize(double nReadXtals) const {
2083   double ruHeaderPayload = 0.;
2084   const unsigned firstEbDcc0 = nEeDccs / 2;
2085   for (unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0) {
2086     //skip barrel:
2087     if (iDcc0 == firstEbDcc0)
2088       iDcc0 += nEbDccs;
2089     ruHeaderPayload += getRuCount(iDcc0) * 8.;
2090   }
2091   return getDccOverhead(EE) * nEeDccs + nReadXtals * getBytesPerCrystal() + ruHeaderPayload;
2092 }
2093 
2094 //This implementation  assumes that int is coded on at least 28-bits,
2095 //which in pratice should be always true.
2096 int EcalSelectiveReadoutValidation::dccZsFIR(const EcalDataFrame& frame,
2097                                              const std::vector<int>& firWeights,
2098                                              int firstFIRSample,
2099                                              bool* saturated) {
2100   const int nFIRTaps = 6;
2101   //FIR filter weights:
2102   const vector<int>& w = firWeights;
2103 
2104   //accumulator used to compute weighted sum of samples
2105   int acc = 0;
2106   bool gain12saturated = false;
2107   const int gain12 = 0x01;
2108   const int lastFIRSample = firstFIRSample + nFIRTaps - 1;
2109   //LogDebug("DccFir") << "DCC FIR operation: ";
2110   int iWeight = 0;
2111   for (int iSample = firstFIRSample - 1; iSample < lastFIRSample; ++iSample, ++iWeight) {
2112     if (iSample >= 0 && iSample < frame.size()) {
2113       EcalMGPASample sample(frame[iSample]);
2114       if (sample.gainId() != gain12)
2115         gain12saturated = true;
2116       LogTrace("DccFir") << (iSample >= firstFIRSample ? "+" : "") << sample.adc() << "*(" << w[iWeight] << ")";
2117       acc += sample.adc() * w[iWeight];
2118     } else {
2119       edm::LogWarning("DccFir") << __FILE__ << ":" << __LINE__
2120                                 << ": Not enough samples in data frame or 'ecalDccZs1stSample' module "
2121                                    "parameter is not valid...";
2122     }
2123   }
2124   LogTrace("DccFir") << "\n";
2125   //discards the 8 LSBs
2126   //(shift operator cannot be used on negative numbers because
2127   // the result depends on compilator implementation)
2128   acc = (acc >= 0) ? (acc >> 8) : -(-acc >> 8);
2129   //ZS passed if weighted sum acc above ZS threshold or if
2130   //one sample has a lower gain than gain 12 (that is gain 12 output
2131   //is saturated)
2132 
2133   LogTrace("DccFir") << "acc: " << acc << "\n"
2134                      << "saturated: " << (gain12saturated ? "yes" : "no") << "\n";
2135 
2136   if (saturated) {
2137     *saturated = gain12saturated;
2138   }
2139 
2140   return gain12saturated ? numeric_limits<int>::max() : acc;
2141 }
2142 
2143 std::vector<int> EcalSelectiveReadoutValidation::getFIRWeights(const std::vector<double>& normalizedWeights) {
2144   const int nFIRTaps = 6;
2145   vector<int> firWeights(nFIRTaps, 0);  //default weight: 0;
2146   const static int maxWeight = 0xEFF;   //weights coded on 11+1 signed bits
2147   for (unsigned i = 0; i < min((size_t)nFIRTaps, normalizedWeights.size()); ++i) {
2148     firWeights[i] = lround(normalizedWeights[i] * (1 << 10));
2149     if (abs(firWeights[i]) > maxWeight) {  //overflow
2150       firWeights[i] = firWeights[i] < 0 ? -maxWeight : maxWeight;
2151     }
2152   }
2153   return firWeights;
2154 }
2155 
2156 void EcalSelectiveReadoutValidation::configFirWeights(const vector<double>& weightsForZsFIR) {
2157   bool notNormalized = false;
2158   bool notInt = false;
2159   for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2160     if (weightsForZsFIR[i] > 1.)
2161       notNormalized = true;
2162     if ((int)weightsForZsFIR[i] != weightsForZsFIR[i])
2163       notInt = true;
2164   }
2165   if (notInt && notNormalized) {
2166     throw cms::Exception("InvalidParameter") << "weigtsForZsFIR paramater values are not valid: they "
2167                                              << "must either be integer and uses the hardware representation "
2168                                              << "of the weights or less or equal than 1 and used the normalized "
2169                                              << "representation.";
2170   }
2171   LogInfo log("DccFir");
2172   if (notNormalized) {
2173     firWeights_ = vector<int>(weightsForZsFIR.size());
2174     for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2175       firWeights_[i] = (int)weightsForZsFIR[i];
2176     }
2177   } else {
2178     firWeights_ = getFIRWeights(weightsForZsFIR);
2179   }
2180 
2181   log << "Input weights for FIR: ";
2182   for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2183     log << weightsForZsFIR[i] << "\t";
2184   }
2185 
2186   double s2 = 0.;
2187   log << "\nActual FIR weights: ";
2188   for (unsigned i = 0; i < firWeights_.size(); ++i) {
2189     log << firWeights_[i] << "\t";
2190     s2 += firWeights_[i] * firWeights_[i];
2191   }
2192 
2193   s2 = sqrt(s2);
2194   log << "\nNormalized FIR weights after hw representation rounding: ";
2195   for (unsigned i = 0; i < firWeights_.size(); ++i) {
2196     log << firWeights_[i] / (double)(1 << 10) << "\t";
2197   }
2198 
2199   log << "\nFirst FIR sample: " << firstFIRSample_;
2200 }
2201 
2202 void EcalSelectiveReadoutValidation::initAsciiFile() {
2203   if (logSrpAlgoErrors_) {
2204     srpAlgoErrorLog_.open(srpAlgoErrorLogFileName_.c_str(), ios::out | ios::trunc);
2205     if (!srpAlgoErrorLog_.good()) {
2206       throw cms::Exception("Output") << "Failed to open the log file '" << srpAlgoErrorLogFileName_
2207                                      << "' for SRP algorithm result check.\n";
2208     }
2209   }
2210 
2211   if (logSrApplicationErrors_) {
2212     srApplicationErrorLog_.open(srApplicationErrorLogFileName_.c_str(), ios::out | ios::trunc);
2213     if (!srApplicationErrorLog_.good()) {
2214       throw cms::Exception("Output") << "Failed to open the log file '" << srApplicationErrorLogFileName_
2215                                      << "' for Selective Readout decision application check.\n";
2216     }
2217   }
2218 }
2219 
2220 //Compares two SR flag sorted collections . Both collections
2221 //are sorted by their key (the detid) and following algorithm is based on
2222 //this feature.
2223 template <class T>  //T must be either an EBSrFlagCollection or an EESrFlagCollection
2224 void EcalSelectiveReadoutValidation::compareSrfColl(const edm::Event& event, T& srfFromData, T& computedSrf) {
2225   typedef typename T::const_iterator SrFlagCollectionConstIt;
2226   typedef typename T::key_type MyRuDetIdType;
2227   SrFlagCollectionConstIt itSrfFromData = srfFromData.begin();
2228   SrFlagCollectionConstIt itComputedSr = computedSrf.begin();
2229 
2230   while (itSrfFromData != srfFromData.end() || itComputedSr != computedSrf.end()) {
2231     MyRuDetIdType inconsistentRu = 0;
2232     bool inconsistent = false;
2233     if (itComputedSr == computedSrf.end() ||
2234         (itSrfFromData != srfFromData.end() && itSrfFromData->id() < itComputedSr->id())) {
2235       //computedSrf is missig a detid found in srfFromData
2236       pair<int, int> ch = dccCh(itSrfFromData->id());
2237       srpAlgoErrorLog_ << event.id() << ": " << itSrfFromData->id() << ", DCC " << ch.first << " ch " << ch.second
2238                        << " found in data (SRF:" << itSrfFromData->flagName()
2239                        << ") but not in the set of SRFs computed from the data TTF.\n";
2240       inconsistentRu = itSrfFromData->id();
2241       inconsistent = true;
2242       ++itSrfFromData;
2243     } else if (itSrfFromData == srfFromData.end() ||
2244                (itComputedSr != computedSrf.end() && itComputedSr->id() < itSrfFromData->id())) {
2245       //ebSrFlags is missing a detid found in computedSrf
2246       pair<int, int> ch = dccCh(itComputedSr->id());
2247       if (logErrForDccs_[ch.first - minDccId_]) {
2248         srpAlgoErrorLog_ << event.id() << ": " << itComputedSr->id() << ", DCC " << ch.first << " ch " << ch.second
2249                          << " not found in data. Computed SRF: " << itComputedSr->flagName() << ".\n";
2250         inconsistentRu = itComputedSr->id();
2251         inconsistent = true;
2252       }
2253       ++itComputedSr;
2254     } else {
2255       //*itSrfFromData and *itComputedSr has same detid
2256       if (itComputedSr->value() != itSrfFromData->value()) {
2257         pair<int, int> ch = dccCh(itSrfFromData->id());
2258         srpAlgoErrorLog_ << event.id() << ", " << itSrfFromData->id() << ", DCC " << ch.first << " ch " << ch.second
2259                          << ", SRF inconsistency: "
2260                          << "from data: " << itSrfFromData->flagName()
2261                          << ", computed from TTF: " << itComputedSr->flagName() << "\n";
2262         inconsistentRu = itComputedSr->id();
2263         inconsistent = true;
2264       }
2265       if (itComputedSr != computedSrf.end())
2266         ++itComputedSr;
2267       if (itSrfFromData != srfFromData.end())
2268         ++itSrfFromData;
2269     }
2270 
2271     if (inconsistent)
2272       fill(meSRFlagsConsistency_, ruGraphX(inconsistentRu), ruGraphY(inconsistentRu));
2273   }
2274 }
2275 
2276 int EcalSelectiveReadoutValidation::dccId(const EcalScDetId& detId) const { return elecMap_->getDCCandSC(detId).first; }
2277 
2278 int EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId& detId) const {
2279   if (detId.ietaAbs() > 17) {
2280     throw cms::Exception("InvalidArgument")
2281         << "Argument of EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId&) "
2282         << "must be a barrel trigger tower Id\n";
2283   }
2284   return dccCh(detId).first;
2285 }
2286 
2287 void EcalSelectiveReadoutValidation::selectFedsForLog() {
2288   logErrForDccs_ = vector<bool>(nDccs_, false);
2289 
2290   for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
2291     int iDcc = dccId(it->id()) - minDccId_;
2292 
2293     logErrForDccs_.at(iDcc) = true;
2294   }
2295 
2296   for (EESrFlagCollection::const_iterator it = eeSrFlags_->begin(); it != eeSrFlags_->end(); ++it) {
2297     int iDcc = dccId(it->id()) - minDccId_;
2298 
2299     logErrForDccs_.at(iDcc) = true;
2300   }
2301 
2302   stringstream buf;
2303   buf << "List of DCCs found in the first processed event: ";
2304   bool first = true;
2305   for (unsigned iDcc = 0; iDcc < nDccs_; ++iDcc) {
2306     if (logErrForDccs_[iDcc]) {
2307       buf << (first ? "" : ", ") << (iDcc + minDccId_);
2308       first = false;
2309     }
2310   }
2311   buf << "\nOnly DCCs from this list will be considered for error logging\n";
2312   srpAlgoErrorLog_ << buf.str();
2313   srApplicationErrorLog_ << buf.str();
2314   LogInfo("EcalSrValid") << buf.str();
2315 }
2316 
2317 template <class T>
2318 void EcalSelectiveReadoutValidation::checkSrApplication(const edm::Event& event, T& srfs) {
2319   typedef typename T::const_iterator SrFlagCollectionConstIt;
2320   typedef typename T::key_type MyRuDetIdType;
2321 
2322   for (SrFlagCollectionConstIt itSrf = srfs.begin(); itSrf != srfs.end(); ++itSrf) {
2323     int flag = itSrf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
2324     pair<int, int> ru = dccCh(itSrf->id());
2325 
2326     if (flag == EcalSrFlag::SRF_FULL) {
2327       if (nPerRu_[ru.first - minDccId_][ru.second - 1] == getCrystalCount(ru.first, ru.second)) {  //no error
2328         fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2329         fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2330       } else if (nPerRu_[ru.first - minDccId_][ru.second - 1] == 0) {  //tower dropped!
2331         fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2332         fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2333         fill(meDroppedFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2334         ++nDroppedFRO_;
2335         srApplicationErrorLog_ << event.id() << ": Flag of RU " << itSrf->id() << " (DCC " << ru.first << " ch "
2336                                << ru.second << ") is 'Full readout' "
2337                                << "while none of its channel was read out\n";
2338       } else {  //tower partially read out
2339         fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2340         fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2341         fill(meIncompleteFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2342         ++nIncompleteFRO_;
2343         srApplicationErrorLog_ << event.id() << ": Flag of RU" << itSrf->id() << " (DCC " << ru.first << " ch "
2344                                << ru.second << ") is 'Full readout' "
2345                                << "while only " << nPerRu_[ru.first - minDccId_][ru.second - 1] << " / "
2346                                << getCrystalCount(ru.first, ru.second) << " channels were read out.\n";
2347       }
2348     }
2349 
2350     if (flag == EcalSrFlag::SRF_ZS1 || flag == EcalSrFlag::SRF_ZS2) {
2351       if (nPerRu_[ru.first - minDccId_][ru.second - 1] == getCrystalCount(ru.first, ru.second)) {
2352         //ZS readout unit whose every channel was read
2353 
2354         fill(meCompleteZSMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()));
2355         fill(meCompleteZSRateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2356 
2357         ++nCompleteZS_;
2358       } else {
2359         fill(meCompleteZSRateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2360       }
2361     }
2362   }
2363 }
2364 
2365 int EcalSelectiveReadoutValidation::getCrystalCount(int iDcc, int iDccCh) {
2366   if (iDcc < minDccId_ || iDcc > maxDccId_) {  //invalid DCC
2367     return 0;
2368   } else if (10 <= iDcc && iDcc <= 45) {  //EB
2369     return 25;
2370   } else {  //EE
2371     int iDccPhi;
2372     if (iDcc < 10)
2373       iDccPhi = iDcc;
2374     else
2375       iDccPhi = iDcc - 45;
2376     switch (iDccPhi * 100 + iDccCh) {
2377       case 110:
2378       case 232:
2379       case 312:
2380       case 412:
2381       case 532:
2382       case 610:
2383       case 830:
2384       case 806:
2385         //inner partials at 12, 3, and 9 o'clock
2386         return 20;
2387       case 134:
2388       case 634:
2389       case 827:
2390       case 803:
2391         return 10;
2392       case 330:
2393       case 430:
2394         return 20;
2395       case 203:
2396       case 503:
2397       case 721:
2398       case 921:
2399         return 21;
2400       default:
2401         return 25;
2402     }
2403   }
2404 }