File indexing completed on 2024-04-06 12:32:40
0001
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <map>
0006 #include <mutex>
0007 #include <string>
0008 #include <vector>
0009
0010 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0011 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0012 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0013 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0014 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0015
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/stream/EDAnalyzer.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 #include "FWCore/Utilities/interface/ESGetToken.h"
0025
0026 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0027 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0028 #include "Geometry/HGCalCommonData/interface/HGCalWaferIndex.h"
0029
0030 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0031 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0032
0033 namespace HGCalValidSimhitCheck {
0034 struct Counters {
0035 Counters() {
0036 badTypes_.clear();
0037 occupancy_.clear();
0038 goodChannels_.clear();
0039 }
0040 CMS_THREAD_GUARD(mtx_) mutable std::map<int, int> badTypes_, occupancy_;
0041 CMS_THREAD_GUARD(mtx_) mutable std::vector<int> goodChannels_;
0042 mutable std::mutex mtx_;
0043 };
0044 }
0045
0046 class HGCalWaferHitCheck : public edm::stream::EDAnalyzer<edm::GlobalCache<HGCalValidSimhitCheck::Counters> > {
0047 public:
0048 explicit HGCalWaferHitCheck(const edm::ParameterSet&, const HGCalValidSimhitCheck::Counters* count);
0049
0050 static std::unique_ptr<HGCalValidSimhitCheck::Counters> initializeGlobalCache(edm::ParameterSet const& iConfig) {
0051 return std::make_unique<HGCalValidSimhitCheck::Counters>();
0052 }
0053
0054 void analyze(edm::Event const&, edm::EventSetup const&) override;
0055 void endStream() override;
0056 static void globalEndJob(const HGCalValidSimhitCheck::Counters* counters);
0057 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0058
0059 protected:
0060 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0061 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0062
0063 private:
0064 template <class T>
0065 void analyzeHits(const std::string&, const T&);
0066
0067
0068 enum inputType { Sim = 1, Digi = 2, Reco = 3 };
0069 const std::string nameDetector_, caloHitSource_;
0070 const edm::InputTag source_;
0071 const int inpType_;
0072 const int verbosity_;
0073 const bool ifNose_;
0074 const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> geomToken_;
0075 const edm::EDGetTokenT<HGCalDigiCollection> digiSource_;
0076 const edm::EDGetTokenT<HGCRecHitCollection> recHitSource_;
0077 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hit_;
0078 const HGCalDDDConstants* hgcons_;
0079 std::map<int, int> badTypes_, occupancy_;
0080 std::vector<int> goodChannels_;
0081 static const int waferMax = 12;
0082 static const int layerMax = 28;
0083 };
0084
0085 HGCalWaferHitCheck::HGCalWaferHitCheck(const edm::ParameterSet& iConfig, const HGCalValidSimhitCheck::Counters*)
0086 : nameDetector_(iConfig.getParameter<std::string>("detectorName")),
0087 caloHitSource_(iConfig.getParameter<std::string>("caloHitSource")),
0088 source_(iConfig.getParameter<edm::InputTag>("source")),
0089 inpType_(iConfig.getParameter<int>("inputType")),
0090 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0091 ifNose_(iConfig.getUntrackedParameter<bool>("ifNose", false)),
0092 geomToken_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0093 edm::ESInputTag{"", nameDetector_})),
0094 digiSource_(consumes<HGCalDigiCollection>(source_)),
0095 recHitSource_(consumes<HGCRecHitCollection>(source_)),
0096 tok_hit_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", caloHitSource_))) {
0097 edm::LogVerbatim("HGCalValidation") << "Validation for input type " << inpType_ << " for " << nameDetector_;
0098 }
0099
0100 void HGCalWaferHitCheck::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0101 edm::ParameterSetDescription desc;
0102 desc.add<std::string>("detectorName", "HGCalEESensitive");
0103 desc.add<std::string>("caloHitSource", "HGCHitsEE");
0104 desc.add<edm::InputTag>("source", edm::InputTag("simHGCalUnsuppressedDigis", "EE"));
0105 desc.add<int>("inputType", 1);
0106 desc.addUntracked<int>("verbosity", 0);
0107 desc.addUntracked<bool>("ifNose", false);
0108 descriptions.add("hgcalWaferHitCheckEE", desc);
0109 }
0110
0111 void HGCalWaferHitCheck::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0112
0113 if (inpType_ <= Sim) {
0114 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainer = iEvent.getHandle(tok_hit_);
0115 if (theCaloHitContainer.isValid()) {
0116 if (verbosity_ > 0)
0117 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << theCaloHitContainer->size() << " SimHits";
0118 analyzeHits(nameDetector_, *(theCaloHitContainer.product()));
0119 } else if (verbosity_ > 0) {
0120 edm::LogVerbatim("HGCalValidation") << "PCaloHitContainer does not exist for " << nameDetector_;
0121 }
0122 } else if (inpType_ == Digi) {
0123 const edm::Handle<HGCalDigiCollection>& theDigiContainer = iEvent.getHandle(digiSource_);
0124 if (theDigiContainer.isValid()) {
0125 if (verbosity_ > 0)
0126 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << theDigiContainer->size() << " Digis";
0127 analyzeHits(nameDetector_, *(theDigiContainer.product()));
0128 } else if (verbosity_ > 0) {
0129 edm::LogVerbatim("HGCalValidation") << "DigiContainer does not exist for " << nameDetector_;
0130 }
0131 } else {
0132 const edm::Handle<HGCRecHitCollection>& theRecHitContainer = iEvent.getHandle(recHitSource_);
0133 if (theRecHitContainer.isValid()) {
0134 if (verbosity_ > 0)
0135 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << theRecHitContainer->size() << " hits";
0136 analyzeHits(nameDetector_, *(theRecHitContainer.product()));
0137 } else if (verbosity_ > 0) {
0138 edm::LogVerbatim("HGCalValidation") << "RecHitContainer does not exist for " << nameDetector_;
0139 }
0140 }
0141 }
0142
0143 template <class T>
0144 void HGCalWaferHitCheck::analyzeHits(std::string const& name, T const& hits) {
0145 for (auto const& hit : hits) {
0146 uint32_t id = hit.id();
0147 int zside, type, layer, waferU, waferV;
0148 if (ifNose_) {
0149 HFNoseDetId detId = HFNoseDetId(id);
0150 waferU = detId.waferU();
0151 waferV = detId.waferV();
0152 type = detId.type();
0153 layer = detId.layer();
0154 zside = detId.zside();
0155 } else {
0156 HGCSiliconDetId detId = HGCSiliconDetId(id);
0157 waferU = detId.waferU();
0158 waferV = detId.waferV();
0159 type = detId.type();
0160 layer = detId.layer();
0161 zside = detId.zside();
0162 }
0163 int index = HGCalWaferIndex::waferIndex(layer, waferU, waferV, false);
0164 int typef = hgcons_->waferType(layer, waferU, waferV, true);
0165 if (zside < 0)
0166 index *= -1;
0167 auto itr = occupancy_.find(index);
0168 if (itr == occupancy_.end())
0169 occupancy_[index] = 1;
0170 else
0171 ++occupancy_[index];
0172 if (type != typef) {
0173 auto ktr = badTypes_.find(index);
0174 if (ktr == badTypes_.end())
0175 badTypes_[index] = 1;
0176 else
0177 ++badTypes_[index];
0178 if (verbosity_ == 1)
0179 edm::LogVerbatim("HGCalValidation")
0180 << "Detector " << name << " zside = " << zside << " layer = " << layer << " type = " << type << ":" << typef
0181 << " wafer = " << waferU << ":" << waferV << " index " << index;
0182 }
0183 if (verbosity_ > 1)
0184 edm::LogVerbatim("HGCalValidation")
0185 << "Detector " << name << " zside = " << zside << " layer = " << layer << " type = " << type << ":" << typef
0186 << " wafer = " << waferU << ":" << waferV;
0187 }
0188 }
0189
0190
0191 void HGCalWaferHitCheck::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0192 hgcons_ = &iSetup.getData(geomToken_);
0193 if (verbosity_ > 0)
0194 edm::LogVerbatim("HGCalValidation") << nameDetector_ << " defined with " << hgcons_->layers(false)
0195 << " Layers with " << (hgcons_->firstLayer() - 1) << " in front";
0196
0197 goodChannels_.clear();
0198 for (int iz = 0; iz < 2; ++iz) {
0199 int zside = iz * 2 - 1;
0200 for (int layer = 1; layer <= layerMax; ++layer) {
0201 for (int waferU = -waferMax; waferU <= waferMax; ++waferU) {
0202 for (int waferV = -waferMax; waferV <= waferMax; ++waferV) {
0203 int index = zside * HGCalWaferIndex::waferIndex(layer, waferU, waferV, false);
0204 if (hgcons_->isValidHex8(layer, waferU, waferV, true))
0205 goodChannels_.emplace_back(index);
0206 }
0207 }
0208 }
0209 }
0210 }
0211
0212 void HGCalWaferHitCheck::endStream() {
0213 std::scoped_lock lock(globalCache()->mtx_);
0214 for (auto [id, count] : occupancy_) {
0215 if (globalCache()->occupancy_.find(id) == globalCache()->occupancy_.end())
0216 globalCache()->occupancy_[id] = count;
0217 else
0218 globalCache()->occupancy_[id] += count;
0219 }
0220 for (auto [id, count] : badTypes_) {
0221 if (globalCache()->badTypes_.find(id) == globalCache()->badTypes_.end())
0222 globalCache()->badTypes_[id] = count;
0223 else
0224 globalCache()->badTypes_[id] += count;
0225 }
0226
0227 globalCache()->goodChannels_ = goodChannels_;
0228 globalCache()->mtx_.unlock();
0229 }
0230
0231 void HGCalWaferHitCheck::globalEndJob(const HGCalValidSimhitCheck::Counters* count) {
0232 int allbad(0), nocc(0);
0233 for (auto const& index : count->goodChannels_) {
0234 int zside = (index < 0) ? -1 : 1;
0235 int layer = HGCalWaferIndex::waferLayer(std::abs(index));
0236 int waferU = HGCalWaferIndex::waferU(std::abs(index));
0237 int waferV = HGCalWaferIndex::waferV(std::abs(index));
0238 int occ = (count->occupancy_.find(index) == count->occupancy_.end()) ? 0 : count->occupancy_[index];
0239 int bad = (count->badTypes_.find(index) == count->badTypes_.end()) ? 0 : count->badTypes_[index];
0240 if (occ == 0)
0241 ++nocc;
0242 if (bad > 0)
0243 ++allbad;
0244 if (occ == 0 || bad > 0) {
0245 edm::LogVerbatim("HGCalValidation") << "ZS:Layer:u:v:index " << zside << ":" << layer << ":" << waferU << ":"
0246 << waferV << ":" << index << " Occ " << occ << " bad " << bad;
0247 }
0248 }
0249 edm::LogVerbatim("HGCalValidation") << "\n\n"
0250 << allbad << " channels with bad types among " << count->goodChannels_.size()
0251 << " channels and " << nocc << " channels with zero occupancy\n\n";
0252 }
0253
0254 #include "FWCore/Framework/interface/MakerMacros.h"
0255
0256 DEFINE_FWK_MODULE(HGCalWaferHitCheck);
