File indexing completed on 2024-04-06 12:15:11
0001 #include <iostream>
0002 #include <sstream>
0003 #include <string>
0004 #include <vector>
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0016 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0017 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0018 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0019 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0020 #include "CoralBase/Exception.h"
0021
0022 class HGCalGeometryTester : public edm::one::EDAnalyzer<> {
0023 public:
0024 explicit HGCalGeometryTester(const edm::ParameterSet&);
0025 ~HGCalGeometryTester() override = default;
0026
0027 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0028
0029 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0030
0031 private:
0032 void doTest(const HGCalGeometry* geom, ForwardSubdetector subdet);
0033 void doTestWafer(const HGCalGeometry* geom, DetId::Detector det);
0034 void doTestScint(const HGCalGeometry* geom, DetId::Detector det);
0035
0036 const std::string name;
0037 const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0038 };
0039
0040 HGCalGeometryTester::HGCalGeometryTester(const edm::ParameterSet& iC)
0041 : name{iC.getParameter<std::string>("Detector")},
0042 geomToken_{esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", name})} {}
0043
0044 void HGCalGeometryTester::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0045 const HGCalGeometry* geom = &(iSetup.getData(geomToken_));
0046 if (geom->topology().waferHexagon6()) {
0047 ForwardSubdetector subdet;
0048 if (name == "HGCalHESiliconSensitive")
0049 subdet = HGCHEF;
0050 else if (name == "HGCalHEScintillatorSensitive")
0051 subdet = HGCHEB;
0052 else
0053 subdet = HGCEE;
0054 edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name << " Detector:Subdetector " << DetId::Forward << ":"
0055 << subdet << " Mode " << geom->topology().dddConstants().geomMode();
0056 doTest(geom, subdet);
0057 } else {
0058 DetId::Detector det;
0059 if (name == "HGCalHESiliconSensitive")
0060 det = DetId::HGCalHSi;
0061 else if (name == "HGCalHEScintillatorSensitive")
0062 det = DetId::HGCalHSc;
0063 else
0064 det = DetId::HGCalEE;
0065 edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name << " Detector " << det << " Mode "
0066 << geom->topology().dddConstants().geomMode();
0067 if (name == "HGCalHEScintillatorSensitive") {
0068 doTestScint(geom, det);
0069 } else {
0070 doTestWafer(geom, det);
0071 }
0072 }
0073 }
0074
0075 void HGCalGeometryTester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0076 edm::ParameterSetDescription desc;
0077 desc.add<std::string>("Detector", "HGCalEESensitive");
0078 descriptions.add("hgcalGeometryTesterEE", desc);
0079 }
0080
0081 void HGCalGeometryTester::doTest(const HGCalGeometry* geom, ForwardSubdetector subdet) {
0082 const std::vector<DetId>& ids = geom->getValidDetIds();
0083 edm::LogVerbatim("HGCalGeomX") << "doTest: " << ids.size() << " valid ids for " << geom->cellElement();
0084
0085 int layers[] = {1, 5, 10};
0086 int zsides[] = {1, -1};
0087 int cells[] = {1, 51, 101};
0088 int wafers[] = {1, 101, 201, 301, 401};
0089 const int ismax(5);
0090 for (int zside : zsides) {
0091 for (int is = 0; is < ismax; ++is) {
0092 int sector = wafers[is];
0093 int type = geom->topology().dddConstants().waferTypeT(sector);
0094 if (type != 1)
0095 type = 0;
0096 for (int layer : layers) {
0097 for (int cell : cells) {
0098 DetId id1;
0099 id1 = static_cast<DetId>(HGCalDetId(subdet, zside, layer, type, sector, cell));
0100 if (geom->topology().valid(id1)) {
0101 auto icell1 = geom->getGeometry(id1);
0102 GlobalPoint global1 = geom->getPosition(id1);
0103 DetId idc1 = geom->getClosestCell(global1);
0104 GlobalPoint global2 = geom->getPosition(idc1);
0105 std::string cherr = (id1.rawId() != idc1.rawId()) ? " ***** ERROR *****" : "";
0106 edm::LogVerbatim("HGCalGeomX")
0107 << "DetId (" << subdet << ":" << zside << ":" << layer << ":" << sector << ":0:" << cell << ") Geom "
0108 << icell1 << " position (" << global1.x() << ", " << global1.y() << ", " << global1.z() << ") ids "
0109 << std::hex << id1.rawId() << ":" << idc1.rawId() << std::dec << ":" << HGCalDetId(id1) << ":"
0110 << HGCalDetId(idc1) << " new position (" << global2.x() << ", " << global2.y() << ", " << global2.z()
0111 << ") parameter[3] = " << icell1->param()[2] << ":" << icell1->param()[2] << cherr;
0112 std::vector<GlobalPoint> corners = geom->getCorners(idc1);
0113 std::ostringstream st1;
0114 st1 << corners.size() << " corners";
0115 for (auto const& cor : corners)
0116 st1 << " [" << cor.x() << "," << cor.y() << "," << cor.z() << "]";
0117 edm::LogVerbatim("HGCalGeomX") << st1.str();
0118 }
0119 }
0120 }
0121 }
0122 }
0123 }
0124
0125 void HGCalGeometryTester::doTestWafer(const HGCalGeometry* geom, DetId::Detector det) {
0126 const std::vector<DetId>& ids = geom->getValidDetIds();
0127 edm::LogVerbatim("HGCalGeomX") << "doTestWafer:: " << ids.size() << " valid ids for " << geom->cellElement();
0128 int layers[] = {1, 5, 10};
0129 int zsides[] = {1, -1};
0130 int cells[] = {1, 4, 7};
0131 int wafers[] = {7, 5, 3, -3, -5, -7};
0132 for (int zside : zsides) {
0133 for (int layer : layers) {
0134 for (int waferU : wafers) {
0135 for (int waferV : wafers) {
0136 int type = geom->topology().dddConstants().getTypeHex(layer, waferU, waferV);
0137 edm::LogVerbatim("HGCalGeomX") << "zside " << zside << " layer " << layer << " wafer " << waferU << ":"
0138 << waferV << " type " << type;
0139 for (int cellU : cells) {
0140 for (int cellV : cells) {
0141 edm::LogVerbatim("HGCalGeomX") << "det " << det << " cell " << cellU << ":" << cellV;
0142 DetId id1 = static_cast<DetId>(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0143 edm::LogVerbatim("HGCalGeomX") << HGCSiliconDetId(id1);
0144 if (geom->topology().valid(id1)) {
0145 auto icell1 = geom->getGeometry(id1);
0146 GlobalPoint global1 = geom->getPosition(id1);
0147 DetId idc1 = geom->getClosestCell(global1);
0148 GlobalPoint global2 = geom->getPosition(idc1);
0149 std::string cherr = (id1.rawId() != idc1.rawId()) ? " ***** ERROR *****" : "";
0150 edm::LogVerbatim("HGCalGeomX")
0151 << "DetId (" << det << ":" << zside << ":" << type << ":" << layer << ":" << waferU << ":" << waferV
0152 << ":" << cellU << ":" << cellV << ") Geom " << icell1 << " position (" << global1.x() << ", "
0153 << global1.y() << ", " << global1.z() << ") ids " << std::hex << id1.rawId() << ":" << idc1.rawId()
0154 << std::dec << ":" << HGCSiliconDetId(id1) << ":" << HGCSiliconDetId(idc1) << " new position ("
0155 << global2.x() << ", " << global2.y() << ", " << global2.z()
0156 << ") parameter[3] = " << icell1->param()[2] << ":" << icell1->param()[2] << cherr;
0157 std::vector<GlobalPoint> corners = geom->getCorners(idc1);
0158 std::ostringstream st1;
0159 st1 << corners.size() << " corners";
0160 for (auto const& cor : corners)
0161 st1 << " [" << cor.x() << "," << cor.y() << "," << cor.z() << "]";
0162 edm::LogVerbatim("HGCalGeomX") << st1.str();
0163 }
0164 }
0165 }
0166 }
0167 }
0168 }
0169 }
0170 }
0171
0172 void HGCalGeometryTester::doTestScint(const HGCalGeometry* geom, DetId::Detector det) {
0173 const std::vector<DetId>& ids = geom->getValidDetIds();
0174 edm::LogVerbatim("HGCalGeomX") << "doTestScint: " << ids.size() << " valid ids for " << geom->cellElement();
0175 int layers[] = {9, 14, 21};
0176 int zsides[] = {1, -1};
0177 int iphis[] = {1, 51, 101, 151, 201};
0178 int ietas[] = {11, 20, 29};
0179 for (int zside : zsides) {
0180 for (int layer : layers) {
0181 int type = geom->topology().dddConstants().getTypeTrap(layer);
0182 for (int ieta : ietas) {
0183 std::pair<int, int> typm = geom->topology().dddConstants().tileType(layer, ieta, 0);
0184 for (int iphi : iphis) {
0185 HGCScintillatorDetId detId(type, layer, zside * ieta, iphi);
0186 if (typm.first >= 0) {
0187 detId.setType(typm.first);
0188 detId.setSiPM(typm.second);
0189 }
0190 DetId id1 = static_cast<DetId>(detId);
0191 if (geom->topology().valid(id1)) {
0192 auto icell1 = geom->getGeometry(id1);
0193 GlobalPoint global1 = geom->getPosition(id1);
0194 DetId idc1 = geom->getClosestCell(global1);
0195 GlobalPoint global2 = geom->getPosition(idc1);
0196 std::string cherr = (id1.rawId() != idc1.rawId()) ? " ***** ERROR *****" : "";
0197 edm::LogVerbatim("HGCalGeomX")
0198 << "DetId (" << det << ":" << zside << ":" << type << ":" << layer << ":" << ieta << ":" << iphi
0199 << ") Geom " << icell1 << " position (" << global1.x() << ", " << global1.y() << ", " << global1.z()
0200 << ":" << global1.perp() << ") ids " << std::hex << id1.rawId() << ":" << idc1.rawId() << std::dec
0201 << ":" << HGCScintillatorDetId(id1) << ":" << HGCScintillatorDetId(idc1) << " new position ("
0202 << global2.x() << ", " << global2.y() << ", " << global2.z() << ":" << global2.perp()
0203 << ") parameter[11] = " << icell1->param()[10] << ":" << icell1->param()[10] << cherr;
0204 std::vector<GlobalPoint> corners = geom->getCorners(idc1);
0205 std::ostringstream st1;
0206 st1 << corners.size() << " corners";
0207 for (auto const& cor : corners)
0208 st1 << " [" << cor.x() << "," << cor.y() << "," << cor.z() << "]";
0209 edm::LogVerbatim("HGCalGeomX") << st1.str();
0210 }
0211 }
0212 }
0213 }
0214 }
0215 for (int layer = geom->topology().dddConstants().firstLayer();
0216 layer <= geom->topology().dddConstants().lastLayer(true);
0217 ++layer) {
0218 HGCScintillatorDetId id(geom->topology().dddConstants().getTypeTrap(layer), layer, 20, 1);
0219 GlobalPoint global1 = geom->getPosition(id);
0220 edm::LogVerbatim("HGCalGeomX") << "Layer " << layer << " DetId " << id << " position (" << global1.x() << ", "
0221 << global1.y() << ", " << global1.z() << ", " << global1.perp() << ")";
0222 }
0223 }
0224
0225
0226 DEFINE_FWK_MODULE(HGCalGeometryTester);