CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HGCalTBAnalyzer.cc
Go to the documentation of this file.
1 // system include files
2 #include <cmath>
3 #include <fstream>
4 #include <iostream>
5 #include <map>
6 #include <memory>
7 #include <string>
8 #include <vector>
9 
10 // user include files
15 
26 
40 
42 
43 // Root objects
44 #include "TH1.h"
45 #include "TH2.h"
46 #include "TProfile.h"
47 #include "TProfile2D.h"
48 #include "TROOT.h"
49 #include "TSystem.h"
50 #include "TTree.h"
51 
52 //#define EDM_ML_DEBUG
53 
54 class HGCalTBAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
55 public:
56  explicit HGCalTBAnalyzer(edm::ParameterSet const&);
57  ~HGCalTBAnalyzer() override;
58 
59  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
60 
61 private:
62  void beginJob() override;
63  void beginRun(edm::Run const&, edm::EventSetup const&) override;
64  void endRun(edm::Run const&, edm::EventSetup const&) override {}
65  void analyze(edm::Event const&, edm::EventSetup const&) override;
66  void analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront);
69  template <class T1>
70  void analyzeDigi(int type, const T1& detId, uint16_t adc);
72  void analyzePassiveHits(edm::Handle<edm::PassiveHitContainer> const& hgcPh, int subdet);
73  static bool sortTime(const std::pair<double, double>& i, const std::pair<double, double>& j);
74 
76  std::unique_ptr<AHCalGeometry> ahcalGeom_;
78  const HGCalGeometry* hgeom_[2];
79  const bool ifEE_, ifFH_, ifBH_, ifBeam_;
81  const bool doTree_, doTreeCell_;
85  const double zFrontEE_, zFrontFH_, zFrontBH_;
86  const int sampleIndex_;
88  std::vector<int> idBeams_;
95 
109 
110  TTree* tree_;
111  TH1D *hSimHitE_[4], *hSimHitT_[4];
112  TH1D *hDigiADC_[3], *hDigiLng_[2];
113  TH1D *hRecHitE_[3], *hSimHitEn_[4], *hBeam_;
114  TH2D *hDigiOcc_[3], *hRecHitOcc_[3];
115  TProfile *hSimHitLng_[3], *hSimHitLng1_[3];
116  TProfile* hSimHitLng2_[3];
117  TProfile *hRecHitLng_[3], *hRecHitLng1_[3];
118  TProfile2D *hSimHitLat_[3], *hRecHitLat_[3];
119  std::vector<TH1D*> hSimHitLayEn1EE_, hSimHitLayEn2EE_;
120  std::vector<TH1D*> hSimHitLayEn1FH_, hSimHitLayEn2FH_;
121  std::vector<TH1D*> hSimHitLayEn1BH_, hSimHitLayEn2BH_;
122  std::vector<TH1D*> hSimHitLayEnBeam_;
123  std::vector<float> simHitLayEn1EE_, simHitLayEn2EE_;
124  std::vector<float> simHitLayEn1FH_, simHitLayEn2FH_;
125  std::vector<float> simHitLayEn1BH_, simHitLayEn2BH_;
126  std::vector<float> simHitLayEnBeam_;
127  std::vector<uint32_t> simHitCellIdEE_, simHitCellIdFH_;
128  std::vector<uint32_t> simHitCellIdBH_, simHitCellIdBeam_;
129  std::vector<float> simHitCellEnEE_, simHitCellEnFH_;
130  std::vector<float> simHitCellEnBH_, simHitCellEnBeam_;
135 
137  std::vector<float> hgcPassiveBeamEnergy_;
139  std::vector<std::string> hgcPassiveBeamName_;
141 
144  int nBeamMC_;
145  std::vector<int> pdgIdBeamMC_;
146  std::vector<float> xBeamMC_, yBeamMC_, zBeamMC_;
147  std::vector<float> pxBeamMC_, pyBeamMC_, pzBeamMC_, pBeamMC_;
148 };
149 
151  : ifEE_(iConfig.getParameter<bool>("useEE")),
152  ifFH_(iConfig.getParameter<bool>("useFH")),
153  ifBH_(iConfig.getParameter<bool>("useBH")),
154  ifBeam_(iConfig.getParameter<bool>("useBeam")),
155  doSimHits_(iConfig.getParameter<bool>("doSimHits")),
156  doDigis_(iConfig.getParameter<bool>("doDigis")),
157  doRecHits_(iConfig.getParameter<bool>("doRecHits")),
158  doTree_(iConfig.getParameter<bool>("doTree")),
159  doTreeCell_(iConfig.getParameter<bool>("doTreeCell")),
160  doPassive_(iConfig.getParameter<bool>("doPassive")),
161  doPassiveEE_(iConfig.getParameter<bool>("doPassiveEE")),
162  doPassiveHE_(iConfig.getParameter<bool>("doPassiveHE")),
163  doPassiveBH_(iConfig.getParameter<bool>("doPassiveBH")),
164  addP_(iConfig.getParameter<bool>("addP")),
165  doBeam_(iConfig.getParameter<bool>("doBeam")),
166  detectorEE_(iConfig.getParameter<std::string>("detectorEE")),
167  detectorFH_(iConfig.getParameter<std::string>("detectorFH")),
168  detectorBH_(iConfig.getParameter<std::string>("detectorBH")),
169  detectorBeam_(iConfig.getParameter<std::string>("detectorBeam")),
170  zFrontEE_(iConfig.getParameter<double>("zFrontEE")),
171  zFrontFH_(iConfig.getParameter<double>("zFrontFH")),
172  zFrontBH_(iConfig.getParameter<double>("zFrontBH")),
173  sampleIndex_(iConfig.getParameter<int>("sampleIndex")),
174  gev2mip200_(iConfig.getUntrackedParameter<double>("gev2mip200", 57.0e-6)),
175  gev2mip300_(iConfig.getUntrackedParameter<double>("gev2mip300", 85.5e-6)),
176  stoc_smear_time_200_(iConfig.getUntrackedParameter<double>("stoc_smear_time_200", 10.24)),
177  stoc_smear_time_300_(iConfig.getUntrackedParameter<double>("stoc_smear_time_300", 15.5)),
178  labelGen_(iConfig.getParameter<edm::InputTag>("generatorSrc")),
179  labelHitEE_(iConfig.getParameter<std::string>("caloHitSrcEE")),
180  labelHitFH_(iConfig.getParameter<std::string>("caloHitSrcFH")),
181  labelHitBH_(iConfig.getParameter<std::string>("caloHitSrcBH")),
182  labelHitBeam_(iConfig.getParameter<std::string>("caloHitSrcBeam")),
183  labelDigiEE_(iConfig.getParameter<edm::InputTag>("digiSrcEE")),
184  labelDigiFH_(iConfig.getParameter<edm::InputTag>("digiSrcFH")),
185  labelDigiBH_(iConfig.getParameter<edm::InputTag>("digiSrcBH")),
186  labelRHitEE_(iConfig.getParameter<edm::InputTag>("recHitSrcEE")),
187  labelRHitFH_(iConfig.getParameter<edm::InputTag>("recHitSrcFH")),
188  labelRHitBH_(iConfig.getParameter<edm::InputTag>("recHitSrcBH")),
189  labelPassiveEE_(iConfig.getParameter<edm::InputTag>("passiveEE")),
190  labelPassiveFH_(iConfig.getParameter<edm::InputTag>("passiveFH")),
191  labelPassiveBH_(iConfig.getParameter<edm::InputTag>("passiveBH")),
192  labelPassiveCMSE_(iConfig.getParameter<edm::InputTag>("passiveCMSE")),
193  labelPassiveBeam_(iConfig.getParameter<edm::InputTag>("passiveBeam")),
194  tok_hepMC_(consumes<edm::HepMCProduct>(labelGen_)),
195  tok_simTk_(consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"))),
196  tok_simVtx_(consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"))),
197  tok_hitsEE_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitEE_))),
198  tok_hitsFH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitFH_))),
199  tok_hitsBH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBH_))),
200  tok_hitsBeam_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBeam_))),
201  tok_digiEE_(consumes<HGCalDigiCollection>(labelDigiEE_)),
202  tok_digiFH_(consumes<HGCalDigiCollection>(labelDigiFH_)),
203  tok_digiBH_(consumes<HGCalDigiCollection>(labelDigiFH_)),
204  tok_hitrEE_(consumes<HGCRecHitCollection>(labelRHitEE_)),
205  tok_hitrFH_(consumes<HGCRecHitCollection>(labelRHitFH_)),
206  tok_hitrBH_(consumes<HGCRecHitCollection>(labelRHitBH_)),
207  tok_hgcPHEE_(consumes<edm::PassiveHitContainer>(labelPassiveEE_)),
208  tok_hgcPHFH_(consumes<edm::PassiveHitContainer>(labelPassiveFH_)),
209  tok_hgcPHBH_(consumes<edm::PassiveHitContainer>(labelPassiveBH_)),
210  tok_hgcPHCMSE_(consumes<edm::PassiveHitContainer>(labelPassiveCMSE_)),
211  tok_hgcPHBeam_(consumes<edm::PassiveHitContainer>(labelPassiveBeam_)),
212  tokDDDEE_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
213  edm::ESInputTag("", detectorEE_))),
214  tokGeomEE_(
215  esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorEE_))),
216  tokDDDFH_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
217  edm::ESInputTag("", detectorFH_))),
218  tokGeomFH_(
219  esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorFH_))) {
220  usesResource("TFileService");
221  ahcalGeom_ = std::make_unique<AHCalGeometry>(iConfig);
222 
223  // now do whatever initialization is needed
229  idBeams_ = (iConfig.getParameter<std::vector<int>>("idBeams"));
230 #ifdef EDM_ML_DEBUG
231  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: SimHits = " << doSimHits_ << " Digis = " << doDigis_ << ":"
232  << sampleIndex_ << " RecHits = " << doRecHits_ << " useDets " << ifEE_ << ":" << ifFH_
233  << ":" << ifBH_ << ":" << ifBeam_ << " zFront " << zFrontEE_ << ":" << zFrontFH_ << ":"
234  << zFrontBH_ << " IdBeam " << idBeams_.size() << ":";
235  for (unsigned int k = 0; k < idBeams_.size(); ++k)
236  edm::LogVerbatim("HGCSim") << " [" << k << "] " << idBeams_[k];
237  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: DoPassive " << doPassive_ << ":" << doPassiveEE_ << ":"
238  << doPassiveHE_ << ":" << doPassiveBH_;
239  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: MIP conversion factors " << gev2mip200_ << ":" << gev2mip300_
240  << " Time smearing " << stoc_smear_time_200_ << ":" << stoc_smear_time_300_ << " AddP "
241  << addP_;
242 #endif
243  if (idBeams_.empty())
244  idBeams_.push_back(1001);
245 
246 #ifdef EDM_ML_DEBUG
247  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: GeneratorSource = " << labelGen_;
248  if (ifEE_)
249  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorEE_ << " with tags " << labelHitEE_ << ", "
250  << labelDigiEE_ << ", " << labelRHitEE_;
251  if (ifFH_)
252  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorFH_ << " with tags " << labelHitFH_ << ", "
253  << labelDigiFH_ << ", " << labelRHitFH_;
254  if (ifBH_)
255  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorBH_ << " with tags " << labelHitBH_ << ", "
256  << labelDigiFH_ << ", " << labelRHitBH_;
257  if (ifBeam_)
258  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorBeam_ << " with tags " << labelHitBeam_;
259 #endif
260 }
261 
263 
266  desc.add<std::string>("detectorEE", "HGCalEESensitive");
267  desc.add<bool>("useEE", true);
268  desc.add<double>("zFrontEE", 0.0);
269  desc.add<std::string>("caloHitSrcEE", "HGCHitsEE");
270  desc.add<edm::InputTag>("digiSrcEE", edm::InputTag("hgcalDigis", "EE"));
271  desc.add<edm::InputTag>("recHitSrcEE", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
272  desc.add<std::string>("detectorFH", "HGCalHESiliconSensitive");
273  desc.add<bool>("useFH", false);
274  desc.add<double>("zFrontFH", 0.0);
275  desc.add<std::string>("caloHitSrcFH", "HGCHitsHEfront");
276  desc.add<edm::InputTag>("digiSrcFH", edm::InputTag("hgcalDigis", "HEfront"));
277  desc.add<edm::InputTag>("recHitSrcFH", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
278  desc.add<std::string>("detectorBH", "AHCal");
279  desc.add<bool>("useBH", false);
280  desc.add<double>("zFrontBH", 0.0);
281  desc.add<std::string>("caloHitSrcBH", "HcalHits");
282  desc.add<edm::InputTag>("digiSrcBH", edm::InputTag("hgcalDigis", "HEback"));
283  desc.add<edm::InputTag>("recHitSrcBH", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
284  desc.add<std::string>("detectorBeam", "HcalTB06BeamDetector");
285  desc.add<bool>("useBeam", false);
286  desc.add<std::string>("caloHitSrcBeam", "HcalTB06BeamHits");
287  std::vector<int> ids = {
288  1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1011, 1012, 1013, 1014, 2001, 2002, 2003, 2004, 2005};
289  desc.add<std::vector<int>>("idBeams", ids);
290  desc.add<edm::InputTag>("generatorSrc", edm::InputTag("generatorSmeared"));
291  desc.add<edm::InputTag>("passiveEE", edm::InputTag("g4SimHits", "HGCalEEPassiveHits"));
292  desc.add<edm::InputTag>("passiveFH", edm::InputTag("g4SimHits", "HGCalHEPassiveHits"));
293  desc.add<edm::InputTag>("passiveBH", edm::InputTag("g4SimHits", "HGCalAHPassiveHits"));
294  desc.add<edm::InputTag>("passiveCMSE", edm::InputTag("g4SimHits", "CMSEPassiveHits"));
295  desc.add<edm::InputTag>("passiveBeam", edm::InputTag("g4SimHits", "HGCalBeamPassiveHits"));
296 
297  desc.add<bool>("doSimHits", true);
298  desc.add<bool>("doDigis", true);
299  desc.add<bool>("doRecHits", true);
300  desc.add<int>("sampleIndex", 0);
301  desc.add<bool>("doTree", true);
302  desc.add<bool>("doTreeCell", true);
303  desc.add<bool>("doPassive", false);
304  desc.add<bool>("doPassiveEE", false);
305  desc.add<bool>("doPassiveHE", false);
306  desc.add<bool>("doPassiveBH", false);
307  desc.add<bool>("addP", false);
308  desc.add<bool>("doBeam", false);
309  desc.addUntracked<double>("gev2mip200", 57.0e-6);
310  desc.addUntracked<double>("gev2mip300", 85.5e-6);
311  desc.addUntracked<double>("stoc_smear_time_200", 10.24);
312  desc.addUntracked<double>("stoc_smear_time_300", 15.5);
313  desc.addUntracked<int>("maxDepth", 12);
314  desc.addUntracked<double>("deltaX", 30.0); // Size of tile along X
315  desc.addUntracked<double>("deltaY", 30.0); // Size of tile along Y
316  desc.addUntracked<double>("deltaZ", 81.0); // Thickness of a single layer
317  desc.addUntracked<double>("zFirst", 17.6); // Position of the center
318  descriptions.add("HGCalTBAnalyzer", desc);
319 }
320 
322  char name[40], title[100];
323  hBeam_ = fs_->make<TH1D>("BeamP", "Beam Momentum", 1000, 0, 1000.0);
324  for (int i = 0; i < 3; ++i) {
325  bool book(ifEE_);
327  if (i == 1) {
328  book = ifFH_;
329  det = detectorFH_;
330  } else if (i == 2) {
331  book = ifBH_;
332  det = detectorBH_;
333  }
334 
335  if (doSimHits_ && book) {
336  sprintf(name, "SimHitEn%s", det.c_str());
337  sprintf(title, "Sim Hit Energy for %s", det.c_str());
338  hSimHitE_[i] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
339  sprintf(name, "SimHitEnX%s", det.c_str());
340  sprintf(title, "Sim Hit Energy for %s", det.c_str());
341  hSimHitEn_[i] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
342  sprintf(name, "SimHitTm%s", det.c_str());
343  sprintf(title, "Sim Hit Timing for %s", det.c_str());
344  hSimHitT_[i] = fs_->make<TH1D>(name, title, 5000, 0., 500.0);
345  sprintf(name, "SimHitLat%s", det.c_str());
346  sprintf(title, "Lateral Shower profile (Sim Hit) for %s", det.c_str());
347  hSimHitLat_[i] = fs_->make<TProfile2D>(name, title, 100, -100., 100., 100, -100., 100.);
348  sprintf(name, "SimHitLng%s", det.c_str());
349  sprintf(title, "Longitudinal Shower profile (Sim Hit) for %s", det.c_str());
350  hSimHitLng_[i] = fs_->make<TProfile>(name, title, 50, 0., 100.);
351  sprintf(name, "SimHitLng1%s", det.c_str());
352  sprintf(title, "Longitudinal Shower profile (Layer) for %s", det.c_str());
353  hSimHitLng1_[i] = fs_->make<TProfile>(name, title, 200, 0., 100.);
354  sprintf(name, "SimHitLng2%s", det.c_str());
355  sprintf(title, "Longitudinal Shower profile (Layer) for %s", det.c_str());
356  hSimHitLng2_[i] = fs_->make<TProfile>(name, title, 200, 0., 100.);
357  }
358 
359  if (doDigis_ && book) {
360  sprintf(name, "DigiADC%s", det.c_str());
361  sprintf(title, "ADC at Digi level for %s", det.c_str());
362  hDigiADC_[i] = fs_->make<TH1D>(name, title, 100, 0., 100.0);
363  sprintf(name, "DigiOcc%s", det.c_str());
364  sprintf(title, "Occupancy (Digi)for %s", det.c_str());
365  hDigiOcc_[i] = fs_->make<TH2D>(name, title, 100, -10., 10., 100, -10., 10.);
366  sprintf(name, "DigiLng%s", det.c_str());
367  sprintf(title, "Longitudinal Shower profile (Digi) for %s", det.c_str());
368  hDigiLng_[i] = fs_->make<TH1D>(name, title, 100, 0., 10.);
369  }
370 
371  if (doRecHits_ && book) {
372  sprintf(name, "RecHitEn%s", det.c_str());
373  sprintf(title, "Rec Hit Energy for %s", det.c_str());
374  hRecHitE_[i] = fs_->make<TH1D>(name, title, 1000, 0., 10.0);
375  sprintf(name, "RecHitOcc%s", det.c_str());
376  sprintf(title, "Occupancy (Rec Hit)for %s", det.c_str());
377  hRecHitOcc_[i] = fs_->make<TH2D>(name, title, 100, -10., 10., 100, -10., 10.);
378  sprintf(name, "RecHitLat%s", det.c_str());
379  sprintf(title, "Lateral Shower profile (Rec Hit) for %s", det.c_str());
380  hRecHitLat_[i] = fs_->make<TProfile2D>(name, title, 100, -10., 10., 100, -10., 10.);
381  sprintf(name, "RecHitLng%s", det.c_str());
382  sprintf(title, "Longitudinal Shower profile (Rec Hit) for %s", det.c_str());
383  hRecHitLng_[i] = fs_->make<TProfile>(name, title, 100, 0., 10.);
384  sprintf(name, "RecHitLng1%s", det.c_str());
385  sprintf(title, "Longitudinal Shower profile vs Layer for %s", det.c_str());
386  hRecHitLng1_[i] = fs_->make<TProfile>(name, title, 120, 0., 60.);
387  }
388  }
389  if (ifBeam_ && doSimHits_) {
390  sprintf(name, "SimHitEn%s", detectorBeam_.c_str());
391  sprintf(title, "Sim Hit Energy for %s", detectorBeam_.c_str());
392  hSimHitE_[3] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
393  sprintf(name, "SimHitEnX%s", detectorBeam_.c_str());
394  sprintf(title, "Sim Hit Energy for %s", detectorBeam_.c_str());
395  hSimHitEn_[3] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
396  sprintf(name, "SimHitTm%s", detectorBeam_.c_str());
397  sprintf(title, "Sim Hit Timing for %s", detectorBeam_.c_str());
398  hSimHitT_[3] = fs_->make<TH1D>(name, title, 5000, 0., 500.0);
399  }
400  if (doSimHits_ && doTree_) {
401  tree_ = fs_->make<TTree>("HGCTB", "SimHitEnergy");
402  tree_->Branch("simHitLayEn1EE", &simHitLayEn1EE_);
403  tree_->Branch("simHitLayEn2EE", &simHitLayEn2EE_);
404  tree_->Branch("simHitLayEn1FH", &simHitLayEn1FH_);
405  tree_->Branch("simHitLayEn2FH", &simHitLayEn2FH_);
406  tree_->Branch("simHitLayEn1BH", &simHitLayEn1BH_);
407  tree_->Branch("simHitLayEn2BH", &simHitLayEn2BH_);
408  tree_->Branch("xBeam", &xBeam_, "xBeam/D");
409  tree_->Branch("yBeam", &yBeam_, "yBeam/D");
410  tree_->Branch("zBeam", &zBeam_, "zBeam/D");
411  tree_->Branch("pBeam", &pBeam_, "pBeam/D");
412  tree_->Branch("thetaBeam", &thetaBeam_, "thetaBeam/D");
413  tree_->Branch("phiBeam", &phiBeam_, "phiBeam/D");
414  if (doBeam_) {
415  tree_->Branch("nBeamMC", &nBeamMC_, "nBeamMC/I");
416  tree_->Branch("pdgIdBeamMC", &pdgIdBeamMC_);
417  tree_->Branch("xBeamMC", &xBeamMC_);
418  tree_->Branch("yBeamMC", &yBeamMC_);
419  tree_->Branch("zBeamMC", &zBeamMC_);
420  tree_->Branch("pxBeamMC", &pxBeamMC_);
421  tree_->Branch("pyBeamMC", &pyBeamMC_);
422  tree_->Branch("pzBeamMC", &pzBeamMC_);
423  tree_->Branch("pBeamMC", &pBeamMC_);
424  }
425  if (doTreeCell_) {
426  tree_->Branch("simHitCellIdEE", &simHitCellIdEE_);
427  tree_->Branch("simHitCellEnEE", &simHitCellEnEE_);
428  tree_->Branch("simHitCellIdFH", &simHitCellIdFH_);
429  tree_->Branch("simHitCellEnFH", &simHitCellEnFH_);
430  tree_->Branch("simHitCellIdBH", &simHitCellIdBH_);
431  tree_->Branch("simHitCellEnBH", &simHitCellEnBH_);
432  tree_->Branch("simHitCellIdBeam", &simHitCellIdBeam_);
433  tree_->Branch("simHitCellEnBeam", &simHitCellEnBeam_);
434 
435  tree_->Branch("simHitCellColBH", &simHitCellColBH_);
436  tree_->Branch("simHitCellRowBH", &simHitCellRowBH_);
437  tree_->Branch("simHitCellLayerBH", &simHitCellLayerBH_);
438  tree_->Branch("simHitCellTimeFirstHitEE", &simHitCellTimeFirstHitEE_);
439  tree_->Branch("simHitCellTimeFirstHitFH", &simHitCellTimeFirstHitFH_);
440  //tree_->Branch("simHitCellTimeFirstHitBH", &simHitCellTimeFirstHitBH_);
441  tree_->Branch("simHitCellTime15MipEE", &simHitCellTime15MipEE_);
442  tree_->Branch("simHitCellTime15MipFH", &simHitCellTime15MipFH_);
443  //tree_->Branch("simHitCellTime15MipBH", &simHitCellTime15MipBH_);
444  tree_->Branch("simHitCellTimeLastHitEE", &simHitCellTimeLastHitEE_);
445  tree_->Branch("simHitCellTimeLastHitFH", &simHitCellTimeLastHitFH_);
446  //tree_->Branch("simHitCellTimeLastHitBH", &simHitCellTimeLastHitBH_);
447  }
448  }
449 
450  if (doPassive_ && doTree_) {
451  if (doPassiveEE_) {
452  tree_->Branch("hgcPassiveEEEnergy", &hgcPassiveEEEnergy_);
453  tree_->Branch("hgcPassiveEEName", &hgcPassiveEEName_);
454  tree_->Branch("hgcPassiveEEID", &hgcPassiveEEID_);
455  }
456  if (doPassiveHE_) {
457  tree_->Branch("hgcPassiveFHEnergy", &hgcPassiveFHEnergy_);
458  tree_->Branch("hgcPassiveFHName", &hgcPassiveFHName_);
459  tree_->Branch("hgcPassiveFHID", &hgcPassiveFHID_);
460  }
461  if (doPassiveBH_) {
462  tree_->Branch("hgcPassiveBHEnergy", &hgcPassiveBHEnergy_);
463  tree_->Branch("hgcPassiveBHName", &hgcPassiveBHName_);
464  tree_->Branch("hgcPassiveBHID", &hgcPassiveBHID_);
465  }
466  tree_->Branch("hgcPassiveCMSEEnergy", &hgcPassiveCMSEEnergy_);
467  tree_->Branch("hgcPassiveCMSEName", &hgcPassiveCMSEName_);
468  tree_->Branch("hgcPassiveCMSEID", &hgcPassiveCMSEID_);
469  tree_->Branch("hgcPassiveBeamEnergy", &hgcPassiveBeamEnergy_);
470  tree_->Branch("hgcPassiveBeamName", &hgcPassiveBeamName_);
471  tree_->Branch("hgcPassiveBeamID", &hgcPassiveBeamID_);
472  }
473 }
474 
476  char name[40], title[100];
477  if (ifEE_) {
478  hgcons_[0] = &iSetup.getData(tokDDDEE_);
479  if (doDigis_ || doRecHits_) {
480  hgeom_[0] = &iSetup.getData(tokGeomEE_);
481  } else {
482  hgeom_[0] = nullptr;
483  }
484  for (unsigned int l = 0; l < hgcons_[0]->layers(false); ++l) {
485  sprintf(name, "SimHitEnA%d%s", l, detectorEE_.c_str());
486  sprintf(title, "Sim Hit Energy in SIM layer %d for %s", l + 1, detectorEE_.c_str());
487  hSimHitLayEn1EE_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
488  if (l % 3 == 0) {
489  sprintf(name, "SimHitEnB%d%s", (l / 3 + 1), detectorEE_.c_str());
490  sprintf(title, "Sim Hit Energy in layer %d for %s", (l / 3 + 1), detectorEE_.c_str());
491  hSimHitLayEn2EE_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
492  }
493  }
494 #ifdef EDM_ML_DEBUG
495  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::" << detectorEE_ << " defined with " << hgcons_[0]->layers(false)
496  << " layers";
497 #endif
498  } else {
499  hgcons_[0] = nullptr;
500  hgeom_[0] = nullptr;
501  }
502 
503  if (ifFH_) {
504  hgcons_[1] = &iSetup.getData(tokDDDFH_);
505  if (doDigis_ || doRecHits_) {
506  hgeom_[1] = &iSetup.getData(tokGeomFH_);
507  } else {
508  hgeom_[1] = nullptr;
509  }
510  for (unsigned int l = 0; l < hgcons_[1]->layers(false); ++l) {
511  sprintf(name, "SimHitEnA%d%s", l, detectorFH_.c_str());
512  sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorFH_.c_str());
513  hSimHitLayEn1FH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
514  if (l % 3 == 0) {
515  sprintf(name, "SimHitEnB%d%s", (l / 3 + 1), detectorFH_.c_str());
516  sprintf(title, "Sim Hit Energy in layer %d for %s", (l / 3 + 1), detectorFH_.c_str());
517  hSimHitLayEn2FH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
518  }
519  }
520 #ifdef EDM_ML_DEBUG
521  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::" << detectorFH_ << " defined with " << hgcons_[1]->layers(false)
522  << " layers";
523 #endif
524  } else {
525  hgcons_[1] = nullptr;
526  hgeom_[1] = nullptr;
527  }
528 
529  if (ifBH_) {
530  for (int l = 0; l < ahcalGeom_->maxDepth(); ++l) {
531  sprintf(name, "SimHitEnA%d%s", l, detectorBH_.c_str());
532  sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorBH_.c_str());
533  hSimHitLayEn1BH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
534  sprintf(name, "SimHitEnB%d%s", l, detectorBH_.c_str());
535  sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorBH_.c_str());
536  hSimHitLayEn2BH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
537  }
538  }
539 
540  if (ifBeam_) {
541  for (unsigned int l = 0; l < idBeams_.size(); ++l) {
542  sprintf(name, "SimHitEna%d%s", l, detectorBeam_.c_str());
543  sprintf(title, "Sim Hit Energy in type %d for %s", idBeams_[l], detectorBeam_.c_str());
544  hSimHitLayEnBeam_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
545  }
546  }
547 }
548 
550  // Generator input
551  const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_hepMC_);
552  if (!evtMC.isValid()) {
553  edm::LogWarning("HGCal") << "no HepMCProduct found";
554  } else {
555  const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
556  unsigned int k(0);
557  HepMC::FourVector pxyz(0, 0, 0, 0);
558  for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
559  ++p, ++k) {
560  edm::LogVerbatim("HGCSim") << "Particle [" << k << "] with p " << (*p)->momentum().rho() << " theta "
561  << (*p)->momentum().theta() << " phi " << (*p)->momentum().phi() << " pxyz ("
562  << (*p)->momentum().px() << ", " << (*p)->momentum().py() << ", "
563  << (*p)->momentum().pz() << ")";
564  if (addP_) {
565  pxyz.setPx(pxyz.px() + (*p)->momentum().px());
566  pxyz.setPy(pxyz.py() + (*p)->momentum().py());
567  pxyz.setPz(pxyz.pz() + (*p)->momentum().pz());
568  pxyz.setE(pxyz.e() + (*p)->momentum().e());
569  } else if (!addP_ && (k == 0)) {
570  pxyz = (*p)->momentum();
571  }
572  }
573  hBeam_->Fill(pxyz.rho());
574  edm::LogVerbatim("HGCSim") << "Particle with p " << pxyz.rho() << " theta " << pxyz.theta() << " phi "
575  << pxyz.phi();
576  }
577 
578  // Now the Simhits
579  if (doSimHits_) {
582  analyzeSimTracks(SimTk, SimVtx);
583 
584  simHitLayEn1EE_.clear();
585  simHitLayEn2EE_.clear();
586  simHitLayEn1FH_.clear();
587  simHitLayEn2FH_.clear();
588  simHitLayEn1BH_.clear();
589  simHitLayEn2BH_.clear();
590  simHitLayEnBeam_.clear();
591  simHitCellIdEE_.clear();
592  simHitCellEnEE_.clear();
593  simHitCellIdFH_.clear();
594  simHitCellEnFH_.clear();
595  simHitCellIdBH_.clear();
596  simHitCellEnBH_.clear();
597  simHitCellIdBeam_.clear();
598  simHitCellEnBeam_.clear();
599  simHitCellColBH_.clear();
600  simHitCellRowBH_.clear();
601  simHitCellLayerBH_.clear();
603  simHitCellTime15MipEE_.clear();
604  simHitCellTimeLastHitEE_.clear();
606  simHitCellTime15MipFH_.clear();
607  simHitCellTimeLastHitFH_.clear();
608  std::vector<PCaloHit> caloHits;
609  if (ifEE_) {
610  simHitLayEn1EE_ = std::vector<float>(hgcons_[0]->layers(false), 0);
611  simHitLayEn2EE_ = std::vector<float>(hgcons_[0]->layers(true), 0);
612  const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsEE_);
613  if (theCaloHitContainers.isValid()) {
614 #ifdef EDM_ML_DEBUG
615  edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorEE_ << " has " << theCaloHitContainers->size()
616  << " hits";
617 #endif
618  caloHits.clear();
619  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
620  analyzeSimHits(0, caloHits, zFrontEE_);
621  } else {
622 #ifdef EDM_ML_DEBUG
623  edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorEE_ << " !!!";
624 #endif
625  }
626  }
627  if (ifFH_) {
628  simHitLayEn1FH_ = std::vector<float>(hgcons_[1]->layers(false), 0);
629  simHitLayEn2FH_ = std::vector<float>(hgcons_[1]->layers(true), 0);
630  const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsFH_);
631  if (theCaloHitContainers.isValid()) {
632 #ifdef EDM_ML_DEBUG
633  edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorFH_ << " has " << theCaloHitContainers->size()
634  << " hits";
635 #endif
636  caloHits.clear();
637  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
638  analyzeSimHits(1, caloHits, zFrontFH_);
639  } else {
640 #ifdef EDM_ML_DEBUG
641  edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorFH_ << " !!!";
642 #endif
643  }
644  }
645  if (ifBH_) {
646  simHitLayEn1BH_ = std::vector<float>(ahcalGeom_->maxDepth(), 0);
647  simHitLayEn2BH_ = std::vector<float>(ahcalGeom_->maxDepth(), 0);
648  const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsBH_);
649  if (theCaloHitContainers.isValid()) {
650 #ifdef EDM_ML_DEBUG
651  edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBH_ << " has " << theCaloHitContainers->size()
652  << " hits";
653 #endif
654  caloHits.clear();
655  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
656  analyzeSimHits(2, caloHits, zFrontBH_);
657  } else {
658 #ifdef EDM_ML_DEBUG
659  edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBH_ << " !!!";
660 #endif
661  }
662  }
663  if (ifBeam_) {
664  simHitLayEnBeam_ = std::vector<float>(idBeams_.size(), 0);
665  const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsBeam_);
666  if (theCaloHitContainers.isValid()) {
667 #ifdef EDM_ML_DEBUG
668  edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBeam_ << " has "
669  << theCaloHitContainers->size() << " hits";
670 #endif
671  caloHits.clear();
672  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
673  analyzeSimHits(3, caloHits, 0.0);
674  } else {
675 #ifdef EDM_ML_DEBUG
676  edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBeam_ << " !!!";
677 #endif
678  }
679  }
680  } // if (doSimHits_)
681 
683  if (doPassive_) {
685  hgcPassiveEEEnergy_.clear();
686  hgcPassiveEEName_.clear();
687  hgcPassiveEEID_.clear();
688  if (doPassiveEE_) {
690  analyzePassiveHits(hgcPHEE, 1);
691  }
693  hgcPassiveFHEnergy_.clear();
694  hgcPassiveFHName_.clear();
695  hgcPassiveFHID_.clear();
696  if (doPassiveHE_) {
698  analyzePassiveHits(hgcPHFH, 2);
699  }
701  hgcPassiveBHEnergy_.clear();
702  hgcPassiveBHName_.clear();
703  hgcPassiveBHID_.clear();
704  if (doPassiveBH_) {
706  analyzePassiveHits(hgcPHBH, 3);
707  }
709  hgcPassiveCMSEEnergy_.clear();
710  hgcPassiveCMSEName_.clear();
711  hgcPassiveCMSEID_.clear();
713  analyzePassiveHits(hgcPHCMSE, 4);
715  hgcPassiveBeamName_.clear();
716  hgcPassiveBeamEnergy_.clear();
717  hgcPassiveBeamID_.clear();
719  analyzePassiveHits(hgcPHBeam, 5);
720  }
721 
722  if ((doSimHits_ || doPassive_) && (doTree_))
723  tree_->Fill();
724 
725  // Now the Digis
726  if (doDigis_) {
727  if (ifEE_) {
728  const edm::Handle<HGCalDigiCollection>& theDigiContainers = iEvent.getHandle(tok_digiEE_);
729  if (theDigiContainers.isValid()) {
730 #ifdef EDM_ML_DEBUG
731  edm::LogVerbatim("HGCSim") << "HGCDigiCintainer for " << detectorEE_ << " with " << theDigiContainers->size()
732  << " element(s)";
733 #endif
734  for (const auto& it : *theDigiContainers) {
735  HGCalDetId detId = (it.id());
736  const HGCSample& hgcSample = it.sample(sampleIndex_);
737  uint16_t adc = hgcSample.data();
738  analyzeDigi(0, detId, adc);
739  }
740  }
741  }
742  if (ifFH_) {
743  const edm::Handle<HGCalDigiCollection>& theDigiContainers = iEvent.getHandle(tok_digiFH_);
744  if (theDigiContainers.isValid()) {
745 #ifdef EDM_ML_DEBUG
746  edm::LogVerbatim("HGCSim") << "HGCDigiContainer for " << detectorFH_ << " with " << theDigiContainers->size()
747  << " element(s)";
748 #endif
749  for (const auto& it : *theDigiContainers) {
750  HGCalDetId detId = (it.id());
751  const HGCSample& hgcSample = it.sample(sampleIndex_);
752  uint16_t adc = hgcSample.data();
753  analyzeDigi(1, detId, adc);
754  }
755  }
756  }
757  }
758 
759  // The Rechits
760  if (doRecHits_) {
761  if (ifEE_) {
762  const edm::Handle<HGCRecHitCollection>& theCaloHitContainers = iEvent.getHandle(tok_hitrEE_);
763  if (theCaloHitContainers.isValid()) {
764 #ifdef EDM_ML_DEBUG
765  edm::LogVerbatim("HGCSim") << "HGCRecHitCollection for " << detectorEE_ << " has "
766  << theCaloHitContainers->size() << " hits";
767 #endif
768  analyzeRecHits(0, theCaloHitContainers);
769  } else {
770 #ifdef EDM_ML_DEBUG
771  edm::LogVerbatim("HGCSim") << "HGCRecHitCollection does not exist for " << detectorEE_ << " !!!";
772 #endif
773  }
774  }
775  if (ifFH_) {
776  const edm::Handle<HGCRecHitCollection>& theCaloHitContainers = iEvent.getHandle(tok_hitrFH_);
777  if (theCaloHitContainers.isValid()) {
778 #ifdef EDM_ML_DEBUG
779  edm::LogVerbatim("HGCSim") << "HGCRecHitCollection for " << detectorFH_ << " has "
780  << theCaloHitContainers->size() << " hits";
781 #endif
782  analyzeRecHits(1, theCaloHitContainers);
783  } else {
784 #ifdef EDM_ML_DEBUG
785  edm::LogVerbatim("HGCSim") << "HGCRecHitCollection does not exist for " << detectorFH_ << " !!!";
786 #endif
787  } // else
788  } // if (ifFH_)
789  } // if (doRecHits_)
790 
791 } // void HGCalTBAnalyzer::analyze
792 
793 void HGCalTBAnalyzer::analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront) {
794  std::map<uint32_t, double> map_hits, map_hitn;
795  std::map<uint32_t, double> map_hittime_firsthit, map_hittime_lasthit, map_hittime_15Mip;
796  std::map<int, double> map_hitDepth, map_hitWafer;
797  std::map<int, std::pair<uint32_t, double>> map_hitLayer, map_hitCell;
798  double entot(0);
799  std::map<uint32_t, double> nhits;
800  std::map<uint32_t, int> ID, Depth;
801  std::map<uint32_t, double> GeV2Mip;
802  std::map<uint32_t, double> StochTermTime;
803  std::vector<int> nSimLayers;
805  std::map<uint32_t, std::vector<std::pair<double, double>>> map_hitTimeEn;
806  //bool debug = true;
807  bool debug = false;
808  for (unsigned int i = 0; i < hits.size(); i++) {
809  double energy = hits[i].energy();
810  double time = hits[i].time();
811  uint32_t id = hits[i].id();
812  entot += energy;
813  int subdet, zside, layer, sector, subsector(0), cell, depth(0), idx(0);
814  if (type == 2) {
815  subdet = HcalDetId(id).subdet();
816  if (subdet != HcalOther)
817  continue;
818  AHCalDetId hid(id);
819  layer = depth = hid.depth();
820  zside = hid.zside();
821  sector = hid.irow();
822  cell = hid.icol();
823  idx = ((hid.irowAbs() * 100) + (hid.icolAbs()));
824  if (debug)
825  edm::LogVerbatim("HGCSim") << "depth, sector, cell " << depth << ":" << sector << ":" << cell;
826  } else if (type == 3) {
827  HcalTestBeamNumbering::unpackIndex(id, subdet, layer, sector, cell);
828  depth = layer;
829  zside = 1;
830  idx = subdet * 1000 + layer;
831  layer = idx;
832  } else {
833  HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, subsector, cell);
834  depth = hgcons_[type]->simToReco(cell, layer, sector, true).second;
835  idx = sector * 1000 + cell;
836 #ifdef EDM_ML_DEBUG
837  std::pair<float, float> xy = hgcons_[type]->locateCell(cell, layer, sector, false);
838  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::detId " << std::hex << id << std::dec << " Layer:Wafer:Cell "
839  << layer << ":" << sector << ":" << cell << " Position " << xy.first << ":"
840  << xy.second << ":" << hgcons_[type]->waferZ(layer, false);
841 #endif
842  }
843 #ifdef EDM_ML_DEBUG
844  edm::LogVerbatim("HGCSim") << "SimHit:Hit[" << i << "] Id " << subdet << ":" << zside << ":" << layer << ":"
845  << sector << ":" << subsector << ":" << cell << ":" << depth << " Energy " << energy
846  << " Time " << time;
847 #endif
848  if (map_hits.count(id) != 0) {
849  map_hits[id] += energy;
850  } else {
851  map_hits[id] = energy;
852  }
853  if (map_hitLayer.count(layer) != 0) {
854  double ee = energy + map_hitLayer[layer].second;
855  map_hitLayer[layer] = std::make_pair(id, ee);
856  } else {
857  map_hitLayer[layer] = std::make_pair(id, energy);
858  }
859  if (map_hitWafer.count(sector) != 0)
860  map_hitWafer[sector] += energy;
861  else
862  map_hitWafer[sector] = energy;
863  if (depth >= 0) {
864  if (map_hitCell.count(idx) != 0) {
865  double ee = energy + map_hitCell[idx].second;
866  map_hitCell[idx] = std::make_pair(id, ee);
867  } else {
868  map_hitCell[idx] = std::make_pair(id, energy);
869  }
870  if (debug) {
871  if (type == 1)
872  edm::LogVerbatim("HGCSim") << "EE, depth is and map_hitDepth[depth] " << depth << " " << map_hitDepth[depth];
873  if (type == 2)
874  edm::LogVerbatim("HGCSim") << "BH, depth is and map_hitDepth[depth] " << depth << " " << map_hitDepth[depth];
875  }
876  if (map_hitDepth.count(depth) != 0) {
877  map_hitDepth[depth] += energy;
878  } else {
879  map_hitDepth[depth] = energy;
880  }
881  uint32_t idn =
882  (type >= 2) ? id : HGCalTestNumbering::packHexagonIndex(subdet, zside, depth, sector, subsector, cell);
883  map_hitTimeEn[idn].push_back(std::make_pair(time, energy));
884  GeV2Mip[idn] = gev2mip300_;
885  StochTermTime[idn] = stoc_smear_time_300_;
886  ID[idn] = id;
887  Depth[idn] = depth;
888  if (map_hitn.count(idn) != 0) {
889  map_hitn[idn] += energy;
890  ++nhits[idn];
891  } else {
892  map_hitn[idn] = energy;
893  nhits[idn] = 1;
894  }
895  }
896  hSimHitT_[type]->Fill(time, energy);
897  }
898 
899  if (type < 2) { //store only for EE and FH
900  edm::LogVerbatim("HGCSim") << "HGCalTAnalyzer:: " << map_hitWafer.size() << " wafers are hit in type " << type;
901  for (auto itr = map_hitWafer.begin(); itr != map_hitWafer.end(); ++itr)
902  edm::LogVerbatim("HGCSim") << "Wafer: " << itr->first << " Deposited Energy " << itr->second;
904  for (const auto& itr : map_hitTimeEn) {
905  uint32_t id = itr.first;
906  int wafer = -99;
908  wafer = HGCalDetId(ID[id]).wafer();
909  double layer = HGCalDetId(id).layer();
910  double thickness = hgcons_[type]->cellThickness(layer, wafer, 0);
911  if (debug)
912  edm::LogVerbatim("HGCSim") << "wafer is : depth (reco) " << wafer << " " << Depth[id]
913  << "\ntype : layer : wafer thickness " << type << " " << layer << " " << thickness
914  << "\nID(sim) and id(reco) " << std::hex << ID[id] << " " << id << std::dec;
915  if (thickness == 300) {
916  GeV2Mip[id] = gev2mip300_;
917  StochTermTime[id] = stoc_smear_time_300_;
918  } else if (thickness == 200) {
919  GeV2Mip[id] = gev2mip200_;
920  StochTermTime[id] = stoc_smear_time_200_;
921  }
922  //first sort
923  std::sort(map_hitTimeEn[id].begin(), map_hitTimeEn[id].end(), sortTime);
925  double threshold = 15.;
926  double totE = 0.;
927  double totEbeforeThreshold = 0.;
928  double timebeforeThreshold = 0.;
929  double timeAtThresohld = 0.;
930  for (unsigned int ihit = 0; ihit < map_hitTimeEn[id].size(); ihit++) {
931  double energy = (map_hitTimeEn[id].at(ihit)).second / GeV2Mip[id];
932  totE += energy;
933  double time = (map_hitTimeEn[id].at(ihit)).first;
934  if (debug)
935  edm::LogVerbatim("HGCSim") << "Tot E till now : time of that E : GeV2Mip[id] is " << totE << " " << time
936  << " " << GeV2Mip[id];
937  if (totE < threshold) {
938  totEbeforeThreshold = totE;
939  timebeforeThreshold = time;
940  } else {
941  timeAtThresohld =
942  (threshold - totEbeforeThreshold) * (time - timebeforeThreshold) / (totE - totEbeforeThreshold) +
943  timebeforeThreshold;
944  map_hittime_15Mip[id] = timeAtThresohld;
945  if (debug)
946  edm::LogVerbatim("HGCSim") << "ihit : energyBefore : timeBefore : energyTot : timeTot : timeAt15MIP "
947  << ihit << " " << totEbeforeThreshold << " " << timebeforeThreshold << " "
948  << totE << " " << time << " " << map_hittime_15Mip[id];
949  break;
950  }
951  }
952  if (!map_hitTimeEn[id].empty()) {
953  map_hittime_firsthit[id] = (map_hitTimeEn[id].at(0)).first;
954  map_hittime_lasthit[id] = (map_hitTimeEn[id].at(map_hitTimeEn[id].size() - 1)).first;
955  if (map_hittime_15Mip[id] < map_hittime_firsthit[id])
956  map_hittime_15Mip[id] = map_hittime_firsthit[id];
957  /*
959  double firsthit_sm = ran3.Gaus(map_hittime_firsthit_vtxCorr[id], StochTermTime[id]);
960  double lasthit_sm = ran3.Gaus(map_hittime_lasthit_vtxCorr[id], StochTermTime[id]);
961  double threshmiphit_sm = ran3.Gaus(map_hittime_15Mip_vtxCorr[id], StochTermTime[id]);
962  */
963  }
964 
965  if (totE < threshold)
966  map_hittime_15Mip[id] = -99;
967  if (debug)
968  edm::LogVerbatim("HGCSim") << "id : first hit time : last hit time " << id << " " << map_hittime_firsthit[id]
969  << " " << map_hittime_lasthit[id] << "\nFinally for this cell, time is "
970  << map_hittime_15Mip[id];
971  }
972  }
973 
974  hSimHitEn_[type]->Fill(entot);
975  for (const auto& itr : map_hits) {
976  hSimHitE_[type]->Fill(itr.second);
977  }
978 
979  if (debug)
980  edm::LogVerbatim("HGCSim") << "Now looping over map_hitLayer";
981  for (const auto& itr : map_hitLayer) {
982  int layer = (type == 2) ? itr.first : (itr.first - 1);
983  double energy = (itr.second).second;
984  double zp(0);
985  if (type < 2)
986  zp = hgcons_[type]->waferZ(layer + 1, false);
987  else if (type == 2)
988  zp = ahcalGeom_->getZ(AHCalDetId((itr.second).first));
989 #ifdef EDM_ML_DEBUG
990  edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer + 1 << " Z " << zp << ":" << zp - zFront << " E " << energy;
991 #endif
992  if (type < 3) {
993  hSimHitLng_[type]->Fill(zp - zFront, energy);
994  hSimHitLng2_[type]->Fill(layer + 1, energy);
995  }
996  if (type == 0) {
997  if (layer < static_cast<int>(hSimHitLayEn1EE_.size())) {
999  hSimHitLayEn1EE_[layer]->Fill(energy);
1000  }
1001  } else if (type == 1) {
1002  if (layer < static_cast<int>(hSimHitLayEn1FH_.size())) {
1004  hSimHitLayEn1FH_[layer]->Fill(energy);
1005  }
1006  } else if (type == 2) {
1007  if (debug)
1008  edm::LogVerbatim("HGCSim") << "layer < hSimHitLayEn1BH_.size()";
1009  if (layer < static_cast<int>(hSimHitLayEn1BH_.size())) {
1011  hSimHitLayEn1BH_[layer]->Fill(energy);
1012  }
1013  } else {
1014  for (unsigned int k = 0; k < idBeams_.size(); ++k) {
1015  if (layer + 1 == idBeams_[k]) {
1017  hSimHitLayEnBeam_[k]->Fill(energy);
1018  break;
1019  }
1020  }
1021  }
1022  }
1023  for (const auto& itr : map_hitDepth) {
1024  int layer = (type == 2) ? itr.first : (itr.first - 1);
1025  double energy = itr.second;
1026 #ifdef EDM_ML_DEBUG
1027  edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer + 1 << " " << energy;
1028 #endif
1029  hSimHitLng1_[type]->Fill(layer + 1, energy);
1030  if (type == 0) {
1031  if (layer < static_cast<int>(hSimHitLayEn2EE_.size())) {
1033  hSimHitLayEn2EE_[layer]->Fill(energy);
1034  }
1035  } else if (type == 1) {
1036  if (layer < static_cast<int>(hSimHitLayEn2FH_.size())) {
1038  hSimHitLayEn2FH_[layer]->Fill(energy);
1039  }
1040  } else if (type == 2) {
1041  if (debug)
1042  edm::LogVerbatim("HGCSim") << "Inside map_hitDepth, layer no. " << layer;
1043  if (layer < static_cast<int>(hSimHitLayEn2BH_.size())) {
1045  hSimHitLayEn2BH_[layer]->Fill(energy);
1046  }
1047  }
1048  }
1049 
1050  if (debug)
1051  edm::LogVerbatim("HGCSim") << "Now looping over map_hitCell";
1052  if (type < 3) {
1053  for (const auto& itr : map_hitCell) {
1054  uint32_t id = ((itr.second).first);
1055  double energy = ((itr.second).second);
1056  std::pair<float, float> xy(0, 0);
1057  double xx(0);
1058  if (type == 2) {
1059  xy = ahcalGeom_->getXY(AHCalDetId(id));
1060  xx = xy.first;
1061  } else {
1062  int subdet, zside, layer, sector, subsector, cell;
1063  HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, subsector, cell);
1064  xy = hgcons_[type]->locateCell(cell, layer, sector, false);
1065  double zp = hgcons_[type]->waferZ(layer, false);
1066  xx = (zp < 0) ? -xy.first : xy.first;
1067  }
1068  hSimHitLat_[type]->Fill(xx, xy.second, energy);
1069  }
1070  }
1071 
1072  for (const auto& itr : map_hitn) {
1073  uint32_t id = itr.first;
1074  double energy = itr.second;
1075  if (type == 0) {
1076  double time_firsthit = map_hittime_firsthit[id];
1077  double time15Mip = map_hittime_15Mip[id];
1078  double time_lasthit = map_hittime_lasthit[id];
1079  simHitCellIdEE_.push_back(id);
1080  simHitCellEnEE_.push_back(energy);
1081  simHitCellTimeFirstHitEE_.push_back(time_firsthit);
1082  simHitCellTime15MipEE_.push_back(time15Mip);
1083  simHitCellTimeLastHitEE_.push_back(time_lasthit);
1084  if (debug && (energy / GeV2Mip[id] < 15) && (map_hittime_15Mip[id] > 0))
1085  edm::LogVerbatim("HGCSim") << "FOUND!!!!rechit energy : Finally for this cell, time is " << energy / GeV2Mip[id]
1086  << " " << map_hittime_15Mip[id];
1087  } else if (type == 1) {
1088  double time_firsthit = map_hittime_firsthit[id];
1089  double time15Mip = map_hittime_15Mip[id];
1090  double time_lasthit = map_hittime_lasthit[id];
1091  simHitCellIdFH_.push_back(id);
1092  simHitCellEnFH_.push_back(energy);
1093  simHitCellTimeFirstHitFH_.push_back(time_firsthit);
1094  simHitCellTime15MipFH_.push_back(time15Mip);
1095  simHitCellTimeLastHitFH_.push_back(time_lasthit);
1096  } else if (type == 2) {
1097  simHitCellIdBH_.push_back(id);
1098  simHitCellEnBH_.push_back(energy);
1099  AHCalDetId hid(id);
1100  int row = hid.irow();
1101  int col = hid.icol();
1102  int layer = hid.depth();
1103  simHitCellColBH_.push_back(col);
1104  simHitCellRowBH_.push_back(row);
1105  simHitCellLayerBH_.push_back(layer);
1106 #ifdef EDM_ML_DEBUG
1107  edm::LogVerbatim("HGCSim") << "ID: " << std::hex << id << std::dec << " Layer: " << layer << " col: " << col
1108  << " row: " << row;
1109 #endif
1110  } else if (type == 3) {
1111  simHitCellIdBeam_.push_back(id);
1112  simHitCellEnBeam_.push_back(energy);
1113  }
1114  }
1115 }
1116 
1118  edm::Handle<edm::SimVertexContainer> const& SimVtx) {
1119  xBeam_ = yBeam_ = zBeam_ = pBeam_ = -9999;
1120  nBeamMC_ = thetaBeam_ = phiBeam_ = -9999;
1121  int nParBeam = 0;
1122  int vertIndex(-1);
1123  if (doBeam_) {
1124  pdgIdBeamMC_.clear();
1125  xBeamMC_.clear();
1126  yBeamMC_.clear();
1127  zBeamMC_.clear();
1128  pxBeamMC_.clear();
1129  pyBeamMC_.clear();
1130  pzBeamMC_.clear();
1131  pBeamMC_.clear();
1132  }
1133  std::vector<float> verX, verY, verZ;
1134  verX.clear();
1135  verY.clear();
1136  verZ.clear();
1137  for (const auto& simVtxItr : *SimVtx) {
1138  verX.push_back(simVtxItr.position().X());
1139  verY.push_back(simVtxItr.position().Y());
1140  verZ.push_back(simVtxItr.position().Z());
1141  }
1142 #ifdef EDM_ML_DEBUG
1143  edm::LogVerbatim("HGCSim") << "Size of track " << SimTk->size();
1144 #endif
1145  HepMC::FourVector pxyz(0, 0, 0, 0);
1146  for (const auto& simTrkItr : *SimTk) {
1147  if (addP_ && !(simTrkItr.noGenpart())) {
1148  pxyz.setPx(pxyz.px() + simTrkItr.momentum().px());
1149  pxyz.setPy(pxyz.py() + simTrkItr.momentum().py());
1150  pxyz.setPz(pxyz.pz() + simTrkItr.momentum().pz());
1151  pxyz.setE(pxyz.e() + simTrkItr.momentum().e());
1152 #ifdef EDM_ML_DEBUG
1153  edm::LogVerbatim("HGCSim") << "Track " << simTrkItr.trackId() << " Vertex " << simTrkItr.vertIndex() << " Type "
1154  << simTrkItr.type() << " Charge " << simTrkItr.charge() << " px "
1155  << simTrkItr.momentum().px() << " py " << simTrkItr.momentum().py() << " pz "
1156  << simTrkItr.momentum().pz() << " P " << simTrkItr.momentum().P() << " GenIndex "
1157  << simTrkItr.genpartIndex();
1158  edm::LogVerbatim("HGCSim") << "Vertex " << simTrkItr.vertIndex()
1159  << " position-> X: " << verX[simTrkItr.vertIndex()]
1160  << " Y: " << verY[simTrkItr.vertIndex()] << " Z: " << verZ[simTrkItr.vertIndex()];
1161 #endif
1162  }
1163  if (doBeam_ && !(simTrkItr.noGenpart())) {
1164  nParBeam++;
1165  pdgIdBeamMC_.push_back(simTrkItr.type());
1166  xBeamMC_.push_back(verX[simTrkItr.vertIndex()]);
1167  yBeamMC_.push_back(verY[simTrkItr.vertIndex()]);
1168  zBeamMC_.push_back(verZ[simTrkItr.vertIndex()]);
1169  pxBeamMC_.push_back(simTrkItr.momentum().px());
1170  pyBeamMC_.push_back(simTrkItr.momentum().py());
1171  pzBeamMC_.push_back(simTrkItr.momentum().pz());
1172  pBeamMC_.push_back(simTrkItr.momentum().P());
1173  } else if (!addP_ && (vertIndex == -1)) {
1174  pxyz = simTrkItr.momentum();
1175  }
1176  if (vertIndex == -1)
1177  vertIndex = simTrkItr.vertIndex();
1178  }
1179  nBeamMC_ = nParBeam;
1180  pBeam_ = pxyz.rho();
1181  thetaBeam_ = pxyz.theta();
1182  phiBeam_ = pxyz.phi();
1183  if (phiBeam_ < 0)
1184  phiBeam_ += (2 * M_PI);
1185  if (vertIndex != -1 && vertIndex < static_cast<int>(SimVtx->size())) {
1186  edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
1187  for (int iv = 0; iv < vertIndex; iv++)
1188  simVtxItr++;
1189  edm::LogVerbatim("HGCSim") << "Vertex " << vertIndex << " position " << simVtxItr->position();
1190  xBeam_ = verX[0];
1191  yBeam_ = verY[0];
1192  zBeam_ = verZ[0];
1193  }
1194 }
1195 
1196 template <class T1>
1197 void HGCalTBAnalyzer::analyzeDigi(int type, const T1& detId, uint16_t adc) {
1198  DetId id1 = DetId(detId.rawId());
1199  GlobalPoint global = hgeom_[type]->getPosition(id1);
1200  hDigiOcc_[type]->Fill(global.x(), global.y());
1201  hDigiLng_[type]->Fill(global.z());
1202  hDigiADC_[type]->Fill(adc);
1203 }
1204 
1206  std::map<int, double> map_hitLayer;
1207  std::map<int, std::pair<DetId, double>> map_hitCell;
1208  for (const auto& it : *hits) {
1209  DetId detId = it.id();
1210  GlobalPoint global = hgeom_[type]->getPosition(detId);
1211  double energy = it.energy();
1212  int layer = HGCalDetId(detId).layer();
1213  int cell = HGCalDetId(detId).cell();
1214 #ifdef EDM_ML_DEBUG
1215  edm::LogVerbatim("HGCSim") << "Layer thickness " << hgcons_[type]->waferTypeL(HGCalDetId(detId).wafer());
1216 #endif
1217  hRecHitOcc_[type]->Fill(global.x(), global.y(), energy);
1218  hRecHitE_[type]->Fill(energy);
1219  if (map_hitLayer.count(layer) != 0) {
1220  map_hitLayer[layer] += energy;
1221  } else {
1222  map_hitLayer[layer] = energy;
1223  }
1224  if (map_hitCell.count(cell) != 0) {
1225  double ee = energy + map_hitCell[cell].second;
1226  map_hitCell[cell] = std::pair<uint32_t, double>(detId, ee);
1227  } else {
1228  map_hitCell[cell] = std::pair<uint32_t, double>(detId, energy);
1229  }
1230 #ifdef EDM_ML_DEBUG
1231  edm::LogVerbatim("HGCSim") << "RecHit: " << layer << " " << global.x() << " " << global.y() << " " << global.z()
1232  << " " << energy;
1233 #endif
1234  }
1235 
1236  for (const auto& itr : map_hitLayer) {
1237  int layer = itr.first;
1238  double energy = itr.second;
1239  double zp = hgcons_[type]->waferZ(layer, true);
1240 #ifdef EDM_ML_DEBUG
1241  edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer << " " << zp << " " << energy;
1242 #endif
1243  hRecHitLng_[type]->Fill(zp, energy);
1244  hRecHitLng1_[type]->Fill(layer, energy);
1245  }
1246 
1247  for (const auto& itr : map_hitCell) {
1248  DetId detId = ((itr.second).first);
1249  double energy = ((itr.second).second);
1250  GlobalPoint global = hgeom_[type]->getPosition(detId);
1251  hRecHitLat_[type]->Fill(global.x(), global.y(), energy);
1252  }
1253 }
1254 
1256  for (const auto& v : *hgcPH) {
1257  double energy = v.energy();
1258  std::string name = v.vname();
1259  unsigned int id = v.id();
1260 #ifdef EDM_ML_DEBUG
1261  double time = v.time();
1262  edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::analyzePassiveHits:Energy:"
1263  << "Time:Name:Id : " << energy << ":" << time << ":" << name << ":" << id;
1264 #endif
1265 
1266  if (subdet == 1) {
1267  hgcPassiveEEEnergy_.push_back(energy);
1268  hgcPassiveEEName_.push_back(name);
1269  hgcPassiveEEID_.push_back(id);
1270  } else if (subdet == 2) {
1271  hgcPassiveFHEnergy_.push_back(energy);
1272  hgcPassiveFHName_.push_back(name);
1273  hgcPassiveFHID_.push_back(id);
1274  } else if (subdet == 3) {
1275  hgcPassiveBHEnergy_.push_back(energy);
1276  hgcPassiveBHName_.push_back(name);
1277  hgcPassiveBHID_.push_back(id);
1278  } else if (subdet == 4) {
1279  hgcPassiveCMSEEnergy_.push_back(energy);
1280  hgcPassiveCMSEName_.push_back(name);
1281  hgcPassiveCMSEID_.push_back(id);
1282  } else if (subdet == 5) {
1283  hgcPassiveBeamEnergy_.push_back(energy);
1284  hgcPassiveBeamName_.push_back(name);
1285  hgcPassiveBeamID_.push_back(id);
1286  }
1287  }
1288 }
1289 
1290 bool HGCalTBAnalyzer::sortTime(const std::pair<double, double>& i, const std::pair<double, double>& j) {
1291  return i.first < j.first;
1292 }
1293 
1294 // define this as a plug-in
std::vector< TH1D * > hSimHitLayEn2EE_
std::vector< float > simHitLayEn2EE_
Log< level::Info, true > LogVerbatim
const bool doTree_
std::vector< float > hgcPassiveBeamEnergy_
const edm::InputTag labelGen_
std::vector< PCaloHit > PCaloHitContainer
const double stoc_smear_time_200_
std::vector< float > simHitLayEn2BH_
std::vector< uint32_t > simHitCellIdFH_
void analyzeSimTracks(edm::Handle< edm::SimTrackContainer > const &SimTk, edm::Handle< edm::SimVertexContainer > const &SimVtx)
std::vector< float > pxBeamMC_
const std::string detectorBH_
int32_t *__restrict__ iv
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
int icolAbs() const
Definition: AHCalDetId.h:38
uint16_t *__restrict__ id
double cellThickness(int layer, int waferU, int waferV) const
const std::string detectorBeam_
TProfile * hSimHitLng1_[3]
TProfile * hSimHitLng2_[3]
TProfile * hRecHitLng_[3]
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
const edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHFH_
std::vector< float > simHitCellTimeFirstHitBH_
std::vector< std::string > hgcPassiveFHName_
std::vector< float > simHitLayEn1EE_
std::vector< float > hgcPassiveBHEnergy_
const edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsBH_
uint32_t ID
Definition: Definitions.h:24
const edm::ESGetToken< HGCalDDDConstants, IdealGeometryRecord > tokDDDFH_
std::vector< float > zBeamMC_
static uint32_t packHexagonIndex(int subdet, int z, int lay, int wafer, int celltyp, int cell)
std::vector< float > pBeamMC_
std::vector< float > xBeamMC_
T y() const
Definition: PV3DBase.h:60
std::vector< int > hgcPassiveBeamID_
edm::Service< TFileService > fs_
HGCalTBAnalyzer(edm::ParameterSet const &)
void beginRun(edm::Run const &, edm::EventSetup const &) override
static bool sortTime(const std::pair< double, double > &i, const std::pair< double, double > &j)
std::vector< float > pyBeamMC_
std::vector< float > simHitCellTimeLastHitFH_
int irow() const
get the row number
Definition: AHCalDetId.cc:29
const edm::EDGetTokenT< HGCRecHitCollection > tok_hitrFH_
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
const edm::InputTag labelPassiveEE_
const edm::InputTag labelDigiEE_
const bool doRecHits_
int zside(DetId const &)
std::vector< std::string > hgcPassiveBHName_
const double zFrontFH_
const edm::ESGetToken< HGCalDDDConstants, IdealGeometryRecord > tokDDDEE_
const edm::InputTag labelRHitBH_
std::pair< float, float > locateCell(int cell, int lay, int type, bool reco) const
std::vector< float > hgcPassiveCMSEEnergy_
const double gev2mip300_
const edm::InputTag labelRHitFH_
std::vector< int > hgcPassiveEEID_
std::vector< float > hgcPassiveEEEnergy_
wrapper for a data word
Definition: HGCSample.h:13
const edm::EDGetTokenT< edm::SimVertexContainer > tok_simVtx_
int depth() const
get the layer number
Definition: AHCalDetId.cc:43
const edm::ESGetToken< HGCalGeometry, IdealGeometryRecord > tokGeomEE_
std::vector< float > simHitCellEnFH_
constexpr std::array< uint8_t, layerIndexSize > layer
TProfile2D * hSimHitLat_[3]
const double zFrontBH_
bool getData(T &iHolder) const
Definition: EventSetup.h:122
U second(std::pair< T, U > const &p)
const std::string labelHitEE_
std::vector< TH1D * > hSimHitLayEn1EE_
const bool doPassiveEE_
void analyzeSimHits(int type, std::vector< PCaloHit > &hits, double zFront)
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Definition: Event.h:563
const bool doBeam_
int icol() const
get the column number
Definition: AHCalDetId.cc:36
std::pair< int, int > simToReco(int cell, int layer, int mod, bool half) const
int iEvent
Definition: GenABIO.cc:224
std::vector< TH1D * > hSimHitLayEn2BH_
unsigned int layers(bool reco) const
std::vector< float > yBeamMC_
const edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHCMSE_
std::vector< int > pdgIdBeamMC_
std::unique_ptr< AHCalGeometry > ahcalGeom_
std::vector< int > idBeams_
std::vector< int > simHitCellColBH_
const HGCalGeometry * hgeom_[2]
const edm::ESGetToken< HGCalGeometry, IdealGeometryRecord > tokGeomFH_
const edm::EDGetTokenT< HGCalDigiCollection > tok_digiFH_
std::vector< float > simHitCellTime15MipBH_
std::vector< std::string > hgcPassiveCMSEName_
const edm::EDGetTokenT< edm::SimTrackContainer > tok_simTk_
const double zFrontEE_
~HGCalTBAnalyzer() override
T z() const
Definition: PV3DBase.h:61
const edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsFH_
const bool ifBeam_
const double stoc_smear_time_300_
void analyze(edm::Event const &, edm::EventSetup const &) override
std::vector< int > simHitCellLayerBH_
std::vector< TH1D * > hSimHitLayEnBeam_
std::vector< int > simHitCellRowBH_
Transition
Definition: Transition.h:12
const bool doPassive_
std::vector< float > simHitCellTimeFirstHitEE_
const HGCalDDDConstants * hgcons_[2]
const edm::EDGetTokenT< HGCRecHitCollection > tok_hitrEE_
int wafer() const
get the wafer #
Definition: HGCalDetId.h:40
int waferTypeL(int wafer) const
std::vector< float > simHitCellTimeFirstHitFH_
const edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsEE_
ParameterDescriptionBase * add(U const &iLabel, T const &value)
void analyzeRecHits(int type, const edm::Handle< HGCRecHitCollection > &hits)
bool isValid() const
Definition: HandleBase.h:70
std::vector< TH1D * > hSimHitLayEn1FH_
std::vector< uint32_t > simHitCellIdBH_
const int sampleIndex_
std::vector< PassiveHit > PassiveHitContainer
Definition: PassiveHit.h:100
void analyzePassiveHits(edm::Handle< edm::PassiveHitContainer > const &hgcPh, int subdet)
double waferZ(int layer, bool reco) const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
#define M_PI
const edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHBeam_
std::vector< float > simHitCellTime15MipEE_
const bool doPassiveBH_
std::vector< std::string > hgcPassiveBeamName_
Basic2DVector< T > xy() const
Definition: DetId.h:17
std::vector< float > hgcPassiveFHEnergy_
std::vector< float > simHitCellTimeLastHitBH_
std::vector< float > simHitCellTimeLastHitEE_
const edm::InputTag labelRHitEE_
std::vector< int > hgcPassiveCMSEID_
const std::string labelHitBH_
TProfile2D * hRecHitLat_[3]
const bool doSimHits_
#define debug
Definition: HDRShower.cc:19
const bool doDigis_
int cell() const
get the absolute value of the cell #&#39;s in x and y
Definition: HGCalDetId.h:37
std::vector< float > simHitCellEnBeam_
const edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHBH_
TProfile * hSimHitLng_[3]
const bool doTreeCell_
const edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsBeam_
std::vector< SimVertex > SimVertexContainer
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
const edm::InputTag labelDigiFH_
std::vector< float > pzBeamMC_
void add(std::string const &label, ParameterSetDescription const &psetDescription)
std::vector< std::string > hgcPassiveEEName_
const edm::InputTag labelPassiveBH_
void endRun(edm::Run const &, edm::EventSetup const &) override
std::vector< float > simHitCellTime15MipFH_
void analyzeDigi(int type, const T1 &detId, uint16_t adc)
const edm::InputTag labelPassiveCMSE_
int irowAbs() const
Definition: AHCalDetId.h:35
std::vector< int > hgcPassiveBHID_
const bool doPassiveHE_
GlobalPoint getPosition(const DetId &id, bool debug=false) const
const edm::EDGetTokenT< HGCalDigiCollection > tok_digiBH_
string end
Definition: dataset.py:937
std::vector< float > simHitLayEn1FH_
std::vector< int > hgcPassiveFHID_
std::vector< float > simHitCellEnBH_
static void unpackIndex(const uint32_t &idx, int &det, int &lay, int &x, int &y)
const edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHEE_
const edm::EDGetTokenT< HGCRecHitCollection > tok_hitrBH_
Log< level::Warning, false > LogWarning
int col
Definition: cuy.py:1009
std::vector< float > simHitLayEn1BH_
std::vector< uint32_t > simHitCellIdEE_
std::vector< TH1D * > hSimHitLayEn1BH_
const double gev2mip200_
const edm::InputTag labelPassiveBeam_
int zside() const
get the z-side of the cell (1/-1)
Definition: AHCalDetId.h:32
std::vector< float > simHitCellEnEE_
T x() const
Definition: PV3DBase.h:59
const edm::InputTag labelDigiBH_
uint16_t data() const
Definition: HGCSample.h:70
const edm::InputTag labelPassiveFH_
std::vector< uint32_t > simHitCellIdBeam_
std::vector< SimTrack > SimTrackContainer
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
std::vector< float > simHitLayEn2FH_
tuple size
Write out results.
const std::string detectorEE_
static void unpackHexagonIndex(const uint32_t &idx, int &subdet, int &z, int &lay, int &wafer, int &celltyp, int &cell)
const std::string labelHitFH_
int layer() const
get the layer #
Definition: HGCalDetId.h:46
TProfile * hRecHitLng1_[3]
Definition: Run.h:45
const edm::EDGetTokenT< HGCalDigiCollection > tok_digiEE_
std::vector< TH1D * > hSimHitLayEn2FH_
uint16_t *__restrict__ uint16_t const *__restrict__ adc
void beginJob() override
const edm::EDGetTokenT< edm::HepMCProduct > tok_hepMC_
const std::string labelHitBeam_
const std::string detectorFH_
std::vector< float > simHitLayEnBeam_