CMS 3D CMS Logo
CaloSimHitStudy.cc
Go to the documentation of this file.
1 #include "FWCore/Framework/interface/Frameworkfwd.h"
2 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
3 
4 #include "FWCore/Framework/interface/Event.h"
5 #include "FWCore/Framework/interface/EventSetup.h"
6 #include "FWCore/Framework/interface/ESHandle.h"
7 #include "FWCore/Framework/interface/MakerMacros.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 #include "CommonTools/UtilAlgos/interface/TFileService.h"
17 
18 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
19 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
20 
21 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
22 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
23 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
24 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
25 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
26 #include "SimG4CMS/Calo/interface/CaloHitID.h"
27 
28 #include <TH1F.h>
29 
30 #include <memory>
31 #include <iostream>
32 #include <fstream>
33 #include <vector>
34 #include <map>
35 #include <string>
36 
37 class CaloSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
38 public:
39  CaloSimHitStudy(const edm::ParameterSet& ps);
40  ~CaloSimHitStudy() override {}
41  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
42 
43 protected:
44  void beginJob() override {}
45  void analyze(edm::Event const&, edm::EventSetup const&) override;
46  void beginRun(edm::Run const&, edm::EventSetup const&) override {}
47  void endRun(edm::Run const&, edm::EventSetup const&) override {}
48 
49  void analyzeHits(std::vector<PCaloHit>&, int);
50  void analyzeHits(edm::Handle<edm::PSimHitContainer>&, int);
51  void analyzeHits(std::vector<PCaloHit>&, std::vector<PCaloHit>&, std::vector<PCaloHit>&);
52 
53 private:
54  std::string g4Label, hitLab[4];
55  edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
56  edm::EDGetTokenT<edm::PCaloHitContainer> toks_calo_[4];
57  edm::EDGetTokenT<edm::PSimHitContainer> toks_track_[3];
58  edm::EDGetTokenT<edm::PSimHitContainer> toks_tkHigh_[6];
59  edm::EDGetTokenT<edm::PSimHitContainer> toks_tkLow_[6];
60  std::string muonLab[3], tkHighLab[6], tkLowLab[6];
61  double tmax_, eMIP_;
62  bool storeRL_, testNumber_;
63 
64  TH1F *hit_[9], *time_[9], *edepEM_[9], *edepHad_[9], *edep_[9];
65  TH1F *etot_[9], *etotg_[9], *timeAll_[9], *hitMu, *hitHigh;
66  TH1F *hitLow, *eneInc_, *etaInc_, *phiInc_, *ptInc_;
67  TH1F *hitTk_[15], *edepTk_[15], *tofTk_[15], *edepC_[9], *edepT_[9];
68 };
69 
70 CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
71  usesResource(TFileService::kSharedResource);
72 
73  tok_evt_ =
74  consumes<edm::HepMCProduct>(edm::InputTag(ps.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")));
75  g4Label = ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits");
76  hitLab[0] = ps.getUntrackedParameter<std::string>("EBCollection", "EcalHitsEB");
77  hitLab[1] = ps.getUntrackedParameter<std::string>("EECollection", "EcalHitsEE");
78  hitLab[2] = ps.getUntrackedParameter<std::string>("ESCollection", "EcalHitsES");
79  hitLab[3] = ps.getUntrackedParameter<std::string>("HCCollection", "HcalHits");
80 
81  double maxEnergy_ = ps.getUntrackedParameter<double>("MaxEnergy", 200.0);
82  tmax_ = ps.getUntrackedParameter<double>("TimeCut", 100.0);
83  eMIP_ = ps.getUntrackedParameter<double>("MIPCut", 0.70);
84  storeRL_ = ps.getUntrackedParameter<bool>("StoreRL", false);
85  testNumber_ = ps.getUntrackedParameter<bool>("TestNumbering", false);
86 
87  muonLab[0] = "MuonRPCHits";
88  muonLab[1] = "MuonCSCHits";
89  muonLab[2] = "MuonDTHits";
90  tkHighLab[0] = "TrackerHitsPixelBarrelHighTof";
91  tkHighLab[1] = "TrackerHitsPixelEndcapHighTof";
92  tkHighLab[2] = "TrackerHitsTECHighTof";
93  tkHighLab[3] = "TrackerHitsTIBHighTof";
94  tkHighLab[4] = "TrackerHitsTIDHighTof";
95  tkHighLab[5] = "TrackerHitsTOBHighTof";
96  tkLowLab[0] = "TrackerHitsPixelBarrelLowTof";
97  tkLowLab[1] = "TrackerHitsPixelEndcapLowTof";
98  tkLowLab[2] = "TrackerHitsTECLowTof";
99  tkLowLab[3] = "TrackerHitsTIBLowTof";
100  tkLowLab[4] = "TrackerHitsTIDLowTof";
101  tkLowLab[5] = "TrackerHitsTOBLowTof";
102 
103  // register for data access
104  for (unsigned i = 0; i != 4; i++)
105  toks_calo_[i] = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hitLab[i]));
106 
107  for (unsigned i = 0; i != 3; i++)
108  toks_track_[i] = consumes<edm::PSimHitContainer>(edm::InputTag(g4Label, muonLab[i]));
109 
110  for (unsigned i = 0; i != 6; i++) {
111  toks_tkHigh_[i] = consumes<edm::PSimHitContainer>(edm::InputTag(g4Label, tkHighLab[i]));
112  toks_tkLow_[i] = consumes<edm::PSimHitContainer>(edm::InputTag(g4Label, tkLowLab[i]));
113  }
114 
115  edm::LogVerbatim("HitStudy") << "Module Label: " << g4Label << " Hits: " << hitLab[0] << ", " << hitLab[1] << ", "
116  << hitLab[2] << ", " << hitLab[3] << " MaxEnergy: " << maxEnergy_
117  << " Tmax: " << tmax_ << " MIP Cut: " << eMIP_;
118 
119  edm::Service<TFileService> tfile;
120 
121  if (!tfile.isAvailable())
122  throw cms::Exception("BadConfig") << "TFileService unavailable: "
123  << "please add it to config file";
124  char name[20], title[120];
125  sprintf(title, "Incident PT (GeV)");
126  ptInc_ = tfile->make<TH1F>("PtInc", title, 1000, 0., maxEnergy_);
127  ptInc_->GetXaxis()->SetTitle(title);
128  ptInc_->GetYaxis()->SetTitle("Events");
129  sprintf(title, "Incident Energy (GeV)");
130  eneInc_ = tfile->make<TH1F>("EneInc", title, 1000, 0., maxEnergy_);
131  eneInc_->GetXaxis()->SetTitle(title);
132  eneInc_->GetYaxis()->SetTitle("Events");
133  sprintf(title, "Incident #eta");
134  etaInc_ = tfile->make<TH1F>("EtaInc", title, 200, -5., 5.);
135  etaInc_->GetXaxis()->SetTitle(title);
136  etaInc_->GetYaxis()->SetTitle("Events");
137  sprintf(title, "Incident #phi");
138  phiInc_ = tfile->make<TH1F>("PhiInc", title, 200, -3.1415926, 3.1415926);
139  phiInc_->GetXaxis()->SetTitle(title);
140  phiInc_->GetYaxis()->SetTitle("Events");
141  std::string dets[9] = {"EB", "EB(APD)", "EB(ATJ)", "EE", "ES", "HB", "HE", "HO", "HF"};
142  for (int i = 0; i < 9; i++) {
143  sprintf(name, "Hit%d", i);
144  sprintf(title, "Number of hits in %s", dets[i].c_str());
145  hit_[i] = tfile->make<TH1F>(name, title, 1000, 0., 20000.);
146  hit_[i]->GetXaxis()->SetTitle(title);
147  hit_[i]->GetYaxis()->SetTitle("Events");
148  sprintf(name, "Time%d", i);
149  sprintf(title, "Time of the hit (ns) in %s", dets[i].c_str());
150  time_[i] = tfile->make<TH1F>(name, title, 1000, 0., 1000.);
151  time_[i]->GetXaxis()->SetTitle(title);
152  time_[i]->GetYaxis()->SetTitle("Hits");
153  sprintf(name, "TimeAll%d", i);
154  sprintf(title, "Hit time (ns) in %s with no check on Track ID", dets[i].c_str());
155  timeAll_[i] = tfile->make<TH1F>(name, title, 1000, 0., 1000.);
156  timeAll_[i]->GetXaxis()->SetTitle(title);
157  timeAll_[i]->GetYaxis()->SetTitle("Hits");
158  double ymax = maxEnergy_;
159  if (i == 1 || i == 2 || i == 4)
160  ymax = 1.0;
161  else if (i > 4 && i < 8)
162  ymax = 10.0;
163  sprintf(name, "Edep%d", i);
164  sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
165  edep_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
166  edep_[i]->GetXaxis()->SetTitle(title);
167  edep_[i]->GetYaxis()->SetTitle("Hits");
168  sprintf(name, "EdepEM%d", i);
169  sprintf(title, "Energy deposit (GeV) by EM particles in %s", dets[i].c_str());
170  edepEM_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
171  edepEM_[i]->GetXaxis()->SetTitle(title);
172  edepEM_[i]->GetYaxis()->SetTitle("Hits");
173  sprintf(name, "EdepHad%d", i);
174  sprintf(title, "Energy deposit (GeV) by hadrons in %s", dets[i].c_str());
175  edepHad_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
176  edepHad_[i]->GetXaxis()->SetTitle(title);
177  edepHad_[i]->GetYaxis()->SetTitle("Hits");
178  sprintf(name, "Etot%d", i);
179  sprintf(title, "Total energy deposit (GeV) in %s", dets[i].c_str());
180  etot_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
181  etot_[i]->GetXaxis()->SetTitle(title);
182  etot_[i]->GetYaxis()->SetTitle("Events");
183  sprintf(name, "EtotG%d", i);
184  sprintf(title, "Total energy deposit (GeV) in %s (t < 100 ns)", dets[i].c_str());
185  etotg_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
186  etotg_[i]->GetXaxis()->SetTitle(title);
187  etotg_[i]->GetYaxis()->SetTitle("Events");
188  }
189  std::string detx[9] = {"EB/EE (MIP)",
190  "HB/HE (MIP)",
191  "HB/HE/HO (MIP)",
192  "EB/EE (no MIP)",
193  "HB/HE (no MIP)",
194  "HB/HE/HO (no MIP)",
195  "EB/EE (All)",
196  "HB/HE (All)",
197  "HB/HE/HO (All)"};
198  for (int i = 0; i < 9; i++) {
199  double ymax = 1.0;
200  if (i == 0 || i == 3 || i == 6)
201  ymax = maxEnergy_;
202  sprintf(name, "EdepCal%d", i);
203  sprintf(title, "Energy deposit in %s", detx[i].c_str());
204  edepC_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
205  edepC_[i]->GetXaxis()->SetTitle(title);
206  edepC_[i]->GetYaxis()->SetTitle("Events");
207  sprintf(name, "EdepCalT%d", i);
208  sprintf(title, "Energy deposit (t < %f ns) in %s", tmax_, detx[i].c_str());
209  edepT_[i] = tfile->make<TH1F>(name, title, 5000, 0., ymax);
210  edepT_[i]->GetXaxis()->SetTitle(title);
211  edepT_[i]->GetYaxis()->SetTitle("Events");
212  }
213  hitLow = tfile->make<TH1F>("HitLow", "Number of hits in Track (Low)", 1000, 0, 10000.);
214  hitLow->GetXaxis()->SetTitle("Number of hits in Track (Low)");
215  hitLow->GetYaxis()->SetTitle("Events");
216  hitHigh = tfile->make<TH1F>("HitHigh", "Number of hits in Track (High)", 1000, 0, 10000.);
217  hitHigh->GetXaxis()->SetTitle("Number of hits in Track (High)");
218  hitHigh->GetYaxis()->SetTitle("Events");
219  hitMu = tfile->make<TH1F>("HitMu", "Number of hits in Track (Muon)", 1000, 0, 5000.);
220  hitMu->GetXaxis()->SetTitle("Number of hits in Muon");
221  hitMu->GetYaxis()->SetTitle("Events");
222  std::string dett[15] = {"Pixel Barrel (High)",
223  "Pixel Endcap (High)",
224  "TEC (High)",
225  "TIB (High)",
226  "TID (High)",
227  "TOB (High)",
228  "Pixel Barrel (Low)",
229  "Pixel Endcap (Low)",
230  "TEC (Low)",
231  "TIB (Low)",
232  "TID (Low)",
233  "TOB (Low)",
234  "RPC",
235  "CSC",
236  "DT"};
237  for (int i = 0; i < 15; i++) {
238  sprintf(name, "HitTk%d", i);
239  sprintf(title, "Number of hits in %s", dett[i].c_str());
240  hitTk_[i] = tfile->make<TH1F>(name, title, 1000, 0., 1000.);
241  hitTk_[i]->GetXaxis()->SetTitle(title);
242  hitTk_[i]->GetYaxis()->SetTitle("Events");
243  sprintf(name, "TimeTk%d", i);
244  sprintf(title, "Time of the hit (ns) in %s", dett[i].c_str());
245  tofTk_[i] = tfile->make<TH1F>(name, title, 1000, 0., 200.);
246  tofTk_[i]->GetXaxis()->SetTitle(title);
247  tofTk_[i]->GetYaxis()->SetTitle("Hits");
248  sprintf(name, "EdepTk%d", i);
249  sprintf(title, "Energy deposit (GeV) in %s", dett[i].c_str());
250  edepTk_[i] = tfile->make<TH1F>(name, title, 5000, 0., 10.);
251  edepTk_[i]->GetXaxis()->SetTitle(title);
252  edepTk_[i]->GetYaxis()->SetTitle("Hits");
253  }
254 }
255 
256 void CaloSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
257  edm::ParameterSetDescription desc;
258  desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
259  desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
260  desc.addUntracked<std::string>("EBCollection", "EcalHitsEB");
261  desc.addUntracked<std::string>("EECollection", "EcalHitsEE");
262  desc.addUntracked<std::string>("ESCollection", "EcalHitsES");
263  desc.addUntracked<std::string>("HCCollection", "HcalHits");
264  desc.addUntracked<double>("MaxEnergy", 200.0);
265  desc.addUntracked<double>("TimeCut", 100.0);
266  desc.addUntracked<double>("MIPCut", 0.70);
267  desc.addUntracked<bool>("StoreRL", false);
268  desc.addUntracked<bool>("TestNumbering", false);
269  descriptions.add("CaloSimHitStudy", desc);
270 }
271 
272 void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const&) {
273  edm::LogVerbatim("HitStudy") << "CaloSimHitStudy:Run = " << e.id().run() << " Event = " << e.id().event();
274 
275  edm::Handle<edm::HepMCProduct> EvtHandle;
276  e.getByToken(tok_evt_, EvtHandle);
277  const HepMC::GenEvent* myGenEvent = EvtHandle->GetEvent();
278 
279  double eInc = 0, etaInc = 0, phiInc = 0;
280  HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
281  if (p != myGenEvent->particles_end()) {
282  eInc = (*p)->momentum().e();
283  etaInc = (*p)->momentum().eta();
284  phiInc = (*p)->momentum().phi();
285  }
286  double ptInc = eInc / std::cosh(etaInc);
287  ptInc_->Fill(ptInc);
288  eneInc_->Fill(eInc);
289  etaInc_->Fill(etaInc);
290  phiInc_->Fill(phiInc);
291 
292  std::vector<PCaloHit> ebHits, eeHits, hcHits;
293  for (int i = 0; i < 4; i++) {
294  bool getHits = false;
295  edm::Handle<edm::PCaloHitContainer> hitsCalo;
296  e.getByToken(toks_calo_[i], hitsCalo);
297  if (hitsCalo.isValid())
298  getHits = true;
299  edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Input flags Hits " << getHits;
300 
301  if (getHits) {
302  std::vector<PCaloHit> caloHits;
303  caloHits.insert(caloHits.end(), hitsCalo->begin(), hitsCalo->end());
304  if (i == 0)
305  ebHits.insert(ebHits.end(), hitsCalo->begin(), hitsCalo->end());
306  else if (i == 1)
307  eeHits.insert(eeHits.end(), hitsCalo->begin(), hitsCalo->end());
308  else if (i == 3)
309  hcHits.insert(hcHits.end(), hitsCalo->begin(), hitsCalo->end());
310  edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Hit buffer " << caloHits.size();
311  analyzeHits(caloHits, i);
312  }
313  }
314  analyzeHits(ebHits, eeHits, hcHits);
315 
316  std::vector<PSimHit> muonHits;
317  edm::Handle<edm::PSimHitContainer> hitsTrack;
318  for (int i = 0; i < 3; i++) {
319  e.getByToken(toks_track_[i], hitsTrack);
320  if (hitsTrack.isValid()) {
321  muonHits.insert(muonHits.end(), hitsTrack->begin(), hitsTrack->end());
322  analyzeHits(hitsTrack, i + 12);
323  }
324  }
325  unsigned int nhmu = muonHits.size();
326  hitMu->Fill(double(nhmu));
327  std::vector<PSimHit> tkHighHits;
328  for (int i = 0; i < 6; i++) {
329  e.getByToken(toks_tkHigh_[i], hitsTrack);
330  if (hitsTrack.isValid()) {
331  tkHighHits.insert(tkHighHits.end(), hitsTrack->begin(), hitsTrack->end());
332  analyzeHits(hitsTrack, i);
333  }
334  }
335  unsigned int nhtkh = tkHighHits.size();
336  hitHigh->Fill(double(nhtkh));
337  std::vector<PSimHit> tkLowHits;
338  for (int i = 0; i < 6; i++) {
339  e.getByToken(toks_tkLow_[i], hitsTrack);
340  if (hitsTrack.isValid()) {
341  tkLowHits.insert(tkLowHits.end(), hitsTrack->begin(), hitsTrack->end());
342  analyzeHits(hitsTrack, i + 6);
343  }
344  }
345  unsigned int nhtkl = tkLowHits.size();
346  hitLow->Fill(double(nhtkl));
347 }
348 
349 void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
350  int nHit = hits.size();
351  int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nEB = 0, nEBAPD = 0, nEBATJ = 0, nEE = 0, nES = 0, nBad = 0;
352  std::map<unsigned int, double> hitMap;
353  std::vector<double> etot(9, 0), etotG(9, 0);
354  for (int i = 0; i < nHit; i++) {
355  double edep = hits[i].energy();
356  double time = hits[i].time();
357  unsigned int id_ = hits[i].id();
358  double edepEM = hits[i].energyEM();
359  double edepHad = hits[i].energyHad();
360  if (indx == 0) {
361  int dep = (hits[i].depth()) & PCaloHit::kEcalDepthIdMask;
362  if (dep == 1)
363  id_ |= 0x20000;
364  else if (dep == 2)
365  id_ |= 0x40000;
366  }
367  std::map<unsigned int, double>::const_iterator it = hitMap.find(id_);
368  if (it == hitMap.end()) {
369  hitMap.insert(std::pair<unsigned int, double>(id_, time));
370  }
371  int idx = -1;
372  if (indx != 3) {
373  if (indx == 0) {
374  if (storeRL_)
375  idx = 0;
376  else
377  idx = ((hits[i].depth()) & PCaloHit::kEcalDepthIdMask);
378  } else
379  idx = indx + 2;
380  time_[idx]->Fill(time);
381  edep_[idx]->Fill(edep);
382  edepEM_[idx]->Fill(edepEM);
383  edepHad_[idx]->Fill(edepHad);
384  if (idx == 0)
385  nEB++;
386  else if (idx == 1)
387  nEBAPD++;
388  else if (idx == 2)
389  nEBATJ++;
390  else if (idx == 3)
391  nEE++;
392  else if (idx == 4)
393  nES++;
394  else
395  nBad++;
396  if (indx >= 0 && indx < 3) {
397  etot[idx] += edep;
398  if (time < 100)
399  etotG[idx] += edep;
400  }
401  } else {
402  int subdet(0);
403  if (testNumber_) {
404  int ieta(0), phi(0), z(0), lay(0), depth(0);
405  HcalTestNumbering::unpackHcalIndex(id_, subdet, z, depth, ieta, phi, lay);
406  } else {
407  HcalDetId id = HcalDetId(id_);
408  subdet = id.subdet();
409  }
410  if (subdet == static_cast<int>(HcalBarrel)) {
411  idx = indx + 2;
412  nHB++;
413  } else if (subdet == static_cast<int>(HcalEndcap)) {
414  idx = indx + 3;
415  nHE++;
416  } else if (subdet == static_cast<int>(HcalOuter)) {
417  idx = indx + 4;
418  nHO++;
419  } else if (subdet == static_cast<int>(HcalForward)) {
420  idx = indx + 5;
421  nHF++;
422  }
423  if (idx > 0) {
424  time_[idx]->Fill(time);
425  edep_[idx]->Fill(edep);
426  edepEM_[idx]->Fill(edepEM);
427  edepHad_[idx]->Fill(edepHad);
428  etot[idx] += edep;
429  if (time < 100)
430  etotG[idx] += edep;
431  } else {
432  nBad++;
433  }
434  }
435  }
436  if (indx < 3) {
437  etot_[indx + 2]->Fill(etot[indx + 2]);
438  etotg_[indx + 2]->Fill(etotG[indx + 2]);
439  if (indx == 0) {
440  etot_[indx]->Fill(etot[indx]);
441  etotg_[indx]->Fill(etotG[indx]);
442  etot_[indx + 1]->Fill(etot[indx + 1]);
443  etotg_[indx + 1]->Fill(etotG[indx + 1]);
444  hit_[indx]->Fill(double(nEB));
445  hit_[indx + 1]->Fill(double(nEBAPD));
446  hit_[indx + 2]->Fill(double(nEBATJ));
447  } else {
448  hit_[indx + 2]->Fill(double(nHit));
449  }
450  } else if (indx == 3) {
451  hit_[5]->Fill(double(nHB));
452  hit_[6]->Fill(double(nHE));
453  hit_[7]->Fill(double(nHO));
454  hit_[8]->Fill(double(nHF));
455  for (int idx = 5; idx < 9; idx++) {
456  etot_[idx]->Fill(etot[idx]);
457  etotg_[idx]->Fill(etotG[idx]);
458  }
459  }
460 
461  edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::analyzeHits: EB " << nEB << ", " << nEBAPD << ", " << nEBATJ
462  << " EE " << nEE << " ES " << nES << " HB " << nHB << " HE " << nHE << " HO " << nHO
463  << " HF " << nHF << " Bad " << nBad << " All " << nHit << " Reduced " << hitMap.size();
464  std::map<unsigned int, double>::const_iterator it = hitMap.begin();
465  for (; it != hitMap.end(); it++) {
466  double time = it->second;
467  unsigned int id_ = (it->first);
468  int idx = -1;
469  if (indx < 3) {
470  if (indx == 0) {
471  if ((id_ & 0x20000) != 0)
472  idx = indx + 1;
473  else if ((id_ & 0x40000) != 0)
474  idx = indx + 1;
475  else
476  idx = indx;
477  } else {
478  idx = indx + 2;
479  }
480  if (idx >= 0 && idx < 5)
481  timeAll_[idx]->Fill(time);
482  } else if (indx == 3) {
483  int idx(-1), subdet(0);
484  if (testNumber_) {
485  int ieta(0), phi(0), z(0), lay(0), depth(0);
486  HcalTestNumbering::unpackHcalIndex(id_, subdet, z, depth, ieta, phi, lay);
487  } else {
488  HcalDetId id = HcalDetId(id_);
489  subdet = id.subdet();
490  }
491  if (subdet == static_cast<int>(HcalBarrel)) {
492  idx = indx + 2;
493  } else if (subdet == static_cast<int>(HcalEndcap)) {
494  idx = indx + 3;
495  } else if (subdet == static_cast<int>(HcalOuter)) {
496  idx = indx + 4;
497  } else if (subdet == static_cast<int>(HcalForward)) {
498  idx = indx + 5;
499  }
500  if (idx > 0) {
501  timeAll_[idx]->Fill(time);
502  }
503  }
504  }
505 }
506 
507 void CaloSimHitStudy::analyzeHits(edm::Handle<edm::PSimHitContainer>& hits, int indx) {
508  int nHit = 0;
509  edm::PSimHitContainer::const_iterator ihit;
510  std::string label(" ");
511  if (indx >= 0 && indx < 6)
512  label = tkHighLab[indx];
513  else if (indx >= 6 && indx < 12)
514  label = tkLowLab[indx - 6];
515  else if (indx >= 12 && indx < 15)
516  label = muonLab[indx - 12];
517  for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
518  edepTk_[indx]->Fill(ihit->energyLoss());
519  tofTk_[indx]->Fill(ihit->timeOfFlight());
520  nHit++;
521  }
522  hitTk_[indx]->Fill(float(nHit));
523  edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::analyzeHits: for " << label << " Index " << indx << " # of Hits "
524  << nHit;
525 }
526 
527 void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& ebHits,
528  std::vector<PCaloHit>& eeHits,
529  std::vector<PCaloHit>& hcHits) {
530  double edepEB = 0, edepEBT = 0;
531  for (unsigned int i = 0; i < ebHits.size(); i++) {
532  double edep = ebHits[i].energy();
533  double time = ebHits[i].time();
534  if (((ebHits[i].depth()) & PCaloHit::kEcalDepthIdMask) == 0) {
535  edepEB += edep;
536  if (time < tmax_)
537  edepEBT += edep;
538  }
539  }
540  double edepEE = 0, edepEET = 0;
541  for (unsigned int i = 0; i < eeHits.size(); i++) {
542  double edep = eeHits[i].energy();
543  double time = eeHits[i].time();
544  edepEE += edep;
545  if (time < tmax_)
546  edepEET += edep;
547  }
548  double edepH = 0, edepHT = 0, edepHO = 0, edepHOT = 0;
549  for (unsigned int i = 0; i < hcHits.size(); i++) {
550  double edep = hcHits[i].energy();
551  double time = hcHits[i].time();
552  int subdet(0);
553  if (testNumber_) {
554  int ieta(0), phi(0), z(0), lay(0), depth(0);
555  HcalTestNumbering::unpackHcalIndex(hcHits[i].id(), subdet, z, depth, ieta, phi, lay);
556  } else {
557  HcalDetId id = HcalDetId(hcHits[i].id());
558  subdet = id.subdet();
559  }
560  if (subdet == static_cast<int>(HcalBarrel) || subdet == static_cast<int>(HcalEndcap)) {
561  edepH += edep;
562  if (time < tmax_)
563  edepHT += edep;
564  } else if (subdet == static_cast<int>(HcalOuter)) {
565  edepHO += edep;
566  if (time < tmax_)
567  edepHOT += edep;
568  }
569  }
570  double edepE = edepEB + edepEE;
571  double edepET = edepEBT + edepEET;
572  double edepHC = edepH + edepHO;
573  double edepHCT = edepHT + edepHOT;
574  edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::energy in EB " << edepEB << " (" << edepEBT << ") from "
575  << ebHits.size() << " hits; "
576  << " energy in EE " << edepEE << " (" << edepEET << ") from " << eeHits.size()
577  << " hits; energy in HC " << edepH << ", " << edepHO << " (" << edepHT << ", " << edepHOT
578  << ") from " << hcHits.size() << " hits";
579 
580  edepC_[6]->Fill(edepE);
581  edepT_[6]->Fill(edepET);
582  edepC_[7]->Fill(edepH);
583  edepT_[7]->Fill(edepHT);
584  edepC_[8]->Fill(edepHC);
585  edepT_[8]->Fill(edepHCT);
586  if (edepE < eMIP_) {
587  edepC_[0]->Fill(edepE);
588  edepC_[1]->Fill(edepH);
589  edepC_[2]->Fill(edepHC);
590  } else {
591  edepC_[3]->Fill(edepE);
592  edepC_[4]->Fill(edepH);
593  edepC_[5]->Fill(edepHC);
594  }
595  if (edepET < eMIP_) {
596  edepT_[0]->Fill(edepET);
597  edepT_[1]->Fill(edepHT);
598  edepT_[2]->Fill(edepHCT);
599  } else {
600  edepT_[3]->Fill(edepET);
601  edepT_[4]->Fill(edepHT);
602  edepT_[5]->Fill(edepHCT);
603  }
604 }
605 
606 //define this as a plug-in
607 DEFINE_FWK_MODULE(CaloSimHitStudy);
ConfigurationDescriptions.h
CaloSimHitStudy::etaInc_
TH1F * etaInc_
Definition: CaloSimHitStudy.cc:66
runGCPTkAlMap.title
string title
Definition: runGCPTkAlMap.py:94
EDAnalyzer.h
mps_fire.i
i
Definition: mps_fire.py:428
MessageLogger.h
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
ESHandle.h
edm::Run
Definition: Run.h:45
edm::EDGetTokenT< edm::HepMCProduct >
CaloSimHitStudy::tmax_
double tmax_
Definition: CaloSimHitStudy.cc:61
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
PSimHitContainer.h
CaloSimHitStudy::g4Label
std::string g4Label
Definition: CaloSimHitStudy.cc:54
CaloSimHitStudy::tok_evt_
edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
Definition: CaloSimHitStudy.cc:55
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
CaloSimHitStudy::toks_tkLow_
edm::EDGetTokenT< edm::PSimHitContainer > toks_tkLow_[6]
Definition: CaloSimHitStudy.cc:59
CaloSimHitStudy::tkLowLab
std::string tkLowLab[6]
Definition: CaloSimHitStudy.cc:60
CaloSimHitStudy::CaloSimHitStudy
CaloSimHitStudy(const edm::ParameterSet &ps)
Definition: CaloSimHitStudy.cc:70
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
CaloHitID.h
HcalTestNumbering.h
CaloSimHitStudy::~CaloSimHitStudy
~CaloSimHitStudy() override
Definition: CaloSimHitStudy.cc:40
edm::one::EDAnalyzer
Definition: EDAnalyzer.h:30
CaloSimHitStudy::hitHigh
TH1F * hitHigh
Definition: CaloSimHitStudy.cc:65
HcalBarrel
Definition: HcalAssistant.h:33
edm::Handle< edm::PSimHitContainer >
HepMC::GenEvent
Definition: hepmc_rootio.cc:9
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
CaloSimHitStudy::edepHad_
TH1F * edepHad_[9]
Definition: CaloSimHitStudy.cc:64
CaloSimHitStudy::hitLow
TH1F * hitLow
Definition: CaloSimHitStudy.cc:66
CaloSimHitStudy::toks_track_
edm::EDGetTokenT< edm::PSimHitContainer > toks_track_[3]
Definition: CaloSimHitStudy.cc:57
MakerMacros.h
CaloSimHitStudy::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: CaloSimHitStudy.cc:272
CaloSimHitStudy::tkHighLab
std::string tkHighLab[6]
Definition: CaloSimHitStudy.cc:60
PSimHit.h
HcalTestNumbering::unpackHcalIndex
static void unpackHcalIndex(const uint32_t &idx, int &det, int &z, int &depth, int &eta, int &phi, int &lay)
Definition: HcalTestNumbering.cc:18
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
L1TOccupancyClient_cfi.ymax
ymax
Definition: L1TOccupancyClient_cfi.py:43
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
CaloSimHitStudy::muonLab
std::string muonLab[3]
Definition: CaloSimHitStudy.cc:60
CaloSimHitStudy::endRun
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: CaloSimHitStudy.cc:47
Service.h
tfile
Definition: tfile.py:1
CaloSimHitStudy::edepT_
TH1F * edepT_[9]
Definition: CaloSimHitStudy.cc:67
DDAxes::z
HcalOuter
Definition: HcalAssistant.h:35
CaloSimHitStudy::timeAll_
TH1F * timeAll_[9]
Definition: CaloSimHitStudy.cc:65
CaloSimHitStudy::hitMu
TH1F * hitMu
Definition: CaloSimHitStudy.cc:65
LEDCalibrationChannels.depth
depth
Definition: LEDCalibrationChannels.py:65
CaloSimHitStudy::beginRun
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: CaloSimHitStudy.cc:46
CaloSimHitStudy::ptInc_
TH1F * ptInc_
Definition: CaloSimHitStudy.cc:66
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TFileService.h
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
CaloSimHitStudy::etotg_
TH1F * etotg_[9]
Definition: CaloSimHitStudy.cc:65
CaloSimHitStudy::edepC_
TH1F * edepC_[9]
Definition: CaloSimHitStudy.cc:67
edm::ParameterSet
Definition: ParameterSet.h:47
CaloSimHitStudy::hitLab
std::string hitLab[4]
Definition: CaloSimHitStudy.cc:54
Event.h
CaloSimHitStudy::hitTk_
TH1F * hitTk_[15]
Definition: CaloSimHitStudy.cc:67
CaloSimHitStudy::edepEM_
TH1F * edepEM_[9]
Definition: CaloSimHitStudy.cc:64
HcalDetId.h
PCaloHit.h
CaloSimHitStudy::etot_
TH1F * etot_[9]
Definition: CaloSimHitStudy.cc:65
HcalDetId
Definition: HcalDetId.h:12
edm::Service< TFileService >
CaloSimHitStudy::edepTk_
TH1F * edepTk_[15]
Definition: CaloSimHitStudy.cc:67
CaloSimHitStudy::eMIP_
double eMIP_
Definition: CaloSimHitStudy.cc:61
CaloSimHitStudy::eneInc_
TH1F * eneInc_
Definition: CaloSimHitStudy.cc:66
edm::EventSetup
Definition: EventSetup.h:57
edm::HepMCProduct::GetEvent
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:34
CaloSimHitStudy::analyzeHits
void analyzeHits(std::vector< PCaloHit > &, int)
Definition: CaloSimHitStudy.cc:349
CaloSimHitStudy::toks_tkHigh_
edm::EDGetTokenT< edm::PSimHitContainer > toks_tkHigh_[6]
Definition: CaloSimHitStudy.cc:58
InputTag.h
compare.tfile
tfile
Definition: compare.py:325
HcalForward
Definition: HcalAssistant.h:36
DDAxes::phi
CaloSimHitStudy::time_
TH1F * time_[9]
Definition: CaloSimHitStudy.cc:64
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
HcalEndcap
Definition: HcalAssistant.h:34
CaloSimHitStudy::beginJob
void beginJob() override
Definition: CaloSimHitStudy.cc:44
Frameworkfwd.h
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
CaloSimHitStudy::edep_
TH1F * edep_[9]
Definition: CaloSimHitStudy.cc:64
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
CaloSimHitStudy::tofTk_
TH1F * tofTk_[15]
Definition: CaloSimHitStudy.cc:67
CaloSimHitStudy::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: CaloSimHitStudy.cc:256
CaloSimHitStudy::storeRL_
bool storeRL_
Definition: CaloSimHitStudy.cc:62
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
Exception.h
PCaloHitContainer.h
CaloSimHitStudy::testNumber_
bool testNumber_
Definition: CaloSimHitStudy.cc:62
cms::Exception
Definition: Exception.h:70
CaloSimHitStudy
Definition: CaloSimHitStudy.cc:37
ParameterSet.h
HepMCProduct.h
CaloSimHitStudy::phiInc_
TH1F * phiInc_
Definition: CaloSimHitStudy.cc:66
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
ntuplemaker.time
time
Definition: ntuplemaker.py:310
edm::Event
Definition: Event.h:73
CaloSimHitStudy::hit_
TH1F * hit_[9]
Definition: CaloSimHitStudy.cc:64
edm::InputTag
Definition: InputTag.h:15
label
const char * label
Definition: PFTauDecayModeTools.cc:11
CaloSimHitStudy::toks_calo_
edm::EDGetTokenT< edm::PCaloHitContainer > toks_calo_[4]
Definition: CaloSimHitStudy.cc:56
PCaloHit::kEcalDepthIdMask
static const int kEcalDepthIdMask
Definition: PCaloHit.h:60
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16