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