File indexing completed on 2024-04-06 12:14:16
0001 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0002 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0003
0004 #include <Math/Transform3D.h>
0005 #include <Math/EulerAngles.h>
0006
0007 #include <algorithm>
0008
0009 typedef CaloCellGeometry::Pt3D Pt3D;
0010 typedef CaloCellGeometry::Pt3DVec Pt3DVec;
0011 typedef CaloCellGeometry::Tr3D Tr3D;
0012
0013 typedef CaloSubdetectorGeometry::CCGFloat CCGFloat;
0014
0015 CaloSubdetectorGeometry::CaloSubdetectorGeometry()
0016 : m_parMgr(nullptr), m_cmgr(nullptr), m_deltaPhi(nullptr), m_deltaEta(nullptr) {}
0017
0018 CaloSubdetectorGeometry::~CaloSubdetectorGeometry() {
0019 delete m_cmgr;
0020 delete m_parMgr;
0021 if (m_deltaPhi)
0022 delete m_deltaPhi.load();
0023 if (m_deltaEta)
0024 delete m_deltaEta.load();
0025 }
0026
0027 void CaloSubdetectorGeometry::addValidID(const DetId& id) {
0028 auto pos = std::lower_bound(m_validIds.begin(), m_validIds.end(), id);
0029 m_validIds.insert(pos, id);
0030 }
0031
0032 const std::vector<DetId>& CaloSubdetectorGeometry::getValidDetIds(DetId::Detector , int ) const {
0033 return m_validIds;
0034 }
0035
0036 std::shared_ptr<const CaloCellGeometry> CaloSubdetectorGeometry::getGeometry(const DetId& id) const {
0037 return cellGeomPtr(CaloGenericDetId(id).denseIndex());
0038 }
0039
0040 bool CaloSubdetectorGeometry::present(const DetId& id) const {
0041 return std::find(m_validIds.begin(), m_validIds.end(), id) != m_validIds.end();
0042 }
0043
0044 DetId CaloSubdetectorGeometry::getClosestCell(const GlobalPoint& r) const {
0045 const CCGFloat eta(r.eta());
0046 const CCGFloat phi(r.phi());
0047 uint32_t index(~0);
0048 CCGFloat closest(1e9);
0049
0050 for (uint32_t i(0); i != m_validIds.size(); ++i) {
0051 std::shared_ptr<const CaloCellGeometry> cell(getGeometry(m_validIds[i]));
0052 if (nullptr != cell) {
0053 const GlobalPoint& p(cell->getPosition());
0054 const CCGFloat eta0(p.eta());
0055 const CCGFloat phi0(p.phi());
0056 const CCGFloat dR2(reco::deltaR2(eta0, phi0, eta, phi));
0057 if (dR2 < closest) {
0058 closest = dR2;
0059 index = i;
0060 }
0061 }
0062 }
0063 return (closest > 0.9e9 || (uint32_t)(~0) == index ? DetId(0) : m_validIds[index]);
0064 }
0065
0066 CaloSubdetectorGeometry::DetIdSet CaloSubdetectorGeometry::getCells(const GlobalPoint& r, double dR) const {
0067 const double dR2(dR * dR);
0068 const double eta(r.eta());
0069 const double phi(r.phi());
0070
0071 DetIdSet dss;
0072
0073 if (0.000001 < dR) {
0074 for (uint32_t i(0); i != m_validIds.size(); ++i) {
0075 std::shared_ptr<const CaloCellGeometry> cell(getGeometry(m_validIds[i]));
0076 if (nullptr != cell) {
0077 const GlobalPoint& p(cell->getPosition());
0078 const CCGFloat eta0(p.eta());
0079 if (fabs(eta - eta0) < dR) {
0080 const CCGFloat phi0(p.phi());
0081 CCGFloat delp(fabs(phi - phi0));
0082 if (delp > M_PI)
0083 delp = 2 * M_PI - delp;
0084 if (delp < dR) {
0085 const CCGFloat dist2(reco::deltaR2(eta0, phi0, eta, phi));
0086 if (dist2 < dR2)
0087 dss.insert(m_validIds[i]);
0088 }
0089 }
0090 }
0091 }
0092 }
0093 return dss;
0094 }
0095
0096 CaloSubdetectorGeometry::CellSet CaloSubdetectorGeometry::getCellSet(const GlobalPoint& r, double dR) const {
0097
0098 DetIdSet ids = getCells(r, dR);
0099 CellSet cells;
0100 cells.reserve(ids.size());
0101 for (auto id : ids)
0102 cells.emplace_back(getGeometry(id));
0103 return cells;
0104 }
0105
0106 void CaloSubdetectorGeometry::allocateCorners(CaloCellGeometry::CornersVec::size_type n) {
0107 assert(nullptr == m_cmgr);
0108 m_cmgr = new CaloCellGeometry::CornersMgr(n * (CaloCellGeometry::k_cornerSize), CaloCellGeometry::k_cornerSize);
0109
0110 m_validIds.reserve(n);
0111 }
0112
0113 void CaloSubdetectorGeometry::allocatePar(ParVec::size_type n, unsigned int m) {
0114 assert(nullptr == m_parMgr);
0115 m_parMgr = new ParMgr(n * m, m);
0116 }
0117
0118 void CaloSubdetectorGeometry::getSummary(CaloSubdetectorGeometry::TrVec& tVec,
0119 CaloSubdetectorGeometry::IVec& iVec,
0120 CaloSubdetectorGeometry::DimVec& dVec,
0121 CaloSubdetectorGeometry::IVec& ) const {
0122 tVec.reserve(m_validIds.size() * numberOfTransformParms());
0123 iVec.reserve(numberOfShapes() == 1 ? 1 : m_validIds.size());
0124 dVec.reserve(numberOfShapes() * numberOfParametersPerShape());
0125
0126 for (const auto& pv : parVecVec()) {
0127 for (float iv : pv) {
0128 dVec.emplace_back(iv);
0129 }
0130 }
0131
0132 for (uint32_t i(0); i != m_validIds.size(); ++i) {
0133 Tr3D tr;
0134 std::shared_ptr<const CaloCellGeometry> ptr(cellGeomPtr(i));
0135 assert(nullptr != ptr);
0136 ptr->getTransform(tr, (Pt3DVec*)nullptr);
0137
0138 if (Tr3D() == tr) {
0139 const GlobalPoint& gp(ptr->getPosition());
0140 tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
0141 }
0142
0143 const CLHEP::Hep3Vector tt(tr.getTranslation());
0144 tVec.emplace_back(tt.x());
0145 tVec.emplace_back(tt.y());
0146 tVec.emplace_back(tt.z());
0147 if (6 == numberOfTransformParms()) {
0148 const CLHEP::HepRotation rr(tr.getRotation());
0149 const ROOT::Math::Transform3D rtr(
0150 rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
0151 ROOT::Math::EulerAngles ea;
0152 rtr.GetRotation(ea);
0153 tVec.emplace_back(ea.Phi());
0154 tVec.emplace_back(ea.Theta());
0155 tVec.emplace_back(ea.Psi());
0156 }
0157
0158 const CCGFloat* par(ptr->param());
0159
0160 unsigned int ishape(9999);
0161 for (unsigned int ivv(0); ivv != parVecVec().size(); ++ivv) {
0162 bool ok(true);
0163 const CCGFloat* pv(&(*parVecVec()[ivv].begin()));
0164 for (unsigned int k(0); k != numberOfParametersPerShape(); ++k) {
0165 ok = ok && (fabs(par[k] - pv[k]) < 1.e-6);
0166 }
0167 if (ok) {
0168 ishape = ivv;
0169 break;
0170 }
0171 }
0172 assert(9999 != ishape);
0173
0174 const unsigned int nn((numberOfShapes() == 1) ? (unsigned int)1 : m_validIds.size());
0175 if (iVec.size() < nn)
0176 iVec.emplace_back(ishape);
0177 }
0178 }
0179
0180 CCGFloat CaloSubdetectorGeometry::deltaPhi(const DetId& detId) const {
0181 const CaloGenericDetId cgId(detId);
0182
0183 if (!m_deltaPhi.load(std::memory_order_acquire)) {
0184 const uint32_t kSize(sizeForDenseIndex(detId));
0185 auto ptr = new std::vector<CCGFloat>(kSize);
0186 for (uint32_t i(0); i != kSize; ++i) {
0187 std::shared_ptr<const CaloCellGeometry> cellPtr(cellGeomPtr(i));
0188 if (nullptr != cellPtr) {
0189 CCGFloat dPhi1(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[1].x()) / 2.,
0190 (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[1].y()) / 2.,
0191 (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[1].z()) / 2.)
0192 .phi() -
0193 GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[3].x()) / 2.,
0194 (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[3].y()) / 2.,
0195 (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[3].z()) / 2.)
0196 .phi()));
0197 CCGFloat dPhi2(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[3].x()) / 2.,
0198 (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[3].y()) / 2.,
0199 (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[3].z()) / 2.)
0200 .phi() -
0201 GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[1].x()) / 2.,
0202 (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[1].y()) / 2.,
0203 (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[1].z()) / 2.)
0204 .phi()));
0205 if (M_PI < dPhi1)
0206 dPhi1 = fabs(dPhi1 - 2. * M_PI);
0207 if (M_PI < dPhi2)
0208 dPhi2 = fabs(dPhi2 - 2. * M_PI);
0209 (*ptr)[i] = dPhi1 > dPhi2 ? dPhi1 : dPhi2;
0210 }
0211 }
0212 std::vector<CCGFloat>* expect = nullptr;
0213 bool exchanged = m_deltaPhi.compare_exchange_strong(expect, ptr, std::memory_order_acq_rel);
0214 if (!exchanged)
0215 delete ptr;
0216 }
0217 return (*m_deltaPhi.load(std::memory_order_acquire))[indexFor(detId)];
0218 }
0219
0220 CCGFloat CaloSubdetectorGeometry::deltaEta(const DetId& detId) const {
0221 if (!m_deltaEta.load(std::memory_order_acquire)) {
0222 const uint32_t kSize(sizeForDenseIndex(detId));
0223 auto ptr = new std::vector<CCGFloat>(kSize);
0224 for (uint32_t i(0); i != kSize; ++i) {
0225 std::shared_ptr<const CaloCellGeometry> cellPtr(cellGeomPtr(i));
0226 if (nullptr != cellPtr) {
0227 const CCGFloat dEta1(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[1].x()) / 2.,
0228 (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[1].y()) / 2.,
0229 (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[1].z()) / 2.)
0230 .eta() -
0231 GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[3].x()) / 2.,
0232 (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[3].y()) / 2.,
0233 (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[3].z()) / 2.)
0234 .eta()));
0235 const CCGFloat dEta2(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[3].x()) / 2.,
0236 (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[3].y()) / 2.,
0237 (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[3].z()) / 2.)
0238 .eta() -
0239 GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[1].x()) / 2.,
0240 (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[1].y()) / 2.,
0241 (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[1].z()) / 2.)
0242 .eta()));
0243 (*ptr)[i] = dEta1 > dEta2 ? dEta1 : dEta2;
0244 }
0245 }
0246 std::vector<CCGFloat>* expect = nullptr;
0247 bool exchanged = m_deltaEta.compare_exchange_strong(expect, ptr, std::memory_order_acq_rel);
0248 if (!exchanged)
0249 delete ptr;
0250 }
0251 return (*m_deltaEta.load(std::memory_order_acquire))[indexFor(detId)];
0252 }
0253
0254 unsigned int CaloSubdetectorGeometry::indexFor(const DetId& id) const { return CaloGenericDetId(id).denseIndex(); }
0255
0256 unsigned int CaloSubdetectorGeometry::sizeForDenseIndex(const DetId& id) const {
0257 return CaloGenericDetId(id).sizeForDenseIndexing();
0258 }
0259
0260 std::shared_ptr<const CaloCellGeometry> CaloSubdetectorGeometry::cellGeomPtr(uint32_t index) const {
0261
0262 auto ptr = getGeometryRawPtr(index);
0263 static const auto do_not_delete = [](const void*) {};
0264 return ptr == nullptr ? nullptr : std::shared_ptr<const CaloCellGeometry>(ptr, do_not_delete);
0265 }