CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HcalSimHitCheck.cc
Go to the documentation of this file.
1 #include "CommonTools/UtilAlgos/interface/TFileService.h"
2 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
3 
4 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
5 #include "FWCore/Framework/interface/Event.h"
6 #include "FWCore/Framework/interface/EventSetup.h"
7 #include "FWCore/Framework/interface/Frameworkfwd.h"
8 
9 #include "FWCore/MessageLogger/interface/MessageLogger.h"
10 #include "FWCore/ParameterSet/interface/ParameterSet.h"
11 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
12 #include "FWCore/Utilities/interface/Exception.h"
13 #include "FWCore/Utilities/interface/InputTag.h"
14 
15 #include "FWCore/ServiceRegistry/interface/Service.h"
16 
17 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
18 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
19 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
20 
21 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
22 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
23 
24 #include "TH1D.h"
25 #include "TH2D.h"
26 #include "TProfile.h"
27 
28 #include <fstream>
29 #include <iostream>
30 #include <map>
31 #include <string>
32 #include <vector>
33 
34 class HcalSimHitCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
35 public:
36  HcalSimHitCheck(const edm::ParameterSet &ps);
37 
38  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
39 
40 protected:
41  void beginRun(edm::Run const &, edm::EventSetup const &) override;
42  void analyze(edm::Event const &, edm::EventSetup const &) override;
43  void endRun(edm::Run const &, edm::EventSetup const &) override {}
44 
45  void analyzeHits(std::vector<PCaloHit> &);
46 
47 private:
48  const HcalDDDRecConstants *hcons_;
49  int maxDepthHB_, maxDepthHE_;
50  int maxDepthHO_, maxDepthHF_;
51  int maxDepth_;
52  edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
53 
54  int iphi_bins;
55  float iphi_min, iphi_max;
56  int ieta_bins_HB;
57  float ieta_min_HB, ieta_max_HB;
58  int ieta_bins_HE;
59  float ieta_min_HE, ieta_max_HE;
60  int ieta_bins_HO;
61  float ieta_min_HO, ieta_max_HO;
62  int ieta_bins_HF;
63  float ieta_min_HF, ieta_max_HF;
64 
65  std::string g4Label, hcalHits, outFile_;
66  int verbose_;
67  bool checkHit_, testNumber_, hep17_;
68 
69  edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
70 
71  TH1D *meAllNHit_, *meBadDetHit_, *meBadSubHit_, *meBadIdHit_;
72  TH1D *meHBNHit_, *meHENHit_, *meHONHit_, *meHFNHit_;
73  TH1D *meDetectHit_, *meSubdetHit_, *meDepthHit_, *meEtaHit_;
74  TH2D *meEtaPhiHit_;
75  std::vector<TH2D *> meEtaPhiHitDepth_;
76  TH1D *mePhiHit_, *mePhiHitb_, *meEnergyHit_, *meTimeHit_, *meTimeWHit_;
77  TH1D *meHBDepHit_, *meHEDepHit_, *meHODepHit_, *meHFDepHit_, *meHFDepHitw_;
78  TH1D *meHBEtaHit_, *meHEEtaHit_, *meHOEtaHit_, *meHFEtaHit_;
79  TH1D *meHBPhiHit_, *meHEPhiHit_, *meHOPhiHit_, *meHFPhiHit_;
80  TH1D *meHBEneHit_, *meHEEneHit_, *meHOEneHit_, *meHFEneHit_;
81  TH2D *meHBEneMap_, *meHEEneMap_, *meHOEneMap_, *meHFEneMap_;
82  TH1D *meHBEneSum_, *meHEEneSum_, *meHOEneSum_, *meHFEneSum_;
83  TProfile *meHBEneSum_vs_ieta_, *meHEEneSum_vs_ieta_, *meHOEneSum_vs_ieta_, *meHFEneSum_vs_ieta_;
84  TH1D *meHBTimHit_, *meHETimHit_, *meHOTimHit_, *meHFTimHit_;
85  TH1D *meHBEneHit2_, *meHEEneHit2_, *meHOEneHit2_, *meHFEneHit2_;
86  TH1D *meHBL10Ene_, *meHEL10Ene_, *meHOL10Ene_, *meHFL10Ene_;
87  TProfile *meHBL10EneP_, *meHEL10EneP_, *meHOL10EneP_, *meHFL10EneP_;
88  TH1D *meHEP17EneHit_, *meHEP17EneHit2_;
89 };
90 
91 HcalSimHitCheck::HcalSimHitCheck(const edm::ParameterSet &ps) {
92  g4Label = ps.getParameter<std::string>("moduleLabel");
93  hcalHits = ps.getParameter<std::string>("HitCollection");
94  outFile_ = ps.getParameter<std::string>("outputFile");
95  verbose_ = ps.getParameter<int>("Verbose");
96  testNumber_ = ps.getParameter<bool>("TestNumber");
97  hep17_ = ps.getParameter<bool>("hep17");
98  checkHit_ = true;
99 
100  tok_hits_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hcalHits));
101  tok_HRNDC_ = esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>();
102 
103  edm::LogVerbatim("HcalSim") << "Module Label: " << g4Label << " Hits: " << hcalHits << " / " << checkHit_
104  << " Output: " << outFile_;
105 }
106 
107 void HcalSimHitCheck::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
108  edm::ParameterSetDescription desc;
109  desc.add<std::string>("moduleLabel", "g4SimHits");
110  desc.add<std::string>("HitCollection", "HcalHits");
111  desc.add<std::string>("outputFile", "hcHit.root");
112  desc.add<int>("Verbose", 0);
113  desc.add<bool>("TestNumber", true);
114  desc.add<bool>("hep17", false);
115  descriptions.add("hcalSimHitCheck", desc);
116 }
117 
118 void HcalSimHitCheck::beginRun(const edm::Run &, const edm::EventSetup &es) {
119  hcons_ = &es.getData(tok_HRNDC_);
120  maxDepthHB_ = hcons_->getMaxDepth(0);
121  maxDepthHE_ = hcons_->getMaxDepth(1);
122  maxDepthHF_ = hcons_->getMaxDepth(2);
123  maxDepthHO_ = hcons_->getMaxDepth(3);
124  maxDepth_ = (maxDepthHB_ > maxDepthHE_ ? maxDepthHB_ : maxDepthHE_);
125  maxDepth_ = (maxDepth_ > maxDepthHF_ ? maxDepth_ : maxDepthHF_);
126  maxDepth_ = (maxDepth_ > maxDepthHO_ ? maxDepth_ : maxDepthHO_);
127 
128  // Get Phi segmentation from geometry, use the max phi number so that all iphi
129  // values are included.
130 
131  int NphiMax = hcons_->getNPhi(0);
132 
133  NphiMax = (hcons_->getNPhi(1) > NphiMax ? hcons_->getNPhi(1) : NphiMax);
134  NphiMax = (hcons_->getNPhi(2) > NphiMax ? hcons_->getNPhi(2) : NphiMax);
135  NphiMax = (hcons_->getNPhi(3) > NphiMax ? hcons_->getNPhi(3) : NphiMax);
136 
137  // Center the iphi bins on the integers
138  iphi_min = 0.5;
139  iphi_max = NphiMax + 0.5;
140  iphi_bins = (int)(iphi_max - iphi_min);
141 
142  int iEtaHBMax = hcons_->getEtaRange(0).second;
143  int iEtaHEMax = std::max(hcons_->getEtaRange(1).second, 1);
144  int iEtaHFMax = hcons_->getEtaRange(2).second;
145  int iEtaHOMax = hcons_->getEtaRange(3).second;
146 
147  // Retain classic behavior, all plots have same ieta range.
148  // Comment out code to allow each subdetector to have its on range
149 
150  int iEtaMax = (iEtaHBMax > iEtaHEMax ? iEtaHBMax : iEtaHEMax);
151  iEtaMax = (iEtaMax > iEtaHFMax ? iEtaMax : iEtaHFMax);
152  iEtaMax = (iEtaMax > iEtaHOMax ? iEtaMax : iEtaHOMax);
153 
154  iEtaHBMax = iEtaMax;
155  iEtaHEMax = iEtaMax;
156  iEtaHFMax = iEtaMax;
157  iEtaHOMax = iEtaMax;
158 
159  // Give an empty bin around the subdet ieta range to make it clear that all
160  // ieta rings have been included
161  ieta_min_HB = -iEtaHBMax - 1.5;
162  ieta_max_HB = iEtaHBMax + 1.5;
163  ieta_bins_HB = (int)(ieta_max_HB - ieta_min_HB);
164 
165  ieta_min_HE = -iEtaHEMax - 1.5;
166  ieta_max_HE = iEtaHEMax + 1.5;
167  ieta_bins_HE = (int)(ieta_max_HE - ieta_min_HE);
168 
169  ieta_min_HF = -iEtaHFMax - 1.5;
170  ieta_max_HF = iEtaHFMax + 1.5;
171  ieta_bins_HF = (int)(ieta_max_HF - ieta_min_HF);
172 
173  ieta_min_HO = -iEtaHOMax - 1.5;
174  ieta_max_HO = iEtaHOMax + 1.5;
175  ieta_bins_HO = (int)(ieta_max_HO - ieta_min_HO);
176 
177  Char_t hname[100];
178  Char_t htitle[100];
179 
180  edm::Service<TFileService> fs;
181 
182  // Histograms for Hits
183  if (checkHit_) {
184  meAllNHit_ = fs->make<TH1D>("Hit01", "Number of Hits in HCal", 20000, 0., 20000.);
185  meBadDetHit_ = fs->make<TH1D>("Hit02", "Hits with wrong Det", 100, 0., 100.);
186  meBadSubHit_ = fs->make<TH1D>("Hit03", "Hits with wrong Subdet", 100, 0., 100.);
187  meBadIdHit_ = fs->make<TH1D>("Hit04", "Hits with wrong ID", 100, 0., 100.);
188  meHBNHit_ = fs->make<TH1D>("Hit05", "Number of Hits in HB", 20000, 0., 20000.);
189  meHENHit_ = fs->make<TH1D>("Hit06", "Number of Hits in HE", 10000, 0., 10000.);
190  meHONHit_ = fs->make<TH1D>("Hit07", "Number of Hits in HO", 10000, 0., 10000.);
191  meHFNHit_ = fs->make<TH1D>("Hit08", "Number of Hits in HF", 10000, 0., 10000.);
192  meHBNHit_->Sumw2();
193  meHENHit_->Sumw2();
194  meHONHit_->Sumw2();
195  meHFNHit_->Sumw2();
196  meDetectHit_ = fs->make<TH1D>("Hit09", "Detector ID", 50, 0., 50.);
197  meSubdetHit_ = fs->make<TH1D>("Hit10", "Subdetectors in HCal", 50, 0., 50.);
198  meDepthHit_ = fs->make<TH1D>("Hit11", "Depths in HCal", 20, 0., 20.);
199  meEtaHit_ = fs->make<TH1D>("Hit12", "Eta in HCal", ieta_bins_HF, ieta_min_HF, ieta_max_HF);
200  meEtaPhiHit_ = fs->make<TH2D>(
201  "Hit12b", "Eta-phi in HCal", ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max);
202  for (int depth = 1; depth <= maxDepth_; depth++) {
203  sprintf(hname, "Hit12bd%d", depth);
204  sprintf(htitle, "Eta-phi in HCal d%d", depth);
205  meEtaPhiHitDepth_.emplace_back(
206  fs->make<TH2D>(hname, htitle, ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max));
207  }
208  // KC: There are different phi segmentation schemes, this plot uses wider
209  // bins to represent the most sparse segmentation
210  mePhiHit_ = fs->make<TH1D>("Hit13", "Phi in HCal (HB,HO)", iphi_bins, iphi_min, iphi_max);
211  mePhiHitb_ = fs->make<TH1D>("Hit13b", "Phi in HCal (HE,HF)", iphi_bins, iphi_min, iphi_max);
212  meEnergyHit_ = fs->make<TH1D>("Hit14", "Energy in HCal", 2000, 0., 20.);
213  meTimeHit_ = fs->make<TH1D>("Hit15", "Time in HCal", 528, 0., 528.);
214  meTimeWHit_ = fs->make<TH1D>("Hit16", "Time in HCal (E wtd)", 528, 0., 528.);
215  meHBDepHit_ = fs->make<TH1D>("Hit17", "Depths in HB", 20, 0., 20.);
216  meHEDepHit_ = fs->make<TH1D>("Hit18", "Depths in HE", 20, 0., 20.);
217  meHODepHit_ = fs->make<TH1D>("Hit19", "Depths in HO", 20, 0., 20.);
218  meHFDepHit_ = fs->make<TH1D>("Hit20", "Depths in HF", 20, 0., 20.);
219  meHBDepHit_->Sumw2();
220  meHEDepHit_->Sumw2();
221  meHODepHit_->Sumw2();
222  meHFDepHit_->Sumw2();
223  meHFDepHitw_ = fs->make<TH1D>("Hit20b", "Depths in HF (p.e. weighted)", 20, 0., 20.);
224  meHBEtaHit_ = fs->make<TH1D>("Hit21", "Eta in HB", ieta_bins_HB, ieta_min_HB, ieta_max_HB);
225  meHEEtaHit_ = fs->make<TH1D>("Hit22", "Eta in HE", ieta_bins_HE, ieta_min_HE, ieta_max_HE);
226  meHOEtaHit_ = fs->make<TH1D>("Hit23", "Eta in HO", ieta_bins_HO, ieta_min_HO, ieta_max_HO);
227  meHFEtaHit_ = fs->make<TH1D>("Hit24", "Eta in HF", ieta_bins_HF, ieta_min_HF, ieta_max_HF);
228  meHBEtaHit_->Sumw2();
229  meHEEtaHit_->Sumw2();
230  meHOEtaHit_->Sumw2();
231  meHFEtaHit_->Sumw2();
232  meHBPhiHit_ = fs->make<TH1D>("Hit25", "Phi in HB", iphi_bins, iphi_min, iphi_max);
233  meHEPhiHit_ = fs->make<TH1D>("Hit26", "Phi in HE", iphi_bins, iphi_min, iphi_max);
234  meHOPhiHit_ = fs->make<TH1D>("Hit27", "Phi in HO", iphi_bins, iphi_min, iphi_max);
235  meHFPhiHit_ = fs->make<TH1D>("Hit28", "Phi in HF", iphi_bins, iphi_min, iphi_max);
236  meHBPhiHit_->Sumw2();
237  meHEPhiHit_->Sumw2();
238  meHOPhiHit_->Sumw2();
239  meHFPhiHit_->Sumw2();
240  meHBEneHit_ = fs->make<TH1D>("Hit29", "Energy in HB", 2000, 0., 20.);
241  meHEEneHit_ = fs->make<TH1D>("Hit30", "Energy in HE", 500, 0., 5.);
242  meHEP17EneHit_ = fs->make<TH1D>("Hit30b", "Energy in HEP17", 500, 0., 5.);
243  meHOEneHit_ = fs->make<TH1D>("Hit31", "Energy in HO", 500, 0., 5.);
244  meHFEneHit_ = fs->make<TH1D>("Hit32", "Energy in HF", 1001, -0.5, 1000.5);
245  meHBEneHit_->Sumw2();
246  meHEEneHit_->Sumw2();
247  meHOEneHit_->Sumw2();
248  meHFEneHit_->Sumw2();
249 
250  // HxEneMap, HxEneSum, HxEneSum_vs_ieta plot the sum of the simhits energy
251  // within a single ieta-iphi tower.
252 
253  meHBEneMap_ =
254  fs->make<TH2D>("HBEneMap", "HBEneMap", ieta_bins_HB, ieta_min_HB, ieta_max_HB, iphi_bins, iphi_min, iphi_max);
255  meHEEneMap_ =
256  fs->make<TH2D>("HEEneMap", "HEEneMap", ieta_bins_HE, ieta_min_HE, ieta_max_HE, iphi_bins, iphi_min, iphi_max);
257  meHOEneMap_ =
258  fs->make<TH2D>("HOEneMap", "HOEneMap", ieta_bins_HO, ieta_min_HO, ieta_max_HO, iphi_bins, iphi_min, iphi_max);
259  meHFEneMap_ =
260  fs->make<TH2D>("HFEneMap", "HFEneMap", ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max);
261 
262  meHBEneSum_ = fs->make<TH1D>("HBEneSum", "HBEneSum", 2000, 0., 20.);
263  meHEEneSum_ = fs->make<TH1D>("HEEneSum", "HEEneSum", 500, 0., 5.);
264  meHOEneSum_ = fs->make<TH1D>("HOEneSum", "HOEneSum", 500, 0., 5.);
265  meHFEneSum_ = fs->make<TH1D>("HFEneSum", "HFEneSum", 1001, -0.5, 1000.5);
266  meHBEneSum_->Sumw2();
267  meHEEneSum_->Sumw2();
268  meHOEneSum_->Sumw2();
269  meHFEneSum_->Sumw2();
270 
271  meHBEneSum_vs_ieta_ = fs->make<TProfile>(
272  "HBEneSum_vs_ieta", "HBEneSum_vs_ieta", ieta_bins_HB, ieta_min_HB, ieta_max_HB, -10.5, 2000.5);
273  meHEEneSum_vs_ieta_ = fs->make<TProfile>(
274  "HEEneSum_vs_ieta", "HEEneSum_vs_ieta", ieta_bins_HE, ieta_min_HE, ieta_max_HE, -10.5, 2000.5);
275  meHOEneSum_vs_ieta_ = fs->make<TProfile>(
276  "HOEneSum_vs_ieta", "HOEneSum_vs_ieta", ieta_bins_HO, ieta_min_HO, ieta_max_HO, -10.5, 2000.5);
277  meHFEneSum_vs_ieta_ = fs->make<TProfile>(
278  "HFEneSum_vs_ieta", "HFEneSum_vs_ieta", ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10.5, 2000.5);
279 
280  meHBTimHit_ = fs->make<TH1D>("Hit33", "Time in HB", 528, 0., 528.);
281  meHETimHit_ = fs->make<TH1D>("Hit34", "Time in HE", 528, 0., 528.);
282  meHOTimHit_ = fs->make<TH1D>("Hit35", "Time in HO", 528, 0., 528.);
283  meHFTimHit_ = fs->make<TH1D>("Hit36", "Time in HF", 528, 0., 528.);
284  meHBTimHit_->Sumw2();
285  meHETimHit_->Sumw2();
286  meHOTimHit_->Sumw2();
287  meHFTimHit_->Sumw2();
288  // These are the zoomed in energy ranges
289  meHBEneHit2_ = fs->make<TH1D>("Hit37", "Energy in HB 2", 100, 0., 0.0001);
290  meHEEneHit2_ = fs->make<TH1D>("Hit38", "Energy in HE 2", 100, 0., 0.0001);
291  meHEP17EneHit2_ = fs->make<TH1D>("Hit38b", "Energy in HEP17 2", 100, 0., 0.0001);
292  meHOEneHit2_ = fs->make<TH1D>("Hit39", "Energy in HO 2", 100, 0., 0.0001);
293  meHFEneHit2_ = fs->make<TH1D>("Hit40", "Energy in HF 2", 100, 0.5, 100.5);
294  meHBL10Ene_ = fs->make<TH1D>("Hit41", "Log10Energy in HB", 140, -10., 4.);
295  meHEL10Ene_ = fs->make<TH1D>("Hit42", "Log10Energy in HE", 140, -10., 4.);
296  meHFL10Ene_ = fs->make<TH1D>("Hit43", "Log10Energy in HF", 50, -1., 4.);
297  meHOL10Ene_ = fs->make<TH1D>("Hit44", "Log10Energy in HO", 140, -10., 4.);
298  meHBL10EneP_ = fs->make<TProfile>("Hit45", "Log10Energy in HB vs Hit contribution", 140, -10., 4., 0., 1.);
299  meHEL10EneP_ = fs->make<TProfile>("Hit46", "Log10Energy in HE vs Hit contribution", 140, -10., 4., 0., 1.);
300  meHFL10EneP_ = fs->make<TProfile>("Hit47", "Log10Energy in HF vs Hit contribution", 140, -10., 4., 0., 1.);
301  meHOL10EneP_ = fs->make<TProfile>("Hit48", "Log10Energy in HO vs Hit contribution", 140, -10., 4., 0., 1.);
302  }
303 }
304 
305 void HcalSimHitCheck::analyze(const edm::Event &e, const edm::EventSetup &) {
306  if (verbose_ > 0)
307  edm::LogVerbatim("HcalSim") << "Run = " << e.id().run() << " Event = " << e.id().event();
308 
309  std::vector<PCaloHit> caloHits;
310  edm::Handle<edm::PCaloHitContainer> hitsHcal;
311 
312  bool getHits = false;
313  if (checkHit_) {
314  e.getByToken(tok_hits_, hitsHcal);
315  if (hitsHcal.isValid())
316  getHits = true;
317  }
318  if (verbose_ > 0)
319  edm::LogVerbatim("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
320 
321  if (getHits) {
322  caloHits.insert(caloHits.end(), hitsHcal->begin(), hitsHcal->end());
323  if (verbose_ > 0)
324  edm::LogVerbatim("HcalSim") << "HcalValidation: Hit buffer " << caloHits.size();
325  analyzeHits(caloHits);
326  }
327 }
328 
329 void HcalSimHitCheck::analyzeHits(std::vector<PCaloHit> &hits) {
330  int nHit = hits.size();
331  int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nBad1 = 0, nBad2 = 0, nBad = 0;
332  std::vector<double> encontHB(140, 0.);
333  std::vector<double> encontHE(140, 0.);
334  std::vector<double> encontHF(140, 0.);
335  std::vector<double> encontHO(140, 0.);
336  double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0;
337 
338  double HBEneMap[ieta_bins_HB][iphi_bins];
339  double HEEneMap[ieta_bins_HE][iphi_bins];
340  double HOEneMap[ieta_bins_HO][iphi_bins];
341  double HFEneMap[ieta_bins_HF][iphi_bins];
342 
343  // Works in ieta_min_Hx is < 0
344  int eta_offset_HB = -(int)ieta_min_HB;
345  int eta_offset_HE = -(int)ieta_min_HE;
346  int eta_offset_HO = -(int)ieta_min_HO;
347  int eta_offset_HF = -(int)ieta_min_HF;
348 
349  for (int i = 0; i < ieta_bins_HB; i++) {
350  for (int j = 0; j < iphi_bins; j++) {
351  HBEneMap[i][j] = 0.;
352  }
353  }
354 
355  for (int i = 0; i < ieta_bins_HE; i++) {
356  for (int j = 0; j < iphi_bins; j++) {
357  HEEneMap[i][j] = 0.;
358  }
359  }
360 
361  for (int i = 0; i < ieta_bins_HO; i++) {
362  for (int j = 0; j < iphi_bins; j++) {
363  HOEneMap[i][j] = 0.;
364  }
365  }
366 
367  for (int i = 0; i < ieta_bins_HF; i++) {
368  for (int j = 0; j < iphi_bins; j++) {
369  HFEneMap[i][j] = 0.;
370  }
371  }
372 
373  for (int i = 0; i < nHit; i++) {
374  double energy = hits[i].energy();
375  double log10en = log10(energy);
376  int log10i = int((log10en + 10.) * 10.);
377  double time = hits[i].time();
378  unsigned int id = hits[i].id();
379  int det, subdet, depth, eta, phi;
380  HcalDetId hid;
381  if (testNumber_)
382  hid = HcalHitRelabeller::relabel(id, hcons_);
383  else
384  hid = HcalDetId(id);
385  det = hid.det();
386  subdet = hid.subdet();
387  depth = hid.depth();
388  eta = hid.ieta();
389  phi = hid.iphi();
390 
391  if (verbose_ > 1)
392  edm::LogVerbatim("HcalSim") << "Hit[" << i << "] ID " << std::hex << id << std::dec << " Det " << det << " Sub "
393  << subdet << " depth " << depth << " Eta " << eta << " Phi " << phi << " E " << energy
394  << " time " << time;
395  if (det == 4) { // Check DetId.h
396  if (subdet == static_cast<int>(HcalBarrel))
397  nHB++;
398  else if (subdet == static_cast<int>(HcalEndcap))
399  nHE++;
400  else if (subdet == static_cast<int>(HcalOuter))
401  nHO++;
402  else if (subdet == static_cast<int>(HcalForward))
403  nHF++;
404  else {
405  nBad++;
406  nBad2++;
407  }
408  } else {
409  nBad++;
410  nBad1++;
411  }
412 
413  meDetectHit_->Fill(double(det));
414  if (det == 4) {
415  meSubdetHit_->Fill(double(subdet));
416  meDepthHit_->Fill(double(depth));
417  meEtaHit_->Fill(double(eta));
418  meEtaPhiHit_->Fill(double(eta), double(phi));
419  meEtaPhiHitDepth_[depth - 1]->Fill(double(eta), double(phi));
420 
421  // We will group the phi plots by HB,HO and HE,HF since these groups share
422  // similar segmentation schemes
423  if (subdet == static_cast<int>(HcalBarrel))
424  mePhiHit_->Fill(double(phi));
425  else if (subdet == static_cast<int>(HcalEndcap))
426  mePhiHitb_->Fill(double(phi));
427  else if (subdet == static_cast<int>(HcalOuter))
428  mePhiHit_->Fill(double(phi));
429  else if (subdet == static_cast<int>(HcalForward))
430  mePhiHitb_->Fill(double(phi));
431 
432  // KC: HF energy is in photoelectrons rather than eV, so it will not be
433  // included in total HCal energy
434  if (subdet != static_cast<int>(HcalForward)) {
435  meEnergyHit_->Fill(energy);
436 
437  // Since the HF energy is a different scale it does not make sense to
438  // include it in the Energy Weighted Plot
439  meTimeWHit_->Fill(double(time), energy);
440  }
441  meTimeHit_->Fill(time);
442 
443  if (subdet == static_cast<int>(HcalBarrel)) {
444  meHBDepHit_->Fill(double(depth));
445  meHBEtaHit_->Fill(double(eta));
446  meHBPhiHit_->Fill(double(phi));
447  meHBEneHit_->Fill(energy);
448  meHBEneHit2_->Fill(energy);
449  meHBTimHit_->Fill(time);
450  meHBL10Ene_->Fill(log10en);
451  if (log10i >= 0 && log10i < 140)
452  encontHB[log10i] += energy;
453  entotHB += energy;
454 
455  HBEneMap[eta + eta_offset_HB][phi - 1] += energy;
456 
457  } else if (subdet == static_cast<int>(HcalEndcap)) {
458  meHEDepHit_->Fill(double(depth));
459  meHEEtaHit_->Fill(double(eta));
460  meHEPhiHit_->Fill(double(phi));
461 
462  bool isHEP17 = (phi >= 63) && (phi <= 66) && (eta > 0);
463  if (hep17_) {
464  if (!isHEP17) {
465  meHEEneHit_->Fill(energy);
466  meHEEneHit2_->Fill(energy);
467  } else {
468  meHEP17EneHit_->Fill(energy);
469  meHEP17EneHit2_->Fill(energy);
470  }
471  } else {
472  meHEEneHit_->Fill(energy);
473  meHEEneHit2_->Fill(energy);
474  }
475 
476  meHETimHit_->Fill(time);
477  meHEL10Ene_->Fill(log10en);
478  if (log10i >= 0 && log10i < 140)
479  encontHE[log10i] += energy;
480  entotHE += energy;
481 
482  HEEneMap[eta + eta_offset_HE][phi - 1] += energy;
483 
484  } else if (subdet == static_cast<int>(HcalOuter)) {
485  meHODepHit_->Fill(double(depth));
486  meHOEtaHit_->Fill(double(eta));
487  meHOPhiHit_->Fill(double(phi));
488  meHOEneHit_->Fill(energy);
489  meHOEneHit2_->Fill(energy);
490  meHOTimHit_->Fill(time);
491  meHOL10Ene_->Fill(log10en);
492  if (log10i >= 0 && log10i < 140)
493  encontHO[log10i] += energy;
494  entotHO += energy;
495 
496  HOEneMap[eta + eta_offset_HO][phi - 1] += energy;
497 
498  } else if (subdet == static_cast<int>(HcalForward)) {
499  meHFDepHit_->Fill(double(depth));
500  meHFDepHitw_->Fill(double(depth), energy);
501  meHFEtaHit_->Fill(double(eta));
502  meHFPhiHit_->Fill(double(phi));
503  meHFEneHit_->Fill(energy);
504  meHFEneHit2_->Fill(energy);
505  meHFTimHit_->Fill(time);
506  meHFL10Ene_->Fill(log10en);
507  if (log10i >= 0 && log10i < 140)
508  encontHF[log10i] += energy;
509  entotHF += energy;
510 
511  HFEneMap[eta + eta_offset_HF][phi - 1] += energy;
512  }
513  }
514  }
515  if (entotHB != 0)
516  for (int i = 0; i < 140; i++)
517  meHBL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHB[i] / entotHB);
518  if (entotHE != 0)
519  for (int i = 0; i < 140; i++)
520  meHEL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHE[i] / entotHE);
521  if (entotHF != 0)
522  for (int i = 0; i < 140; i++)
523  meHFL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHF[i] / entotHF);
524  if (entotHO != 0)
525  for (int i = 0; i < 140; i++)
526  meHOL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHO[i] / entotHO);
527 
528  meAllNHit_->Fill(double(nHit));
529  meBadDetHit_->Fill(double(nBad1));
530  meBadSubHit_->Fill(double(nBad2));
531  meBadIdHit_->Fill(double(nBad));
532  meHBNHit_->Fill(double(nHB));
533  meHENHit_->Fill(double(nHE));
534  meHONHit_->Fill(double(nHO));
535  meHFNHit_->Fill(double(nHF));
536 
537  for (int i = 0; i < ieta_bins_HB; i++) {
538  for (int j = 0; j < iphi_bins; j++) {
539  if (HBEneMap[i][j] != 0) {
540  meHBEneSum_->Fill(HBEneMap[i][j]);
541  meHBEneSum_vs_ieta_->Fill((i - eta_offset_HB), HBEneMap[i][j]);
542  meHBEneMap_->Fill((i - eta_offset_HB), j + 1, HBEneMap[i][j]);
543  }
544  }
545  }
546 
547  for (int i = 0; i < ieta_bins_HE; i++) {
548  for (int j = 0; j < iphi_bins; j++) {
549  if (HEEneMap[i][j] != 0) {
550  meHEEneSum_->Fill(HEEneMap[i][j]);
551  meHEEneSum_vs_ieta_->Fill((i - eta_offset_HE), HEEneMap[i][j]);
552  meHEEneMap_->Fill((i - eta_offset_HE), j + 1, HEEneMap[i][j]);
553  }
554  }
555  }
556 
557  for (int i = 0; i < ieta_bins_HO; i++) {
558  for (int j = 0; j < iphi_bins; j++) {
559  if (HOEneMap[i][j] != 0) {
560  meHOEneSum_->Fill(HOEneMap[i][j]);
561  meHOEneSum_vs_ieta_->Fill((i - eta_offset_HO), HOEneMap[i][j]);
562  meHOEneMap_->Fill((i - eta_offset_HO), j + 1, HOEneMap[i][j]);
563  }
564  }
565  }
566 
567  for (int i = 0; i < ieta_bins_HF; i++) {
568  for (int j = 0; j < iphi_bins; j++) {
569  if (HFEneMap[i][j] != 0) {
570  meHFEneSum_->Fill(HFEneMap[i][j]);
571  meHFEneSum_vs_ieta_->Fill((i - eta_offset_HF), HFEneMap[i][j]);
572  meHFEneMap_->Fill((i - eta_offset_HF), j + 1, HFEneMap[i][j]);
573  }
574  }
575  }
576 
577  if (verbose_ > 0)
578  edm::LogVerbatim("HcalSim") << "HcalSimHitCheck::analyzeHits: HB " << nHB << " HE " << nHE << " HO " << nHO
579  << " HF " << nHF << " Bad " << nBad << " All " << nHit;
580 }
581 
582 #include "FWCore/Framework/interface/MakerMacros.h"
583 DEFINE_FWK_MODULE(HcalSimHitCheck);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:38
HcalSimHitCheck::meHFEneSum_vs_ieta_
TProfile * meHFEneSum_vs_ieta_
Definition: HcalSimHitCheck.cc:83
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:40
HcalSimHitCheck::hcons_
const HcalDDDRecConstants * hcons_
Definition: HcalSimHitCheck.cc:48
mps_fire.i
i
Definition: mps_fire.py:428
HcalSimHitCheck::outFile_
std::string outFile_
Definition: HcalSimHitCheck.cc:65
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
HcalSimHitCheck::tok_HRNDC_
edm::ESGetToken< HcalDDDRecConstants, HcalRecNumberingRecord > tok_HRNDC_
Definition: HcalSimHitCheck.cc:52
HcalSimHitCheck::meHBEneSum_vs_ieta_
TProfile * meHBEneSum_vs_ieta_
Definition: HcalSimHitCheck.cc:83
HcalSimHitCheck::g4Label
std::string g4Label
Definition: HcalSimHitCheck.cc:65
HcalSimHitCheck::meHOEneSum_vs_ieta_
TProfile * meHOEneSum_vs_ieta_
Definition: HcalSimHitCheck.cc:83
HcalSimHitCheck::meEtaPhiHitDepth_
std::vector< TH2D * > meEtaPhiHitDepth_
Definition: HcalSimHitCheck.cc:75
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
HcalSimHitCheck::analyzeHits
void analyzeHits(std::vector< PCaloHit > &)
Definition: HcalSimHitCheck.cc:329
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:128
HcalDetId::subdet
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
HcalDetId::iphi
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157
HcalSimHitCheck::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: HcalSimHitCheck.cc:305
HcalDetId::ieta
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
HcalSimHitCheck::meHOL10EneP_
TProfile * meHOL10EneP_
Definition: HcalSimHitCheck.cc:87
edm::ParameterSetDescription::add
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:95
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
HcalSimHitCheck::beginRun
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: HcalSimHitCheck.cc:118
HcalSimHitCheck::meHEEneSum_vs_ieta_
TProfile * meHEEneSum_vs_ieta_
Definition: HcalSimHitCheck.cc:83
HLT_FULL_cff.depth
tuple depth
Definition: HLT_FULL_cff.py:15256
HcalSimHitCheck::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HcalSimHitCheck.cc:107
HcalDDDRecConstants::getMaxDepth
int getMaxDepth(const int &type) const
Definition: HcalDDDRecConstants.h:88
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:86956
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
HcalSimHitCheck::endRun
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: HcalSimHitCheck.cc:43
HcalSimHitCheck::tok_hits_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hits_
Definition: HcalSimHitCheck.cc:69
HcalSimHitCheck::HcalSimHitCheck
HcalSimHitCheck(const edm::ParameterSet &ps)
Definition: HcalSimHitCheck.cc:91
HcalSimHitCheck::meHFL10EneP_
TProfile * meHFL10EneP_
Definition: HcalSimHitCheck.cc:87
HcalSimHitCheck::meHEL10EneP_
TProfile * meHEL10EneP_
Definition: HcalSimHitCheck.cc:87
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
HcalSimHitCheck::hcalHits
std::string hcalHits
Definition: HcalSimHitCheck.cc:65
HcalSimHitCheck::meHBL10EneP_
TProfile * meHBL10EneP_
Definition: HcalSimHitCheck.cc:87
HcalHitRelabeller::relabel
DetId relabel(const uint32_t testId) const
Definition: HcalHitRelabeller.cc:49
HcalDetId::depth
constexpr int depth() const
get the tower depth
Definition: HcalDetId.h:164
HcalDDDRecConstants::getEtaRange
std::pair< int, int > getEtaRange(const int &i) const
Definition: HcalDDDRecConstants.h:74
edm::Run
Definition: Run.h:45
HcalDDDRecConstants::getNPhi
int getNPhi(const int &type) const
Definition: HcalDDDRecConstants.h:93
DetId::det
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
Generated for CMSSW Reference Manual by doxygen 1.8.5