Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-28 23:36:04

0001 // system include files
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <map>
0006 #include <memory>
0007 #include <string>
0008 #include <vector>
0009 
0010 // user include files
0011 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0012 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0013 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0014 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0015 
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ServiceRegistry/interface/Service.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026 
0027 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0028 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0031 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0032 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0033 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0034 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0035 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
0036 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0037 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0038 #include "SimG4CMS/HGCalTestBeam/interface/AHCalDetId.h"
0039 #include "SimG4CMS/HGCalTestBeam/interface/AHCalGeometry.h"
0040 
0041 #include "SimDataFormats/CaloHit/interface/PassiveHit.h"
0042 
0043 // Root objects
0044 #include "TH1.h"
0045 #include "TH2.h"
0046 #include "TProfile.h"
0047 #include "TProfile2D.h"
0048 #include "TROOT.h"
0049 #include "TSystem.h"
0050 #include "TTree.h"
0051 
0052 //#define EDM_ML_DEBUG
0053 
0054 class HGCalTBAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0055 public:
0056   explicit HGCalTBAnalyzer(edm::ParameterSet const&);
0057   ~HGCalTBAnalyzer() override;
0058 
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060 
0061 private:
0062   void beginJob() override;
0063   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0064   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0065   void analyze(edm::Event const&, edm::EventSetup const&) override;
0066   void analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront);
0067   void analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
0068                         edm::Handle<edm::SimVertexContainer> const& SimVtx);
0069   template <class T1>
0070   void analyzeDigi(int type, const T1& detId, uint16_t adc);
0071   void analyzeRecHits(int type, const edm::Handle<HGCRecHitCollection>& hits);
0072   void analyzePassiveHits(edm::Handle<edm::PassiveHitContainer> const& hgcPh, int subdet);
0073   static bool sortTime(const std::pair<double, double>& i, const std::pair<double, double>& j);
0074 
0075   edm::Service<TFileService> fs_;
0076   std::unique_ptr<AHCalGeometry> ahcalGeom_;
0077   const HGCalDDDConstants* hgcons_[2];
0078   const HGCalGeometry* hgeom_[2];
0079   const bool ifEE_, ifFH_, ifBH_, ifBeam_;
0080   const bool doSimHits_, doDigis_, doRecHits_;
0081   const bool doTree_, doTreeCell_;
0082   const bool doPassive_, doPassiveEE_, doPassiveHE_, doPassiveBH_, addP_, doBeam_;
0083   const std::string detectorEE_, detectorFH_;
0084   const std::string detectorBH_, detectorBeam_;
0085   const double zFrontEE_, zFrontFH_, zFrontBH_;
0086   const int sampleIndex_;
0087   const double gev2mip200_, gev2mip300_, stoc_smear_time_200_, stoc_smear_time_300_;
0088   std::vector<int> idBeams_;
0089   const edm::InputTag labelGen_;
0090   const std::string labelHitEE_, labelHitFH_, labelHitBH_, labelHitBeam_;
0091   const edm::InputTag labelDigiEE_, labelDigiFH_, labelDigiBH_;
0092   const edm::InputTag labelRHitEE_, labelRHitFH_, labelRHitBH_;
0093   const edm::InputTag labelPassiveEE_, labelPassiveFH_, labelPassiveBH_;
0094   const edm::InputTag labelPassiveCMSE_, labelPassiveBeam_;
0095 
0096   const edm::EDGetTokenT<edm::HepMCProduct> tok_hepMC_;
0097   const edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0098   const edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0099   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsEE_, tok_hitsFH_;
0100   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsBH_, tok_hitsBeam_;
0101   const edm::EDGetTokenT<HGCalDigiCollection> tok_digiEE_, tok_digiFH_, tok_digiBH_;
0102   const edm::EDGetTokenT<HGCRecHitCollection> tok_hitrEE_, tok_hitrFH_, tok_hitrBH_;
0103   const edm::EDGetTokenT<edm::PassiveHitContainer> tok_hgcPHEE_, tok_hgcPHFH_;
0104   const edm::EDGetTokenT<edm::PassiveHitContainer> tok_hgcPHBH_, tok_hgcPHCMSE_, tok_hgcPHBeam_;
0105   const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tokDDDEE_;
0106   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> tokGeomEE_;
0107   const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tokDDDFH_;
0108   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> tokGeomFH_;
0109 
0110   TTree* tree_;
0111   TH1D *hSimHitE_[4], *hSimHitT_[4];
0112   TH1D *hDigiADC_[3], *hDigiLng_[2];
0113   TH1D *hRecHitE_[3], *hSimHitEn_[4], *hBeam_;
0114   TH2D *hDigiOcc_[3], *hRecHitOcc_[3];
0115   TProfile *hSimHitLng_[3], *hSimHitLng1_[3];
0116   TProfile* hSimHitLng2_[3];
0117   TProfile *hRecHitLng_[3], *hRecHitLng1_[3];
0118   TProfile2D *hSimHitLat_[3], *hRecHitLat_[3];
0119   std::vector<TH1D*> hSimHitLayEn1EE_, hSimHitLayEn2EE_;
0120   std::vector<TH1D*> hSimHitLayEn1FH_, hSimHitLayEn2FH_;
0121   std::vector<TH1D*> hSimHitLayEn1BH_, hSimHitLayEn2BH_;
0122   std::vector<TH1D*> hSimHitLayEnBeam_;
0123   std::vector<float> simHitLayEn1EE_, simHitLayEn2EE_;
0124   std::vector<float> simHitLayEn1FH_, simHitLayEn2FH_;
0125   std::vector<float> simHitLayEn1BH_, simHitLayEn2BH_;
0126   std::vector<float> simHitLayEnBeam_;
0127   std::vector<uint32_t> simHitCellIdEE_, simHitCellIdFH_;
0128   std::vector<uint32_t> simHitCellIdBH_, simHitCellIdBeam_;
0129   std::vector<float> simHitCellEnEE_, simHitCellEnFH_;
0130   std::vector<float> simHitCellEnBH_, simHitCellEnBeam_;
0131   std::vector<int> simHitCellColBH_, simHitCellRowBH_, simHitCellLayerBH_;
0132   std::vector<float> simHitCellTimeFirstHitEE_, simHitCellTimeFirstHitFH_, simHitCellTimeFirstHitBH_;
0133   std::vector<float> simHitCellTime15MipEE_, simHitCellTime15MipFH_, simHitCellTime15MipBH_;
0134   std::vector<float> simHitCellTimeLastHitEE_, simHitCellTimeLastHitFH_, simHitCellTimeLastHitBH_;
0135 
0136   std::vector<float> hgcPassiveEEEnergy_, hgcPassiveFHEnergy_, hgcPassiveBHEnergy_, hgcPassiveCMSEEnergy_;
0137   std::vector<float> hgcPassiveBeamEnergy_;
0138   std::vector<std::string> hgcPassiveEEName_, hgcPassiveFHName_, hgcPassiveBHName_, hgcPassiveCMSEName_;
0139   std::vector<std::string> hgcPassiveBeamName_;
0140   std::vector<int> hgcPassiveEEID_, hgcPassiveFHID_, hgcPassiveBHID_, hgcPassiveCMSEID_, hgcPassiveBeamID_;
0141 
0142   double xBeam_, yBeam_, zBeam_, pBeam_;
0143   double thetaBeam_, phiBeam_;
0144   int nBeamMC_;
0145   std::vector<int> pdgIdBeamMC_;
0146   std::vector<float> xBeamMC_, yBeamMC_, zBeamMC_;
0147   std::vector<float> pxBeamMC_, pyBeamMC_, pzBeamMC_, pBeamMC_;
0148 };
0149 
0150 HGCalTBAnalyzer::HGCalTBAnalyzer(const edm::ParameterSet& iConfig)
0151     : ifEE_(iConfig.getParameter<bool>("useEE")),
0152       ifFH_(iConfig.getParameter<bool>("useFH")),
0153       ifBH_(iConfig.getParameter<bool>("useBH")),
0154       ifBeam_(iConfig.getParameter<bool>("useBeam")),
0155       doSimHits_(iConfig.getParameter<bool>("doSimHits")),
0156       doDigis_(iConfig.getParameter<bool>("doDigis")),
0157       doRecHits_(iConfig.getParameter<bool>("doRecHits")),
0158       doTree_(iConfig.getParameter<bool>("doTree")),
0159       doTreeCell_(iConfig.getParameter<bool>("doTreeCell")),
0160       doPassive_(iConfig.getParameter<bool>("doPassive")),
0161       doPassiveEE_(iConfig.getParameter<bool>("doPassiveEE")),
0162       doPassiveHE_(iConfig.getParameter<bool>("doPassiveHE")),
0163       doPassiveBH_(iConfig.getParameter<bool>("doPassiveBH")),
0164       addP_(iConfig.getParameter<bool>("addP")),
0165       doBeam_(iConfig.getParameter<bool>("doBeam")),
0166       detectorEE_(iConfig.getParameter<std::string>("detectorEE")),
0167       detectorFH_(iConfig.getParameter<std::string>("detectorFH")),
0168       detectorBH_(iConfig.getParameter<std::string>("detectorBH")),
0169       detectorBeam_(iConfig.getParameter<std::string>("detectorBeam")),
0170       zFrontEE_(iConfig.getParameter<double>("zFrontEE")),
0171       zFrontFH_(iConfig.getParameter<double>("zFrontFH")),
0172       zFrontBH_(iConfig.getParameter<double>("zFrontBH")),
0173       sampleIndex_(iConfig.getParameter<int>("sampleIndex")),
0174       gev2mip200_(iConfig.getUntrackedParameter<double>("gev2mip200", 57.0e-6)),
0175       gev2mip300_(iConfig.getUntrackedParameter<double>("gev2mip300", 85.5e-6)),
0176       stoc_smear_time_200_(iConfig.getUntrackedParameter<double>("stoc_smear_time_200", 10.24)),
0177       stoc_smear_time_300_(iConfig.getUntrackedParameter<double>("stoc_smear_time_300", 15.5)),
0178       labelGen_(iConfig.getParameter<edm::InputTag>("generatorSrc")),
0179       labelHitEE_(iConfig.getParameter<std::string>("caloHitSrcEE")),
0180       labelHitFH_(iConfig.getParameter<std::string>("caloHitSrcFH")),
0181       labelHitBH_(iConfig.getParameter<std::string>("caloHitSrcBH")),
0182       labelHitBeam_(iConfig.getParameter<std::string>("caloHitSrcBeam")),
0183       labelDigiEE_(iConfig.getParameter<edm::InputTag>("digiSrcEE")),
0184       labelDigiFH_(iConfig.getParameter<edm::InputTag>("digiSrcFH")),
0185       labelDigiBH_(iConfig.getParameter<edm::InputTag>("digiSrcBH")),
0186       labelRHitEE_(iConfig.getParameter<edm::InputTag>("recHitSrcEE")),
0187       labelRHitFH_(iConfig.getParameter<edm::InputTag>("recHitSrcFH")),
0188       labelRHitBH_(iConfig.getParameter<edm::InputTag>("recHitSrcBH")),
0189       labelPassiveEE_(iConfig.getParameter<edm::InputTag>("passiveEE")),
0190       labelPassiveFH_(iConfig.getParameter<edm::InputTag>("passiveFH")),
0191       labelPassiveBH_(iConfig.getParameter<edm::InputTag>("passiveBH")),
0192       labelPassiveCMSE_(iConfig.getParameter<edm::InputTag>("passiveCMSE")),
0193       labelPassiveBeam_(iConfig.getParameter<edm::InputTag>("passiveBeam")),
0194       tok_hepMC_(consumes<edm::HepMCProduct>(labelGen_)),
0195       tok_simTk_(consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"))),
0196       tok_simVtx_(consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"))),
0197       tok_hitsEE_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitEE_))),
0198       tok_hitsFH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitFH_))),
0199       tok_hitsBH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBH_))),
0200       tok_hitsBeam_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBeam_))),
0201       tok_digiEE_(consumes<HGCalDigiCollection>(labelDigiEE_)),
0202       tok_digiFH_(consumes<HGCalDigiCollection>(labelDigiFH_)),
0203       tok_digiBH_(consumes<HGCalDigiCollection>(labelDigiFH_)),
0204       tok_hitrEE_(consumes<HGCRecHitCollection>(labelRHitEE_)),
0205       tok_hitrFH_(consumes<HGCRecHitCollection>(labelRHitFH_)),
0206       tok_hitrBH_(consumes<HGCRecHitCollection>(labelRHitBH_)),
0207       tok_hgcPHEE_(consumes<edm::PassiveHitContainer>(labelPassiveEE_)),
0208       tok_hgcPHFH_(consumes<edm::PassiveHitContainer>(labelPassiveFH_)),
0209       tok_hgcPHBH_(consumes<edm::PassiveHitContainer>(labelPassiveBH_)),
0210       tok_hgcPHCMSE_(consumes<edm::PassiveHitContainer>(labelPassiveCMSE_)),
0211       tok_hgcPHBeam_(consumes<edm::PassiveHitContainer>(labelPassiveBeam_)),
0212       tokDDDEE_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0213           edm::ESInputTag("", detectorEE_))),
0214       tokGeomEE_(
0215           esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorEE_))),
0216       tokDDDFH_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0217           edm::ESInputTag("", detectorFH_))),
0218       tokGeomFH_(
0219           esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag("", detectorFH_))) {
0220   usesResource("TFileService");
0221   ahcalGeom_ = std::make_unique<AHCalGeometry>(iConfig);
0222 
0223   // now do whatever initialization is needed
0224   ///TIME SMEARING
0225   ///Const = 50ps = 0.05ns
0226   ///Stochastic is 4ns/Q(fC)
0227   ////For 300 um, 1MIP = 3.84fC, 200um, it is 3.84*2./3.=2.56fC
0228   ////This stochastic for 300um = 4*3.84/E_MIP=15.5ns/E_MIP and 200um = 4*2.56/E_MIPs=10.24ns/E_MIP
0229   idBeams_ = (iConfig.getParameter<std::vector<int>>("idBeams"));
0230 #ifdef EDM_ML_DEBUG
0231   edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: SimHits = " << doSimHits_ << " Digis = " << doDigis_ << ":"
0232                              << sampleIndex_ << " RecHits = " << doRecHits_ << " useDets " << ifEE_ << ":" << ifFH_
0233                              << ":" << ifBH_ << ":" << ifBeam_ << " zFront " << zFrontEE_ << ":" << zFrontFH_ << ":"
0234                              << zFrontBH_ << " IdBeam " << idBeams_.size() << ":";
0235   for (unsigned int k = 0; k < idBeams_.size(); ++k)
0236     edm::LogVerbatim("HGCSim") << " [" << k << "] " << idBeams_[k];
0237   edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: DoPassive " << doPassive_ << ":" << doPassiveEE_ << ":"
0238                              << doPassiveHE_ << ":" << doPassiveBH_;
0239   edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: MIP conversion factors " << gev2mip200_ << ":" << gev2mip300_
0240                              << " Time smearing " << stoc_smear_time_200_ << ":" << stoc_smear_time_300_ << " AddP "
0241                              << addP_;
0242 #endif
0243   if (idBeams_.empty())
0244     idBeams_.push_back(1001);
0245 
0246 #ifdef EDM_ML_DEBUG
0247   edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: GeneratorSource = " << labelGen_;
0248   if (ifEE_)
0249     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorEE_ << " with tags " << labelHitEE_ << ", "
0250                                << labelDigiEE_ << ", " << labelRHitEE_;
0251   if (ifFH_)
0252     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorFH_ << " with tags " << labelHitFH_ << ", "
0253                                << labelDigiFH_ << ", " << labelRHitFH_;
0254   if (ifBH_)
0255     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorBH_ << " with tags " << labelHitBH_ << ", "
0256                                << labelDigiFH_ << ", " << labelRHitBH_;
0257   if (ifBeam_)
0258     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer:: Detector " << detectorBeam_ << " with tags " << labelHitBeam_;
0259 #endif
0260 }
0261 
0262 HGCalTBAnalyzer::~HGCalTBAnalyzer() {}
0263 
0264 void HGCalTBAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0265   edm::ParameterSetDescription desc;
0266   desc.add<std::string>("detectorEE", "HGCalEESensitive");
0267   desc.add<bool>("useEE", true);
0268   desc.add<double>("zFrontEE", 0.0);
0269   desc.add<std::string>("caloHitSrcEE", "HGCHitsEE");
0270   desc.add<edm::InputTag>("digiSrcEE", edm::InputTag("hgcalDigis", "EE"));
0271   desc.add<edm::InputTag>("recHitSrcEE", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0272   desc.add<std::string>("detectorFH", "HGCalHESiliconSensitive");
0273   desc.add<bool>("useFH", false);
0274   desc.add<double>("zFrontFH", 0.0);
0275   desc.add<std::string>("caloHitSrcFH", "HGCHitsHEfront");
0276   desc.add<edm::InputTag>("digiSrcFH", edm::InputTag("hgcalDigis", "HEfront"));
0277   desc.add<edm::InputTag>("recHitSrcFH", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0278   desc.add<std::string>("detectorBH", "AHCal");
0279   desc.add<bool>("useBH", false);
0280   desc.add<double>("zFrontBH", 0.0);
0281   desc.add<std::string>("caloHitSrcBH", "HcalHits");
0282   desc.add<edm::InputTag>("digiSrcBH", edm::InputTag("hgcalDigis", "HEback"));
0283   desc.add<edm::InputTag>("recHitSrcBH", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0284   desc.add<std::string>("detectorBeam", "HcalTB06BeamDetector");
0285   desc.add<bool>("useBeam", false);
0286   desc.add<std::string>("caloHitSrcBeam", "HcalTB06BeamHits");
0287   std::vector<int> ids = {
0288       1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1011, 1012, 1013, 1014, 2001, 2002, 2003, 2004, 2005};
0289   desc.add<std::vector<int>>("idBeams", ids);
0290   desc.add<edm::InputTag>("generatorSrc", edm::InputTag("generatorSmeared"));
0291   desc.add<edm::InputTag>("passiveEE", edm::InputTag("g4SimHits", "HGCalEEPassiveHits"));
0292   desc.add<edm::InputTag>("passiveFH", edm::InputTag("g4SimHits", "HGCalHEPassiveHits"));
0293   desc.add<edm::InputTag>("passiveBH", edm::InputTag("g4SimHits", "HGCalAHPassiveHits"));
0294   desc.add<edm::InputTag>("passiveCMSE", edm::InputTag("g4SimHits", "CMSEPassiveHits"));
0295   desc.add<edm::InputTag>("passiveBeam", edm::InputTag("g4SimHits", "HGCalBeamPassiveHits"));
0296 
0297   desc.add<bool>("doSimHits", true);
0298   desc.add<bool>("doDigis", true);
0299   desc.add<bool>("doRecHits", true);
0300   desc.add<int>("sampleIndex", 0);
0301   desc.add<bool>("doTree", true);
0302   desc.add<bool>("doTreeCell", true);
0303   desc.add<bool>("doPassive", false);
0304   desc.add<bool>("doPassiveEE", false);
0305   desc.add<bool>("doPassiveHE", false);
0306   desc.add<bool>("doPassiveBH", false);
0307   desc.add<bool>("addP", false);
0308   desc.add<bool>("doBeam", false);
0309   desc.addUntracked<double>("gev2mip200", 57.0e-6);
0310   desc.addUntracked<double>("gev2mip300", 85.5e-6);
0311   desc.addUntracked<double>("stoc_smear_time_200", 10.24);
0312   desc.addUntracked<double>("stoc_smear_time_300", 15.5);
0313   desc.addUntracked<int>("maxDepth", 12);
0314   desc.addUntracked<double>("deltaX", 30.0);  // Size of tile along X
0315   desc.addUntracked<double>("deltaY", 30.0);  // Size of tile along Y
0316   desc.addUntracked<double>("deltaZ", 81.0);  // Thickness of a single layer
0317   desc.addUntracked<double>("zFirst", 17.6);  // Position of the center
0318   descriptions.add("HGCalTBAnalyzer", desc);
0319 }
0320 
0321 void HGCalTBAnalyzer::beginJob() {
0322   char name[40], title[100];
0323   hBeam_ = fs_->make<TH1D>("BeamP", "Beam Momentum", 1000, 0, 1000.0);
0324   for (int i = 0; i < 3; ++i) {
0325     bool book(ifEE_);
0326     std::string det(detectorEE_);
0327     if (i == 1) {
0328       book = ifFH_;
0329       det = detectorFH_;
0330     } else if (i == 2) {
0331       book = ifBH_;
0332       det = detectorBH_;
0333     }
0334 
0335     if (doSimHits_ && book) {
0336       sprintf(name, "SimHitEn%s", det.c_str());
0337       sprintf(title, "Sim Hit Energy for %s", det.c_str());
0338       hSimHitE_[i] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0339       sprintf(name, "SimHitEnX%s", det.c_str());
0340       sprintf(title, "Sim Hit Energy for %s", det.c_str());
0341       hSimHitEn_[i] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0342       sprintf(name, "SimHitTm%s", det.c_str());
0343       sprintf(title, "Sim Hit Timing for %s", det.c_str());
0344       hSimHitT_[i] = fs_->make<TH1D>(name, title, 5000, 0., 500.0);
0345       sprintf(name, "SimHitLat%s", det.c_str());
0346       sprintf(title, "Lateral Shower profile (Sim Hit) for %s", det.c_str());
0347       hSimHitLat_[i] = fs_->make<TProfile2D>(name, title, 100, -100., 100., 100, -100., 100.);
0348       sprintf(name, "SimHitLng%s", det.c_str());
0349       sprintf(title, "Longitudinal Shower profile (Sim Hit) for %s", det.c_str());
0350       hSimHitLng_[i] = fs_->make<TProfile>(name, title, 50, 0., 100.);
0351       sprintf(name, "SimHitLng1%s", det.c_str());
0352       sprintf(title, "Longitudinal Shower profile (Layer) for %s", det.c_str());
0353       hSimHitLng1_[i] = fs_->make<TProfile>(name, title, 200, 0., 100.);
0354       sprintf(name, "SimHitLng2%s", det.c_str());
0355       sprintf(title, "Longitudinal Shower profile (Layer) for %s", det.c_str());
0356       hSimHitLng2_[i] = fs_->make<TProfile>(name, title, 200, 0., 100.);
0357     }
0358 
0359     if (doDigis_ && book) {
0360       sprintf(name, "DigiADC%s", det.c_str());
0361       sprintf(title, "ADC at Digi level for %s", det.c_str());
0362       hDigiADC_[i] = fs_->make<TH1D>(name, title, 100, 0., 100.0);
0363       sprintf(name, "DigiOcc%s", det.c_str());
0364       sprintf(title, "Occupancy (Digi)for %s", det.c_str());
0365       hDigiOcc_[i] = fs_->make<TH2D>(name, title, 100, -10., 10., 100, -10., 10.);
0366       sprintf(name, "DigiLng%s", det.c_str());
0367       sprintf(title, "Longitudinal Shower profile (Digi) for %s", det.c_str());
0368       hDigiLng_[i] = fs_->make<TH1D>(name, title, 100, 0., 10.);
0369     }
0370 
0371     if (doRecHits_ && book) {
0372       sprintf(name, "RecHitEn%s", det.c_str());
0373       sprintf(title, "Rec Hit Energy for %s", det.c_str());
0374       hRecHitE_[i] = fs_->make<TH1D>(name, title, 1000, 0., 10.0);
0375       sprintf(name, "RecHitOcc%s", det.c_str());
0376       sprintf(title, "Occupancy (Rec Hit)for %s", det.c_str());
0377       hRecHitOcc_[i] = fs_->make<TH2D>(name, title, 100, -10., 10., 100, -10., 10.);
0378       sprintf(name, "RecHitLat%s", det.c_str());
0379       sprintf(title, "Lateral Shower profile (Rec Hit) for %s", det.c_str());
0380       hRecHitLat_[i] = fs_->make<TProfile2D>(name, title, 100, -10., 10., 100, -10., 10.);
0381       sprintf(name, "RecHitLng%s", det.c_str());
0382       sprintf(title, "Longitudinal Shower profile (Rec Hit) for %s", det.c_str());
0383       hRecHitLng_[i] = fs_->make<TProfile>(name, title, 100, 0., 10.);
0384       sprintf(name, "RecHitLng1%s", det.c_str());
0385       sprintf(title, "Longitudinal Shower profile vs Layer for %s", det.c_str());
0386       hRecHitLng1_[i] = fs_->make<TProfile>(name, title, 120, 0., 60.);
0387     }
0388   }
0389   if (ifBeam_ && doSimHits_) {
0390     sprintf(name, "SimHitEn%s", detectorBeam_.c_str());
0391     sprintf(title, "Sim Hit Energy for %s", detectorBeam_.c_str());
0392     hSimHitE_[3] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0393     sprintf(name, "SimHitEnX%s", detectorBeam_.c_str());
0394     sprintf(title, "Sim Hit Energy for %s", detectorBeam_.c_str());
0395     hSimHitEn_[3] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0396     sprintf(name, "SimHitTm%s", detectorBeam_.c_str());
0397     sprintf(title, "Sim Hit Timing for %s", detectorBeam_.c_str());
0398     hSimHitT_[3] = fs_->make<TH1D>(name, title, 5000, 0., 500.0);
0399   }
0400   if (doSimHits_ && doTree_) {
0401     tree_ = fs_->make<TTree>("HGCTB", "SimHitEnergy");
0402     tree_->Branch("simHitLayEn1EE", &simHitLayEn1EE_);
0403     tree_->Branch("simHitLayEn2EE", &simHitLayEn2EE_);
0404     tree_->Branch("simHitLayEn1FH", &simHitLayEn1FH_);
0405     tree_->Branch("simHitLayEn2FH", &simHitLayEn2FH_);
0406     tree_->Branch("simHitLayEn1BH", &simHitLayEn1BH_);
0407     tree_->Branch("simHitLayEn2BH", &simHitLayEn2BH_);
0408     tree_->Branch("xBeam", &xBeam_, "xBeam/D");
0409     tree_->Branch("yBeam", &yBeam_, "yBeam/D");
0410     tree_->Branch("zBeam", &zBeam_, "zBeam/D");
0411     tree_->Branch("pBeam", &pBeam_, "pBeam/D");
0412     tree_->Branch("thetaBeam", &thetaBeam_, "thetaBeam/D");
0413     tree_->Branch("phiBeam", &phiBeam_, "phiBeam/D");
0414     if (doBeam_) {
0415       tree_->Branch("nBeamMC", &nBeamMC_, "nBeamMC/I");
0416       tree_->Branch("pdgIdBeamMC", &pdgIdBeamMC_);
0417       tree_->Branch("xBeamMC", &xBeamMC_);
0418       tree_->Branch("yBeamMC", &yBeamMC_);
0419       tree_->Branch("zBeamMC", &zBeamMC_);
0420       tree_->Branch("pxBeamMC", &pxBeamMC_);
0421       tree_->Branch("pyBeamMC", &pyBeamMC_);
0422       tree_->Branch("pzBeamMC", &pzBeamMC_);
0423       tree_->Branch("pBeamMC", &pBeamMC_);
0424     }
0425     if (doTreeCell_) {
0426       tree_->Branch("simHitCellIdEE", &simHitCellIdEE_);
0427       tree_->Branch("simHitCellEnEE", &simHitCellEnEE_);
0428       tree_->Branch("simHitCellIdFH", &simHitCellIdFH_);
0429       tree_->Branch("simHitCellEnFH", &simHitCellEnFH_);
0430       tree_->Branch("simHitCellIdBH", &simHitCellIdBH_);
0431       tree_->Branch("simHitCellEnBH", &simHitCellEnBH_);
0432       tree_->Branch("simHitCellIdBeam", &simHitCellIdBeam_);
0433       tree_->Branch("simHitCellEnBeam", &simHitCellEnBeam_);
0434 
0435       tree_->Branch("simHitCellColBH", &simHitCellColBH_);
0436       tree_->Branch("simHitCellRowBH", &simHitCellRowBH_);
0437       tree_->Branch("simHitCellLayerBH", &simHitCellLayerBH_);
0438       tree_->Branch("simHitCellTimeFirstHitEE", &simHitCellTimeFirstHitEE_);
0439       tree_->Branch("simHitCellTimeFirstHitFH", &simHitCellTimeFirstHitFH_);
0440       //tree_->Branch("simHitCellTimeFirstHitBH", &simHitCellTimeFirstHitBH_);
0441       tree_->Branch("simHitCellTime15MipEE", &simHitCellTime15MipEE_);
0442       tree_->Branch("simHitCellTime15MipFH", &simHitCellTime15MipFH_);
0443       //tree_->Branch("simHitCellTime15MipBH", &simHitCellTime15MipBH_);
0444       tree_->Branch("simHitCellTimeLastHitEE", &simHitCellTimeLastHitEE_);
0445       tree_->Branch("simHitCellTimeLastHitFH", &simHitCellTimeLastHitFH_);
0446       //tree_->Branch("simHitCellTimeLastHitBH", &simHitCellTimeLastHitBH_);
0447     }
0448   }
0449 
0450   if (doPassive_ && doTree_) {
0451     if (doPassiveEE_) {
0452       tree_->Branch("hgcPassiveEEEnergy", &hgcPassiveEEEnergy_);
0453       tree_->Branch("hgcPassiveEEName", &hgcPassiveEEName_);
0454       tree_->Branch("hgcPassiveEEID", &hgcPassiveEEID_);
0455     }
0456     if (doPassiveHE_) {
0457       tree_->Branch("hgcPassiveFHEnergy", &hgcPassiveFHEnergy_);
0458       tree_->Branch("hgcPassiveFHName", &hgcPassiveFHName_);
0459       tree_->Branch("hgcPassiveFHID", &hgcPassiveFHID_);
0460     }
0461     if (doPassiveBH_) {
0462       tree_->Branch("hgcPassiveBHEnergy", &hgcPassiveBHEnergy_);
0463       tree_->Branch("hgcPassiveBHName", &hgcPassiveBHName_);
0464       tree_->Branch("hgcPassiveBHID", &hgcPassiveBHID_);
0465     }
0466     tree_->Branch("hgcPassiveCMSEEnergy", &hgcPassiveCMSEEnergy_);
0467     tree_->Branch("hgcPassiveCMSEName", &hgcPassiveCMSEName_);
0468     tree_->Branch("hgcPassiveCMSEID", &hgcPassiveCMSEID_);
0469     tree_->Branch("hgcPassiveBeamEnergy", &hgcPassiveBeamEnergy_);
0470     tree_->Branch("hgcPassiveBeamName", &hgcPassiveBeamName_);
0471     tree_->Branch("hgcPassiveBeamID", &hgcPassiveBeamID_);
0472   }
0473 }
0474 
0475 void HGCalTBAnalyzer::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0476   char name[40], title[100];
0477   if (ifEE_) {
0478     hgcons_[0] = &iSetup.getData(tokDDDEE_);
0479     if (doDigis_ || doRecHits_) {
0480       hgeom_[0] = &iSetup.getData(tokGeomEE_);
0481     } else {
0482       hgeom_[0] = nullptr;
0483     }
0484     for (unsigned int l = 0; l < hgcons_[0]->layers(false); ++l) {
0485       sprintf(name, "SimHitEnA%d%s", l, detectorEE_.c_str());
0486       sprintf(title, "Sim Hit Energy in SIM layer %d for %s", l + 1, detectorEE_.c_str());
0487       hSimHitLayEn1EE_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0488       if (l % 3 == 0) {
0489         sprintf(name, "SimHitEnB%d%s", (l / 3 + 1), detectorEE_.c_str());
0490         sprintf(title, "Sim Hit Energy in layer %d for %s", (l / 3 + 1), detectorEE_.c_str());
0491         hSimHitLayEn2EE_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0492       }
0493     }
0494 #ifdef EDM_ML_DEBUG
0495     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::" << detectorEE_ << " defined with " << hgcons_[0]->layers(false)
0496                                << " layers";
0497 #endif
0498   } else {
0499     hgcons_[0] = nullptr;
0500     hgeom_[0] = nullptr;
0501   }
0502 
0503   if (ifFH_) {
0504     hgcons_[1] = &iSetup.getData(tokDDDFH_);
0505     if (doDigis_ || doRecHits_) {
0506       hgeom_[1] = &iSetup.getData(tokGeomFH_);
0507     } else {
0508       hgeom_[1] = nullptr;
0509     }
0510     for (unsigned int l = 0; l < hgcons_[1]->layers(false); ++l) {
0511       sprintf(name, "SimHitEnA%d%s", l, detectorFH_.c_str());
0512       sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorFH_.c_str());
0513       hSimHitLayEn1FH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0514       if (l % 3 == 0) {
0515         sprintf(name, "SimHitEnB%d%s", (l / 3 + 1), detectorFH_.c_str());
0516         sprintf(title, "Sim Hit Energy in layer %d for %s", (l / 3 + 1), detectorFH_.c_str());
0517         hSimHitLayEn2FH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0518       }
0519     }
0520 #ifdef EDM_ML_DEBUG
0521     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::" << detectorFH_ << " defined with " << hgcons_[1]->layers(false)
0522                                << " layers";
0523 #endif
0524   } else {
0525     hgcons_[1] = nullptr;
0526     hgeom_[1] = nullptr;
0527   }
0528 
0529   if (ifBH_) {
0530     for (int l = 0; l < ahcalGeom_->maxDepth(); ++l) {
0531       sprintf(name, "SimHitEnA%d%s", l, detectorBH_.c_str());
0532       sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorBH_.c_str());
0533       hSimHitLayEn1BH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0534       sprintf(name, "SimHitEnB%d%s", l, detectorBH_.c_str());
0535       sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorBH_.c_str());
0536       hSimHitLayEn2BH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0537     }
0538   }
0539 
0540   if (ifBeam_) {
0541     for (unsigned int l = 0; l < idBeams_.size(); ++l) {
0542       sprintf(name, "SimHitEna%d%s", l, detectorBeam_.c_str());
0543       sprintf(title, "Sim Hit Energy in type %d for %s", idBeams_[l], detectorBeam_.c_str());
0544       hSimHitLayEnBeam_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0545     }
0546   }
0547 }
0548 
0549 void HGCalTBAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0550   // Generator input
0551   const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_hepMC_);
0552   if (!evtMC.isValid()) {
0553     edm::LogWarning("HGCal") << "no HepMCProduct found";
0554   } else {
0555     const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0556     unsigned int k(0);
0557     HepMC::FourVector pxyz(0, 0, 0, 0);
0558     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0559          ++p, ++k) {
0560       edm::LogVerbatim("HGCSim") << "Particle [" << k << "] with p " << (*p)->momentum().rho() << " theta "
0561                                  << (*p)->momentum().theta() << " phi " << (*p)->momentum().phi() << " pxyz ("
0562                                  << (*p)->momentum().px() << ", " << (*p)->momentum().py() << ", "
0563                                  << (*p)->momentum().pz() << ")";
0564       if (addP_) {
0565         pxyz.setPx(pxyz.px() + (*p)->momentum().px());
0566         pxyz.setPy(pxyz.py() + (*p)->momentum().py());
0567         pxyz.setPz(pxyz.pz() + (*p)->momentum().pz());
0568         pxyz.setE(pxyz.e() + (*p)->momentum().e());
0569       } else if (!addP_ && (k == 0)) {
0570         pxyz = (*p)->momentum();
0571       }
0572     }
0573     hBeam_->Fill(pxyz.rho());
0574     edm::LogVerbatim("HGCSim") << "Particle with p " << pxyz.rho() << " theta " << pxyz.theta() << " phi "
0575                                << pxyz.phi();
0576   }
0577 
0578   // Now the Simhits
0579   if (doSimHits_) {
0580     const edm::Handle<edm::SimTrackContainer>& SimTk = iEvent.getHandle(tok_simTk_);
0581     const edm::Handle<edm::SimVertexContainer>& SimVtx = iEvent.getHandle(tok_simVtx_);
0582     analyzeSimTracks(SimTk, SimVtx);
0583 
0584     simHitLayEn1EE_.clear();
0585     simHitLayEn2EE_.clear();
0586     simHitLayEn1FH_.clear();
0587     simHitLayEn2FH_.clear();
0588     simHitLayEn1BH_.clear();
0589     simHitLayEn2BH_.clear();
0590     simHitLayEnBeam_.clear();
0591     simHitCellIdEE_.clear();
0592     simHitCellEnEE_.clear();
0593     simHitCellIdFH_.clear();
0594     simHitCellEnFH_.clear();
0595     simHitCellIdBH_.clear();
0596     simHitCellEnBH_.clear();
0597     simHitCellIdBeam_.clear();
0598     simHitCellEnBeam_.clear();
0599     simHitCellColBH_.clear();
0600     simHitCellRowBH_.clear();
0601     simHitCellLayerBH_.clear();
0602     simHitCellTimeFirstHitEE_.clear();
0603     simHitCellTime15MipEE_.clear();
0604     simHitCellTimeLastHitEE_.clear();
0605     simHitCellTimeFirstHitFH_.clear();
0606     simHitCellTime15MipFH_.clear();
0607     simHitCellTimeLastHitFH_.clear();
0608     std::vector<PCaloHit> caloHits;
0609     if (ifEE_) {
0610       simHitLayEn1EE_ = std::vector<float>(hgcons_[0]->layers(false), 0);
0611       simHitLayEn2EE_ = std::vector<float>(hgcons_[0]->layers(true), 0);
0612       const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsEE_);
0613       if (theCaloHitContainers.isValid()) {
0614 #ifdef EDM_ML_DEBUG
0615         edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorEE_ << " has " << theCaloHitContainers->size()
0616                                    << " hits";
0617 #endif
0618         caloHits.clear();
0619         caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0620         analyzeSimHits(0, caloHits, zFrontEE_);
0621       } else {
0622 #ifdef EDM_ML_DEBUG
0623         edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorEE_ << " !!!";
0624 #endif
0625       }
0626     }
0627     if (ifFH_) {
0628       simHitLayEn1FH_ = std::vector<float>(hgcons_[1]->layers(false), 0);
0629       simHitLayEn2FH_ = std::vector<float>(hgcons_[1]->layers(true), 0);
0630       const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsFH_);
0631       if (theCaloHitContainers.isValid()) {
0632 #ifdef EDM_ML_DEBUG
0633         edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorFH_ << " has " << theCaloHitContainers->size()
0634                                    << " hits";
0635 #endif
0636         caloHits.clear();
0637         caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0638         analyzeSimHits(1, caloHits, zFrontFH_);
0639       } else {
0640 #ifdef EDM_ML_DEBUG
0641         edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorFH_ << " !!!";
0642 #endif
0643       }
0644     }
0645     if (ifBH_) {
0646       simHitLayEn1BH_ = std::vector<float>(ahcalGeom_->maxDepth(), 0);
0647       simHitLayEn2BH_ = std::vector<float>(ahcalGeom_->maxDepth(), 0);
0648       const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsBH_);
0649       if (theCaloHitContainers.isValid()) {
0650 #ifdef EDM_ML_DEBUG
0651         edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBH_ << " has " << theCaloHitContainers->size()
0652                                    << " hits";
0653 #endif
0654         caloHits.clear();
0655         caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0656         analyzeSimHits(2, caloHits, zFrontBH_);
0657       } else {
0658 #ifdef EDM_ML_DEBUG
0659         edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBH_ << " !!!";
0660 #endif
0661       }
0662     }
0663     if (ifBeam_) {
0664       simHitLayEnBeam_ = std::vector<float>(idBeams_.size(), 0);
0665       const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsBeam_);
0666       if (theCaloHitContainers.isValid()) {
0667 #ifdef EDM_ML_DEBUG
0668         edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBeam_ << " has "
0669                                    << theCaloHitContainers->size() << " hits";
0670 #endif
0671         caloHits.clear();
0672         caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0673         analyzeSimHits(3, caloHits, 0.0);
0674       } else {
0675 #ifdef EDM_ML_DEBUG
0676         edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBeam_ << " !!!";
0677 #endif
0678       }
0679     }
0680   }  // if (doSimHits_)
0681 
0682   ////Store the info about the Passive hits
0683   if (doPassive_) {
0684     /// EE
0685     hgcPassiveEEEnergy_.clear();
0686     hgcPassiveEEName_.clear();
0687     hgcPassiveEEID_.clear();
0688     if (doPassiveEE_) {
0689       const edm::Handle<edm::PassiveHitContainer>& hgcPHEE = iEvent.getHandle(tok_hgcPHEE_);
0690       analyzePassiveHits(hgcPHEE, 1);
0691     }
0692     /// FH
0693     hgcPassiveFHEnergy_.clear();
0694     hgcPassiveFHName_.clear();
0695     hgcPassiveFHID_.clear();
0696     if (doPassiveHE_) {
0697       const edm::Handle<edm::PassiveHitContainer>& hgcPHFH = iEvent.getHandle(tok_hgcPHFH_);
0698       analyzePassiveHits(hgcPHFH, 2);
0699     }
0700     /// BH
0701     hgcPassiveBHEnergy_.clear();
0702     hgcPassiveBHName_.clear();
0703     hgcPassiveBHID_.clear();
0704     if (doPassiveBH_) {
0705       const edm::Handle<edm::PassiveHitContainer>& hgcPHBH = iEvent.getHandle(tok_hgcPHBH_);
0706       analyzePassiveHits(hgcPHBH, 3);
0707     }
0708     /// CMSE
0709     hgcPassiveCMSEEnergy_.clear();
0710     hgcPassiveCMSEName_.clear();
0711     hgcPassiveCMSEID_.clear();
0712     const edm::Handle<edm::PassiveHitContainer>& hgcPHCMSE = iEvent.getHandle(tok_hgcPHCMSE_);
0713     analyzePassiveHits(hgcPHCMSE, 4);
0714     ///Beam
0715     hgcPassiveBeamName_.clear();
0716     hgcPassiveBeamEnergy_.clear();
0717     hgcPassiveBeamID_.clear();
0718     const edm::Handle<edm::PassiveHitContainer>& hgcPHBeam = iEvent.getHandle(tok_hgcPHBeam_);
0719     analyzePassiveHits(hgcPHBeam, 5);
0720   }
0721 
0722   if ((doSimHits_ || doPassive_) && (doTree_))
0723     tree_->Fill();
0724 
0725   // Now the Digis
0726   if (doDigis_) {
0727     if (ifEE_) {
0728       const edm::Handle<HGCalDigiCollection>& theDigiContainers = iEvent.getHandle(tok_digiEE_);
0729       if (theDigiContainers.isValid()) {
0730 #ifdef EDM_ML_DEBUG
0731         edm::LogVerbatim("HGCSim") << "HGCDigiCintainer for " << detectorEE_ << " with " << theDigiContainers->size()
0732                                    << " element(s)";
0733 #endif
0734         for (const auto& it : *theDigiContainers) {
0735           HGCalDetId detId = (it.id());
0736           const HGCSample& hgcSample = it.sample(sampleIndex_);
0737           uint16_t adc = hgcSample.data();
0738           analyzeDigi(0, detId, adc);
0739         }
0740       }
0741     }
0742     if (ifFH_) {
0743       const edm::Handle<HGCalDigiCollection>& theDigiContainers = iEvent.getHandle(tok_digiFH_);
0744       if (theDigiContainers.isValid()) {
0745 #ifdef EDM_ML_DEBUG
0746         edm::LogVerbatim("HGCSim") << "HGCDigiContainer for " << detectorFH_ << " with " << theDigiContainers->size()
0747                                    << " element(s)";
0748 #endif
0749         for (const auto& it : *theDigiContainers) {
0750           HGCalDetId detId = (it.id());
0751           const HGCSample& hgcSample = it.sample(sampleIndex_);
0752           uint16_t adc = hgcSample.data();
0753           analyzeDigi(1, detId, adc);
0754         }
0755       }
0756     }
0757   }
0758 
0759   // The Rechits
0760   if (doRecHits_) {
0761     if (ifEE_) {
0762       const edm::Handle<HGCRecHitCollection>& theCaloHitContainers = iEvent.getHandle(tok_hitrEE_);
0763       if (theCaloHitContainers.isValid()) {
0764 #ifdef EDM_ML_DEBUG
0765         edm::LogVerbatim("HGCSim") << "HGCRecHitCollection for " << detectorEE_ << " has "
0766                                    << theCaloHitContainers->size() << " hits";
0767 #endif
0768         analyzeRecHits(0, theCaloHitContainers);
0769       } else {
0770 #ifdef EDM_ML_DEBUG
0771         edm::LogVerbatim("HGCSim") << "HGCRecHitCollection does not exist for " << detectorEE_ << " !!!";
0772 #endif
0773       }
0774     }
0775     if (ifFH_) {
0776       const edm::Handle<HGCRecHitCollection>& theCaloHitContainers = iEvent.getHandle(tok_hitrFH_);
0777       if (theCaloHitContainers.isValid()) {
0778 #ifdef EDM_ML_DEBUG
0779         edm::LogVerbatim("HGCSim") << "HGCRecHitCollection for " << detectorFH_ << " has "
0780                                    << theCaloHitContainers->size() << " hits";
0781 #endif
0782         analyzeRecHits(1, theCaloHitContainers);
0783       } else {
0784 #ifdef EDM_ML_DEBUG
0785         edm::LogVerbatim("HGCSim") << "HGCRecHitCollection does not exist for " << detectorFH_ << " !!!";
0786 #endif
0787       }  // else
0788     }    // if (ifFH_)
0789   }      // if (doRecHits_)
0790 
0791 }  // void HGCalTBAnalyzer::analyze
0792 
0793 void HGCalTBAnalyzer::analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront) {
0794   std::map<uint32_t, double> map_hits, map_hitn;
0795   std::map<uint32_t, double> map_hittime_firsthit, map_hittime_lasthit, map_hittime_15Mip;
0796   std::map<int, double> map_hitDepth, map_hitWafer;
0797   std::map<int, std::pair<uint32_t, double>> map_hitLayer, map_hitCell;
0798   double entot(0);
0799   std::map<uint32_t, double> nhits;
0800   std::map<uint32_t, int> ID, Depth;
0801   std::map<uint32_t, double> GeV2Mip;
0802   std::map<uint32_t, double> StochTermTime;
0803   std::vector<int> nSimLayers;
0804   ///storing time of each pcalohit in each cell ///SJ
0805   std::map<uint32_t, std::vector<std::pair<double, double>>> map_hitTimeEn;
0806   //bool debug = true;
0807   bool debug = false;
0808   for (unsigned int i = 0; i < hits.size(); i++) {
0809     double energy = hits[i].energy();
0810     double time = hits[i].time();
0811     uint32_t id = hits[i].id();
0812     entot += energy;
0813     int subdet, zside, layer, sector, subsector(0), cell, depth(0), idx(0);
0814     if (type == 2) {
0815       subdet = HcalDetId(id).subdet();
0816       if (subdet != HcalOther)
0817         continue;
0818       AHCalDetId hid(id);
0819       layer = depth = hid.depth();
0820       zside = hid.zside();
0821       sector = hid.irow();
0822       cell = hid.icol();
0823       idx = ((hid.irowAbs() * 100) + (hid.icolAbs()));
0824       if (debug)
0825         edm::LogVerbatim("HGCSim") << "depth, sector, cell " << depth << ":" << sector << ":" << cell;
0826     } else if (type == 3) {
0827       HcalTestBeamNumbering::unpackIndex(id, subdet, layer, sector, cell);
0828       depth = layer;
0829       zside = 1;
0830       idx = subdet * 1000 + layer;
0831       layer = idx;
0832     } else {
0833       HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, subsector, cell);
0834       depth = hgcons_[type]->simToReco(cell, layer, sector, true).second;
0835       idx = sector * 1000 + cell;
0836 #ifdef EDM_ML_DEBUG
0837       std::pair<float, float> xy = hgcons_[type]->locateCell(cell, layer, sector, false);
0838       edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::detId " << std::hex << id << std::dec << " Layer:Wafer:Cell "
0839                                  << layer << ":" << sector << ":" << cell << " Position " << xy.first << ":"
0840                                  << xy.second << ":" << hgcons_[type]->waferZ(layer, false);
0841 #endif
0842     }
0843 #ifdef EDM_ML_DEBUG
0844     edm::LogVerbatim("HGCSim") << "SimHit:Hit[" << i << "] Id " << subdet << ":" << zside << ":" << layer << ":"
0845                                << sector << ":" << subsector << ":" << cell << ":" << depth << " Energy " << energy
0846                                << " Time " << time;
0847 #endif
0848     if (map_hits.count(id) != 0) {
0849       map_hits[id] += energy;
0850     } else {
0851       map_hits[id] = energy;
0852     }
0853     if (map_hitLayer.count(layer) != 0) {
0854       double ee = energy + map_hitLayer[layer].second;
0855       map_hitLayer[layer] = std::make_pair(id, ee);
0856     } else {
0857       map_hitLayer[layer] = std::make_pair(id, energy);
0858     }
0859     if (map_hitWafer.count(sector) != 0)
0860       map_hitWafer[sector] += energy;
0861     else
0862       map_hitWafer[sector] = energy;
0863     if (depth >= 0) {
0864       if (map_hitCell.count(idx) != 0) {
0865         double ee = energy + map_hitCell[idx].second;
0866         map_hitCell[idx] = std::make_pair(id, ee);
0867       } else {
0868         map_hitCell[idx] = std::make_pair(id, energy);
0869       }
0870       if (debug) {
0871         if (type == 1)
0872           edm::LogVerbatim("HGCSim") << "EE, depth is and map_hitDepth[depth] " << depth << " " << map_hitDepth[depth];
0873         if (type == 2)
0874           edm::LogVerbatim("HGCSim") << "BH, depth is and map_hitDepth[depth] " << depth << " " << map_hitDepth[depth];
0875       }
0876       if (map_hitDepth.count(depth) != 0) {
0877         map_hitDepth[depth] += energy;
0878       } else {
0879         map_hitDepth[depth] = energy;
0880       }
0881       uint32_t idn =
0882           (type >= 2) ? id : HGCalTestNumbering::packHexagonIndex(subdet, zside, depth, sector, subsector, cell);
0883       map_hitTimeEn[idn].push_back(std::make_pair(time, energy));
0884       GeV2Mip[idn] = gev2mip300_;
0885       StochTermTime[idn] = stoc_smear_time_300_;
0886       ID[idn] = id;
0887       Depth[idn] = depth;
0888       if (map_hitn.count(idn) != 0) {
0889         map_hitn[idn] += energy;
0890         ++nhits[idn];
0891       } else {
0892         map_hitn[idn] = energy;
0893         nhits[idn] = 1;
0894       }
0895     }
0896     hSimHitT_[type]->Fill(time, energy);
0897   }
0898 
0899   if (type < 2) {  //store only for EE and FH
0900     edm::LogVerbatim("HGCSim") << "HGCalTAnalyzer:: " << map_hitWafer.size() << " wafers are hit in type " << type;
0901     for (auto itr = map_hitWafer.begin(); itr != map_hitWafer.end(); ++itr)
0902       edm::LogVerbatim("HGCSim") << "Wafer: " << itr->first << " Deposited Energy " << itr->second;
0903     ///now sort the vector of each cell hits
0904     for (const auto& itr : map_hitTimeEn) {
0905       uint32_t id = itr.first;
0906       int wafer = -99;
0907       ///id: reco and ID[id]: sim ID
0908       wafer = HGCalDetId(ID[id]).wafer();
0909       double layer = HGCalDetId(id).layer();
0910       double thickness = hgcons_[type]->cellThickness(layer, wafer, 0);
0911       if (debug)
0912         edm::LogVerbatim("HGCSim") << "wafer is : depth (reco) " << wafer << " " << Depth[id]
0913                                    << "\ntype : layer : wafer thickness " << type << " " << layer << " " << thickness
0914                                    << "\nID(sim) and id(reco) " << std::hex << ID[id] << " " << id << std::dec;
0915       if (thickness == 300) {
0916         GeV2Mip[id] = gev2mip300_;
0917         StochTermTime[id] = stoc_smear_time_300_;
0918       } else if (thickness == 200) {
0919         GeV2Mip[id] = gev2mip200_;
0920         StochTermTime[id] = stoc_smear_time_200_;
0921       }
0922       //first sort
0923       std::sort(map_hitTimeEn[id].begin(), map_hitTimeEn[id].end(), sortTime);
0924       ///once it is sorted now start adding the energy
0925       double threshold = 15.;
0926       double totE = 0.;
0927       double totEbeforeThreshold = 0.;
0928       double timebeforeThreshold = 0.;
0929       double timeAtThresohld = 0.;
0930       for (unsigned int ihit = 0; ihit < map_hitTimeEn[id].size(); ihit++) {
0931         double energy = (map_hitTimeEn[id].at(ihit)).second / GeV2Mip[id];
0932         totE += energy;
0933         double time = (map_hitTimeEn[id].at(ihit)).first;
0934         if (debug)
0935           edm::LogVerbatim("HGCSim") << "Tot E till now : time of that E : GeV2Mip[id] is " << totE << " " << time
0936                                      << " " << GeV2Mip[id];
0937         if (totE < threshold) {
0938           totEbeforeThreshold = totE;
0939           timebeforeThreshold = time;
0940         } else {
0941           timeAtThresohld =
0942               (threshold - totEbeforeThreshold) * (time - timebeforeThreshold) / (totE - totEbeforeThreshold) +
0943               timebeforeThreshold;
0944           map_hittime_15Mip[id] = timeAtThresohld;
0945           if (debug)
0946             edm::LogVerbatim("HGCSim") << "ihit : energyBefore : timeBefore : energyTot : timeTot : timeAt15MIP  "
0947                                        << ihit << " " << totEbeforeThreshold << " " << timebeforeThreshold << " "
0948                                        << totE << " " << time << " " << map_hittime_15Mip[id];
0949           break;
0950         }
0951       }
0952       if (!map_hitTimeEn[id].empty()) {
0953         map_hittime_firsthit[id] = (map_hitTimeEn[id].at(0)).first;
0954         map_hittime_lasthit[id] = (map_hitTimeEn[id].at(map_hitTimeEn[id].size() - 1)).first;
0955         if (map_hittime_15Mip[id] < map_hittime_firsthit[id])
0956           map_hittime_15Mip[id] = map_hittime_firsthit[id];
0957         /*
0958         ///Put the smeared values here after correcting for the vertex time
0959         double firsthit_sm = ran3.Gaus(map_hittime_firsthit_vtxCorr[id], StochTermTime[id]);
0960         double lasthit_sm = ran3.Gaus(map_hittime_lasthit_vtxCorr[id], StochTermTime[id]);
0961         double threshmiphit_sm = ran3.Gaus(map_hittime_15Mip_vtxCorr[id], StochTermTime[id]);
0962         */
0963       }
0964 
0965       if (totE < threshold)
0966         map_hittime_15Mip[id] = -99;
0967       if (debug)
0968         edm::LogVerbatim("HGCSim") << "id : first hit time :  last hit time " << id << " " << map_hittime_firsthit[id]
0969                                    << " " << map_hittime_lasthit[id] << "\nFinally for this cell, time is "
0970                                    << map_hittime_15Mip[id];
0971     }
0972   }
0973 
0974   hSimHitEn_[type]->Fill(entot);
0975   for (const auto& itr : map_hits) {
0976     hSimHitE_[type]->Fill(itr.second);
0977   }
0978 
0979   if (debug)
0980     edm::LogVerbatim("HGCSim") << "Now looping over map_hitLayer";
0981   for (const auto& itr : map_hitLayer) {
0982     int layer = (type == 2) ? itr.first : (itr.first - 1);
0983     double energy = (itr.second).second;
0984     double zp(0);
0985     if (type < 2)
0986       zp = hgcons_[type]->waferZ(layer + 1, false);
0987     else if (type == 2)
0988       zp = ahcalGeom_->getZ(AHCalDetId((itr.second).first));
0989 #ifdef EDM_ML_DEBUG
0990     edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer + 1 << " Z " << zp << ":" << zp - zFront << " E " << energy;
0991 #endif
0992     if (type < 3) {
0993       hSimHitLng_[type]->Fill(zp - zFront, energy);
0994       hSimHitLng2_[type]->Fill(layer + 1, energy);
0995     }
0996     if (type == 0) {
0997       if (layer < static_cast<int>(hSimHitLayEn1EE_.size())) {
0998         simHitLayEn1EE_[layer] = energy;
0999         hSimHitLayEn1EE_[layer]->Fill(energy);
1000       }
1001     } else if (type == 1) {
1002       if (layer < static_cast<int>(hSimHitLayEn1FH_.size())) {
1003         simHitLayEn1FH_[layer] = energy;
1004         hSimHitLayEn1FH_[layer]->Fill(energy);
1005       }
1006     } else if (type == 2) {
1007       if (debug)
1008         edm::LogVerbatim("HGCSim") << "layer < hSimHitLayEn1BH_.size()";
1009       if (layer < static_cast<int>(hSimHitLayEn1BH_.size())) {
1010         simHitLayEn1BH_[layer] = energy;
1011         hSimHitLayEn1BH_[layer]->Fill(energy);
1012       }
1013     } else {
1014       for (unsigned int k = 0; k < idBeams_.size(); ++k) {
1015         if (layer + 1 == idBeams_[k]) {
1016           simHitLayEnBeam_[k] = energy;
1017           hSimHitLayEnBeam_[k]->Fill(energy);
1018           break;
1019         }
1020       }
1021     }
1022   }
1023   for (const auto& itr : map_hitDepth) {
1024     int layer = (type == 2) ? itr.first : (itr.first - 1);
1025     double energy = itr.second;
1026 #ifdef EDM_ML_DEBUG
1027     edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer + 1 << " " << energy;
1028 #endif
1029     hSimHitLng1_[type]->Fill(layer + 1, energy);
1030     if (type == 0) {
1031       if (layer < static_cast<int>(hSimHitLayEn2EE_.size())) {
1032         simHitLayEn2EE_[layer] = energy;
1033         hSimHitLayEn2EE_[layer]->Fill(energy);
1034       }
1035     } else if (type == 1) {
1036       if (layer < static_cast<int>(hSimHitLayEn2FH_.size())) {
1037         simHitLayEn2FH_[layer] = energy;
1038         hSimHitLayEn2FH_[layer]->Fill(energy);
1039       }
1040     } else if (type == 2) {
1041       if (debug)
1042         edm::LogVerbatim("HGCSim") << "Inside map_hitDepth, layer no. " << layer;
1043       if (layer < static_cast<int>(hSimHitLayEn2BH_.size())) {
1044         simHitLayEn2BH_[layer] = energy;
1045         hSimHitLayEn2BH_[layer]->Fill(energy);
1046       }
1047     }
1048   }
1049 
1050   if (debug)
1051     edm::LogVerbatim("HGCSim") << "Now looping over map_hitCell";
1052   if (type < 3) {
1053     for (const auto& itr : map_hitCell) {
1054       uint32_t id = ((itr.second).first);
1055       double energy = ((itr.second).second);
1056       std::pair<float, float> xy(0, 0);
1057       double xx(0);
1058       if (type == 2) {
1059         xy = ahcalGeom_->getXY(AHCalDetId(id));
1060         xx = xy.first;
1061       } else {
1062         int subdet, zside, layer, sector, subsector, cell;
1063         HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, subsector, cell);
1064         xy = hgcons_[type]->locateCell(cell, layer, sector, false);
1065         double zp = hgcons_[type]->waferZ(layer, false);
1066         xx = (zp < 0) ? -xy.first : xy.first;
1067       }
1068       hSimHitLat_[type]->Fill(xx, xy.second, energy);
1069     }
1070   }
1071 
1072   for (const auto& itr : map_hitn) {
1073     uint32_t id = itr.first;
1074     double energy = itr.second;
1075     if (type == 0) {
1076       double time_firsthit = map_hittime_firsthit[id];
1077       double time15Mip = map_hittime_15Mip[id];
1078       double time_lasthit = map_hittime_lasthit[id];
1079       simHitCellIdEE_.push_back(id);
1080       simHitCellEnEE_.push_back(energy);
1081       simHitCellTimeFirstHitEE_.push_back(time_firsthit);
1082       simHitCellTime15MipEE_.push_back(time15Mip);
1083       simHitCellTimeLastHitEE_.push_back(time_lasthit);
1084       if (debug && (energy / GeV2Mip[id] < 15) && (map_hittime_15Mip[id] > 0))
1085         edm::LogVerbatim("HGCSim") << "FOUND!!!!rechit energy : Finally for this cell, time is " << energy / GeV2Mip[id]
1086                                    << " " << map_hittime_15Mip[id];
1087     } else if (type == 1) {
1088       double time_firsthit = map_hittime_firsthit[id];
1089       double time15Mip = map_hittime_15Mip[id];
1090       double time_lasthit = map_hittime_lasthit[id];
1091       simHitCellIdFH_.push_back(id);
1092       simHitCellEnFH_.push_back(energy);
1093       simHitCellTimeFirstHitFH_.push_back(time_firsthit);
1094       simHitCellTime15MipFH_.push_back(time15Mip);
1095       simHitCellTimeLastHitFH_.push_back(time_lasthit);
1096     } else if (type == 2) {
1097       simHitCellIdBH_.push_back(id);
1098       simHitCellEnBH_.push_back(energy);
1099       AHCalDetId hid(id);
1100       int row = hid.irow();
1101       int col = hid.icol();
1102       int layer = hid.depth();
1103       simHitCellColBH_.push_back(col);
1104       simHitCellRowBH_.push_back(row);
1105       simHitCellLayerBH_.push_back(layer);
1106 #ifdef EDM_ML_DEBUG
1107       edm::LogVerbatim("HGCSim") << "ID: " << std::hex << id << std::dec << " Layer: " << layer << " col: " << col
1108                                  << " row: " << row;
1109 #endif
1110     } else if (type == 3) {
1111       simHitCellIdBeam_.push_back(id);
1112       simHitCellEnBeam_.push_back(energy);
1113     }
1114   }
1115 }
1116 
1117 void HGCalTBAnalyzer::analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
1118                                        edm::Handle<edm::SimVertexContainer> const& SimVtx) {
1119   xBeam_ = yBeam_ = zBeam_ = pBeam_ = -9999;
1120   nBeamMC_ = thetaBeam_ = phiBeam_ = -9999;
1121   int nParBeam = 0;
1122   int vertIndex(-1);
1123   if (doBeam_) {
1124     pdgIdBeamMC_.clear();
1125     xBeamMC_.clear();
1126     yBeamMC_.clear();
1127     zBeamMC_.clear();
1128     pxBeamMC_.clear();
1129     pyBeamMC_.clear();
1130     pzBeamMC_.clear();
1131     pBeamMC_.clear();
1132   }
1133   std::vector<float> verX, verY, verZ;
1134   verX.clear();
1135   verY.clear();
1136   verZ.clear();
1137   for (const auto& simVtxItr : *SimVtx) {
1138     verX.push_back(simVtxItr.position().X());
1139     verY.push_back(simVtxItr.position().Y());
1140     verZ.push_back(simVtxItr.position().Z());
1141   }
1142 #ifdef EDM_ML_DEBUG
1143   edm::LogVerbatim("HGCSim") << "Size of track " << SimTk->size();
1144 #endif
1145   HepMC::FourVector pxyz(0, 0, 0, 0);
1146   for (const auto& simTrkItr : *SimTk) {
1147     if (addP_ && !(simTrkItr.noGenpart())) {
1148       pxyz.setPx(pxyz.px() + simTrkItr.momentum().px());
1149       pxyz.setPy(pxyz.py() + simTrkItr.momentum().py());
1150       pxyz.setPz(pxyz.pz() + simTrkItr.momentum().pz());
1151       pxyz.setE(pxyz.e() + simTrkItr.momentum().e());
1152 #ifdef EDM_ML_DEBUG
1153       edm::LogVerbatim("HGCSim") << "Track " << simTrkItr.trackId() << " Vertex " << simTrkItr.vertIndex() << " Type "
1154                                  << simTrkItr.type() << " Charge " << simTrkItr.charge() << " px "
1155                                  << simTrkItr.momentum().px() << " py " << simTrkItr.momentum().py() << " pz "
1156                                  << simTrkItr.momentum().pz() << " P " << simTrkItr.momentum().P() << " GenIndex "
1157                                  << simTrkItr.genpartIndex();
1158       edm::LogVerbatim("HGCSim") << "Vertex " << simTrkItr.vertIndex()
1159                                  << " position-> X: " << verX[simTrkItr.vertIndex()]
1160                                  << " Y: " << verY[simTrkItr.vertIndex()] << " Z: " << verZ[simTrkItr.vertIndex()];
1161 #endif
1162     }
1163     if (doBeam_ && !(simTrkItr.noGenpart())) {
1164       nParBeam++;
1165       pdgIdBeamMC_.push_back(simTrkItr.type());
1166       xBeamMC_.push_back(verX[simTrkItr.vertIndex()]);
1167       yBeamMC_.push_back(verY[simTrkItr.vertIndex()]);
1168       zBeamMC_.push_back(verZ[simTrkItr.vertIndex()]);
1169       pxBeamMC_.push_back(simTrkItr.momentum().px());
1170       pyBeamMC_.push_back(simTrkItr.momentum().py());
1171       pzBeamMC_.push_back(simTrkItr.momentum().pz());
1172       pBeamMC_.push_back(simTrkItr.momentum().P());
1173     } else if (!addP_ && (vertIndex == -1)) {
1174       pxyz = simTrkItr.momentum();
1175     }
1176     if (vertIndex == -1)
1177       vertIndex = simTrkItr.vertIndex();
1178   }
1179   nBeamMC_ = nParBeam;
1180   pBeam_ = pxyz.rho();
1181   thetaBeam_ = pxyz.theta();
1182   phiBeam_ = pxyz.phi();
1183   if (phiBeam_ < 0)
1184     phiBeam_ += (2 * M_PI);
1185   if (vertIndex != -1 && vertIndex < static_cast<int>(SimVtx->size())) {
1186     edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
1187     for (int iv = 0; iv < vertIndex; iv++)
1188       simVtxItr++;
1189     edm::LogVerbatim("HGCSim") << "Vertex " << vertIndex << " position " << simVtxItr->position();
1190     xBeam_ = verX[0];
1191     yBeam_ = verY[0];
1192     zBeam_ = verZ[0];
1193   }
1194 }
1195 
1196 template <class T1>
1197 void HGCalTBAnalyzer::analyzeDigi(int type, const T1& detId, uint16_t adc) {
1198   DetId id1 = DetId(detId.rawId());
1199   GlobalPoint global = hgeom_[type]->getPosition(id1);
1200   hDigiOcc_[type]->Fill(global.x(), global.y());
1201   hDigiLng_[type]->Fill(global.z());
1202   hDigiADC_[type]->Fill(adc);
1203 }
1204 
1205 void HGCalTBAnalyzer::analyzeRecHits(int type, const edm::Handle<HGCRecHitCollection>& hits) {
1206   std::map<int, double> map_hitLayer;
1207   std::map<int, std::pair<DetId, double>> map_hitCell;
1208   for (const auto& it : *hits) {
1209     DetId detId = it.id();
1210     GlobalPoint global = hgeom_[type]->getPosition(detId);
1211     double energy = it.energy();
1212     int layer = HGCalDetId(detId).layer();
1213     int cell = HGCalDetId(detId).cell();
1214 #ifdef EDM_ML_DEBUG
1215     edm::LogVerbatim("HGCSim") << "Layer thickness " << hgcons_[type]->waferTypeL(HGCalDetId(detId).wafer());
1216 #endif
1217     hRecHitOcc_[type]->Fill(global.x(), global.y(), energy);
1218     hRecHitE_[type]->Fill(energy);
1219     if (map_hitLayer.count(layer) != 0) {
1220       map_hitLayer[layer] += energy;
1221     } else {
1222       map_hitLayer[layer] = energy;
1223     }
1224     if (map_hitCell.count(cell) != 0) {
1225       double ee = energy + map_hitCell[cell].second;
1226       map_hitCell[cell] = std::pair<uint32_t, double>(detId, ee);
1227     } else {
1228       map_hitCell[cell] = std::pair<uint32_t, double>(detId, energy);
1229     }
1230 #ifdef EDM_ML_DEBUG
1231     edm::LogVerbatim("HGCSim") << "RecHit: " << layer << " " << global.x() << " " << global.y() << " " << global.z()
1232                                << " " << energy;
1233 #endif
1234   }
1235 
1236   for (const auto& itr : map_hitLayer) {
1237     int layer = itr.first;
1238     double energy = itr.second;
1239     double zp = hgcons_[type]->waferZ(layer, true);
1240 #ifdef EDM_ML_DEBUG
1241     edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer << " " << zp << " " << energy;
1242 #endif
1243     hRecHitLng_[type]->Fill(zp, energy);
1244     hRecHitLng1_[type]->Fill(layer, energy);
1245   }
1246 
1247   for (const auto& itr : map_hitCell) {
1248     DetId detId = ((itr.second).first);
1249     double energy = ((itr.second).second);
1250     GlobalPoint global = hgeom_[type]->getPosition(detId);
1251     hRecHitLat_[type]->Fill(global.x(), global.y(), energy);
1252   }
1253 }
1254 
1255 void HGCalTBAnalyzer::analyzePassiveHits(edm::Handle<edm::PassiveHitContainer> const& hgcPH, int subdet) {
1256   for (const auto& v : *hgcPH) {
1257     double energy = v.energy();
1258     std::string name = v.vname();
1259     unsigned int id = v.id();
1260 #ifdef EDM_ML_DEBUG
1261     double time = v.time();
1262     edm::LogVerbatim("HGCSim") << "HGCalTBAnalyzer::analyzePassiveHits:Energy:"
1263                                << "Time:Name:Id : " << energy << ":" << time << ":" << name << ":" << id;
1264 #endif
1265 
1266     if (subdet == 1) {
1267       hgcPassiveEEEnergy_.push_back(energy);
1268       hgcPassiveEEName_.push_back(name);
1269       hgcPassiveEEID_.push_back(id);
1270     } else if (subdet == 2) {
1271       hgcPassiveFHEnergy_.push_back(energy);
1272       hgcPassiveFHName_.push_back(name);
1273       hgcPassiveFHID_.push_back(id);
1274     } else if (subdet == 3) {
1275       hgcPassiveBHEnergy_.push_back(energy);
1276       hgcPassiveBHName_.push_back(name);
1277       hgcPassiveBHID_.push_back(id);
1278     } else if (subdet == 4) {
1279       hgcPassiveCMSEEnergy_.push_back(energy);
1280       hgcPassiveCMSEName_.push_back(name);
1281       hgcPassiveCMSEID_.push_back(id);
1282     } else if (subdet == 5) {
1283       hgcPassiveBeamEnergy_.push_back(energy);
1284       hgcPassiveBeamName_.push_back(name);
1285       hgcPassiveBeamID_.push_back(id);
1286     }
1287   }
1288 }
1289 
1290 bool HGCalTBAnalyzer::sortTime(const std::pair<double, double>& i, const std::pair<double, double>& j) {
1291   return i.first < j.first;
1292 }
1293 
1294 // define this as a plug-in
1295 DEFINE_FWK_MODULE(HGCalTBAnalyzer);