CMS 3D CMS Logo

HcalDigisValidation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HcalDigisValidation
4 // Class: HcalDigisValidation
5 //
13 //
14 // Original Author: Ali Fahim,22 R-013,+41227672649,
15 // Created: Wed Mar 23 11:42:34 CET 2011
16 //
17 //
18 
25 
27 
28  using namespace std;
29 
30  subdet_ = iConfig.getUntrackedParameter<std::string > ("subdetector", "all");
31  outputFile_ = iConfig.getUntrackedParameter<std::string > ("outputFile", "");
32 // inputLabel_ = iConfig.getParameter<std::string > ("digiLabel");
33  inputTag_ = iConfig.getParameter<edm::InputTag > ("digiTag");
34  QIE10inputTag_ = iConfig.getParameter<edm::InputTag > ("QIE10digiTag");
35  QIE11inputTag_ = iConfig.getParameter<edm::InputTag > ("QIE11digiTag");
36  emulTPsTag_ = iConfig.getParameter<edm::InputTag > ("emulTPs");
37  dataTPsTag_ = iConfig.getParameter<edm::InputTag > ("dataTPs");
38  mc_ = iConfig.getUntrackedParameter<std::string > ("mc", "no");
39  mode_ = iConfig.getUntrackedParameter<std::string > ("mode", "multi");
40  dirName_ = iConfig.getUntrackedParameter<std::string > ("dirName", "HcalDigisV/HcalDigiTask");
41  testNumber_= iConfig.getParameter<bool>("TestNumber");
42  hep17_ = iConfig.getParameter<bool>("hep17");
43 
44  // register for data access
45  if (iConfig.exists("simHits"))
46  {
47  tok_mc_ = consumes<edm::PCaloHitContainer>(iConfig.getUntrackedParameter<edm::InputTag>("simHits"));
48  }
49  tok_hbhe_ = consumes< HBHEDigiCollection >(inputTag_);
50  tok_ho_ = consumes< HODigiCollection >(inputTag_);
51  tok_hf_ = consumes< HFDigiCollection >(inputTag_);
52  tok_emulTPs_ = consumes<HcalTrigPrimDigiCollection>(emulTPsTag_);
53  if(dataTPsTag_==edm::InputTag("")) skipDataTPs = true;
54  else {
55  skipDataTPs = false;
56  tok_dataTPs_ = consumes<HcalTrigPrimDigiCollection>(dataTPsTag_);
57  }
58 
59  tok_qie10_hf_ = consumes< QIE10DigiCollection >(QIE10inputTag_);
60  tok_qie11_hbhe_ = consumes< QIE11DigiCollection >(QIE11inputTag_);
61 
62  nevent1 = 0;
63  nevent2 = 0;
64  nevent3 = 0;
65  nevent4 = 0;
66  nevtot = 0;
67 
68  msm_ = new std::map<std::string, MonitorElement*>();
69 
70  if (outputFile_.size() != 0) edm::LogInfo("OutputInfo") << " Hcal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
71  else edm::LogInfo("OutputInfo") << " Hcal Digi Task histograms will NOT be saved";
72 
73 }
74 
75 
77  delete msm_;
78 }
79 
80 
82 
84  es.get<HcalRecNumberingRecord>().get( pHRNDC );
85  hcons = &(*pHRNDC);
86 
88 
89  maxDepth_[1] = hcons->getMaxDepth(0); // HB
90  maxDepth_[2] = hcons->getMaxDepth(1); // HE
91  maxDepth_[3] = hcons->getMaxDepth(3); // HO
92  maxDepth_[4] = hcons->getMaxDepth(2); // HF
93  maxDepth_[0] = (maxDepth_[1] > maxDepth_[2] ? maxDepth_[1] : maxDepth_[2]);
94  maxDepth_[0] = (maxDepth_[0] > maxDepth_[3] ? maxDepth_[0] : maxDepth_[3]);
95  maxDepth_[0] = (maxDepth_[0] > maxDepth_[4] ? maxDepth_[0] : maxDepth_[4]); // any of HB/HE/HO/HF
96 
97  es.get<CaloGeometryRecord>().get(geometry);
98  const CaloGeometry* geo = geometry.product();
103 
104  nChannels_[1] = gHB->getHxSize(1);
105  nChannels_[2] = gHE->getHxSize(2);
106  nChannels_[3] = gHO->getHxSize(3);
107  nChannels_[4] = gHF->getHxSize(4);
108 
109  nChannels_[0] = nChannels_[1] + nChannels_[2] + nChannels_[3] + nChannels_[4];
110 
111  //std::cout << "Channels HB:" << nChannels_[1] << " HE:" << nChannels_[2] << " HO:" << nChannels_[3] << " HF:" << nChannels_[4] << std::endl;
112 
113 }
114 
116 {
117 
119 
120  // book
121  book1D(ib,"nevtot", 1, 0, 1);
122  int bnoise = 0;
123  int bmc = 0;
124  if (subdet_ == "noise") bnoise = 1;
125  if (mc_ == "yes") bmc = 1;
126  if (subdet_ == "noise" || subdet_ == "all") {
127  booking(ib,"HB", bnoise, bmc);
128  booking(ib,"HO", bnoise, bmc);
129  booking(ib,"HF", bnoise, bmc);
130  booking(ib,"HE", bnoise, bmc);
131  } else {
132  booking(ib,subdet_, 0, bmc);
133  }
134 
135  if(skipDataTPs) return;
136 
137  HistLim tp_hl_et(260, -10, 250);
138  HistLim tp_hl_ntp(640, -20, 3180);
139  HistLim tp_hl_ntp_sub(404, -20, 2000);
140  HistLim tp_hl_ieta(85, -42.5, 42.5);
141  HistLim tp_hl_iphi(74,-0.5,73.5);
142 
143 
144  book1D(ib,"HcalDigiTask_tp_et", tp_hl_et);
145  book1D(ib,"HcalDigiTask_tp_et_v0", tp_hl_et);
146  book1D(ib,"HcalDigiTask_tp_et_v1", tp_hl_et);
147  book1D(ib,"HcalDigiTask_tp_et_HB", tp_hl_et);
148  book1D(ib,"HcalDigiTask_tp_et_HE", tp_hl_et);
149  book1D(ib,"HcalDigiTask_tp_et_HF", tp_hl_et);
150  book1D(ib,"HcalDigiTask_tp_et_HF_v0", tp_hl_et);
151  book1D(ib,"HcalDigiTask_tp_et_HF_v1", tp_hl_et);
152  book1D(ib,"HcalDigiTask_tp_ntp", tp_hl_ntp);
153  book1D(ib,"HcalDigiTask_tp_ntp_v0", tp_hl_ntp);
154  book1D(ib,"HcalDigiTask_tp_ntp_v1", tp_hl_ntp);
155  book1D(ib,"HcalDigiTask_tp_ntp_HB", tp_hl_ntp_sub);
156  book1D(ib,"HcalDigiTask_tp_ntp_HE", tp_hl_ntp_sub);
157  book1D(ib,"HcalDigiTask_tp_ntp_HF", tp_hl_ntp_sub);
158  book1D(ib,"HcalDigiTask_tp_ntp_HF_v0", tp_hl_ntp_sub);
159  book1D(ib,"HcalDigiTask_tp_ntp_HF_v1", tp_hl_ntp_sub);
160  book1D(ib,"HcalDigiTask_tp_ntp_ieta", tp_hl_ieta);
161  book1D(ib,"HcalDigiTask_tp_ntp_ieta_v0", tp_hl_ieta);
162  book1D(ib,"HcalDigiTask_tp_ntp_ieta_v1", tp_hl_ieta);
163  book1D(ib,"HcalDigiTask_tp_ntp_iphi", tp_hl_iphi);
164  book1D(ib,"HcalDigiTask_tp_ntp_iphi_v0", tp_hl_iphi);
165  book1D(ib,"HcalDigiTask_tp_ntp_iphi_v1", tp_hl_iphi);
166  book1D(ib,"HcalDigiTask_tp_ntp_10_ieta", tp_hl_ieta);
167  book1D(ib,"HcalDigiTask_tp_ntp_10_ieta_v0", tp_hl_ieta);
168  book1D(ib,"HcalDigiTask_tp_ntp_10_ieta_v1", tp_hl_ieta);
169  book2D(ib,"HcalDigiTask_tp_et_ieta", tp_hl_ieta, tp_hl_et);
170  book2D(ib,"HcalDigiTask_tp_et_ieta_v0", tp_hl_ieta, tp_hl_et);
171  book2D(ib,"HcalDigiTask_tp_et_ieta_v1", tp_hl_ieta, tp_hl_et);
172  book2D(ib,"HcalDigiTask_tp_ieta_iphi", tp_hl_ieta, tp_hl_iphi);
173  book2D(ib,"HcalDigiTask_tp_ieta_iphi_v0", tp_hl_ieta, tp_hl_iphi);
174  book2D(ib,"HcalDigiTask_tp_ieta_iphi_v1", tp_hl_ieta, tp_hl_iphi);
175  bookPf(ib,"HcalDigiTask_tp_ave_et_ieta", tp_hl_ieta, tp_hl_et, " ");
176  bookPf(ib,"HcalDigiTask_tp_ave_et_ieta_v0", tp_hl_ieta, tp_hl_et, " ");
177  bookPf(ib,"HcalDigiTask_tp_ave_et_ieta_v1", tp_hl_ieta, tp_hl_et, " ");
178 
179 }
180 
181 void HcalDigisValidation::booking(DQMStore::IBooker &ib, const std::string bsubdet, int bnoise, int bmc) {
182 
183  // Adjust/Optimize binning (JR Dittmann, 16-JUL-2015)
184 
185  HistLim Ndigis(2600, 0., 2600.);
186  HistLim ndigis(520, -20., 1020.);
187  HistLim sime(200, 0., 1.0);
188 
189  HistLim digiAmp(360, -100., 7100.);
190  HistLim digiAmpWide(360, -10000., 710000.);
191  HistLim ratio(2000, -100., 3900.);
192  HistLim sumAmp(100, -500., 1500.);
193 
194  HistLim nbin(10, 0., 10.);
195 
196  HistLim pedestal(80, -1.0, 15.);
197  HistLim pedestalfC(400, -10., 30.);
198 
199  HistLim frac(80, -0.20, 1.40);
200 
201  HistLim pedLim(80, 0., 8.);
202  HistLim pedWidthLim(100, 0., 2.);
203 
204  HistLim gainLim(120, 0., 0.6);
205  HistLim gainWidthLim(160, 0., 0.32);
206 
207  HistLim ietaLim(85, -42.5, 42.5);
208  HistLim iphiLim(74, -0.5, 73.5);
209 
210  HistLim depthLim(15,-0.5,14.5);
211 
212  if (bsubdet == "HB") {
213  Ndigis = HistLim( ((int)(nChannels_[1]/100) + 1)*100, 0., (float)((int)(nChannels_[1]/100) + 1)*100);
214  } else if (bsubdet == "HE") {
215  sime = HistLim(200, 0., 1.0);
216  Ndigis = HistLim( ((int)(nChannels_[2]/100) + 1)*100, 0., (float)((int)(nChannels_[2]/100) + 1)*100);
217  } else if (bsubdet == "HF") {
218  sime = HistLim(100, 0., 100.);
219  pedLim = HistLim(100, 0., 20.);
220  pedWidthLim = HistLim(100, 0., 5.);
221  frac = HistLim(400, -4.00, 4.00);
222  Ndigis = HistLim( ((int)(nChannels_[4]/100) + 1)*100, 0., (float)((int)(nChannels_[4]/100) + 1)*100);
223  } else if (bsubdet == "HO") {
224  sime = HistLim(200, 0., 1.0);
225  gainLim = HistLim(160, 0., 1.6);
226  Ndigis = HistLim( ((int)(nChannels_[3]/100) + 1)*100, 0., (float)((int)(nChannels_[3]/100) + 1)*100);
227  }
228 
229  int isubdet=0;
230  if (bsubdet=="HB") isubdet=1;
231  else if (bsubdet=="HE") isubdet=2;
232  else if (bsubdet=="HO") isubdet=3;
233  else if (bsubdet=="HF") isubdet=4;
234  else edm::LogWarning("HcalDigisValidation") << "HcalDigisValidation Warning: not HB/HE/HF/HO " << bsubdet << std::endl;
235 
236  Char_t histo[100];
237  const char * sub = bsubdet.c_str();
238  if (bnoise == 0) {
239  // number of digis in each subdetector
240  sprintf(histo, "HcalDigiTask_Ndigis_%s", sub);
241  book1D(ib, histo, Ndigis);
242 
243  // maps of occupancies
244  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
245  sprintf(histo, "HcalDigiTask_ieta_iphi_occupancy_map_depth%d_%s", depth, sub);
246  book2D(ib, histo, ietaLim, iphiLim);
247  }
248 
249  //Depths
250  sprintf(histo, "HcalDigiTask_depths_%s",sub);
251  book1D(ib,histo, depthLim);
252 
253  // occupancies vs ieta
254  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
255  sprintf(histo, "HcalDigiTask_occupancy_vs_ieta_depth%d_%s", depth, sub);
256  book1D(ib, histo, ietaLim);
257  }
258 
259 
260  // just 1D of all cells' amplitudes
261  sprintf(histo, "HcalDigiTask_sum_all_amplitudes_%s", sub);
262  book1D(ib, histo, sumAmp);
263 
264  sprintf(histo, "HcalDigiTask_number_of_amplitudes_above_10fC_%s", sub);
265  book1D(ib, histo, ndigis);
266 
267  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
268  sprintf(histo, "HcalDigiTask_ADC0_adc_depth%d_%s", depth, sub);
269  book1D(ib, histo, pedestal);
270  }
271 
272  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
273  sprintf(histo, "HcalDigiTask_ADC0_fC_depth%d_%s", depth, sub);
274  book1D(ib, histo, pedestalfC);
275  }
276 
277  sprintf(histo, "HcalDigiTask_signal_amplitude_%s", sub);
278  book1D(ib, histo, digiAmp);
279 
280  if(hep17_ && bsubdet=="HE"){
281  sprintf(histo, "HcalDigiTask_signal_amplitude_HEP17");
282  book1D(ib, histo, digiAmpWide);
283  }
284  //
285  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
286  sprintf(histo, "HcalDigiTask_signal_amplitude_depth%d_%s", depth, sub);
287  book1D(ib, histo, digiAmp);
288  if(hep17_ && bsubdet=="HE"){
289  sprintf(histo, "HcalDigiTask_signal_amplitude_depth%d_HEP17", depth);
290  book1D(ib, histo, digiAmpWide);
291  }
292  }
293 
294  sprintf(histo, "HcalDigiTask_signal_amplitude_vs_bin_all_depths_%s", sub);
295  book2D(ib, histo, nbin, digiAmp);
296  if(hep17_ && bsubdet=="HE"){
297  sprintf(histo, "HcalDigiTask_signal_amplitude_vs_bin_all_depths_HEP17");
298  book2D(ib, histo, nbin, digiAmpWide);
299  }
300 
301  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
302  sprintf(histo, "HcalDigiTask_all_amplitudes_vs_bin_1D_depth%d_%s", depth, sub);
303  book1D(ib, histo, nbin);
304  if(hep17_ && bsubdet=="HE"){
305  sprintf(histo, "HcalDigiTask_all_amplitudes_vs_bin_1D_depth%d_HEP17", depth);
306  book1D(ib, histo, nbin);
307  }
308  }
309 
310  sprintf(histo, "HcalDigiTask_bin_5_frac_%s", sub);
311  book1D(ib, histo, frac);
312  sprintf(histo, "HcalDigiTask_bin_6_7_frac_%s", sub);
313  book1D(ib, histo, frac);
314 
315  if (bmc == 1) {
316  sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_%s", sub);
317  book2D(ib, histo, sime, digiAmp);
318  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
319  sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_depth%d_%s", depth, sub);
320  book2D(ib, histo, sime, digiAmp);
321  }
322 
323  sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_profile_%s", sub);
324  bookPf(ib, histo, sime, digiAmp);
325  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
326  sprintf(histo, "HcalDigiTask_amplitude_vs_simhits_profile_depth%d_%s", depth, sub);
327  bookPf(ib, histo, sime, digiAmp);
328  }
329 
330  sprintf(histo, "HcalDigiTask_ratio_amplitude_vs_simhits_%s", sub);
331  book1D(ib, histo, ratio);
332  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
333  sprintf(histo, "HcalDigiTask_ratio_amplitude_vs_simhits_depth%d_%s", depth, sub);
334  book1D(ib, histo, ratio);
335  }
336  }//mc only
337 
338  } else { // noise only
339 
340  // EVENT "1" distributions of all cells properties
341 
342  //KH
343  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
344  sprintf(histo, "HcalDigiTask_gain_capId0_Depth%d_%s", depth, sub);
345  book1D(ib, histo, gainLim);
346  sprintf(histo, "HcalDigiTask_gain_capId1_Depth%d_%s", depth, sub);
347  book1D(ib, histo, gainLim);
348  sprintf(histo, "HcalDigiTask_gain_capId2_Depth%d_%s", depth, sub);
349  book1D(ib, histo, gainLim);
350  sprintf(histo, "HcalDigiTask_gain_capId3_Depth%d_%s", depth, sub);
351  book1D(ib, histo, gainLim);
352 
353  sprintf(histo, "HcalDigiTask_gainWidth_capId0_Depth%d_%s", depth, sub);
354  book1D(ib, histo, gainWidthLim);
355  sprintf(histo, "HcalDigiTask_gainWidth_capId1_Depth%d_%s", depth, sub);
356  book1D(ib, histo, gainWidthLim);
357  sprintf(histo, "HcalDigiTask_gainWidth_capId2_Depth%d_%s", depth, sub);
358  book1D(ib, histo, gainWidthLim);
359  sprintf(histo, "HcalDigiTask_gainWidth_capId3_Depth%d_%s", depth, sub);
360  book1D(ib, histo, gainWidthLim);
361 
362  sprintf(histo, "HcalDigiTask_pedestal_capId0_Depth%d_%s", depth, sub);
363  book1D(ib, histo, pedLim);
364  sprintf(histo, "HcalDigiTask_pedestal_capId1_Depth%d_%s", depth, sub);
365  book1D(ib, histo, pedLim);
366  sprintf(histo, "HcalDigiTask_pedestal_capId2_Depth%d_%s", depth, sub);
367  book1D(ib, histo, pedLim);
368  sprintf(histo, "HcalDigiTask_pedestal_capId3_Depth%d_%s", depth, sub);
369  book1D(ib, histo, pedLim);
370 
371  sprintf(histo, "HcalDigiTask_pedestal_width_capId0_Depth%d_%s", depth, sub);
372  book1D(ib, histo, pedWidthLim);
373  sprintf(histo, "HcalDigiTask_pedestal_width_capId1_Depth%d_%s", depth, sub);
374  book1D(ib, histo, pedWidthLim);
375  sprintf(histo, "HcalDigiTask_pedestal_width_capId2_Depth%d_%s", depth, sub);
376  book1D(ib, histo, pedWidthLim);
377  sprintf(histo, "HcalDigiTask_pedestal_width_capId3_Depth%d_%s", depth, sub);
378  book1D(ib, histo, pedWidthLim);
379 
380  }
381 
382  //KH
383  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
384  sprintf(histo, "HcalDigiTask_gainMap_Depth%d_%s", depth, sub);
385  book2D(ib, histo, ietaLim, iphiLim);
386  sprintf(histo, "HcalDigiTask_pwidthMap_Depth%d_%s", depth, sub);
387  book2D(ib, histo, ietaLim, iphiLim);
388  }
389 
390  } //end of noise-only
391 }//book
392 
394  using namespace edm;
395  using namespace std;
396 
397  iSetup.get<HcalDbRecord > ().get(conditions);
398 
399  //TP Code
401  iSetup.get<CaloTPGRecord>().get(decoder);
402 
404  iSetup.get<CaloGeometryRecord>().get(tp_geometry);
405 
406  iSetup.get<HcalRecNumberingRecord>().get(htopo);
407 
408  //Get all handles
410  iEvent.getByToken(tok_emulTPs_, emulTPs);
411 
413  if(!skipDataTPs) iEvent.getByToken(tok_dataTPs_, dataTPs);
414  //iEvent.getByLabel("hcalDigis", dataTPs);
415 
416  //~TP Code
417 
418  if (subdet_ != "all") {
419  noise_ = 0;
420  if (subdet_ == "HB") reco<HBHEDataFrame > (iEvent, iSetup, tok_hbhe_);
421  if (subdet_ == "HE"){
422  reco<HBHEDataFrame > (iEvent, iSetup, tok_hbhe_);
423  reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
424  }
425  if (subdet_ == "HO") reco<HODataFrame > (iEvent, iSetup, tok_ho_);
426  if (subdet_ == "HF"){
427  reco<HFDataFrame > (iEvent, iSetup, tok_hf_);
428  reco<QIE10DataFrame>(iEvent, iSetup, tok_qie10_hf_);
429  }
430 
431  if (subdet_ == "noise") {
432  noise_ = 1;
433  subdet_ = "HB";
434  reco<HBHEDataFrame > (iEvent, iSetup, tok_hbhe_);
435  subdet_ = "HE";
436  reco<HBHEDataFrame > (iEvent, iSetup, tok_hbhe_);
437  reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
438  subdet_ = "HO";
439  reco<HODataFrame > (iEvent, iSetup, tok_ho_);
440  subdet_ = "HF";
441  reco<HFDataFrame > (iEvent, iSetup, tok_hf_);
442  reco<QIE10DataFrame>(iEvent, iSetup, tok_qie10_hf_);
443  subdet_ = "noise";
444  }
445  }// all subdetectors
446  else {
447  noise_ = 0;
448 
449  subdet_ = "HB";
450  reco<HBHEDataFrame > (iEvent, iSetup, tok_hbhe_);
451  subdet_ = "HE";
452  reco<HBHEDataFrame > (iEvent, iSetup, tok_hbhe_);
453  reco<QIE11DataFrame>(iEvent, iSetup, tok_qie11_hbhe_);
454  subdet_ = "HO";
455  reco<HODataFrame > (iEvent, iSetup, tok_ho_);
456  subdet_ = "HF";
457  reco<HFDataFrame > (iEvent, iSetup, tok_hf_);
458  reco<QIE10DataFrame>(iEvent, iSetup, tok_qie10_hf_);
459  subdet_ = "all";
460  }
461 
462  fill1D("nevtot", 0);
463  nevtot++;
464 
465  //TP Code
466  //Counters
467  int c = 0, cv0 = 0, cv1 = 0, chb = 0, che = 0, chf = 0, chfv0 = 0, chfv1 = 0;
468 
469  if(skipDataTPs) return;
470 
471  for (HcalTrigPrimDigiCollection::const_iterator itr = dataTPs->begin(); itr != dataTPs->end(); ++itr) {
472  int ieta = itr->id().ieta();
473  int iphi = itr->id().iphi();
474 
475  HcalSubdetector subdet = (HcalSubdetector) 0;
476 
477 
478  if ( abs(ieta) <= 16 )
479  subdet = HcalSubdetector::HcalBarrel ;
480  else if ( abs(ieta) < tp_geometry->firstHFTower(itr->id().version()) )
481  subdet = HcalSubdetector::HcalEndcap ;
482  else if ( abs(ieta) <= 42 )
484 
485  //Right now, the only case where version matters is in HF
486  //If the subdetector is not HF, set version to -1
487  int tpVersion = (subdet == HcalSubdetector::HcalForward ? itr->id().version() : -1);
488 
489  float en = decoder->hcaletValue(itr->id(), itr->t0());
490 
491  if (en < 0.00001) continue;
492 
493  //Plot the variables
494  //Retain classic behavior (include all tps)
495  //Additional plots that only include HF 3x2 or HF 1x1
496 
497  //Classic
498  fill1D("HcalDigiTask_tp_et",en);
499  fill2D("HcalDigiTask_tp_et_ieta",ieta,en);
500  fill2D("HcalDigiTask_tp_ieta_iphi",ieta,iphi);
501  fillPf("HcalDigiTask_tp_ave_et_ieta",ieta,en);
502  fill1D("HcalDigiTask_tp_ntp_ieta",ieta);
503  fill1D("HcalDigiTask_tp_ntp_iphi",iphi);
504  if ( en > 10. ) fill1D("HcalDigiTask_tp_ntp_10_ieta",ieta);
505 
506  //3x2 Trig Primitives (tpVersion == 0)
507  if( subdet != HcalSubdetector::HcalForward || tpVersion==0){
508  fill1D("HcalDigiTask_tp_et_v0",en);
509  fill2D("HcalDigiTask_tp_et_ieta_v0",ieta,en);
510  fill2D("HcalDigiTask_tp_ieta_iphi_v0",ieta,iphi);
511  fillPf("HcalDigiTask_tp_ave_et_ieta_v0",ieta,en);
512  fill1D("HcalDigiTask_tp_ntp_ieta_v0",ieta);
513  fill1D("HcalDigiTask_tp_ntp_iphi_v0",iphi);
514  if ( en > 10. ) fill1D("HcalDigiTask_tp_ntp_10_ieta_v0",ieta);
515  }
516 
517  //1x1 Trig Primitives (tpVersion == 1)
518  if( subdet != HcalSubdetector::HcalForward || tpVersion==1){
519  fill1D("HcalDigiTask_tp_et_v1",en);
520  fill2D("HcalDigiTask_tp_et_ieta_v1",ieta,en);
521  fill2D("HcalDigiTask_tp_ieta_iphi_v1",ieta,iphi);
522  fillPf("HcalDigiTask_tp_ave_et_ieta_v1",ieta,en);
523  fill1D("HcalDigiTask_tp_ntp_ieta_v1",ieta);
524  fill1D("HcalDigiTask_tp_ntp_iphi_v1",iphi);
525  if ( en > 10. ) fill1D("HcalDigiTask_tp_ntp_10_ieta_v1",ieta);
526  }
527 
528  ++c;
529  if ( subdet == HcalSubdetector::HcalBarrel ) {
530  fill1D("HcalDigiTask_tp_et_HB",en);
531  ++chb;
532  ++cv0;
533  ++cv1;
534  }
535  if ( subdet == HcalSubdetector::HcalEndcap ) {
536  fill1D("HcalDigiTask_tp_et_HE",en);
537  ++che;
538  ++cv0;
539  ++cv1;
540  }
541  if ( subdet == HcalSubdetector::HcalForward ) {
542  fill1D("HcalDigiTask_tp_et_HF",en);
543  ++chf;
544 
545  if(tpVersion == 0){
546  fill1D("HcalDigiTask_tp_et_HF_v0",en);
547  ++chfv0;
548  ++cv0;
549  }
550 
551  if(tpVersion == 1){
552  fill1D("HcalDigiTask_tp_et_HF_v1",en);
553  ++chfv1;
554  ++cv1;
555  }
556 
557  }
558 
559  }//end data TP collection
560 
561  fill1D("HcalDigiTask_tp_ntp",c);
562  fill1D("HcalDigiTask_tp_ntp_v0",cv0);
563  fill1D("HcalDigiTask_tp_ntp_v1",cv1);
564  fill1D("HcalDigiTask_tp_ntp_HB",chb);
565  fill1D("HcalDigiTask_tp_ntp_HE",che);
566  fill1D("HcalDigiTask_tp_ntp_HF",chf);
567  fill1D("HcalDigiTask_tp_ntp_HF_v0",chfv0);
568  fill1D("HcalDigiTask_tp_ntp_HF_v1",chfv1);
569 
570  //~TP Code
571 }
572 
573 template<class Digi> void HcalDigisValidation::reco(const edm::Event& iEvent, const edm::EventSetup& iSetup, const edm::EDGetTokenT<edm::SortedCollection<Digi> > & tok) {
574 
575 
576  // HistLim =============================================================
577 
578  std::string strtmp;
579 
580  // ======================================================================
581  using namespace edm;
584 
585  // ADC2fC
587  CaloSamples tool;
588  iEvent.getByToken(tok, digiCollection);
589  if (!digiCollection.isValid()) return;
590  int isubdet = 0;
591  if (subdet_ == "HB") isubdet = 1;
592  if (subdet_ == "HE") isubdet = 2;
593  if (subdet_ == "HO") isubdet = 3;
594  if (subdet_ == "HF") isubdet = 4;
595 
596  if (isubdet == 1) nevent1++;
597  if (isubdet == 2) nevent2++;
598  if (isubdet == 3) nevent3++;
599  if (isubdet == 4) nevent4++;
600 
601  int indigis = 0;
602  // amplitude for signal cell at diff. depths
603  std::vector<double> v_ampl_c(maxDepth_[isubdet]+1,0);
604 
605  // is set to 1 if "seed" SimHit is found
606  int seedSimHit = 0;
607 
608  int ieta_Sim = 9999;
609  int iphi_Sim = 9999;
610  double emax_Sim = -9999.;
611 
612 
613  // SimHits MC only
614  if (mc_ == "yes") {
616  iEvent.getByToken(tok_mc_, hcalHits);
617  const edm::PCaloHitContainer * simhitResult = hcalHits.product();
618 
619  if (isubdet != 0 && noise_ == 0) { // signal only SimHits
620 
621  for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end(); ++simhits) {
622 
623  unsigned int id_ = simhits->id();
624  int sub, ieta, iphi;
625  HcalDetId hid;
627  else hid = HcalDetId(id_);
628  sub = hid.subdet();
629  ieta = hid.ieta();
630  iphi = hid.iphi();
631 
632  double en = simhits->energy();
633 
634  if (en > emax_Sim && sub == isubdet) {
635  emax_Sim = en;
636  ieta_Sim = ieta;
637  iphi_Sim = iphi;
638  // to limit "seed" SimHit energy in case of "multi" event
639  if (mode_ == "multi" &&
640  ((sub == 4 && en < 100. && en > 1.)
641  || ((sub != 4) && en < 1. && en > 0.02))) {
642  seedSimHit = 1;
643  break;
644  }
645  }
646 
647  } // end of SimHits cycle
648 
649 
650  // found highest-energy SimHit for single-particle
651  if (mode_ != "multi" && emax_Sim > 0.) seedSimHit = 1;
652  } // end of SimHits
653  }// end of mc_ == "yes"
654 
655  // CYCLE OVER CELLS ========================================================
656  int Ndig = 0;
657 
658  for (digiItr = digiCollection->begin(); digiItr != digiCollection->end(); digiItr++) {
659 
660  HcalDetId cell(digiItr->id());
661  int depth = cell.depth();
662  int iphi = cell.iphi();
663  int ieta = cell.ieta();
664  int sub = cell.subdet();
665 
666  if(depth > maxDepth_[isubdet] && sub == isubdet){
667  edm::LogWarning("HcalDetId") << "HcalDetID presents conflicting information. Depth: " << depth << ", iphi: " << iphi << ", ieta: " << ieta << ". Max depth from geometry is: " << maxDepth_[isubdet] << ". TestNumber = " << testNumber_;
668  continue;
669  }
670 
671  // amplitude for signal cell at diff. depths
672  std::vector<double> v_ampl(maxDepth_[isubdet]+1,0);
673 
674  // Gains, pedestals (once !) and only for "noise" case
675  if (((nevent1 == 1 && isubdet == 1) ||
676  (nevent2 == 1 && isubdet == 2) ||
677  (nevent3 == 1 && isubdet == 3) ||
678  (nevent4 == 1 && isubdet == 4)) && noise_ == 1 && sub == isubdet) {
679 
680  HcalGenericDetId hcalGenDetId(digiItr->id());
681  const HcalPedestal* pedestal = conditions->getPedestal(hcalGenDetId);
682  const HcalGain* gain = conditions->getGain(hcalGenDetId);
683  const HcalGainWidth* gainWidth = conditions->getGainWidth(hcalGenDetId);
684  const HcalPedestalWidth* pedWidth = conditions-> getPedestalWidth(hcalGenDetId);
685 
686  for (int i = 0; i < 4; i++) {
687  fill1D("HcalDigiTask_gain_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gain->getValue(i));
688  fill1D("HcalDigiTask_gainWidth_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gainWidth->getValue(i));
689  fill1D("HcalDigiTask_pedestal_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, pedestal->getValue(i));
690  fill1D("HcalDigiTask_pedestal_width_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, pedWidth->getWidth(i));
691  }
692 
693  fill2D("HcalDigiTask_gainMap_Depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi), gain->getValue(0));
694  fill2D("HcalDigiTask_pwidthMap_Depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi), pedWidth->getWidth(0));
695 
696  }// end of event #1
697 
698  if (sub == isubdet) Ndig++; // subdet number of digi
699 
700  // No-noise case, only single subdet selected ===========================
701 
702  if (sub == isubdet && noise_ == 0) {
703 
704  HcalCalibrations calibrations = conditions->getHcalCalibrations(cell);
705 
706  const HcalQIECoder* channelCoder = conditions->getHcalCoder(cell);
707  const HcalQIEShape* shape = conditions->getHcalShape(channelCoder);
708  HcalCoderDb coder(*channelCoder, *shape);
709  coder.adc2fC(*digiItr, tool);
710 
711  double noiseADC = (*digiItr)[0].adc();
712  double noisefC = tool[0];
713  // noise evaluations from "pre-samples"
714  fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
715  fill1D("HcalDigiTask_ADC0_fC_depth" + str(depth) + "_" + subdet_, noisefC);
716 
717 
718  // OCCUPANCY maps fill
719  fill2D("HcalDigiTask_ieta_iphi_occupancy_map_depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi));
720 
721  fill1D("HcalDigiTask_depths_"+subdet_,double(depth));
722 
723  // Cycle on time slices
724  // - for each Digi
725  // - for one Digi with max SimHits E in subdet
726 
727 
728  int closen = 0; // =1 if 1) seedSimHit = 1 and 2) the cell is the same
729  if (ieta == ieta_Sim && iphi == iphi_Sim) closen = seedSimHit;
730 
731  for (int ii = 0; ii < tool.size(); ii++) {
732  int capid = (*digiItr)[ii].capid();
733  // single ts amplitude
734  double val = (tool[ii] - calibrations.pedestal(capid));
735 
736  if (val > 100.) {
737  fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
738  strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_" + subdet_;
739  fill1D(strtmp, double(ii), val);
740  }
741 
742  if (closen == 1) {
743  strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_" + subdet_;
744  fill2D(strtmp, double(ii), val);
745  }
746 
747 
748  // HB/HE/HO
749  if (isubdet != 4 && ii >= 4 && ii <= 7) {
750  v_ampl[0] += val;
751  v_ampl[depth] += val;
752 
753  if (closen == 1) {
754  v_ampl_c[0] += val;
755  v_ampl_c[depth] += val;
756  }
757  }
758 
759  // HF
760  if (isubdet == 4 && ii >= 2 && ii <= 4) {
761  v_ampl[0] += val;
762  v_ampl[depth] += val;
763  if (closen == 1) {
764  v_ampl_c[0] += val;
765  v_ampl_c[depth] += val;
766  }
767  }
768  }
769  // end of time bucket sample
770 
771  // maps of sum of amplitudes (sum lin.digis(4,5,6,7) - ped) all depths
772  // just 1D of all cells' amplitudes
773  strtmp = "HcalDigiTask_sum_all_amplitudes_" + subdet_;
774  fill1D(strtmp, v_ampl[0]);
775 
776  std::vector<int> v_ampl_sub(v_ampl.begin() + 1, v_ampl.end()); // remove element 0, which is the sum of any depth
777  double ampl_max = *std::max_element(v_ampl_sub.begin(),v_ampl_sub.end());
778  if (ampl_max>10.) indigis++;
779  //KH if (ampl1 > 10. || ampl2 > 10. || ampl3 > 10. || ampl4 > 10.) indigis++;
780 
781  // fraction 5,6 bins if ampl. is big.
782  if (v_ampl[1] > 30. && depth == 1 && closen == 1 && isubdet != 4) {
783  double fBin5 = tool[4] - calibrations.pedestal((*digiItr)[4].capid());
784  double fBin67 = tool[5] + tool[6]
785  - calibrations.pedestal((*digiItr)[5].capid())
786  - calibrations.pedestal((*digiItr)[6].capid());
787 
788  fBin5 /= v_ampl[1];
789  fBin67 /= v_ampl[1];
790 
791  strtmp = "HcalDigiTask_bin_5_frac_" + subdet_;
792  fill1D(strtmp, fBin5);
793  strtmp = "HcalDigiTask_bin_6_7_frac_" + subdet_;
794  fill1D(strtmp, fBin67);
795 
796  }
797 
798  //Special for HF
799  if (isubdet == 4 && v_ampl[1] > 30. && depth == 1) {
800  double fBin5 = tool[2] - calibrations.pedestal((*digiItr)[2].capid());
801  double fBin67 = tool[3] + tool[4]
802  - calibrations.pedestal((*digiItr)[3].capid())
803  - calibrations.pedestal((*digiItr)[4].capid());
804  fBin5 /= v_ampl[1];
805  fBin67 /= v_ampl[1];
806  strtmp = "HcalDigiTask_bin_5_frac_" + subdet_;
807  fill1D(strtmp, fBin5);
808  strtmp = "HcalDigiTask_bin_6_7_frac_" + subdet_;
809  fill1D(strtmp, fBin67);
810  }
811 
812  strtmp = "HcalDigiTask_signal_amplitude_" + subdet_;
813  fill1D(strtmp, v_ampl[0]);
814  strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_" + subdet_;
815  fill1D(strtmp, v_ampl[depth]);
816  }
817  } // End of CYCLE OVER CELLS =============================================
818 
819  if (isubdet != 0 && noise_ == 0) { // signal only, once per event
820  strtmp = "HcalDigiTask_number_of_amplitudes_above_10fC_" + subdet_;
821  fill1D(strtmp, indigis);
822 
823  // SimHits once again !!!
824  double eps = 1.e-3;
825  std::vector<double> v_ehits(maxDepth_[isubdet]+1,0);
826 
827  if (mc_ == "yes") {
829  iEvent.getByToken(tok_mc_, hcalHits);
830  const edm::PCaloHitContainer * simhitResult = hcalHits.product();
831  for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end(); ++simhits) {
832 
833  unsigned int id_ = simhits->id();
834  int sub, depth, ieta, iphi;
835  HcalDetId hid;
837  else hid = HcalDetId(id_);
838  sub = hid.subdet();
839  depth = hid.depth();
840  ieta = hid.ieta();
841  iphi = hid.iphi();
842 
843  if(depth > maxDepth_[isubdet] && sub == isubdet){
844  edm::LogWarning("HcalDetId") << "HcalDetID(SimHit) presents conflicting information. Depth: " << depth << ", iphi: " << iphi << ", ieta: " << ieta << ". Max depth from geometry is: " << maxDepth_[isubdet] << ". TestNumber = " << testNumber_;
845  continue;
846  }
847 
848 
849  // take cell already found to be max energy in a particular subdet
850  if (sub == isubdet && ieta == ieta_Sim && iphi == iphi_Sim) {
851  double en = simhits->energy();
852 
853  v_ehits[0] += en;
854  v_ehits[depth] += en;
855 
856  }
857  } // simhit loop
858 
859  strtmp = "HcalDigiTask_amplitude_vs_simhits_" + subdet_;
860  if (v_ehits[0] > eps) fill2D(strtmp, v_ehits[0], v_ampl_c[0]);
861  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
862  strtmp = "HcalDigiTask_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
863  if (v_ehits[depth] > eps) fill2D(strtmp, v_ehits[depth], v_ampl_c[depth]);
864  }
865 
866  strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_" + subdet_;
867  if (v_ehits[0] > eps) fillPf(strtmp, v_ehits[0], v_ampl_c[0]);
868  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
869  strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
870  if (v_ehits[depth] > eps) fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
871  }
872 
873  strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_" + subdet_;
874  if (v_ehits[0] > eps) fill1D(strtmp, v_ampl_c[0] / v_ehits[0]);
875  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
876  strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
877  if (v_ehits[depth] > eps) fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
878  strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
879  if (v_ehits[depth] > eps) fill1D(strtmp, v_ampl_c[depth] / v_ehits[depth]);
880  }
881 
882 
883  } // end of if(mc_ == "yes")
884 
885  strtmp = "HcalDigiTask_Ndigis_" + subdet_;
886  fill1D(strtmp, double(Ndig));
887 
888  } // end of if( subdet != 0 && noise_ == 0) { // signal only
889 
890 }
891 template<class dataFrameType> void HcalDigisValidation::reco(const edm::Event& iEvent, const edm::EventSetup& iSetup, const edm::EDGetTokenT<HcalDataFrameContainer<dataFrameType> > & tok) {
892 
893 
894  // HistLim =============================================================
895 
896  std::string strtmp;
897 
898  // ======================================================================
899  using namespace edm;
901  //typename HcalDataFrameContainer<dataFrameType>::const_iterator digiItr;
902 
903 
904  // ADC2fC
906  CaloSamples tool;
907  iEvent.getByToken(tok, digiHandle);
908  if (!digiHandle.isValid()) return;
910  int isubdet = 0;
911  if (subdet_ == "HB") isubdet = 1;
912  if (subdet_ == "HE") isubdet = 2;
913  if (subdet_ == "HO") isubdet = 3;
914  if (subdet_ == "HF") isubdet = 4;
915 
916  if (isubdet == 1) nevent1++;
917  if (isubdet == 2) nevent2++;
918  if (isubdet == 3) nevent3++;
919  if (isubdet == 4) nevent4++;
920 
921  int indigis = 0;
922  // amplitude for signal cell at diff. depths
923  std::vector<double> v_ampl_c(maxDepth_[isubdet]+1,0);
924 
925  // is set to 1 if "seed" SimHit is found
926  int seedSimHit = 0;
927 
928  int ieta_Sim = 9999;
929  int iphi_Sim = 9999;
930  double emax_Sim = -9999.;
931 
932 
933  // SimHits MC only
934  if (mc_ == "yes") {
936  iEvent.getByToken(tok_mc_, hcalHits);
937  const edm::PCaloHitContainer * simhitResult = hcalHits.product();
938 
939  if (isubdet != 0 && noise_ == 0) { // signal only SimHits
940 
941  for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end(); ++simhits) {
942 
943  unsigned int id_ = simhits->id();
944  int sub, ieta, iphi;
945  HcalDetId hid;
947  else hid = HcalDetId(id_);
948  sub = hid.subdet();
949  ieta = hid.ieta();
950  iphi = hid.iphi();
951 
952  double en = simhits->energy();
953 
954  if (en > emax_Sim && sub == isubdet) {
955  emax_Sim = en;
956  ieta_Sim = ieta;
957  iphi_Sim = iphi;
958  // to limit "seed" SimHit energy in case of "multi" event
959  if (mode_ == "multi" &&
960  ((sub == 4 && en < 100. && en > 1.)
961  || ((sub != 4) && en < 1. && en > 0.02))) {
962  seedSimHit = 1;
963  break;
964  }
965  }
966 
967  } // end of SimHits cycle
968 
969 
970  // found highest-energy SimHit for single-particle
971  if (mode_ != "multi" && emax_Sim > 0.) seedSimHit = 1;
972  } // end of SimHits
973  }// end of mc_ == "yes"
974 
975  // CYCLE OVER CELLS ========================================================
976  int Ndig = 0;
977 
978  for (typename HcalDataFrameContainer<dataFrameType>::const_iterator digiItr = digiCollection->begin(); digiItr != digiCollection->end(); digiItr++) {
979 
980  dataFrameType dataFrame = *digiItr;
981 
982  HcalDetId cell(digiItr->id());
983  int depth = cell.depth();
984  int iphi = cell.iphi();
985  int ieta = cell.ieta();
986  int sub = cell.subdet();
987 
988  //Is this in HEP17
989  bool isHEP17 = (iphi>=63)&&(iphi<=66)&&(ieta>0)&&(sub == 2);
990 
991 
992  if(depth > maxDepth_[isubdet] && sub == isubdet){
993  edm::LogWarning("HcalDetId") << "HcalDetID presents conflicting information. Depth: " << depth << ", iphi: " << iphi << ", ieta: " << ieta << ". Max depth from geometry is: " << maxDepth_[isubdet] << ". TestNumber = " << testNumber_;
994  continue;
995  }
996 
997  // amplitude for signal cell at diff. depths
998  std::vector<double> v_ampl(maxDepth_[isubdet]+1,0);
999 
1000  // Gains, pedestals (once !) and only for "noise" case
1001  if (((nevent1 == 1 && isubdet == 1) ||
1002  (nevent2 == 1 && isubdet == 2) ||
1003  (nevent3 == 1 && isubdet == 3) ||
1004  (nevent4 == 1 && isubdet == 4)) && noise_ == 1 && sub == isubdet) {
1005 
1006  HcalGenericDetId hcalGenDetId(digiItr->id());
1007  const HcalPedestal* pedestal = conditions->getPedestal(hcalGenDetId);
1008  const HcalGain* gain = conditions->getGain(hcalGenDetId);
1009  const HcalGainWidth* gainWidth = conditions->getGainWidth(hcalGenDetId);
1010  const HcalPedestalWidth* pedWidth = conditions-> getPedestalWidth(hcalGenDetId);
1011 
1012  for (int i = 0; i < 4; i++) {
1013  fill1D("HcalDigiTask_gain_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gain->getValue(i));
1014  fill1D("HcalDigiTask_gainWidth_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, gainWidth->getValue(i));
1015  fill1D("HcalDigiTask_pedestal_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, pedestal->getValue(i));
1016  fill1D("HcalDigiTask_pedestal_width_capId" + str(i) + "_Depth" + str(depth) + "_" + subdet_, pedWidth->getWidth(i));
1017  }
1018 
1019  fill2D("HcalDigiTask_gainMap_Depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi), gain->getValue(0));
1020  fill2D("HcalDigiTask_pwidthMap_Depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi), pedWidth->getWidth(0));
1021 
1022  }// end of event #1
1023  //std::cout << "==== End of event noise block in cell cycle" << std::endl;
1024 
1025  if (sub == isubdet) Ndig++; // subdet number of digi
1026 
1027  // No-noise case, only single subdet selected ===========================
1028 
1029  if (sub == isubdet && noise_ == 0) {
1030 
1031  HcalCalibrations calibrations = conditions->getHcalCalibrations(cell);
1032 
1033  const HcalQIECoder* channelCoder = conditions->getHcalCoder(cell);
1034  const HcalQIEShape* shape = conditions->getHcalShape(channelCoder);
1035  HcalCoderDb coder(*channelCoder, *shape);
1036  coder.adc2fC(dataFrame, tool);
1037 
1038  double noiseADC = (dataFrame)[0].adc();
1039  double noisefC = tool[0];
1040  // noise evaluations from "pre-samples"
1041  fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
1042  fill1D("HcalDigiTask_ADC0_fC_depth" + str(depth) + "_" + subdet_, noisefC);
1043 
1044 
1045  // OCCUPANCY maps fill
1046  fill2D("HcalDigiTask_ieta_iphi_occupancy_map_depth" + str(depth) + "_" + subdet_, double(ieta), double(iphi));
1047 
1048  fill1D("HcalDigiTask_depths_"+subdet_,double(depth));
1049 
1050  // Cycle on time slices
1051  // - for each Digi
1052  // - for one Digi with max SimHits E in subdet
1053 
1054 
1055  int closen = 0; // =1 if 1) seedSimHit = 1 and 2) the cell is the same
1056  if (ieta == ieta_Sim && iphi == iphi_Sim) closen = seedSimHit;
1057 
1058  for (int ii = 0; ii < tool.size(); ii++) {
1059  int capid = (dataFrame)[ii].capid();
1060  // single ts amplitude
1061  double val = (tool[ii] - calibrations.pedestal(capid));
1062 
1063  if (val > 100.) {
1064  fill1D("HcalDigiTask_ADC0_adc_depth" + str(depth) + "_" + subdet_, noiseADC);
1065  if(hep17_){
1066  if(!isHEP17){
1067  strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_" + subdet_;
1068  fill1D(strtmp, double(ii), val);
1069  } else {
1070  strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_HEP17";
1071  fill1D(strtmp, double(ii), val);
1072  }
1073  } else {
1074  strtmp = "HcalDigiTask_all_amplitudes_vs_bin_1D_depth" + str(depth) + "_" + subdet_;
1075  fill1D(strtmp, double(ii), val);
1076  }
1077  }
1078 
1079  if (closen == 1) {
1080  if(hep17_){
1081  if(!isHEP17){
1082  strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_" + subdet_;
1083  fill2D(strtmp, double(ii), val);
1084  }else{
1085  strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_HEP17";
1086  fill2D(strtmp, double(ii), val);
1087  }
1088  } else {
1089  strtmp = "HcalDigiTask_signal_amplitude_vs_bin_all_depths_" + subdet_;
1090  fill2D(strtmp, double(ii), val);
1091  }
1092  }
1093 
1094 
1095  // HB/HE/HO
1096  if (isubdet != 4 && ii >= 4 && ii <= 7) {
1097  v_ampl[0] += val;
1098  v_ampl[depth] += val;
1099 
1100  if (closen == 1) {
1101  v_ampl_c[0] += val;
1102  v_ampl_c[depth] += val;
1103  }
1104  }
1105 
1106  // HF
1107  if (isubdet == 4 && ii >= 2 && ii <= 4) {
1108  v_ampl[0] += val;
1109  v_ampl[depth] += val;
1110  if (closen == 1) {
1111  v_ampl_c[0] += val;
1112  v_ampl_c[depth] += val;
1113  }
1114  }
1115  }
1116  // end of time bucket sample
1117 
1118  // just 1D of all cells' amplitudes
1119  strtmp = "HcalDigiTask_sum_all_amplitudes_" + subdet_;
1120  fill1D(strtmp, v_ampl[0]);
1121 
1122  std::vector<int> v_ampl_sub(v_ampl.begin() + 1, v_ampl.end()); // remove element 0, which is the sum of any depth
1123  double ampl_max = *std::max_element(v_ampl_sub.begin(),v_ampl_sub.end());
1124  if (ampl_max>10.) indigis++;
1125  //KH if (ampl1 > 10. || ampl2 > 10. || ampl3 > 10. || ampl4 > 10.) indigis++;
1126 
1127  // fraction 5,6 bins if ampl. is big.
1128  if (v_ampl[1] > 30. && depth == 1 && closen == 1 && isubdet != 4) {
1129  double fBin5 = tool[4] - calibrations.pedestal((dataFrame)[4].capid());
1130  double fBin67 = tool[5] + tool[6]
1131  - calibrations.pedestal((dataFrame)[5].capid())
1132  - calibrations.pedestal((dataFrame)[6].capid());
1133 
1134  fBin5 /= v_ampl[1];
1135  fBin67 /= v_ampl[1];
1136 
1137  strtmp = "HcalDigiTask_bin_5_frac_" + subdet_;
1138  fill1D(strtmp, fBin5);
1139  strtmp = "HcalDigiTask_bin_6_7_frac_" + subdet_;
1140  fill1D(strtmp, fBin67);
1141 
1142  }
1143 
1144  //Special for HF
1145  if (isubdet == 4 && v_ampl[1] > 30. && depth == 1) {
1146  double fBin5 = tool[2] - calibrations.pedestal((dataFrame)[2].capid());
1147  double fBin67 = tool[3] + tool[4]
1148  - calibrations.pedestal((dataFrame)[3].capid())
1149  - calibrations.pedestal((dataFrame)[4].capid());
1150  fBin5 /= v_ampl[1];
1151  fBin67 /= v_ampl[1];
1152  strtmp = "HcalDigiTask_bin_5_frac_" + subdet_;
1153  fill1D(strtmp, fBin5);
1154  strtmp = "HcalDigiTask_bin_6_7_frac_" + subdet_;
1155  fill1D(strtmp, fBin67);
1156  }
1157 
1158 
1159  if(hep17_){
1160  if(!isHEP17){
1161  strtmp = "HcalDigiTask_signal_amplitude_" + subdet_;
1162  fill1D(strtmp, v_ampl[0]);
1163  strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_" + subdet_;
1164  fill1D(strtmp, v_ampl[depth]);
1165  } else {
1166  strtmp = "HcalDigiTask_signal_amplitude_HEP17";
1167  fill1D(strtmp, v_ampl[0]);
1168  strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_HEP17";
1169  fill1D(strtmp, v_ampl[depth]);
1170  }
1171  }else{
1172  strtmp = "HcalDigiTask_signal_amplitude_" + subdet_;
1173  fill1D(strtmp, v_ampl[0]);
1174  strtmp = "HcalDigiTask_signal_amplitude_depth" + str(depth) + "_" + subdet_;
1175  fill1D(strtmp, v_ampl[depth]);
1176  }
1177  }
1178  } // End of CYCLE OVER CELLS =============================================
1179 
1180  if (isubdet != 0 && noise_ == 0) { // signal only, once per event
1181  strtmp = "HcalDigiTask_number_of_amplitudes_above_10fC_" + subdet_;
1182  fill1D(strtmp, indigis);
1183 
1184  // SimHits once again !!!
1185  double eps = 1.e-3;
1186  std::vector<double> v_ehits(maxDepth_[isubdet]+1,0);
1187 
1188  if (mc_ == "yes") {
1190  iEvent.getByToken(tok_mc_, hcalHits);
1191  const edm::PCaloHitContainer * simhitResult = hcalHits.product();
1192  for (std::vector<PCaloHit>::const_iterator simhits = simhitResult->begin(); simhits != simhitResult->end(); ++simhits) {
1193 
1194  unsigned int id_ = simhits->id();
1195  int sub, depth, ieta, iphi;
1196  HcalDetId hid;
1198  else hid = HcalDetId(id_);
1199  sub = hid.subdet();
1200  depth = hid.depth();
1201  ieta = hid.ieta();
1202  iphi = hid.iphi();
1203 
1204  if(depth > maxDepth_[isubdet] && sub == isubdet){
1205  edm::LogWarning("HcalDetId") << "HcalDetID(SimHit) presents conflicting information. Depth: " << depth << ", iphi: " << iphi << ", ieta: " << ieta << ". Max depth from geometry is: " << maxDepth_[isubdet] << ". TestNumber = " << testNumber_;
1206  continue;
1207  }
1208 
1209 
1210  // take cell already found to be max energy in a particular subdet
1211  if (sub == isubdet && ieta == ieta_Sim && iphi == iphi_Sim) {
1212  double en = simhits->energy();
1213 
1214  v_ehits[0] += en;
1215  v_ehits[depth] += en;
1216 
1217  }
1218  } // simhit loop
1219 
1220  strtmp = "HcalDigiTask_amplitude_vs_simhits_" + subdet_;
1221  if (v_ehits[0] > eps) fill2D(strtmp, v_ehits[0], v_ampl_c[0]);
1222  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
1223  strtmp = "HcalDigiTask_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
1224  if (v_ehits[depth] > eps) fill2D(strtmp, v_ehits[depth], v_ampl_c[depth]);
1225  }
1226 
1227  strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_" + subdet_;
1228  if (v_ehits[0] > eps) fillPf(strtmp, v_ehits[0], v_ampl_c[0]);
1229  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
1230  strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
1231  if (v_ehits[depth] > eps) fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
1232  }
1233 
1234  strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_" + subdet_;
1235  if (v_ehits[0] > eps) fill1D(strtmp, v_ampl_c[0] / v_ehits[0]);
1236  for (int depth = 1; depth <= maxDepth_[isubdet]; depth++) {
1237  strtmp = "HcalDigiTask_amplitude_vs_simhits_profile_depth" + str(depth) + "_" + subdet_;
1238  if (v_ehits[depth] > eps) fillPf(strtmp, v_ehits[depth], v_ampl_c[depth]);
1239  strtmp = "HcalDigiTask_ratio_amplitude_vs_simhits_depth" + str(depth) + "_" + subdet_;
1240  if (v_ehits[depth] > eps) fill1D(strtmp, v_ampl_c[depth] / v_ehits[depth]);
1241  }
1242 
1243 
1244  } // end of if(mc_ == "yes")
1245 
1246  strtmp = "HcalDigiTask_Ndigis_" + subdet_;
1247  fill1D(strtmp, double(Ndig));
1248 
1249  } // end of if( subdet != 0 && noise_ == 0) { // signal only
1250 
1251 } //HcalDataFrameContainer
1252 
1254  if (!msm_->count(name)) (*msm_)[name] = ib.book1D(name.c_str(), name.c_str(), n, min, max);
1255 }
1256 
1258  if (!msm_->count(name)) (*msm_)[name] = ib.book1D(name.c_str(), name.c_str(), limX.n, limX.min, limX.max);
1259 }
1260 
1262  msm_->find(name)->second->Fill(X, weight);
1263 }
1264 
1266  if (!msm_->count(name)) (*msm_)[name] = ib.book2D(name.c_str(), name.c_str(), limX.n, limX.min, limX.max, limY.n, limY.min, limY.max);
1267 }
1268 
1269 void HcalDigisValidation::fill2D(std::string name, double X, double Y, double weight) {
1270  msm_->find(name)->second->Fill(X, Y, weight);
1271 }
1272 
1274  if (!msm_->count(name)) (*msm_)[name] = ib.bookProfile(name.c_str(), name.c_str(), limX.n, limX.min, limX.max, limY.n, limY.min, limY.max);
1275 }
1276 
1277 void HcalDigisValidation::bookPf(DQMStore::IBooker &ib, std::string name, const HistLim& limX, const HistLim& limY, const char *option) {
1278  if (!msm_->count(name)) (*msm_)[name] = ib.bookProfile(name.c_str(), name.c_str(), limX.n, limX.min, limX.max, limY.n, limY.min, limY.max, option);
1279 }
1280 
1282  msm_->find(name)->second->Fill(X, Y);
1283 }
1284 
1286  if (!msm_->count(name)) return NULL;
1287  else return msm_->find(name)->second;
1288 }
1289 
1291  std::stringstream out;
1292  out << x;
1293  return out.str();
1294 }
1295 
1296 
1297 //define this as a plug-in
1299 
1300 
int adc(sample_type sample)
get the ADC sample (12 bits)
virtual void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &)
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:45
const HcalGainWidth * getGainWidth(const HcalGenericDetId &fId) const
T getUntrackedParameter(std::string const &, T const &) const
std::vector< PCaloHit > PCaloHitContainer
boost::transform_iterator< IterHelp, boost::counting_iterator< int > > const_iterator
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:49
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
unsigned int getHxSize(const int type) const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
const HcalTopology * htopology
std::vector< HcalTriggerPrimitiveDigi >::const_iterator const_iterator
bool exists(std::string const &parameterName) const
checks if a parameter exists
#define X(str)
Definition: MuonsGrabber.cc:48
#define NULL
Definition: scimark2.h:8
Definition: weight.py:1
void reco(const edm::Event &iEvent, const edm::EventSetup &iSetup, const edm::EDGetTokenT< edm::SortedCollection< Digi > > &tok)
double pedestal(int fCapId) const
get pedestal for capid=0..3
const_iterator begin() const
float getValue(int fCapId) const
get value for capId = 0..3
Definition: HcalGain.h:22
virtual void analyze(const edm::Event &, const edm::EventSetup &)
float getValue(int fCapId) const
get value for capId = 0..3
Definition: HcalGainWidth.h:21
edm::EDGetTokenT< QIE11DigiCollection > tok_qie11_hbhe_
edm::EDGetTokenT< QIE10DigiCollection > tok_qie10_hf_
std::map< std::string, MonitorElement * > * msm_
int getMaxDepth(const int type) const
int depth() const
get the tower depth
Definition: HcalDetId.cc:108
int iEvent
Definition: GenABIO.cc:230
edm::EDGetTokenT< HBHEDigiCollection > tok_hbhe_
MonitorElement * monitor(std::string name)
virtual void adc2fC(const HBHEDataFrame &df, CaloSamples &lf) const
Definition: HcalCoderDb.cc:68
void fillPf(std::string name, double X, double Y)
void fill2D(std::string name, double X, double Y, double weight=1)
edm::ESHandle< HcalTopology > htopo
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
HcalSubdetector
Definition: HcalAssistant.h:31
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
HcalDigisValidation(const edm::ParameterSet &)
void bookPf(DQMStore::IBooker &ib, std::string name, const HistLim &limX, const HistLim &limY)
T min(T a, T b)
Definition: MathUtil.h:58
bool isValid() const
Definition: HandleBase.h:74
const HcalDDDRecConstants * hcons
ii
Definition: cuy.py:588
virtual void dqmBeginRun(const edm::Run &run, const edm::EventSetup &c)
edm::InputTag QIE10inputTag_
void book2D(DQMStore::IBooker &ib, std::string name, const HistLim &limX, const HistLim &limY)
const_iterator end() const
int iphi() const
get the cell iphi
Definition: HcalDetId.cc:103
edm::EDGetTokenT< HFDigiCollection > tok_hf_
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:277
T const * product() const
Definition: Handle.h:81
void fill1D(std::string name, double X, double weight=1)
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
int size() const
get the size
Definition: CaloSamples.h:24
edm::EDGetTokenT< edm::PCaloHitContainer > tok_mc_
edm::EDGetTokenT< HODigiCollection > tok_ho_
const T & get() const
Definition: EventSetup.h:56
void booking(DQMStore::IBooker &ib, std::string subdetopt, int bnoise, int bmc)
const HcalGain * getGain(const HcalGenericDetId &fId) const
float getWidth(int fCapId) const
get width (sqrt(sigma_i_i)) for capId = 0..3
const HcalQIECoder * getHcalCoder(const HcalGenericDetId &fId) const
edm::EDGetTokenT< HcalTrigPrimDigiCollection > tok_emulTPs_
edm::InputTag QIE11inputTag_
const HcalQIEShape * getHcalShape(const HcalGenericDetId &fId) const
std::string str(int x)
const_iterator end() const
HLT enums.
edm::EDGetTokenT< HcalTrigPrimDigiCollection > tok_dataTPs_
edm::ESHandle< HcalDbService > conditions
virtual double hcaletValue(const int &ieta, const int &iphi, const int &version, const int &compressedValue) const =0
void book1D(DQMStore::IBooker &ib, std::string name, int n, double min, double max)
DetId relabel(const uint32_t testId) const
const HcalCalibrations & getHcalCalibrations(const HcalGenericDetId &fId) const
const HcalPedestal * getPedestal(const HcalGenericDetId &fId) const
const_iterator begin() const
Definition: Run.h:42
ib
Definition: cuy.py:660
int firstHFTower(int version) const