File indexing completed on 2023-03-17 13:03:38
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 uint32_t lines(0);
0086 const std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0087 while (fInput.getline(buffer, 80)) {
0088 ++lines;
0089 std::vector<std::string> items = HGCalGeomUtils::splitString(std::string(buffer));
0090 if (items.size() == 8) {
0091 DetId::Detector det = static_cast<DetId::Detector>(std::atoi(items[0].c_str()));
0092 auto itr = std::find(dets.begin(), dets.end(), det);
0093 if (itr != dets.end()) {
0094 uint32_t pos = static_cast<uint32_t>(itr - dets.begin());
0095 DetId id(0);
0096 if ((det == DetId::HGCalEE) || (det == DetId::HGCalHSi)) {
0097 int type = std::atoi(items[1].c_str());
0098 int zside = std::atoi(items[2].c_str());
0099 int layer = std::atoi(items[3].c_str());
0100 int waferU = std::atoi(items[4].c_str());
0101 int waferV = std::atoi(items[5].c_str());
0102 int cellU = std::atoi(items[6].c_str());
0103 int cellV = std::atoi(items[7].c_str());
0104 id = static_cast<DetId>(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0105 } else if (det == DetId::HGCalHSc) {
0106 bool trig = (std::atoi(items[1].c_str()) > 0);
0107 int type = std::atoi(items[2].c_str());
0108 int zside = std::atoi(items[3].c_str());
0109 int sipm = std::atoi(items[4].c_str());
0110 int layer = std::atoi(items[5].c_str());
0111 int ring = std::atoi(items[6].c_str());
0112 int iphi = std::atoi(items[7].c_str());
0113 id = static_cast<DetId>(HGCScintillatorDetId(type, layer, zside * ring, iphi, trig, sipm));
0114 }
0115 if (id.rawId() != 0)
0116 detIds_.emplace_back(id, pos);
0117 }
0118 }
0119 }
0120 fInput.close();
0121 }
0122 }
0123 edm::LogVerbatim("HGCGeom") << "Reads " << detIds_.size() << " ID's from " << fileName_;
0124 for (unsigned int k = 0; k < detIds_.size(); ++k) {
0125 if ((detIds_[k].first).det() == DetId::HGCalHSc)
0126 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << HGCScintillatorDetId(detIds_[k].first) << " from DDConstant "
0127 << (detIds_[k].second);
0128 else
0129 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << HGCSiliconDetId(detIds_[k].first) << " from DDConstant "
0130 << (detIds_[k].second);
0131 }
0132 }
0133
0134 void HGCalValidityTester::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0135 std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0136 edm::ParameterSetDescription desc;
0137 desc.add<std::vector<std::string>>("nameDetectors", names);
0138 desc.add<std::string>("fileName", "missD88.txt");
0139 descriptions.add("hgcalValidityTester", desc);
0140 }
0141
0142
0143 void HGCalValidityTester::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
0144
0145 std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0146 std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0147 std::map<DetId::Detector, uint32_t> detMap;
0148 for (uint32_t i = 0; i < nameDetectors_.size(); i++) {
0149 edm::LogVerbatim("HGCGeom") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i << ":"
0150 << nameDetectors_[i];
0151 const edm::ESHandle<HGCalDDDConstants> &hgcCons = iSetup.getHandle(tok_hgcal_[i]);
0152 if (hgcCons.isValid()) {
0153 hgcCons_.push_back(hgcCons.product());
0154 } else {
0155 edm::LogWarning("HGCGeom") << "Cannot initiate HGCalDDDConstants for " << nameDetectors_[i] << std::endl;
0156 }
0157 auto ii = std::find(names.begin(), names.end(), nameDetectors_[i]);
0158 if (ii != names.end()) {
0159 uint32_t k = static_cast<uint32_t>(ii - names.begin());
0160 detMap[dets[k]] = i;
0161 }
0162 }
0163 edm::LogVerbatim("HGCGeom") << "Loaded HGCalDDConstants for " << detMap.size() << " detectors";
0164
0165 for (auto itr = detMap.begin(); itr != detMap.end(); ++itr)
0166 edm::LogVerbatim("HGCGeom") << "[" << itr->second << "]: " << nameDetectors_[itr->second] << " for Detector "
0167 << itr->first;
0168
0169 for (unsigned int k = 0; k < detIds_.size(); ++k) {
0170 std::ostringstream st1;
0171 std::pair<float, float> xy;
0172 bool valid(false);
0173 int layer(0), zside(1);
0174 const HGCalDDDConstants *cons = hgcCons_[detMap[(detIds_[k].first).det()]];
0175 if ((detIds_[k].first).det() == DetId::HGCalHSc) {
0176 HGCScintillatorDetId id(detIds_[k].first);
0177 layer = id.layer();
0178 zside = id.zside();
0179 xy = cons->locateCellTrap(zside, layer, id.ring(), id.iphi(), true, false);
0180 double z = zside * (cons->waferZ(layer, true));
0181 valid = cons->isValidTrap(zside, layer, id.ring(), id.iphi());
0182 auto cell = cons->assignCellTrap(zside * xy.first, xy.second, z, layer, true);
0183 HGCScintillatorDetId newId(cell[2], layer, zside * cell[0], cell[1], id.trigger(), id.sipm());
0184 st1 << "Old: " << id << " New: " << newId << " OK " << (id.rawId() == newId.rawId());
0185 } else {
0186 HGCSiliconDetId id(detIds_[k].first);
0187 layer = id.layer();
0188 zside = id.zside();
0189 xy = cons->locateCell(zside, layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), true, true, false, false);
0190 valid = cons->isValidHex8(layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), false);
0191 auto cell = cons->assignCellHex(zside * xy.first, xy.second, zside, layer, true, false, false);
0192 HGCSiliconDetId newId(id.det(), id.zside(), cell[2], id.layer(), cell[0], cell[1], cell[3], cell[4]);
0193 st1 << "Old: " << id << " New: " << newId << " OK " << (id.rawId() == newId.rawId());
0194 }
0195 double r = std::sqrt(xy.first * xy.first + xy.second * xy.second);
0196 double z = zside * (cons->waferZ(layer, true));
0197 auto range = cons->getRangeR(layer, true);
0198 edm::LogVerbatim("HGCalMiss") << "Hit[" << k << "] " << st1.str() << " Position (" << xy.first << ", " << xy.second
0199 << ", " << z << ") Valid " << valid << " R " << r << " (" << range.first << ":"
0200 << range.second << ")";
0201 }
0202 }
0203
0204
0205 DEFINE_FWK_MODULE(HGCalValidityTester);
