Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-13 23:03:34

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