File indexing completed on 2024-07-16 22:52:43
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0023 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0024 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0025 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0026
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0029 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0030 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0031 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0032 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0033 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0034
0035 #include "FWCore/Framework/interface/Frameworkfwd.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0039 #include "FWCore/Framework/interface/ESHandle.h"
0040 #include "FWCore/ServiceRegistry/interface/Service.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "FWCore/Utilities/interface/Exception.h"
0043 #include "FWCore/Utilities/interface/EDGetToken.h"
0044 #include "FWCore/Utilities/interface/transform.h"
0045
0046 #include "SimG4CMS/Calo/interface/HGCNumberingScheme.h"
0047 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0049 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0050
0051 #include <TH1.h>
0052 #include <TH2.h>
0053 #include <TTree.h>
0054 #include <TVector3.h>
0055 #include <TSystem.h>
0056 #include <TFile.h>
0057
0058 #include <cmath>
0059 #include <memory>
0060 #include <iostream>
0061 #include <string>
0062 #include <vector>
0063
0064 class HGCHitValidation : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0065 public:
0066 explicit HGCHitValidation(const edm::ParameterSet &);
0067 ~HGCHitValidation() override = default;
0068 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0069
0070 private:
0071 typedef std::tuple<float, float, float, float> HGCHitTuple;
0072
0073 void beginJob() override;
0074 void endJob() override {}
0075 void beginRun(edm::Run const &, edm::EventSetup const &) override;
0076 void analyze(edm::Event const &, edm::EventSetup const &) override;
0077 void endRun(edm::Run const &, edm::EventSetup const &) override {}
0078 virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0079 virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0080 void analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const &simHits,
0081 int idet,
0082 TH1F *,
0083 std::map<unsigned int, HGCHitTuple> &);
0084 template <class T1>
0085 void analyzeHGCalRecHit(T1 const &theHits, std::map<unsigned int, HGCHitTuple> const &hitRefs);
0086
0087 private:
0088
0089 const std::vector<std::string> geometrySource_;
0090 const std::vector<int> ietaExcludeBH_;
0091 const bool makeTree_;
0092 const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_hgcal_;
0093 const std::vector<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>> tok_hgcalg_;
0094 std::vector<const HGCalDDDConstants *> hgcCons_;
0095 std::vector<const HGCalGeometry *> hgcGeometry_;
0096
0097 const edm::InputTag eeSimHitSource, fhSimHitSource, bhSimHitSource;
0098 const edm::EDGetTokenT<std::vector<PCaloHit>> eeSimHitToken_;
0099 const edm::EDGetTokenT<std::vector<PCaloHit>> fhSimHitToken_;
0100 const edm::EDGetTokenT<std::vector<PCaloHit>> bhSimHitToken_;
0101 const edm::InputTag eeRecHitSource, fhRecHitSource, bhRecHitSource;
0102 const edm::EDGetTokenT<HGCeeRecHitCollection> eeRecHitToken_;
0103 const edm::EDGetTokenT<HGChefRecHitCollection> fhRecHitToken_;
0104 const edm::EDGetTokenT<HGChebRecHitCollection> bhRecHitToken_;
0105
0106 TTree *hgcHits_;
0107 std::vector<float> *heeRecX_, *heeRecY_, *heeRecZ_, *heeRecEnergy_;
0108 std::vector<float> *hefRecX_, *hefRecY_, *hefRecZ_, *hefRecEnergy_;
0109 std::vector<float> *hebRecX_, *hebRecY_, *hebRecZ_, *hebRecEnergy_;
0110 std::vector<float> *heeSimX_, *heeSimY_, *heeSimZ_, *heeSimEnergy_;
0111 std::vector<float> *hefSimX_, *hefSimY_, *hefSimZ_, *hefSimEnergy_;
0112 std::vector<float> *hebSimX_, *hebSimY_, *hebSimZ_, *hebSimEnergy_;
0113 std::vector<float> *hebSimEta_, *hebRecEta_, *hebSimPhi_, *hebRecPhi_;
0114 std::vector<unsigned int> *heeDetID_, *hefDetID_, *hebDetID_;
0115
0116 TH2F *heedzVsZ_, *heedyVsY_, *heedxVsX_;
0117 TH2F *hefdzVsZ_, *hefdyVsY_, *hefdxVsX_;
0118 TH2F *hebdzVsZ_, *hebdPhiVsPhi_, *hebdEtaVsEta_;
0119 TH2F *heeRecVsSimZ_, *heeRecVsSimY_, *heeRecVsSimX_;
0120 TH2F *hefRecVsSimZ_, *hefRecVsSimY_, *hefRecVsSimX_;
0121 TH2F *hebRecVsSimZ_, *hebRecVsSimY_, *hebRecVsSimX_;
0122 TH2F *heeEnSimRec_, *hefEnSimRec_, *hebEnSimRec_;
0123 TH1F *hebEnRec_, *hebEnSim_, *hefEnRec_;
0124 TH1F *hefEnSim_, *heeEnRec_, *heeEnSim_;
0125 };
0126
0127 HGCHitValidation::HGCHitValidation(const edm::ParameterSet &cfg)
0128 : geometrySource_(cfg.getUntrackedParameter<std::vector<std::string>>("geometrySource")),
0129 ietaExcludeBH_(cfg.getParameter<std::vector<int>>("ietaExcludeBH")),
0130 makeTree_(cfg.getUntrackedParameter<bool>("makeTree", true)),
0131 tok_hgcal_{
0132 edm::vector_transform(geometrySource_,
0133 [this](const std::string &name) {
0134 return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0135 edm::ESInputTag{"", name});
0136 })},
0137 tok_hgcalg_{edm::vector_transform(
0138 geometrySource_,
0139 [this](const std::string &name) {
0140 return esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0141 })},
0142 eeSimHitSource(cfg.getParameter<edm::InputTag>("eeSimHitSource")),
0143 fhSimHitSource(cfg.getParameter<edm::InputTag>("fhSimHitSource")),
0144 bhSimHitSource(cfg.getParameter<edm::InputTag>("bhSimHitSource")),
0145 eeSimHitToken_(consumes<std::vector<PCaloHit>>(eeSimHitSource)),
0146 fhSimHitToken_(consumes<std::vector<PCaloHit>>(fhSimHitSource)),
0147 bhSimHitToken_(consumes<std::vector<PCaloHit>>(bhSimHitSource)),
0148 eeRecHitSource(cfg.getParameter<edm::InputTag>("eeRecHitSource")),
0149 fhRecHitSource(cfg.getParameter<edm::InputTag>("fhRecHitSource")),
0150 bhRecHitSource(cfg.getParameter<edm::InputTag>("bhRecHitSource")),
0151 eeRecHitToken_(consumes<HGCeeRecHitCollection>(eeRecHitSource)),
0152 fhRecHitToken_(consumes<HGChefRecHitCollection>(fhRecHitSource)),
0153 bhRecHitToken_(consumes<HGChebRecHitCollection>(bhRecHitSource)) {
0154 usesResource(TFileService::kSharedResource);
0155
0156 hgcHits_ = nullptr;
0157 heeRecX_ = heeRecY_ = heeRecZ_ = heeRecEnergy_ = nullptr;
0158 hefRecX_ = hefRecY_ = hefRecZ_ = hefRecEnergy_ = nullptr;
0159 hebRecX_ = hebRecY_ = hebRecZ_ = hebRecEnergy_ = nullptr;
0160 heeSimX_ = heeSimY_ = heeSimZ_ = heeSimEnergy_ = nullptr;
0161 hefSimX_ = hefSimY_ = hefSimZ_ = hefSimEnergy_ = nullptr;
0162 hebSimX_ = hebSimY_ = hebSimZ_ = hebSimEnergy_ = nullptr;
0163 hebSimEta_ = hebRecEta_ = hebSimPhi_ = hebRecPhi_ = nullptr;
0164 heeDetID_ = hefDetID_ = hebDetID_ = nullptr;
0165 heedzVsZ_ = heedyVsY_ = heedxVsX_ = nullptr;
0166 hefdzVsZ_ = hefdyVsY_ = hefdxVsX_ = nullptr;
0167 hebdzVsZ_ = hebdPhiVsPhi_ = hebdEtaVsEta_ = nullptr;
0168 heeRecVsSimZ_ = heeRecVsSimY_ = heeRecVsSimX_ = nullptr;
0169 hefRecVsSimZ_ = hefRecVsSimY_ = hefRecVsSimX_ = nullptr;
0170 hebRecVsSimZ_ = hebRecVsSimY_ = hebRecVsSimX_ = nullptr;
0171 heeEnSimRec_ = hefEnSimRec_ = hebEnSimRec_ = nullptr;
0172 hebEnRec_ = hebEnSim_ = hefEnRec_ = nullptr;
0173 hefEnSim_ = heeEnRec_ = heeEnSim_ = nullptr;
0174
0175 edm::LogVerbatim("HGCalValid") << "MakeTree Flag set to " << makeTree_ << " and use " << geometrySource_.size()
0176 << " Geometry sources";
0177 for (auto const &s : geometrySource_)
0178 edm::LogVerbatim("HGCalValid") << " " << s;
0179 edm::LogVerbatim("HGCalValid") << "SimHit labels: " << eeSimHitSource << " " << fhSimHitSource << " "
0180 << bhSimHitSource;
0181 edm::LogVerbatim("HGCalValid") << "RecHit labels: " << eeRecHitSource << " " << fhRecHitSource << " "
0182 << bhRecHitSource;
0183 edm::LogVerbatim("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots";
0184 for (unsigned int k = 0; k < ietaExcludeBH_.size(); ++k)
0185 edm::LogVerbatim("HGCalValid") << " [" << k << "] " << ietaExcludeBH_[k];
0186 }
0187
0188 void HGCHitValidation::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0189 std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0190 std::vector<int> etas;
0191 edm::ParameterSetDescription desc;
0192 desc.addUntracked<bool>("makeTree", true);
0193 desc.addUntracked<std::vector<std::string>>("geometrySource", names);
0194 desc.add<edm::InputTag>("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE"));
0195 desc.add<edm::InputTag>("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront"));
0196 desc.add<edm::InputTag>("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback"));
0197 desc.add<edm::InputTag>("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0198 desc.add<edm::InputTag>("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0199 desc.add<edm::InputTag>("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0200 desc.add<std::vector<int>>("ietaExcludeBH", etas);
0201 descriptions.add("hgcHitAnalysis", desc);
0202 }
0203
0204 void HGCHitValidation::beginJob() {
0205
0206 edm::Service<TFileService> fs;
0207 if (makeTree_) {
0208 hgcHits_ = fs->make<TTree>("hgcHits", "Hit Collection");
0209 hgcHits_->Branch("heeRecX", &heeRecX_);
0210 hgcHits_->Branch("heeRecY", &heeRecY_);
0211 hgcHits_->Branch("heeRecZ", &heeRecZ_);
0212 hgcHits_->Branch("heeRecEnergy", &heeRecEnergy_);
0213 hgcHits_->Branch("hefRecX", &hefRecX_);
0214 hgcHits_->Branch("hefRecY", &hefRecY_);
0215 hgcHits_->Branch("hefRecZ", &hefRecZ_);
0216 hgcHits_->Branch("hefRecEnergy", &hefRecEnergy_);
0217 hgcHits_->Branch("hebRecX", &hebRecX_);
0218 hgcHits_->Branch("hebRecY", &hebRecY_);
0219 hgcHits_->Branch("hebRecZ", &hebRecZ_);
0220 hgcHits_->Branch("hebRecEta", &hebRecEta_);
0221 hgcHits_->Branch("hebRecPhi", &hebRecPhi_);
0222 hgcHits_->Branch("hebRecEnergy", &hebRecEnergy_);
0223
0224 hgcHits_->Branch("heeSimX", &heeSimX_);
0225 hgcHits_->Branch("heeSimY", &heeSimY_);
0226 hgcHits_->Branch("heeSimZ", &heeSimZ_);
0227 hgcHits_->Branch("heeSimEnergy", &heeSimEnergy_);
0228 hgcHits_->Branch("hefSimX", &hefSimX_);
0229 hgcHits_->Branch("hefSimY", &hefSimY_);
0230 hgcHits_->Branch("hefSimZ", &hefSimZ_);
0231 hgcHits_->Branch("hefSimEnergy", &hefSimEnergy_);
0232 hgcHits_->Branch("hebSimX", &hebSimX_);
0233 hgcHits_->Branch("hebSimY", &hebSimY_);
0234 hgcHits_->Branch("hebSimZ", &hebSimZ_);
0235 hgcHits_->Branch("hebSimEta", &hebSimEta_);
0236 hgcHits_->Branch("hebSimPhi", &hebSimPhi_);
0237 hgcHits_->Branch("hebSimEnergy", &hebSimEnergy_);
0238
0239 hgcHits_->Branch("heeDetID", &heeDetID_);
0240 hgcHits_->Branch("hefDetID", &hefDetID_);
0241 hgcHits_->Branch("hebDetID", &hebDetID_);
0242 } else {
0243 heedzVsZ_ = fs->make<TH2F>("heedzVsZ", "", 7200, -360, 360, 100, -0.1, 0.1);
0244 heedyVsY_ = fs->make<TH2F>("heedyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0245 heedxVsX_ = fs->make<TH2F>("heedxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0246 heeRecVsSimZ_ = fs->make<TH2F>("heeRecVsSimZ", "", 7200, -360, 360, 7200, -360, 360);
0247 heeRecVsSimY_ = fs->make<TH2F>("heeRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0248 heeRecVsSimX_ = fs->make<TH2F>("heeRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0249 hefdzVsZ_ = fs->make<TH2F>("hefdzVsZ", "", 8200, -410, 410, 100, -0.1, 0.1);
0250 hefdyVsY_ = fs->make<TH2F>("hefdyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0251 hefdxVsX_ = fs->make<TH2F>("hefdxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0252 hefRecVsSimZ_ = fs->make<TH2F>("hefRecVsSimZ", "", 8200, -410, 410, 8200, -410, 410);
0253 hefRecVsSimY_ = fs->make<TH2F>("hefRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0254 hefRecVsSimX_ = fs->make<TH2F>("hefRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0255 hebdzVsZ_ = fs->make<TH2F>("hebdzVsZ", "", 1080, -540, 540, 100, -1.0, 1.0);
0256 hebdPhiVsPhi_ = fs->make<TH2F>("hebdPhiVsPhi", "", M_PI * 100, -0.5, M_PI + 0.5, 200, -0.2, 0.2);
0257 hebdEtaVsEta_ = fs->make<TH2F>("hebdEtaVsEta", "", 1000, -5, 5, 200, -0.1, 0.1);
0258 hebRecVsSimZ_ = fs->make<TH2F>("hebRecVsSimZ", "", 1080, -540, 540, 1080, -540, 540);
0259 hebRecVsSimY_ = fs->make<TH2F>("hebRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0260 hebRecVsSimX_ = fs->make<TH2F>("hebRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0261 heeEnRec_ = fs->make<TH1F>("heeEnRec", "", 1000, 0, 10);
0262 heeEnSim_ = fs->make<TH1F>("heeEnSim", "", 1000, 0, 0.01);
0263 heeEnSimRec_ = fs->make<TH2F>("heeEnSimRec", "", 200, 0, 0.002, 200, 0, 0.2);
0264 hefEnRec_ = fs->make<TH1F>("hefEnRec", "", 1000, 0, 10);
0265 hefEnSim_ = fs->make<TH1F>("hefEnSim", "", 1000, 0, 0.01);
0266 hefEnSimRec_ = fs->make<TH2F>("hefEnSimRec", "", 200, 0, 0.001, 200, 0, 0.5);
0267 hebEnRec_ = fs->make<TH1F>("hebEnRec", "", 1000, 0, 15);
0268 hebEnSim_ = fs->make<TH1F>("hebEnSim", "", 1000, 0, 0.01);
0269 hebEnSimRec_ = fs->make<TH2F>("hebEnSimRec", "", 200, 0, 0.02, 200, 0, 4);
0270 }
0271 }
0272
0273 void HGCHitValidation::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
0274
0275 for (size_t i = 0; i < geometrySource_.size(); i++) {
0276 edm::LogVerbatim("HGCalValid") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i;
0277 edm::ESHandle<HGCalDDDConstants> hgcCons = iSetup.getHandle(tok_hgcal_[i]);
0278 if (hgcCons.isValid()) {
0279 hgcCons_.push_back(hgcCons.product());
0280 } else {
0281 edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl;
0282 }
0283 edm::ESHandle<HGCalGeometry> hgcGeom = iSetup.getHandle(tok_hgcalg_[i]);
0284 if (hgcGeom.isValid()) {
0285 hgcGeometry_.push_back(hgcGeom.product());
0286 } else {
0287 edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for " << geometrySource_[i] << std::endl;
0288 }
0289 }
0290 }
0291
0292 void HGCHitValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0293 std::map<unsigned int, HGCHitTuple> eeHitRefs, fhHitRefs, bhHitRefs;
0294
0295
0296 const edm::Handle<std::vector<PCaloHit>> &eeSimHits = iEvent.getHandle(eeSimHitToken_);
0297
0298 if (eeSimHits.isValid()) {
0299 analyzeHGCalSimHit(eeSimHits, 0, heeEnSim_, eeHitRefs);
0300 for (std::map<unsigned int, HGCHitTuple>::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) {
0301 int idx = std::distance(eeHitRefs.begin(), itr);
0302 edm::LogVerbatim("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0303 << std::get<0>(itr->second) << "; Position"
0304 << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", "
0305 << std::get<3>(itr->second) << ")";
0306 }
0307 } else {
0308 edm::LogWarning("HGCalValid") << "No EE SimHit Found " << std::endl;
0309 }
0310
0311
0312 const edm::Handle<std::vector<PCaloHit>> &fhSimHits = iEvent.getHandle(fhSimHitToken_);
0313 if (fhSimHits.isValid()) {
0314 analyzeHGCalSimHit(fhSimHits, 1, hefEnSim_, fhHitRefs);
0315 for (std::map<unsigned int, HGCHitTuple>::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) {
0316 int idx = std::distance(fhHitRefs.begin(), itr);
0317 edm::LogVerbatim("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0318 << std::get<0>(itr->second) << "; Position"
0319 << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", "
0320 << std::get<3>(itr->second) << ")";
0321 }
0322 } else {
0323 edm::LogWarning("HGCalValid") << "No FH SimHit Found " << std::endl;
0324 }
0325
0326
0327 const edm::Handle<std::vector<PCaloHit>> &bhSimHits = iEvent.getHandle(bhSimHitToken_);
0328 if (bhSimHits.isValid()) {
0329 analyzeHGCalSimHit(bhSimHits, 2, hebEnSim_, bhHitRefs);
0330 for (std::map<unsigned int, HGCHitTuple>::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) {
0331 int idx = std::distance(bhHitRefs.begin(), itr);
0332 edm::LogVerbatim("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0333 << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second) << ", "
0334 << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ")";
0335 }
0336 } else {
0337 edm::LogWarning("HGCalValid") << "No BH SimHit Found " << std::endl;
0338 }
0339
0340
0341 const edm::Handle<HGCeeRecHitCollection> &eeRecHit = iEvent.getHandle(eeRecHitToken_);
0342 if (eeRecHit.isValid()) {
0343 const HGCeeRecHitCollection *theHits = (eeRecHit.product());
0344 for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0345 double energy = it->energy();
0346 if (!makeTree_)
0347 heeEnRec_->Fill(energy);
0348 std::map<unsigned int, HGCHitTuple>::const_iterator itr = eeHitRefs.find(it->id().rawId());
0349 if (itr != eeHitRefs.end()) {
0350 GlobalPoint xyz = hgcGeometry_[0]->getPosition(it->id());
0351 if (makeTree_) {
0352 heeRecX_->push_back(xyz.x());
0353 heeRecY_->push_back(xyz.y());
0354 heeRecZ_->push_back(xyz.z());
0355 heeRecEnergy_->push_back(energy);
0356 heeSimX_->push_back(std::get<1>(itr->second));
0357 heeSimY_->push_back(std::get<2>(itr->second));
0358 heeSimZ_->push_back(std::get<3>(itr->second));
0359 heeSimEnergy_->push_back(std::get<0>(itr->second));
0360 heeDetID_->push_back(itr->first);
0361 } else {
0362 heeRecVsSimX_->Fill(std::get<1>(itr->second), xyz.x());
0363 heeRecVsSimY_->Fill(std::get<2>(itr->second), xyz.y());
0364 heeRecVsSimZ_->Fill(std::get<3>(itr->second), xyz.z());
0365 heedxVsX_->Fill(std::get<1>(itr->second), (xyz.x() - std::get<1>(itr->second)));
0366 heedyVsY_->Fill(std::get<2>(itr->second), (xyz.y() - std::get<2>(itr->second)));
0367 heedzVsZ_->Fill(std::get<3>(itr->second), (xyz.z() - std::get<3>(itr->second)));
0368 heeEnSimRec_->Fill(std::get<0>(itr->second), energy);
0369 }
0370 edm::LogVerbatim("HGCalValid") << "EEHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0371 << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ", "
0372 << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ") Rec ("
0373 << energy << ", " << xyz.x() << ", " << xyz.y() << ", " << xyz.z();
0374 }
0375 }
0376 } else {
0377 edm::LogWarning("HGCalValid") << "No EE RecHit Found " << std::endl;
0378 }
0379
0380
0381 const edm::Handle<HGChefRecHitCollection> &fhRecHit = iEvent.getHandle(fhRecHitToken_);
0382 if (fhRecHit.isValid()) {
0383 const HGChefRecHitCollection *theHits = (fhRecHit.product());
0384 for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0385 double energy = it->energy();
0386 if (!makeTree_)
0387 hefEnRec_->Fill(energy);
0388 std::map<unsigned int, HGCHitTuple>::const_iterator itr = fhHitRefs.find(it->id().rawId());
0389 if (itr != fhHitRefs.end()) {
0390 GlobalPoint xyz = hgcGeometry_[1]->getPosition(it->id());
0391 if (makeTree_) {
0392 hefRecX_->push_back(xyz.x());
0393 hefRecY_->push_back(xyz.y());
0394 hefRecZ_->push_back(xyz.z());
0395 hefRecEnergy_->push_back(energy);
0396 hefSimX_->push_back(std::get<1>(itr->second));
0397 hefSimY_->push_back(std::get<2>(itr->second));
0398 hefSimZ_->push_back(std::get<3>(itr->second));
0399 hefSimEnergy_->push_back(std::get<0>(itr->second));
0400 hefDetID_->push_back(itr->first);
0401 } else {
0402 hefRecVsSimX_->Fill(std::get<1>(itr->second), xyz.x());
0403 hefRecVsSimY_->Fill(std::get<2>(itr->second), xyz.y());
0404 hefRecVsSimZ_->Fill(std::get<3>(itr->second), xyz.z());
0405 hefdxVsX_->Fill(std::get<1>(itr->second), (xyz.x() - std::get<1>(itr->second)));
0406 hefdyVsY_->Fill(std::get<2>(itr->second), (xyz.y() - std::get<2>(itr->second)));
0407 hefdzVsZ_->Fill(std::get<3>(itr->second), (xyz.z() - std::get<3>(itr->second)));
0408 hefEnSimRec_->Fill(std::get<0>(itr->second), energy);
0409 }
0410 edm::LogVerbatim("HGCalValid") << "FHHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0411 << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ", "
0412 << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ") Rec ("
0413 << energy << "," << xyz.x() << ", " << xyz.y() << ", " << xyz.z();
0414 }
0415 }
0416 } else {
0417 edm::LogWarning("HGCalValid") << "No FH RecHit Found " << std::endl;
0418 }
0419
0420
0421 const edm::Handle<HGChebRecHitCollection> &bhRecHit = iEvent.getHandle(bhRecHitToken_);
0422 if (bhRecHit.isValid()) {
0423 const HGChebRecHitCollection *theHits = (bhRecHit.product());
0424 analyzeHGCalRecHit(theHits, bhHitRefs);
0425 } else {
0426 edm::LogWarning("HGCalValid") << "No BH RecHit Found " << std::endl;
0427 }
0428
0429 if (makeTree_) {
0430 hgcHits_->Fill();
0431
0432 heeRecX_->clear();
0433 heeRecY_->clear();
0434 heeRecZ_->clear();
0435 hefRecX_->clear();
0436 hefRecY_->clear();
0437 hefRecZ_->clear();
0438 hebRecX_->clear();
0439 hebRecY_->clear();
0440 hebRecZ_->clear();
0441 heeRecEnergy_->clear();
0442 hefRecEnergy_->clear();
0443 hebRecEnergy_->clear();
0444 heeSimX_->clear();
0445 heeSimY_->clear();
0446 heeSimZ_->clear();
0447 hefSimX_->clear();
0448 hefSimY_->clear();
0449 hefSimZ_->clear();
0450 hebSimX_->clear();
0451 hebSimY_->clear();
0452 hebSimZ_->clear();
0453 heeSimEnergy_->clear();
0454 hefSimEnergy_->clear();
0455 hebSimEnergy_->clear();
0456 hebSimEta_->clear();
0457 hebRecEta_->clear();
0458 hebSimPhi_->clear();
0459 hebRecPhi_->clear();
0460 heeDetID_->clear();
0461 hefDetID_->clear();
0462 hebDetID_->clear();
0463 }
0464 }
0465
0466 void HGCHitValidation::analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const &simHits,
0467 int idet,
0468 TH1F *hist,
0469 std::map<unsigned int, HGCHitTuple> &hitRefs) {
0470 const HGCalTopology &hTopo = hgcGeometry_[idet]->topology();
0471 for (auto const &simHit : *simHits) {
0472 unsigned int id = simHit.id();
0473 std::pair<float, float> xy;
0474 bool ok(true);
0475 int subdet(0), zside, layer, wafer, celltype, cell, wafer2(0), cell2(0);
0476 if (hgcCons_[idet]->waferHexagon8()) {
0477 HGCSiliconDetId detId = HGCSiliconDetId(id);
0478 subdet = (int)(detId.det());
0479 cell = detId.cellU();
0480 cell2 = detId.cellV();
0481 wafer = detId.waferU();
0482 wafer2 = detId.waferV();
0483 celltype = detId.type();
0484 layer = detId.layer();
0485 zside = detId.zside();
0486 xy = hgcCons_[idet]->locateCell(zside, layer, wafer, wafer2, cell, cell2, false, true, false, false, false);
0487 } else if (hgcCons_[idet]->tileTrapezoid()) {
0488 HGCScintillatorDetId detId = HGCScintillatorDetId(id);
0489 subdet = (int)(detId.det());
0490 cell = detId.ietaAbs();
0491 wafer = detId.iphi();
0492 celltype = detId.type();
0493 layer = detId.layer();
0494 zside = detId.zside();
0495 xy = hgcCons_[idet]->locateCellTrap(zside, layer, wafer, cell, false, false);
0496 } else {
0497 HGCalTestNumbering::unpackHexagonIndex(simHit.id(), subdet, zside, layer, wafer, celltype, cell);
0498 xy = hgcCons_[idet]->locateCell(cell, layer, wafer, false);
0499
0500
0501 std::pair<int, int> recoLayerCell = hgcCons_[idet]->simToReco(cell, layer, wafer, hTopo.detectorType());
0502 ok = !(recoLayerCell.second < 0 || recoLayerCell.first < 0);
0503 id = HGCalDetId((ForwardSubdetector)(subdet), zside, layer, celltype, wafer, cell).rawId();
0504 }
0505
0506 edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id << std::dec << " (" << subdet << ":" << zside << ":"
0507 << layer << ":" << celltype << ":" << wafer << ":" << wafer2 << ":" << cell << ":"
0508 << cell2 << ") Flag " << ok;
0509
0510 if (ok) {
0511 float zp = hgcCons_[idet]->waferZ(layer, false);
0512 if (zside < 0)
0513 zp = -zp;
0514 float xp = (zside < 0) ? -xy.first / 10 : xy.first / 10;
0515 float yp = xy.second / 10.0;
0516 float energy = simHit.energy();
0517
0518 float energySum(energy);
0519 if (hitRefs.count(id) != 0)
0520 energySum += std::get<0>(hitRefs[id]);
0521 hitRefs[id] = std::make_tuple(energySum, xp, yp, zp);
0522 if (hist != nullptr)
0523 hist->Fill(energy);
0524 edm::LogVerbatim("HGCalValid") << "Position (" << xp << ", " << yp << ", " << zp << ") "
0525 << " Energy " << simHit.energy() << ":" << energySum;
0526 }
0527 }
0528 }
0529
0530 template <class T1>
0531 void HGCHitValidation::analyzeHGCalRecHit(T1 const &theHits, std::map<unsigned int, HGCHitTuple> const &hitRefs) {
0532 for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0533 DetId id = it->id();
0534 double energy = it->energy();
0535 if (!makeTree_)
0536 hebEnRec_->Fill(energy);
0537 GlobalPoint xyz = hgcGeometry_[2]->getPosition(id);
0538
0539 std::map<unsigned int, HGCHitTuple>::const_iterator itr = hitRefs.find(id.rawId());
0540 if (itr != hitRefs.end()) {
0541 float ang3 = xyz.phi().value();
0542 double fac = sinh(std::get<1>(itr->second));
0543 double pT = std::get<3>(itr->second) / fac;
0544 double xp = pT * cos(std::get<2>(itr->second));
0545 double yp = pT * sin(std::get<2>(itr->second));
0546 if (makeTree_) {
0547 hebRecX_->push_back(xyz.x());
0548 hebRecY_->push_back(xyz.y());
0549 hebRecZ_->push_back(xyz.z());
0550 hebRecEnergy_->push_back(energy);
0551 hebSimX_->push_back(xp);
0552 hebSimY_->push_back(yp);
0553 hebSimZ_->push_back(std::get<3>(itr->second));
0554 hebSimEnergy_->push_back(std::get<0>(itr->second));
0555 hebSimEta_->push_back(std::get<1>(itr->second));
0556 hebRecEta_->push_back(xyz.eta());
0557 hebSimPhi_->push_back(std::get<2>(itr->second));
0558 hebRecPhi_->push_back(ang3);
0559 hebDetID_->push_back(itr->first);
0560 } else {
0561 hebRecVsSimX_->Fill(xp, xyz.x());
0562 hebRecVsSimY_->Fill(yp, xyz.y());
0563 hebRecVsSimZ_->Fill(std::get<3>(itr->second), xyz.z());
0564 hebdEtaVsEta_->Fill(std::get<1>(itr->second), (xyz.eta() - std::get<1>(itr->second)));
0565 hebdPhiVsPhi_->Fill(std::get<2>(itr->second), (ang3 - std::get<2>(itr->second)));
0566 hebdzVsZ_->Fill(std::get<3>(itr->second), (xyz.z() - std::get<3>(itr->second)));
0567 hebEnSimRec_->Fill(std::get<0>(itr->second), energy);
0568 }
0569 edm::LogVerbatim("HGCalValid") << "BHHit: " << std::hex << id.rawId() << std::dec << " Sim ("
0570 << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ", "
0571 << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ") Rec ("
0572 << energy << ", " << xyz.eta() << ", " << ang3 << ", " << xyz.z() << ")";
0573 }
0574 }
0575 }
0576
0577
0578 DEFINE_FWK_MODULE(HGCHitValidation);