CMS 3D CMS Logo

EcalSelectiveReadoutValidation.cc
Go to the documentation of this file.
1 /*
2  * \file EcalSelectiveReadoutValidation.cc
3  *
4  *
5  */
6 
8 
9 #include "ecalDccMap.h"
10 
17 
25 
27 
28 using namespace cms;
29 using namespace edm;
30 using namespace std;
31 
32 const double EcalSelectiveReadoutValidation::rad2deg = 45. / atan(1.);
33 
34 const int EcalSelectiveReadoutValidation::nDccRus_[nDccs_] = {
35  //EE- DCCs:
36  // 1 2 3 4 5 6 7 8 9
37  34,
38  32,
39  33,
40  33,
41  32,
42  34,
43  33,
44  34,
45  33,
46  //EB- DCCs:
47  // 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
48  68,
49  68,
50  68,
51  68,
52  68,
53  68,
54  68,
55  68,
56  68,
57  68,
58  68,
59  68,
60  68,
61  68,
62  68,
63  68,
64  68,
65  68,
66  //EB+ DCCs:
67  // 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
68  68,
69  68,
70  68,
71  68,
72  68,
73  68,
74  68,
75  68,
76  68,
77  68,
78  68,
79  68,
80  68,
81  68,
82  68,
83  68,
84  68,
85  68,
86  //EE+ DCCs:
87  // 46 47 48 49 50 51 52 53 54
88  32,
89  33,
90  33,
91  32,
92  34,
93  33,
94  34,
95  33,
96  34};
97 
99  : geoToken(esConsumes()),
100  ecalmapping(esConsumes<edm::Transition::BeginRun>()),
101  hTriggerTowerMap(esConsumes<edm::Transition::BeginRun>()),
102  physHandle(esConsumes()),
103  lutGrpHandle(esConsumes()),
104  lutMapHandle(esConsumes()),
105  collNotFoundWarn_(ps.getUntrackedParameter<bool>("warnIfCollectionNotFound", true)),
106  ebDigis_(ps.getParameter<edm::InputTag>("EbDigiCollection"), false, collNotFoundWarn_),
107  eeDigis_(ps.getParameter<edm::InputTag>("EeDigiCollection"), false, collNotFoundWarn_),
108  ebNoZsDigis_(ps.getParameter<edm::InputTag>("EbUnsuppressedDigiCollection"), false, false /*collNotFoundWarn_*/),
109  eeNoZsDigis_(ps.getParameter<edm::InputTag>("EeUnsuppressedDigiCollection"), false, false /*collNotFoundWarn_*/),
110  ebSrFlags_(ps.getParameter<edm::InputTag>("EbSrFlagCollection"), false, collNotFoundWarn_),
111  eeSrFlags_(ps.getParameter<edm::InputTag>("EeSrFlagCollection"), false, collNotFoundWarn_),
112  ebComputedSrFlags_(
113  ps.getParameter<edm::InputTag>("EbSrFlagFromTTCollection"), false, false /*collNotFoundWarn_*/),
114  eeComputedSrFlags_(
115  ps.getParameter<edm::InputTag>("EeSrFlagFromTTCollection"), false, false /*collNotFoundWarn_*/),
116  ebSimHits_(ps.getParameter<edm::InputTag>("EbSimHitCollection"), false, false /*collNotFoundWarn_*/),
117  eeSimHits_(ps.getParameter<edm::InputTag>("EeSimHitCollection"), false, false /*collNotFoundWarn_*/),
118  tps_(ps.getParameter<edm::InputTag>("TrigPrimCollection"), false, collNotFoundWarn_),
119  ebRecHits_(ps.getParameter<edm::InputTag>("EbRecHitCollection"), false, false /*collNotFoundWarn_*/),
120  eeRecHits_(ps.getParameter<edm::InputTag>("EeRecHitCollection"), false, false /*collNotFoundWarn_*/),
121  fedRaw_(ps.getParameter<edm::InputTag>("FEDRawCollection"), false, false /*collNotFoundWarn_*/),
122  tmax(0),
123  tmin(numeric_limits<int64_t>::max()),
124  l1aOfTmin(0),
125  l1aOfTmax(0),
126  triggerTowerMap_(nullptr),
127  localReco_(ps.getParameter<bool>("LocalReco")),
128  weights_(ps.getParameter<vector<double>>("weights")),
129  tpInGeV_(ps.getParameter<bool>("tpInGeV")),
130  firstFIRSample_(ps.getParameter<int>("ecalDccZs1stSample")),
131  useEventRate_(ps.getParameter<bool>("useEventRate")),
132  logErrForDccs_(nDccs_, false),
133  ievt_(0),
134  allHists_(false),
135  histDir_(ps.getParameter<string>("histDir")),
136  withEeSimHit_(false),
137  withEbSimHit_(false) {
139  ebDigis_.setToken(collector);
140  eeDigis_.setToken(collector);
141  ebNoZsDigis_.setToken(collector);
142  eeNoZsDigis_.setToken(collector);
143  ebSrFlags_.setToken(collector);
144  eeSrFlags_.setToken(collector);
145  ebComputedSrFlags_.setToken(collector);
146  eeComputedSrFlags_.setToken(collector);
147  ebSimHits_.setToken(collector);
148  eeSimHits_.setToken(collector);
149  tps_.setToken(collector);
150  ebRecHits_.setToken(collector);
151  eeRecHits_.setToken(collector);
152  fedRaw_.setToken(collector);
153 
154  double ebZsThr = ps.getParameter<double>("ebZsThrADCCount");
155  double eeZsThr = ps.getParameter<double>("eeZsThrADCCount");
156 
157  ebZsThr_ = lround(ebZsThr * 4);
158  eeZsThr_ = lround(eeZsThr * 4);
159 
160  //File to log SRP algorithem inconsistency
161  srpAlgoErrorLogFileName_ = ps.getUntrackedParameter<string>("srpAlgoErrorLogFile", "");
163 
164  //File to log SRP decision application inconsistency
165  srApplicationErrorLogFileName_ = ps.getUntrackedParameter<string>("srApplicationErrorLogFile", "");
167 
168  //FIR ZS weights
169  configFirWeights(ps.getParameter<vector<double>>("dccWeights"));
170 
171  // DQM ROOT output
172  outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
173 
174  if (!outputFile_.empty()) {
175  LogInfo("OutputInfo") << " Ecal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
176  } else {
177  LogInfo("OutputInfo") << " Ecal Digi Task histograms will NOT be saved";
178  }
179 
180  // verbosity switch
181  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
182 
183  // get hold of back-end interface
184 
185  vector<string> hists(ps.getUntrackedParameter<vector<string>>("histograms", vector<string>(1, "all")));
186 
187  for (vector<string>::iterator it = hists.begin(); it != hists.end(); ++it)
188  histList_.insert(*it);
189  if (histList_.find("all") != histList_.end())
190  allHists_ = true;
191 }
192 
194  const int32_t bx = event.bunchCrossing();
195  if (bx < 1 || bx > 3564)
196  return;
197 
198  int64_t t = event.bunchCrossing() + (event.orbitNumber() - 1) * 3564;
199 
200  if (t < tmin) {
201  tmin = t;
202  l1aOfTmin = event.id().event();
203  }
204 
205  if (t > tmax) {
206  tmax = t;
207  l1aOfTmax = event.id().event();
208  }
209 }
210 
212  LogDebug("EcalSrValid") << __FILE__ << ":" << __LINE__ << ": "
213  << "Tmax = " << tmax << " x 25ns; Tmin = " << tmin << " x 25ns; L1A(Tmax) = " << l1aOfTmax
214  << "; L1A(Tmin) = " << l1aOfTmin << "\n";
215  return (double)(l1aOfTmax - l1aOfTmin) / ((tmax - tmin) * 25e-9);
216 }
217 
220 
221  //retrieves event products:
223 
224  withEeSimHit_ = (!eeSimHits_->empty());
225  withEbSimHit_ = (!ebSimHits_->empty());
226 
227  if (ievt_ < 10) {
228  edm::LogInfo("EcalSrValid") << "Size of TP collection: " << tps_->size() << std::endl
229  << "Size of EB SRF collection read from data: " << ebSrFlags_->size() << std::endl
230  << "Size of EB SRF collection computed from data TTFs: " << ebComputedSrFlags_->size()
231  << std::endl
232  << "Size of EE SRF collection read from data: " << eeSrFlags_->size() << std::endl
233  << "Size of EE SRF collection computed from data TTFs: " << eeComputedSrFlags_->size()
234  << std::endl;
235  }
236 
237  if (ievt_ == 0) {
238  selectFedsForLog(); //note: must be called after readAllCollection
239  }
240 
241  //computes Et sum trigger tower crystals:
243 
244  //Data Volume
246 
247  //EB digis
248  //must be called after analyzeDataVolume because it uses
249  //isRuComplete_ array that this method fills
250  analyzeEB(event, es);
251 
252  //EE digis
253  //must be called after analyzeDataVolume because it uses
254  //isRuComplete_ array that this method fills
255  analyzeEE(event, es);
256 
260 
264 
268 
269  //TP
270  analyzeTP(event, es);
271 
272  if (!ebComputedSrFlags_->empty()) {
274  }
275  if (!eeComputedSrFlags_->empty()) {
277  }
278  nDroppedFRO_ = 0;
279  nIncompleteFRO_ = 0;
280  nCompleteZS_ = 0;
286  ++ievt_;
287 }
288 
290  bool eventError = false;
291  nEeZsErrors_ = 0;
292  nEeZsErrorsType1_ = 0;
293 
297  std::unique_ptr<std::array<std::array<std::array<energiesEe_t, nEeY>, nEeX>, nEndcaps>> eeEnergies =
298  std::make_unique<std::array<std::array<std::array<energiesEe_t, nEeY>, nEeX>, nEndcaps>>();
299 
300  for (int iZ0 = 0; iZ0 < nEndcaps; ++iZ0) {
301  for (int iX0 = 0; iX0 < nEeX; ++iX0) {
302  for (int iY0 = 0; iY0 < nEeY; ++iY0) {
303  (*eeEnergies)[iZ0][iX0][iY0].noZsRecE = -numeric_limits<double>::max();
304  (*eeEnergies)[iZ0][iX0][iY0].recE = -numeric_limits<double>::max();
305  (*eeEnergies)[iZ0][iX0][iY0].simE = 0; //must be set to zero.
306  (*eeEnergies)[iZ0][iX0][iY0].simHit = 0;
307  (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
308  }
309  }
310  }
311 
312  // gets the endcap geometry:
313  auto geoHandle = es.getHandle(geoToken);
314  const CaloSubdetectorGeometry* geometry_p = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
315  //CaloSubdetectorGeometry const& geometry = *geometry_p;
316 
317  //EE unsupressed digis:
318  for (unsigned int digis = 0; digis < eeNoZsDigis_->size(); ++digis) {
319  EEDataFrame frame = (*eeNoZsDigis_)[digis];
320  int iX0 = iXY2cIndex(frame.id().ix());
321  int iY0 = iXY2cIndex(frame.id().iy());
322  int iZ0 = frame.id().zside() > 0 ? 1 : 0;
323 
324  if (iX0 < 0 || iX0 >= nEeX) {
325  edm::LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
326  << "[0," << nEeX - 1 << "]\n";
327  }
328  if (iY0 < 0 || iY0 >= nEeY) {
329  edm::LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
330  << "[0," << nEeY - 1 << "]\n";
331  }
332  // cout << "EE no ZS energy computation..." ;
333  (*eeEnergies)[iZ0][iX0][iY0].noZsRecE = frame2Energy(frame);
334 
335  (*eeEnergies)[iZ0][iX0][iY0].gain12 = true;
336  for (int i = 0; i < frame.size(); ++i) {
337  const int gain12Code = 0x1;
338  if (frame[i].gainId() != gain12Code)
339  (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
340  }
341 
342  const GlobalPoint xtalPos = geometry_p->getGeometry(frame.id())->getPosition();
343 
344  (*eeEnergies)[iZ0][iX0][iY0].phi = rad2deg * ((double)xtalPos.phi());
345  (*eeEnergies)[iZ0][iX0][iY0].eta = xtalPos.eta();
346  }
347 
348  //EE rec hits:
349  if (!localReco_) {
350  for (RecHitCollection::const_iterator it = eeRecHits_->begin(); it != eeRecHits_->end(); ++it) {
351  const RecHit& hit = *it;
352  int iX0 = iXY2cIndex(static_cast<const EEDetId&>(hit.id()).ix());
353  int iY0 = iXY2cIndex(static_cast<const EEDetId&>(hit.id()).iy());
354  int iZ0 = static_cast<const EEDetId&>(hit.id()).zside() > 0 ? 1 : 0;
355 
356  if (iX0 < 0 || iX0 >= nEeX) {
357  LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
358  << "[0," << nEeX - 1 << "]\n";
359  }
360  if (iY0 < 0 || iY0 >= nEeY) {
361  LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
362  << "[0," << nEeY - 1 << "]\n";
363  }
364  // cout << "EE no ZS energy computation..." ;
365  (*eeEnergies)[iZ0][iX0][iY0].recE = hit.energy();
366  }
367  }
368 
369  //EE sim hits:
370  for (vector<PCaloHit>::const_iterator it = eeSimHits_->begin(); it != eeSimHits_->end(); ++it) {
371  const PCaloHit& simHit = *it;
372  EEDetId detId(simHit.id());
373  int iX = detId.ix();
374  int iX0 = iXY2cIndex(iX);
375  int iY = detId.iy();
376  int iY0 = iXY2cIndex(iY);
377  int iZ0 = detId.zside() > 0 ? 1 : 0;
378  (*eeEnergies)[iZ0][iX0][iY0].simE += simHit.energy();
379  ++(*eeEnergies)[iZ0][iX0][iY0].simHit;
380  }
381 
382  bool EEcrystalShot[nEeX][nEeY][2];
383  pair<int, int> EExtalCoor[nEeX][nEeY][2];
384 
385  for (int iEeZ = 0; iEeZ < 2; ++iEeZ) {
386  for (int iEeX = 0; iEeX < nEeX; ++iEeX) {
387  for (int iEeY = 0; iEeY < nEeY; ++iEeY) {
388  EEcrystalShot[iEeX][iEeY][iEeZ] = false;
389  EExtalCoor[iEeX][iEeY][iEeZ] = make_pair(0, 0);
390  }
391  }
392  }
393 
394  //EE suppressed digis
396  const EEDataFrame& frame = *it;
397  int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
398  int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
399  int iZ0 = static_cast<const EEDetId&>(frame.id()).zside() > 0 ? 1 : 0;
400  if (iX0 < 0 || iX0 >= nEeX) {
401  LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
402  << "[0," << nEeX - 1 << "]\n";
403  }
404  if (iY0 < 0 || iY0 >= nEeY) {
405  LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
406  << "[0," << nEeY - 1 << "]\n";
407  }
408 
409  if (!EEcrystalShot[iX0][iY0][iZ0]) {
410  EEcrystalShot[iX0][iY0][iZ0] = true;
411  EExtalCoor[iX0][iY0][iZ0] = make_pair(xtalGraphX(frame.id()), xtalGraphY(frame.id()));
412  } else {
413  cout << "Error: several digi for same crystal!";
414  abort();
415  }
416 
417  if (localReco_) {
418  (*eeEnergies)[iZ0][iX0][iY0].recE = frame2Energy(frame);
419  }
420 
421  (*eeEnergies)[iZ0][iX0][iY0].gain12 = true;
422  for (int i = 0; i < frame.size(); ++i) {
423  const int gain12Code = 0x1;
424  if (frame[i].gainId() != gain12Code) {
425  (*eeEnergies)[iZ0][iX0][iY0].gain12 = false;
426  }
427  }
428 
430 
431  bool highInterest = false;
432 
433  if (srf == eeSrFlags_->end())
434  continue;
435 
436  if (srf != eeSrFlags_->end()) {
437  highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK) == EcalSrFlag::SRF_FULL);
438  }
439 
440  if (highInterest) {
442  } else {
443  int v = dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr);
444  fill(meEeLiZsFir_, v);
445  if (v < eeZsThr_) {
446  eventError = true;
447  ++nEeZsErrors_;
448  pair<int, int> ru = dccCh(frame.id());
449  if (isRuComplete_[ru.first - 1][ru.second - 1])
451  if (nEeZsErrors_ < 3) {
452  srApplicationErrorLog_ << event.id() << ", "
453  << "RU " << frame.id() << ", "
454  << "DCC " << ru.first << " Ch : " << ru.second << ": "
455  << "LI channel under ZS threshold.\n";
456  }
457  if (nEeZsErrors_ == 3) {
458  srApplicationErrorLog_ << event.id() << ": "
459  << "more ZS errors for this event...\n";
460  }
461  }
462  }
463  } //next ZS digi.
464 
465  for (int iEeZ = 0; iEeZ < 2; ++iEeZ) {
466  for (int iEeX = 0; iEeX < nEeX; ++iEeX) {
467  for (int iEeY = 0; iEeY < nEeY; ++iEeY) {
468  fill(meChOcc_,
469  EExtalCoor[iEeX][iEeY][iEeZ].first,
470  EExtalCoor[iEeX][iEeY][iEeZ].second,
471  EEcrystalShot[iEeX][iEeY][iEeZ] ? 1 : 0);
472  }
473  }
474  }
475 
476  for (int iZ0 = 0; iZ0 < nEndcaps; ++iZ0) {
477  for (int iX0 = 0; iX0 < nEeX; ++iX0) {
478  for (int iY0 = 0; iY0 < nEeY; ++iY0) {
479  double recE = (*eeEnergies)[iZ0][iX0][iY0].recE;
480  if (recE == -numeric_limits<double>::max())
481  continue; //not a crystal or ZS
482  fill(meEeRecE_, (*eeEnergies)[iZ0][iX0][iY0].recE);
483 
484  fill(meEeEMean_, ievt_ + 1, (*eeEnergies)[iZ0][iX0][iY0].recE);
485 
486  if (withEeSimHit_) {
487  if (!(*eeEnergies)[iZ0][iX0][iY0].simHit) { //noise only crystal channel
488  fill(meEeNoise_, (*eeEnergies)[iZ0][iX0][iY0].noZsRecE);
489  } else {
490  fill(meEeSimE_, (*eeEnergies)[iZ0][iX0][iY0].simE);
491  fill(meEeRecEHitXtal_, (*eeEnergies)[iZ0][iX0][iY0].recE);
492  }
493  fill(meEeRecVsSimE_, (*eeEnergies)[iZ0][iX0][iY0].simE, (*eeEnergies)[iZ0][iX0][iY0].recE);
494  fill(meEeNoZsRecVsSimE_, (*eeEnergies)[iZ0][iX0][iY0].simE, (*eeEnergies)[iZ0][iX0][iY0].noZsRecE);
495  }
496  }
497  }
498  }
499 
500  int EEZs1RuCount[2][20][20];
501  int EEFullRuCount[2][20][20];
502  int EEForcedRuCount[2][20][20];
503  for (int iZ(0); iZ < 2; iZ++) {
504  for (int iX(0); iX < 20; iX++) {
505  for (int iY(0); iY < 20; iY++) {
506  EEZs1RuCount[iZ][iX][iY] = 0;
507  EEFullRuCount[iZ][iX][iY] = 0;
508  EEForcedRuCount[iZ][iX][iY] = 0;
509  }
510  }
511  }
512 
513  nEeFROCnt_ = 0;
514  char eeSrfMark[2][20][20];
515  bzero(eeSrfMark, sizeof(eeSrfMark));
516  //Filling RU histo
517  for (EESrFlagCollection::const_iterator it = eeSrFlags_->begin(); it != eeSrFlags_->end(); ++it) {
518  const EESrFlag& srf = *it;
519  // srf.id() is EcalScDetId; 1 <= ix <= 20 1 <= iy <= 20
520  int iX = srf.id().ix();
521  int iY = srf.id().iy();
522  int zside = srf.id().zside(); //-1 for EE-, +1 for EE+
523  if (iX < 1 || iY > 100)
524  throw cms::Exception("EcalSelectiveReadoutValidation")
525  << "Found an endcap SRF with an invalid det ID: " << srf.id() << ".\n";
526  ++eeSrfMark[zside > 0 ? 1 : 0][iX - 1][iY - 1];
527  if (eeSrfMark[zside > 0 ? 1 : 0][iX - 1][iY - 1] > 1)
528  throw cms::Exception("EcalSelectiveReadoutValidation") << "Duplicate SRF for supercrystal " << srf.id() << ".\n";
529  int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
530  if (flag == EcalSrFlag::SRF_ZS1) {
531  EEZs1RuCount[zside > 0 ? 1 : 0][iX - 1][iY - 1] += 1;
532  }
533 
534  if (flag == EcalSrFlag::SRF_FULL) {
535  EEFullRuCount[zside > 0 ? 1 : 0][iX - 1][iY - 1] += 1;
536  ++nEeFROCnt_;
537  }
538 
539  if (srf.value() & EcalSrFlag::SRF_FORCED_MASK) {
540  EEForcedRuCount[zside > 0 ? 1 : 0][iX - 1][iY - 1] += 1;
541  }
542  }
543  for (int iZ(0); iZ < 2; iZ++) {
544  int xOffset(iZ == 0 ? -40 : 20);
545  for (int iX(0); iX < 20; iX++) {
546  for (int iY(0); iY < 20; iY++) {
547  int GraphX = (iX + 1) + xOffset;
548  int GraphY = (iY + 1);
549  fill(meZs1Ru_, GraphX, GraphY, EEZs1RuCount[iZ][iX][iY]);
550  fill(meFullRoRu_, GraphX, GraphY, EEFullRuCount[iZ][iX][iY]);
551  fill(meForcedRu_, GraphX, GraphY, EEForcedRuCount[iZ][iX][iY]);
552  }
553  }
554  }
555 
556  if (eventError)
557  srApplicationErrorLog_ << event.id() << ": " << nEeZsErrors_
558  << " ZS-flagged EE channels under "
559  "the ZS threshold, whose "
560  << nEeZsErrorsType1_ << " in a complete RU.\n";
561 } //end of analyzeEE
562 
564  bool eventError = false;
565  nEbZsErrors_ = 0;
566  nEbZsErrorsType1_ = 0;
567  vector<pair<int, int>> xtalEtaPhi;
568 
569  xtalEtaPhi.reserve(nEbPhi * nEbEta);
570 
574  std::unique_ptr<std::array<std::array<energiesEb_t, nEbPhi>, nEbEta>> ebEnergies =
575  std::make_unique<std::array<std::array<energiesEb_t, nEbPhi>, nEbEta>>();
576 
577  for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
578  for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0) {
579  (*ebEnergies)[iEta0][iPhi0].noZsRecE = -numeric_limits<double>::max();
580  (*ebEnergies)[iEta0][iPhi0].recE = -numeric_limits<double>::max();
581  (*ebEnergies)[iEta0][iPhi0].simE = 0; //must be zero.
582  (*ebEnergies)[iEta0][iPhi0].simHit = 0;
583  (*ebEnergies)[iEta0][iPhi0].gain12 = false;
584  xtalEtaPhi.push_back(pair<int, int>(iEta0, iPhi0));
585  }
586  }
587 
588  // get the barrel geometry:
589  auto geoHandle = es.getHandle(geoToken);
590  const CaloSubdetectorGeometry* geometry_p = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
591  //CaloSubdetectorGeometry const& geometry = *geometry_p;
592 
593  //EB unsuppressed digis:
595  const EBDataFrame& frame = *it;
596  int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(frame.id()).ieta());
597  int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(frame.id()).iphi());
598  if (iEta0 < 0 || iEta0 >= nEbEta) {
599  stringstream s;
600  s << "EcalSelectiveReadoutValidation: "
601  << "iEta0 (= " << iEta0 << ") is out of range ("
602  << "[0," << nEbEta - 1 << "]\n";
603  throw cms::Exception(s.str());
604  }
605  if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
606  stringstream s;
607  s << "EcalSelectiveReadoutValidation: "
608  << "iPhi0 (= " << iPhi0 << ") is out of range ("
609  << "[0," << nEbPhi - 1 << "]\n";
610  throw cms::Exception(s.str());
611  }
612 
613  (*ebEnergies)[iEta0][iPhi0].noZsRecE = frame2Energy(frame);
614  (*ebEnergies)[iEta0][iPhi0].gain12 = true;
615  for (int i = 0; i < frame.size(); ++i) {
616  const int gain12Code = 0x1;
617  if (frame[i].gainId() != gain12Code)
618  (*ebEnergies)[iEta0][iPhi0].gain12 = false;
619  }
620 
621  const GlobalPoint xtalPos = geometry_p->getGeometry(frame.id())->getPosition();
622 
623  (*ebEnergies)[iEta0][iPhi0].phi = rad2deg * ((double)xtalPos.phi());
624  (*ebEnergies)[iEta0][iPhi0].eta = xtalPos.eta();
625  } //next non-zs digi
626 
627  //EB sim hits
628  for (vector<PCaloHit>::const_iterator it = ebSimHits_->begin(); it != ebSimHits_->end(); ++it) {
629  const PCaloHit& simHit = *it;
630  EBDetId detId(simHit.id());
631  int iEta = detId.ieta();
632  int iEta0 = iEta2cIndex(iEta);
633  int iPhi = detId.iphi();
634  int iPhi0 = iPhi2cIndex(iPhi);
635  (*ebEnergies)[iEta0][iPhi0].simE += simHit.energy();
636  ++(*ebEnergies)[iEta0][iPhi0].simHit;
637  }
638 
639  bool crystalShot[nEbEta][nEbPhi];
640  pair<int, int> EBxtalCoor[nEbEta][nEbPhi];
641 
642  for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
643  for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0) {
644  crystalShot[iEta0][iPhi0] = false;
645  EBxtalCoor[iEta0][iPhi0] = make_pair(0, 0);
646  }
647  }
648 
650  const EBDataFrame& frame = *it;
651  int iEta = static_cast<const EBDetId&>(frame.id()).ieta();
652  int iPhi = static_cast<const EBDetId&>(frame.id()).iphi();
653  int iEta0 = iEta2cIndex(iEta);
654  int iPhi0 = iPhi2cIndex(iPhi);
655  if (iEta0 < 0 || iEta0 >= nEbEta) {
656  throw(cms::Exception("EcalSelectiveReadoutValidation") << "iEta0 (= " << iEta0 << ") is out of range ("
657  << "[0," << nEbEta - 1 << "]");
658  }
659  if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
660  throw(cms::Exception("EcalSelectiveReadoutValidation") << "iPhi0 (= " << iPhi0 << ") is out of range ("
661  << "[0," << nEbPhi - 1 << "]");
662  }
663  assert(iEta0 >= 0 && iEta0 < nEbEta);
664  assert(iPhi0 >= 0 && iPhi0 < nEbPhi);
665  if (!crystalShot[iEta0][iPhi0]) {
666  crystalShot[iEta0][iPhi0] = true;
667  EBxtalCoor[iEta0][iPhi0] = make_pair(xtalGraphX(frame.id()), xtalGraphY(frame.id()));
668  } else {
669  cout << "Error: several digi for same crystal!";
670  abort();
671  }
672  if (localReco_) {
673  (*ebEnergies)[iEta0][iPhi0].recE = frame2Energy(frame);
674  }
675 
676  (*ebEnergies)[iEta0][iPhi0].gain12 = true;
677  for (int i = 0; i < frame.size(); ++i) {
678  const int gain12Code = 0x1;
679  if (frame[i].gainId() != gain12Code) {
680  (*ebEnergies)[iEta0][iPhi0].gain12 = false;
681  }
682  }
683 
685 
686  bool highInterest = false;
687 
688  // if(srf == ebSrFlags_->end()){
689  // throw cms::Exception("EcalSelectiveReadoutValidation")
690  // << __FILE__ << ":" << __LINE__ << ": SR flag not found";
691  //}
692 
693  if (srf != ebSrFlags_->end()) {
694  highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK) == EcalSrFlag::SRF_FULL);
695  }
696 
697  if (highInterest) {
699  } else {
700  int v = dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr);
701  fill(meEbLiZsFir_, v);
702  if (v < ebZsThr_) {
703  eventError = true;
704  ++nEbZsErrors_;
705  pair<int, int> ru = dccCh(frame.id());
706  if (isRuComplete_[ru.first - 1][ru.second - 1])
708  if (nEbZsErrors_ < 3) {
709  srApplicationErrorLog_ << event.id() << ", "
710  << "RU " << frame.id() << ", "
711  << "DCC " << ru.first << " Ch : " << ru.second << ": "
712  << "LI channel under ZS threshold.\n";
713  }
714  if (nEbZsErrors_ == 3) {
715  srApplicationErrorLog_ << event.id() << ": "
716  << "more ZS errors for this event...\n";
717  }
718  }
719  }
720  } //next EB digi
721 
722  for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
723  for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0)
724  fill(meChOcc_,
725  EBxtalCoor[iEta0][iPhi0].first,
726  EBxtalCoor[iEta0][iPhi0].second,
727  crystalShot[iEta0][iPhi0] ? 1. : 0.);
728  }
729 
730  if (!localReco_) {
731  for (RecHitCollection::const_iterator it = ebRecHits_->begin(); it != ebRecHits_->end(); ++it) {
732  const RecHit& hit = *it;
733  int iEta = static_cast<const EBDetId&>(hit.id()).ieta();
734  int iPhi = static_cast<const EBDetId&>(hit.id()).iphi();
735  int iEta0 = iEta2cIndex(iEta);
736  int iPhi0 = iPhi2cIndex(iPhi);
737  if (iEta0 < 0 || iEta0 >= nEbEta) {
738  LogError("EcalSrValid") << "iEta0 (= " << iEta0 << ") is out of range ("
739  << "[0," << nEbEta - 1 << "]\n";
740  }
741  if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
742  LogError("EcalSrValid") << "iPhi0 (= " << iPhi0 << ") is out of range ("
743  << "[0," << nEbPhi - 1 << "]\n";
744  }
745  (*ebEnergies)[iEta0][iPhi0].recE = hit.energy();
746  }
747  }
748 
749  for (unsigned int i = 0; i < xtalEtaPhi.size(); ++i) {
750  int iEta0 = xtalEtaPhi[i].first;
751  int iPhi0 = xtalEtaPhi[i].second;
752  energiesEb_t& energies = (*ebEnergies)[iEta0][iPhi0];
753 
754  double recE = energies.recE;
755  if (recE != -numeric_limits<double>::max()) { //not zero suppressed
756  fill(meEbRecE_, (*ebEnergies)[iEta0][iPhi0].recE);
757  fill(meEbEMean_, ievt_ + 1, recE);
758  } //not zero suppressed
759 
760  if (withEbSimHit_) {
761  if (!energies.simHit) { //noise only crystal channel
762  fill(meEbNoise_, energies.noZsRecE);
763  } else {
764  fill(meEbSimE_, energies.simE);
765  fill(meEbRecEHitXtal_, energies.recE);
766  }
767  fill(meEbRecVsSimE_, energies.simE, energies.recE);
768  fill(meEbNoZsRecVsSimE_, energies.simE, energies.noZsRecE);
769  }
770  }
771 
772  int EBZs1RuCount[2][17][72];
773  int EBFullRuCount[2][17][72];
774  int EBForcedRuCount[2][17][72];
775  std::pair<int, int> EBtowerCoor[2][17][72];
776  for (int iZ(0); iZ < 2; iZ++) {
777  for (int iEta(0); iEta < 17; iEta++) {
778  for (int iPhi(0); iPhi < 72; iPhi++) {
779  EBZs1RuCount[iZ][iEta][iPhi] = 0;
780  EBFullRuCount[iZ][iEta][iPhi] = 0;
781  EBForcedRuCount[iZ][iEta][iPhi] = 0;
782  }
783  }
784  }
785 
786  //SRF
787  nEbFROCnt_ = 0;
788  char ebSrfMark[2][17][72];
789  bzero(ebSrfMark, sizeof(ebSrfMark));
790  // int idbg = 0;
791  for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
792  const EBSrFlag& srf = *it;
793  int iEtaAbs = srf.id().ietaAbs();
794  int iPhi = srf.id().iphi();
795  int iZ = srf.id().zside();
796 
797  // cout << "--> " << ++idbg << iEtaAbs << " " << iPhi << " " << iZ
798  // << " " << srf.id() << "\n";
799 
800  if (iEtaAbs < 1 || iEtaAbs > 17 || iPhi < 1 || iPhi > 72)
801  throw cms::Exception("EcalSelectiveReadoutValidation")
802  << "Found a barrel SRF with an invalid det ID: " << srf.id() << ".\n";
803  ++ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1];
804  if (ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] > 1)
805  throw cms::Exception("EcalSelectiveReadoutValidation") << "Duplicate SRF for RU " << srf.id() << ".\n";
806 
807  EBtowerCoor[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] = std::pair<int, int>(srf.id().ieta(), srf.id().iphi());
808 
809  int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
810  if (flag == EcalSrFlag::SRF_ZS1) {
811  EBZs1RuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
812  }
813  if (flag == EcalSrFlag::SRF_FULL) {
814  EBFullRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
815  ++nEbFROCnt_;
816  }
817  if (srf.value() & EcalSrFlag::SRF_FORCED_MASK) {
818  EBForcedRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
819  }
820  }
821  for (int iZ(0); iZ < 2; iZ++) {
822  for (int iEta(0); iEta < 17; iEta++) {
823  for (int iPhi(0); iPhi < 72; iPhi++) {
824  float x(EBtowerCoor[iZ][iEta][iPhi].first);
825  float y(EBtowerCoor[iZ][iEta][iPhi].second);
826  fill(meZs1Ru_, x, y, EBZs1RuCount[iZ][iEta][iPhi]);
827  fill(meFullRoRu_, x, y, EBFullRuCount[iZ][iEta][iPhi]);
828  fill(meForcedRu_, x, y, EBForcedRuCount[iZ][iEta][iPhi]);
829  }
830  }
831  }
832 
833  if (eventError)
834  srApplicationErrorLog_ << event.id() << ": " << nEbZsErrors_
835  << " ZS-flagged EB channels under "
836  "the ZS threshold, whose "
837  << nEbZsErrorsType1_ << " in a complete RU.\n";
838 }
839 
841 
843  // endcap mapping
845 
846  //electronics map
848 
849  initAsciiFile();
850 }
851 
854 }
855 
857  edm::Run const&,
858  edm::EventSetup const&) {
859  ibooker.setCurrentFolder("EcalDigisV/SelectiveReadout");
860 
861  {
862  auto scope = DQMStore::IBooker::UseRunScope(ibooker);
863  meL1aRate_ = bookFloat(ibooker, "l1aRate_");
864  }
865 
866  meDccVol_ = bookProfile(ibooker,
867  "hDccVol", //"EcalDccEventSizeComputed",
868  "ECAL DCC event fragment size;Dcc id; "
869  "<Event size> (kB)",
870  nDccs_,
871  .5,
872  .5 + nDccs_);
873 
874  meDccLiVol_ = bookProfile(ibooker,
875  "hDccLiVol",
876  "LI channel payload per DCC;Dcc id; "
877  "<Event size> (kB)",
878  nDccs_,
879  .5,
880  .5 + nDccs_);
881 
882  meDccHiVol_ = bookProfile(ibooker,
883  "hDccHiVol",
884  "HI channel payload per DCC;Dcc id; "
885  "<Event size> (kB)",
886  nDccs_,
887  .5,
888  .5 + nDccs_);
889 
890  meDccVolFromData_ = bookProfile(ibooker,
891  "hDccVolFromData", //"EcalDccEventSize",
892  "ECAL DCC event fragment size;Dcc id; "
893  "<Event size> (kB)",
894  nDccs_,
895  .5,
896  .5 + nDccs_);
897 
898  meVolBLI_ = book1D(ibooker,
899  "hVolBLI", // "EBLowInterestPayload",
900  "ECAL Barrel low interest crystal data payload;"
901  "Event size (kB);Nevts",
902  100,
903  0.,
904  200.);
905 
906  meVolELI_ = book1D(ibooker,
907  "hVolELI", //"EELowInterestPayload",
908  "Endcap low interest crystal data payload;"
909  "Event size (kB);Nevts",
910  100,
911  0.,
912  200.);
913 
914  meVolLI_ = book1D(ibooker,
915  "hVolLI", //"EcalLowInterestPayload",
916  "ECAL low interest crystal data payload;"
917  "Event size (kB);Nevts",
918  100,
919  0.,
920  200.);
921 
922  meVolBHI_ = book1D(ibooker,
923  "hVolBHI", //"EBHighInterestPayload",
924  "Barrel high interest crystal data payload;"
925  "Event size (kB);Nevts",
926  100,
927  0.,
928  200.);
929 
930  meVolEHI_ = book1D(ibooker,
931  "hVolEHI", //"EEHighInterestPayload",
932  "Endcap high interest crystal data payload;"
933  "Event size (kB);Nevts",
934  100,
935  0.,
936  200.);
937 
938  meVolHI_ = book1D(ibooker,
939  "hVolHI", //"EcalHighInterestPayload",
940  "ECAL high interest crystal data payload;"
941  "Event size (kB);Nevts",
942  100,
943  0.,
944  200.);
945 
946  meVolB_ = book1D(ibooker,
947  "hVolB", //"EBEventSize",
948  "Barrel data volume;Event size (kB);Nevts",
949  100,
950  0.,
951  200.);
952 
953  meVolE_ = book1D(ibooker,
954  "hVolE", //"EEEventSize",
955  "Endcap data volume;Event size (kB);Nevts",
956  100,
957  0.,
958  200.);
959 
960  meVol_ = book1D(ibooker,
961  "hVol", //"EcalEventSize",
962  "ECAL data volume;Event size (kB);Nevts",
963  100,
964  0.,
965  200.);
966 
967  meChOcc_ = bookProfile2D(ibooker,
968  "h2ChOcc", //"EcalChannelOccupancy",
969  "ECAL crystal channel occupancy after zero suppression;"
970  "iX -200 / iEta / iX + 100;"
971  "iY / iPhi (starting from -10^{o}!);"
972  "Event count rate",
973  401,
974  -200.5,
975  200.5,
976  360,
977  .5,
978  360.5);
979 
980  //TP
981  string tpUnit;
982  if (tpInGeV_)
983  tpUnit = string("GeV");
984  else
985  tpUnit = string("TP hw unit");
986  string title;
987  title = string("Trigger primitive TT E_{T};E_{T} ") + tpUnit + string(";Event Count");
988  meTp_ = bookProfile(ibooker,
989  "hTp", //"EcalTriggerPrimitiveEt",
990  title,
991  (tpInGeV_ ? 100 : 40),
992  0.,
993  (tpInGeV_ ? 10. : 40.));
994 
995  meTtf_ = bookProfile(ibooker,
996  "hTtf", //"EcalTriggerTowerFlag",
997  "Trigger primitive TT flag;Flag number;Event count",
998  8,
999  -.5,
1000  7.5);
1001 
1002  title = string("Trigger tower flag vs TP;E_{T}(TT) (") + tpUnit + string(");Flag number");
1003  meTtfVsTp_ = book2D(ibooker, "h2TtfVsTp", title, 100, 0., (tpInGeV_ ? 10. : 40.), 8, -.5, 7.5);
1004 
1005  meTtfVsEtSum_ = book2D(ibooker,
1006  "h2TtfVsEtSum",
1007  "Trigger tower flag vs #sumE_{T};"
1008  "E_{T}(TT) (GeV);"
1009  "TTF",
1010  100,
1011  0.,
1012  10.,
1013  8,
1014  -.5,
1015  7.5);
1016  title = string(
1017  "Trigger primitive Et (TP) vs #sumE_{T};"
1018  "E_{T} (sum) (GeV);"
1019  "E_{T} (TP) (") +
1020  tpUnit + string(")");
1021 
1022  meTpVsEtSum_ = book2D(ibooker, "h2TpVsEtSum", title, 100, 0., 10., 100, 0., (tpInGeV_ ? 10. : 40.));
1023 
1024  title = string(
1025  "Trigger primitive E_{T};"
1026  "iEta;"
1027  "iPhi;"
1028  "E_{T} (TP) (") +
1029  tpUnit + string(")");
1030  meTpMap_ = bookProfile2D(ibooker, "h2Tp", title, 57, -28.5, 28.5, 72, .5, 72.5);
1031 
1032  //SRF
1033  meFullRoRu_ = book2D(ibooker,
1034  "h2FRORu", //"EcalFullReadoutSRFlagMap",
1035  "Full Read-out readout unit;"
1036  "iX - 40 / iEta / iX + 20;"
1037  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1038  "Event count",
1039  80,
1040  -39.5,
1041  40.5,
1042  72,
1043  .5,
1044  72.5);
1045 
1046  meFullRoCnt_ = book1D(ibooker,
1047  "hFROCnt",
1048  "Number of Full-readout-flagged readout units;"
1049  "FRO RU count;Event count",
1050  300,
1051  -.5,
1052  299.5);
1053 
1054  meEbFullRoCnt_ = book1D(ibooker,
1055  "hEbFROCnt",
1056  "Number of EB Full-readout-flagged readout units;"
1057  "FRO RU count;Event count",
1058  200,
1059  -.5,
1060  199.5);
1061 
1062  meEeFullRoCnt_ = book1D(ibooker,
1063  "hEeFROCnt",
1064  "Number of EE Full-readout-flagged readout units;"
1065  "FRO RU count;Event count",
1066  200,
1067  -.5,
1068  199.5);
1069 
1070  meZs1Ru_ = book2D(ibooker,
1071  "h2Zs1Ru", //"EbZeroSupp1SRFlagMap",
1072  "Readout unit with ZS-thr-1 flag;"
1073  "iX - 40 / iEta / iX + 20;"
1074  "iY0 / iPhi0 (iPhi = 1 at phi = 0 rad);"
1075  "Event count",
1076  80,
1077  -39.5,
1078  40.5,
1079  72,
1080  .5,
1081  72.5);
1082 
1083  meForcedRu_ = book2D(ibooker,
1084  "h2ForcedRu", //"EcalReadoutUnitForcedBitMap",
1085  "ECAL readout unit with forced bit of SR flag on;"
1086  "iX - 40 / iEta / iX + 20;"
1087  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1088  "Event count",
1089  80,
1090  -39.5,
1091  40.5,
1092  72,
1093  .5,
1094  72.5);
1095 
1096  meLiTtf_ = bookProfile2D(ibooker,
1097  "h2LiTtf", //"EcalLowInterestTriggerTowerFlagMap",
1098  "Low interest trigger tower flags;"
1099  "iEta;"
1100  "iPhi;"
1101  "Event count",
1102  57,
1103  -28.5,
1104  28.5,
1105  72,
1106  .5,
1107  72.5);
1108 
1109  meMiTtf_ = bookProfile2D(ibooker,
1110  "h2MiTtf", //"EcalMidInterestTriggerTowerFlagMap",
1111  "Mid interest trigger tower flags;"
1112  "iEta;"
1113  "iPhi;"
1114  "Event count",
1115  57,
1116  -28.5,
1117  28.5,
1118  72,
1119  .5,
1120  72.5);
1121 
1122  meHiTtf_ = bookProfile2D(ibooker,
1123  "h2HiTtf", //"EcalHighInterestTriggerTowerFlagMap",
1124  "High interest trigger tower flags;"
1125  "iEta;"
1126  "iPhi;"
1127  "Event count",
1128  57,
1129  -28.5,
1130  28.5,
1131  72,
1132  .5,
1133  72.5);
1134 
1135  meForcedTtf_ = book2D(ibooker,
1136  "h2ForcedTtf", //"EcalTtfForcedBitMap",
1137  "Trigger tower flags with forced bit set;"
1138  "iEta;"
1139  "iPhi;"
1140  "Event count",
1141  57,
1142  -28.5,
1143  28.5,
1144  72,
1145  .5,
1146  72.5);
1147 
1148  const float ebMinNoise = -1.;
1149  const float ebMaxNoise = 1.;
1150 
1151  const float eeMinNoise = -1.;
1152  const float eeMaxNoise = 1.;
1153 
1154  const float ebMinE = ebMinNoise;
1155  const float ebMaxE = ebMaxNoise;
1156 
1157  const float eeMinE = eeMinNoise;
1158  const float eeMaxE = eeMaxNoise;
1159 
1160  const int evtMax = 500;
1161 
1162  meEbRecE_ = book1D(ibooker, "hEbRecE", "Crystal reconstructed energy;E (GeV);Event count", 100, ebMinE, ebMaxE);
1163 
1164  meEbEMean_ = bookProfile(ibooker, "hEbEMean", "EE <E_hit>;event #;<E_hit> (GeV)", evtMax, .5, evtMax + .5);
1165 
1166  meEbNoise_ = book1D(ibooker,
1167  "hEbNoise",
1168  "Crystal noise "
1169  "(rec E of crystal without deposited energy)"
1170  ";Rec E (GeV);Event count",
1171  100,
1172  ebMinNoise,
1173  ebMaxNoise);
1174 
1175  meEbLiZsFir_ = book1D(ibooker,
1176  "zsEbLiFIRemu",
1177  "Emulated ouput of ZS FIR filter for EB "
1178  "low interest crystals;"
1179  "ADC count*4;"
1180  "Event count",
1181  60,
1182  -30,
1183  30);
1184 
1185  meEbHiZsFir_ = book1D(ibooker,
1186  "zsEbHiFIRemu",
1187  "Emulated ouput of ZS FIR filter for EB "
1188  "high interest crystals;"
1189  "ADC count*4;"
1190  "Event count",
1191  60,
1192  -30,
1193  30);
1194 
1195  //TODO: Fill this histogram...
1196  // meEbIncompleteRUZsFir_ = book1D(ibooker, "zsEbIncompleteRUFIRemu",
1197  // "Emulated ouput of ZS FIR filter for EB "
1198  // "incomplete FRO-flagged RU;"
1199  // "ADC count*4;"
1200  // "Event count",
1201  // 60, -30, 30);
1202 
1203  meEbSimE_ = book1D(ibooker, "hEbSimE", "EB hit crystal simulated energy", 100, ebMinE, ebMaxE);
1204 
1205  meEbRecEHitXtal_ = book1D(ibooker, "hEbRecEHitXtal", "EB rec energy of hit crystals", 100, ebMinE, ebMaxE);
1206 
1207  meEbRecVsSimE_ = book2D(ibooker,
1208  "hEbRecVsSimE",
1209  "Crystal simulated vs reconstructed energy;"
1210  "Esim (GeV);Erec GeV);Event count",
1211  100,
1212  ebMinE,
1213  ebMaxE,
1214  100,
1215  ebMinE,
1216  ebMaxE);
1217 
1218  meEbNoZsRecVsSimE_ = book2D(ibooker,
1219  "hEbNoZsRecVsSimE",
1220  "Crystal no-zs simulated vs reconstructed "
1221  "energy;"
1222  "Esim (GeV);Erec GeV);Event count",
1223  100,
1224  ebMinE,
1225  ebMaxE,
1226  100,
1227  ebMinE,
1228  ebMaxE);
1229 
1230  meEeRecE_ = book1D(ibooker,
1231  "hEeRecE",
1232  "EE crystal reconstructed energy;E (GeV);"
1233  "Event count",
1234  100,
1235  eeMinE,
1236  eeMaxE);
1237 
1238  meEeEMean_ = bookProfile(ibooker, "hEeEMean", "<E_{EE hit}>;event;<E_{hit}> (GeV)", evtMax, .5, evtMax + .5);
1239 
1240  meEeNoise_ = book1D(ibooker,
1241  "hEeNoise",
1242  "EE crystal noise "
1243  "(rec E of crystal without deposited energy);"
1244  "E (GeV);Event count",
1245  200,
1246  eeMinNoise,
1247  eeMaxNoise);
1248 
1249  meEeLiZsFir_ = book1D(ibooker,
1250  "zsEeLiFIRemu",
1251  "Emulated ouput of ZS FIR filter for EE "
1252  "low interest crystals;"
1253  "ADC count*4;"
1254  "Event count",
1255  60,
1256  -30,
1257  30);
1258 
1259  meEeHiZsFir_ = book1D(ibooker,
1260  "zsEeHiFIRemu",
1261  "Emulated ouput of ZS FIR filter for EE "
1262  "high interest crystals;"
1263  "ADC count*4;"
1264  "Event count",
1265  60,
1266  -30,
1267  30);
1268 
1269  //TODO: Fill this histogram...
1270  // meEeIncompleteRUZsFir_ = book1D(ibooker, "zsEeIncompleteRUFIRemu",
1271  // "Emulated ouput of ZS FIR filter for EE "
1272  // "incomplete FRO-flagged RU;"
1273  // "ADC count*4;"
1274  // "Event count",
1275  // 60, -30, 30);
1276 
1277  meEeSimE_ = book1D(ibooker, "hEeSimE", "EE hit crystal simulated energy", 100, eeMinE, eeMaxE);
1278 
1279  meEeRecEHitXtal_ = book1D(ibooker, "hEeRecEHitXtal", "EE rec energy of hit crystals", 100, eeMinE, eeMaxE);
1280 
1281  meEeRecVsSimE_ = book2D(ibooker,
1282  "hEeRecVsSimE",
1283  "EE crystal simulated vs reconstructed energy;"
1284  "Esim (GeV);Erec GeV);Event count",
1285  100,
1286  eeMinE,
1287  eeMaxE,
1288  100,
1289  eeMinE,
1290  eeMaxE);
1291 
1292  meEeNoZsRecVsSimE_ = book2D(ibooker,
1293  "hEeNoZsRecVsSimE",
1294  "EE crystal no-zs simulated vs "
1295  "reconstructed "
1296  "energy;Esim (GeV);Erec GeV);Event count",
1297  100,
1298  eeMinE,
1299  eeMaxE,
1300  100,
1301  eeMinE,
1302  eeMaxE);
1303 
1304  meSRFlagsConsistency_ = book2D(ibooker,
1305  "hSRAlgoErrorMap",
1306  "TTFlags and SR Flags mismatch;"
1307  "iX - 40 / iEta / iX + 20;"
1308  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1309  "Event count",
1310  80,
1311  -39.5,
1312  40.5,
1313  72,
1314  .5,
1315  72.5);
1316 
1317  //Readout Units histos (interest/Ncrystals)
1318  meIncompleteFROMap_ = book2D(ibooker,
1319  "hIncompleteFROMap",
1320  "Incomplete full-readout-flagged readout units;"
1321  "iX - 40 / iEta / iX + 20;"
1322  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1323  "Event count",
1324  80,
1325  -39.5,
1326  40.5,
1327  72,
1328  .5,
1329  72.5);
1330 
1331  meIncompleteFROCnt_ = book1D(ibooker,
1332  "hIncompleteFROCnt",
1333  "Number of incomplete full-readout-flagged "
1334  "readout units;"
1335  "Number of RUs;Event count;",
1336  200,
1337  -.5,
1338  199.5);
1339 
1341  "hIncompleteFRORateMap",
1342  "Incomplete full-readout-flagged readout units;"
1343  "iX - 40 / iEta / iX + 20;"
1344  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1345  "Incomplete error rate",
1346  80,
1347  -39.5,
1348  40.5,
1349  72,
1350  .5,
1351  72.5);
1352 
1353  meDroppedFROMap_ = book2D(ibooker,
1354  "hDroppedFROMap",
1355  "Dropped full-readout-flagged readout units;"
1356  "iX - 40 / iEta / iX + 20;"
1357  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1358  "Event count",
1359  80,
1360  -39.5,
1361  40.5,
1362  72,
1363  .5,
1364  72.5);
1365 
1366  meDroppedFROCnt_ = book1D(ibooker,
1367  "hDroppedFROCnt",
1368  "Number of dropped full-readout-flagged "
1369  "RU count;RU count;Event count",
1370  200,
1371  -.5,
1372  199.5);
1373 
1374  meCompleteZSCnt_ = book1D(ibooker,
1375  "hCompleteZsCnt",
1376  "Number of zero-suppressed-flagged RU fully "
1377  "readout;"
1378  "RU count;Event count",
1379  200,
1380  -.5,
1381  199.5);
1382 
1383  stringstream buf;
1384  buf << "Number of LI EB channels below the " << ebZsThr_ / 4.
1385  << " ADC count ZS threshold;"
1386  "Channel count;Event count",
1387  meEbZsErrCnt_ = book1D(ibooker, "hEbZsErrCnt", buf.str(), 200, -.5, 199.5);
1388 
1389  buf.str("");
1390  buf << "Number of LI EE channels below the " << eeZsThr_ / 4.
1391  << " ADC count ZS theshold;"
1392  "Channel count;Event count",
1393  meEeZsErrCnt_ = book1D(ibooker, "hEeZsErrCnt", buf.str(), 200, -.5, 199.5);
1394 
1395  meZsErrCnt_ = book1D(ibooker,
1396  "hZsErrCnt",
1397  "Number of LI channels below the ZS threshold;"
1398  "Channel count;Event count",
1399  200,
1400  -.5,
1401  199.5);
1402 
1403  meEbZsErrType1Cnt_ = book1D(ibooker,
1404  "hEbZsErrType1Cnt",
1405  "Number of EB channels below the ZS "
1406  "threshold in a LI but fully readout RU;"
1407  "Channel count;Event count;",
1408  200,
1409  -.5,
1410  199.5);
1411 
1412  meEeZsErrType1Cnt_ = book1D(ibooker,
1413  "hEeZsErrType1Cnt",
1414  "Number EE channels below the ZS threshold"
1415  " in a LI but fully readout RU;"
1416  "Channel count;Event count",
1417  200,
1418  -.5,
1419  199.5);
1420 
1421  meZsErrType1Cnt_ = book1D(ibooker,
1422  "hZsErrType1Cnt",
1423  "Number of LI channels below the ZS threshold "
1424  "in a LI but fully readout RU;"
1425  "Channel count;Event count",
1426  200,
1427  -.5,
1428  199.5);
1429 
1431  "hDroppedFRORateMap",
1432  "Dropped full-readout-flagged readout units"
1433  "iX - 40 / iEta / iX + 20;"
1434  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1435  "Dropping rate",
1436  80,
1437  -39.5,
1438  40.5,
1439  72,
1440  .5,
1441  72.5);
1442 
1443  meCompleteZSMap_ = book2D(ibooker,
1444  "hCompleteZSMap",
1445  "Complete zero-suppressed-flagged readout units;"
1446  "iX - 40 / iEta / iX + 20;"
1447  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1448  "Event count",
1449  80,
1450  -39.5,
1451  40.5,
1452  72,
1453  .5,
1454  72.5);
1455 
1457  "hCompleteZSRate",
1458  "Complete zero-suppressed-flagged readout units;"
1459  "iX - 40 / iEta / iX + 20;"
1460  "iY / iPhi (iPhi = 1 at phi = 0 rad);"
1461  "Completeness rate",
1462  80,
1463  -39.5,
1464  40.5,
1465  72,
1466  .5,
1467  72.5);
1468 
1469  //print list of available histograms (must be called after
1470  //the bookXX methods):
1472 
1473  //check the histList parameter:
1474  stringstream s;
1475  for (set<string>::iterator it = histList_.begin(); it != histList_.end(); ++it) {
1476  if (*it != string("all") && availableHistList_.find(*it) == availableHistList_.end()) {
1477  s << (s.str().empty() ? "" : ", ") << *it;
1478  }
1479  }
1480  if (!s.str().empty()) {
1481  LogWarning("Configuration") << "Parameter 'histList' contains some unknown histogram(s). "
1482  "Check spelling. Following name were not found: "
1483  << s.str();
1484  }
1485 }
1486 
1488  int TTFlagCount[8];
1489  int LiTTFlagCount[nTtEta][nTtPhi];
1490  int MiTTFlagCount[nTtEta][nTtPhi];
1491  int HiTTFlagCount[nTtEta][nTtPhi];
1492  for (int iTTFlag(0); iTTFlag < 8; iTTFlag++) {
1493  TTFlagCount[iTTFlag] = 0;
1494  }
1495  for (int iTtEta(0); iTtEta < nTtEta; iTtEta++) {
1496  for (int iTtPhi(0); iTtPhi < nTtPhi; iTtPhi++) {
1497  LiTTFlagCount[iTtEta][iTtPhi] = 0;
1498  MiTTFlagCount[iTtEta][iTtPhi] = 0;
1499  HiTTFlagCount[iTtEta][iTtPhi] = 0;
1500  }
1501  }
1502  int tpEtCount[100];
1503  for (int iEt(0); iEt < 100; iEt++) {
1504  tpEtCount[iEt] = 0;
1505  }
1506 
1507  const EcalTPGPhysicsConstMap& physMap = es.getData(physHandle).getMap();
1508 
1509  const EcalTPGGroups::EcalTPGGroupsMap& lutGrpMap = es.getData(lutGrpHandle).getMap();
1510 
1511  const EcalTPGLutIdMap::EcalTPGLutMap& lutMap = es.getData(lutMapHandle).getMap();
1512 
1514  double lsb10bitsEB(ebItr == physMap.end() ? 0. : ebItr->second.EtSat / 1024.);
1516  double lsb10bitsEE(eeItr == physMap.end() ? 0. : eeItr->second.EtSat / 1024.);
1517 
1518  for (EcalTrigPrimDigiCollection::const_iterator it = tps_->begin(); it != tps_->end(); ++it) {
1519  double tpEt;
1520  if (tpInGeV_) {
1521  EcalTrigTowerDetId const& towerId(it->id());
1522  unsigned int ADC = it->compressedEt();
1523 
1524  double lsb10bits(0.);
1525  if (towerId.subDet() == EcalBarrel)
1526  lsb10bits = lsb10bitsEB;
1527  else if (towerId.subDet() == EcalEndcap)
1528  lsb10bits = lsb10bitsEE;
1529 
1530  int tpg10bits = 0;
1531  EcalTPGGroups::EcalTPGGroupsMapItr itgrp = lutGrpMap.find(towerId.rawId());
1532  uint32_t lutGrp = 999;
1533  if (itgrp != lutGrpMap.end())
1534  lutGrp = itgrp->second;
1535 
1536  EcalTPGLutIdMap::EcalTPGLutMapItr itLut = lutMap.find(lutGrp);
1537  if (itLut != lutMap.end()) {
1538  const unsigned int* lut = (itLut->second).getLut();
1539  for (unsigned int i = 0; i < 1024; i++)
1540  if (ADC == (0xff & lut[i])) {
1541  tpg10bits = i;
1542  break;
1543  }
1544  }
1545 
1546  tpEt = lsb10bits * tpg10bits;
1547  } else {
1548  tpEt = it->compressedEt();
1549  }
1550  int iEta = it->id().ieta();
1551  int iEta0 = iTtEta2cIndex(iEta);
1552  int iPhi = it->id().iphi();
1553  int iPhi0 = iTtPhi2cIndex(iPhi);
1554  double etSum = ttEtSums[iEta0][iPhi0];
1555 
1556  int iE = meTp_->getTProfile()->FindFixBin(tpEt);
1557  if ((iE >= 0) && (iE < 100)) {
1558  ++tpEtCount[iE];
1559  } else {
1560  // FindFixBin might return an overflow bin (outside tpEtCount range).
1561  // To prevent a memory overflow / segfault, these values are ignored.
1562  //std::cout << "EcalSelectiveReadoutValidation: Invalid iE value: " << iE << std::endl;
1563  }
1564 
1565  fill(meTpVsEtSum_, etSum, tpEt);
1566  ++TTFlagCount[it->ttFlag()];
1567  if ((it->ttFlag() & 0x3) == 0) {
1568  LiTTFlagCount[iEta0][iPhi0] += 1;
1569  } else if ((it->ttFlag() & 0x3) == 1) {
1570  MiTTFlagCount[iEta0][iPhi0] += 1;
1571  } else if ((it->ttFlag() & 0x3) == 3) {
1572  HiTTFlagCount[iEta0][iPhi0] += 1;
1573  }
1574  if ((it->ttFlag() & 0x4)) {
1575  fill(meForcedTtf_, iEta, iPhi);
1576  }
1577 
1578  fill(meTtfVsTp_, tpEt, it->ttFlag());
1579  fill(meTtfVsEtSum_, etSum, it->ttFlag());
1580  fill(meTpMap_, iEta, iPhi, tpEt, 1.);
1581  }
1582 
1583  for (int ittflag(0); ittflag < 8; ittflag++) {
1584  fill(meTtf_, ittflag, TTFlagCount[ittflag]);
1585  }
1586  for (int iTtEta(0); iTtEta < nTtEta; iTtEta++) {
1587  for (int iTtPhi(0); iTtPhi < nTtPhi; iTtPhi++) {
1588  fill(meLiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), LiTTFlagCount[iTtEta][iTtPhi]);
1589  fill(meMiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), MiTTFlagCount[iTtEta][iTtPhi]);
1590  fill(meHiTtf_, cIndex2iTtEta(iTtEta), cIndex2iTtPhi(iTtPhi), HiTTFlagCount[iTtEta][iTtPhi]);
1591  }
1592  }
1593  if (tpInGeV_) {
1594  for (int iE(0); iE < 100; iE++) {
1595  fill(meTp_, iE, tpEtCount[iE]);
1596  }
1597  } else {
1598  for (int iE(0); iE < 40; iE++) {
1599  fill(meTp_, iE, tpEtCount[iE]);
1600  }
1601  }
1602 }
1603 
1605  anaDigiInit();
1606 
1607  //Complete RU, i.e. RU actually fully readout
1608  for (int iDcc = minDccId_; iDcc <= maxDccId_; ++iDcc) {
1609  for (int iCh = 1; iCh < nDccRus_[iDcc - minDccId_]; ++iCh) {
1610  isRuComplete_[iDcc - minDccId_][iCh - 1] = (nPerRu_[iDcc - minDccId_][iCh - 1] == getCrystalCount(iDcc, iCh));
1611  }
1612  }
1613 
1614  //Barrel
1615  for (unsigned int digis = 0; digis < ebDigis_->size(); ++digis) {
1616  EBDataFrame ebdf = (*ebDigis_)[digis];
1617  anaDigi(ebdf, *ebSrFlags_);
1618  }
1619 
1620  // Endcap
1621  for (unsigned int digis = 0; digis < eeDigis_->size(); ++digis) {
1622  EEDataFrame eedf = (*eeDigis_)[digis];
1623  anaDigi(eedf, *eeSrFlags_);
1624  }
1625 
1626  //histos
1627  for (unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0) {
1628  fill(meDccVol_, iDcc0 + 1, getDccEventSize(iDcc0, nPerDcc_[iDcc0]) / kByte_);
1629  fill(meDccLiVol_, iDcc0 + 1, getDccSrDependentPayload(iDcc0, nLiRuPerDcc_[iDcc0], nLiPerDcc_[iDcc0]) / kByte_);
1630  fill(meDccHiVol_, iDcc0 + 1, getDccSrDependentPayload(iDcc0, nHiRuPerDcc_[iDcc0], nHiPerDcc_[iDcc0]) / kByte_);
1631  const FEDRawDataCollection& raw = *fedRaw_;
1632  fill(meDccVolFromData_, iDcc0 + 1, ((double)raw.FEDData(601 + iDcc0).size()) / kByte_);
1633  }
1634 
1635  //low interesest channels:
1636  double a = nEbLI_ * getBytesPerCrystal() / kByte_; //getEbEventSize(nEbLI_)/kByte_;
1637  fill(meVolBLI_, a);
1638  double b = nEeLI_ * getBytesPerCrystal() / kByte_; //getEeEventSize(nEeLI_)/kByte_;
1639  fill(meVolELI_, b);
1640  fill(meVolLI_, a + b);
1641 
1642  //high interest chanels:
1643  a = nEbHI_ * getBytesPerCrystal() / kByte_; //getEbEventSize(nEbHI_)/kByte_;
1644  fill(meVolBHI_, a);
1645  b = nEeHI_ * getBytesPerCrystal() / kByte_; //getEeEventSize(nEeHI_)/kByte_;
1646  fill(meVolEHI_, b);
1647  fill(meVolHI_, a + b);
1648 
1649  //any-interest channels:
1650  a = getEbEventSize(nEb_) / kByte_;
1651  fill(meVolB_, a);
1652  b = getEeEventSize(nEe_) / kByte_;
1653  fill(meVolE_, b);
1654  fill(meVol_, a + b);
1655 }
1656 
1657 template <class T, class U>
1658 void EcalSelectiveReadoutValidation::anaDigi(const T& frame, const U& srFlagColl) {
1659  const DetId& xtalId = frame.id();
1660  typedef typename U::key_type RuDetId;
1661  const RuDetId& ruId = readOutUnitOf(frame.id());
1662  typename U::const_iterator srf = srFlagColl.find(ruId);
1663 
1664  bool highInterest = false;
1665  int flag = 0;
1666 
1667  if (srf != srFlagColl.end()) {
1668  flag = srf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
1669 
1670  highInterest = (flag == EcalSrFlag::SRF_FULL);
1671  }
1672 
1673  bool barrel = (xtalId.subdetId() == EcalBarrel);
1674 
1675  pair<int, int> ch = dccCh(xtalId);
1676 
1677  if (barrel) {
1678  ++nEb_;
1679  if (highInterest) {
1680  ++nEbHI_;
1681  } else { //low interest
1682  ++nEbLI_;
1683  }
1684  int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(xtalId).ieta());
1685  int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(xtalId).iphi());
1686  if (!ebRuActive_[iEta0 / ebTtEdge][iPhi0 / ebTtEdge]) {
1687  ++nRuPerDcc_[ch.first - minDccId_];
1688  if (highInterest) {
1689  ++nHiRuPerDcc_[ch.first - minDccId_];
1690  } else {
1691  ++nLiRuPerDcc_[ch.first - minDccId_];
1692  }
1693 
1694  ebRuActive_[iEta0 / ebTtEdge][iPhi0 / ebTtEdge] = true;
1695  }
1696  } else { //endcap
1697  ++nEe_;
1698  if (highInterest) {
1699  ++nEeHI_;
1700  } else { //low interest
1701  ++nEeLI_;
1702  }
1703  int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
1704  int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
1705  int iZ0 = static_cast<const EEDetId&>(frame.id()).zside() > 0 ? 1 : 0;
1706 
1707  if (!eeRuActive_[iZ0][iX0 / scEdge][iY0 / scEdge]) {
1708  ++nRuPerDcc_[ch.first - minDccId_];
1709  if (highInterest) {
1710  ++nHiRuPerDcc_[ch.first - minDccId_];
1711  } else {
1712  ++nLiRuPerDcc_[ch.first - minDccId_];
1713  }
1714 
1715  eeRuActive_[iZ0][iX0 / scEdge][iY0 / scEdge] = true;
1716  }
1717  }
1718 
1719  if (ch.second < 1 || ch.second > 68) {
1720  throw cms::Exception("EcalSelectiveReadoutValidation")
1721  << "Error in DCC channel retrieval for crystal with detId " << xtalId.rawId()
1722  << "DCC channel out of allowed range [1..68]\n";
1723  }
1724  ++nPerDcc_[ch.first - minDccId_];
1725  ++nPerRu_[ch.first - minDccId_][ch.second - 1];
1726  if (highInterest) {
1727  ++nHiPerDcc_[ch.first - minDccId_];
1728  } else { //low interest channel
1729  ++nLiPerDcc_[ch.first - minDccId_];
1730  }
1731 }
1732 
1734  nEb_ = 0;
1735  nEe_ = 0;
1736  nEeLI_ = 0;
1737  nEeHI_ = 0;
1738  nEbLI_ = 0;
1739  nEbHI_ = 0;
1740  bzero(nPerDcc_, sizeof(nPerDcc_));
1741  bzero(nLiPerDcc_, sizeof(nLiPerDcc_));
1742  bzero(nHiPerDcc_, sizeof(nHiPerDcc_));
1743  bzero(nRuPerDcc_, sizeof(nRuPerDcc_));
1744  bzero(ebRuActive_, sizeof(ebRuActive_));
1745  bzero(eeRuActive_, sizeof(eeRuActive_));
1746  bzero(nPerRu_, sizeof(nPerRu_));
1747  bzero(nLiRuPerDcc_, sizeof(nLiRuPerDcc_));
1748  bzero(nHiRuPerDcc_, sizeof(nHiRuPerDcc_));
1749 }
1750 
1752  static std::atomic<bool> firstCall{true};
1753  bool expected = true;
1754  if (firstCall.compare_exchange_strong(expected, false)) {
1755  stringstream buf;
1756  buf << "Weights:";
1757  for (unsigned i = 0; i < weights_.size(); ++i) {
1758  buf << "\t" << weights_[i];
1759  }
1760  edm::LogInfo("EcalSrValid") << buf.str() << "\n";
1761  firstCall = false;
1762  }
1763  double adc2GeV = 0.;
1764 
1765  if (typeid(EBDataFrame) == typeid(frame)) { //barrel APD
1766  adc2GeV = .035;
1767  } else if (typeid(EEDataFrame) == typeid(frame)) { //endcap VPT
1768  adc2GeV = 0.06;
1769  } else {
1770  assert(false);
1771  }
1772 
1773  double acc = 0;
1774 
1775  const int n = min(frame.size(), (int)weights_.size());
1776 
1777  double gainInv[] = {12., 1., 6., 12.};
1778 
1779  for (int i = 0; i < n; ++i) {
1780  acc += weights_[i] * frame[i].adc() * gainInv[frame[i].gainId()] * adc2GeV;
1781  }
1782  return acc;
1783 }
1784 
1785 int EcalSelectiveReadoutValidation::getRuCount(int iDcc0) const { return nRuPerDcc_[iDcc0]; }
1786 
1787 pair<int, int> EcalSelectiveReadoutValidation::dccCh(const DetId& detId) const {
1788  if (detId.det() != DetId::Ecal) {
1789  throw cms::Exception("InvalidParameter") << "Wrong type of DetId passed to the "
1790  "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
1791  "An ECAL DetId was expected.\n";
1792  }
1793 
1794  DetId xtalId;
1795  switch (detId.subdetId()) {
1796  case EcalTriggerTower: //Trigger tower
1797  {
1798  const EcalTrigTowerDetId tt = detId;
1799  //pick up one crystal of the trigger tower: they are however all readout by
1800  //the same DCC channel in the barrel.
1801  //Arithmetic is easier on the "c" indices:
1802  const int iTtPhi0 = iTtPhi2cIndex(tt.iphi());
1803  const int iTtEta0 = iTtEta2cIndex(tt.ieta());
1804  const int oneXtalPhi0 = iTtPhi0 * 5;
1805  const int oneXtalEta0 = (iTtEta0 - nOneEeTtEta) * 5;
1806 
1807  xtalId = EBDetId(cIndex2iEta(oneXtalEta0), cIndex2iPhi(oneXtalPhi0));
1808  } break;
1809  case EcalEndcap:
1810  if (detId.rawId() & 0x8000) { //Supercrystal
1812  } else { //EE crystal
1813  xtalId = detId;
1814  }
1815  break;
1816  case EcalBarrel: //EB crystal
1817  xtalId = detId;
1818  break;
1819  default:
1820  throw cms::Exception("InvalidParameter")
1821  << "Wrong type of DetId passed to the method "
1822  "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
1823  "A valid EcalTriggerTower, EcalBarrel or EcalEndcap DetId was expected. "
1824  "detid = "
1825  << xtalId.rawId() << ".\n";
1826  }
1827 
1828  const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
1829 
1830  pair<int, int> result;
1831  result.first = EcalElecId.dccId();
1832 
1833  if (result.first < minDccId_ || result.second > maxDccId_) {
1834  throw cms::Exception("OutOfRange") << "Got an invalid DCC ID, DCCID = " << result.first << " for DetId 0x" << hex
1835  << detId.rawId() << " and 0x" << xtalId.rawId() << dec << "\n";
1836  }
1837 
1838  result.second = EcalElecId.towerId();
1839 
1840  if (result.second < 1 || result.second > 68) {
1841  throw cms::Exception("OutOfRange") << "Got an invalid DCC channel ID, DCC_CH = " << result.second << " for DetId 0x"
1842  << hex << detId.rawId() << " and 0x" << xtalId.rawId() << dec << "\n";
1843  }
1844 
1845  return result;
1846 }
1847 
1849  return triggerTowerMap_->towerOf(xtalId);
1850 }
1851 
1853  const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
1854  int iDCC = EcalElecId.dccId();
1855  int iDccChan = EcalElecId.towerId();
1856  const bool ignoreSingle = true;
1857  const vector<EcalScDetId> id = elecMap_->getEcalScDetId(iDCC, iDccChan, ignoreSingle);
1858  return !id.empty() ? id[0] : EcalScDetId();
1859 }
1860 
1862  const EBDigiCollection& ebDigis,
1863  const EEDigiCollection& eeDigis) {
1864  //ecal geometry:
1865  const CaloSubdetectorGeometry* eeGeometry = nullptr;
1866  const CaloSubdetectorGeometry* ebGeometry = nullptr;
1867  if (eeGeometry == nullptr || ebGeometry == nullptr) {
1868  es.getData(geoToken);
1869  auto geoHandle = es.getHandle(geoToken);
1870  eeGeometry = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
1871  ebGeometry = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
1872  }
1873 
1874  //init etSum array:
1875  for (int iEta0 = 0; iEta0 < nTtEta; ++iEta0) {
1876  for (int iPhi0 = 0; iPhi0 < nTtPhi; ++iPhi0) {
1877  ttEtSums[iEta0][iPhi0] = 0.;
1878  }
1879  }
1880 
1882  const EBDataFrame& frame = *it;
1884 
1885  const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
1886  const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
1887  double theta = ebGeometry->getGeometry(frame.id())->getPosition().theta();
1888  double e = frame2EnergyForTp(frame);
1889  if ((frame2EnergyForTp(frame, -1) < e) && (frame2EnergyForTp(frame, 1) < e)) {
1890  ttEtSums[iTtEta0][iTtPhi0] += e * sin(theta);
1891  }
1892  }
1893 
1894  for (EEDigiCollection::const_iterator it = eeDigis.begin(); it != eeDigis.end(); ++it) {
1895  const EEDataFrame& frame = *it;
1897  const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
1898  const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
1899 
1900  double theta = eeGeometry->getGeometry(frame.id())->getPosition().theta();
1901  double e = frame2EnergyForTp(frame);
1902  if ((frame2EnergyForTp(frame, -1) < e) && (frame2EnergyForTp(frame, 1) < e)) {
1903  ttEtSums[iTtEta0][iTtPhi0] += e * sin(theta);
1904  }
1905  }
1906 
1907  //dealing with pseudo-TT in two inner EE eta-ring:
1908  int innerTTEtas[] = {0, 1, 54, 55};
1909  for (unsigned iRing = 0; iRing < sizeof(innerTTEtas) / sizeof(innerTTEtas[0]); ++iRing) {
1910  int iTtEta0 = innerTTEtas[iRing];
1911  //this detector eta-section is divided in only 36 phi bins
1912  //For this eta regions,
1913  //current tower eta numbering scheme is inconsistent. For geometry
1914  //version 133:
1915  //- TT are numbered from 0 to 72 for 36 bins
1916  //- some TT have an even index, some an odd index
1917  //For geometry version 125, there are 72 phi bins.
1918  //The code below should handle both geometry definition.
1919  //If there are 72 input trigger primitives for each inner eta-ring,
1920  //then the average of the trigger primitive of the two pseudo-TT of
1921  //a pair (nEta, nEta+1) is taken as Et of both pseudo TTs.
1922  //If there are only 36 input TTs for each inner eta ring, then half
1923  //of the present primitive of a pseudo TT pair is used as Et of both
1924  //pseudo TTs.
1925 
1926  for (unsigned iTtPhi0 = 0; iTtPhi0 < nTtPhi - 1; iTtPhi0 += 2) {
1927  double et = .5 * (ttEtSums[iTtEta0][iTtPhi0] + ttEtSums[iTtEta0][iTtPhi0 + 1]);
1928  //divides the TT into 2 phi bins in order to match with 72 phi-bins SRP
1929  //scheme or average the Et on the two pseudo TTs if the TT is already
1930  //divided into two trigger primitives.
1931  ttEtSums[iTtEta0][iTtPhi0] = et;
1932  ttEtSums[iTtEta0][iTtPhi0 + 1] = et;
1933  }
1934  }
1935 }
1936 
1937 template <class T>
1939  //we have to start by 0 in order to handle offset=-1
1940  //(however Fenix FIR has AFAK only 5 taps)
1941  double weights[] = {0., -1 / 3., -1 / 3., -1 / 3., 0., 1.};
1942 
1943  double adc2GeV = 0.;
1944  if (typeid(frame) == typeid(EBDataFrame)) {
1945  adc2GeV = 0.035;
1946  } else if (typeid(frame) == typeid(EEDataFrame)) {
1947  adc2GeV = 0.060;
1948  } else { //T is an invalid type!
1949  //TODO: replace message by a cms exception
1950  throw cms::Exception("Severe Error") << __FILE__ << ":" << __LINE__ << ": "
1951  << "this is a bug. Please report it.\n";
1952  }
1953 
1954  double acc = 0;
1955 
1956  const int n = min<int>(frame.size(), sizeof(weights) / sizeof(weights[0]));
1957 
1958  double gainInv[] = {12., 1., 6., 12};
1959 
1960  for (int i = offset; i < n; ++i) {
1961  int iframe = i + offset;
1962  if (iframe >= 0 && iframe < frame.size()) {
1963  acc += weights[i] * frame[iframe].adc() * gainInv[frame[iframe].gainId()] * adc2GeV;
1964  }
1965  }
1966  //cout << "\n";
1967  return acc;
1968 }
1969 
1971  const std::string& name) {
1972  if (!registerHist(name, ""))
1973  return nullptr; //this histo is disabled
1974  MonitorElement* result = ibook.bookFloat(name);
1975  if (result == nullptr) {
1976  throw cms::Exception("DQM") << "Failed to book integer DQM monitor element" << name;
1977  }
1978  return result;
1979 }
1980 
1982  DQMStore::IBooker& ibook, const std::string& name, const std::string& title, int nbins, double xmin, double xmax) {
1983  if (!registerHist(name, title))
1984  return nullptr; //this histo is disabled
1986  if (result == nullptr) {
1987  throw cms::Exception("Histo") << "Failed to book histogram " << name;
1988  }
1989  return result;
1990 }
1991 
1993  const std::string& name,
1994  const std::string& title,
1995  int nxbins,
1996  double xmin,
1997  double xmax,
1998  int nybins,
1999  double ymin,
2000  double ymax) {
2001  if (!registerHist(name, title))
2002  return nullptr; //this histo is disabled
2003  MonitorElement* result = ibook.book2D(name, title, nxbins, xmin, xmax, nybins, ymin, ymax);
2004  if (result == nullptr) {
2005  throw cms::Exception("Histo") << "Failed to book histogram " << name;
2006  }
2007  return result;
2008 }
2009 
2011  DQMStore::IBooker& ibook, const std::string& name, const std::string& title, int nbins, double xmin, double xmax) {
2012  if (!registerHist(name, title))
2013  return nullptr; //this histo is disabled
2014  MonitorElement* result = ibook.bookProfile(name, title, nbins, xmin, xmax, 0, 0, 0);
2015  if (result == nullptr) {
2016  throw cms::Exception("Histo") << "Failed to book histogram " << name;
2017  }
2018  return result;
2019 }
2020 
2022  const std::string& name,
2023  const std::string& title,
2024  int nbinx,
2025  double xmin,
2026  double xmax,
2027  int nbiny,
2028  double ymin,
2029  double ymax,
2030  const char* option) {
2031  if (!registerHist(name, title))
2032  return nullptr; //this histo is disabled
2033  MonitorElement* result = ibook.bookProfile2D(name, title, nbinx, xmin, xmax, nbiny, ymin, ymax, 0, 0, 0, option);
2034  if (result == nullptr) {
2035  throw cms::Exception("Histo") << "Failed to book histogram " << name;
2036  }
2037  return result;
2038 }
2039 
2041  availableHistList_.insert(pair<string, string>(name, title));
2042  return allHists_ || histList_.find(name) != histList_.end();
2043 }
2044 
2046  ebRecHits_.read(event);
2047  eeRecHits_.read(event);
2048  ebDigis_.read(event);
2049  eeDigis_.read(event);
2052  ebSrFlags_.read(event);
2053  eeSrFlags_.read(event);
2054  ebComputedSrFlags_.read(event);
2055  eeComputedSrFlags_.read(event);
2058  tps_.read(event);
2059  fedRaw_.read(event);
2060 }
2061 
2063  LogInfo log("HistoList");
2064  log << "Avalailable histograms (DQM monitor elements): \n";
2065  for (map<string, string>::iterator it = availableHistList_.begin(); it != availableHistList_.end(); ++it) {
2066  log << it->first << ": " << it->second << "\n";
2067  }
2068  log << "\nTo include an histogram add its name in the vstring parameter "
2069  "'histograms' of the EcalSelectiveReadoutValidation module\n";
2070 }
2071 
2072 double EcalSelectiveReadoutValidation::getEbEventSize(double nReadXtals) const {
2073  double ruHeaderPayload = 0.;
2074  const int firstEbDcc0 = nEeDccs / 2;
2075  for (int iDcc0 = firstEbDcc0; iDcc0 < firstEbDcc0 + nEbDccs; ++iDcc0) {
2076  ruHeaderPayload += getRuCount(iDcc0) * 8.;
2077  }
2078 
2079  return getDccOverhead(EB) * nEbDccs + nReadXtals * getBytesPerCrystal() + ruHeaderPayload;
2080 }
2081 
2082 double EcalSelectiveReadoutValidation::getEeEventSize(double nReadXtals) const {
2083  double ruHeaderPayload = 0.;
2084  const unsigned firstEbDcc0 = nEeDccs / 2;
2085  for (unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0) {
2086  //skip barrel:
2087  if (iDcc0 == firstEbDcc0)
2088  iDcc0 += nEbDccs;
2089  ruHeaderPayload += getRuCount(iDcc0) * 8.;
2090  }
2091  return getDccOverhead(EE) * nEeDccs + nReadXtals * getBytesPerCrystal() + ruHeaderPayload;
2092 }
2093 
2094 //This implementation assumes that int is coded on at least 28-bits,
2095 //which in pratice should be always true.
2097  const std::vector<int>& firWeights,
2098  int firstFIRSample,
2099  bool* saturated) {
2100  const int nFIRTaps = 6;
2101  //FIR filter weights:
2102  const vector<int>& w = firWeights;
2103 
2104  //accumulator used to compute weighted sum of samples
2105  int acc = 0;
2106  bool gain12saturated = false;
2107  const int gain12 = 0x01;
2108  const int lastFIRSample = firstFIRSample + nFIRTaps - 1;
2109  //LogDebug("DccFir") << "DCC FIR operation: ";
2110  int iWeight = 0;
2111  for (int iSample = firstFIRSample - 1; iSample < lastFIRSample; ++iSample, ++iWeight) {
2112  if (iSample >= 0 && iSample < frame.size()) {
2113  EcalMGPASample sample(frame[iSample]);
2114  if (sample.gainId() != gain12)
2115  gain12saturated = true;
2116  LogTrace("DccFir") << (iSample >= firstFIRSample ? "+" : "") << sample.adc() << "*(" << w[iWeight] << ")";
2117  acc += sample.adc() * w[iWeight];
2118  } else {
2119  edm::LogWarning("DccFir") << __FILE__ << ":" << __LINE__
2120  << ": Not enough samples in data frame or 'ecalDccZs1stSample' module "
2121  "parameter is not valid...";
2122  }
2123  }
2124  LogTrace("DccFir") << "\n";
2125  //discards the 8 LSBs
2126  //(shift operator cannot be used on negative numbers because
2127  // the result depends on compilator implementation)
2128  acc = (acc >= 0) ? (acc >> 8) : -(-acc >> 8);
2129  //ZS passed if weighted sum acc above ZS threshold or if
2130  //one sample has a lower gain than gain 12 (that is gain 12 output
2131  //is saturated)
2132 
2133  LogTrace("DccFir") << "acc: " << acc << "\n"
2134  << "saturated: " << (gain12saturated ? "yes" : "no") << "\n";
2135 
2136  if (saturated) {
2137  *saturated = gain12saturated;
2138  }
2139 
2140  return gain12saturated ? numeric_limits<int>::max() : acc;
2141 }
2142 
2143 std::vector<int> EcalSelectiveReadoutValidation::getFIRWeights(const std::vector<double>& normalizedWeights) {
2144  const int nFIRTaps = 6;
2145  vector<int> firWeights(nFIRTaps, 0); //default weight: 0;
2146  const static int maxWeight = 0xEFF; //weights coded on 11+1 signed bits
2147  for (unsigned i = 0; i < min((size_t)nFIRTaps, normalizedWeights.size()); ++i) {
2148  firWeights[i] = lround(normalizedWeights[i] * (1 << 10));
2149  if (abs(firWeights[i]) > maxWeight) { //overflow
2150  firWeights[i] = firWeights[i] < 0 ? -maxWeight : maxWeight;
2151  }
2152  }
2153  return firWeights;
2154 }
2155 
2156 void EcalSelectiveReadoutValidation::configFirWeights(const vector<double>& weightsForZsFIR) {
2157  bool notNormalized = false;
2158  bool notInt = false;
2159  for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2160  if (weightsForZsFIR[i] > 1.)
2161  notNormalized = true;
2162  if ((int)weightsForZsFIR[i] != weightsForZsFIR[i])
2163  notInt = true;
2164  }
2165  if (notInt && notNormalized) {
2166  throw cms::Exception("InvalidParameter") << "weigtsForZsFIR paramater values are not valid: they "
2167  << "must either be integer and uses the hardware representation "
2168  << "of the weights or less or equal than 1 and used the normalized "
2169  << "representation.";
2170  }
2171  LogInfo log("DccFir");
2172  if (notNormalized) {
2173  firWeights_ = vector<int>(weightsForZsFIR.size());
2174  for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2175  firWeights_[i] = (int)weightsForZsFIR[i];
2176  }
2177  } else {
2178  firWeights_ = getFIRWeights(weightsForZsFIR);
2179  }
2180 
2181  log << "Input weights for FIR: ";
2182  for (unsigned i = 0; i < weightsForZsFIR.size(); ++i) {
2183  log << weightsForZsFIR[i] << "\t";
2184  }
2185 
2186  double s2 = 0.;
2187  log << "\nActual FIR weights: ";
2188  for (unsigned i = 0; i < firWeights_.size(); ++i) {
2189  log << firWeights_[i] << "\t";
2190  s2 += firWeights_[i] * firWeights_[i];
2191  }
2192 
2193  s2 = sqrt(s2);
2194  log << "\nNormalized FIR weights after hw representation rounding: ";
2195  for (unsigned i = 0; i < firWeights_.size(); ++i) {
2196  log << firWeights_[i] / (double)(1 << 10) << "\t";
2197  }
2198 
2199  log << "\nFirst FIR sample: " << firstFIRSample_;
2200 }
2201 
2203  if (logSrpAlgoErrors_) {
2205  if (!srpAlgoErrorLog_.good()) {
2206  throw cms::Exception("Output") << "Failed to open the log file '" << srpAlgoErrorLogFileName_
2207  << "' for SRP algorithm result check.\n";
2208  }
2209  }
2210 
2213  if (!srApplicationErrorLog_.good()) {
2214  throw cms::Exception("Output") << "Failed to open the log file '" << srApplicationErrorLogFileName_
2215  << "' for Selective Readout decision application check.\n";
2216  }
2217  }
2218 }
2219 
2220 //Compares two SR flag sorted collections . Both collections
2221 //are sorted by their key (the detid) and following algorithm is based on
2222 //this feature.
2223 template <class T> //T must be either an EBSrFlagCollection or an EESrFlagCollection
2224 void EcalSelectiveReadoutValidation::compareSrfColl(const edm::Event& event, T& srfFromData, T& computedSrf) {
2225  typedef typename T::const_iterator SrFlagCollectionConstIt;
2226  typedef typename T::key_type MyRuDetIdType;
2227  SrFlagCollectionConstIt itSrfFromData = srfFromData.begin();
2228  SrFlagCollectionConstIt itComputedSr = computedSrf.begin();
2229 
2230  while (itSrfFromData != srfFromData.end() || itComputedSr != computedSrf.end()) {
2231  MyRuDetIdType inconsistentRu = 0;
2232  bool inconsistent = false;
2233  if (itComputedSr == computedSrf.end() ||
2234  (itSrfFromData != srfFromData.end() && itSrfFromData->id() < itComputedSr->id())) {
2235  //computedSrf is missig a detid found in srfFromData
2236  pair<int, int> ch = dccCh(itSrfFromData->id());
2237  srpAlgoErrorLog_ << event.id() << ": " << itSrfFromData->id() << ", DCC " << ch.first << " ch " << ch.second
2238  << " found in data (SRF:" << itSrfFromData->flagName()
2239  << ") but not in the set of SRFs computed from the data TTF.\n";
2240  inconsistentRu = itSrfFromData->id();
2241  inconsistent = true;
2242  ++itSrfFromData;
2243  } else if (itSrfFromData == srfFromData.end() ||
2244  (itComputedSr != computedSrf.end() && itComputedSr->id() < itSrfFromData->id())) {
2245  //ebSrFlags is missing a detid found in computedSrf
2246  pair<int, int> ch = dccCh(itComputedSr->id());
2247  if (logErrForDccs_[ch.first - minDccId_]) {
2248  srpAlgoErrorLog_ << event.id() << ": " << itComputedSr->id() << ", DCC " << ch.first << " ch " << ch.second
2249  << " not found in data. Computed SRF: " << itComputedSr->flagName() << ".\n";
2250  inconsistentRu = itComputedSr->id();
2251  inconsistent = true;
2252  }
2253  ++itComputedSr;
2254  } else {
2255  //*itSrfFromData and *itComputedSr has same detid
2256  if (itComputedSr->value() != itSrfFromData->value()) {
2257  pair<int, int> ch = dccCh(itSrfFromData->id());
2258  srpAlgoErrorLog_ << event.id() << ", " << itSrfFromData->id() << ", DCC " << ch.first << " ch " << ch.second
2259  << ", SRF inconsistency: "
2260  << "from data: " << itSrfFromData->flagName()
2261  << ", computed from TTF: " << itComputedSr->flagName() << "\n";
2262  inconsistentRu = itComputedSr->id();
2263  inconsistent = true;
2264  }
2265  if (itComputedSr != computedSrf.end())
2266  ++itComputedSr;
2267  if (itSrfFromData != srfFromData.end())
2268  ++itSrfFromData;
2269  }
2270 
2271  if (inconsistent)
2272  fill(meSRFlagsConsistency_, ruGraphX(inconsistentRu), ruGraphY(inconsistentRu));
2273  }
2274 }
2275 
2277 
2279  if (detId.ietaAbs() > 17) {
2280  throw cms::Exception("InvalidArgument")
2281  << "Argument of EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId&) "
2282  << "must be a barrel trigger tower Id\n";
2283  }
2284  return dccCh(detId).first;
2285 }
2286 
2288  logErrForDccs_ = vector<bool>(nDccs_, false);
2289 
2290  for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
2291  int iDcc = dccId(it->id()) - minDccId_;
2292 
2293  logErrForDccs_.at(iDcc) = true;
2294  }
2295 
2296  for (EESrFlagCollection::const_iterator it = eeSrFlags_->begin(); it != eeSrFlags_->end(); ++it) {
2297  int iDcc = dccId(it->id()) - minDccId_;
2298 
2299  logErrForDccs_.at(iDcc) = true;
2300  }
2301 
2302  stringstream buf;
2303  buf << "List of DCCs found in the first processed event: ";
2304  bool first = true;
2305  for (unsigned iDcc = 0; iDcc < nDccs_; ++iDcc) {
2306  if (logErrForDccs_[iDcc]) {
2307  buf << (first ? "" : ", ") << (iDcc + minDccId_);
2308  first = false;
2309  }
2310  }
2311  buf << "\nOnly DCCs from this list will be considered for error logging\n";
2312  srpAlgoErrorLog_ << buf.str();
2313  srApplicationErrorLog_ << buf.str();
2314  LogInfo("EcalSrValid") << buf.str();
2315 }
2316 
2317 template <class T>
2319  typedef typename T::const_iterator SrFlagCollectionConstIt;
2320  typedef typename T::key_type MyRuDetIdType;
2321 
2322  for (SrFlagCollectionConstIt itSrf = srfs.begin(); itSrf != srfs.end(); ++itSrf) {
2323  int flag = itSrf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
2324  pair<int, int> ru = dccCh(itSrf->id());
2325 
2326  if (flag == EcalSrFlag::SRF_FULL) {
2327  if (nPerRu_[ru.first - minDccId_][ru.second - 1] == getCrystalCount(ru.first, ru.second)) { //no error
2328  fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2329  fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2330  } else if (nPerRu_[ru.first - minDccId_][ru.second - 1] == 0) { //tower dropped!
2331  fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2332  fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2333  fill(meDroppedFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2334  ++nDroppedFRO_;
2335  srApplicationErrorLog_ << event.id() << ": Flag of RU " << itSrf->id() << " (DCC " << ru.first << " ch "
2336  << ru.second << ") is 'Full readout' "
2337  << "while none of its channel was read out\n";
2338  } else { //tower partially read out
2339  fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2340  fill(meDroppedFRORateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2341  fill(meIncompleteFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2342  ++nIncompleteFRO_;
2343  srApplicationErrorLog_ << event.id() << ": Flag of RU" << itSrf->id() << " (DCC " << ru.first << " ch "
2344  << ru.second << ") is 'Full readout' "
2345  << "while only " << nPerRu_[ru.first - minDccId_][ru.second - 1] << " / "
2346  << getCrystalCount(ru.first, ru.second) << " channels were read out.\n";
2347  }
2348  }
2349 
2351  if (nPerRu_[ru.first - minDccId_][ru.second - 1] == getCrystalCount(ru.first, ru.second)) {
2352  //ZS readout unit whose every channel was read
2353 
2354  fill(meCompleteZSMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()));
2355  fill(meCompleteZSRateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
2356 
2357  ++nCompleteZS_;
2358  } else {
2359  fill(meCompleteZSRateMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
2360  }
2361  }
2362  }
2363 }
2364 
2366  if (iDcc < minDccId_ || iDcc > maxDccId_) { //invalid DCC
2367  return 0;
2368  } else if (10 <= iDcc && iDcc <= 45) { //EB
2369  return 25;
2370  } else { //EE
2371  int iDccPhi;
2372  if (iDcc < 10)
2373  iDccPhi = iDcc;
2374  else
2375  iDccPhi = iDcc - 45;
2376  switch (iDccPhi * 100 + iDccCh) {
2377  case 110:
2378  case 232:
2379  case 312:
2380  case 412:
2381  case 532:
2382  case 610:
2383  case 830:
2384  case 806:
2385  //inner partials at 12, 3, and 9 o'clock
2386  return 20;
2387  case 134:
2388  case 634:
2389  case 827:
2390  case 803:
2391  return 10;
2392  case 330:
2393  case 430:
2394  return 20;
2395  case 203:
2396  case 503:
2397  case 721:
2398  case 921:
2399  return 21;
2400  default:
2401  return 25;
2402  }
2403  }
2404 }
void configFirWeights(const std::vector< double > &weightsForZsFIR)
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
MonitorElement * bookFloat(TString const &name, FUNC onbooking=NOOP())
Definition: DQMStore.h:80
CollHandle< EBDigiCollection > ebNoZsDigis_
std::ofstream srApplicationErrorLog_
Output ascii file for unconsistency between Xtals and RU Flags.
std::pair< int, int > getDCCandSC(EcalScDetId id) const
MonitorElement * bookProfile2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, double lowZ, double highZ, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:485
static const int nEeX
EE crystal grid size along X.
int zside() const
Definition: EcalScDetId.h:64
EcalTrigTowerDetId readOutUnitOf(const EBDetId &xtalId) const
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
int nIncompleteFRO_
Counter of FRO-flagged RU only partial data.
int ieta() const
get the tower ieta
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
void dqmBeginRun(edm::Run const &, edm::EventSetup const &) override
static const int SRF_ZS2
Definition: EcalSrFlag.h:21
T w() const
static const int nEbPhi
number of crystals along Phi in EB
int ruGraphY(const EcalScDetId &id) const
std::map< uint32_t, uint32_t >::const_iterator EcalTPGGroupsMapItr
Definition: EcalTPGGroups.h:20
Ecal readout channel identification [32:20] Unused (so far) [19:13] DCC id [12:6] tower [5:3] strip [...
static const int SRF_FORCED_MASK
Definition: EcalSrFlag.h:29
static const double rad2deg
Conversion factor from radian to degree.
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
double getDccSrDependentPayload(int iDcc0, double nReadRus, double nReadXtals) const
int simHit
energy reconstructed from zero-suppressed digi
CollHandle< EcalTrigPrimDigiCollection > tps_
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
T eta() const
Definition: PV3DBase.h:73
int dccId() const
get the DCC (Ecal Local DCC value not global one) id
std::pair< int, int > dccCh(const DetId &xtalId) const
void analyzeEB(const edm::Event &event, const edm::EventSetup &es)
edm::ESGetToken< EcalTPGLutIdMap, EcalTPGLutIdMapRcd > lutMapHandle
std::vector< EcalRecHit >::const_iterator const_iterator
int nEeZsErrors_
Counter of EE ZS errors (LI channel below ZS threshold)
size_t size() const
Lenght of the data buffer in bytes.
Definition: FEDRawData.h:48
CollHandle< FEDRawDataCollection > fedRaw_
void dqmEndRun(const edm::Run &r, const edm::EventSetup &c) override
void analyzeEE(const edm::Event &event, const edm::EventSetup &es)
int zside(DetId const &)
static const int nEeDccs
number of DCCs for EE
Log< level::Error, false > LogError
assert(be >=bs)
double frame2Energy(const EcalDataFrame &frame) const
unsigned ttId(DetId const &, EcalElectronicsMapping const *)
const EcalTrigTowerDetId & id() const override
Definition: EBSrFlag.h:36
double frame2EnergyForTp(const T &frame, int offset=0) const
bool registerHist(const std::string &name, const std::string &title)
void setTtEtSums(const edm::EventSetup &es, const EBDigiCollection &ebDigis, const EEDigiCollection &eeDigis)
#define LogTrace(id)
EcalTrigTowerDetId towerOf(const DetId &id) const
Get the tower id for this det id (or null if not known)
const EcalScDetId & id() const override
Definition: EESrFlag.h:37
static const int nEbEta
number of crystals along Eta in EB
T getUntrackedParameter(std::string const &, T const &) const
void Fill(long long x)
U second(std::pair< T, U > const &p)
edm::ESGetToken< EcalTrigTowerConstituentsMap, IdealGeometryRecord > hTriggerTowerMap
CollHandle< EESrFlagCollection > eeComputedSrFlags_
static const int nTtPhi
Number of Trigger Towers along Phi.
Definition: TTTypes.h:54
CollHandle< EEDigiCollection > eeDigis_
int zside() const
get the z-side of the tower (1/-1)
void checkSrApplication(const edm::Event &event, T &srfs)
std::string outputFile_
Output file for histograms.
CollHandle< RecHitCollection > eeRecHits_
unsigned towerId(DetId const &, EcalElectronicsMapping const *)
int ruGraphX(const EcalScDetId &id) const
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:408
MonitorElement * bookProfile2D(DQMStore::IBooker &, const std::string &name, const std::string &title, int nbinx, double xmin, double xmax, int nbiny, double ymin, double ymax, const char *option="")
double getDccEventSize(int iDcc0, double nReadXtals) const
void analyzeDataVolume(const edm::Event &e, const edm::EventSetup &es)
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
int towerId() const
get the tower id
T sqrt(T t)
Definition: SSEVec.h:23
edm::ESGetToken< EcalTPGLutGroup, EcalTPGLutGroupRcd > lutGrpHandle
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > geoToken
void analyze(edm::Event const &e, edm::EventSetup const &c) override
Analyzes the event.
void anaDigi(const T &frame, const U &srFlagColl)
CollHandle< EBSrFlagCollection > ebSrFlags_
edm::ESGetToken< EcalTPGPhysicsConst, EcalTPGPhysicsConstRcd > physHandle
void updateL1aRate(const edm::Event &event)
static const int SRF_FULL
Definition: EcalSrFlag.h:24
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int nCompleteZS_
Counter of ZS-flagged RU fully read out.
Transition
Definition: Transition.h:12
int value() const
Definition: EcalSrFlag.h:44
double getDccOverhead(subdet_t subdet) const
CollHandle< std::vector< PCaloHit > > eeSimHits_
static int dccZsFIR(const EcalDataFrame &frame, const std::vector< int > &firWeights, int firstFIRSample, bool *saturated=nullptr)
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
int nEeFROCnt_
Counter of EE FRO-flagged RUs.
void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
int iy() const
Definition: EcalScDetId.h:76
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
virtual TProfile * getTProfile()
void read(const edm::Event &event)
Definition: CollHandle.h:58
static const int SRF_ZS1
Definition: EcalSrFlag.h:18
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
const FEDRawData & FEDData(int fedid) const
retrieve data for fed
MonitorElement * book2D(DQMStore::IBooker &, const std::string &name, const std::string &title, int nxbins, double xmin, double xmax, int nybins, double ymin, double ymax)
CollHandle< EESrFlagCollection > eeSrFlags_
static const int nDccRus_[nDccs_]
number of RUs for each DCC
double getEbEventSize(double nReadXtals) const
const_iterator end() const
unsigned int id
Namespace of DDCMS conversion namespace.
static const int nOneEeTtEta
Number of Trigger Towers in an endcap along Eta.
int dccId(const EcalScDetId &detId) const
int nEbFROCnt_
Counter of EB FRO-flagged RUs.
static const int nEbDccs
number of DCCs for EB
Log< level::Info, false > LogInfo
int ietaAbs() const
get the absolute value of the tower ieta
CollHandle< std::vector< PCaloHit > > ebSimHits_
int nEbZsErrors_
Counter of EB ZS errors (LI channel below ZS threshold)
static const int scEdge
Number of crystals along a supercrystal edge.
ALPAKA_FN_ACC ALPAKA_FN_INLINE uint32_t ix(uint32_t id)
double recE
energy reconstructed from unsuppressed digi
static const double tmax[3]
Definition: DetId.h:17
static std::vector< int > getFIRWeights(const std::vector< double > &normalizedWeights)
static const int ebTtEdge
Number of crystals along an EB TT.
MonitorElement * book1D(DQMStore::IBooker &, const std::string &name, const std::string &title, int nbins, double xmin, double xmax)
EcalSelectiveReadoutValidation(const edm::ParameterSet &ps)
Constructor.
constexpr int gainId(sample_type sample)
get the gainId (2 bits)
static const int kByte_
number of bytes in 1 kByte:
const_iterator begin() const
The iterator returned can not safely be used across threads.
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
bool ebRuActive_[nEbEta/ebTtEdge][nEbPhi/ebTtEdge]
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:221
edm::ESGetToken< EcalElectronicsMapping, EcalMappingRcd > ecalmapping
MonitorElement * bookProfile(DQMStore::IBooker &, const std::string &name, const std::string &title, int nbins, double xmin, double xmax)
double b
Definition: hdecay.h:120
CollHandle< EEDigiCollection > eeNoZsDigis_
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
std::map< std::string, std::string > availableHistList_
std::map< uint32_t, EcalTPGLut > EcalTPGLutMap
std::vector< EcalScDetId > getEcalScDetId(int DCCid, int DCC_Channel, bool ignoreSingleCrystal=true) const
HLT enums.
EcalElectronicsId getElectronicsId(const DetId &id) const
Get the electronics id for this det id.
static const int nTtEta
Number of Trigger Towers along Eta.
double a
Definition: hdecay.h:121
ALPAKA_FN_ACC ALPAKA_FN_INLINE uint32_t iy(uint32_t id)
int nDroppedFRO_
Counter of FRO-flagged RU dropped from data.
void setToken(edm::ConsumesCollector &collector)
Definition: CollHandle.h:48
void fill(MonitorElement *me, float x)
std::map< uint32_t, EcalTPGPhysicsConst::Item > EcalTPGPhysicsConstMap
void compareSrfColl(const edm::Event &event, T &srfFromData, T &computedSrf)
Log< level::Warning, false > LogWarning
std::ofstream srpAlgoErrorLog_
Output ascii file for unconsistency on SR flags.
static const int nEeY
EE crystal grid size along Y.
CollHandle< EBDigiCollection > ebDigis_
std::map< uint32_t, EcalTPGPhysicsConst::Item >::const_iterator EcalTPGPhysicsConstMapIterator
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
int iphi() const
get the tower iphi
UseScope< MonitorElementData::Scope::RUN > UseRunScope
Definition: DQMStore.h:550
double getEeEventSize(double nReadXtals) const
long double T
Geom::Theta< T > theta() const
static const unsigned nDccs_
Total number of DCCs.
MonitorElement * bookFloat(DQMStore::IBooker &, const std::string &name)
static const int nEndcaps
number of endcaps
std::map< uint32_t, uint32_t > EcalTPGGroupsMap
Definition: EcalTPGGroups.h:19
std::map< uint32_t, EcalTPGLut >::const_iterator EcalTPGLutMapItr
const EcalTrigTowerConstituentsMap * triggerTowerMap_
void analyzeTP(const edm::Event &event, const edm::EventSetup &es)
const EcalElectronicsMapping * elecMap_
CollHandle< RecHitCollection > ebRecHits_
Definition: event.py:1
Definition: Run.h:45
bool eeRuActive_[nEndcaps][nEeX/scEdge][nEeY/scEdge]
CollHandle< EBSrFlagCollection > ebComputedSrFlags_
#define LogDebug(id)
int ix() const
Definition: EcalScDetId.h:70