File indexing completed on 2024-09-07 04:38:09
0001
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <map>
0006 #include <memory>
0007 #include <string>
0008 #include <vector>
0009
0010
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/HGCalTBCommonData/interface/HGCalTBDDDConstants.h"
0028 #include "Geometry/HGCalGeometry/interface/HGCalTBGeometry.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
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
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 HGCalTBDDDConstants* hgcons_[2];
0078 const HGCalTBGeometry* 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<HGCalTBDDDConstants, IdealGeometryRecord> tokDDDEE_;
0106 const edm::ESGetToken<HGCalTBGeometry, IdealGeometryRecord> tokGeomEE_;
0107 const edm::ESGetToken<HGCalTBDDDConstants, IdealGeometryRecord> tokDDDFH_;
0108 const edm::ESGetToken<HGCalTBGeometry, 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<HGCalTBDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0213 edm::ESInputTag("", detectorEE_))),
0214 tokGeomEE_(esConsumes<HGCalTBGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
0215 edm::ESInputTag("", detectorEE_))),
0216 tokDDDFH_(esConsumes<HGCalTBDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0217 edm::ESInputTag("", detectorFH_))),
0218 tokGeomFH_(esConsumes<HGCalTBGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
0219 edm::ESInputTag("", detectorFH_))) {
0220 usesResource("TFileService");
0221 ahcalGeom_ = std::make_unique<AHCalGeometry>(iConfig);
0222
0223
0224
0225
0226
0227
0228
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);
0315 desc.addUntracked<double>("deltaY", 30.0);
0316 desc.addUntracked<double>("deltaZ", 81.0);
0317 desc.addUntracked<double>("zFirst", 17.6);
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
0441 tree_->Branch("simHitCellTime15MipEE", &simHitCellTime15MipEE_);
0442 tree_->Branch("simHitCellTime15MipFH", &simHitCellTime15MipFH_);
0443
0444 tree_->Branch("simHitCellTimeLastHitEE", &simHitCellTimeLastHitEE_);
0445 tree_->Branch("simHitCellTimeLastHitFH", &simHitCellTimeLastHitFH_);
0446
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
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
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 }
0681
0682
0683 if (doPassive_) {
0684
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
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
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
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
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
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
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 }
0788 }
0789 }
0790
0791 }
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
0805 std::map<uint32_t, std::vector<std::pair<double, double>>> map_hitTimeEn;
0806
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) {
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
0904 for (const auto& itr : map_hitTimeEn) {
0905 uint32_t id = itr.first;
0906 int wafer = -99;
0907
0908 wafer = HGCalDetId(ID[id]).wafer();
0909 double layer = HGCalDetId(id).layer();
0910 double thickness = hgcons_[type]->cellThickness(layer, wafer);
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
0923 std::sort(map_hitTimeEn[id].begin(), map_hitTimeEn[id].end(), sortTime);
0924
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
0959
0960
0961
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
1295 DEFINE_FWK_MODULE(HGCalTBAnalyzer);
