File indexing completed on 2024-07-16 22:52:24
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include <fstream>
0022 #include <iostream>
0023 #include <sstream>
0024 #include <string>
0025 #include <vector>
0026
0027
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "FWCore/ParameterSet/interface/FileInPath.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0037 #include "FWCore/Utilities/interface/transform.h"
0038
0039 #include "DataFormats/DetId/interface/DetId.h"
0040 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0041 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0042 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0043 #include "Geometry/HGCalCommonData/interface/HGCalGeomUtils.h"
0044 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0045
0046 class HGCalValidityTester : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0047 public:
0048 explicit HGCalValidityTester(const edm::ParameterSet &);
0049 ~HGCalValidityTester() override = default;
0050 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0051
0052 void beginJob() override {}
0053 void beginRun(edm::Run const &, edm::EventSetup const &) override;
0054 void analyze(edm::Event const &iEvent, edm::EventSetup const &) override {}
0055 void endRun(edm::Run const &, edm::EventSetup const &) override {}
0056 void endJob() override {}
0057
0058 private:
0059 const std::vector<std::string> nameDetectors_;
0060 const std::string fileName_;
0061 const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_hgcal_;
0062 std::vector<const HGCalDDDConstants *> hgcCons_;
0063 std::vector<std::pair<DetId, uint32_t>> detIds_;
0064 };
0065
0066 HGCalValidityTester::HGCalValidityTester(const edm::ParameterSet &iC)
0067 : nameDetectors_(iC.getParameter<std::vector<std::string>>("nameDetectors")),
0068 fileName_(iC.getParameter<std::string>("fileName")),
0069 tok_hgcal_{edm::vector_transform(nameDetectors_, [this](const std::string &name) {
0070 return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0071 })} {
0072 std::ostringstream st1;
0073 for (const auto &name : nameDetectors_)
0074 st1 << " : " << name;
0075 edm::LogVerbatim("HGCGeom") << "Test validity of cells for " << nameDetectors_.size() << " detectors" << st1.str()
0076 << " with inputs from " << fileName_;
0077 if (!fileName_.empty()) {
0078 edm::FileInPath filetmp("Geometry/HGCalCommonData/data/" + fileName_);
0079 std::string fileName = filetmp.fullPath();
0080 std::ifstream fInput(fileName.c_str());
0081 if (!fInput.good()) {
0082 edm::LogVerbatim("HGCGeom") << "Cannot open file " << fileName;
0083 } else {
0084 char buffer[80];
0085 const std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0086 while (fInput.getline(buffer, 80)) {
0087 std::vector<std::string> items = HGCalGeomUtils::splitString(std::string(buffer));
0088 if (items.size() == 8) {
0089 DetId::Detector det = static_cast<DetId::Detector>(std::atoi(items[0].c_str()));
0090 auto itr = std::find(dets.begin(), dets.end(), det);
0091 if (itr != dets.end()) {
0092 uint32_t pos = static_cast<uint32_t>(itr - dets.begin());
0093 DetId id(0);
0094 if ((det == DetId::HGCalEE) || (det == DetId::HGCalHSi)) {
0095 int type = std::atoi(items[1].c_str());
0096 int zside = std::atoi(items[2].c_str());
0097 int layer = std::atoi(items[3].c_str());
0098 int waferU = std::atoi(items[4].c_str());
0099 int waferV = std::atoi(items[5].c_str());
0100 int cellU = std::atoi(items[6].c_str());
0101 int cellV = std::atoi(items[7].c_str());
0102 id = static_cast<DetId>(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0103 } else if (det == DetId::HGCalHSc) {
0104 bool trig = (std::atoi(items[1].c_str()) > 0);
0105 int type = std::atoi(items[2].c_str());
0106 int zside = std::atoi(items[3].c_str());
0107 int sipm = std::atoi(items[4].c_str());
0108 int layer = std::atoi(items[5].c_str());
0109 int ring = std::atoi(items[6].c_str());
0110 int iphi = std::atoi(items[7].c_str());
0111 id = static_cast<DetId>(HGCScintillatorDetId(type, layer, zside * ring, iphi, trig, sipm));
0112 }
0113 if (id.rawId() != 0)
0114 detIds_.emplace_back(id, pos);
0115 }
0116 }
0117 }
0118 fInput.close();
0119 }
0120 }
0121 edm::LogVerbatim("HGCGeom") << "Reads " << detIds_.size() << " ID's from " << fileName_;
0122 for (unsigned int k = 0; k < detIds_.size(); ++k) {
0123 if ((detIds_[k].first).det() == DetId::HGCalHSc)
0124 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << HGCScintillatorDetId(detIds_[k].first) << " from DDConstant "
0125 << (detIds_[k].second);
0126 else
0127 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << HGCSiliconDetId(detIds_[k].first) << " from DDConstant "
0128 << (detIds_[k].second);
0129 }
0130 }
0131
0132 void HGCalValidityTester::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0133 std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0134 edm::ParameterSetDescription desc;
0135 desc.add<std::vector<std::string>>("nameDetectors", names);
0136 desc.add<std::string>("fileName", "missD88.txt");
0137 descriptions.add("hgcalValidityTester", desc);
0138 }
0139
0140
0141 void HGCalValidityTester::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
0142
0143 std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0144 std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0145 std::map<DetId::Detector, uint32_t> detMap;
0146 for (uint32_t i = 0; i < nameDetectors_.size(); i++) {
0147 edm::LogVerbatim("HGCGeom") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i << ":"
0148 << nameDetectors_[i];
0149 const edm::ESHandle<HGCalDDDConstants> &hgcCons = iSetup.getHandle(tok_hgcal_[i]);
0150 if (hgcCons.isValid()) {
0151 hgcCons_.push_back(hgcCons.product());
0152 } else {
0153 edm::LogWarning("HGCGeom") << "Cannot initiate HGCalDDDConstants for " << nameDetectors_[i] << std::endl;
0154 }
0155 auto ii = std::find(names.begin(), names.end(), nameDetectors_[i]);
0156 if (ii != names.end()) {
0157 uint32_t k = static_cast<uint32_t>(ii - names.begin());
0158 detMap[dets[k]] = i;
0159 }
0160 }
0161 edm::LogVerbatim("HGCGeom") << "Loaded HGCalDDConstants for " << detMap.size() << " detectors";
0162
0163 for (auto itr = detMap.begin(); itr != detMap.end(); ++itr)
0164 edm::LogVerbatim("HGCGeom") << "[" << itr->second << "]: " << nameDetectors_[itr->second] << " for Detector "
0165 << itr->first;
0166
0167 for (unsigned int k = 0; k < detIds_.size(); ++k) {
0168 std::ostringstream st1;
0169 std::pair<float, float> xy;
0170 bool valid(false);
0171 int layer(0), zside(1);
0172 const HGCalDDDConstants *cons = hgcCons_[detMap[(detIds_[k].first).det()]];
0173 if ((detIds_[k].first).det() == DetId::HGCalHSc) {
0174 HGCScintillatorDetId id(detIds_[k].first);
0175 layer = id.layer();
0176 zside = id.zside();
0177 xy = cons->locateCellTrap(zside, layer, id.ring(), id.iphi(), true, false);
0178 double z = zside * (cons->waferZ(layer, true));
0179 valid = cons->isValidTrap(zside, layer, id.ring(), id.iphi());
0180 auto cell = cons->assignCellTrap(zside * xy.first, xy.second, z, layer, true);
0181 HGCScintillatorDetId newId(cell[2], layer, zside * cell[0], cell[1], id.trigger(), id.sipm());
0182 st1 << "Old: " << id << " New: " << newId << " OK " << (id.rawId() == newId.rawId());
0183 } else {
0184 HGCSiliconDetId id(detIds_[k].first);
0185 layer = id.layer();
0186 zside = id.zside();
0187 xy = cons->locateCell(
0188 zside, layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), true, true, false, false, false);
0189 valid = cons->isValidHex8(layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), false);
0190 auto cell = cons->assignCellHex(zside * xy.first, xy.second, zside, layer, true, false, false);
0191 HGCSiliconDetId newId(id.det(), id.zside(), cell[2], id.layer(), cell[0], cell[1], cell[3], cell[4]);
0192 st1 << "Old: " << id << " New: " << newId << " OK " << (id.rawId() == newId.rawId());
0193 }
0194 double r = std::sqrt(xy.first * xy.first + xy.second * xy.second);
0195 double z = zside * (cons->waferZ(layer, true));
0196 auto range = cons->getRangeR(layer, true);
0197 edm::LogVerbatim("HGCalMiss") << "Hit[" << k << "] " << st1.str() << " Position (" << xy.first << ", " << xy.second
0198 << ", " << z << ") Valid " << valid << " R " << r << " (" << range.first << ":"
0199 << range.second << ")";
0200 }
0201 }
0202
0203
0204 DEFINE_FWK_MODULE(HGCalValidityTester);
