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
395  for (EEDigiCollection::const_iterator it = eeDigis_->begin(); it != eeDigis_->end(); ++it) {
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:
594  for (EBDigiCollection::const_iterator it = ebNoZsDigis_->begin(); it != ebNoZsDigis_->end(); ++it) {
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 
649  int nEbDigi = 0;
650 
651  for (EBDigiCollection::const_iterator it = ebDigis_->begin(); it != ebDigis_->end(); ++it) {
652  ++nEbDigi;
653  const EBDataFrame& frame = *it;
654  int iEta = static_cast<const EBDetId&>(frame.id()).ieta();
655  int iPhi = static_cast<const EBDetId&>(frame.id()).iphi();
656  int iEta0 = iEta2cIndex(iEta);
657  int iPhi0 = iPhi2cIndex(iPhi);
658  if (iEta0 < 0 || iEta0 >= nEbEta) {
659  throw(cms::Exception("EcalSelectiveReadoutValidation") << "iEta0 (= " << iEta0 << ") is out of range ("
660  << "[0," << nEbEta - 1 << "]");
661  }
662  if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
663  throw(cms::Exception("EcalSelectiveReadoutValidation") << "iPhi0 (= " << iPhi0 << ") is out of range ("
664  << "[0," << nEbPhi - 1 << "]");
665  }
666  assert(iEta0 >= 0 && iEta0 < nEbEta);
667  assert(iPhi0 >= 0 && iPhi0 < nEbPhi);
668  if (!crystalShot[iEta0][iPhi0]) {
669  crystalShot[iEta0][iPhi0] = true;
670  EBxtalCoor[iEta0][iPhi0] = make_pair(xtalGraphX(frame.id()), xtalGraphY(frame.id()));
671  } else {
672  cout << "Error: several digi for same crystal!";
673  abort();
674  }
675  if (localReco_) {
676  (*ebEnergies)[iEta0][iPhi0].recE = frame2Energy(frame);
677  }
678 
679  (*ebEnergies)[iEta0][iPhi0].gain12 = true;
680  for (int i = 0; i < frame.size(); ++i) {
681  const int gain12Code = 0x1;
682  if (frame[i].gainId() != gain12Code) {
683  (*ebEnergies)[iEta0][iPhi0].gain12 = false;
684  }
685  }
686 
688 
689  bool highInterest = false;
690 
691  // if(srf == ebSrFlags_->end()){
692  // throw cms::Exception("EcalSelectiveReadoutValidation")
693  // << __FILE__ << ":" << __LINE__ << ": SR flag not found";
694  //}
695 
696  if (srf != ebSrFlags_->end()) {
697  highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK) == EcalSrFlag::SRF_FULL);
698  }
699 
700  if (highInterest) {
702  } else {
703  int v = dccZsFIR(frame, firWeights_, firstFIRSample_, nullptr);
704  fill(meEbLiZsFir_, v);
705  if (v < ebZsThr_) {
706  eventError = true;
707  ++nEbZsErrors_;
708  pair<int, int> ru = dccCh(frame.id());
709  if (isRuComplete_[ru.first - 1][ru.second - 1])
711  if (nEbZsErrors_ < 3) {
712  srApplicationErrorLog_ << event.id() << ", "
713  << "RU " << frame.id() << ", "
714  << "DCC " << ru.first << " Ch : " << ru.second << ": "
715  << "LI channel under ZS threshold.\n";
716  }
717  if (nEbZsErrors_ == 3) {
718  srApplicationErrorLog_ << event.id() << ": "
719  << "more ZS errors for this event...\n";
720  }
721  }
722  }
723  } //next EB digi
724 
725  for (int iEta0 = 0; iEta0 < nEbEta; ++iEta0) {
726  for (int iPhi0 = 0; iPhi0 < nEbPhi; ++iPhi0)
727  fill(meChOcc_,
728  EBxtalCoor[iEta0][iPhi0].first,
729  EBxtalCoor[iEta0][iPhi0].second,
730  crystalShot[iEta0][iPhi0] ? 1. : 0.);
731  }
732 
733  if (!localReco_) {
734  for (RecHitCollection::const_iterator it = ebRecHits_->begin(); it != ebRecHits_->end(); ++it) {
735  ++nEbDigi;
736  const RecHit& hit = *it;
737  int iEta = static_cast<const EBDetId&>(hit.id()).ieta();
738  int iPhi = static_cast<const EBDetId&>(hit.id()).iphi();
739  int iEta0 = iEta2cIndex(iEta);
740  int iPhi0 = iPhi2cIndex(iPhi);
741  if (iEta0 < 0 || iEta0 >= nEbEta) {
742  LogError("EcalSrValid") << "iEta0 (= " << iEta0 << ") is out of range ("
743  << "[0," << nEbEta - 1 << "]\n";
744  }
745  if (iPhi0 < 0 || iPhi0 >= nEbPhi) {
746  LogError("EcalSrValid") << "iPhi0 (= " << iPhi0 << ") is out of range ("
747  << "[0," << nEbPhi - 1 << "]\n";
748  }
749  (*ebEnergies)[iEta0][iPhi0].recE = hit.energy();
750  }
751  }
752 
753  for (unsigned int i = 0; i < xtalEtaPhi.size(); ++i) {
754  int iEta0 = xtalEtaPhi[i].first;
755  int iPhi0 = xtalEtaPhi[i].second;
756  energiesEb_t& energies = (*ebEnergies)[iEta0][iPhi0];
757 
758  double recE = energies.recE;
759  if (recE != -numeric_limits<double>::max()) { //not zero suppressed
760  fill(meEbRecE_, (*ebEnergies)[iEta0][iPhi0].recE);
761  fill(meEbEMean_, ievt_ + 1, recE);
762  } //not zero suppressed
763 
764  if (withEbSimHit_) {
765  if (!energies.simHit) { //noise only crystal channel
766  fill(meEbNoise_, energies.noZsRecE);
767  } else {
768  fill(meEbSimE_, energies.simE);
769  fill(meEbRecEHitXtal_, energies.recE);
770  }
771  fill(meEbRecVsSimE_, energies.simE, energies.recE);
772  fill(meEbNoZsRecVsSimE_, energies.simE, energies.noZsRecE);
773  }
774  }
775 
776  int EBZs1RuCount[2][17][72];
777  int EBFullRuCount[2][17][72];
778  int EBForcedRuCount[2][17][72];
779  std::pair<int, int> EBtowerCoor[2][17][72];
780  for (int iZ(0); iZ < 2; iZ++) {
781  for (int iEta(0); iEta < 17; iEta++) {
782  for (int iPhi(0); iPhi < 72; iPhi++) {
783  EBZs1RuCount[iZ][iEta][iPhi] = 0;
784  EBFullRuCount[iZ][iEta][iPhi] = 0;
785  EBForcedRuCount[iZ][iEta][iPhi] = 0;
786  }
787  }
788  }
789 
790  //SRF
791  nEbFROCnt_ = 0;
792  char ebSrfMark[2][17][72];
793  bzero(ebSrfMark, sizeof(ebSrfMark));
794  // int idbg = 0;
795  for (EBSrFlagCollection::const_iterator it = ebSrFlags_->begin(); it != ebSrFlags_->end(); ++it) {
796  const EBSrFlag& srf = *it;
797  int iEtaAbs = srf.id().ietaAbs();
798  int iPhi = srf.id().iphi();
799  int iZ = srf.id().zside();
800 
801  // cout << "--> " << ++idbg << iEtaAbs << " " << iPhi << " " << iZ
802  // << " " << srf.id() << "\n";
803 
804  if (iEtaAbs < 1 || iEtaAbs > 17 || iPhi < 1 || iPhi > 72)
805  throw cms::Exception("EcalSelectiveReadoutValidation")
806  << "Found a barrel SRF with an invalid det ID: " << srf.id() << ".\n";
807  ++ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1];
808  if (ebSrfMark[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] > 1)
809  throw cms::Exception("EcalSelectiveReadoutValidation") << "Duplicate SRF for RU " << srf.id() << ".\n";
810 
811  EBtowerCoor[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] = std::pair<int, int>(srf.id().ieta(), srf.id().iphi());
812 
813  int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
814  if (flag == EcalSrFlag::SRF_ZS1) {
815  EBZs1RuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
816  }
817  if (flag == EcalSrFlag::SRF_FULL) {
818  EBFullRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
819  ++nEbFROCnt_;
820  }
821  if (srf.value() & EcalSrFlag::SRF_FORCED_MASK) {
822  EBForcedRuCount[iZ > 0 ? 1 : 0][iEtaAbs - 1][iPhi - 1] += 1;
823  }
824  }
825  for (int iZ(0); iZ < 2; iZ++) {
826  for (int iEta(0); iEta < 17; iEta++) {
827  for (int iPhi(0); iPhi < 72; iPhi++) {
828  float x(EBtowerCoor[iZ][iEta][iPhi].first);
829  float y(EBtowerCoor[iZ][iEta][iPhi].second);
830  fill(meZs1Ru_, x, y, EBZs1RuCount[iZ][iEta][iPhi]);
831  fill(meFullRoRu_, x, y, EBFullRuCount[iZ][iEta][iPhi]);
832  fill(meForcedRu_, x, y, EBForcedRuCount[iZ][iEta][iPhi]);
833  }
834  }
835  }
836 
837  if (eventError)
838  srApplicationErrorLog_ << event.id() << ": " << nEbZsErrors_
839  << " ZS-flagged EB channels under "
840  "the ZS threshold, whose "
841  << nEbZsErrorsType1_ << " in a complete RU.\n";
842 }
843 
845 
847  // endcap mapping
849 
850  //electronics map
852 
853  initAsciiFile();
854 }
855 
858 }
859 
861  edm::Run const&,
862  edm::EventSetup const&) {
863  ibooker.setCurrentFolder("EcalDigisV/SelectiveReadout");
864 
865  {
866  auto scope = DQMStore::IBooker::UseRunScope(ibooker);
867  meL1aRate_ = bookFloat(ibooker, "l1aRate_");
868  }
869 
870  meDccVol_ = bookProfile(ibooker,
871  "hDccVol", //"EcalDccEventSizeComputed",
872  "ECAL DCC event fragment size;Dcc id; "
873  "<Event size> (kB)",
874  nDccs_,
875  .5,
876  .5 + nDccs_);
877 
878  meDccLiVol_ = bookProfile(ibooker,
879  "hDccLiVol",
880  "LI channel payload per DCC;Dcc id; "
881  "<Event size> (kB)",
882  nDccs_,
883  .5,
884  .5 + nDccs_);
885 
886  meDccHiVol_ = bookProfile(ibooker,
887  "hDccHiVol",
888  "HI channel payload per DCC;Dcc id; "
889  "<Event size> (kB)",
890  nDccs_,
891  .5,
892  .5 + nDccs_);
893 
894  meDccVolFromData_ = bookProfile(ibooker,
895  "hDccVolFromData", //"EcalDccEventSize",
896  "ECAL DCC event fragment size;Dcc id; "
897  "<Event size> (kB)",
898  nDccs_,
899  .5,
900  .5 + nDccs_);
901 
902  meVolBLI_ = book1D(ibooker,
903  "hVolBLI", // "EBLowInterestPayload",
904  "ECAL Barrel low interest crystal data payload;"
905  "Event size (kB);Nevts",
906  100,
907  0.,
908  200.);
909 
910  meVolELI_ = book1D(ibooker,
911  "hVolELI", //"EELowInterestPayload",
912  "Endcap low interest crystal data payload;"
913  "Event size (kB);Nevts",
914  100,
915  0.,
916  200.);
917 
918  meVolLI_ = book1D(ibooker,
919  "hVolLI", //"EcalLowInterestPayload",
920  "ECAL low interest crystal data payload;"
921  "Event size (kB);Nevts",
922  100,
923  0.,
924  200.);
925 
926  meVolBHI_ = book1D(ibooker,
927  "hVolBHI", //"EBHighInterestPayload",
928  "Barrel high interest crystal data payload;"
929  "Event size (kB);Nevts",
930  100,
931  0.,
932  200.);
933 
934  meVolEHI_ = book1D(ibooker,
935  "hVolEHI", //"EEHighInterestPayload",
936  "Endcap high interest crystal data payload;"
937  "Event size (kB);Nevts",
938  100,
939  0.,
940  200.);
941 
942  meVolHI_ = book1D(ibooker,
943  "hVolHI", //"EcalHighInterestPayload",
944  "ECAL high interest crystal data payload;"
945  "Event size (kB);Nevts",
946  100,
947  0.,
948  200.);
949 
950  meVolB_ = book1D(ibooker,
951  "hVolB", //"EBEventSize",
952  "Barrel data volume;Event size (kB);Nevts",
953  100,
954  0.,
955  200.);
956 
957  meVolE_ = book1D(ibooker,
958  "hVolE", //"EEEventSize",
959  "Endcap data volume;Event size (kB);Nevts",
960  100,
961  0.,
962  200.);
963 
964  meVol_ = book1D(ibooker,
965  "hVol", //"EcalEventSize",
966  "ECAL data volume;Event size (kB);Nevts",
967  100,
968  0.,
969  200.);
970 
971  meChOcc_ = bookProfile2D(ibooker,
972  "h2ChOcc", //"EcalChannelOccupancy",
973  "ECAL crystal channel occupancy after zero suppression;"
974  "iX -200 / iEta / iX + 100;"
975  "iY / iPhi (starting from -10^{o}!);"
976  "Event count rate",
977  401,
978  -200.5,
979  200.5,
980  360,
981  .5,
982  360.5);
983 
984  //TP
985  string tpUnit;
986  if (tpInGeV_)
987  tpUnit = string("GeV");
988  else
989  tpUnit = string("TP hw unit");
990  string title;
991  title = string("Trigger primitive TT E_{T};E_{T} ") + tpUnit + string(";Event Count");
992  meTp_ = bookProfile(ibooker,
993  "hTp", //"EcalTriggerPrimitiveEt",
994  title,
995  (tpInGeV_ ? 100 : 40),
996  0.,
997  (tpInGeV_ ? 10. : 40.));
998 
999  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 
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 
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 
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):
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 
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 
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 
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 
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 
1853  return triggerTowerMap_->towerOf(xtalId);
1854 }
1855 
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 
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;
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;
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>
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 
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 
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
1990  if (result == nullptr) {
1991  throw cms::Exception("Histo") << "Failed to book histogram " << name;
1992  }
1993  return result;
1994 }
1995 
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 
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 
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 
2045  availableHistList_.insert(pair<string, string>(name, title));
2046  return allHists_ || histList_.find(name) != histList_.end();
2047 }
2048 
2050  ebRecHits_.read(event);
2051  eeRecHits_.read(event);
2052  ebDigis_.read(event);
2053  eeDigis_.read(event);
2056  ebSrFlags_.read(event);
2057  eeSrFlags_.read(event);
2058  ebComputedSrFlags_.read(event);
2059  eeComputedSrFlags_.read(event);
2062  tps_.read(event);
2063  fedRaw_.read(event);
2064 }
2065 
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.
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 
2207  if (logSrpAlgoErrors_) {
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 
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 
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 
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>
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 
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 
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 }
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:303
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:476
static const int nEeX
EE crystal grid size along X.
int zside() const
Definition: EcalScDetId.h:64
EcalTrigTowerDetId readOutUnitOf(const EBDetId &xtalId) const
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:32
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
int ix() const
Definition: EEDetId.h:77
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:45
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
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
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.
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:399
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:19
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.
bool getData(T &iHolder) const
Definition: EventSetup.h:122
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:151
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.
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)
constexpr int16_t xOffset
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:212
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:118
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:119
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:541
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