Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:04:14

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