CMS 3D CMS Logo

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