File indexing completed on 2023-03-17 13:03:42
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
0047 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0048 const CaloGeometry* geom_;
0049 const int mode_;
0050 int layerEE_, layerFH_, layerBH_;
0051 int layerEE1000_, layerFH1000_, layerBH1000_;
0052 };
0053
0054 HGCalTestRecHitTool::HGCalTestRecHitTool(const edm::ParameterSet& iC)
0055 : geomToken_{esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})},
0056 geom_(nullptr),
0057 mode_(iC.getParameter<int>("Mode")) {
0058 layerEE_ = layerFH_ = layerBH_ = 0;
0059 layerEE1000_ = layerFH1000_ = layerBH1000_ = 0;
0060 }
0061
0062 HGCalTestRecHitTool::~HGCalTestRecHitTool() {}
0063
0064 void HGCalTestRecHitTool::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0065 if (auto pG = iSetup.getHandle(geomToken_)) {
0066 geom_ = pG.product();
0067 auto geomEE = ((mode_ == 0) ? static_cast<const HGCalGeometry*>(
0068 geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HGCEE))
0069 : static_cast<const HGCalGeometry*>(
0070 geom_->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty)));
0071 layerEE_ = (geomEE->topology().dddConstants()).layers(true);
0072 layerEE1000_ = (geomEE->topology().dddConstants()).getLayer(10000., true);
0073 auto geomFH = ((mode_ == 0) ? static_cast<const HGCalGeometry*>(
0074 geom_->getSubdetectorGeometry(DetId::Forward, ForwardSubdetector::HGCHEF))
0075 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0076 DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty)));
0077 layerFH_ = (geomFH->topology().dddConstants()).layers(true);
0078 layerFH1000_ = (geomFH->topology().dddConstants()).getLayer(10000., true);
0079 if (mode_ == 0) {
0080 auto geomBH =
0081 static_cast<const HcalGeometry*>(geom_->getSubdetectorGeometry(DetId::Hcal, HcalSubdetector::HcalEndcap));
0082 layerBH_ = (geomBH->topology().dddConstants())->getMaxDepth(1);
0083 } else {
0084 auto geomBH = static_cast<const HGCalGeometry*>(
0085 geom_->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0086 layerBH_ = (geomBH->topology().dddConstants()).layers(true);
0087 layerBH1000_ = (geomBH->topology().dddConstants()).getLayer(10000., true);
0088 }
0089 edm::LogVerbatim("HGCalGeom") << "Layers " << layerEE_ << ":" << layerFH_ << ":" << layerBH_
0090 << "\nLayer # at 1000 cm " << layerEE1000_ << ":" << layerFH1000_ << ":"
0091 << layerBH1000_;
0092 for (int layer = 1; layer <= layerEE_; ++layer)
0093 edm::LogVerbatim("HGCalGeom") << "EE Layer " << layer << " Wafers "
0094 << (geomEE->topology().dddConstants()).wafers(layer, 0) << ":"
0095 << (geomEE->topology().dddConstants()).wafers(layer, 1) << ":"
0096 << (geomEE->topology().dddConstants()).wafers(layer, 2) << std::endl;
0097 for (int layer = 1; layer <= layerFH_; ++layer)
0098 edm::LogVerbatim("HGCalGeom") << "FH Layer " << layer << " Wafers "
0099 << (geomFH->topology().dddConstants()).wafers(layer, 0) << ":"
0100 << (geomFH->topology().dddConstants()).wafers(layer, 1) << ":"
0101 << (geomFH->topology().dddConstants()).wafers(layer, 2) << std::endl;
0102 int nlayer = ((mode_ == 0) ? (layerEE_ + layerFH_ + layerBH_) : (layerEE_ + layerFH_));
0103 retrieveLayerPositions(nlayer);
0104 } else {
0105 edm::LogVerbatim("HGCalGeom") << "Cannot get valid CaloGeometry Object" << std::endl;
0106 }
0107 }
0108
0109 template <typename GEOM>
0110 void HGCalTestRecHitTool::check_geom(const GEOM* geom) const {
0111 if (nullptr == geom) {
0112 throw cms::Exception("HGCalTestRecHitTools") << "Geometry not provided yet to HGCalTestRecHitTools!";
0113 }
0114 }
0115
0116 template <typename DDD>
0117 void HGCalTestRecHitTool::check_ddd(const DDD* ddd) const {
0118 if (nullptr == ddd) {
0119 throw cms::Exception("HGCalTestRecHitTools") << "DDDConstants not accessibl to HGCalTestRecHitTools!";
0120 }
0121 }
0122
0123 std::vector<double> HGCalTestRecHitTool::retrieveLayerPositions(unsigned layers) {
0124 DetId id;
0125 std::vector<double> layerPositions(layers);
0126 for (int ilayer = 1; ilayer <= (int)(layers); ++ilayer) {
0127 int lay(ilayer), type(0);
0128 if (ilayer <= layerEE_) {
0129 id = ((mode_ == 0) ? static_cast<DetId>(HGCalDetId(ForwardSubdetector::HGCEE, 1, lay, 1, 2, 1))
0130 : static_cast<DetId>(HGCSiliconDetId(DetId::HGCalEE, 1, 0, lay, 3, 3, 1, 1)));
0131 } else if (ilayer > layerEE_ && ilayer <= (layerEE_ + layerFH_)) {
0132 lay = ilayer - layerEE_;
0133 id = ((mode_ == 0) ? static_cast<DetId>(HGCalDetId(ForwardSubdetector::HGCHEF, 1, lay, 1, 2, 1))
0134 : static_cast<DetId>(HGCSiliconDetId(DetId::HGCalHSi, 1, 0, lay, 3, 3, 1, 1)));
0135 type = 1;
0136 } else {
0137 lay = ilayer - (layerEE_ + layerFH_);
0138 id = HcalDetId(HcalSubdetector::HcalEndcap, 50, 100, lay);
0139 type = 2;
0140 }
0141 const GlobalPoint pos = getPosition(id);
0142 if (id.det() == DetId::Hcal) {
0143 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HcalDetId(id);
0144 } else if (id.det() == DetId::Forward) {
0145 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HGCalDetId(id);
0146 } else if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
0147 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HGCSiliconDetId(id);
0148 } else {
0149 edm::LogVerbatim("HGCalGeom") << "GEOM layer " << ilayer << " ID " << HGCScintillatorDetId(id);
0150 }
0151 edm::LogVerbatim("HGCalGeom") << " Z " << pos.z() << ":" << getLayerZ(id) << ":" << getLayerZ(type, lay)
0152 << std::endl;
0153 layerPositions[ilayer - 1] = pos.z();
0154 }
0155 return layerPositions;
0156 }
0157
0158 const HcalDDDRecConstants* HGCalTestRecHitTool::get_ddd(const HcalDetId& detid) const {
0159 auto geom = geom_->getSubdetectorGeometry(detid);
0160 check_geom(geom);
0161 const HcalGeometry* hc = static_cast<const HcalGeometry*>(geom);
0162 const HcalDDDRecConstants* ddd = hc->topology().dddConstants();
0163 check_ddd(ddd);
0164 return ddd;
0165 }
0166
0167 const HGCalDDDConstants* HGCalTestRecHitTool::get_ddd(const DetId& detid) const {
0168 auto geom = geom_->getSubdetectorGeometry(detid);
0169 check_geom(geom);
0170 const HGCalGeometry* hg = static_cast<const HGCalGeometry*>(geom);
0171 const HGCalDDDConstants* ddd = &(hg->topology().dddConstants());
0172 check_ddd(ddd);
0173 return ddd;
0174 }
0175
0176 GlobalPoint HGCalTestRecHitTool::getPosition(const DetId& id) const {
0177 auto geom = geom_->getSubdetectorGeometry(id);
0178 check_geom(geom);
0179 GlobalPoint position;
0180 if (id.det() == DetId::Hcal) {
0181 position = geom->getGeometry(id)->getPosition();
0182 } else {
0183 position = (dynamic_cast<const HGCalGeometry*>(geom))->getPosition(id);
0184 }
0185 return position;
0186 }
0187
0188 double HGCalTestRecHitTool::getLayerZ(const DetId& id) const {
0189 double zpos(0);
0190 if (id.det() == DetId::Hcal) {
0191 auto geom = geom_->getSubdetectorGeometry(id);
0192 check_geom(geom);
0193 zpos = (static_cast<const HcalGeometry*>(geom))->getGeometry(id)->getPosition().z();
0194 } else {
0195 const HGCalDDDConstants* ddd = get_ddd(id);
0196 int layer = ((id.det() == DetId::Forward) ? HGCalDetId(id).layer()
0197 : ((id.det() != DetId::HGCalHSc) ? HGCSiliconDetId(id).layer()
0198 : HGCScintillatorDetId(id).layer()));
0199 zpos = ddd->waferZ(layer, true);
0200 }
0201 return zpos;
0202 }
0203
0204 double HGCalTestRecHitTool::getLayerZ(int type, int layer) const {
0205 double zpos(0);
0206 if (type == 2) {
0207 auto geom =
0208 static_cast<const HcalGeometry*>(geom_->getSubdetectorGeometry(DetId::Hcal, HcalSubdetector::HcalEndcap));
0209 check_geom(geom);
0210 auto ddd = (geom->topology().dddConstants());
0211 check_ddd(ddd);
0212 std::pair<int, int> etar = ddd->getEtaRange(1);
0213 zpos = ddd->getRZ(2, etar.second, layer);
0214 } else {
0215 auto geom = ((type == 1) ? ((mode_ == 0) ? static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0216 DetId::Forward, ForwardSubdetector::HGCHEF))
0217 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0218 DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty)))
0219 : ((mode_ == 0) ? static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0220 DetId::Forward, ForwardSubdetector::HGCEE))
0221 : static_cast<const HGCalGeometry*>(geom_->getSubdetectorGeometry(
0222 DetId::HGCalEE, ForwardSubdetector::ForwardEmpty))));
0223 check_geom(geom);
0224 auto ddd = &(geom->topology().dddConstants());
0225 check_ddd(ddd);
0226 zpos = ddd->waferZ(layer, true);
0227 }
0228 return zpos;
0229 }
0230
0231
0232 DEFINE_FWK_MODULE(HGCalTestRecHitTool);
