File indexing completed on 2023-10-25 10:06:38
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0016 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0017 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0018 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0019 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0020
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0023 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0024 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/Exception.h"
0033 #include "FWCore/Utilities/interface/EDGetToken.h"
0034 #include "FWCore/Utilities/interface/transform.h"
0035 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0036 #include "DQMServices/Core/interface/DQMStore.h"
0037
0038 #include "SimG4CMS/Calo/interface/HGCNumberingScheme.h"
0039 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0040 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0041 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0042
0043 #include <cmath>
0044 #include <memory>
0045 #include <iostream>
0046 #include <string>
0047 #include <vector>
0048
0049
0050
0051 class HGCalHitValidation : public DQMEDAnalyzer {
0052 public:
0053 explicit HGCalHitValidation(const edm::ParameterSet&);
0054 ~HGCalHitValidation() override = default;
0055 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0056
0057 protected:
0058 typedef std::tuple<float, GlobalPoint> HGCHitTuple;
0059
0060 void dqmBeginRun(edm::Run const&, edm::EventSetup const&) override;
0061 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0062 void analyze(const edm::Event&, const edm::EventSetup&) override;
0063 void analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const& simHits,
0064 int idet,
0065 MonitorElement* hist,
0066 std::map<unsigned int, HGCHitTuple>&);
0067
0068 private:
0069
0070 std::vector<const HGCalDDDConstants*> hgcCons_;
0071 std::vector<const HGCalGeometry*> hgcGeometry_;
0072 const std::vector<std::string> geometrySource_;
0073 const std::vector<int> ietaExcludeBH_;
0074 const edm::EDGetTokenT<std::vector<PCaloHit>> eeSimHitToken_;
0075 const edm::EDGetTokenT<std::vector<PCaloHit>> fhSimHitToken_;
0076 const edm::EDGetTokenT<std::vector<PCaloHit>> bhSimHitToken_;
0077 const edm::EDGetTokenT<HGCeeRecHitCollection> eeRecHitToken_;
0078 const edm::EDGetTokenT<HGChefRecHitCollection> fhRecHitToken_;
0079 const edm::EDGetTokenT<HGChebRecHitCollection> bhRecHitToken_;
0080 const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_ddd_;
0081 const std::vector<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>> tok_geom_;
0082
0083
0084 MonitorElement *heedzVsZ, *heedyVsY, *heedxVsX;
0085 MonitorElement *hefdzVsZ, *hefdyVsY, *hefdxVsX;
0086 MonitorElement *hebdzVsZ, *hebdPhiVsPhi, *hebdEtaVsEta;
0087
0088 MonitorElement *heeRecVsSimZ, *heeRecVsSimY, *heeRecVsSimX;
0089 MonitorElement *hefRecVsSimZ, *hefRecVsSimY, *hefRecVsSimX;
0090 MonitorElement *hebRecVsSimZ, *hebRecVsSimY, *hebRecVsSimX;
0091 MonitorElement *heeEnSimRec, *hefEnSimRec, *hebEnSimRec;
0092
0093 MonitorElement *hebEnRec, *hebEnSim;
0094 MonitorElement *hefEnRec, *hefEnSim;
0095 MonitorElement *heeEnRec, *heeEnSim;
0096 };
0097
0098 HGCalHitValidation::HGCalHitValidation(const edm::ParameterSet& cfg)
0099 : geometrySource_(cfg.getParameter<std::vector<std::string>>("geometrySource")),
0100 ietaExcludeBH_(cfg.getParameter<std::vector<int>>("ietaExcludeBH")),
0101 eeSimHitToken_(consumes<std::vector<PCaloHit>>(cfg.getParameter<edm::InputTag>("eeSimHitSource"))),
0102 fhSimHitToken_(consumes<std::vector<PCaloHit>>(cfg.getParameter<edm::InputTag>("fhSimHitSource"))),
0103 bhSimHitToken_(consumes<std::vector<PCaloHit>>(cfg.getParameter<edm::InputTag>("bhSimHitSource"))),
0104 eeRecHitToken_(consumes<HGCeeRecHitCollection>(cfg.getParameter<edm::InputTag>("eeRecHitSource"))),
0105 fhRecHitToken_(consumes<HGChefRecHitCollection>(cfg.getParameter<edm::InputTag>("fhRecHitSource"))),
0106 bhRecHitToken_(consumes<HGChebRecHitCollection>(cfg.getParameter<edm::InputTag>("bhRecHitSource"))),
0107 tok_ddd_{
0108 edm::vector_transform(geometrySource_,
0109 [this](const std::string& name) {
0110 return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0111 edm::ESInputTag{"", name});
0112 })},
0113 tok_geom_{edm::vector_transform(geometrySource_, [this](const std::string& name) {
0114 return esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0115 })} {
0116 #ifdef EDM_ML_DEBUG
0117 edm::LogInfo("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots";
0118 for (unsigned int k = 0; k < ietaExcludeBH_.size(); ++k)
0119 edm::LogInfo("HGCalValid") << " [" << k << "] " << ietaExcludeBH_[k];
0120 #endif
0121 }
0122
0123 void HGCalHitValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0124 edm::ParameterSetDescription desc;
0125 std::vector<std::string> source = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0126 desc.add<std::vector<std::string>>("geometrySource", source);
0127 desc.add<edm::InputTag>("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE"));
0128 desc.add<edm::InputTag>("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront"));
0129 desc.add<edm::InputTag>("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback"));
0130 desc.add<edm::InputTag>("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0131 desc.add<edm::InputTag>("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0132 desc.add<edm::InputTag>("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0133 std::vector<int> dummy;
0134 desc.add<std::vector<int>>("ietaExcludeBH", dummy);
0135 descriptions.add("hgcalHitValidation", desc);
0136 }
0137
0138 void HGCalHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&, edm::EventSetup const&) {
0139 iB.setCurrentFolder("HGCAL/HGCalSimHitsV/HitValidation");
0140
0141
0142 heedzVsZ = iB.book2D("heedzVsZ", "", 720, -360, 360, 100, -0.1, 0.1);
0143 heedyVsY = iB.book2D("heedyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0144 heedxVsX = iB.book2D("heedxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0145 heeRecVsSimZ = iB.book2D("heeRecVsSimZ", "", 720, -360, 360, 720, -360, 360);
0146 heeRecVsSimY = iB.book2D("heeRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0147 heeRecVsSimX = iB.book2D("heeRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0148
0149 hefdzVsZ = iB.book2D("hefdzVsZ", "", 820, -410, 410, 100, -0.1, 0.1);
0150 hefdyVsY = iB.book2D("hefdyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0151 hefdxVsX = iB.book2D("hefdxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0152 hefRecVsSimZ = iB.book2D("hefRecVsSimZ", "", 820, -410, 410, 820, -410, 410);
0153 hefRecVsSimY = iB.book2D("hefRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0154 hefRecVsSimX = iB.book2D("hefRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0155
0156 hebdzVsZ = iB.book2D("hebdzVsZ", "", 1080, -540, 540, 100, -1.0, 1.0);
0157 hebdPhiVsPhi = iB.book2D("hebdPhiVsPhi", "", M_PI * 100, -0.5, M_PI + 0.5, 200, -0.2, 0.2);
0158 hebdEtaVsEta = iB.book2D("hebdEtaVsEta", "", 1000, -5, 5, 200, -0.1, 0.1);
0159 hebRecVsSimZ = iB.book2D("hebRecVsSimZ", "", 1080, -540, 540, 1080, -540, 540);
0160 hebRecVsSimY = iB.book2D("hebRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0161 hebRecVsSimX = iB.book2D("hebRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0162
0163 heeEnRec = iB.book1D("heeEnRec", "", 1000, 0, 10);
0164 heeEnSim = iB.book1D("heeEnSim", "", 1000, 0, 0.01);
0165 heeEnSimRec = iB.book2D("heeEnSimRec", "", 200, 0, 0.002, 200, 0, 0.2);
0166
0167 hefEnRec = iB.book1D("hefEnRec", "", 1000, 0, 10);
0168 hefEnSim = iB.book1D("hefEnSim", "", 1000, 0, 0.01);
0169 hefEnSimRec = iB.book2D("hefEnSimRec", "", 200, 0, 0.001, 200, 0, 0.5);
0170
0171 hebEnRec = iB.book1D("hebEnRec", "", 1000, 0, 15);
0172 hebEnSim = iB.book1D("hebEnSim", "", 1000, 0, 0.01);
0173 hebEnSimRec = iB.book2D("hebEnSimRec", "", 200, 0, 0.02, 200, 0, 4);
0174 }
0175
0176 void HGCalHitValidation::dqmBeginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0177
0178 for (size_t i = 0; i < geometrySource_.size(); i++) {
0179 const edm::ESHandle<HGCalDDDConstants>& hgcCons = iSetup.getHandle(tok_ddd_[i]);
0180 if (hgcCons.isValid()) {
0181 hgcCons_.push_back(hgcCons.product());
0182 } else {
0183 edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl;
0184 }
0185 const edm::ESHandle<HGCalGeometry>& hgcGeom = iSetup.getHandle(tok_geom_[i]);
0186 if (hgcGeom.isValid()) {
0187 hgcGeometry_.push_back(hgcGeom.product());
0188 } else {
0189 edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for " << geometrySource_[i] << std::endl;
0190 }
0191 }
0192 }
0193
0194 void HGCalHitValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0195 std::map<unsigned int, HGCHitTuple> eeHitRefs, fhHitRefs, bhHitRefs;
0196
0197
0198 const edm::Handle<std::vector<PCaloHit>>& eeSimHits = iEvent.getHandle(eeSimHitToken_);
0199
0200 if (eeSimHits.isValid()) {
0201 analyzeHGCalSimHit(eeSimHits, 0, heeEnSim, eeHitRefs);
0202 #ifdef EDM_ML_DEBUG
0203 for (std::map<unsigned int, HGCHitTuple>::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) {
0204 int idx = std::distance(eeHitRefs.begin(), itr);
0205 edm::LogInfo("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0206 << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second).x() << ", "
0207 << std::get<1>(itr->second).y() << ", " << std::get<1>(itr->second).z() << ")";
0208 }
0209 #endif
0210 } else {
0211 edm::LogVerbatim("HGCalValid") << "No EE SimHit Found ";
0212 }
0213
0214
0215 const edm::Handle<std::vector<PCaloHit>>& fhSimHits = iEvent.getHandle(fhSimHitToken_);
0216 if (fhSimHits.isValid()) {
0217 analyzeHGCalSimHit(fhSimHits, 1, hefEnSim, fhHitRefs);
0218 #ifdef EDM_ML_DEBUG
0219 for (std::map<unsigned int, HGCHitTuple>::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) {
0220 int idx = std::distance(fhHitRefs.begin(), itr);
0221 edm::LogInfo("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0222 << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second).x() << ", "
0223 << std::get<1>(itr->second).y() << ", " << std::get<1>(itr->second).z() << ")";
0224 }
0225 #endif
0226 } else {
0227 edm::LogVerbatim("HGCalValid") << "No FH SimHit Found ";
0228 }
0229
0230
0231 const edm::Handle<std::vector<PCaloHit>>& bhSimHits = iEvent.getHandle(bhSimHitToken_);
0232 if (bhSimHits.isValid()) {
0233 analyzeHGCalSimHit(bhSimHits, 2, hebEnSim, bhHitRefs);
0234 #ifdef EDM_ML_DEBUG
0235 for (std::map<unsigned int, HGCHitTuple>::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) {
0236 int idx = std::distance(bhHitRefs.begin(), itr);
0237 edm::LogInfo("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0238 << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second).x() << ", "
0239 << std::get<1>(itr->second).y() << ", " << std::get<1>(itr->second).z() << ")";
0240 }
0241 #endif
0242 } else {
0243 edm::LogVerbatim("HGCalValid") << "No BH SimHit Found ";
0244 }
0245
0246
0247 const edm::Handle<HGCeeRecHitCollection>& eeRecHit = iEvent.getHandle(eeRecHitToken_);
0248 if (eeRecHit.isValid()) {
0249 const HGCeeRecHitCollection* theHits = (eeRecHit.product());
0250 for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0251 double energy = it->energy();
0252 heeEnRec->Fill(energy);
0253 std::map<unsigned int, HGCHitTuple>::const_iterator itr = eeHitRefs.find(it->id().rawId());
0254 if (itr != eeHitRefs.end()) {
0255 GlobalPoint xyz = hgcGeometry_[0]->getPosition(it->id());
0256 heeRecVsSimX->Fill(std::get<1>(itr->second).x(), xyz.x());
0257 heeRecVsSimY->Fill(std::get<1>(itr->second).y(), xyz.y());
0258 heeRecVsSimZ->Fill(std::get<1>(itr->second).z(), xyz.z());
0259 heedxVsX->Fill(std::get<1>(itr->second).x(), (xyz.x() - std::get<1>(itr->second).x()));
0260 heedyVsY->Fill(std::get<1>(itr->second).y(), (xyz.y() - std::get<1>(itr->second).y()));
0261 heedzVsZ->Fill(std::get<1>(itr->second).z(), (xyz.z() - std::get<1>(itr->second).z()));
0262 heeEnSimRec->Fill(std::get<0>(itr->second), energy);
0263 #ifdef EDM_ML_DEBUG
0264 edm::LogInfo("HGCalValid") << "EEHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0265 << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ") Rec ("
0266 << energy << ", " << xyz.x() << ", " << xyz.y() << ", " << xyz.z() << ")";
0267 #endif
0268 }
0269 }
0270 } else {
0271 edm::LogVerbatim("HGCalValid") << "No EE RecHit Found ";
0272 }
0273
0274
0275 const edm::Handle<HGChefRecHitCollection>& fhRecHit = iEvent.getHandle(fhRecHitToken_);
0276 if (fhRecHit.isValid()) {
0277 const HGChefRecHitCollection* theHits = (fhRecHit.product());
0278 for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0279 double energy = it->energy();
0280 hefEnRec->Fill(energy);
0281 std::map<unsigned int, HGCHitTuple>::const_iterator itr = fhHitRefs.find(it->id().rawId());
0282 if (itr != fhHitRefs.end()) {
0283 GlobalPoint xyz = hgcGeometry_[1]->getPosition(it->id());
0284
0285 hefRecVsSimX->Fill(std::get<1>(itr->second).x(), xyz.x());
0286 hefRecVsSimY->Fill(std::get<1>(itr->second).y(), xyz.y());
0287 hefRecVsSimZ->Fill(std::get<1>(itr->second).z(), xyz.z());
0288 hefdxVsX->Fill(std::get<1>(itr->second).x(), (xyz.x() - std::get<1>(itr->second).x()));
0289 hefdyVsY->Fill(std::get<1>(itr->second).y(), (xyz.y() - std::get<1>(itr->second).y()));
0290 hefdzVsZ->Fill(std::get<1>(itr->second).z(), (xyz.z() - std::get<1>(itr->second).z()));
0291 hefEnSimRec->Fill(std::get<0>(itr->second), energy);
0292 #ifdef EDM_ML_DEBUG
0293 edm::LogInfo("HGCalValid") << "FHHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0294 << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ") Rec ("
0295 << energy << "," << xyz.x() << ", " << xyz.y() << ", " << xyz.z() << ")";
0296 #endif
0297 }
0298 }
0299 } else {
0300 edm::LogVerbatim("HGCalValid") << "No FH RecHit Found ";
0301 }
0302
0303
0304 const edm::Handle<HGChebRecHitCollection>& bhRecHit = iEvent.getHandle(bhRecHitToken_);
0305 if (bhRecHit.isValid()) {
0306 const HGChebRecHitCollection* theHits = (bhRecHit.product());
0307 for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0308 double energy = it->energy();
0309 hebEnRec->Fill(energy);
0310 std::map<unsigned int, HGCHitTuple>::const_iterator itr = bhHitRefs.find(it->id().rawId());
0311 GlobalPoint xyz = hgcGeometry_[2]->getPosition(it->id());
0312 if (itr != bhHitRefs.end()) {
0313 float ang3 = xyz.phi().value();
0314 float ang3s = std::get<1>(itr->second).phi().value();
0315 hebRecVsSimX->Fill(std::get<1>(itr->second).x(), xyz.x());
0316 hebRecVsSimY->Fill(std::get<1>(itr->second).y(), xyz.y());
0317 hebRecVsSimZ->Fill(std::get<1>(itr->second).z(), xyz.z());
0318 hebdEtaVsEta->Fill(std::get<1>(itr->second).eta(), (xyz.eta() - std::get<1>(itr->second).eta()));
0319 hebdPhiVsPhi->Fill(std::get<1>(itr->second).phi(), (ang3 - ang3s));
0320 hebdzVsZ->Fill(std::get<1>(itr->second).z(), (xyz.z() - std::get<1>(itr->second).z()));
0321 hebEnSimRec->Fill(std::get<0>(itr->second), energy);
0322
0323 #ifdef EDM_ML_DEBUG
0324 edm::LogInfo("HGCalValid") << "BHHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0325 << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ") Rec ("
0326 << energy << ", " << xyz.x() << ", " << xyz.y() << ", " << xyz.z() << ")\n";
0327 #endif
0328 }
0329 }
0330 } else {
0331 edm::LogVerbatim("HGCalValid") << "No BH RecHit Found ";
0332 }
0333 }
0334
0335 void HGCalHitValidation::analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const& simHits,
0336 int idet,
0337 MonitorElement* hist,
0338 std::map<unsigned int, HGCHitTuple>& hitRefs) {
0339 for (std::vector<PCaloHit>::const_iterator simHit = simHits->begin(); simHit != simHits->end(); ++simHit) {
0340 DetId id(simHit->id());
0341 GlobalPoint p = hgcGeometry_[idet]->getPosition(id);
0342 float energy = simHit->energy();
0343
0344 float energySum(energy);
0345 if (hitRefs.count(id.rawId()) != 0)
0346 energySum += std::get<0>(hitRefs[id.rawId()]);
0347 hitRefs[id.rawId()] = std::make_tuple(energySum, p);
0348 hist->Fill(energy);
0349 }
0350 }
0351
0352
0353 #include "FWCore/Framework/interface/MakerMacros.h"
0354 DEFINE_FWK_MODULE(HGCalHitValidation);
