File indexing completed on 2024-04-06 12:15:12
0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004
0005 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0006 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0008 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0010 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020
0021 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0022 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0025 #include "CoralBase/Exception.h"
0026
0027 class HGCalTestNeighbor : public edm::one::EDAnalyzer<> {
0028 public:
0029 explicit HGCalTestNeighbor(const edm::ParameterSet&);
0030 ~HGCalTestNeighbor() override = default;
0031
0032 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0033 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0034
0035 private:
0036 void doTest(const HGCalGeometry* geom, ForwardSubdetector subdet, const MagneticField* bField);
0037 void doTestWafer(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField);
0038 void doTestScint(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField);
0039
0040 const std::string name_;
0041 const std::vector<double> px_, py_, pz_;
0042 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> fieldToken_;
0043 const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0044 };
0045
0046 HGCalTestNeighbor::HGCalTestNeighbor(const edm::ParameterSet& iC)
0047 : name_(iC.getParameter<std::string>("detector")),
0048 px_(iC.getParameter<std::vector<double> >("pX")),
0049 py_(iC.getParameter<std::vector<double> >("pY")),
0050 pz_(iC.getParameter<std::vector<double> >("pZ")),
0051 fieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>(edm::ESInputTag{})),
0052 geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", name_})) {}
0053
0054 void HGCalTestNeighbor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0055 edm::ParameterSetDescription desc;
0056 desc.add<std::string>("detector", "HGCalEESensitive");
0057 std::vector<double> pxy = {2.0, 5.0, 10.0, 100.0};
0058 std::vector<double> pz = {10.0, 5.0, 20.0, 50.0};
0059 desc.add<std::vector<double> >("pX", pxy);
0060 desc.add<std::vector<double> >("pY", pxy);
0061 desc.add<std::vector<double> >("pZ", pz);
0062 descriptions.add("hgcalEETestNeighbor", desc);
0063 }
0064
0065 void HGCalTestNeighbor::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0066 const auto& bFieldR = iSetup.getData(fieldToken_);
0067 const MagneticField* bField = &bFieldR;
0068
0069 const auto& geomR = iSetup.getData(geomToken_);
0070 const HGCalGeometry* geom = &geomR;
0071 if (geom->topology().waferHexagon6()) {
0072 ForwardSubdetector subdet;
0073 if (name_ == "HGCalHESiliconSensitive")
0074 subdet = HGCHEF;
0075 else if (name_ == "HGCalHEScintillatorSensitive")
0076 subdet = HGCHEB;
0077 else
0078 subdet = HGCEE;
0079 edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name_ << " Detector:Subdetector " << DetId::Forward << ":"
0080 << subdet;
0081 doTest(geom, subdet, bField);
0082 } else {
0083 DetId::Detector det;
0084 if (name_ == "HGCalHESiliconSensitive")
0085 det = DetId::HGCalHSi;
0086 else if (name_ == "HGCalHEScintillatorSensitive")
0087 det = DetId::HGCalHSc;
0088 else
0089 det = DetId::HGCalEE;
0090 edm::LogVerbatim("HGCalGeomX") << "Perform test for " << name_ << " Detector " << det;
0091 if (name_ == "HGCalHEScintillatorSensitive") {
0092 doTestScint(geom, det, bField);
0093 } else {
0094 doTestWafer(geom, det, bField);
0095 }
0096 }
0097 }
0098
0099 void HGCalTestNeighbor::doTest(const HGCalGeometry* geom, ForwardSubdetector subdet, const MagneticField* bField) {
0100 const std::vector<DetId>& ids = geom->getValidDetIds();
0101 edm::LogVerbatim("HGCalGeomX") << "doTest: " << ids.size() << " valid ids for " << geom->cellElement();
0102
0103 int layers[] = {1, 5, 10};
0104 int zsides[] = {1, -1};
0105 int cells[] = {1, 51, 101};
0106 int wafers[] = {1, 101, 201, 301, 401};
0107 const int ismax(5);
0108 for (int zside : zsides) {
0109 for (int is = 0; is < ismax; ++is) {
0110 int sector = wafers[is];
0111 int type = geom->topology().dddConstants().waferTypeT(sector);
0112 if (type != 1)
0113 type = 0;
0114 for (int layer : layers) {
0115 for (int cell : cells) {
0116 DetId id1;
0117 id1 = (DetId)(HGCalDetId(subdet, zside, layer, type, sector, cell));
0118 if (geom->topology().valid(id1)) {
0119 auto icell1 = geom->getGeometry(id1);
0120 GlobalPoint global1 = geom->getPosition(id1);
0121 for (unsigned int k = 0; k < px_.size(); ++k) {
0122 GlobalVector p(px_[k], py_[k], zside * pz_[k]);
0123 DetId id2 = geom->neighborZ(id1, p);
0124 DetId id3 = geom->neighborZ(id1, bField, 1, p);
0125 edm::LogVerbatim("HGCalGeomX")
0126 << "DetId" << HGCalDetId(id1) << " :" << global1 << " p" << p << " ID2" << HGCalDetId(id2) << " :"
0127 << geom->getPosition(id2) << " ID3" << HGCalDetId(id3) << " :" << geom->getPosition(id3);
0128 }
0129 }
0130 }
0131 }
0132 }
0133 }
0134 }
0135
0136 void HGCalTestNeighbor::doTestWafer(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField) {
0137 const std::vector<DetId>& ids = geom->getValidDetIds();
0138 edm::LogVerbatim("HGCalGeomX") << "doTestWafer:: " << ids.size() << " valid ids for " << geom->cellElement();
0139 int layers[] = {1, 5, 10};
0140 int zsides[] = {1, -1};
0141 int cells[] = {1, 4, 7};
0142 int wafers[] = {7, 5, 3, -3, -5, -7};
0143 for (int zside : zsides) {
0144 for (int layer : layers) {
0145 for (int waferU : wafers) {
0146 for (int waferV : wafers) {
0147 int type = geom->topology().dddConstants().getTypeHex(layer, waferU, waferV);
0148 for (int cellU : cells) {
0149 for (int cellV : cells) {
0150 DetId id1 = (DetId)(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0151 if (geom->topology().valid(id1)) {
0152 auto icell1 = geom->getGeometry(id1);
0153 GlobalPoint global1 = geom->getPosition(id1);
0154 for (unsigned int k = 0; k < px_.size(); ++k) {
0155 GlobalVector p(px_[k], py_[k], zside * pz_[k]);
0156 DetId id2 = geom->neighborZ(id1, p);
0157 DetId id3 = geom->neighborZ(id1, bField, 1, p);
0158 edm::LogVerbatim("HGCalGeomX") << "DetId" << HGCSiliconDetId(id1) << " :" << global1 << " p" << p
0159 << " ID2" << HGCSiliconDetId(id2) << " :" << geom->getPosition(id2)
0160 << " ID3" << HGCSiliconDetId(id3) << " :" << geom->getPosition(id3);
0161 }
0162 }
0163 }
0164 }
0165 }
0166 }
0167 }
0168 }
0169 }
0170
0171 void HGCalTestNeighbor::doTestScint(const HGCalGeometry* geom, DetId::Detector det, const MagneticField* bField) {
0172 const std::vector<DetId>& ids = geom->getValidDetIds();
0173 edm::LogVerbatim("HGCalGeomX") << "doTestScint: " << ids.size() << " valid ids for " << geom->cellElement();
0174 int layers[] = {9, 15, 22};
0175 int zsides[] = {1, -1};
0176 int iphis[] = {1, 51, 101, 151, 201};
0177 int ietas[] = {11, 20, 29};
0178 for (int zside : zsides) {
0179 for (int layer : layers) {
0180 int type = geom->topology().dddConstants().getTypeTrap(layer);
0181 for (int ieta : ietas) {
0182 for (int iphi : iphis) {
0183 DetId id1 = (DetId)(HGCScintillatorDetId(type, layer, zside * ieta, iphi));
0184 if (geom->topology().valid(id1)) {
0185 auto icell1 = geom->getGeometry(id1);
0186 GlobalPoint global1 = geom->getPosition(id1);
0187 for (unsigned int k = 0; k < px_.size(); ++k) {
0188 GlobalVector p(px_[k], py_[k], zside * pz_[k]);
0189 DetId id2 = geom->neighborZ(id1, p);
0190 DetId id3 = geom->neighborZ(id1, bField, 1, p);
0191 edm::LogVerbatim("HGCalGeomX") << "DetId" << HGCScintillatorDetId(id1) << " :" << global1 << " p" << p
0192 << " ID2" << HGCScintillatorDetId(id2) << " :" << geom->getPosition(id2)
0193 << " ID3" << HGCScintillatorDetId(id3) << " :" << geom->getPosition(id3);
0194 }
0195 }
0196 }
0197 }
0198 }
0199 }
0200 }
0201
0202
0203 DEFINE_FWK_MODULE(HGCalTestNeighbor);