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