CMS 3D CMS Logo

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