CMS 3D CMS Logo

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