File indexing completed on 2024-10-25 23:57:07
0001 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0002 #include "Geometry/HcalTowerAlgo/interface/HcalDDDGeometry.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include <algorithm>
0006 #include <mutex>
0007
0008
0009
0010 static std::mutex s_fillLock;
0011
0012 HcalDDDGeometry::HcalDDDGeometry(const HcalTopology& topo)
0013 : topo_(topo),
0014 etaMax_(0),
0015 m_hbCellVec(topo.getHBSize()),
0016 m_heCellVec(topo.getHESize()),
0017 m_hoCellVec(topo.getHOSize()),
0018 m_hfCellVec(topo.getHFSize()),
0019 m_filledDetIds(false) {}
0020
0021 HcalDDDGeometry::~HcalDDDGeometry() {}
0022
0023 void HcalDDDGeometry::fillDetIds() const {
0024 std::lock_guard<std::mutex> guard(s_fillLock);
0025 if (m_filledDetIds) {
0026
0027 return;
0028 }
0029 const std::vector<DetId>& baseIds(CaloSubdetectorGeometry::getValidDetIds());
0030 for (unsigned int i(0); i != baseIds.size(); ++i) {
0031 const DetId id(baseIds[i]);
0032 if (id.subdetId() == HcalBarrel) {
0033 m_hbIds.emplace_back(id);
0034 } else {
0035 if (id.subdetId() == HcalEndcap) {
0036 m_heIds.emplace_back(id);
0037 } else {
0038 if (id.subdetId() == HcalOuter) {
0039 m_hoIds.emplace_back(id);
0040 } else {
0041 if (id.subdetId() == HcalForward) {
0042 m_hfIds.emplace_back(id);
0043 }
0044 }
0045 }
0046 }
0047 }
0048 std::sort(m_hbIds.begin(), m_hbIds.end());
0049 std::sort(m_heIds.begin(), m_heIds.end());
0050 std::sort(m_hoIds.begin(), m_hoIds.end());
0051 std::sort(m_hfIds.begin(), m_hfIds.end());
0052
0053 m_emptyIds.resize(0);
0054 m_filledDetIds = true;
0055 }
0056
0057 std::vector<DetId> const& HcalDDDGeometry::getValidDetIds(DetId::Detector det, int subdet) const {
0058 if (0 != subdet && not m_filledDetIds)
0059 fillDetIds();
0060 return (0 == subdet
0061 ? CaloSubdetectorGeometry::getValidDetIds()
0062 : (HcalBarrel == subdet
0063 ? m_hbIds
0064 : (HcalEndcap == subdet
0065 ? m_heIds
0066 : (HcalOuter == subdet ? m_hoIds : (HcalForward == subdet ? m_hfIds : m_emptyIds)))));
0067 }
0068
0069 DetId HcalDDDGeometry::getClosestCell(const GlobalPoint& r) const {
0070 constexpr double twopi = M_PI + M_PI;
0071
0072
0073
0074 double abseta = fabs(r.eta());
0075 double phi = r.phi();
0076 if (phi < 0)
0077 phi += twopi;
0078 double radius = r.mag();
0079 double z = fabs(r.z());
0080
0081 #ifdef EDM_ML_DEBUG
0082 constexpr double deg = M_PI / 180.;
0083 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometry::getClosestCell for eta " << r.eta() << " phi " << phi / deg << " z "
0084 << r.z() << " radius " << radius;
0085 #endif
0086 HcalDetId bestId;
0087 if (abseta <= etaMax_) {
0088 for (const auto& hcalCell : hcalCells_) {
0089 if (abseta >= hcalCell.etaMin() && abseta <= hcalCell.etaMax()) {
0090 HcalSubdetector bc = hcalCell.detType();
0091 int etaring = hcalCell.etaBin();
0092 int phibin = 0;
0093 if (hcalCell.unitPhi() == 4) {
0094
0095
0096
0097
0098 phibin =
0099 static_cast<int>((phi + hcalCell.phiOffset() + 0.5 * hcalCell.phiBinWidth()) / hcalCell.phiBinWidth());
0100 if (phibin == 0)
0101 phibin = hcalCell.nPhiBins();
0102 phibin = phibin * 4 - 1;
0103 } else {
0104 phibin = static_cast<int>((phi + hcalCell.phiOffset()) / hcalCell.phiBinWidth()) + 1;
0105
0106 phibin = (phibin - 1) * (hcalCell.unitPhi()) + 1;
0107 }
0108
0109 int dbin = 1;
0110 int etabin = (r.z() > 0) ? etaring : -etaring;
0111 if (bc == HcalForward) {
0112 bestId = HcalDetId(bc, etabin, phibin, dbin);
0113 break;
0114 } else {
0115 double rz = z;
0116 if (hcalCell.depthType())
0117 rz = radius;
0118 if (rz < hcalCell.depthMax()) {
0119 dbin = hcalCell.depthSegment();
0120 bestId = HcalDetId(bc, etabin, phibin, dbin);
0121 break;
0122 }
0123 }
0124 }
0125 }
0126 }
0127
0128 #ifdef EDM_ML_DEBUG
0129 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometry::getClosestCell " << bestId;
0130 #endif
0131 return bestId;
0132 }
0133
0134 int HcalDDDGeometry::insertCell(std::vector<HcalCellType> const& cells) {
0135 hcalCells_.insert(hcalCells_.end(), cells.begin(), cells.end());
0136 int num = static_cast<int>(hcalCells_.size());
0137 for (const auto& cell : cells) {
0138 if (cell.etaMax() > etaMax_)
0139 etaMax_ = cell.etaMax();
0140 }
0141 #ifdef EDM_ML_DEBUG
0142 edm::LogVerbatim("HCalGeom") << "HcalDDDGeometry::insertCell " << cells.size() << " cells inserted == Total " << num
0143 << " EtaMax = " << etaMax_;
0144 #endif
0145 return num;
0146 }
0147
0148 void HcalDDDGeometry::newCellImpl(
0149 const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0150 assert(detId.det() == DetId::Hcal);
0151
0152 const unsigned int din(topo_.detId2denseId(detId));
0153
0154 HcalDetId hId(detId);
0155
0156 if (hId.subdet() == HcalBarrel) {
0157 m_hbCellVec[din] = IdealObliquePrism(f1, cornersMgr(), parm);
0158 } else {
0159 if (hId.subdet() == HcalEndcap) {
0160 const unsigned int index(din - m_hbCellVec.size());
0161 m_heCellVec[index] = IdealObliquePrism(f1, cornersMgr(), parm);
0162 } else {
0163 if (hId.subdet() == HcalOuter) {
0164 const unsigned int index(din - m_hbCellVec.size() - m_heCellVec.size());
0165 m_hoCellVec[index] = IdealObliquePrism(f1, cornersMgr(), parm);
0166 } else {
0167 const unsigned int index(din - m_hbCellVec.size() - m_heCellVec.size() - m_hoCellVec.size());
0168 m_hfCellVec[index] =
0169 IdealZPrism(f1, cornersMgr(), parm, hId.depth() == 1 ? IdealZPrism::EM : IdealZPrism::HADR);
0170 }
0171 }
0172 }
0173 }
0174
0175 void HcalDDDGeometry::newCell(
0176 const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0177 newCellImpl(f1, f2, f3, parm, detId);
0178 addValidID(detId);
0179 }
0180
0181 void HcalDDDGeometry::newCellFast(
0182 const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0183 newCellImpl(f1, f2, f3, parm, detId);
0184 m_validIds.emplace_back(detId);
0185 }
0186
0187 CaloCellGeometryPtr HcalDDDGeometry::getGeometryRawPtr(uint32_t din) const {
0188
0189 const CaloCellGeometry* cell(nullptr);
0190 if (m_hbCellVec.size() > din) {
0191 cell = (&m_hbCellVec[din]);
0192 } else if (m_hbCellVec.size() + m_heCellVec.size() > din) {
0193 const unsigned int ind(din - m_hbCellVec.size());
0194 cell = (&m_heCellVec[ind]);
0195 } else if (m_hbCellVec.size() + m_heCellVec.size() + m_hoCellVec.size() > din) {
0196 const unsigned int ind(din - m_hbCellVec.size() - m_heCellVec.size());
0197 cell = (&m_hoCellVec[ind]);
0198 } else if (m_hbCellVec.size() + m_heCellVec.size() + m_hoCellVec.size() + m_hfCellVec.size() > din) {
0199 const unsigned int ind(din - m_hbCellVec.size() - m_heCellVec.size() - m_hoCellVec.size());
0200 cell = (&m_hfCellVec[ind]);
0201 }
0202
0203 return CaloCellGeometryPtr((nullptr == cell || nullptr == cell->param()) ? nullptr : cell);
0204 }
0205
0206 void HcalDDDGeometry::increaseReserve(unsigned int extra) { m_validIds.reserve(m_validIds.size() + extra); }
0207
0208 void HcalDDDGeometry::sortValidIds() { std::sort(m_validIds.begin(), m_validIds.end()); }