From cd0cc38a40b6e4349721da1ef48a55921b0a2b73 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Thu, 22 Sep 2022 05:42:02 +0200 Subject: [PATCH 1/8] Add a validation code to look into missed RecHits --- .../HGCalValidation/test/HGCMissingRecHit.cc | 357 ++++++++++++++++++ .../test/python/runHGCMissingRecHit_cfg.py | 65 ++++ 2 files changed, 422 insertions(+) create mode 100644 Validation/HGCalValidation/test/HGCMissingRecHit.cc create mode 100644 Validation/HGCalValidation/test/python/runHGCMissingRecHit_cfg.py diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc new file mode 100644 index 0000000000000..c4b9df0d69eb0 --- /dev/null +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -0,0 +1,357 @@ +/// -*- C++ -*- +// +// Package: HGCMissingRecHit +// Class: HGCMissingRecHit +// +/**\class HGCMissingRecHit HGCMissingRecHit.cc Validation/HGCalValidation/test/HGCMissingRecHit.cc + + Description: [one line class summary] + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: "Sunanda Banerjee" +// Created: Tue September 20 17:55:26 CDT 2022 +// $Id$ +// +// + +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h" +#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHit.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/transform.h" + +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" +#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" + +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" + +#include +#include +#include +#include +#include +#include + +#include + +class HGCMissingRecHit : public edm::one::EDAnalyzer { +public: + explicit HGCMissingRecHit(const edm::ParameterSet &); + ~HGCMissingRecHit() override = default; + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + typedef std::tuple HGCHitTuple; + + void beginJob() override; + void endJob() override {} + void beginRun(edm::Run const &, edm::EventSetup const &) override; + void analyze(edm::Event const &, edm::EventSetup const &) override; + void endRun(edm::Run const &, edm::EventSetup const &) override {} + virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {} + virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {} + void analyzeHGCalSimHit(edm::Handle> const &simHits, + int idet, + std::map &); + template + void analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs); + +private: + //HGC Geometry + const std::vector geometrySource_, detectors_; + const std::vector ietaExcludeBH_; + const std::vector> tok_hgcal_; + const std::vector> tok_hgcalg_; + std::vector hgcCons_; + std::vector hgcGeometry_; + + const edm::InputTag eeSimHitSource, fhSimHitSource, bhSimHitSource; + const edm::EDGetTokenT> eeSimHitToken_; + const edm::EDGetTokenT> fhSimHitToken_; + const edm::EDGetTokenT> bhSimHitToken_; + const edm::InputTag eeRecHitSource, fhRecHitSource, bhRecHitSource; + const edm::EDGetTokenT eeRecHitToken_; + const edm::EDGetTokenT fhRecHitToken_; + const edm::EDGetTokenT bhRecHitToken_; + + std::vector goodHitsE_, missedHitsE_; +}; + +HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) + : geometrySource_(cfg.getParameter>("geometrySource")), + detectors_(cfg.getParameter>("detectors")), + ietaExcludeBH_(cfg.getParameter>("ietaExcludeBH")), + tok_hgcal_{ + edm::vector_transform(geometrySource_, + [this](const std::string &name) { + return esConsumes( + edm::ESInputTag{"", name}); + })}, + tok_hgcalg_{edm::vector_transform( + geometrySource_, + [this](const std::string &name) { + return esConsumes(edm::ESInputTag{"", name}); + })}, + eeSimHitSource(cfg.getParameter("eeSimHitSource")), + fhSimHitSource(cfg.getParameter("fhSimHitSource")), + bhSimHitSource(cfg.getParameter("bhSimHitSource")), + eeSimHitToken_(consumes>(eeSimHitSource)), + fhSimHitToken_(consumes>(fhSimHitSource)), + bhSimHitToken_(consumes>(bhSimHitSource)), + eeRecHitSource(cfg.getParameter("eeRecHitSource")), + fhRecHitSource(cfg.getParameter("fhRecHitSource")), + bhRecHitSource(cfg.getParameter("bhRecHitSource")), + eeRecHitToken_(consumes(eeRecHitSource)), + fhRecHitToken_(consumes(fhRecHitSource)), + bhRecHitToken_(consumes(bhRecHitSource)) { + usesResource(TFileService::kSharedResource); + + edm::LogVerbatim("HGCalValid") << "Use " << geometrySource_.size() << " Geometry sources"; + for (unsigned int k = 0; k < geometrySource_.size(); k++) + edm::LogVerbatim("HGCalValid") << " " << detectors_[k] << ":" << geometrySource_[k]; + edm::LogVerbatim("HGCalValid") << "SimHit labels: " << eeSimHitSource << " " << fhSimHitSource << " " + << bhSimHitSource; + edm::LogVerbatim("HGCalValid") << "RecHit labels: " << eeRecHitSource << " " << fhRecHitSource << " " + << bhRecHitSource; + edm::LogVerbatim("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots"; + for (unsigned int k = 0; k < ietaExcludeBH_.size(); ++k) + edm::LogVerbatim("HGCalValid") << " [" << k << "] " << ietaExcludeBH_[k]; +} + +void HGCMissingRecHit::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + std::vector sources = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"}; + std::vector names = {"EE", "HE Silicon", "HE Scintillator"}; + std::vector etas; + edm::ParameterSetDescription desc; + desc.add>("geometrySource", sources); + desc.add>("detectors", names); + desc.add("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE")); + desc.add("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront")); + desc.add("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback")); + desc.add("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits")); + desc.add("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits")); + desc.add("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); + desc.add>("ietaExcludeBH", etas); + descriptions.add("hgcMissingRecHit", desc); +} + +void HGCMissingRecHit::beginJob() { + //initiating fileservice + edm::Service fs; + for (unsigned int k = 0; k < detectors_.size(); ++k) { + char name[50], title[100]; + sprintf(name, "GoodE%s", geometrySource_[k].c_str()); + sprintf(title, "SimHit energy for good hits in %s", detectors_[k].c_str()); + goodHitsE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); + goodHitsE_.back()->Sumw2(); + sprintf(name, "MissE%s", geometrySource_[k].c_str()); + sprintf(title, "SimHit energy for missed hits in %s", detectors_[k].c_str()); + missedHitsE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); + missedHitsE_.back()->Sumw2(); + } +} + +void HGCMissingRecHit::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) { + //initiating hgc Geometry + for (size_t i = 0; i < geometrySource_.size(); i++) { + edm::LogVerbatim("HGCalValid") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i; + const edm::ESHandle& hgcCons = iSetup.getHandle(tok_hgcal_[i]); + if (hgcCons.isValid()) { + hgcCons_.push_back(hgcCons.product()); + } else { + edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl; + } + const edm::ESHandle& hgcGeom = iSetup.getHandle(tok_hgcalg_[i]); + if (hgcGeom.isValid()) { + hgcGeometry_.push_back(hgcGeom.product()); + } else { + edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for " << geometrySource_[i] << std::endl; + } + } +} + +void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) { + std::map eeHitRefs, fhHitRefs, bhHitRefs; + + //Accesing ee simhits + const edm::Handle> &eeSimHits = iEvent.getHandle(eeSimHitToken_); + + if (eeSimHits.isValid()) { + analyzeHGCalSimHit(eeSimHits, 0, eeHitRefs); + for (std::map::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) { + int idx = std::distance(eeHitRefs.begin(), itr); + edm::LogVerbatim("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " + << std::get<0>(itr->second) << "; Position" + << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", " + << std::get<3>(itr->second) << ")"; + } + } else { + edm::LogWarning("HGCalValid") << "No EE SimHit Found " << std::endl; + } + + //Accesing fh simhits + const edm::Handle> &fhSimHits = iEvent.getHandle(fhSimHitToken_); + if (fhSimHits.isValid()) { + analyzeHGCalSimHit(fhSimHits, 1, fhHitRefs); + for (std::map::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) { + int idx = std::distance(fhHitRefs.begin(), itr); + edm::LogVerbatim("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " + << std::get<0>(itr->second) << "; Position" + << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", " + << std::get<3>(itr->second) << ")"; + } + } else { + edm::LogWarning("HGCalValid") << "No FH SimHit Found " << std::endl; + } + + //Accessing bh simhits + const edm::Handle> &bhSimHits = iEvent.getHandle(bhSimHitToken_); + if (bhSimHits.isValid()) { + analyzeHGCalSimHit(bhSimHits, 2, bhHitRefs); + for (std::map::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) { + int idx = std::distance(bhHitRefs.begin(), itr); + edm::LogVerbatim("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " + << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second) << ", " + << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ")"; + } + } else { + edm::LogWarning("HGCalValid") << "No BH SimHit Found " << std::endl; + } + + //accessing EE Rechit information + const edm::Handle &eeRecHit = iEvent.getHandle(eeRecHitToken_); + if (eeRecHit.isValid()) { + const HGCeeRecHitCollection *theHits = (eeRecHit.product()); + analyzeHGCalRecHit(theHits, 0, eeHitRefs); + } else { + edm::LogWarning("HGCalValid") << "No EE RecHit Found " << std::endl; + } + + //accessing FH Rechit information + const edm::Handle &fhRecHit = iEvent.getHandle(fhRecHitToken_); + if (fhRecHit.isValid()) { + const HGChefRecHitCollection *theHits = (fhRecHit.product()); + analyzeHGCalRecHit(theHits, 1, fhHitRefs); + } else { + edm::LogWarning("HGCalValid") << "No FH RecHit Found " << std::endl; + } + + //accessing BH Rechit information + const edm::Handle &bhRecHit = iEvent.getHandle(bhRecHitToken_); + if (bhRecHit.isValid()) { + const HGChebRecHitCollection *theHits = (bhRecHit.product()); + analyzeHGCalRecHit(theHits, 2, bhHitRefs); + } else { + edm::LogWarning("HGCalValid") << "No BH RecHit Found " << std::endl; + } +} + +void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> const &simHits, + int idet, + std::map &hitRefs) { + const bool debug(false); + for (auto const &simHit : *simHits) { + unsigned int id = simHit.id(); + std::pair xy; + bool ok(true); + int subdet(0), zside(0), layer(0), wafer(0), celltype(0), cell(0), wafer2(0), cell2(0); + if ((hgcCons_[idet]->waferHexagon8()) && ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi))) { + HGCSiliconDetId detId = HGCSiliconDetId(id); + subdet = static_cast(detId.det()); + cell = detId.cellU(); + cell2 = detId.cellV(); + wafer = detId.waferU(); + wafer2 = detId.waferV(); + celltype = detId.type(); + layer = detId.layer(); + zside = detId.zside(); + xy = hgcCons_[idet]->locateCell(layer, wafer, wafer2, cell, cell2, true, true, debug); + } else if ((hgcCons_[idet]->tileTrapezoid()) && (DetId(id).det() == DetId::HGCalHSc)) { + HGCScintillatorDetId detId = HGCScintillatorDetId(id); + subdet = static_cast(detId.det()); + wafer = detId.ietaAbs(); + cell = detId.iphi(); + celltype = detId.type(); + layer = detId.layer(); + zside = detId.zside(); + xy = hgcCons_[idet]->locateCellTrap(layer, wafer, cell, false, debug); + edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " LocateCellTrap i/p " << layer << ":" << wafer << ":" << cell << " o/p " << xy.first << ":" << xy.second; + } else { + // This is an invalid cell + ok = false; + std::ostringstream st1; + if (DetId(id).det() == DetId::HGCalHSc) { + st1 << HGCScintillatorDetId(id); + } else if ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi)) { + st1 << HGCSiliconDetId(id); + } else { + st1 << "Not a Standard One"; + } + edm::LogVerbatim("HGCalError") << "Hit " << std::hex << id << std::dec << " " << st1.str() << " in the wrong collection for detector " << idet << ":" << geometrySource_[idet] << " ***** ERROR *****"; + } + + edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id << std::dec << " (" << subdet << ":" << zside << ":" + << layer << ":" << celltype << ":" << wafer << ":" << wafer2 << ":" << cell << ":" + << cell2 << ") Flag " << ok; + + if (ok) { + float zp = hgcCons_[idet]->waferZ(layer, false); + if (zside < 0) + zp = -zp; + float xp = (zside < 0) ? -xy.first / 10 : xy.first / 10; + float yp = xy.second / 10.0; + float energy = simHit.energy(); + + float energySum(energy); + if (hitRefs.count(id) != 0) + energySum += std::get<0>(hitRefs[id]); + hitRefs[id] = std::make_tuple(energySum, xp, yp, zp); + edm::LogVerbatim("HGCalValid") << "Position (" << xp << ", " << yp << ", " << zp << ") " + << " Energy " << simHit.energy() << ":" << energySum; + } + } +} + +template +void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs) { + std::vector ids; + for (auto it = theHits->begin(); it != theHits->end(); ++it) + ids.emplace_back((it->id().rawId())); + for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { + auto itr = std::find(ids.begin(), ids.end(), it->first); + if (itr == ids.end()) { + missedHitsE_[idet]->Fill(std::get<0>(it->second)); + std::ostringstream st1; + if (DetId(it->first).det() == DetId::HGCalHSc) + st1 << HGCScintillatorDetId(it->first); + else + st1 << HGCSiliconDetId(it->first); + edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() << " SimHit (E = " << std::get<0>(it->second) << ", X = " << std::get<1>(it->second) << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) << ") is missing in the RecHit collection"; + } else { + goodHitsE_[idet]->Fill(std::get<0>(it->second)); + } + } +} + +//define this as a plug-in +DEFINE_FWK_MODULE(HGCMissingRecHit); diff --git a/Validation/HGCalValidation/test/python/runHGCMissingRecHit_cfg.py b/Validation/HGCalValidation/test/python/runHGCMissingRecHit_cfg.py new file mode 100644 index 0000000000000..f583857902ed6 --- /dev/null +++ b/Validation/HGCalValidation/test/python/runHGCMissingRecHit_cfg.py @@ -0,0 +1,65 @@ +############################################################################### +# Way to use this: +# cmsRun runHGCMissingRecHit_cfg.py geometry=D88 +# +# Options for geometry D88, D92, D93 +# +############################################################################### +import FWCore.ParameterSet.Config as cms +import os, sys, imp, re +import FWCore.ParameterSet.VarParsing as VarParsing + +#################################################################### +### SETUP OPTIONS +options = VarParsing.VarParsing('standard') +options.register('geometry', + "D88", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "geometry of operations: D88, D92, D93") + +### get and parse the command line arguments +options.parseArguments() + +print(options) + +#################################################################### +# Use the options +from Configuration.Eras.Era_Phase2C11M9_cff import Phase2C11M9 +process = cms.Process('HGCHitAnalysis',Phase2C11M9) + +geomFile = "Configuration.Geometry.GeometryExtended2026" + options.geometry + "Reco_cff" +inFile = "file:step3" + options.geometry + ".root" +outFile = "missedRecHit" + options.geometry + ".root" + +print("Geometry file: ", geomFile) +print("Input file: ", inFile) +print("Output file: ", outFile) + +process.load(geomFile) +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic_T21', '') +process.MessageLogger.cerr.FwkReport.reportEvery = 1 +if hasattr(process,'MessageLogger'): + process.MessageLogger.HGCalMiss=dict() + process.MessageLogger.HGCalError=dict() +# process.MessageLogger.HGCalGeom=dict() + +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring(inFile) +) + +process.load('Validation.HGCalValidation.hgcMissingRecHit_cfi') +process.TFileService = cms.Service("TFileService", + fileName = cms.string(outFile)) + +process.p = cms.Path(process.hgcMissingRecHit) + + From cd54e22c17dab47e554bf323c53d24d7d5615da8 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Thu, 22 Sep 2022 06:16:27 +0200 Subject: [PATCH 2/8] Code check --- .../HGCalValidation/test/HGCMissingRecHit.cc | 39 ++++++++++++------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index c4b9df0d69eb0..a1f817811f3e5 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -94,7 +94,7 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer fhRecHitToken_; const edm::EDGetTokenT bhRecHitToken_; - std::vector goodHitsE_, missedHitsE_; + std::vector goodHitsE_, missedHitsE_; }; HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) @@ -175,13 +175,13 @@ void HGCMissingRecHit::beginRun(edm::Run const &iRun, edm::EventSetup const &iSe //initiating hgc Geometry for (size_t i = 0; i < geometrySource_.size(); i++) { edm::LogVerbatim("HGCalValid") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i; - const edm::ESHandle& hgcCons = iSetup.getHandle(tok_hgcal_[i]); + const edm::ESHandle &hgcCons = iSetup.getHandle(tok_hgcal_[i]); if (hgcCons.isValid()) { hgcCons_.push_back(hgcCons.product()); } else { edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl; } - const edm::ESHandle& hgcGeom = iSetup.getHandle(tok_hgcalg_[i]); + const edm::ESHandle &hgcGeom = iSetup.getHandle(tok_hgcalg_[i]); if (hgcGeom.isValid()) { hgcGeometry_.push_back(hgcGeom.product()); } else { @@ -275,7 +275,8 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con std::pair xy; bool ok(true); int subdet(0), zside(0), layer(0), wafer(0), celltype(0), cell(0), wafer2(0), cell2(0); - if ((hgcCons_[idet]->waferHexagon8()) && ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi))) { + if ((hgcCons_[idet]->waferHexagon8()) && + ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi))) { HGCSiliconDetId detId = HGCSiliconDetId(id); subdet = static_cast(detId.det()); cell = detId.cellU(); @@ -295,20 +296,23 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con layer = detId.layer(); zside = detId.zside(); xy = hgcCons_[idet]->locateCellTrap(layer, wafer, cell, false, debug); - edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " LocateCellTrap i/p " << layer << ":" << wafer << ":" << cell << " o/p " << xy.first << ":" << xy.second; + edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " LocateCellTrap i/p " << layer << ":" + << wafer << ":" << cell << " o/p " << xy.first << ":" << xy.second; } else { // This is an invalid cell ok = false; std::ostringstream st1; if (DetId(id).det() == DetId::HGCalHSc) { - st1 << HGCScintillatorDetId(id); + st1 << HGCScintillatorDetId(id); } else if ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi)) { - st1 << HGCSiliconDetId(id); + st1 << HGCSiliconDetId(id); } else { - st1 << "Not a Standard One"; + st1 << "Not a Standard One"; } - edm::LogVerbatim("HGCalError") << "Hit " << std::hex << id << std::dec << " " << st1.str() << " in the wrong collection for detector " << idet << ":" << geometrySource_[idet] << " ***** ERROR *****"; - } + edm::LogVerbatim("HGCalError") << "Hit " << std::hex << id << std::dec << " " << st1.str() + << " in the wrong collection for detector " << idet << ":" << geometrySource_[idet] + << " ***** ERROR *****"; + } edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id << std::dec << " (" << subdet << ":" << zside << ":" << layer << ":" << celltype << ":" << wafer << ":" << wafer2 << ":" << cell << ":" @@ -333,9 +337,11 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con } template -void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs) { +void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, + int idet, + std::map const &hitRefs) { std::vector ids; - for (auto it = theHits->begin(); it != theHits->end(); ++it) + for (auto it = theHits->begin(); it != theHits->end(); ++it) ids.emplace_back((it->id().rawId())); for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { auto itr = std::find(ids.begin(), ids.end(), it->first); @@ -343,10 +349,13 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, int idet, std::map< missedHitsE_[idet]->Fill(std::get<0>(it->second)); std::ostringstream st1; if (DetId(it->first).det() == DetId::HGCalHSc) - st1 << HGCScintillatorDetId(it->first); + st1 << HGCScintillatorDetId(it->first); else - st1 << HGCSiliconDetId(it->first); - edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() << " SimHit (E = " << std::get<0>(it->second) << ", X = " << std::get<1>(it->second) << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) << ") is missing in the RecHit collection"; + st1 << HGCSiliconDetId(it->first); + edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() + << " SimHit (E = " << std::get<0>(it->second) << ", X = " << std::get<1>(it->second) + << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) + << ") is missing in the RecHit collection"; } else { goodHitsE_[idet]->Fill(std::get<0>(it->second)); } From bf9f354c521fde102970372a36087c719d2f0318 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Thu, 22 Sep 2022 14:28:38 +0200 Subject: [PATCH 3/8] Extend for comparison with Digis as well --- .../HGCalValidation/test/HGCMissingRecHit.cc | 105 +++++++++++++++--- 1 file changed, 92 insertions(+), 13 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index a1f817811f3e5..a6b067717d9a9 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -23,6 +23,7 @@ #include "DataFormats/DetId/interface/DetId.h" #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h" #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h" +#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h" #include "DataFormats/HGCRecHit/interface/HGCRecHit.h" #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" @@ -74,6 +75,8 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer &); template + void analyzeHGCalDigi(T1 const &theHits, int idet, std::map const &hitRefs); + template void analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs); private: @@ -89,12 +92,17 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer> eeSimHitToken_; const edm::EDGetTokenT> fhSimHitToken_; const edm::EDGetTokenT> bhSimHitToken_; + const edm::InputTag eeDigiSource, fhDigiSource, bhDigiSource; + const edm::EDGetTokenT eeDigiToken_; + const edm::EDGetTokenT fhDigiToken_; + const edm::EDGetTokenT bhDigiToken_; const edm::InputTag eeRecHitSource, fhRecHitSource, bhRecHitSource; const edm::EDGetTokenT eeRecHitToken_; const edm::EDGetTokenT fhRecHitToken_; const edm::EDGetTokenT bhRecHitToken_; - std::vector goodHitsE_, missedHitsE_; + std::vector goodHitsDE_, missedHitsDE_; + std::vector goodHitsRE_, missedHitsRE_; }; HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) @@ -118,6 +126,12 @@ HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) eeSimHitToken_(consumes>(eeSimHitSource)), fhSimHitToken_(consumes>(fhSimHitSource)), bhSimHitToken_(consumes>(bhSimHitSource)), + eeDigiSource(cfg.getParameter("eeDigiSource")), + fhDigiSource(cfg.getParameter("fhDigiSource")), + bhDigiSource(cfg.getParameter("bhDigiSource")), + eeDigiToken_(consumes(eeDigiSource)), + fhDigiToken_(consumes(fhDigiSource)), + bhDigiToken_(consumes(bhDigiSource)), eeRecHitSource(cfg.getParameter("eeRecHitSource")), fhRecHitSource(cfg.getParameter("fhRecHitSource")), bhRecHitSource(cfg.getParameter("bhRecHitSource")), @@ -131,6 +145,8 @@ HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) edm::LogVerbatim("HGCalValid") << " " << detectors_[k] << ":" << geometrySource_[k]; edm::LogVerbatim("HGCalValid") << "SimHit labels: " << eeSimHitSource << " " << fhSimHitSource << " " << bhSimHitSource; + edm::LogVerbatim("HGCalValid") << "Digi labels: " << eeDigiSource << " " << fhDigiSource << " " + << bhDigiSource; edm::LogVerbatim("HGCalValid") << "RecHit labels: " << eeRecHitSource << " " << fhRecHitSource << " " << bhRecHitSource; edm::LogVerbatim("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots"; @@ -148,9 +164,12 @@ void HGCMissingRecHit::fillDescriptions(edm::ConfigurationDescriptions &descript desc.add("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE")); desc.add("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront")); desc.add("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback")); + desc.add("eeDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "EE")); + desc.add("fhDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "HEfront")); + desc.add("bhDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "HEback")); desc.add("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits")); desc.add("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits")); - desc.add("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); + desc.add("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); desc.add>("ietaExcludeBH", etas); descriptions.add("hgcMissingRecHit", desc); } @@ -160,14 +179,22 @@ void HGCMissingRecHit::beginJob() { edm::Service fs; for (unsigned int k = 0; k < detectors_.size(); ++k) { char name[50], title[100]; - sprintf(name, "GoodE%s", geometrySource_[k].c_str()); - sprintf(title, "SimHit energy for good hits in %s", detectors_[k].c_str()); - goodHitsE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); - goodHitsE_.back()->Sumw2(); - sprintf(name, "MissE%s", geometrySource_[k].c_str()); - sprintf(title, "SimHit energy for missed hits in %s", detectors_[k].c_str()); - missedHitsE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); - missedHitsE_.back()->Sumw2(); + sprintf(name, "GoodDE%s", geometrySource_[k].c_str()); + sprintf(title, "SimHit energy present among Digis in %s", detectors_[k].c_str()); + goodHitsDE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); + goodHitsDE_.back()->Sumw2(); + sprintf(name, "MissDE%s", geometrySource_[k].c_str()); + sprintf(title, "SimHit energy absent among Digis in %s", detectors_[k].c_str()); + missedHitsDE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); + missedHitsDE_.back()->Sumw2(); + sprintf(name, "GoodRE%s", geometrySource_[k].c_str()); + sprintf(title, "SimHit energy present among RecHits in %s", detectors_[k].c_str()); + goodHitsRE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); + goodHitsRE_.back()->Sumw2(); + sprintf(name, "MissRE%s", geometrySource_[k].c_str()); + sprintf(title, "SimHit energy absent among RecHits in %s", detectors_[k].c_str()); + missedHitsRE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); + missedHitsRE_.back()->Sumw2(); } } @@ -195,7 +222,6 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & //Accesing ee simhits const edm::Handle> &eeSimHits = iEvent.getHandle(eeSimHitToken_); - if (eeSimHits.isValid()) { analyzeHGCalSimHit(eeSimHits, 0, eeHitRefs); for (std::map::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) { @@ -238,6 +264,33 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & edm::LogWarning("HGCalValid") << "No BH SimHit Found " << std::endl; } + //accessing EE Digi information + const edm::Handle &eeDigi = iEvent.getHandle(eeDigiToken_); + if (eeDigi.isValid()) { + const HGCalDigiCollection *theHits = (eeDigi.product()); + analyzeHGCalDigi(theHits, 0, eeHitRefs); + } else { + edm::LogWarning("HGCalValid") << "No EE Digi Found " << std::endl; + } + + //accessing FH Digi information + const edm::Handle &fhDigi = iEvent.getHandle(fhDigiToken_); + if (fhDigi.isValid()) { + const HGCalDigiCollection *theHits = (fhDigi.product()); + analyzeHGCalDigi(theHits, 1, fhHitRefs); + } else { + edm::LogWarning("HGCalValid") << "No FH Digi Found " << std::endl; + } + + //accessing BH Digi information + const edm::Handle &bhDigi = iEvent.getHandle(bhDigiToken_); + if (bhDigi.isValid()) { + const HGCalDigiCollection *theHits = (bhDigi.product()); + analyzeHGCalDigi(theHits, 2, bhHitRefs); + } else { + edm::LogWarning("HGCalValid") << "No BH Digi Found " << std::endl; + } + //accessing EE Rechit information const edm::Handle &eeRecHit = iEvent.getHandle(eeRecHitToken_); if (eeRecHit.isValid()) { @@ -336,6 +389,32 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con } } +template +void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, + int idet, + std::map const &hitRefs) { + std::vector ids; + for (auto it = theHits->begin(); it != theHits->end(); ++it) + ids.emplace_back((it->id().rawId())); + for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { + auto itr = std::find(ids.begin(), ids.end(), it->first); + if (itr == ids.end()) { + missedHitsDE_[idet]->Fill(std::get<0>(it->second)); + std::ostringstream st1; + if (DetId(it->first).det() == DetId::HGCalHSc) + st1 << HGCScintillatorDetId(it->first); + else + st1 << HGCSiliconDetId(it->first); + edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() + << " SimHit (E = " << std::get<0>(it->second) << ", X = " << std::get<1>(it->second) + << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) + << ") is missing in the Digi collection"; + } else { + goodHitsDE_[idet]->Fill(std::get<0>(it->second)); + } + } +} + template void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, int idet, @@ -346,7 +425,7 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { auto itr = std::find(ids.begin(), ids.end(), it->first); if (itr == ids.end()) { - missedHitsE_[idet]->Fill(std::get<0>(it->second)); + missedHitsRE_[idet]->Fill(std::get<0>(it->second)); std::ostringstream st1; if (DetId(it->first).det() == DetId::HGCalHSc) st1 << HGCScintillatorDetId(it->first); @@ -357,7 +436,7 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) << ") is missing in the RecHit collection"; } else { - goodHitsE_[idet]->Fill(std::get<0>(it->second)); + goodHitsRE_[idet]->Fill(std::get<0>(it->second)); } } } From cbad48cbc95a55830b352d9b90cbe068602add17 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Thu, 22 Sep 2022 14:45:20 +0200 Subject: [PATCH 4/8] Code check --- Validation/HGCalValidation/test/HGCMissingRecHit.cc | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index a6b067717d9a9..ab6201697a5f6 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -145,8 +145,7 @@ HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) edm::LogVerbatim("HGCalValid") << " " << detectors_[k] << ":" << geometrySource_[k]; edm::LogVerbatim("HGCalValid") << "SimHit labels: " << eeSimHitSource << " " << fhSimHitSource << " " << bhSimHitSource; - edm::LogVerbatim("HGCalValid") << "Digi labels: " << eeDigiSource << " " << fhDigiSource << " " - << bhDigiSource; + edm::LogVerbatim("HGCalValid") << "Digi labels: " << eeDigiSource << " " << fhDigiSource << " " << bhDigiSource; edm::LogVerbatim("HGCalValid") << "RecHit labels: " << eeRecHitSource << " " << fhRecHitSource << " " << bhRecHitSource; edm::LogVerbatim("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots"; @@ -169,7 +168,7 @@ void HGCMissingRecHit::fillDescriptions(edm::ConfigurationDescriptions &descript desc.add("bhDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "HEback")); desc.add("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits")); desc.add("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits")); - desc.add("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); + desc.add("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits")); desc.add>("ietaExcludeBH", etas); descriptions.add("hgcMissingRecHit", desc); } @@ -391,8 +390,8 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con template void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, - int idet, - std::map const &hitRefs) { + int idet, + std::map const &hitRefs) { std::vector ids; for (auto it = theHits->begin(); it != theHits->end(); ++it) ids.emplace_back((it->id().rawId())); From 1b1bcbea6a2be74a22030551be6bad815a563d5f Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 23 Sep 2022 05:27:37 +0200 Subject: [PATCH 5/8] Introduce eta plots --- .../HGCalValidation/test/HGCMissingRecHit.cc | 103 ++++++++++-------- 1 file changed, 57 insertions(+), 46 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index ab6201697a5f6..9f88b120b6f59 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -23,6 +23,7 @@ #include "DataFormats/DetId/interface/DetId.h" #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h" #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h" #include "DataFormats/HGCRecHit/interface/HGCRecHit.h" #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" @@ -78,7 +79,7 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer const &hitRefs); template void analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs); - + double getEta(const HGCHitTuple& hitRef); private: //HGC Geometry const std::vector geometrySource_, detectors_; @@ -101,8 +102,8 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer fhRecHitToken_; const edm::EDGetTokenT bhRecHitToken_; - std::vector goodHitsDE_, missedHitsDE_; - std::vector goodHitsRE_, missedHitsRE_; + std::vector goodHitsDE_, missedHitsDE_, goodHitsDT_, missedHitsDT_; + std::vector goodHitsRE_, missedHitsRE_, goodHitsRT_, missedHitsRT_; }; HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg) @@ -186,6 +187,14 @@ void HGCMissingRecHit::beginJob() { sprintf(title, "SimHit energy absent among Digis in %s", detectors_[k].c_str()); missedHitsDE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); missedHitsDE_.back()->Sumw2(); + sprintf(name, "GoodDT%s", geometrySource_[k].c_str()); + sprintf(title, "|#eta| of SimHits present among Digis in %s", detectors_[k].c_str()); + goodHitsDT_.emplace_back(fs->make(name, title, 320, 2.7, 3.1)); + goodHitsDT_.back()->Sumw2(); + sprintf(name, "MissDT%s", geometrySource_[k].c_str()); + sprintf(title, "|#eta| of SimHits absent among Digis in %s", detectors_[k].c_str()); + missedHitsDT_.emplace_back(fs->make(name, title, 320, 2.7, 3.1)); + missedHitsDT_.back()->Sumw2(); sprintf(name, "GoodRE%s", geometrySource_[k].c_str()); sprintf(title, "SimHit energy present among RecHits in %s", detectors_[k].c_str()); goodHitsRE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); @@ -194,6 +203,14 @@ void HGCMissingRecHit::beginJob() { sprintf(title, "SimHit energy absent among RecHits in %s", detectors_[k].c_str()); missedHitsRE_.emplace_back(fs->make(name, title, 1000, 0, 0.01)); missedHitsRE_.back()->Sumw2(); + sprintf(name, "GoodRT%s", geometrySource_[k].c_str()); + sprintf(title, "|#eta| of SimHits present among RecHits in %s", detectors_[k].c_str()); + goodHitsRT_.emplace_back(fs->make(name, title, 320, 2.7, 3.1)); + goodHitsRT_.back()->Sumw2(); + sprintf(name, "MissRT%s", geometrySource_[k].c_str()); + sprintf(title, "|#eta| of SimHits absent among RecHits in %s", detectors_[k].c_str()); + missedHitsRT_.emplace_back(fs->make(name, title, 320, 2.7, 3.1)); + missedHitsRT_.back()->Sumw2(); } } @@ -321,61 +338,38 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> const &simHits, int idet, std::map &hitRefs) { - const bool debug(false); for (auto const &simHit : *simHits) { - unsigned int id = simHit.id(); - std::pair xy; + DetId id(simHit.id()); bool ok(true); - int subdet(0), zside(0), layer(0), wafer(0), celltype(0), cell(0), wafer2(0), cell2(0); + GlobalPoint p; + std::ostringstream st1; + if (id.det() == DetId::HGCalHSc) { + st1 << HGCScintillatorDetId(id); + } else if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) { + st1 << HGCSiliconDetId(id); + } else { + st1 << "Not a Standard One"; + } if ((hgcCons_[idet]->waferHexagon8()) && - ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi))) { - HGCSiliconDetId detId = HGCSiliconDetId(id); - subdet = static_cast(detId.det()); - cell = detId.cellU(); - cell2 = detId.cellV(); - wafer = detId.waferU(); - wafer2 = detId.waferV(); - celltype = detId.type(); - layer = detId.layer(); - zside = detId.zside(); - xy = hgcCons_[idet]->locateCell(layer, wafer, wafer2, cell, cell2, true, true, debug); - } else if ((hgcCons_[idet]->tileTrapezoid()) && (DetId(id).det() == DetId::HGCalHSc)) { - HGCScintillatorDetId detId = HGCScintillatorDetId(id); - subdet = static_cast(detId.det()); - wafer = detId.ietaAbs(); - cell = detId.iphi(); - celltype = detId.type(); - layer = detId.layer(); - zside = detId.zside(); - xy = hgcCons_[idet]->locateCellTrap(layer, wafer, cell, false, debug); - edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " LocateCellTrap i/p " << layer << ":" - << wafer << ":" << cell << " o/p " << xy.first << ":" << xy.second; + ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi))) { + p = hgcGeometry_[idet]->getPosition(id); + } else if ((hgcCons_[idet]->tileTrapezoid()) && (id.det() == DetId::HGCalHSc)) { + p = hgcGeometry_[idet]->getPosition(id); + edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " position (" << p.x() << ", " << p.y() << ", " << p.z() << ")"; } else { // This is an invalid cell ok = false; - std::ostringstream st1; - if (DetId(id).det() == DetId::HGCalHSc) { - st1 << HGCScintillatorDetId(id); - } else if ((DetId(id).det() == DetId::HGCalEE) || (DetId(id).det() == DetId::HGCalHSi)) { - st1 << HGCSiliconDetId(id); - } else { - st1 << "Not a Standard One"; - } - edm::LogVerbatim("HGCalError") << "Hit " << std::hex << id << std::dec << " " << st1.str() + edm::LogVerbatim("HGCalError") << "Hit " << std::hex << id.rawId() << std::dec << " " << st1.str() << " in the wrong collection for detector " << idet << ":" << geometrySource_[idet] << " ***** ERROR *****"; } - edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id << std::dec << " (" << subdet << ":" << zside << ":" - << layer << ":" << celltype << ":" << wafer << ":" << wafer2 << ":" << cell << ":" - << cell2 << ") Flag " << ok; + edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id.rawId() << std::dec << " " << st1.str() << " Flag " << ok; if (ok) { - float zp = hgcCons_[idet]->waferZ(layer, false); - if (zside < 0) - zp = -zp; - float xp = (zside < 0) ? -xy.first / 10 : xy.first / 10; - float yp = xy.second / 10.0; + float xp = p.x(); + float yp = p.y(); + float zp = p.z(); float energy = simHit.energy(); float energySum(energy); @@ -396,9 +390,11 @@ void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, for (auto it = theHits->begin(); it != theHits->end(); ++it) ids.emplace_back((it->id().rawId())); for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { + double eta = getEta(it->second); auto itr = std::find(ids.begin(), ids.end(), it->first); if (itr == ids.end()) { missedHitsDE_[idet]->Fill(std::get<0>(it->second)); + missedHitsDT_[idet]->Fill(eta); std::ostringstream st1; if (DetId(it->first).det() == DetId::HGCalHSc) st1 << HGCScintillatorDetId(it->first); @@ -410,6 +406,7 @@ void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, << ") is missing in the Digi collection"; } else { goodHitsDE_[idet]->Fill(std::get<0>(it->second)); + goodHitsDT_[idet]->Fill(eta); } } } @@ -422,9 +419,11 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, for (auto it = theHits->begin(); it != theHits->end(); ++it) ids.emplace_back((it->id().rawId())); for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { + double eta = getEta(it->second); auto itr = std::find(ids.begin(), ids.end(), it->first); if (itr == ids.end()) { missedHitsRE_[idet]->Fill(std::get<0>(it->second)); + missedHitsRT_[idet]->Fill(eta); std::ostringstream st1; if (DetId(it->first).det() == DetId::HGCalHSc) st1 << HGCScintillatorDetId(it->first); @@ -436,9 +435,21 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, << ") is missing in the RecHit collection"; } else { goodHitsRE_[idet]->Fill(std::get<0>(it->second)); + goodHitsRT_[idet]->Fill(eta); } } } +double HGCMissingRecHit::getEta(const HGCHitTuple& hitRef) { + double x = std::get<1>(hitRef); + double y = std::get<2>(hitRef); + double r = std::sqrt(x * x + y * y); + double z = std::abs(std::get<3>(hitRef)); + double theta = std::atan(r / z); + double eta = (z > 0) ? -std::log(std::tan(0.5 * theta)) : 100.0; + edm::LogVerbatim("HGCalValid") << " x:y:z:r:theta:eta " << x << ":" << y << ":" << z << ":" << r << ":" << theta << ":" << eta; + return eta; +} + //define this as a plug-in DEFINE_FWK_MODULE(HGCMissingRecHit); From fb602d2df094e8fa2b442c48864dd29fe775ed4e Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 23 Sep 2022 06:21:19 +0200 Subject: [PATCH 6/8] Code check --- .../HGCalValidation/test/HGCMissingRecHit.cc | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index 9f88b120b6f59..a78235985902e 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -79,7 +79,8 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer const &hitRefs); template void analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs); - double getEta(const HGCHitTuple& hitRef); + double getEta(const HGCHitTuple &hitRef); + private: //HGC Geometry const std::vector geometrySource_, detectors_; @@ -350,12 +351,12 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con } else { st1 << "Not a Standard One"; } - if ((hgcCons_[idet]->waferHexagon8()) && - ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi))) { + if ((hgcCons_[idet]->waferHexagon8()) && ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi))) { p = hgcGeometry_[idet]->getPosition(id); } else if ((hgcCons_[idet]->tileTrapezoid()) && (id.det() == DetId::HGCalHSc)) { p = hgcGeometry_[idet]->getPosition(id); - edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " position (" << p.x() << ", " << p.y() << ", " << p.z() << ")"; + edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " position (" << p.x() << ", " << p.y() + << ", " << p.z() << ")"; } else { // This is an invalid cell ok = false; @@ -364,7 +365,8 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con << " ***** ERROR *****"; } - edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id.rawId() << std::dec << " " << st1.str() << " Flag " << ok; + edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id.rawId() << std::dec << " " << st1.str() << " Flag " + << ok; if (ok) { float xp = p.x(); @@ -440,14 +442,15 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, } } -double HGCMissingRecHit::getEta(const HGCHitTuple& hitRef) { +double HGCMissingRecHit::getEta(const HGCHitTuple &hitRef) { double x = std::get<1>(hitRef); double y = std::get<2>(hitRef); double r = std::sqrt(x * x + y * y); double z = std::abs(std::get<3>(hitRef)); double theta = std::atan(r / z); double eta = (z > 0) ? -std::log(std::tan(0.5 * theta)) : 100.0; - edm::LogVerbatim("HGCalValid") << " x:y:z:r:theta:eta " << x << ":" << y << ":" << z << ":" << r << ":" << theta << ":" << eta; + edm::LogVerbatim("HGCalValid") << " x:y:z:r:theta:eta " << x << ":" << y << ":" << z << ":" << r << ":" << theta + << ":" << eta; return eta; } From c2b923cf40c3c6b2ba0e554594e509ac518603d0 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 23 Sep 2022 11:18:59 +0200 Subject: [PATCH 7/8] Take Marco's comments --- .../HGCalValidation/test/HGCMissingRecHit.cc | 47 +++++-------------- 1 file changed, 12 insertions(+), 35 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index a78235985902e..e85ab87530fa9 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -63,7 +63,7 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer HGCHitTuple; + typedef std::tuple HGCHitTuple; void beginJob() override; void endJob() override {} @@ -79,7 +79,6 @@ class HGCMissingRecHit : public edm::one::EDAnalyzer const &hitRefs); template void analyzeHGCalRecHit(T1 const &theHits, int idet, std::map const &hitRefs); - double getEta(const HGCHitTuple &hitRef); private: //HGC Geometry @@ -244,9 +243,8 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & for (std::map::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) { int idx = std::distance(eeHitRefs.begin(), itr); edm::LogVerbatim("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " - << std::get<0>(itr->second) << "; Position" - << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", " - << std::get<3>(itr->second) << ")"; + << std::get<0>(itr->second) << "; Position = " + << std::get<1>(itr->second) << ")"; } } else { edm::LogWarning("HGCalValid") << "No EE SimHit Found " << std::endl; @@ -259,9 +257,8 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & for (std::map::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) { int idx = std::distance(fhHitRefs.begin(), itr); edm::LogVerbatim("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " - << std::get<0>(itr->second) << "; Position" - << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", " - << std::get<3>(itr->second) << ")"; + << std::get<0>(itr->second) << "; Position = " + << std::get<1>(itr->second) << ")"; } } else { edm::LogWarning("HGCalValid") << "No FH SimHit Found " << std::endl; @@ -274,8 +271,7 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & for (std::map::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) { int idx = std::distance(bhHitRefs.begin(), itr); edm::LogVerbatim("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " - << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second) << ", " - << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ")"; + << std::get<0>(itr->second) << "; Position = (" << std::get<1>(itr->second) << ")"; } } else { edm::LogWarning("HGCalValid") << "No BH SimHit Found " << std::endl; @@ -369,17 +365,13 @@ void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle> con << ok; if (ok) { - float xp = p.x(); - float yp = p.y(); - float zp = p.z(); float energy = simHit.energy(); float energySum(energy); if (hitRefs.count(id) != 0) energySum += std::get<0>(hitRefs[id]); - hitRefs[id] = std::make_tuple(energySum, xp, yp, zp); - edm::LogVerbatim("HGCalValid") << "Position (" << xp << ", " << yp << ", " << zp << ") " - << " Energy " << simHit.energy() << ":" << energySum; + hitRefs[id] = std::make_tuple(energySum, p); + edm::LogVerbatim("HGCalValid") << "Position = " << p << " Energy " << simHit.energy() << ":" << energySum; } } } @@ -392,7 +384,7 @@ void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, for (auto it = theHits->begin(); it != theHits->end(); ++it) ids.emplace_back((it->id().rawId())); for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { - double eta = getEta(it->second); + double eta = std::get<1>(it->second).eta(); auto itr = std::find(ids.begin(), ids.end(), it->first); if (itr == ids.end()) { missedHitsDE_[idet]->Fill(std::get<0>(it->second)); @@ -403,9 +395,7 @@ void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, else st1 << HGCSiliconDetId(it->first); edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() - << " SimHit (E = " << std::get<0>(it->second) << ", X = " << std::get<1>(it->second) - << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) - << ") is missing in the Digi collection"; + << " SimHit (E = " << std::get<0>(it->second) << ", Position = " << std::get<1>(it->second) << ") is missing in the Digi collection"; } else { goodHitsDE_[idet]->Fill(std::get<0>(it->second)); goodHitsDT_[idet]->Fill(eta); @@ -421,7 +411,7 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, for (auto it = theHits->begin(); it != theHits->end(); ++it) ids.emplace_back((it->id().rawId())); for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) { - double eta = getEta(it->second); + double eta = std::get<1>(it->second).eta(); auto itr = std::find(ids.begin(), ids.end(), it->first); if (itr == ids.end()) { missedHitsRE_[idet]->Fill(std::get<0>(it->second)); @@ -432,8 +422,7 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, else st1 << HGCSiliconDetId(it->first); edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() - << " SimHit (E = " << std::get<0>(it->second) << ", X = " << std::get<1>(it->second) - << ", Y = " << std::get<2>(it->second) << ", Z = " << std::get<3>(it->second) + << " SimHit (E = " << std::get<0>(it->second) << ", Position = " << std::get<1>(it->second) << ") is missing in the RecHit collection"; } else { goodHitsRE_[idet]->Fill(std::get<0>(it->second)); @@ -442,17 +431,5 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, } } -double HGCMissingRecHit::getEta(const HGCHitTuple &hitRef) { - double x = std::get<1>(hitRef); - double y = std::get<2>(hitRef); - double r = std::sqrt(x * x + y * y); - double z = std::abs(std::get<3>(hitRef)); - double theta = std::atan(r / z); - double eta = (z > 0) ? -std::log(std::tan(0.5 * theta)) : 100.0; - edm::LogVerbatim("HGCalValid") << " x:y:z:r:theta:eta " << x << ":" << y << ":" << z << ":" << r << ":" << theta - << ":" << eta; - return eta; -} - //define this as a plug-in DEFINE_FWK_MODULE(HGCMissingRecHit); From b66d829b23601d983cdd133933663d94bcd9efd8 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Fri, 23 Sep 2022 11:46:13 +0200 Subject: [PATCH 8/8] Code check --- Validation/HGCalValidation/test/HGCMissingRecHit.cc | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Validation/HGCalValidation/test/HGCMissingRecHit.cc b/Validation/HGCalValidation/test/HGCMissingRecHit.cc index e85ab87530fa9..879890ea5b2df 100644 --- a/Validation/HGCalValidation/test/HGCMissingRecHit.cc +++ b/Validation/HGCalValidation/test/HGCMissingRecHit.cc @@ -243,8 +243,7 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & for (std::map::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) { int idx = std::distance(eeHitRefs.begin(), itr); edm::LogVerbatim("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " - << std::get<0>(itr->second) << "; Position = " - << std::get<1>(itr->second) << ")"; + << std::get<0>(itr->second) << "; Position = " << std::get<1>(itr->second) << ")"; } } else { edm::LogWarning("HGCalValid") << "No EE SimHit Found " << std::endl; @@ -257,8 +256,7 @@ void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup & for (std::map::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) { int idx = std::distance(fhHitRefs.begin(), itr); edm::LogVerbatim("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy " - << std::get<0>(itr->second) << "; Position = " - << std::get<1>(itr->second) << ")"; + << std::get<0>(itr->second) << "; Position = " << std::get<1>(itr->second) << ")"; } } else { edm::LogWarning("HGCalValid") << "No FH SimHit Found " << std::endl; @@ -395,7 +393,9 @@ void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits, else st1 << HGCSiliconDetId(it->first); edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() - << " SimHit (E = " << std::get<0>(it->second) << ", Position = " << std::get<1>(it->second) << ") is missing in the Digi collection"; + << " SimHit (E = " << std::get<0>(it->second) + << ", Position = " << std::get<1>(it->second) + << ") is missing in the Digi collection"; } else { goodHitsDE_[idet]->Fill(std::get<0>(it->second)); goodHitsDT_[idet]->Fill(eta); @@ -422,7 +422,8 @@ void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits, else st1 << HGCSiliconDetId(it->first); edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str() - << " SimHit (E = " << std::get<0>(it->second) << ", Position = " << std::get<1>(it->second) + << " SimHit (E = " << std::get<0>(it->second) + << ", Position = " << std::get<1>(it->second) << ") is missing in the RecHit collection"; } else { goodHitsRE_[idet]->Fill(std::get<0>(it->second));