File indexing completed on 2025-06-09 22:19:19
0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "DataFormats/DetId/interface/DetId.h"
0015 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0017 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0018 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0019 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0020 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0021 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0022 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0023 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0026 #include "CoralBase/Exception.h"
0027
0028 class HGCalTestRecHitTool : public edm::one::EDAnalyzer<> {
0029 public:
0030 explicit HGCalTestRecHitTool(const edm::ParameterSet&);
0031 ~HGCalTestRecHitTool() override;
0032
0033 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0034
0035 private:
0036 std::vector<double> retrieveLayerPositions(unsigned int);
0037 template <typename GEOM>
0038 void check_geom(const GEOM* geom) const;
0039 template <typename DDD>
0040 void check_ddd(const DDD* ddd) const;
0041 const HGCalDDDConstants* get_ddd(const DetId& detid) const;
0042 const HcalDDDRecConstants* get_ddd(const HcalDetId& detid) const;
0043 double getLayerZ(DetId const&) const;
0044 double getLayerZ(int type, int layer) const;
0045 GlobalPoint getPosition(DetId const&) const;
0046 int getScintMaxIphi(const DetId& id) const;
0047 bool isScintillatorFine(const DetId& id) const;
0048
0049 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0050 const CaloGeometry* geom_;
0051 const int mode_;
0052 int layerEE_, layerFH_, layerBH_;
0053 int eeOffset_, fhOffset_, bhOffset_;
0054 int layerEE1000_, layerFH1000_, layerBH1000_;
0055 };
0056
0057 HGCalTestRecHitTool::HGCalTestRecHitTool(const edm::ParameterSet& iC)
0058 : geomToken_{esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})},
0059 geom_(nullptr),
0060 mode_(iC.getParameter<int>("Mode")) {
0061 layerEE_ = layerFH_ = layerBH_ = 0;
0062 eeOffset_ = fhOffset_ = bhOffset_ = 0;
0063 layerEE1000_ = layerFH1000_ = layerBH1000_ = 0;
0064 edm::LogVerbatim("HGCalGeom") << "Instantiate HGCalTestRecHitTool with mode " << mode_;
0065 }
0066
0067 HGCalTestRecHitTool::~HGCalTestRecHitTool() {}
0068
0069 void HGCalTestRecHitTool::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0070 if (auto pG = iSetup.getHandle(geomToken_)) {
0071 geom_ = pG.product();
0072 auto geomEE = ((mode_ == 0) ? static_cast<const HGCalGeometry*>(
0073 geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HGCEE))
0074 : static_cast<const HGCalGeometry*>(
0075 geom_->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty)));
0076 layerEE_ = (geomEE->topology().dddConstants()).layers(true);
0077 eeOffset_ = (geomEE->topology().dddConstants()).getLayerOffset();
0078 layerEE1000_ = (geomEE->topology().dddConstants()).getLayer(10000., true);
0079 edm::LogVerbatim("HGCalGeom") << "EE::Layers " << layerEE_ << " Offset " << eeOffset_ << " Layer # at 1000 cm "
0080 << layerEE1000_;
0081 auto geomFH = ((mode_ == 0) ? static_cast<const HGCalGeometry*>(
0082 geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HGCHEF))
0083 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0084 DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty)));
0085 layerFH_ = (geomFH->topology().dddConstants()).layers(true);
0086 fhOffset_ = (geomFH->topology().dddConstants()).getLayerOffset();
0087 layerFH1000_ = (geomFH->topology().dddConstants()).getLayer(10000., true);
0088 edm::LogVerbatim("HGCalGeom") << "FH::Layers " << layerFH_ << " Offsets " << fhOffset_ << " Layer # at 1000 cm "
0089 << layerFH1000_;
0090 if (mode_ == 0) {
0091 auto geomBH =
0092 static_cast<const HcalGeometry*>(geom_->getSubdetectorGeometry(DetId::Hcal, HcalSubdetector::HcalEndcap));
0093 layerBH_ = (geomBH->topology().dddConstants())->getMaxDepth(1);
0094 } else {
0095 auto geomBH = static_cast<const HGCalGeometry*>(
0096 geom_->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0097 layerBH_ = (geomBH->topology().dddConstants()).layers(true);
0098 bhOffset_ = (geomBH->topology().dddConstants()).getLayerOffset();
0099 layerBH1000_ = (geomBH->topology().dddConstants()).getLayer(10000., true);
0100 }
0101 edm::LogVerbatim("HGCalGeom") << "BH::Layers " << layerBH_ << " nOffsets " << bhOffset_ << " Layer # at 1000 cm "
0102 << layerBH1000_;
0103 for (int layer = 1; layer <= layerEE_; ++layer)
0104 edm::LogVerbatim("HGCalGeom") << "EE Layer " << layer << " Wafers "
0105 << (geomEE->topology().dddConstants()).wafers(layer, 0) << ":"
0106 << (geomEE->topology().dddConstants()).wafers(layer, 1) << ":"
0107 << (geomEE->topology().dddConstants()).wafers(layer, 2);
0108 for (int layer = 1; layer <= layerFH_; ++layer)
0109 edm::LogVerbatim("HGCalGeom") << "FH Layer " << layer << " Wafers "
0110 << (geomFH->topology().dddConstants()).wafers(layer, 0) << ":"
0111 << (geomFH->topology().dddConstants()).wafers(layer, 1) << ":"
0112 << (geomFH->topology().dddConstants()).wafers(layer, 2);
0113 if (mode_ != 0) {
0114 auto geomBH = static_cast<const HGCalGeometry*>(
0115 geom_->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0116 int firstL = (geomBH->topology().dddConstants()).firstLayer();
0117 edm::LogVerbatim("HGCalGeom") << "BH First Layer " << firstL << " Total " << layerBH_;
0118 for (int lay = 1; lay <= layerBH_; ++lay) {
0119 int layer = firstL + lay - 1;
0120 int ring = (geomBH->topology().dddConstants()).tileRings(layer).first + 1;
0121 auto typm = (geomBH->topology().dddConstants()).tileType(layer, ring, 1);
0122 HGCScintillatorDetId id(typm.first, layer, ring, 1, false, typm.second, 0);
0123 edm::LogVerbatim("HGCalGeom") << "BH Layer " << layer << " Ring " << ring << " Max phi "
0124 << getScintMaxIphi(DetId(id)) << " Fine " << isScintillatorFine(DetId(id));
0125 }
0126 }
0127 int nlayer = ((mode_ == 0) ? (layerEE_ + layerFH_ + layerBH_) : (layerEE_ + layerFH_));
0128 retrieveLayerPositions(nlayer);
0129 } else {
0130 edm::LogVerbatim("HGCalGeom") << "Cannot get valid CaloGeometry Object" << std::endl;
0131 }
0132 }
0133
0134 template <typename GEOM>
0135 void HGCalTestRecHitTool::check_geom(const GEOM* geom) const {
0136 if (nullptr == geom) {
0137 throw cms::Exception("HGCalTestRecHitTools") << "Geometry not provided yet to HGCalTestRecHitTools!";
0138 }
0139 }
0140
0141 template <typename DDD>
0142 void HGCalTestRecHitTool::check_ddd(const DDD* ddd) const {
0143 if (nullptr == ddd) {
0144 throw cms::Exception("HGCalTestRecHitTools") << "DDDConstants not accessibl to HGCalTestRecHitTools!";
0145 }
0146 }
0147
0148 std::vector<double> HGCalTestRecHitTool::retrieveLayerPositions(unsigned layers) {
0149 DetId id;
0150 std::vector<double> layerPositions(layers);
0151 for (int ilayer = 1; ilayer <= (int)(layers); ++ilayer) {
0152 int lay(ilayer), type(0);
0153 if (ilayer <= layerEE_) {
0154 id = ((mode_ == 0) ? static_cast<DetId>(HGCalDetId(ForwardSubdetector::HGCEE, 1, lay, 1, 2, 1))
0155 : static_cast<DetId>(HGCSiliconDetId(DetId::HGCalEE, 1, 0, lay, 3, 3, 1, 1)));
0156 } else if (ilayer > layerEE_ && ilayer <= (layerEE_ + layerFH_)) {
0157 lay = ilayer - layerEE_;
0158 id = ((mode_ == 0) ? static_cast<DetId>(HGCalDetId(ForwardSubdetector::HGCHEF, 1, lay, 1, 2, 1))
0159 : static_cast<DetId>(HGCSiliconDetId(DetId::HGCalHSi, 1, 0, lay, 3, 3, 1, 1)));
0160 type = 1;
0161 } else {
0162 lay = ilayer - (layerEE_ + layerFH_);
0163 id = HcalDetId(HcalSubdetector::HcalEndcap, 50, 100, lay);
0164 type = 2;
0165 }
0166 const GlobalPoint pos = getPosition(id);
0167 if (id.det() == DetId::Hcal) {
0168 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HcalDetId(id);
0169 } else if (id.det() == DetId::Forward) {
0170 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HGCalDetId(id);
0171 } else if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
0172 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HGCSiliconDetId(id);
0173 } else {
0174 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HGCScintillatorDetId(id);
0175 }
0176 edm::LogVerbatim("HGCalGeom") << " Z " << pos.z() << ":" << getLayerZ(id) << ":" << getLayerZ(type, lay)
0177 << std::endl;
0178 layerPositions[ilayer - 1] = pos.z();
0179 }
0180 return layerPositions;
0181 }
0182
0183 const HcalDDDRecConstants* HGCalTestRecHitTool::get_ddd(const HcalDetId& detid) const {
0184 auto geom = geom_->getSubdetectorGeometry(detid);
0185 check_geom(geom);
0186 const HcalGeometry* hc = static_cast<const HcalGeometry*>(geom);
0187 const HcalDDDRecConstants* ddd = hc->topology().dddConstants();
0188 check_ddd(ddd);
0189 return ddd;
0190 }
0191
0192 const HGCalDDDConstants* HGCalTestRecHitTool::get_ddd(const DetId& detid) const {
0193 auto geom = geom_->getSubdetectorGeometry(detid);
0194 check_geom(geom);
0195 const HGCalGeometry* hg = static_cast<const HGCalGeometry*>(geom);
0196 const HGCalDDDConstants* ddd = &(hg->topology().dddConstants());
0197 check_ddd(ddd);
0198 return ddd;
0199 }
0200
0201 GlobalPoint HGCalTestRecHitTool::getPosition(const DetId& id) const {
0202 auto geom = geom_->getSubdetectorGeometry(id);
0203 check_geom(geom);
0204 GlobalPoint position;
0205 if (id.det() == DetId::Hcal) {
0206 position = geom->getGeometry(id)->getPosition();
0207 } else {
0208 position = (dynamic_cast<const HGCalGeometry*>(geom))->getPosition(id);
0209 }
0210 return position;
0211 }
0212
0213 double HGCalTestRecHitTool::getLayerZ(const DetId& id) const {
0214 double zpos(0);
0215 if (id.det() == DetId::Hcal) {
0216 auto geom = geom_->getSubdetectorGeometry(id);
0217 check_geom(geom);
0218 zpos = (static_cast<const HcalGeometry*>(geom))->getGeometry(id)->getPosition().z();
0219 } else {
0220 const HGCalDDDConstants* ddd = get_ddd(id);
0221 int layer = ((id.det() == DetId::Forward) ? HGCalDetId(id).layer()
0222 : ((id.det() != DetId::HGCalHSc) ? HGCSiliconDetId(id).layer()
0223 : HGCScintillatorDetId(id).layer()));
0224 zpos = ddd->waferZ(layer, true);
0225 }
0226 return zpos;
0227 }
0228
0229 double HGCalTestRecHitTool::getLayerZ(int type, int layer) const {
0230 double zpos(0);
0231 if (type == 2) {
0232 auto geom =
0233 static_cast<const HcalGeometry*>(geom_->getSubdetectorGeometry(DetId::Hcal, HcalSubdetector::HcalEndcap));
0234 check_geom(geom);
0235 auto ddd = (geom->topology().dddConstants());
0236 check_ddd(ddd);
0237 std::pair<int, int> etar = ddd->getEtaRange(1);
0238 zpos = ddd->getRZ(2, etar.second, layer);
0239 } else {
0240 auto geom = ((type == 1) ? ((mode_ == 0) ? static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0241 DetId::Forward, ForwardSubdetector::HGCHEF))
0242 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0243 DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty)))
0244 : ((mode_ == 0) ? static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0245 DetId::Forward, ForwardSubdetector::HGCEE))
0246 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0247 DetId::HGCalEE, ForwardSubdetector::ForwardEmpty))));
0248 check_geom(geom);
0249 auto ddd = &(geom->topology().dddConstants());
0250 check_ddd(ddd);
0251 zpos = ddd->waferZ(layer, true);
0252 }
0253 return zpos;
0254 }
0255
0256 int HGCalTestRecHitTool::getScintMaxIphi(const DetId& id) const {
0257 if (id.det() == DetId::HGCalHSc) {
0258 int layer = HGCScintillatorDetId(id).layer();
0259 auto hg = static_cast<const HGCalGeometry*>(
0260 geom_->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0261 return hg->topology().dddConstants().maxCells(layer, true);
0262 } else {
0263 return 0;
0264 }
0265 }
0266
0267 bool HGCalTestRecHitTool::isScintillatorFine(const DetId& id) const {
0268 edm::LogVerbatim("HGCalGeom") << "isScintillatorFine " << id.det() << ":" << DetId::HGCalHSc << ":" << bhOffset_;
0269 if (id.det() == DetId::HGCalHSc) {
0270 int layer = HGCScintillatorDetId(id).layer();
0271 auto hg = static_cast<const HGCalGeometry*>(
0272 geom_->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0273 return hg->topology().dddConstants().scintFine(layer);
0274 } else {
0275 return false;
0276 }
0277 }
0278
0279
0280 DEFINE_FWK_MODULE(HGCalTestRecHitTool);
