File indexing completed on 2024-04-06 12:14:51
0001 #include "Geometry/HcalTowerAlgo/interface/HcalDDDGeometryLoader.h"
0002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0003 #include "Geometry/HcalCommonData/interface/HcalParametersFromDD.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
0006 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include <string>
0010
0011 typedef CaloCellGeometry::CCGFloat CCGFloat;
0012
0013
0014
0015 HcalDDDGeometryLoader::HcalDDDGeometryLoader(const HcalDDDRecConstants* hcons) : hcalConstants_(hcons) {
0016 isBH_ = hcalConstants_->isBH();
0017 }
0018
0019 HcalDDDGeometryLoader::~HcalDDDGeometryLoader() {}
0020
0021 HcalDDDGeometryLoader::ReturnType HcalDDDGeometryLoader::load(const HcalTopology& topo,
0022 DetId::Detector det,
0023 int subdet) {
0024 HcalSubdetector hsub = static_cast<HcalSubdetector>(subdet);
0025
0026 HcalDDDGeometry* geom(new HcalDDDGeometry(topo));
0027
0028 if (geom->cornersMgr() == nullptr) {
0029 const unsigned int count(hcalConstants_->numberOfCells(HcalBarrel) + hcalConstants_->numberOfCells(HcalEndcap) +
0030 2 * hcalConstants_->numberOfCells(HcalForward) + hcalConstants_->numberOfCells(HcalOuter));
0031 geom->allocateCorners(count);
0032 }
0033
0034
0035
0036 if (geom->parMgr() == nullptr)
0037 geom->allocatePar(500, 3);
0038
0039 fill(hsub, geom);
0040
0041 geom->sortValidIds();
0042 return geom;
0043 }
0044
0045 HcalDDDGeometryLoader::ReturnType HcalDDDGeometryLoader::load(const HcalTopology& topo) {
0046 HcalDDDGeometry* geom(new HcalDDDGeometry(topo));
0047
0048 if (geom->cornersMgr() == nullptr) {
0049 const unsigned int count(hcalConstants_->numberOfCells(HcalBarrel) + hcalConstants_->numberOfCells(HcalEndcap) +
0050 2 * hcalConstants_->numberOfCells(HcalForward) + hcalConstants_->numberOfCells(HcalOuter));
0051 geom->allocateCorners(count);
0052 }
0053 if (geom->parMgr() == nullptr)
0054 geom->allocatePar(500, 3);
0055
0056 fill(HcalBarrel, geom);
0057 fill(HcalEndcap, geom);
0058 fill(HcalForward, geom);
0059 fill(HcalOuter, geom);
0060
0061 geom->sortValidIds();
0062 return geom;
0063 }
0064
0065 void HcalDDDGeometryLoader::fill(HcalSubdetector subdet, HcalDDDGeometry* geom) {
0066
0067 std::vector<HcalCellType> hcalCells = hcalConstants_->HcalCellTypes(subdet);
0068 geom->insertCell(hcalCells);
0069 #ifdef EDM_ML_DEBUG
0070 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::fill gets " << hcalCells.size() << " cells for subdetector "
0071 << subdet;
0072 #endif
0073
0074
0075 std::vector<HcalDetId> hcalIds;
0076 for (auto& hcalCell : hcalCells) {
0077 int etaRing = hcalCell.etaBin();
0078 int iside = hcalCell.zside();
0079 int depthBin = hcalCell.depthSegment();
0080 double dphi = hcalCell.phiBinWidth();
0081 std::vector<std::pair<int, double> > phis = hcalCell.phis();
0082 #ifdef EDM_ML_DEBUG
0083 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader: Subdet " << subdet << " side " << iside << " eta "
0084 << etaRing << " depth " << depthBin << " with " << phis.size() << "modules:";
0085 size_t i(0);
0086 #endif
0087 geom->increaseReserve(phis.size());
0088 for (auto& phi : phis) {
0089 #ifdef EDM_ML_DEBUG
0090 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::fill Cell " << i << " eta " << iside * etaRing << " phi "
0091 << phi.first << "(" << phi.second / CLHEP::deg << ", " << dphi / CLHEP::deg
0092 << ") depth " << depthBin;
0093 i++;
0094 #endif
0095 HcalDetId id(subdet, iside * etaRing, phi.first, depthBin);
0096 hcalIds.emplace_back(id);
0097 makeCell(id, hcalCell, phi.second, dphi, geom);
0098 }
0099 }
0100
0101 edm::LogInfo("HCalGeom") << "Number of HCAL DetIds made for " << subdet << " is " << hcalIds.size();
0102 }
0103
0104 void HcalDDDGeometryLoader::makeCell(
0105 const HcalDetId& detId, const HcalCellType& hcalCell, double phi, double dphi, HcalDDDGeometry* geom) const {
0106
0107 double eta1 = hcalCell.etaMin();
0108 double eta2 = hcalCell.etaMax();
0109 HcalSubdetector subdet = detId.subdet();
0110 double eta = 0.5 * (eta1 + eta2) * detId.zside();
0111 double deta = (eta2 - eta1);
0112 double theta = 2.0 * atan(exp(-eta));
0113
0114
0115 bool rzType = hcalCell.depthType();
0116 bool isBarrel = (subdet == HcalBarrel || subdet == HcalOuter);
0117
0118 double z, r, thickness;
0119 #ifdef EDM_ML_DEBUG
0120 double r0, r1, r2;
0121 #endif
0122 if (rzType) {
0123 r = hcalCell.depthMin();
0124 if (isBarrel) {
0125 z = r * sinh(eta);
0126 thickness = (hcalCell.depthMax() - r) * cosh(eta);
0127 #ifdef EDM_ML_DEBUG
0128 r1 = r;
0129 r2 = hcalCell.depthMax();
0130 r0 = 0.5 * (r1 + r2);
0131 #endif
0132 } else {
0133 z = r * sinh(eta2);
0134 thickness = 2. * hcalCell.halfSize();
0135 r = z / sinh(std::abs(eta));
0136 #ifdef EDM_ML_DEBUG
0137 r0 = z / sinh(std::abs(eta));
0138 r1 = z / sinh(std::abs(eta) + 0.5 * deta);
0139 r2 = z / sinh(std::abs(eta) - 0.5 * deta);
0140 #endif
0141 }
0142 #ifdef EDM_ML_DEBUG
0143 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::makeCell SubDet " << subdet << " eta = " << eta
0144 << " theta = " << theta << " r = " << r << " thickness = " << thickness << " r0-r2 ("
0145 << r0 << ":" << r1 << ":" << r2 << ")";
0146 #endif
0147 } else {
0148 z = hcalCell.depthMin();
0149 thickness = hcalCell.depthMax() - z;
0150 if (isBH_)
0151 z += (0.5 * thickness);
0152 z *= detId.zside();
0153 r = z * tan(theta);
0154 thickness /= std::abs(cos(theta));
0155 #ifdef EDM_ML_DEBUG
0156 r0 = z / sinh(std::abs(eta));
0157 r1 = z / sinh(std::abs(eta) + 0.5 * deta);
0158 r2 = z / sinh(std::abs(eta) - 0.5 * deta);
0159 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::makeCell SubDet " << subdet << " eta = " << eta
0160 << " theta = " << theta << " z = " << z << " r = " << r << " thickness = " << thickness
0161 << " r0-r2 (" << r0 << ":" << r1 << ":" << r2 << ")";
0162 #endif
0163 }
0164
0165 double x = r * cos(phi);
0166 double y = r * sin(phi);
0167 GlobalPoint point(x, y, z);
0168
0169 #ifdef EDM_ML_DEBUG
0170 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::makeCell for " << detId << " Point " << point
0171 << " deta = " << deta << " dphi = " << dphi << " thickness = " << thickness
0172 << " isBarrel = " << isBarrel << " " << rzType;
0173 #endif
0174
0175 std::vector<CCGFloat> hp;
0176 hp.reserve(3);
0177
0178 if (subdet == HcalForward) {
0179 hp.emplace_back(deta / 2.);
0180 hp.emplace_back(dphi / 2.);
0181 hp.emplace_back(thickness / 2.);
0182 } else {
0183 const double sign(isBarrel ? 1 : -1);
0184 hp.emplace_back(deta / 2.);
0185 hp.emplace_back(dphi / 2.);
0186 hp.emplace_back(sign * thickness / 2.);
0187 }
0188 geom->newCellFast(point, point, point, CaloCellGeometry::getParmPtr(hp, geom->parMgr(), geom->parVecVec()), detId);
0189 }