File indexing completed on 2024-04-06 12:32:27
0001 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0003
0004 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016
0017 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0018 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0019 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0020
0021 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0022 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0023
0024 #include "TH1D.h"
0025 #include "TH2D.h"
0026 #include "TProfile.h"
0027
0028 #include <fstream>
0029 #include <iostream>
0030 #include <map>
0031 #include <string>
0032 #include <vector>
0033
0034 class HcalSimHitCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0035 public:
0036 HcalSimHitCheck(const edm::ParameterSet &ps);
0037
0038 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0039
0040 protected:
0041 void beginRun(edm::Run const &, edm::EventSetup const &) override;
0042 void analyze(edm::Event const &, edm::EventSetup const &) override;
0043 void endRun(edm::Run const &, edm::EventSetup const &) override {}
0044
0045 void analyzeHits(std::vector<PCaloHit> &);
0046
0047 private:
0048 const HcalDDDRecConstants *hcons_;
0049 int maxDepthHB_, maxDepthHE_;
0050 int maxDepthHO_, maxDepthHF_;
0051 int maxDepth_;
0052
0053 int iphi_bins;
0054 float iphi_min, iphi_max;
0055 int ieta_bins_HB;
0056 float ieta_min_HB, ieta_max_HB;
0057 int ieta_bins_HE;
0058 float ieta_min_HE, ieta_max_HE;
0059 int ieta_bins_HO;
0060 float ieta_min_HO, ieta_max_HO;
0061 int ieta_bins_HF;
0062 float ieta_min_HF, ieta_max_HF;
0063
0064 const std::string g4Label, hcalHits, outFile_;
0065 const int verbose_;
0066 const bool checkHit_, testNumber_, hep17_;
0067
0068 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0069 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
0070
0071 TH1D *meAllNHit_, *meBadDetHit_, *meBadSubHit_, *meBadIdHit_;
0072 TH1D *meHBNHit_, *meHENHit_, *meHONHit_, *meHFNHit_;
0073 TH1D *meDetectHit_, *meSubdetHit_, *meDepthHit_, *meEtaHit_;
0074 TH2D *meEtaPhiHit_;
0075 std::vector<TH2D *> meEtaPhiHitDepth_;
0076 TH1D *mePhiHit_, *mePhiHitb_, *meEnergyHit_, *meTimeHit_, *meTimeWHit_;
0077 TH1D *meHBDepHit_, *meHEDepHit_, *meHODepHit_, *meHFDepHit_, *meHFDepHitw_;
0078 TH1D *meHBEtaHit_, *meHEEtaHit_, *meHOEtaHit_, *meHFEtaHit_;
0079 TH1D *meHBPhiHit_, *meHEPhiHit_, *meHOPhiHit_, *meHFPhiHit_;
0080 TH1D *meHBEneHit_, *meHEEneHit_, *meHOEneHit_, *meHFEneHit_;
0081 TH2D *meHBEneMap_, *meHEEneMap_, *meHOEneMap_, *meHFEneMap_;
0082 TH1D *meHBEneSum_, *meHEEneSum_, *meHOEneSum_, *meHFEneSum_;
0083 TProfile *meHBEneSum_vs_ieta_, *meHEEneSum_vs_ieta_, *meHOEneSum_vs_ieta_, *meHFEneSum_vs_ieta_;
0084 TH1D *meHBTimHit_, *meHETimHit_, *meHOTimHit_, *meHFTimHit_;
0085 TH1D *meHBEneHit2_, *meHEEneHit2_, *meHOEneHit2_, *meHFEneHit2_;
0086 TH1D *meHBL10Ene_, *meHEL10Ene_, *meHOL10Ene_, *meHFL10Ene_;
0087 TProfile *meHBL10EneP_, *meHEL10EneP_, *meHOL10EneP_, *meHFL10EneP_;
0088 TH1D *meHEP17EneHit_, *meHEP17EneHit2_;
0089 };
0090
0091 HcalSimHitCheck::HcalSimHitCheck(const edm::ParameterSet &ps)
0092 : g4Label(ps.getParameter<std::string>("moduleLabel")),
0093 hcalHits(ps.getParameter<std::string>("HitCollection")),
0094 outFile_(ps.getParameter<std::string>("outputFile")),
0095 verbose_(ps.getParameter<int>("Verbose")),
0096 checkHit_(true),
0097 testNumber_(ps.getParameter<bool>("TestNumber")),
0098 hep17_(ps.getParameter<bool>("hep17")),
0099 tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hcalHits))),
0100 tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {
0101 edm::LogVerbatim("HcalSim") << "Module Label: " << g4Label << " Hits: " << hcalHits << " / " << checkHit_
0102 << " Output: " << outFile_;
0103 }
0104
0105 void HcalSimHitCheck::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0106 edm::ParameterSetDescription desc;
0107 desc.add<std::string>("moduleLabel", "g4SimHits");
0108 desc.add<std::string>("HitCollection", "HcalHits");
0109 desc.add<std::string>("outputFile", "hcHit.root");
0110 desc.add<int>("Verbose", 0);
0111 desc.add<bool>("TestNumber", true);
0112 desc.add<bool>("hep17", false);
0113 descriptions.add("hcalSimHitCheck", desc);
0114 }
0115
0116 void HcalSimHitCheck::beginRun(const edm::Run &, const edm::EventSetup &es) {
0117 hcons_ = &es.getData(tok_HRNDC_);
0118 maxDepthHB_ = hcons_->getMaxDepth(0);
0119 maxDepthHE_ = hcons_->getMaxDepth(1);
0120 maxDepthHF_ = hcons_->getMaxDepth(2);
0121 maxDepthHO_ = hcons_->getMaxDepth(3);
0122 maxDepth_ = (maxDepthHB_ > maxDepthHE_ ? maxDepthHB_ : maxDepthHE_);
0123 maxDepth_ = (maxDepth_ > maxDepthHF_ ? maxDepth_ : maxDepthHF_);
0124 maxDepth_ = (maxDepth_ > maxDepthHO_ ? maxDepth_ : maxDepthHO_);
0125
0126
0127
0128
0129 int NphiMax = hcons_->getNPhi(0);
0130
0131 NphiMax = (hcons_->getNPhi(1) > NphiMax ? hcons_->getNPhi(1) : NphiMax);
0132 NphiMax = (hcons_->getNPhi(2) > NphiMax ? hcons_->getNPhi(2) : NphiMax);
0133 NphiMax = (hcons_->getNPhi(3) > NphiMax ? hcons_->getNPhi(3) : NphiMax);
0134
0135
0136 iphi_min = 0.5;
0137 iphi_max = NphiMax + 0.5;
0138 iphi_bins = (int)(iphi_max - iphi_min);
0139
0140 int iEtaHBMax = hcons_->getEtaRange(0).second;
0141 int iEtaHEMax = std::max(hcons_->getEtaRange(1).second, 1);
0142 int iEtaHFMax = hcons_->getEtaRange(2).second;
0143 int iEtaHOMax = hcons_->getEtaRange(3).second;
0144
0145
0146
0147
0148 int iEtaMax = (iEtaHBMax > iEtaHEMax ? iEtaHBMax : iEtaHEMax);
0149 iEtaMax = (iEtaMax > iEtaHFMax ? iEtaMax : iEtaHFMax);
0150 iEtaMax = (iEtaMax > iEtaHOMax ? iEtaMax : iEtaHOMax);
0151
0152 iEtaHBMax = iEtaMax;
0153 iEtaHEMax = iEtaMax;
0154 iEtaHFMax = iEtaMax;
0155 iEtaHOMax = iEtaMax;
0156
0157
0158
0159 ieta_min_HB = -iEtaHBMax - 1.5;
0160 ieta_max_HB = iEtaHBMax + 1.5;
0161 ieta_bins_HB = (int)(ieta_max_HB - ieta_min_HB);
0162
0163 ieta_min_HE = -iEtaHEMax - 1.5;
0164 ieta_max_HE = iEtaHEMax + 1.5;
0165 ieta_bins_HE = (int)(ieta_max_HE - ieta_min_HE);
0166
0167 ieta_min_HF = -iEtaHFMax - 1.5;
0168 ieta_max_HF = iEtaHFMax + 1.5;
0169 ieta_bins_HF = (int)(ieta_max_HF - ieta_min_HF);
0170
0171 ieta_min_HO = -iEtaHOMax - 1.5;
0172 ieta_max_HO = iEtaHOMax + 1.5;
0173 ieta_bins_HO = (int)(ieta_max_HO - ieta_min_HO);
0174
0175 Char_t hname[100];
0176 Char_t htitle[100];
0177
0178 edm::Service<TFileService> fs;
0179
0180
0181 if (checkHit_) {
0182 meAllNHit_ = fs->make<TH1D>("Hit01", "Number of Hits in HCal", 20000, 0., 20000.);
0183 meBadDetHit_ = fs->make<TH1D>("Hit02", "Hits with wrong Det", 100, 0., 100.);
0184 meBadSubHit_ = fs->make<TH1D>("Hit03", "Hits with wrong Subdet", 100, 0., 100.);
0185 meBadIdHit_ = fs->make<TH1D>("Hit04", "Hits with wrong ID", 100, 0., 100.);
0186 meHBNHit_ = fs->make<TH1D>("Hit05", "Number of Hits in HB", 20000, 0., 20000.);
0187 meHENHit_ = fs->make<TH1D>("Hit06", "Number of Hits in HE", 10000, 0., 10000.);
0188 meHONHit_ = fs->make<TH1D>("Hit07", "Number of Hits in HO", 10000, 0., 10000.);
0189 meHFNHit_ = fs->make<TH1D>("Hit08", "Number of Hits in HF", 10000, 0., 10000.);
0190 meHBNHit_->Sumw2();
0191 meHENHit_->Sumw2();
0192 meHONHit_->Sumw2();
0193 meHFNHit_->Sumw2();
0194 meDetectHit_ = fs->make<TH1D>("Hit09", "Detector ID", 50, 0., 50.);
0195 meSubdetHit_ = fs->make<TH1D>("Hit10", "Subdetectors in HCal", 50, 0., 50.);
0196 meDepthHit_ = fs->make<TH1D>("Hit11", "Depths in HCal", 20, 0., 20.);
0197 meEtaHit_ = fs->make<TH1D>("Hit12", "Eta in HCal", ieta_bins_HF, ieta_min_HF, ieta_max_HF);
0198 meEtaPhiHit_ = fs->make<TH2D>(
0199 "Hit12b", "Eta-phi in HCal", ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max);
0200 for (int depth = 1; depth <= maxDepth_; depth++) {
0201 sprintf(hname, "Hit12bd%d", depth);
0202 sprintf(htitle, "Eta-phi in HCal d%d", depth);
0203 meEtaPhiHitDepth_.emplace_back(
0204 fs->make<TH2D>(hname, htitle, ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max));
0205 }
0206
0207
0208 mePhiHit_ = fs->make<TH1D>("Hit13", "Phi in HCal (HB,HO)", iphi_bins, iphi_min, iphi_max);
0209 mePhiHitb_ = fs->make<TH1D>("Hit13b", "Phi in HCal (HE,HF)", iphi_bins, iphi_min, iphi_max);
0210 meEnergyHit_ = fs->make<TH1D>("Hit14", "Energy in HCal", 2000, 0., 20.);
0211 meTimeHit_ = fs->make<TH1D>("Hit15", "Time in HCal", 528, 0., 528.);
0212 meTimeWHit_ = fs->make<TH1D>("Hit16", "Time in HCal (E wtd)", 528, 0., 528.);
0213 meHBDepHit_ = fs->make<TH1D>("Hit17", "Depths in HB", 20, 0., 20.);
0214 meHEDepHit_ = fs->make<TH1D>("Hit18", "Depths in HE", 20, 0., 20.);
0215 meHODepHit_ = fs->make<TH1D>("Hit19", "Depths in HO", 20, 0., 20.);
0216 meHFDepHit_ = fs->make<TH1D>("Hit20", "Depths in HF", 20, 0., 20.);
0217 meHBDepHit_->Sumw2();
0218 meHEDepHit_->Sumw2();
0219 meHODepHit_->Sumw2();
0220 meHFDepHit_->Sumw2();
0221 meHFDepHitw_ = fs->make<TH1D>("Hit20b", "Depths in HF (p.e. weighted)", 20, 0., 20.);
0222 meHBEtaHit_ = fs->make<TH1D>("Hit21", "Eta in HB", ieta_bins_HB, ieta_min_HB, ieta_max_HB);
0223 meHEEtaHit_ = fs->make<TH1D>("Hit22", "Eta in HE", ieta_bins_HE, ieta_min_HE, ieta_max_HE);
0224 meHOEtaHit_ = fs->make<TH1D>("Hit23", "Eta in HO", ieta_bins_HO, ieta_min_HO, ieta_max_HO);
0225 meHFEtaHit_ = fs->make<TH1D>("Hit24", "Eta in HF", ieta_bins_HF, ieta_min_HF, ieta_max_HF);
0226 meHBEtaHit_->Sumw2();
0227 meHEEtaHit_->Sumw2();
0228 meHOEtaHit_->Sumw2();
0229 meHFEtaHit_->Sumw2();
0230 meHBPhiHit_ = fs->make<TH1D>("Hit25", "Phi in HB", iphi_bins, iphi_min, iphi_max);
0231 meHEPhiHit_ = fs->make<TH1D>("Hit26", "Phi in HE", iphi_bins, iphi_min, iphi_max);
0232 meHOPhiHit_ = fs->make<TH1D>("Hit27", "Phi in HO", iphi_bins, iphi_min, iphi_max);
0233 meHFPhiHit_ = fs->make<TH1D>("Hit28", "Phi in HF", iphi_bins, iphi_min, iphi_max);
0234 meHBPhiHit_->Sumw2();
0235 meHEPhiHit_->Sumw2();
0236 meHOPhiHit_->Sumw2();
0237 meHFPhiHit_->Sumw2();
0238 meHBEneHit_ = fs->make<TH1D>("Hit29", "Energy in HB", 2000, 0., 20.);
0239 meHEEneHit_ = fs->make<TH1D>("Hit30", "Energy in HE", 500, 0., 5.);
0240 meHEP17EneHit_ = fs->make<TH1D>("Hit30b", "Energy in HEP17", 500, 0., 5.);
0241 meHOEneHit_ = fs->make<TH1D>("Hit31", "Energy in HO", 500, 0., 5.);
0242 meHFEneHit_ = fs->make<TH1D>("Hit32", "Energy in HF", 1001, -0.5, 1000.5);
0243 meHBEneHit_->Sumw2();
0244 meHEEneHit_->Sumw2();
0245 meHOEneHit_->Sumw2();
0246 meHFEneHit_->Sumw2();
0247
0248
0249
0250
0251 meHBEneMap_ =
0252 fs->make<TH2D>("HBEneMap", "HBEneMap", ieta_bins_HB, ieta_min_HB, ieta_max_HB, iphi_bins, iphi_min, iphi_max);
0253 meHEEneMap_ =
0254 fs->make<TH2D>("HEEneMap", "HEEneMap", ieta_bins_HE, ieta_min_HE, ieta_max_HE, iphi_bins, iphi_min, iphi_max);
0255 meHOEneMap_ =
0256 fs->make<TH2D>("HOEneMap", "HOEneMap", ieta_bins_HO, ieta_min_HO, ieta_max_HO, iphi_bins, iphi_min, iphi_max);
0257 meHFEneMap_ =
0258 fs->make<TH2D>("HFEneMap", "HFEneMap", ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max);
0259
0260 meHBEneSum_ = fs->make<TH1D>("HBEneSum", "HBEneSum", 2000, 0., 20.);
0261 meHEEneSum_ = fs->make<TH1D>("HEEneSum", "HEEneSum", 500, 0., 5.);
0262 meHOEneSum_ = fs->make<TH1D>("HOEneSum", "HOEneSum", 500, 0., 5.);
0263 meHFEneSum_ = fs->make<TH1D>("HFEneSum", "HFEneSum", 1001, -0.5, 1000.5);
0264 meHBEneSum_->Sumw2();
0265 meHEEneSum_->Sumw2();
0266 meHOEneSum_->Sumw2();
0267 meHFEneSum_->Sumw2();
0268
0269 meHBEneSum_vs_ieta_ = fs->make<TProfile>(
0270 "HBEneSum_vs_ieta", "HBEneSum_vs_ieta", ieta_bins_HB, ieta_min_HB, ieta_max_HB, -10.5, 2000.5);
0271 meHEEneSum_vs_ieta_ = fs->make<TProfile>(
0272 "HEEneSum_vs_ieta", "HEEneSum_vs_ieta", ieta_bins_HE, ieta_min_HE, ieta_max_HE, -10.5, 2000.5);
0273 meHOEneSum_vs_ieta_ = fs->make<TProfile>(
0274 "HOEneSum_vs_ieta", "HOEneSum_vs_ieta", ieta_bins_HO, ieta_min_HO, ieta_max_HO, -10.5, 2000.5);
0275 meHFEneSum_vs_ieta_ = fs->make<TProfile>(
0276 "HFEneSum_vs_ieta", "HFEneSum_vs_ieta", ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10.5, 2000.5);
0277
0278 meHBTimHit_ = fs->make<TH1D>("Hit33", "Time in HB", 528, 0., 528.);
0279 meHETimHit_ = fs->make<TH1D>("Hit34", "Time in HE", 528, 0., 528.);
0280 meHOTimHit_ = fs->make<TH1D>("Hit35", "Time in HO", 528, 0., 528.);
0281 meHFTimHit_ = fs->make<TH1D>("Hit36", "Time in HF", 528, 0., 528.);
0282 meHBTimHit_->Sumw2();
0283 meHETimHit_->Sumw2();
0284 meHOTimHit_->Sumw2();
0285 meHFTimHit_->Sumw2();
0286
0287 meHBEneHit2_ = fs->make<TH1D>("Hit37", "Energy in HB 2", 100, 0., 0.0001);
0288 meHEEneHit2_ = fs->make<TH1D>("Hit38", "Energy in HE 2", 100, 0., 0.0001);
0289 meHEP17EneHit2_ = fs->make<TH1D>("Hit38b", "Energy in HEP17 2", 100, 0., 0.0001);
0290 meHOEneHit2_ = fs->make<TH1D>("Hit39", "Energy in HO 2", 100, 0., 0.0001);
0291 meHFEneHit2_ = fs->make<TH1D>("Hit40", "Energy in HF 2", 100, 0.5, 100.5);
0292 meHBL10Ene_ = fs->make<TH1D>("Hit41", "Log10Energy in HB", 140, -10., 4.);
0293 meHEL10Ene_ = fs->make<TH1D>("Hit42", "Log10Energy in HE", 140, -10., 4.);
0294 meHFL10Ene_ = fs->make<TH1D>("Hit43", "Log10Energy in HF", 50, -1., 4.);
0295 meHOL10Ene_ = fs->make<TH1D>("Hit44", "Log10Energy in HO", 140, -10., 4.);
0296 meHBL10EneP_ = fs->make<TProfile>("Hit45", "Log10Energy in HB vs Hit contribution", 140, -10., 4., 0., 1.);
0297 meHEL10EneP_ = fs->make<TProfile>("Hit46", "Log10Energy in HE vs Hit contribution", 140, -10., 4., 0., 1.);
0298 meHFL10EneP_ = fs->make<TProfile>("Hit47", "Log10Energy in HF vs Hit contribution", 140, -10., 4., 0., 1.);
0299 meHOL10EneP_ = fs->make<TProfile>("Hit48", "Log10Energy in HO vs Hit contribution", 140, -10., 4., 0., 1.);
0300 }
0301 }
0302
0303 void HcalSimHitCheck::analyze(const edm::Event &e, const edm::EventSetup &) {
0304 if (verbose_ > 0)
0305 edm::LogVerbatim("HcalSim") << "Run = " << e.id().run() << " Event = " << e.id().event();
0306
0307 bool getHits = false;
0308 if (checkHit_) {
0309 const edm::Handle<edm::PCaloHitContainer> &hitsHcal = e.getHandle(tok_hits_);
0310 if (hitsHcal.isValid()) {
0311 getHits = true;
0312 if (verbose_ > 0)
0313 edm::LogVerbatim("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
0314 std::vector<PCaloHit> caloHits;
0315 caloHits.insert(caloHits.end(), hitsHcal->begin(), hitsHcal->end());
0316 if (verbose_ > 0)
0317 edm::LogVerbatim("HcalSim") << "HcalValidation: Hit buffer " << caloHits.size();
0318 analyzeHits(caloHits);
0319 } else {
0320 if (verbose_ > 0)
0321 edm::LogVerbatim("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
0322 }
0323 }
0324 }
0325
0326 void HcalSimHitCheck::analyzeHits(std::vector<PCaloHit> &hits) {
0327 int nHit = hits.size();
0328 int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nBad1 = 0, nBad2 = 0, nBad = 0;
0329 std::vector<double> encontHB(140, 0.);
0330 std::vector<double> encontHE(140, 0.);
0331 std::vector<double> encontHF(140, 0.);
0332 std::vector<double> encontHO(140, 0.);
0333 double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0;
0334
0335 double HBEneMap[ieta_bins_HB][iphi_bins];
0336 double HEEneMap[ieta_bins_HE][iphi_bins];
0337 double HOEneMap[ieta_bins_HO][iphi_bins];
0338 double HFEneMap[ieta_bins_HF][iphi_bins];
0339
0340
0341 int eta_offset_HB = -(int)ieta_min_HB;
0342 int eta_offset_HE = -(int)ieta_min_HE;
0343 int eta_offset_HO = -(int)ieta_min_HO;
0344 int eta_offset_HF = -(int)ieta_min_HF;
0345
0346 for (int i = 0; i < ieta_bins_HB; i++) {
0347 for (int j = 0; j < iphi_bins; j++) {
0348 HBEneMap[i][j] = 0.;
0349 }
0350 }
0351
0352 for (int i = 0; i < ieta_bins_HE; i++) {
0353 for (int j = 0; j < iphi_bins; j++) {
0354 HEEneMap[i][j] = 0.;
0355 }
0356 }
0357
0358 for (int i = 0; i < ieta_bins_HO; i++) {
0359 for (int j = 0; j < iphi_bins; j++) {
0360 HOEneMap[i][j] = 0.;
0361 }
0362 }
0363
0364 for (int i = 0; i < ieta_bins_HF; i++) {
0365 for (int j = 0; j < iphi_bins; j++) {
0366 HFEneMap[i][j] = 0.;
0367 }
0368 }
0369
0370 for (int i = 0; i < nHit; i++) {
0371 double energy = hits[i].energy();
0372 double log10en = log10(energy);
0373 int log10i = int((log10en + 10.) * 10.);
0374 double time = hits[i].time();
0375 unsigned int id = hits[i].id();
0376 int det, subdet, depth, eta, phi;
0377 HcalDetId hid;
0378 if (testNumber_)
0379 hid = HcalHitRelabeller::relabel(id, hcons_);
0380 else
0381 hid = HcalDetId(id);
0382 det = hid.det();
0383 subdet = hid.subdet();
0384 depth = hid.depth();
0385 eta = hid.ieta();
0386 phi = hid.iphi();
0387
0388 if (verbose_ > 1)
0389 edm::LogVerbatim("HcalSim") << "Hit[" << i << "] ID " << std::hex << id << std::dec << " Det " << det << " Sub "
0390 << subdet << " depth " << depth << " Eta " << eta << " Phi " << phi << " E " << energy
0391 << " time " << time;
0392 if (det == 4) {
0393 if (subdet == static_cast<int>(HcalBarrel))
0394 nHB++;
0395 else if (subdet == static_cast<int>(HcalEndcap))
0396 nHE++;
0397 else if (subdet == static_cast<int>(HcalOuter))
0398 nHO++;
0399 else if (subdet == static_cast<int>(HcalForward))
0400 nHF++;
0401 else {
0402 nBad++;
0403 nBad2++;
0404 }
0405 } else {
0406 nBad++;
0407 nBad1++;
0408 }
0409
0410 meDetectHit_->Fill(double(det));
0411 if (det == 4) {
0412 meSubdetHit_->Fill(double(subdet));
0413 meDepthHit_->Fill(double(depth));
0414 meEtaHit_->Fill(double(eta));
0415 meEtaPhiHit_->Fill(double(eta), double(phi));
0416 meEtaPhiHitDepth_[depth - 1]->Fill(double(eta), double(phi));
0417
0418
0419
0420 if (subdet == static_cast<int>(HcalBarrel))
0421 mePhiHit_->Fill(double(phi));
0422 else if (subdet == static_cast<int>(HcalEndcap))
0423 mePhiHitb_->Fill(double(phi));
0424 else if (subdet == static_cast<int>(HcalOuter))
0425 mePhiHit_->Fill(double(phi));
0426 else if (subdet == static_cast<int>(HcalForward))
0427 mePhiHitb_->Fill(double(phi));
0428
0429
0430
0431 if (subdet != static_cast<int>(HcalForward)) {
0432 meEnergyHit_->Fill(energy);
0433
0434
0435
0436 meTimeWHit_->Fill(double(time), energy);
0437 }
0438 meTimeHit_->Fill(time);
0439
0440 if (subdet == static_cast<int>(HcalBarrel)) {
0441 meHBDepHit_->Fill(double(depth));
0442 meHBEtaHit_->Fill(double(eta));
0443 meHBPhiHit_->Fill(double(phi));
0444 meHBEneHit_->Fill(energy);
0445 meHBEneHit2_->Fill(energy);
0446 meHBTimHit_->Fill(time);
0447 meHBL10Ene_->Fill(log10en);
0448 if (log10i >= 0 && log10i < 140)
0449 encontHB[log10i] += energy;
0450 entotHB += energy;
0451
0452 HBEneMap[eta + eta_offset_HB][phi - 1] += energy;
0453
0454 } else if (subdet == static_cast<int>(HcalEndcap)) {
0455 meHEDepHit_->Fill(double(depth));
0456 meHEEtaHit_->Fill(double(eta));
0457 meHEPhiHit_->Fill(double(phi));
0458
0459 bool isHEP17 = (phi >= 63) && (phi <= 66) && (eta > 0);
0460 if (hep17_) {
0461 if (!isHEP17) {
0462 meHEEneHit_->Fill(energy);
0463 meHEEneHit2_->Fill(energy);
0464 } else {
0465 meHEP17EneHit_->Fill(energy);
0466 meHEP17EneHit2_->Fill(energy);
0467 }
0468 } else {
0469 meHEEneHit_->Fill(energy);
0470 meHEEneHit2_->Fill(energy);
0471 }
0472
0473 meHETimHit_->Fill(time);
0474 meHEL10Ene_->Fill(log10en);
0475 if (log10i >= 0 && log10i < 140)
0476 encontHE[log10i] += energy;
0477 entotHE += energy;
0478
0479 HEEneMap[eta + eta_offset_HE][phi - 1] += energy;
0480
0481 } else if (subdet == static_cast<int>(HcalOuter)) {
0482 meHODepHit_->Fill(double(depth));
0483 meHOEtaHit_->Fill(double(eta));
0484 meHOPhiHit_->Fill(double(phi));
0485 meHOEneHit_->Fill(energy);
0486 meHOEneHit2_->Fill(energy);
0487 meHOTimHit_->Fill(time);
0488 meHOL10Ene_->Fill(log10en);
0489 if (log10i >= 0 && log10i < 140)
0490 encontHO[log10i] += energy;
0491 entotHO += energy;
0492
0493 HOEneMap[eta + eta_offset_HO][phi - 1] += energy;
0494
0495 } else if (subdet == static_cast<int>(HcalForward)) {
0496 meHFDepHit_->Fill(double(depth));
0497 meHFDepHitw_->Fill(double(depth), energy);
0498 meHFEtaHit_->Fill(double(eta));
0499 meHFPhiHit_->Fill(double(phi));
0500 meHFEneHit_->Fill(energy);
0501 meHFEneHit2_->Fill(energy);
0502 meHFTimHit_->Fill(time);
0503 meHFL10Ene_->Fill(log10en);
0504 if (log10i >= 0 && log10i < 140)
0505 encontHF[log10i] += energy;
0506 entotHF += energy;
0507
0508 HFEneMap[eta + eta_offset_HF][phi - 1] += energy;
0509 }
0510 }
0511 }
0512 if (entotHB != 0)
0513 for (int i = 0; i < 140; i++)
0514 meHBL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHB[i] / entotHB);
0515 if (entotHE != 0)
0516 for (int i = 0; i < 140; i++)
0517 meHEL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHE[i] / entotHE);
0518 if (entotHF != 0)
0519 for (int i = 0; i < 140; i++)
0520 meHFL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHF[i] / entotHF);
0521 if (entotHO != 0)
0522 for (int i = 0; i < 140; i++)
0523 meHOL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHO[i] / entotHO);
0524
0525 meAllNHit_->Fill(double(nHit));
0526 meBadDetHit_->Fill(double(nBad1));
0527 meBadSubHit_->Fill(double(nBad2));
0528 meBadIdHit_->Fill(double(nBad));
0529 meHBNHit_->Fill(double(nHB));
0530 meHENHit_->Fill(double(nHE));
0531 meHONHit_->Fill(double(nHO));
0532 meHFNHit_->Fill(double(nHF));
0533
0534 for (int i = 0; i < ieta_bins_HB; i++) {
0535 for (int j = 0; j < iphi_bins; j++) {
0536 if (HBEneMap[i][j] != 0) {
0537 meHBEneSum_->Fill(HBEneMap[i][j]);
0538 meHBEneSum_vs_ieta_->Fill((i - eta_offset_HB), HBEneMap[i][j]);
0539 meHBEneMap_->Fill((i - eta_offset_HB), j + 1, HBEneMap[i][j]);
0540 }
0541 }
0542 }
0543
0544 for (int i = 0; i < ieta_bins_HE; i++) {
0545 for (int j = 0; j < iphi_bins; j++) {
0546 if (HEEneMap[i][j] != 0) {
0547 meHEEneSum_->Fill(HEEneMap[i][j]);
0548 meHEEneSum_vs_ieta_->Fill((i - eta_offset_HE), HEEneMap[i][j]);
0549 meHEEneMap_->Fill((i - eta_offset_HE), j + 1, HEEneMap[i][j]);
0550 }
0551 }
0552 }
0553
0554 for (int i = 0; i < ieta_bins_HO; i++) {
0555 for (int j = 0; j < iphi_bins; j++) {
0556 if (HOEneMap[i][j] != 0) {
0557 meHOEneSum_->Fill(HOEneMap[i][j]);
0558 meHOEneSum_vs_ieta_->Fill((i - eta_offset_HO), HOEneMap[i][j]);
0559 meHOEneMap_->Fill((i - eta_offset_HO), j + 1, HOEneMap[i][j]);
0560 }
0561 }
0562 }
0563
0564 for (int i = 0; i < ieta_bins_HF; i++) {
0565 for (int j = 0; j < iphi_bins; j++) {
0566 if (HFEneMap[i][j] != 0) {
0567 meHFEneSum_->Fill(HFEneMap[i][j]);
0568 meHFEneSum_vs_ieta_->Fill((i - eta_offset_HF), HFEneMap[i][j]);
0569 meHFEneMap_->Fill((i - eta_offset_HF), j + 1, HFEneMap[i][j]);
0570 }
0571 }
0572 }
0573
0574 if (verbose_ > 0)
0575 edm::LogVerbatim("HcalSim") << "HcalSimHitCheck::analyzeHits: HB " << nHB << " HE " << nHE << " HO " << nHO
0576 << " HF " << nHF << " Bad " << nBad << " All " << nHit;
0577 }
0578
0579 #include "FWCore/Framework/interface/MakerMacros.h"
0580 DEFINE_FWK_MODULE(HcalSimHitCheck);