CMS 3D CMS Logo

HOSimHitStudy.cc
Go to the documentation of this file.
3 
7 
13 
16 
22 
25 
26 #include <CLHEP/Units/GlobalPhysicalConstants.h>
27 #include <CLHEP/Units/GlobalSystemOfUnits.h>
28 
29 #include <TH1F.h>
30 #include <TProfile.h>
31 #include <TProfile2D.h>
32 
33 #include <fstream>
34 #include <iostream>
35 #include <iomanip>
36 #include <memory>
37 #include <map>
38 #include <string>
39 #include <vector>
40 
41 class HOSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
42 public:
44  ~HOSimHitStudy() override = default;
45  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
46 
47 protected:
48  void beginJob() override;
49  void beginRun(edm::Run const &, edm::EventSetup const &) override {}
50  void endJob() override {}
51  void endRun(edm::Run const &, edm::EventSetup const &) override {}
52  void analyze(const edm::Event &e, const edm::EventSetup &c) override;
53 
54  void analyzeHits();
55 
56 private:
58  const std::vector<std::string> hitLab_;
60  const double tcut_;
61  const bool scheme_, print_;
63  const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>> toks_calo_;
64  std::vector<PCaloHit> ecalHits, hcalHits;
65  TH1F *hit_[3], *time_[3], *edepTW_[3], *edepTWT_[3];
66  TH1F *edep_[3], *hitTow_[3], *eneInc_, *etaInc_, *phiInc_;
67  TH1F *edepT_[3], *eEB_, *eEBHB_, *eEBHBHO_, *eEBHBHOT_;
68  TH1F *edepZon_[3], *edepZonT_[3], *eEBT_, *eEBHBT_;
69  TH1F *eHOE17_[15], *eHOE18_[15], *eHOE_[15];
70  TH1F *eHOE17T_[15], *eHOE18T_[15], *eHOET_[15];
71  TH1F *eHOEta17_[15], *eHOEta18_[15], *eHOEta_[15];
72  TH1F *eHOEta17T_[15], *eHOEta18T_[15], *eHOEtaT_[15];
73  TH1F *nHOE1_[15], *nHOE1T_[15], *nHOEta1_[15], *nHOEta1T_[15];
74  TProfile *eHO1_, *eHO1T_, *eHO17_, *eHO17T_, *eHO18_, *eHO18T_;
75  TProfile *nHO1_, *nHO1T_, *nHOE2_[15], *nHOE2T_[15];
76  TProfile *nHOEta2_[15], *nHOEta2T_[15];
77  TProfile2D *eHO2_, *eHO2T_, *nHO2_, *nHO2T_;
78  double eInc, etaInc, phiInc;
79 };
80 
82  : g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
83  hitLab_(ps.getParameter<std::vector<std::string>>("HitCollection")),
84  maxEnergy_(ps.getUntrackedParameter<double>("MaxEnergy", 200.0)),
85  scaleEB_(ps.getUntrackedParameter<double>("ScaleEB", 1.0)),
86  scaleHB_(ps.getUntrackedParameter<double>("ScaleHB", 100.0)),
87  scaleHO_(ps.getUntrackedParameter<double>("ScaleHO", 2.0)),
88  tcut_(ps.getUntrackedParameter<double>("TimeCut", 100.0)),
89  scheme_(ps.getUntrackedParameter<bool>("TestNumbering", false)),
90  print_(ps.getUntrackedParameter<bool>("PrintExcessEnergy", true)),
91  tok_evt_(consumes<edm::HepMCProduct>(
92  edm::InputTag(ps.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")))),
93  toks_calo_{edm::vector_transform(hitLab_, [this](const std::string &name) {
94  return consumes<edm::PCaloHitContainer>(edm::InputTag{g4Label_, name});
95  })} {
96  usesResource(TFileService::kSharedResource);
97 
98  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Module Label: " << g4Label_ << " Hits: " << hitLab_[0] << ", "
99  << hitLab_[1] << " MaxEnergy: " << maxEnergy_ << " Scale factor for EB " << scaleEB_
100  << ", for HB " << scaleHB_ << " and for HO " << scaleHO_ << " time Cut " << tcut_;
101 }
102 
105  std::vector<std::string> labels = {"EcalHitsEB", "HcalHits"};
106  desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
107  desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
108  desc.addUntracked<std::vector<std::string>>("HitCollection", labels);
109  desc.addUntracked<double>("MaxEnergy", 50.0);
110  desc.addUntracked<double>("ScaleEB", 1.02);
111  desc.addUntracked<double>("ScaleHB", 104.4);
112  desc.addUntracked<double>("ScaleHO", 2.33);
113  desc.addUntracked<double>("TimeCut", 100.0);
114  desc.addUntracked<bool>("PrintExcessEnergy", true);
115  desc.addUntracked<bool>("TestNumbering", false);
116  descriptions.add("HOSimHitStudy", desc);
117 }
118 
121 
122  if (!tfile.isAvailable())
123  throw cms::Exception("BadConfig") << "TFileService unavailable: "
124  << "please add it to config file";
125  std::string dets[3] = {"EB", "HB", "HO"};
126  char name[60], title[100];
127  double ymax = maxEnergy_;
128  sprintf(title, "Incident Energy (GeV)");
129  eneInc_ = tfile->make<TH1F>("EneInc", title, 1000, 0., ymax);
130  eneInc_->GetXaxis()->SetTitle(title);
131  eneInc_->GetYaxis()->SetTitle("Events");
132  sprintf(title, "Incident #eta");
133  etaInc_ = tfile->make<TH1F>("EtaInc", title, 200, -5., 5.);
134  etaInc_->GetXaxis()->SetTitle(title);
135  etaInc_->GetYaxis()->SetTitle("Events");
136  sprintf(title, "Incident #phi");
137  phiInc_ = tfile->make<TH1F>("PhiInc", title, 200, -3.1415926, 3.1415926);
138  phiInc_->GetXaxis()->SetTitle(title);
139  phiInc_->GetYaxis()->SetTitle("Events");
140  int itcut = (int)(tcut_);
141  for (int i = 0; i < 3; i++) {
142  sprintf(name, "Hit%d", i);
143  sprintf(title, "Number of hits in %s", dets[i].c_str());
144  hit_[i] = tfile->make<TH1F>(name, title, 1000, 0., 20000.);
145  hit_[i]->GetXaxis()->SetTitle(title);
146  hit_[i]->GetYaxis()->SetTitle("Events");
147  sprintf(name, "Time%d", i);
148  sprintf(title, "Time of the hit (ns) in %s", dets[i].c_str());
149  time_[i] = tfile->make<TH1F>(name, title, 1200, 0., 1200.);
150  time_[i]->GetXaxis()->SetTitle(title);
151  time_[i]->GetYaxis()->SetTitle("Hits");
152  if (i > 0)
153  ymax = 1.0;
154  else
155  ymax = 50.0;
156  sprintf(name, "Edep%d", i);
157  sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
158  edep_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
159  edep_[i]->GetXaxis()->SetTitle(title);
160  edep_[i]->GetYaxis()->SetTitle("Hits");
161  sprintf(name, "EdepT%d", i);
162  sprintf(title, "Energy deposit (GeV) in %s for t < %d ns", dets[i].c_str(), itcut);
163  edepT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
164  edepT_[i]->GetXaxis()->SetTitle(title);
165  edepT_[i]->GetYaxis()->SetTitle("Hits");
166  sprintf(name, "HitTow%d", i);
167  sprintf(title, "Number of towers with hits in %s", dets[i].c_str());
168  hitTow_[i] = tfile->make<TH1F>(name, title, 1000, 0., 20000.);
169  hitTow_[i]->GetXaxis()->SetTitle(title);
170  hitTow_[i]->GetYaxis()->SetTitle("Events");
171  if (i > 0)
172  ymax = 0.05 * maxEnergy_;
173  else
174  ymax = maxEnergy_;
175  sprintf(name, "EdepTW%d", i);
176  sprintf(title, "Energy deposit (GeV) in %s Tower", dets[i].c_str());
177  edepTW_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
178  edepTW_[i]->GetXaxis()->SetTitle(title);
179  edepTW_[i]->GetYaxis()->SetTitle("Towers");
180  sprintf(name, "EdepTWT%d", i);
181  sprintf(title, "Energy deposit (GeV) in %s Tower for t < %d ns", dets[i].c_str(), itcut);
182  edepTWT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
183  edepTWT_[i]->GetXaxis()->SetTitle(title);
184  edepTWT_[i]->GetYaxis()->SetTitle("Towers");
185  sprintf(name, "EdepZone%d", i);
186  sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
187  edepZon_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
188  edepZon_[i]->GetXaxis()->SetTitle(title);
189  edepZon_[i]->GetYaxis()->SetTitle("Events");
190  sprintf(name, "EdepZoneT%d", i);
191  sprintf(title, "Energy deposit (GeV) in %s for t < %d ns", dets[i].c_str(), itcut);
192  edepZonT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
193  edepZonT_[i]->GetXaxis()->SetTitle(title);
194  edepZonT_[i]->GetYaxis()->SetTitle("Events");
195  }
196  sprintf(title, "Energy Measured in EB (GeV)");
197  eEB_ = tfile->make<TH1F>("EEB", title, 5000, 0., maxEnergy_);
198  eEB_->GetXaxis()->SetTitle(title);
199  eEB_->GetYaxis()->SetTitle("Events");
200  sprintf(title, "Energy Measured in EB+HB (GeV)");
201  eEBHB_ = tfile->make<TH1F>("EEBHB", title, 5000, 0., maxEnergy_);
202  eEBHB_->GetXaxis()->SetTitle(title);
203  eEBHB_->GetYaxis()->SetTitle("Events");
204  sprintf(title, "Energy Measured in EB+HB+HO (GeV)");
205  eEBHBHO_ = tfile->make<TH1F>("EEBHBHO", title, 5000, 0., maxEnergy_);
206  eEBHBHO_->GetXaxis()->SetTitle(title);
207  eEBHBHO_->GetYaxis()->SetTitle("Events");
208  sprintf(title, "Energy Measured in EB (GeV) for t < %d ns", itcut);
209  eEBT_ = tfile->make<TH1F>("EEBT", title, 5000, 0., maxEnergy_);
210  eEBT_->GetXaxis()->SetTitle(title);
211  eEBT_->GetYaxis()->SetTitle("Events");
212  sprintf(title, "Energy Measured in EB+HB (GeV) for t < %d ns", itcut);
213  eEBHBT_ = tfile->make<TH1F>("EEBHBT", title, 5000, 0., maxEnergy_);
214  eEBHBT_->GetXaxis()->SetTitle(title);
215  eEBHBT_->GetYaxis()->SetTitle("Events");
216  sprintf(title, "Energy Measured in EB+HB+HO (GeV) for t < %d ns", itcut);
217  eEBHBHOT_ = tfile->make<TH1F>("EEBHBHOT", title, 5000, 0., maxEnergy_);
218  eEBHBHOT_->GetXaxis()->SetTitle(title);
219  eEBHBHOT_->GetYaxis()->SetTitle("Events");
220  sprintf(title, "SimHit energy in HO");
221  eHO1_ = tfile->make<TProfile>("EHO1", title, 30, -1.305, 1.305);
222  eHO1_->GetXaxis()->SetTitle(title);
223  eHO1_->GetYaxis()->SetTitle("Events");
224  eHO2_ = tfile->make<TProfile2D>("EHO2", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
225  eHO2_->GetXaxis()->SetTitle(title);
226  eHO2_->GetYaxis()->SetTitle("Events");
227  sprintf(title, "SimHit energy in HO Layer 17");
228  eHO17_ = tfile->make<TProfile>("EHO17", title, 30, -1.305, 1.305);
229  eHO17_->GetXaxis()->SetTitle(title);
230  eHO17_->GetYaxis()->SetTitle("Events");
231  sprintf(title, "SimHit energy in HO Layer 18");
232  eHO18_ = tfile->make<TProfile>("EHO18", title, 30, -1.305, 1.305);
233  eHO18_->GetXaxis()->SetTitle(title);
234  eHO18_->GetYaxis()->SetTitle("Events");
235  sprintf(title, "SimHit energy in HO for t < %d ns", itcut);
236  eHO1T_ = tfile->make<TProfile>("EHO1T", title, 30, -1.305, 1.305);
237  eHO1T_->GetXaxis()->SetTitle(title);
238  eHO1T_->GetYaxis()->SetTitle("Events");
239  eHO2T_ = tfile->make<TProfile2D>("EHO2T", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
240  eHO2T_->GetXaxis()->SetTitle(title);
241  eHO2T_->GetYaxis()->SetTitle("Events");
242  sprintf(title, "SimHit energy in HO Layer 17 for t < %d ns", itcut);
243  eHO17T_ = tfile->make<TProfile>("EHO17T", title, 30, -1.305, 1.305);
244  eHO17T_->GetXaxis()->SetTitle(title);
245  eHO17T_->GetYaxis()->SetTitle("Events");
246  sprintf(title, "SimHit energy in HO Layer 18 for t < %d ns", itcut);
247  eHO18T_ = tfile->make<TProfile>("EHO18T", title, 30, -1.305, 1.305);
248  eHO18T_->GetXaxis()->SetTitle(title);
249  eHO18T_->GetYaxis()->SetTitle("Events");
250  sprintf(title, "Number of layers hit in HO");
251  nHO1_ = tfile->make<TProfile>("NHO1", title, 30, -1.305, 1.305);
252  nHO1_->GetXaxis()->SetTitle(title);
253  nHO1_->GetYaxis()->SetTitle("Events");
254  nHO2_ = tfile->make<TProfile2D>("NHO2", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
255  nHO2_->GetXaxis()->SetTitle(title);
256  nHO2_->GetYaxis()->SetTitle("Events");
257  sprintf(title, "Number of layers hit in HO for t < %d ns", itcut);
258  nHO1T_ = tfile->make<TProfile>("NHO1T", title, 30, -1.305, 1.305);
259  nHO1T_->GetXaxis()->SetTitle(title);
260  nHO1T_->GetYaxis()->SetTitle("Events");
261  nHO2T_ = tfile->make<TProfile2D>("NHO2T", title, 30, -1.305, 1.305, 72, -3.1415926, 3.1415926);
262  nHO2T_->GetXaxis()->SetTitle(title);
263  nHO2T_->GetYaxis()->SetTitle("Events");
264  for (int i = 0; i < 15; i++) {
265  sprintf(name, "EHOE%d", i + 1);
266  sprintf(title, "SimHit energy in HO (Beam in #eta=%d bin)", i + 1);
267  eHOE_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
268  eHOE_[i]->GetXaxis()->SetTitle(title);
269  eHOE_[i]->GetYaxis()->SetTitle("Events");
270  sprintf(name, "EHOE17%d", i + 1);
271  sprintf(title, "SimHit energy in Layer 17 (Beam in #eta=%d bin)", i + 1);
272  eHOE17_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
273  eHOE17_[i]->GetXaxis()->SetTitle(title);
274  eHOE17_[i]->GetYaxis()->SetTitle("Events");
275  sprintf(name, "EHOE18%d", i + 1);
276  sprintf(title, "SimHit energy in Layer 18 (Beam in #eta=%d bin)", i + 1);
277  eHOE18_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
278  eHOE18_[i]->GetXaxis()->SetTitle(title);
279  eHOE18_[i]->GetYaxis()->SetTitle("Events");
280  sprintf(name, "EHOE%dT", i + 1);
281  sprintf(title, "SimHit energy in HO (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
282  eHOET_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
283  eHOET_[i]->GetXaxis()->SetTitle(title);
284  eHOET_[i]->GetYaxis()->SetTitle("Events");
285  sprintf(name, "EHOE17%dT", i + 1);
286  sprintf(title, "SimHit energy in Layer 17 (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
287  eHOE17T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
288  eHOE17T_[i]->GetXaxis()->SetTitle(title);
289  eHOE17T_[i]->GetYaxis()->SetTitle("Events");
290  sprintf(name, "EHOE18%dT", i + 1);
291  sprintf(title, "SimHit energy in Layer 18 (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
292  eHOE18T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
293  eHOE18T_[i]->GetXaxis()->SetTitle(title);
294  eHOE18T_[i]->GetYaxis()->SetTitle("Events");
295  sprintf(name, "EHOEta%d", i + 1);
296  sprintf(title, "SimHit energy in HO #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
297  eHOEta_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
298  eHOEta_[i]->GetXaxis()->SetTitle(title);
299  eHOEta_[i]->GetYaxis()->SetTitle("Events");
300  sprintf(name, "EHOEta17%d", i + 1);
301  sprintf(title, "SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
302  eHOEta17_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
303  eHOEta17_[i]->GetXaxis()->SetTitle(title);
304  eHOEta17_[i]->GetYaxis()->SetTitle("Events");
305  sprintf(name, "EHOEta18%d", i + 1);
306  sprintf(title, "SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
307  eHOEta18_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
308  eHOEta18_[i]->GetXaxis()->SetTitle(title);
309  eHOEta18_[i]->GetYaxis()->SetTitle("Events");
310  sprintf(name, "EHOEta%dT", i + 1);
311  sprintf(title, "SimHit energy in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
312  eHOEtaT_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
313  eHOEtaT_[i]->GetXaxis()->SetTitle(title);
314  eHOEtaT_[i]->GetYaxis()->SetTitle("Events");
315  sprintf(name, "EHOEta17%dT", i + 1);
316  sprintf(title, "SimHit energy in Layer 17 #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
317  eHOEta17T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
318  eHOEta17T_[i]->GetXaxis()->SetTitle(title);
319  eHOEta17T_[i]->GetYaxis()->SetTitle("Events");
320  sprintf(name, "EHOEta18%dT", i + 1);
321  sprintf(title, "SimHit energy in Layer 18 #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
322  eHOEta18T_[i] = tfile->make<TH1F>(name, title, 1000, 0., 0.25);
323  eHOEta18T_[i]->GetXaxis()->SetTitle(title);
324  eHOEta18T_[i]->GetYaxis()->SetTitle("Events");
325  sprintf(name, "NHOE1%d", i + 1);
326  sprintf(title, "Number of layers hit in HO (Beam in #eta=%d bin)", i + 1);
327  nHOE1_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
328  nHOE1_[i]->GetXaxis()->SetTitle(title);
329  nHOE1_[i]->GetYaxis()->SetTitle("Events");
330  sprintf(name, "NHOE2%d", i + 1);
331  nHOE2_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
332  nHOE2_[i]->GetXaxis()->SetTitle(title);
333  nHOE2_[i]->GetYaxis()->SetTitle("Events");
334  sprintf(name, "NHOE1%dT", i + 1);
335  sprintf(title, "Number of layers hit in HO (Beam in #eta=%d bin, t < %d ns)", i + 1, itcut);
336  nHOE1T_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
337  nHOE1T_[i]->GetXaxis()->SetTitle(title);
338  nHOE1T_[i]->GetYaxis()->SetTitle("Events");
339  sprintf(name, "NHOE2%dT", i + 1);
340  nHOE2T_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
341  nHOE2T_[i]->GetXaxis()->SetTitle(title);
342  nHOE2T_[i]->GetYaxis()->SetTitle("Events");
343  sprintf(name, "NHOEta1%d", i + 1);
344  sprintf(title, "Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin)", i + 1, i + 1);
345  nHOEta1_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
346  nHOEta1_[i]->GetXaxis()->SetTitle(title);
347  nHOEta1_[i]->GetYaxis()->SetTitle("Events");
348  sprintf(name, "NHOEta2%d", i + 1);
349  nHOEta2_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
350  nHOEta2_[i]->GetXaxis()->SetTitle(title);
351  nHOEta2_[i]->GetYaxis()->SetTitle("Events");
352  sprintf(name, "NHOEta1%dT", i + 1);
353  sprintf(title, "Number of layers hit in HO #eta bin %d (Beam in #eta=%d bin, t < %d ns)", i + 1, i + 1, itcut);
354  nHOEta1T_[i] = tfile->make<TH1F>(name, title, 20, 0, 20.);
355  nHOEta1T_[i]->GetXaxis()->SetTitle(title);
356  nHOEta1T_[i]->GetYaxis()->SetTitle("Events");
357  sprintf(name, "NHOEta2%dT", i + 1);
358  nHOEta2T_[i] = tfile->make<TProfile>(name, title, 72, -3.1415926, 3.1415926);
359  nHOEta2T_[i]->GetXaxis()->SetTitle(title);
360  nHOEta2T_[i]->GetYaxis()->SetTitle("Events");
361  }
362 }
363 
365  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Run = " << e.id().run() << " Event = " << e.id().event();
366 
367  const edm::Handle<edm::HepMCProduct> &EvtHandle = e.getHandle(tok_evt_);
368  const HepMC::GenEvent *myGenEvent = EvtHandle->GetEvent();
369 
370  eInc = etaInc = phiInc = 0;
371  HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
372  if (p != myGenEvent->particles_end()) {
373  eInc = (*p)->momentum().e();
374  etaInc = (*p)->momentum().eta();
375  phiInc = (*p)->momentum().phi();
376  }
377 
378  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Energy = " << eInc << " Eta = " << etaInc
379  << " Phi = " << phiInc / CLHEP::deg;
380 
381  for (int i = 0; i < 2; i++) {
382  bool getHits = false;
383  if (i == 0)
384  ecalHits.clear();
385  else
386  hcalHits.clear();
387  const edm::Handle<edm::PCaloHitContainer> &hitsCalo = e.getHandle(toks_calo_[i]);
388  if (hitsCalo.isValid())
389  getHits = true;
390  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Input flag " << hitLab_[i] << " getHits flag " << getHits;
391 
392  if (getHits) {
393  unsigned int isiz;
394  if (i == 0) {
395  ecalHits.insert(ecalHits.end(), hitsCalo->begin(), hitsCalo->end());
396  isiz = ecalHits.size();
397  } else {
398  hcalHits.insert(hcalHits.end(), hitsCalo->begin(), hitsCalo->end());
399  isiz = hcalHits.size();
400  }
401  edm::LogVerbatim("HitStudy") << "HOSimHitStudy:: Hit buffer for " << hitLab_[i] << " has " << isiz << " hits";
402  }
403  }
404  analyzeHits();
405 }
406 
408  //initialize
409  int nhit[3];
410  double etot[3], etotT[3];
411  std::vector<unsigned int> ebID, hbID, hoID;
412  std::vector<double> ebEtow, hbEtow, hoEtow;
413  std::vector<double> ebEtowT, hbEtowT, hoEtowT;
414  for (int k = 0; k < 3; k++) {
415  nhit[k] = 0;
416  etot[k] = 0;
417  etotT[k] = 0;
418  }
419  eneInc_->Fill(eInc);
420  etaInc_->Fill(etaInc);
421  phiInc_->Fill(phiInc);
422 
423  double eHO17 = 0, eHO18 = 0, eHO17T = 0, eHO18T = 0;
424  double eHOE17[15], eHOE18[15], eHOE17T[15], eHOE18T[15];
425  for (int k = 0; k < 15; k++) {
426  eHOE17[k] = eHOE18[k] = eHOE17T[k] = eHOE18T[k] = 0;
427  }
428  // Loop over containers
429  for (int k = 0; k < 2; k++) {
430  int nHit;
431  if (k == 0) {
432  nHit = ecalHits.size();
433  } else {
434  nHit = hcalHits.size();
435  }
436  for (int i = 0; i < nHit; i++) {
437  double edep, time;
438  int indx = -1;
439  unsigned int id_;
440  if (k == 0) {
441  indx = 0;
442  edep = ecalHits[i].energy();
443  time = ecalHits[i].time();
444  id_ = ecalHits[i].id();
445  } else {
446  edep = hcalHits[i].energy();
447  time = hcalHits[i].time();
448  id_ = hcalHits[i].id();
449  int subdet, zside, depth, eta, phi, lay;
450  if (scheme_) {
452  } else {
453  HcalDetId id = HcalDetId(id_);
454  subdet = id.subdet();
455  zside = id.zside();
456  depth = id.depth();
457  eta = id.ietaAbs();
458  phi = id.iphi();
459  lay = -1;
460  }
461  edm::LogVerbatim("HitStudy") << "HOSimHitStudy:: Hit " << k << " Subdet:" << subdet << " zside:" << zside
462  << " depth:" << depth << " layer:" << lay << " eta:" << eta << " phi:" << phi;
463  if (subdet == static_cast<int>(HcalBarrel))
464  indx = 1;
465  else if (subdet == static_cast<int>(HcalOuter)) {
466  indx = 2;
467  if (lay == 18) {
468  eHO17 += edep;
469  if (time < tcut_)
470  eHO17T += edep;
471  if (eta >= 0 && eta < 15) {
472  eHOE17[eta] += edep;
473  if (time < tcut_)
474  eHOE17T[eta] += edep;
475  }
476  } else {
477  eHO18 += edep;
478  if (time < tcut_)
479  eHO18T += edep;
480  if (eta >= 0 && eta < 15) {
481  eHOE18[eta] += edep;
482  if (time < tcut_)
483  eHOE18T[eta] += edep;
484  }
485  }
486  }
487  }
488  if (indx >= 0) {
489  double edepT = edep;
490  time_[indx]->Fill(time);
491  edep_[indx]->Fill(edep);
492  etot[indx] += edep;
493  if (time < tcut_) {
494  etotT[indx] += edep;
495  edepT_[indx]->Fill(edep);
496  edepT = 0;
497  }
498  nhit[indx]++;
499  if (indx == 0) {
500  bool ok = false;
501  for (unsigned int j = 0; j < ebID.size(); j++) {
502  if (id_ == ebID[j]) {
503  ebEtow[j] += edep;
504  ebEtowT[j] += edepT;
505  ok = true;
506  break;
507  }
508  }
509  if (!ok) {
510  ebID.push_back(id_);
511  ebEtow.push_back(edep);
512  ebEtowT.push_back(edepT);
513  }
514  } else if (indx == 1) {
515  bool ok = false;
516  for (unsigned int j = 0; j < hbID.size(); j++) {
517  if (id_ == hbID[j]) {
518  hbEtow[j] += edep;
519  hbEtowT[j] += edepT;
520  ok = true;
521  break;
522  }
523  }
524  if (!ok) {
525  hbID.push_back(id_);
526  hbEtow.push_back(edep);
527  hbEtowT.push_back(edepT);
528  }
529  } else {
530  bool ok = false;
531  for (unsigned int j = 0; j < hoID.size(); j++) {
532  if (id_ == hoID[j]) {
533  hoEtow[j] += edep;
534  hoEtowT[j] += edepT;
535  ok = true;
536  break;
537  }
538  }
539  if (!ok) {
540  hoID.push_back(id_);
541  hoEtow.push_back(edep);
542  hoEtowT.push_back(edepT);
543  }
544  }
545  }
546  }
547  }
548 
549  // Now for towers and total energy deposits
550  for (int k = 0; k < 3; k++) {
551  hit_[k]->Fill(double(nhit[k]));
552  edepZon_[k]->Fill(etot[k]);
553  edepZonT_[k]->Fill(etotT[k]);
554  }
555  hitTow_[0]->Fill(double(ebEtow.size()));
556  for (unsigned int i = 0; i < ebEtow.size(); i++) {
557  edepTW_[0]->Fill(ebEtow[i]);
558  edepTWT_[0]->Fill(ebEtowT[i]);
559  }
560  hitTow_[1]->Fill(double(hbEtow.size()));
561  for (unsigned int i = 0; i < hbEtow.size(); i++) {
562  edepTW_[1]->Fill(hbEtow[i]);
563  edepTWT_[1]->Fill(hbEtowT[i]);
564  }
565  hitTow_[2]->Fill(double(hoEtow.size()));
566  for (unsigned int i = 0; i < hoEtow.size(); i++) {
567  edepTW_[2]->Fill(hoEtow[i]);
568  edepTWT_[2]->Fill(hoEtowT[i]);
569  }
570  double eEB = scaleEB_ * etot[0];
571  double eEBHB = eEB + scaleHB_ * etot[1];
572  double eEBHBHO = eEBHB + scaleHB_ * scaleHO_ * etot[2];
573  eEB_->Fill(eEB);
574  eEBHB_->Fill(eEBHB);
575  eEBHBHO_->Fill(eEBHBHO);
576  double eEBT = scaleEB_ * etotT[0];
577  double eEBHBT = eEBT + scaleHB_ * etotT[1];
578  double eEBHBHOT = eEBHBT + scaleHB_ * scaleHO_ * etotT[2];
579  eEBT_->Fill(eEBT);
580  eEBHBT_->Fill(eEBHBT);
581  eEBHBHOT_->Fill(eEBHBHOT);
582  eHO1_->Fill(etaInc, eHO17 + eHO18);
583  eHO2_->Fill(etaInc, phiInc, eHO17 + eHO18);
584  eHO17_->Fill(etaInc, eHO17);
585  eHO18_->Fill(etaInc, eHO18);
586  eHO1T_->Fill(etaInc, eHO17T + eHO18T);
587  eHO2T_->Fill(etaInc, phiInc, eHO17T + eHO18T);
588  eHO17T_->Fill(etaInc, eHO17T);
589  eHO18T_->Fill(etaInc, eHO18T);
590  int nHO = 0, nHOT = 0;
591  if (eHO17 > 0)
592  nHO++;
593  if (eHO17T > 0)
594  nHOT++;
595  if (eHO18 > 0)
596  nHO++;
597  if (eHO18T > 0)
598  nHOT++;
599  nHO1_->Fill(etaInc, (double)(nHO));
600  nHO2_->Fill(etaInc, phiInc, (double)(nHO));
601  nHO1T_->Fill(etaInc, (double)(nHOT));
602  nHO2T_->Fill(etaInc, phiInc, (double)(nHOT));
603  int ieta = 15;
604  for (int k = 0; k < 15; ++k) {
605  if (std::abs(etaInc) < 0.087 * (k + 1)) {
606  ieta = k;
607  break;
608  }
609  }
610  if (ieta >= 0 && ieta < 15) {
611  eHOE_[ieta]->Fill(eHO17 + eHO18);
612  eHOE17_[ieta]->Fill(eHO17);
613  eHOE18_[ieta]->Fill(eHO18);
614  eHOET_[ieta]->Fill(eHO17T + eHO18T);
615  eHOE17T_[ieta]->Fill(eHO17T);
616  eHOE18T_[ieta]->Fill(eHO18T);
617  eHOEta_[ieta]->Fill(eHOE17[ieta] + eHOE18[ieta]);
618  eHOEta17_[ieta]->Fill(eHOE17[ieta]);
619  eHOEta18_[ieta]->Fill(eHOE18[ieta]);
620  nHOE1_[ieta]->Fill((double)(nHO));
621  nHOE2_[ieta]->Fill(phiInc, (double)(nHO));
622  nHOE1T_[ieta]->Fill((double)(nHOT));
623  nHOE2T_[ieta]->Fill(phiInc, (double)(nHOT));
624  int nHOE = 0, nHOET = 0;
625  if (eHOE17[ieta] > 0)
626  nHOE++;
627  if (eHOE17T[ieta] > 0)
628  nHOET++;
629  if (eHOE18[ieta] > 0)
630  nHOE++;
631  if (eHOE18T[ieta] > 0)
632  nHOET++;
633  nHOEta1_[ieta]->Fill((double)(nHOE));
634  nHOEta2_[ieta]->Fill(phiInc, (double)(nHOE));
635  nHOEta1T_[ieta]->Fill((double)(nHOET));
636  nHOEta2T_[ieta]->Fill(phiInc, (double)(nHOET));
637  }
638 
639  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::analyzeHits: Hits in EB " << nhit[0] << " in " << ebEtow.size()
640  << " towers with total E " << etot[0] << "|" << etotT[0]
641  << "\n Hits in HB " << nhit[1] << " in " << hbEtow.size()
642  << " towers with total E " << etot[1] << "|" << etotT[1]
643  << "\n Hits in HO " << nhit[2] << " in " << hoEtow.size()
644  << " towers with total E " << etot[2] << "|" << etotT[2]
645  << "\n Energy in EB " << eEB << "|" << eEBT << " with HB "
646  << eEBHB << "|" << eEBHBT << " and with HO " << eEBHBHO << "|" << eEBHBHOT
647  << "\n E in HO layers " << eHO17 << "|" << eHO17T << " "
648  << eHO18 << "|" << eHO18T << " number of HO hits " << nHO << "|" << nHOT;
649 
650  if (eEBHBHO > 0.75 * maxEnergy_ && print_) {
651  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::Event with excessive energy: EB = " << eEB << " EB+HB = " << eEBHB
652  << " EB+HB+HO = " << eEBHBHO;
653  const std::string Dets[3] = {"EB", "HB", "HO"};
654  for (int k = 0; k < 2; k++) {
655  int nHit;
656  if (k == 0) {
657  nHit = ecalHits.size();
658  } else {
659  nHit = hcalHits.size();
660  }
661  for (int i = 0; i < nHit; i++) {
662  double edep, time;
663  int indx = -1;
664  unsigned int id_;
665  int ieta, iphi, depth = 0;
666  if (k == 0) {
667  indx = 0;
668  edep = ecalHits[i].energy();
669  time = ecalHits[i].time();
670  id_ = ecalHits[i].id();
671  EBDetId id = EBDetId(id_);
672  ieta = id.ieta();
673  iphi = id.iphi();
674  } else {
675  indx = -1;
676  edep = hcalHits[i].energy();
677  time = hcalHits[i].time();
678  id_ = hcalHits[i].id();
679  HcalDetId id = HcalDetId(id_);
680  int subdet = id.subdet();
681  if (subdet == static_cast<int>(HcalBarrel))
682  indx = 1;
683  else if (subdet == static_cast<int>(HcalOuter))
684  indx = 2;
685  ieta = id.ieta();
686  iphi = id.iphi();
687  depth = id.depth();
688  }
689  if (indx >= 0) {
690  edm::LogVerbatim("HitStudy") << "HOSimHitStudy::" << Dets[indx] << " " << i << std::hex << id_ << std::dec
691  << " (" << ieta << "|" << iphi << "|" << depth << ") " << std::setw(8) << edep
692  << " " << std::setw(8) << time;
693  }
694  }
695  }
696  }
697 }
698 
699 //define this as a plug-in
static const std::string kSharedResource
Definition: TFileService.h:76
Log< level::Info, true > LogVerbatim
void beginRun(edm::Run const &, edm::EventSetup const &) override
TProfile * nHOE2_[15]
TProfile * nHOE2T_[15]
void analyze(const edm::Event &e, const edm::EventSetup &c) override
std::ostream & print_(std::ostream &os, value_type const &hash)
Definition: Hash.cc:92
TProfile * eHO1_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const double scaleHO_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
TProfile * eHO18_
TProfile2D * eHO2_
TH1F * eHOE18_[15]
const double tcut_
TH1F * eHOE18T_[15]
TProfile2D * eHO2T_
int zside(DetId const &)
TH1F * eHOE17T_[15]
void endRun(edm::Run const &, edm::EventSetup const &) override
TH1F * time_[3]
TH1F * eHOEta18_[15]
auto vector_transform(std::vector< InputType > const &input, Function predicate) -> std::vector< typename std::remove_cv< typename std::remove_reference< decltype(predicate(input.front()))>::type >::type >
Definition: transform.h:11
const std::vector< std::string > hitLab_
TH1F * eHOEta17T_[15]
const double scaleHB_
TH1F * edepZonT_[3]
TH1F * nHOE1T_[15]
const bool scheme_
TH1F * edepTWT_[3]
const double maxEnergy_
TH1F * eHOE17_[15]
TH1F * eHOEta18T_[15]
TH1F * edep_[3]
Definition: tfile.py:1
void endJob() override
TProfile * eHO17_
TH1F * nHOE1_[15]
TH1F * edepTW_[3]
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const bool print_
TProfile * nHOEta2_[15]
TProfile * eHO18T_
TProfile * eHO17T_
TProfile * nHO1T_
TProfile2D * nHO2_
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
TH1F * eHOE_[15]
TH1F * eHOET_[15]
TH1F * nHOEta1T_[15]
TProfile * nHOEta2T_[15]
void beginJob() override
TH1F * hitTow_[3]
const std::vector< edm::EDGetTokenT< edm::PCaloHitContainer > > toks_calo_
TProfile2D * nHO2T_
~HOSimHitStudy() override=default
TH1F * edepZon_[3]
void add(std::string const &label, ParameterSetDescription const &psetDescription)
TH1F * eHOEta_[15]
TProfile * nHO1_
bool isValid() const
Definition: HandleBase.h:70
std::vector< PCaloHit > hcalHits
HLT enums.
std::vector< PCaloHit > ecalHits
TH1F * edepT_[3]
TProfile * eHO1T_
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
TH1F * hit_[3]
const std::string g4Label_
TH1F * eHOEta17_[15]
HOSimHitStudy(const edm::ParameterSet &ps)
const edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
const double scaleEB_
Definition: Run.h:45
TH1F * nHOEta1_[15]
TH1F * eHOEtaT_[15]