File indexing completed on 2024-04-06 12:15:09
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 "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0038 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0039 #include "Geometry/HGCalCommonData/interface/HGCalGeomUtils.h"
0040 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0041
0042 class HGCalScintIDTester : public edm::one::EDAnalyzer<> {
0043 public:
0044 explicit HGCalScintIDTester(const edm::ParameterSet&);
0045 ~HGCalScintIDTester() override = default;
0046 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048 void beginJob() override {}
0049 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0050 void endJob() override {}
0051
0052 private:
0053 const std::string nameSense_, errorFile_;
0054 const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> dddToken_;
0055 std::vector<HGCScintillatorDetId> detIds_;
0056 std::vector<std::pair<double, double>> posXY_;
0057 };
0058
0059 HGCalScintIDTester::HGCalScintIDTester(const edm::ParameterSet& iC)
0060 : nameSense_(iC.getParameter<std::string>("nameSense")),
0061 errorFile_(iC.getParameter<std::string>("fileName")),
0062 dddToken_(esConsumes<HGCalDDDConstants, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})) {
0063 edm::LogVerbatim("HGCalGeom") << "Test HGCScintillator DetID for " << nameSense_ << " of positions from the file "
0064 << errorFile_;
0065
0066 edm::FileInPath filetmp("Geometry/HGCalCommonData/data/" + errorFile_);
0067 std::string fileName = filetmp.fullPath();
0068 std::ifstream fInput(fileName.c_str());
0069 if (!fInput.good()) {
0070 edm::LogVerbatim("HGCalSim") << "Cannot open file " << fileName;
0071 } else {
0072 char buffer[80];
0073 int kount(0);
0074 while (fInput.getline(buffer, 80)) {
0075 std::vector<std::string> items = HGCalGeomUtils::splitString(std::string(buffer));
0076 ++kount;
0077 if (items.size() > 8) {
0078 bool trig = (std::atoi(items[1].c_str()) > 0);
0079 int zp = std::atoi(items[2].c_str());
0080 int type = std::atoi(items[3].c_str());
0081 int sipm = std::atoi(items[4].c_str());
0082 int layer = std::atoi(items[5].c_str());
0083 int ring = std::atoi(items[6].c_str());
0084 int iphi = std::atoi(items[7].c_str());
0085 double xx = std::atof(items[8].c_str());
0086 double yy = std::atof(items[9].c_str());
0087 HGCScintillatorDetId id(type, layer, zp * ring, iphi, trig, sipm);
0088 detIds_.emplace_back(id);
0089 posXY_.emplace_back(std::make_pair(xx, yy));
0090 }
0091 }
0092 fInput.close();
0093 edm::LogVerbatim("HGCalGeom") << "Reads a total of " << detIds_.size() << ":" << posXY_.size() << " entries out of "
0094 << kount << "\n";
0095 for (unsigned int k = 0; k < detIds_.size(); ++k)
0096 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << detIds_[k] << " (" << posXY_[k].first << ", "
0097 << posXY_[k].second << ")";
0098 }
0099 }
0100
0101 void HGCalScintIDTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0102 edm::ParameterSetDescription desc;
0103 desc.add<std::string>("nameSense", "HGCalHEScintillatorSensitive");
0104 desc.add<std::string>("fileName", "errorScintD88.txt");
0105 descriptions.add("hgcalScintIDTester", desc);
0106 }
0107
0108
0109 void HGCalScintIDTester::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110 const HGCalDDDConstants& hgdc = iSetup.getData(dddToken_);
0111 edm::LogVerbatim("HGCalGeom") << "\nStart testing " << nameSense_ << std::endl;
0112
0113 for (unsigned int k = 0; k < detIds_.size(); ++k) {
0114 std::ostringstream st1;
0115 st1 << "Hit[" << k << "] " << detIds_[k];
0116 double xx = posXY_[k].first;
0117 double yy = posXY_[k].second;
0118 int layer = detIds_[k].layer();
0119 int zside = detIds_[k].zside();
0120 double zz = zside * hgdc.waferZ(layer, false);
0121 std::array<int, 3> idx = hgdc.assignCellTrap(xx, yy, zz, layer, false);
0122 HGCScintillatorDetId id(idx[2], layer, (zside * idx[0]), idx[1], false, 0);
0123 std::pair<int, int> typm = hgdc.tileType(layer, idx[0], 0);
0124 if (typm.first >= 0) {
0125 id.setType(typm.first);
0126 id.setSiPM(typm.second);
0127 }
0128 if (id.rawId() != detIds_[k].rawId())
0129 st1 << " non-matching DetId: new ID " << id;
0130 auto xy = hgdc.locateCell(id, false);
0131 double xx0 = (id.zside() > 0) ? xy.first : -xy.first;
0132 double yy0 = xy.second;
0133 double dx = xx0 - xx;
0134 double dy = yy0 - yy;
0135 double diff = std::sqrt(dx * dx + dy * dy);
0136 st1 << " input position: (" << xx << ", " << yy << "); position from ID (" << xx0 << ", " << yy0 << ") distance "
0137 << diff;
0138 constexpr double tol = 3.0;
0139 if (diff > tol)
0140 st1 << " ***** CheckID *****";
0141 bool valid1 = hgdc.isValidTrap(detIds_[k].zside(), detIds_[k].layer(), detIds_[k].ring(), detIds_[k].iphi());
0142 bool valid2 = hgdc.isValidTrap(id.zside(), id.layer(), id.ring(), id.iphi());
0143 st1 << " Validity flag: " << valid1 << ":" << valid2;
0144 if ((!valid1) || (!valid2))
0145 st1 << " +++++ Validity Check +++++ ";
0146 edm::LogVerbatim("HGCalGeom") << st1.str();
0147 }
0148 }
0149
0150
0151 DEFINE_FWK_MODULE(HGCalScintIDTester);
