CMS 3D CMS Logo
HGCalTBAnalyzer.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 #include <iostream>
4 #include <fstream>
5 #include <vector>
6 #include <map>
7 #include <string>
8 
9 // user include files
10 #include "CommonTools/UtilAlgos/interface/TFileService.h"
11 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
12 #include "DataFormats/ForwardDetId/interface/HGCEEDetId.h"
13 #include "DataFormats/ForwardDetId/interface/HGCHEDetId.h"
14 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
15 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
16 
17 #include "FWCore/Framework/interface/Frameworkfwd.h"
18 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
19 #include "FWCore/Framework/interface/ESHandle.h"
20 #include "FWCore/Framework/interface/Event.h"
21 #include "FWCore/Framework/interface/EventSetup.h"
22 #include "FWCore/Framework/interface/MakerMacros.h"
23 #include "FWCore/ParameterSet/interface/ParameterSet.h"
24 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
25 #include "FWCore/ServiceRegistry/interface/Service.h"
26 #include "FWCore/Utilities/interface/InputTag.h"
27 
28 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
29 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
30 #include "Geometry/Records/interface/CaloGeometryRecord.h"
31 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
32 #include "Geometry/Records/interface/IdealGeometryRecord.h"
33 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
34 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
35 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
36 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
37 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
38 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
39 #include "SimG4CMS/HGCalTestBeam/interface/AHCalDetId.h"
40 
41 #include "SimDataFormats/CaloHit/interface/PassiveHit.h"
42 
43 // Root objects
44 #include "TROOT.h"
45 #include "TSystem.h"
46 #include "TFile.h"
47 #include "TH1.h"
48 #include "TH2.h"
49 #include "TProfile.h"
50 #include "TProfile2D.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 
57 public:
58  explicit HGCalTBAnalyzer(edm::ParameterSet const&);
59  ~HGCalTBAnalyzer();
60 
61  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
62 
63 private:
64  virtual void beginJob() override ;
65  virtual void beginRun(edm::Run const&, edm::EventSetup const&) override;
66  virtual void endRun(edm::Run const&, edm::EventSetup const&) override {}
67  virtual void analyze(edm::Event const&, edm::EventSetup const&) override;
68  void analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront);
69  void analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
70  edm::Handle<edm::SimVertexContainer> const& SimVtx);
71  template<class T1> void analyzeDigi(int type, const T1& detId, uint16_t adc);
72  void analyzeRecHits(int type, edm::Handle<HGCRecHitCollection> & hits);
73  void analyzePassiveHits (edm::Handle<edm::PassiveHitContainer> const& hgcPh,
74  int subdet);
75 
76  edm::Service<TFileService> fs_;
77  const HGCalDDDConstants *hgcons_[2];
78  const HGCalGeometry *hgeom_[2];
79  bool ifEE_, ifFH_, ifBH_, ifBeam_;
80  bool doTree_, doTreeCell_;
81  bool doSimHits_, doDigis_, doRecHits_;
82  bool doPassive_;
83  std::string detectorEE_, detectorFH_;
84  std::string detectorBH_, detectorBeam_;
85  double zFrontEE_, zFrontFH_, zFrontBH_;
86  int sampleIndex_;
87  std::vector<int> idBeams_;
88  edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsEE_, tok_hitsFH_;
89  edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsBH_, tok_hitsBeam_;
90  edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
91  edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
92  edm::EDGetToken tok_digiEE_, tok_digiFH_, tok_digiBH_;
93  edm::EDGetToken tok_hitrEE_, tok_hitrFH_, tok_hitrBH_;
94  edm::EDGetTokenT<edm::HepMCProduct> tok_hepMC_;
95  edm::EDGetTokenT<edm::PassiveHitContainer> tok_hgcPHEE_, tok_hgcPHFH_;
96  edm::EDGetTokenT<edm::PassiveHitContainer> tok_hgcPHBH_, tok_hgcPHCMSE_;
97 
98  TTree *tree_;
99  TH1D *hSimHitE_[4], *hSimHitT_[4];
100  TH1D *hDigiADC_[3], *hDigiLng_[2];
101  TH1D *hRecHitE_[3], *hSimHitEn_[4], *hBeam_;
102  TH2D *hDigiOcc_[3], *hRecHitOcc_[3];
103  TProfile *hSimHitLng_[3], *hSimHitLng1_[3];
104  TProfile *hSimHitLng2_[3];
105  TProfile *hRecHitLng_[3], *hRecHitLng1_[3];
106  TProfile2D *hSimHitLat_[3], *hRecHitLat_[3];
107  std::vector<TH1D*> hSimHitLayEn1EE_, hSimHitLayEn2EE_;
108  std::vector<TH1D*> hSimHitLayEn1FH_, hSimHitLayEn2FH_;
109  std::vector<TH1D*> hSimHitLayEn1BH_, hSimHitLayEn2BH_;
110  std::vector<TH1D*> hSimHitLayEnBeam_;
111  std::vector<float> simHitLayEn1EE_, simHitLayEn2EE_;
112  std::vector<float> simHitLayEn1FH_, simHitLayEn2FH_;
113  std::vector<float> simHitLayEn1BH_, simHitLayEn2BH_;
114  std::vector<float> simHitLayEnBeam_;
115  std::vector<uint32_t> simHitCellIdEE_, simHitCellIdFH_;
116  std::vector<uint32_t> simHitCellIdBH_, simHitCellIdBeam_;
117  std::vector<float> simHitCellEnEE_, simHitCellEnFH_;
118  std::vector<float> simHitCellEnBH_, simHitCellEnBeam_;
119 
120  std::vector<float> hgcPassiveEEEnergy_, hgcPassiveFHEnergy_, hgcPassiveBHEnergy_, hgcPassiveCMSEEnergy_;
121  std::vector<std::string> hgcPassiveEEName_, hgcPassiveFHName_, hgcPassiveBHName_, hgcPassiveCMSEName_;
122  std::vector<int> hgcPassiveEEID_, hgcPassiveFHID_, hgcPassiveBHID_, hgcPassiveCMSEID_;
123 
124  double xBeam_, yBeam_, zBeam_, pBeam_;
125 };
126 
127 HGCalTBAnalyzer::HGCalTBAnalyzer(const edm::ParameterSet& iConfig) {
128 
129  usesResource("TFileService");
130 
131  //now do whatever initialization is needed
132  detectorEE_ = iConfig.getParameter<std::string>("DetectorEE");
133  detectorFH_ = iConfig.getParameter<std::string>("DetectorFH");
134  detectorBH_ = iConfig.getParameter<std::string>("DetectorBH");
135  detectorBeam_= iConfig.getParameter<std::string>("DetectorBeam");
136  ifEE_ = iConfig.getParameter<bool>("UseEE");
137  ifFH_ = iConfig.getParameter<bool>("UseFH");
138  ifBH_ = iConfig.getParameter<bool>("UseBH");
139  ifBeam_ = iConfig.getParameter<bool>("UseBeam");
140  zFrontEE_ = iConfig.getParameter<double>("ZFrontEE");
141  zFrontFH_ = iConfig.getParameter<double>("ZFrontFH");
142  zFrontBH_ = iConfig.getParameter<double>("ZFrontBH");
143  idBeams_ = iConfig.getParameter<std::vector<int>>("IDBeams");
144  doSimHits_ = iConfig.getParameter<bool>("DoSimHits");
145  doDigis_ = iConfig.getParameter<bool>("DoDigis");
146  sampleIndex_ = iConfig.getParameter<int>("SampleIndex");
147  doRecHits_ = iConfig.getParameter<bool>("DoRecHits");
148  doTree_ = iConfig.getUntrackedParameter<bool>("DoTree",false);
149  doTreeCell_ = iConfig.getUntrackedParameter<bool>("DoTreeCell",false);
150  doPassive_ = iConfig.getUntrackedParameter<bool>("DoPassive",false);
151 
152 #ifdef EDM_ML_DEBUG
153  std::cout << "HGCalTBAnalyzer:: SimHits = " << doSimHits_ << " Digis = "
154  << doDigis_ << ":" << sampleIndex_ << " RecHits = " << doRecHits_
155  << " useDets " << ifEE_ << ":" << ifFH_ << ":" << ifBH_ << ":"
156  << ifBeam_ << " zFront " << zFrontEE_ << ":" << zFrontFH_ << ":"
157  << zFrontBH_ << " IdBeam " << idBeams_.size() << ":";
158  for (auto id : idBeams_) std::cout << " " << id;
159  std::cout << std::endl;
160 #endif
161  if (idBeams_.size() == 0) idBeams_.push_back(1001);
162 
163  edm::InputTag tmp0 = iConfig.getParameter<edm::InputTag>("GeneratorSrc");
164  tok_hepMC_ = consumes<edm::HepMCProduct>(tmp0);
165 
166 #ifdef EDM_ML_DEBUG
167  std::cout << "HGCalTBAnalyzer:: GeneratorSource = " << tmp0 << std::endl;
168 #endif
169  std::string tmp1 = iConfig.getParameter<std::string>("CaloHitSrcEE");
170  tok_hitsEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits",tmp1));
171  tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
172  tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
173  edm::InputTag tmp2 = iConfig.getParameter<edm::InputTag>("DigiSrcEE");
174  tok_digiEE_ = consumes<HGCEEDigiCollection>(tmp2);
175  edm::InputTag tmp3 = iConfig.getParameter<edm::InputTag>("RecHitSrcEE");
176  tok_hitrEE_ = consumes<HGCRecHitCollection>(tmp3);
177 #ifdef EDM_ML_DEBUG
178  if (ifEE_) {
179  std::cout << "HGCalTBAnalyzer:: Detector " << detectorEE_ << " with tags "
180  << tmp1 << ", " << tmp2 << ", " << tmp3 << std::endl;
181  }
182 #endif
183  tmp1 = iConfig.getParameter<std::string>("CaloHitSrcFH");
184  tok_hitsFH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits",tmp1));
185  tmp2 = iConfig.getParameter<edm::InputTag>("DigiSrcFH");
186  tok_digiFH_ = consumes<HGCHEDigiCollection>(tmp2);
187  tmp3 = iConfig.getParameter<edm::InputTag>("RecHitSrcFH");
188  tok_hitrFH_ = consumes<HGCRecHitCollection>(tmp3);
189 #ifdef EDM_ML_DEBUG
190  if (ifFH_) {
191  std::cout << "HGCalTBAnalyzer:: Detector " << detectorFH_ << " with tags "
192  << tmp1 << ", " << tmp2 << ", " << tmp3 << std::endl;
193  }
194 #endif
195  tmp1 = iConfig.getParameter<std::string>("CaloHitSrcBH");
196  tok_hitsBH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits",tmp1));
197  tmp2 = iConfig.getParameter<edm::InputTag>("DigiSrcBH");
198  tok_digiBH_ = consumes<HGCBHDigiCollection>(tmp2);
199  tmp3 = iConfig.getParameter<edm::InputTag>("RecHitSrcBH");
200  tok_hitrBH_ = consumes<HGCRecHitCollection>(tmp3);
201 
203  edm::InputTag tmp = iConfig.getParameter<edm::InputTag>("HGCPassiveEE");
204  tok_hgcPHEE_ = consumes<edm::PassiveHitContainer>(tmp);
205 
206  tmp = iConfig.getParameter<edm::InputTag>("HGCPassiveFH");
207  tok_hgcPHFH_ = consumes<edm::PassiveHitContainer>(tmp);
208 
209  tmp = iConfig.getParameter<edm::InputTag>("HGCPassiveBH");
210  tok_hgcPHBH_ = consumes<edm::PassiveHitContainer>(tmp);
211 
212  tmp = iConfig.getParameter<edm::InputTag>("HGCPassiveCMSE");
213  tok_hgcPHCMSE_ = consumes<edm::PassiveHitContainer>(tmp);
214 
215 #ifdef EDM_ML_DEBUG
216  if (ifBH_) {
217  std::cout << "HGCalTBAnalyzer:: Detector " << detectorBH_ << " with tags "
218  << tmp1 << ", " << tmp2 << ", " << tmp3 << std::endl;
219  }
220 #endif
221  tmp1 = iConfig.getParameter<std::string>("CaloHitSrcBeam");
222  tok_hitsBeam_= consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits",tmp1));
223 #ifdef EDM_ML_DEBUG
224  if (ifBeam_) {
225  std::cout << "HGCalTBAnalyzer:: Detector " << detectorBeam_
226  << " with tags " << tmp1 << std::endl;
227  }
228 #endif
229 }
230 
231 HGCalTBAnalyzer::~HGCalTBAnalyzer() {}
232 
233 void HGCalTBAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
234  edm::ParameterSetDescription desc;
235  desc.setUnknown();
236  desc.add<std::string>("DetectorEE","HGCalEESensitive");
237  desc.add<bool>("UseEE",true);
238  desc.add<double>("ZFrontEE",0.0);
239  desc.add<std::string>("CaloHitSrcEE","HGCHitsEE");
240  desc.add<edm::InputTag>("DigiSrcEE",edm::InputTag("mix","HGCDigisEE"));
241  desc.add<edm::InputTag>("RecHitSrcEE",edm::InputTag("HGCalRecHit","HGCEERecHits"));
242  desc.add<std::string>("DetectorFH","HGCalHESiliconSensitive");
243  desc.add<bool>("UseFH",false);
244  desc.add<double>("ZFrontFH",0.0);
245  desc.add<std::string>("CaloHitSrcFH","HGCHitsHEfront");
246  desc.add<edm::InputTag>("DigiSrcFH",edm::InputTag("mix","HGCDigisHEfront"));
247  desc.add<edm::InputTag>("RecHitSrcFH",edm::InputTag("HGCalRecHit","HGCHEFRecHits"));
248  desc.add<std::string>("DetectorBH","AHCal");
249  desc.add<bool>("UseBH",false);
250  desc.add<double>("ZFrontBH",0.0);
251  desc.add<std::string>("CaloHitSrcBH","HcalHits");
252  desc.add<edm::InputTag>("DigiSrcBH",edm::InputTag("mix","HGCDigisHEback"));
253  desc.add<edm::InputTag>("RecHitSrcBH",edm::InputTag("HGCalRecHit","HGCHEBRecHits"));
254  desc.add<std::string>("DetectorBeam","HcalTB06BeamDetector");
255  desc.add<bool>("UseBeam",false);
256  desc.add<std::string>("CaloHitSrcBeam","HcalTB06BeamHits");
257  std::vector<int> ids = {1000,1001,1002,1003,1004,1005,1006,1007,1008,1011,1012,1013,1014,2001,2002,2003,2004,2005};
258  desc.add<std::vector<int>>("IDBeams",ids);
259  desc.add<edm::InputTag>("GeneratorSrc",edm::InputTag("generatorSmeared"));
260  desc.add<edm::InputTag>("HGCPassiveEE",edm::InputTag("g4SimHits","HGCalEEPassiveHits"));
261  desc.add<edm::InputTag>("HGCPassiveFH",edm::InputTag("g4SimHits","HGCalHEPassiveHits"));
262  desc.add<edm::InputTag>("HGCPassiveBH",edm::InputTag("g4SimHits","HGCalAHPassiveHits"));
263  desc.add<edm::InputTag>("HGCPassiveCMSE",edm::InputTag("g4SimHits","CMSEPassiveHits"));
264 
265  desc.add<bool>("DoSimHits",true);
266  desc.add<bool>("DoDigis",true);
267  desc.add<bool>("DoRecHits",true);
268  desc.add<int>("SampleIndex",0);
269  desc.addUntracked<bool>("DoTree",true);
270  desc.addUntracked<bool>("DoTreeCell",true);
271  desc.addUntracked<bool>("DoPassive",false);
272 
273  descriptions.add("HGCalTBAnalyzer",desc);
274 }
275 
276 void HGCalTBAnalyzer::beginJob() {
277 
278  char name[40], title[100];
279  hBeam_ = fs_->make<TH1D>("BeamP", "Beam Momentum", 1000, 0, 1000.0);
280  for (int i=0; i<3; ++i) {
281  bool book(ifEE_);
282  std::string det(detectorEE_);
283  if (i == 1) {
284  book = ifFH_;
285  det = detectorFH_;
286  } else if (i == 2) {
287  book = ifBH_;
288  det = detectorBH_;
289  }
290 
291  if (doSimHits_ && book) {
292  sprintf (name, "SimHitEn%s", det.c_str());
293  sprintf (title,"Sim Hit Energy for %s", det.c_str());
294  hSimHitE_[i] = fs_->make<TH1D>(name,title,100000,0.,0.2);
295  sprintf (name, "SimHitEnX%s", det.c_str());
296  sprintf (title,"Sim Hit Energy for %s", det.c_str());
297  hSimHitEn_[i] = fs_->make<TH1D>(name,title,100000,0.,0.2);
298  sprintf (name, "SimHitTm%s", det.c_str());
299  sprintf (title,"Sim Hit Timing for %s", det.c_str());
300  hSimHitT_[i] = fs_->make<TH1D>(name,title,5000,0.,500.0);
301  sprintf (name, "SimHitLat%s", det.c_str());
302  sprintf (title,"Lateral Shower profile (Sim Hit) for %s", det.c_str());
303  hSimHitLat_[i] = fs_->make<TProfile2D>(name,title,100,-100.,100.,100,-100.,100.);
304  sprintf (name, "SimHitLng%s", det.c_str());
305  sprintf (title,"Longitudinal Shower profile (Sim Hit) for %s",det.c_str());
306  hSimHitLng_[i] = fs_->make<TProfile>(name,title,50,0.,100.);
307  sprintf (name, "SimHitLng1%s", det.c_str());
308  sprintf (title,"Longitudinal Shower profile (Layer) for %s",det.c_str());
309  hSimHitLng1_[i] = fs_->make<TProfile>(name,title,200,0.,100.);
310  sprintf (name, "SimHitLng2%s", det.c_str());
311  sprintf (title,"Longitudinal Shower profile (Layer) for %s",det.c_str());
312  hSimHitLng2_[i] = fs_->make<TProfile>(name,title,200,0.,100.);
313  }
314 
315  if (doDigis_ && book) {
316  sprintf (name, "DigiADC%s", det.c_str());
317  sprintf (title,"ADC at Digi level for %s", det.c_str());
318  hDigiADC_[i] = fs_->make<TH1D>(name,title,100,0.,100.0);
319  sprintf (name, "DigiOcc%s", det.c_str());
320  sprintf (title,"Occupancy (Digi)for %s", det.c_str());
321  hDigiOcc_[i] = fs_->make<TH2D>(name,title,100,-10.,10.,100,-10.,10.);
322  sprintf (name, "DigiLng%s", det.c_str());
323  sprintf (title,"Longitudinal Shower profile (Digi) for %s",det.c_str());
324  hDigiLng_[i] = fs_->make<TH1D>(name,title,100,0.,10.);
325  }
326 
327  if (doRecHits_ && book) {
328  sprintf (name, "RecHitEn%s", det.c_str());
329  sprintf (title,"Rec Hit Energy for %s", det.c_str());
330  hRecHitE_[i] = fs_->make<TH1D>(name,title,1000,0.,10.0);
331  sprintf (name, "RecHitOcc%s", det.c_str());
332  sprintf (title,"Occupancy (Rec Hit)for %s", det.c_str());
333  hRecHitOcc_[i] = fs_->make<TH2D>(name,title,100,-10.,10.,100,-10.,10.);
334  sprintf (name, "RecHitLat%s", det.c_str());
335  sprintf (title,"Lateral Shower profile (Rec Hit) for %s", det.c_str());
336  hRecHitLat_[i] = fs_->make<TProfile2D>(name,title,100,-10.,10.,100,-10.,10.);
337  sprintf (name, "RecHitLng%s", det.c_str());
338  sprintf (title,"Longitudinal Shower profile (Rec Hit) for %s",det.c_str());
339  hRecHitLng_[i] = fs_->make<TProfile>(name,title,100,0.,10.);
340  sprintf (name, "RecHitLng1%s", det.c_str());
341  sprintf (title,"Longitudinal Shower profile vs Layer for %s",det.c_str());
342  hRecHitLng1_[i] = fs_->make<TProfile>(name,title,120,0.,60.);
343  }
344  }
345  if (ifBeam_ && doSimHits_) {
346  sprintf (name, "SimHitEn%s", detectorBeam_.c_str());
347  sprintf (title,"Sim Hit Energy for %s", detectorBeam_.c_str());
348  hSimHitE_[3] = fs_->make<TH1D>(name,title,100000,0.,0.2);
349  sprintf (name, "SimHitEnX%s", detectorBeam_.c_str());
350  sprintf (title,"Sim Hit Energy for %s", detectorBeam_.c_str());
351  hSimHitEn_[3] = fs_->make<TH1D>(name,title,100000,0.,0.2);
352  sprintf (name, "SimHitTm%s", detectorBeam_.c_str());
353  sprintf (title,"Sim Hit Timing for %s", detectorBeam_.c_str());
354  hSimHitT_[3] = fs_->make<TH1D>(name,title,5000,0.,500.0);
355  }
356  if (doSimHits_ && doTree_) {
357  tree_ = fs_->make<TTree>("HGCTB","SimHitEnergy");
358  tree_->Branch("simHitLayEn1EE", &simHitLayEn1EE_);
359  tree_->Branch("simHitLayEn2EE", &simHitLayEn2EE_);
360  tree_->Branch("simHitLayEn1FH", &simHitLayEn1FH_);
361  tree_->Branch("simHitLayEn2FH", &simHitLayEn2FH_);
362  tree_->Branch("simHitLayEn1BH", &simHitLayEn1BH_);
363  tree_->Branch("simHitLayEn2BH", &simHitLayEn2BH_);
364  tree_->Branch("xBeam", &xBeam_, "xBeam/D");
365  tree_->Branch("yBeam", &yBeam_, "yBeam/D");
366  tree_->Branch("zBeam", &zBeam_, "zBeam/D");
367  tree_->Branch("pBeam", &pBeam_, "pBeam/D");
368  if (doTreeCell_) {
369  tree_->Branch("simHitCellIdEE", &simHitCellIdEE_);
370  tree_->Branch("simHitCellEnEE", &simHitCellEnEE_);
371  tree_->Branch("simHitCellIdFH", &simHitCellIdFH_);
372  tree_->Branch("simHitCellEnFH", &simHitCellEnFH_);
373  tree_->Branch("simHitCellIdBH", &simHitCellIdBH_);
374  tree_->Branch("simHitCellEnBH", &simHitCellEnBH_);
375  tree_->Branch("simHitCellIdBeam", &simHitCellIdBeam_);
376  tree_->Branch("simHitCellEnBeam", &simHitCellEnBeam_);
377  }
378  }
379 
380  if (doPassive_ && doTree_) {
381  tree_->Branch("hgcPassiveEEEnergy", &hgcPassiveEEEnergy_);
382  tree_->Branch("hgcPassiveEEName", &hgcPassiveEEName_);
383  tree_->Branch("hgcPassiveEEID", &hgcPassiveEEID_);
384  tree_->Branch("hgcPassiveFHEnergy", &hgcPassiveFHEnergy_);
385  tree_->Branch("hgcPassiveFHName", &hgcPassiveFHName_);
386  tree_->Branch("hgcPassiveFHID", &hgcPassiveFHID_);
387  tree_->Branch("hgcPassiveBHEnergy", &hgcPassiveBHEnergy_);
388  tree_->Branch("hgcPassiveBHName", &hgcPassiveBHName_);
389  tree_->Branch("hgcPassiveBHID", &hgcPassiveBHID_);
390  tree_->Branch("hgcPassiveCMSEEnergy", &hgcPassiveCMSEEnergy_);
391  tree_->Branch("hgcPassiveCMSEName", &hgcPassiveCMSEName_);
392  tree_->Branch("hgcPassiveCMSEID", &hgcPassiveCMSEID_);
393  }
394 }
395 
396 void HGCalTBAnalyzer::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
397 
398  char name[40], title[100];
399  if (ifEE_) {
400  edm::ESHandle<HGCalDDDConstants> pHGDC;
401  iSetup.get<IdealGeometryRecord>().get(detectorEE_, pHGDC);
402  hgcons_[0] = &(*pHGDC);
403  if (doDigis_ || doRecHits_) {
404  edm::ESHandle<HGCalGeometry> geom;
405  iSetup.get<IdealGeometryRecord>().get(detectorEE_, geom);
406  hgeom_[0] = geom.product();
407  } else {
408  hgeom_[0] = 0;
409  }
410  for (unsigned int l=0; l<hgcons_[0]->layers(false); ++l) {
411  sprintf (name, "SimHitEnA%d%s", l, detectorEE_.c_str());
412  sprintf (title,"Sim Hit Energy in SIM layer %d for %s",l+1,
413  detectorEE_.c_str());
414  hSimHitLayEn1EE_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
415  if (l%3 == 0) {
416  sprintf (name, "SimHitEnB%d%s", (l/3+1), detectorEE_.c_str());
417  sprintf (title,"Sim Hit Energy in layer %d for %s",(l/3+1),
418  detectorEE_.c_str());
419  hSimHitLayEn2EE_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
420  }
421  }
422 #ifdef EDM_ML_DEBUG
423  std::cout << "HGCalTBAnalyzer::" << detectorEE_ << " defined with "
424  << hgcons_[0]->layers(false) << " layers" << std::endl;
425 #endif
426  } else {
427  hgcons_[0] = 0;
428  hgeom_[0] = 0;
429  }
430 
431  if (ifFH_) {
432  edm::ESHandle<HGCalDDDConstants> pHGDC;
433  iSetup.get<IdealGeometryRecord>().get(detectorFH_, pHGDC);
434  hgcons_[1] = &(*pHGDC);
435  if (doDigis_ || doRecHits_) {
436  edm::ESHandle<HGCalGeometry> geom;
437  iSetup.get<IdealGeometryRecord>().get(detectorFH_, geom);
438  hgeom_[1] = geom.product();
439  } else {
440  hgeom_[1] = 0;
441  }
442  for (unsigned int l=0; l<hgcons_[1]->layers(false); ++l) {
443  sprintf (name, "SimHitEnA%d%s", l, detectorFH_.c_str());
444  sprintf (title,"Sim Hit Energy in layer %d for %s",l+1,
445  detectorFH_.c_str());
446  hSimHitLayEn1FH_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
447  if (l%3 == 0) {
448  sprintf (name, "SimHitEnB%d%s", (l/3+1), detectorFH_.c_str());
449  sprintf (title,"Sim Hit Energy in layer %d for %s",(l/3+1),
450  detectorFH_.c_str());
451  hSimHitLayEn2FH_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
452  }
453  }
454 #ifdef EDM_ML_DEBUG
455  std::cout << "HGCalTBAnalyzer::" << detectorFH_ << " defined with "
456  << hgcons_[1]->layers(false) << " layers" << std::endl;
457 #endif
458  } else {
459  hgcons_[1] = 0;
460  hgeom_[1] = 0;
461  }
462 
463  if (ifBH_) {
464  for (int l=0; l<AHCalDetId::MaxDepth; ++l) {
465  sprintf (name, "SimHitEnA%d%s", l, detectorBH_.c_str());
466  sprintf (title,"Sim Hit Energy in layer %d for %s",l+1,
467  detectorBH_.c_str());
468  hSimHitLayEn1BH_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
469  sprintf (name, "SimHitEnB%d%s", l, detectorBH_.c_str());
470  sprintf (title,"Sim Hit Energy in layer %d for %s",l+1,
471  detectorBH_.c_str());
472  hSimHitLayEn2BH_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
473  }
474  }
475 
476  if (ifBeam_) {
477  for (unsigned int l=0; l<idBeams_.size(); ++l) {
478  sprintf (name, "SimHitEna%d%s", l, detectorBeam_.c_str());
479  sprintf (title, "Sim Hit Energy in type %d for %s",idBeams_[l],
480  detectorBeam_.c_str());
481  hSimHitLayEnBeam_.push_back(fs_->make<TH1D>(name,title,100000,0.,0.2));
482  }
483  }
484 }
485 
486 void HGCalTBAnalyzer::analyze(const edm::Event& iEvent,
487  const edm::EventSetup& iSetup) {
488 
489  //Generator input
490  edm::Handle<edm::HepMCProduct> evtMC;
491  iEvent.getByToken(tok_hepMC_,evtMC);
492  if (!evtMC.isValid()) {
493  edm::LogWarning("HGCal") << "no HepMCProduct found";
494  } else {
495  const HepMC::GenEvent * myGenEvent = evtMC->GetEvent();
496  unsigned int k(0);
497  for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
498  p != myGenEvent->particles_end(); ++p, ++k) {
499  if (k == 0) hBeam_->Fill((*p)->momentum().rho());
500 #ifdef EDM_ML_DEBUG
501  std::cout << "Particle[" << k << "] with p " << (*p)->momentum().rho()
502  << " theta " << (*p)->momentum().theta() << " phi "
503  << (*p)->momentum().phi() << std::endl;
504 #endif
505  }
506  }
507 
508  //Now the Simhits
509  if (doSimHits_) {
510  edm::Handle<edm::SimTrackContainer> SimTk;
511  iEvent.getByToken(tok_simTk_, SimTk);
512  edm::Handle<edm::SimVertexContainer> SimVtx;
513  iEvent.getByToken(tok_simVtx_, SimVtx);
514  analyzeSimTracks(SimTk, SimVtx);
515 
516  simHitLayEn1EE_.clear(); simHitLayEn2EE_.clear();
517  simHitLayEn1FH_.clear(); simHitLayEn2FH_.clear();
518  simHitLayEn1BH_.clear(); simHitLayEn2BH_.clear();
519  simHitLayEnBeam_.clear();
520  simHitCellIdEE_.clear(); simHitCellEnEE_.clear();
521  simHitCellIdFH_.clear(); simHitCellEnFH_.clear();
522  simHitCellIdBH_.clear(); simHitCellEnBH_.clear();
523  simHitCellIdBeam_.clear(); simHitCellEnBeam_.clear();
524  edm::Handle<edm::PCaloHitContainer> theCaloHitContainers;
525  std::vector<PCaloHit> caloHits;
526  if (ifEE_) {
527  simHitLayEn1EE_ = std::vector<float>(hgcons_[0]->layers(false),0);
528  simHitLayEn2EE_ = std::vector<float>(hgcons_[0]->layers(true),0);
529  iEvent.getByToken(tok_hitsEE_, theCaloHitContainers);
530  if (theCaloHitContainers.isValid()) {
531 #ifdef EDM_ML_DEBUG
532  std::cout << "PcalohitContainer for " << detectorEE_ << " has "
533  << theCaloHitContainers->size() << " hits" << std::endl;
534 #endif
535  caloHits.clear();
536  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(),
537  theCaloHitContainers->end());
538  analyzeSimHits(0, caloHits, zFrontEE_);
539  } else {
540 #ifdef EDM_ML_DEBUG
541  std::cout << "PCaloHitContainer does not exist for " << detectorEE_
542  << " !!!" << std::endl;
543 #endif
544  }
545  }
546  if (ifFH_) {
547  simHitLayEn1FH_ = std::vector<float>(hgcons_[1]->layers(false),0);
548  simHitLayEn2FH_ = std::vector<float>(hgcons_[1]->layers(true),0);
549  iEvent.getByToken(tok_hitsFH_, theCaloHitContainers);
550  if (theCaloHitContainers.isValid()) {
551 #ifdef EDM_ML_DEBUG
552  std::cout << "PcalohitContainer for " << detectorFH_ << " has "
553  << theCaloHitContainers->size() << " hits" << std::endl;
554 #endif
555  caloHits.clear();
556  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(),
557  theCaloHitContainers->end());
558  analyzeSimHits(1, caloHits, zFrontFH_);
559  } else {
560 #ifdef EDM_ML_DEBUG
561  std::cout << "PCaloHitContainer does not exist for " << detectorFH_
562  << " !!!" << std::endl;
563 #endif
564  }
565  }
566  if (ifBH_) {
567  simHitLayEn1BH_ = std::vector<float>(AHCalDetId::MaxDepth,0);
568  simHitLayEn2BH_ = std::vector<float>(AHCalDetId::MaxDepth,0);
569  iEvent.getByToken(tok_hitsBH_, theCaloHitContainers);
570  if (theCaloHitContainers.isValid()) {
571 #ifdef EDM_ML_DEBUG
572  std::cout << "PcalohitContainer for " << detectorBH_ << " has "
573  << theCaloHitContainers->size() << " hits" << std::endl;
574 #endif
575  caloHits.clear();
576  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(),
577  theCaloHitContainers->end());
578  analyzeSimHits(2, caloHits, zFrontBH_);
579  } else {
580 #ifdef EDM_ML_DEBUG
581  std::cout << "PCaloHitContainer does not exist for " << detectorBH_
582  << " !!!" << std::endl;
583 #endif
584  }
585  }
586  if (ifBeam_) {
587  simHitLayEnBeam_ = std::vector<float>(idBeams_.size(),0);
588  iEvent.getByToken(tok_hitsBeam_, theCaloHitContainers);
589  if (theCaloHitContainers.isValid()) {
590 #ifdef EDM_ML_DEBUG
591  std::cout << "PcalohitContainer for " << detectorBeam_ << " has "
592  << theCaloHitContainers->size() << " hits" << std::endl;
593 #endif
594  caloHits.clear();
595  caloHits.insert(caloHits.end(), theCaloHitContainers->begin(),
596  theCaloHitContainers->end());
597  analyzeSimHits(3, caloHits, 0.0);
598  } else {
599 #ifdef EDM_ML_DEBUG
600  std::cout << "PCaloHitContainer does not exist for " << detectorBeam_
601  << " !!!" << std::endl;
602 #endif
603  }
604  }
605  //if (doTree_) tree_->Fill();
606  }//if (doSimHits_)
607 
609  if (doPassive_) {
610  hgcPassiveEEEnergy_.clear(); hgcPassiveFHEnergy_.clear();
611  hgcPassiveBHEnergy_.clear(); hgcPassiveCMSEEnergy_.clear();
612  hgcPassiveEEName_.clear(); hgcPassiveFHName_.clear();
613  hgcPassiveBHName_.clear(); hgcPassiveCMSEName_.clear();
614  hgcPassiveEEID_.clear(); hgcPassiveFHID_.clear();
615  hgcPassiveBHID_.clear(); hgcPassiveCMSEID_.clear();
616 
618  edm::Handle<edm::PassiveHitContainer> hgcPHEE;
619  iEvent.getByToken(tok_hgcPHEE_,hgcPHEE);
620  analyzePassiveHits(hgcPHEE, 1);
621 
623  edm::Handle<edm::PassiveHitContainer> hgcPHFH;
624  iEvent.getByToken(tok_hgcPHFH_,hgcPHFH);
625  analyzePassiveHits(hgcPHFH, 2);
626 
628  edm::Handle<edm::PassiveHitContainer> hgcPHBH;
629  iEvent.getByToken(tok_hgcPHBH_,hgcPHBH);
630  analyzePassiveHits(hgcPHBH, 3);
631 
633  edm::Handle<edm::PassiveHitContainer> hgcPHCMSE;
634  iEvent.getByToken(tok_hgcPHCMSE_,hgcPHCMSE);
635  analyzePassiveHits(hgcPHCMSE, 4);
636  }
637 
638  if ((doSimHits_ || doPassive_) && (doTree_)) tree_->Fill();
639 
640  //Now the Digis
641  if (doDigis_) {
642  if (ifEE_) {
643  edm::Handle<HGCEEDigiCollection> theDigiContainers;
644  iEvent.getByToken(tok_digiEE_, theDigiContainers);
645  if (theDigiContainers.isValid()) {
646 #ifdef EDM_ML_DEBUG
647  std::cout << "HGCDigiCintainer for " << detectorEE_ << " with "
648  << theDigiContainers->size() << " element(s)" << std::endl;
649 #endif
650  for (auto it : *theDigiContainers) {
651  HGCEEDetId detId = (it.id());
652  HGCSample hgcSample = it.sample(sampleIndex_);
653  uint16_t adc = hgcSample.data();
654  analyzeDigi(0, detId, adc);
655  }
656  }
657  }
658  if (ifFH_) {
659  edm::Handle<HGCHEDigiCollection> theDigiContainers;
660  iEvent.getByToken(tok_digiFH_, theDigiContainers);
661  if (theDigiContainers.isValid()) {
662 #ifdef EDM_ML_DEBUG
663  std::cout << "HGCDigiContainer for " << detectorFH_ << " with "
664  << theDigiContainers->size() << " element(s)" << std::endl;
665 #endif
666  for (auto it : *theDigiContainers) {
667  HGCHEDetId detId = (it.id());
668  HGCSample hgcSample = it.sample(sampleIndex_);
669  uint16_t adc = hgcSample.data();
670  analyzeDigi(1, detId, adc);
671  }
672  }
673  }
674  }
675 
676  //The Rechits
677  if (doRecHits_) {
678  edm::Handle<HGCRecHitCollection> theCaloHitContainers;
679  if (ifEE_) {
680  iEvent.getByToken(tok_hitrEE_, theCaloHitContainers);
681  if (theCaloHitContainers.isValid()) {
682 #ifdef EDM_ML_DEBUG
683  std::cout << "HGCRecHitCollection for " << detectorEE_ << " has "
684  << theCaloHitContainers->size() << " hits" << std::endl;
685 #endif
686  analyzeRecHits(0, theCaloHitContainers);
687  } else {
688 #ifdef EDM_ML_DEBUG
689  std::cout << "HGCRecHitCollection does not exist for " << detectorEE_
690  << " !!!" << std::endl;
691 #endif
692  }
693  }
694  if (ifFH_) {
695  iEvent.getByToken(tok_hitrFH_, theCaloHitContainers);
696  if (theCaloHitContainers.isValid()) {
697 #ifdef EDM_ML_DEBUG
698  std::cout << "HGCRecHitCollection for " << detectorFH_ << " has "
699  << theCaloHitContainers->size() << " hits" << std::endl;
700 #endif
701  analyzeRecHits(1, theCaloHitContainers);
702  } else {
703 #ifdef EDM_ML_DEBUG
704  std::cout << "HGCRecHitCollection does not exist for " << detectorFH_
705  << " !!!" << std::endl;
706 #endif
707  }//else
708  }//if (ifFH_)
709  }//if (doRecHits_)
710 
711 }//void HGCalTBAnalyzer::analyze
712 
713 void HGCalTBAnalyzer::analyzeSimHits (int type, std::vector<PCaloHit>& hits,
714  double zFront) {
715 
716  std::map<uint32_t,double> map_hits, map_hitn;
717  std::map<int,double> map_hitDepth;
718  std::map<int,std::pair<uint32_t,double> > map_hitLayer, map_hitCell;
719  double entot(0);
720  for (unsigned int i=0; i<hits.size(); i++) {
721  double energy = hits[i].energy();
722  double time = hits[i].time();
723  uint32_t id = hits[i].id();
724  entot += energy;
725  int subdet, zside, layer, sector, subsector(0), cell, depth(0), idx(0);
726  if (type == 2) {
727  subdet = HcalDetId(id).subdet();
728  if (subdet != HcalOther) continue;
729  AHCalDetId hid(id);
730  layer = depth = hid.depth();
731  zside = hid.zside();
732  sector = hid.irow();
733  cell = hid.icol();
734  idx = ((hid.irowAbs()*100) + (hid.icolAbs()));
735  } else if (type == 3) {
736  HcalTestBeamNumbering::unpackIndex(id, subdet, layer, sector, cell);
737  depth = layer; zside = 1;
738  idx = subdet*1000 + layer;
739  layer = idx;
740  } else {
741  HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector,
742  subsector, cell);
743  depth = hgcons_[type]->simToReco(cell,layer,sector,true).second;
744  idx = sector*1000+cell;
745  }
746 #ifdef EDM_ML_DEBUG
747  std::cout << "SimHit:Hit[" << i << "] Id " << subdet << ":" << zside << ":"
748  << layer << ":" << sector << ":" << subsector << ":" << cell
749  << ":" << depth << " Energy " << energy << " Time " << time
750  << std::endl;
751 #endif
752  if (map_hits.count(id) != 0) {
753  map_hits[id] += energy;
754  } else {
755  map_hits[id] = energy;
756  }
757  if (map_hitLayer.count(layer) != 0) {
758  double ee = energy + map_hitLayer[layer].second;
759  map_hitLayer[layer] = std::pair<uint32_t,double>(id,ee);
760  } else {
761  map_hitLayer[layer] = std::pair<uint32_t,double>(id,energy);
762  }
763  if (depth >= 0) {
764  if (map_hitCell.count(idx) != 0) {
765  double ee = energy + map_hitCell[idx].second;
766  map_hitCell[idx] = std::pair<uint32_t,double>(id,ee);
767  } else {
768  map_hitCell[idx] = std::pair<uint32_t,double>(id,energy);
769  }
770  if (map_hitDepth.count(depth) != 0) {
771  map_hitDepth[depth] += energy;
772  } else {
773  map_hitDepth[depth] = energy;
774  }
775  uint32_t idn = (type >= 2) ? id :
776  HGCalTestNumbering::packHexagonIndex(subdet, zside, depth, sector,
777  subsector, cell);
778  if (map_hitn.count(idn) != 0) {
779  map_hitn[idn] += energy;
780  } else {
781  map_hitn[idn] = energy;
782  }
783  }
784  hSimHitT_[type]->Fill(time,energy);
785  }
786 
787  hSimHitEn_[type]->Fill(entot);
788  for (auto itr : map_hits) {
789  hSimHitE_[type]->Fill(itr.second);
790  }
791 
792  for (auto itr : map_hitLayer) {
793  int layer = itr.first - 1;
794  double energy = (itr.second).second;
795  double zp(0);
796  if (type < 2) zp = hgcons_[type]->waferZ(layer+1,false);
797  else if (type == 2) zp = AHCalDetId((itr.second).first).getZ();
798 #ifdef EDM_ML_DEBUG
799  std::cout << "SimHit:Layer " << layer+1 << " Z " << zp << ":" << zp-zFront
800  << " E " << energy << std::endl;
801 #endif
802  if (type < 3) {
803  hSimHitLng_[type]->Fill(zp-zFront,energy);
804  hSimHitLng2_[type]->Fill(layer+1,energy);
805  }
806  if (type == 0) {
807  if (layer < (int)(hSimHitLayEn1EE_.size())) {
808  simHitLayEn1EE_[layer] = energy;
809  hSimHitLayEn1EE_[layer]->Fill(energy);
810  }
811  } else if (type == 1) {
812  if (layer < (int)(hSimHitLayEn1FH_.size())) {
813  simHitLayEn1FH_[layer] = energy;
814  hSimHitLayEn1FH_[layer]->Fill(energy);
815  }
816  } else if (type == 2) {
817  if (layer < (int)(hSimHitLayEn1BH_.size())) {
818  simHitLayEn1BH_[layer] = energy;
819  hSimHitLayEn1BH_[layer]->Fill(energy);
820  }
821  } else {
822  for (unsigned int k=0; k<idBeams_.size(); ++k) {
823  if (layer+1 == idBeams_[k]) {
824  simHitLayEnBeam_[k] = energy;
825  hSimHitLayEnBeam_[k]->Fill(energy);
826  break;
827  }
828  }
829  }
830  }
831  for (auto itr : map_hitDepth) {
832  int layer = itr.first - 1;
833  double energy = itr.second;
834 #ifdef EDM_ML_DEBUG
835  std::cout << "SimHit:Layer " << layer+1 << " " << energy << std::endl;
836 #endif
837  hSimHitLng1_[type]->Fill(layer+1,energy);
838  if (type == 0) {
839  if (layer < (int)(hSimHitLayEn2EE_.size())) {
840  simHitLayEn2EE_[layer] = energy;
841  hSimHitLayEn2EE_[layer]->Fill(energy);
842  }
843  } else if (type == 1) {
844  if (layer < (int)(hSimHitLayEn2FH_.size())) {
845  simHitLayEn2FH_[layer] = energy;
846  hSimHitLayEn2FH_[layer]->Fill(energy);
847  }
848  } else if (type == 2) {
849  if (layer < (int)(hSimHitLayEn2BH_.size())) {
850  simHitLayEn2BH_[layer] = energy;
851  hSimHitLayEn2BH_[layer]->Fill(energy);
852  }
853  }
854  }
855 
856  if (type < 3) {
857  for (auto itr : map_hitCell) {
858  uint32_t id = ((itr.second).first);
859  double energy = ((itr.second).second);
860  std::pair<float,float> xy(0,0);
861  double xx(0);
862  if (type == 2) {
863  xy = AHCalDetId(id).getXY();
864  xx = xy.first;
865  } else {
866  int subdet, zside, layer, sector, subsector, cell;
867  HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector,
868  subsector, cell);
869  xy = hgcons_[type]->locateCell(cell,layer,sector,false);
870  double zp = hgcons_[type]->waferZ(layer,false);
871  xx = (zp < 0) ? -xy.first : xy.first;
872  }
873  hSimHitLat_[type]->Fill(xx,xy.second,energy);
874  }
875  }
876 
877  for (auto itr : map_hitn) {
878  uint32_t id = itr.first;
879  double energy = itr.second;
880  if (type == 0) {
881  simHitCellIdEE_.push_back(id); simHitCellEnEE_.push_back(energy);
882  } else if (type == 1) {
883  simHitCellIdFH_.push_back(id); simHitCellEnFH_.push_back(energy);
884  } else if (type == 2) {
885  simHitCellIdBH_.push_back(id); simHitCellEnBH_.push_back(energy);
886  } else if (type == 3) {
887  simHitCellIdBeam_.push_back(id); simHitCellEnBeam_.push_back(energy);
888  }
889  }
890 }
891 
892 void HGCalTBAnalyzer::analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
893  edm::Handle<edm::SimVertexContainer> const& SimVtx) {
894 
895  xBeam_ = yBeam_ = zBeam_ = pBeam_ = -1000000;
896  int vertIndex(-1);
897  for (auto simTrkItr : *SimTk) {
898 #ifdef EDM_ML_DEBUG
899  std::cout << "Track " << simTrkItr.trackId() << " Vertex "
900  << simTrkItr.vertIndex() << " Type " << simTrkItr.type()
901  << " Charge " << simTrkItr.charge() << " momentum "
902  << simTrkItr.momentum() << " " << simTrkItr.momentum().P()
903  << std::endl;
904 #endif
905  if (vertIndex == -1) {
906  vertIndex = simTrkItr.vertIndex();
907  pBeam_ = simTrkItr.momentum().P();
908  }
909  }
910  if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
911  edm::SimVertexContainer::const_iterator simVtxItr= SimVtx->begin();
912  for (int iv=0; iv<vertIndex; iv++) simVtxItr++;
913 #ifdef EDM_ML_DEBUG
914  std::cout << "Vertex " << vertIndex << " position "
915  << simVtxItr->position() << std::endl;
916 #endif
917  xBeam_ = simVtxItr->position().X();
918  yBeam_ = simVtxItr->position().Y();
919  zBeam_ = simVtxItr->position().Z();
920  }
921 
922 }
923 
924 template<class T1>
925 void HGCalTBAnalyzer::analyzeDigi (int type, const T1& detId, uint16_t adc) {
926 
927  DetId id1 = DetId(detId.rawId());
928  GlobalPoint global = hgeom_[type]->getPosition(id1);
929  hDigiOcc_[type]->Fill(global.x(),global.y());
930  hDigiLng_[type]->Fill(global.z());
931  hDigiADC_[type]->Fill(adc);
932 }
933 
934 void HGCalTBAnalyzer::analyzeRecHits (int type,
935  edm::Handle<HGCRecHitCollection>& hits) {
936 
937  std::map<int,double> map_hitLayer;
938  std::map<int,std::pair<DetId,double> > map_hitCell;
939  for (auto it : *hits) {
940  DetId detId = it.id();
941  GlobalPoint global = hgeom_[type]->getPosition(detId);
942  double energy = it.energy();
943  int layer = HGCalDetId(detId).layer();
944  int cell = HGCalDetId(detId).cell();
945  hRecHitOcc_[type]->Fill(global.x(),global.y(),energy);
946  hRecHitE_[type]->Fill(energy);
947  if (map_hitLayer.count(layer) != 0) {
948  map_hitLayer[layer] += energy;
949  } else {
950  map_hitLayer[layer] = energy;
951  }
952  if (map_hitCell.count(cell) != 0) {
953  double ee = energy + map_hitCell[cell].second;
954  map_hitCell[cell] = std::pair<uint32_t,double>(detId,ee);
955  } else {
956  map_hitCell[cell] = std::pair<uint32_t,double>(detId,energy);
957  }
958 #ifdef EDM_ML_DEBUG
959  std::cout << "RecHit: " << layer << " " << global.x() << " " << global.y()
960  << " " << global.z() << " " << energy << std::endl;
961 #endif
962  }
963 
964  for (auto itr : map_hitLayer) {
965  int layer = itr.first;
966  double energy = itr.second;
967  double zp = hgcons_[type]->waferZ(layer,true);
968 #ifdef EDM_ML_DEBUG
969  std::cout << "SimHit:Layer " << layer << " " << zp << " " << energy
970  << std::endl;
971 #endif
972  hRecHitLng_[type]->Fill(zp,energy);
973  hRecHitLng1_[type]->Fill(layer,energy);
974  }
975 
976  for (auto itr : map_hitCell) {
977  DetId detId = ((itr.second).first);
978  double energy = ((itr.second).second);
979  GlobalPoint global = hgeom_[type]->getPosition(detId);
980  hRecHitLat_[type]->Fill(global.x(),global.y(),energy);
981  }
982 }
983 
984 void HGCalTBAnalyzer::analyzePassiveHits (edm::Handle<edm::PassiveHitContainer>const& hgcPH, int subdet) {
985 
986  for (auto v : *hgcPH) {
987  double energy = v.energy();
988  std::string name = v.vname();
989  unsigned int id = v.id();
990 #ifdef EDM_ML_DEBUG
991  double time = v.time();
992  std::cout << "HGCalTBAnalyzer::analyzePassiveHits:Energy:Time:Name:Id : "
993  << energy << ":" << time << ":" << name << ":" << id
994  << std::endl;
995 #endif
996 
997  if (subdet==1) {
998  hgcPassiveEEEnergy_.push_back(energy);
999  hgcPassiveEEName_.push_back(name);
1000  hgcPassiveEEID_.push_back(id);
1001  } else if (subdet==2) {
1002  hgcPassiveFHEnergy_.push_back(energy);
1003  hgcPassiveFHName_.push_back(name);
1004  hgcPassiveFHID_.push_back(id);
1005  } else if (subdet==3) {
1006  hgcPassiveBHEnergy_.push_back(energy);
1007  hgcPassiveBHName_.push_back(name);
1008  hgcPassiveBHID_.push_back(id);
1009  } else if (subdet==4) {
1010  hgcPassiveCMSEEnergy_.push_back(energy);
1011  hgcPassiveCMSEName_.push_back(name);
1012  hgcPassiveCMSEID_.push_back(id);
1013  }
1014  }
1015 }
1016 
1017 //define this as a plug-in
1018 DEFINE_FWK_MODULE(HGCalTBAnalyzer);
HGCalTBAnalyzer::hSimHitLayEn2EE_
std::vector< TH1D * > hSimHitLayEn2EE_
Definition: HGCalTBAnalyzer.cc:107
ecalMGPA::adc
int adc(sample_type sample)
get the ADC sample (12 bits)
Definition: EcalMGPASample.h:11
HGCalTBAnalyzer::detectorEE_
std::string detectorEE_
Definition: HGCalTBAnalyzer.cc:83
HGCalTBAnalyzer::simHitLayEn2EE_
std::vector< float > simHitLayEn2EE_
Definition: HGCalTBAnalyzer.cc:111
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
HGCSample::data
uint32_t data() const
Definition: HGCSample.h:60
HGCalTBAnalyzer::simHitLayEn2BH_
std::vector< float > simHitLayEn2BH_
Definition: HGCalTBAnalyzer.cc:113
HGCalTBAnalyzer::simHitCellIdFH_
std::vector< uint32_t > simHitCellIdFH_
Definition: HGCalTBAnalyzer.cc:115
mps_fire.i
i
Definition: mps_fire.py:156
HGCalTBAnalyzer::detectorBH_
std::string detectorBH_
Definition: HGCalTBAnalyzer.cc:84
HGCalTBAnalyzer::analyzeSimTracks
void analyzeSimTracks(edm::Handle< edm::SimTrackContainer > const &SimTk, edm::Handle< edm::SimVertexContainer > const &SimVtx)
Definition: HGCalTBAnalyzer.cc:892
HcalDetId::subdet
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:49
edm::ParameterSetDescription::addUntracked
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:94
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
AHCalDetId::icolAbs
int icolAbs() const
Definition: AHCalDetId.h:39
HGCalTBAnalyzer::tok_hitrBH_
edm::EDGetToken tok_hitrBH_
Definition: HGCalTBAnalyzer.cc:93
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
HGCalTBAnalyzer::hSimHitLng1_
TProfile * hSimHitLng1_[3]
Definition: HGCalTBAnalyzer.cc:103
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
HGCalTBAnalyzer::hSimHitLng2_
TProfile * hSimHitLng2_[3]
Definition: HGCalTBAnalyzer.cc:104
HGCalTBAnalyzer::hRecHitLng_
TProfile * hRecHitLng_[3]
Definition: HGCalTBAnalyzer.cc:105
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
HGCalTBAnalyzer::hgcPassiveFHName_
std::vector< std::string > hgcPassiveFHName_
Definition: HGCalTBAnalyzer.cc:121
HGCalTBAnalyzer::simHitLayEn1EE_
std::vector< float > simHitLayEn1EE_
Definition: HGCalTBAnalyzer.cc:111
HGCalTBAnalyzer::hgcPassiveBHEnergy_
std::vector< float > hgcPassiveBHEnergy_
Definition: HGCalTBAnalyzer.cc:120
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
AHCalDetId::MaxDepth
static const int MaxDepth
get the layer number
Definition: AHCalDetId.h:41
HGCalTestNumbering::packHexagonIndex
static uint32_t packHexagonIndex(int subdet, int z, int lay, int wafer, int celltyp, int cell)
Definition: HGCalTestNumbering.cc:23
HGCalTBAnalyzer::detectorFH_
std::string detectorFH_
Definition: HGCalTBAnalyzer.cc:83
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
HGCalTBAnalyzer::fs_
edm::Service< TFileService > fs_
Definition: HGCalTBAnalyzer.cc:76
HGCalTBAnalyzer::HGCalTBAnalyzer
HGCalTBAnalyzer(edm::ParameterSet const &)
Definition: HGCalTBAnalyzer.cc:127
HGCalTBAnalyzer::beginRun
virtual void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: HGCalTBAnalyzer.cc:396
HGCalTBAnalyzer::tok_digiFH_
edm::EDGetToken tok_digiFH_
Definition: HGCalTBAnalyzer.cc:92
HGCalTBAnalyzer::tok_hitsFH_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsFH_
Definition: HGCalTBAnalyzer.cc:88
AHCalDetId::getZ
double getZ() const
Definition: AHCalDetId.cc:57
AHCalDetId::irow
int irow() const
get the row number
Definition: AHCalDetId.cc:33
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
HGCalGeometry::getPosition
GlobalPoint getPosition(const DetId &id) const
Definition: HGCalGeometry.cc:144
HGCalTBAnalyzer::hgcPassiveBHName_
std::vector< std::string > hgcPassiveBHName_
Definition: HGCalTBAnalyzer.cc:121
HGCalDDDConstants::locateCell
std::pair< float, float > locateCell(int cell, int lay, int type, bool reco) const
Definition: HGCalDDDConstants.cc:348
HGCalTBAnalyzer::hgcPassiveCMSEEnergy_
std::vector< float > hgcPassiveCMSEEnergy_
Definition: HGCalTBAnalyzer.cc:120
HGCalTBAnalyzer::hgcPassiveEEID_
std::vector< int > hgcPassiveEEID_
Definition: HGCalTBAnalyzer.cc:122
HGCalTBAnalyzer::hgcPassiveEEEnergy_
std::vector< float > hgcPassiveEEEnergy_
Definition: HGCalTBAnalyzer.cc:120
HGCSample
wrapper for a data word
Definition: HGCSample.h:13
AHCalDetId::depth
int depth() const
Definition: AHCalDetId.cc:45
HGCalTBAnalyzer::simHitCellEnFH_
std::vector< float > simHitCellEnFH_
Definition: HGCalTBAnalyzer.cc:117
HGCalTBAnalyzer::hSimHitLat_
TProfile2D * hSimHitLat_[3]
Definition: HGCalTBAnalyzer.cc:106
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:250
HGCalTBAnalyzer::hSimHitLayEn1EE_
std::vector< TH1D * > hSimHitLayEn1EE_
Definition: HGCalTBAnalyzer.cc:107
HGCalTBAnalyzer::tok_hitsBH_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsBH_
Definition: HGCalTBAnalyzer.cc:89
HGCalTBAnalyzer::analyzeSimHits
void analyzeSimHits(int type, std::vector< PCaloHit > &hits, double zFront)
Definition: HGCalTBAnalyzer.cc:713
AHCalDetId::icol
int icol() const
get the column number
Definition: AHCalDetId.cc:39
HGCalDDDConstants::simToReco
std::pair< int, int > simToReco(int cell, int layer, int mod, bool half) const
Definition: HGCalDDDConstants.cc:599
iEvent
int iEvent
Definition: GenABIO.cc:230
HGCalTBAnalyzer::hSimHitLayEn2BH_
std::vector< TH1D * > hSimHitLayEn2BH_
Definition: HGCalTBAnalyzer.cc:109
particleFlowClusterECALTimeSelected_cfi.depth
depth
Definition: particleFlowClusterECALTimeSelected_cfi.py:8
HGCalDDDConstants::layers
unsigned int layers(bool reco) const
Definition: HGCalDDDConstants.cc:340
HGCalTBAnalyzer::idBeams_
std::vector< int > idBeams_
Definition: HGCalTBAnalyzer.cc:87
HGCalTBAnalyzer::analyzeRecHits
void analyzeRecHits(int type, edm::Handle< HGCRecHitCollection > &hits)
Definition: HGCalTBAnalyzer.cc:934
HGCalTBAnalyzer::hgeom_
const HGCalGeometry * hgeom_[2]
Definition: HGCalTBAnalyzer.cc:78
HGCalTBAnalyzer::hgcPassiveCMSEName_
std::vector< std::string > hgcPassiveCMSEName_
Definition: HGCalTBAnalyzer.cc:121
AHCalDetId::getXY
std::pair< double, double > getXY() const
get the local coordinate in the plane and along depth
Definition: AHCalDetId.cc:49
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
HGCalTBAnalyzer::tok_simTk_
edm::EDGetTokenT< edm::SimTrackContainer > tok_simTk_
Definition: HGCalTBAnalyzer.cc:90
HGCalTBAnalyzer::analyze
virtual void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: HGCalTBAnalyzer.cc:486
HGCalTBAnalyzer::hSimHitLayEnBeam_
std::vector< TH1D * > hSimHitLayEnBeam_
Definition: HGCalTBAnalyzer.cc:110
HGCalTBAnalyzer::hgcons_
const HGCalDDDConstants * hgcons_[2]
Definition: HGCalTBAnalyzer.cc:77
HGCalTBAnalyzer::tok_hitrFH_
edm::EDGetToken tok_hitrFH_
Definition: HGCalTBAnalyzer.cc:93
HGCalTBAnalyzer::tok_hitsBeam_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsBeam_
Definition: HGCalTBAnalyzer.cc:89
HGCalTBAnalyzer::tok_hepMC_
edm::EDGetTokenT< edm::HepMCProduct > tok_hepMC_
Definition: HGCalTBAnalyzer.cc:94
edm::ParameterSetDescription::add
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:89
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:74
HGCalTBAnalyzer::hSimHitLayEn1FH_
std::vector< TH1D * > hSimHitLayEn1FH_
Definition: HGCalTBAnalyzer.cc:108
HGCalTBAnalyzer::simHitCellIdBH_
std::vector< uint32_t > simHitCellIdBH_
Definition: HGCalTBAnalyzer.cc:116
HGCalTBAnalyzer::analyzePassiveHits
void analyzePassiveHits(edm::Handle< edm::PassiveHitContainer > const &hgcPh, int subdet)
Definition: HGCalTBAnalyzer.cc:984
HGCalDDDConstants::waferZ
double waferZ(int layer, bool reco) const
Definition: HGCalDDDConstants.cc:727
HGCalTBAnalyzer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HGCalTBAnalyzer.cc:233
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
DetId
Definition: DetId.h:18
HcalDigiParam_cfi.zside
zside
Definition: HcalDigiParam_cfi.py:5
HGCalTBAnalyzer::hgcPassiveFHEnergy_
std::vector< float > hgcPassiveFHEnergy_
Definition: HGCalTBAnalyzer.cc:120
HGCalTBAnalyzer::tok_hgcPHEE_
edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHEE_
Definition: HGCalTBAnalyzer.cc:95
HGCalTBAnalyzer::hgcPassiveCMSEID_
std::vector< int > hgcPassiveCMSEID_
Definition: HGCalTBAnalyzer.cc:122
HGCalTBAnalyzer::hRecHitLat_
TProfile2D * hRecHitLat_[3]
Definition: HGCalTBAnalyzer.cc:106
edm::HepMCProduct::GetEvent
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:38
HGCalTBAnalyzer::tok_hgcPHBH_
edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHBH_
Definition: HGCalTBAnalyzer.cc:96
HGCalDetId::cell
int cell() const
get the absolute value of the cell #&#39;s in x and y
Definition: HGCalDetId.h:39
HGCalTBAnalyzer::simHitCellEnBeam_
std::vector< float > simHitCellEnBeam_
Definition: HGCalTBAnalyzer.cc:118
HGCalTBAnalyzer::hSimHitLng_
TProfile * hSimHitLng_[3]
Definition: HGCalTBAnalyzer.cc:103
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
training_settings.idx
idx
Definition: training_settings.py:16
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:63
HGCalTBAnalyzer::hgcPassiveEEName_
std::vector< std::string > hgcPassiveEEName_
Definition: HGCalTBAnalyzer.cc:121
HGCalTBAnalyzer::detectorBeam_
std::string detectorBeam_
Definition: HGCalTBAnalyzer.cc:84
HGCalTBAnalyzer::analyzeDigi
void analyzeDigi(int type, const T1 &detId, uint16_t adc)
Definition: HGCalTBAnalyzer.cc:925
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
fftjetcommon_cfi.title
title
Definition: fftjetcommon_cfi.py:32
AHCalDetId::irowAbs
int irowAbs() const
Definition: AHCalDetId.h:35
edm::SortedCollection::size
size_type size() const
Definition: SortedCollection.h:236
HGCalTBAnalyzer::hgcPassiveBHID_
std::vector< int > hgcPassiveBHID_
Definition: HGCalTBAnalyzer.cc:122
HGCalTBAnalyzer::endRun
virtual void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: HGCalTBAnalyzer.cc:66
HGCalTBAnalyzer::tok_hitsEE_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hitsEE_
Definition: HGCalTBAnalyzer.cc:88
HGCalTBAnalyzer::simHitLayEn1FH_
std::vector< float > simHitLayEn1FH_
Definition: HGCalTBAnalyzer.cc:112
HGCalTBAnalyzer::tok_digiBH_
edm::EDGetToken tok_digiBH_
Definition: HGCalTBAnalyzer.cc:92
HGCalTBAnalyzer::hgcPassiveFHID_
std::vector< int > hgcPassiveFHID_
Definition: HGCalTBAnalyzer.cc:122
HGCalTBAnalyzer::simHitCellEnBH_
std::vector< float > simHitCellEnBH_
Definition: HGCalTBAnalyzer.cc:118
HcalTestBeamNumbering::unpackIndex
static void unpackIndex(const uint32_t &idx, int &det, int &lay, int &x, int &y)
Definition: HcalTestBeamNumbering.cc:28
gather_cfg.cout
cout
Definition: gather_cfg.py:145
HGCalTBAnalyzer::tok_hgcPHFH_
edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHFH_
Definition: HGCalTBAnalyzer.cc:95
HGCalTBAnalyzer::simHitLayEn1BH_
std::vector< float > simHitLayEn1BH_
Definition: HGCalTBAnalyzer.cc:113
HGCalTBAnalyzer::simHitCellIdEE_
std::vector< uint32_t > simHitCellIdEE_
Definition: HGCalTBAnalyzer.cc:115
HGCalTBAnalyzer::hSimHitLayEn1BH_
std::vector< TH1D * > hSimHitLayEn1BH_
Definition: HGCalTBAnalyzer.cc:109
AHCalDetId::zside
int zside() const
get the z-side of the cell (1/-1)
Definition: AHCalDetId.h:32
HGCalTBAnalyzer::simHitCellEnEE_
std::vector< float > simHitCellEnEE_
Definition: HGCalTBAnalyzer.cc:117
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
HGCalTBAnalyzer::simHitCellIdBeam_
std::vector< uint32_t > simHitCellIdBeam_
Definition: HGCalTBAnalyzer.cc:116
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
HGCalTBAnalyzer::tok_hitrEE_
edm::EDGetToken tok_hitrEE_
Definition: HGCalTBAnalyzer.cc:93
HGCalTBAnalyzer::simHitLayEn2FH_
std::vector< float > simHitLayEn2FH_
Definition: HGCalTBAnalyzer.cc:112
HGCalTBAnalyzer::tok_simVtx_
edm::EDGetTokenT< edm::SimVertexContainer > tok_simVtx_
Definition: HGCalTBAnalyzer.cc:91
HGCalTestNumbering::unpackHexagonIndex
static void unpackHexagonIndex(const uint32_t &idx, int &subdet, int &z, int &lay, int &wafer, int &celltyp, int &cell)
Definition: HGCalTestNumbering.cc:53
hcalTTPDigis_cfi.id
id
Definition: hcalTTPDigis_cfi.py:7
HGCalDetId::layer
int layer() const
get the layer #
Definition: HGCalDetId.h:48
HGCalTBAnalyzer::hRecHitLng1_
TProfile * hRecHitLng1_[3]
Definition: HGCalTBAnalyzer.cc:105
edm::Run
Definition: Run.h:42
HGCalTBAnalyzer::hSimHitLayEn2FH_
std::vector< TH1D * > hSimHitLayEn2FH_
Definition: HGCalTBAnalyzer.cc:108
HGCalTBAnalyzer::beginJob
virtual void beginJob() override
Definition: HGCalTBAnalyzer.cc:276
HGCalTBAnalyzer::tok_digiEE_
edm::EDGetToken tok_digiEE_
Definition: HGCalTBAnalyzer.cc:92
HGCalTBAnalyzer::tok_hgcPHCMSE_
edm::EDGetTokenT< edm::PassiveHitContainer > tok_hgcPHCMSE_
Definition: HGCalTBAnalyzer.cc:96
HGCalTBAnalyzer::simHitLayEnBeam_
std::vector< float > simHitLayEnBeam_
Definition: HGCalTBAnalyzer.cc:114
Generated for CMSSW Reference Manual by doxygen 1.8.11