CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:38
TFileService::kSharedResource
static const std::string kSharedResource
Definition: TFileService.h:76
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
CaloSimHitStudy::toks_tkHigh_
edm::EDGetTokenT< edm::PSimHitContainer > toks_tkHigh_[6]
Definition: CaloSimHitStudy.cc:57
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:40
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
mps_fire.i
i
Definition: mps_fire.py:428
CaloSimHitStudy::g4Label
std::string g4Label
Definition: CaloSimHitStudy.cc:53
edm::ParameterSetDescription::addUntracked
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:100
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
CaloSimHitStudy::analyzeHits
void analyzeHits(std::vector< PCaloHit > &, int)
Definition: CaloSimHitStudy.cc:348
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
SiStrip_OfflineMonitoring_cff.ymax
tuple ymax
Definition: SiStrip_OfflineMonitoring_cff.py:53
CaloSimHitStudy::tok_evt_
edm::EDGetTokenT< edm::HepMCProduct > tok_evt_
Definition: CaloSimHitStudy.cc:54
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
CaloSimHitStudy::tkLowLab
std::string tkLowLab[6]
Definition: CaloSimHitStudy.cc:59
CaloSimHitStudy::beginRun
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: CaloSimHitStudy.cc:45
label
char const * label
Definition: PFTauDecayModeTools.cc:11
CaloSimHitStudy::toks_tkLow_
edm::EDGetTokenT< edm::PSimHitContainer > toks_tkLow_[6]
Definition: CaloSimHitStudy.cc:58
edm::Service::isAvailable
bool isAvailable() const
Definition: Service.h:40
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
CaloSimHitStudy::edepTk_
TH1F * edepTk_[15]
Definition: CaloSimHitStudy.cc:66
CaloSimHitStudy::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: CaloSimHitStudy.cc:255
CaloSimHitStudy::CaloSimHitStudy
CaloSimHitStudy(const edm::ParameterSet &ps)
Definition: CaloSimHitStudy.cc:69
CaloSimHitStudy::tkHighLab
std::string tkHighLab[6]
Definition: CaloSimHitStudy.cc:59
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
HLT_FULL_cff.depth
tuple depth
Definition: HLT_FULL_cff.py:15256
CaloSimHitStudy::beginJob
void beginJob() override
Definition: CaloSimHitStudy.cc:43
CaloSimHitStudy::toks_track_
edm::EDGetTokenT< edm::PSimHitContainer > toks_track_[3]
Definition: CaloSimHitStudy.cc:56
CaloSimHitStudy::muonLab
std::string muonLab[3]
Definition: CaloSimHitStudy.cc:59
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
CaloSimHitStudy::toks_calo_
edm::EDGetTokenT< edm::PCaloHitContainer > toks_calo_[4]
Definition: CaloSimHitStudy.cc:55
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
CaloSimHitStudy::~CaloSimHitStudy
~CaloSimHitStudy() override
Definition: CaloSimHitStudy.cc:39
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:86956
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
CaloSimHitStudy::endRun
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: CaloSimHitStudy.cc:46
CaloSimHitStudy::hitLab
std::string hitLab[4]
Definition: CaloSimHitStudy.cc:53
compare.tfile
tuple tfile
Definition: compare.py:324
CaloSimHitStudy::analyze
void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: CaloSimHitStudy.cc:271
edm::Run
Definition: Run.h:45
PCaloHit::kEcalDepthIdMask
static const int kEcalDepthIdMask
Definition: PCaloHit.h:60
Generated for CMSSW Reference Manual by doxygen 1.8.5