CaloGeometry/src/CaloSubdetectorGeometry.cc

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 /*det*/, int /*subdet*/) const {
0033   return m_validIds;
0034 }
0035
0036 CaloSubdetectorGeometry::CellMayOwnPtr CaloSubdetectorGeometry::getGeometry(const DetId& id) const {
0037   return CellMayOwnPtr(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     auto 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       auto 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   // stupid implementation not to be really used...
0098   DetIdSet ids = getCells(r, dR);
0099   CellSet cells;
0100   cells.reserve(ids.size());
0101   static const auto do_not_delete = [](const void*) {};
0102   for (auto id : ids)
0103     cells.emplace_back(std::shared_ptr<CaloCellGeometry const>(getGeometry(id).get(), do_not_delete));
0104   return cells;
0105 }
0106
0107 void CaloSubdetectorGeometry::allocateCorners(CaloCellGeometry::CornersVec::size_type n) {
0108   assert(nullptr == m_cmgr);
0109   m_cmgr = new CaloCellGeometry::CornersMgr(n * (CaloCellGeometry::k_cornerSize), CaloCellGeometry::k_cornerSize);
0110
0111   m_validIds.reserve(n);
0112 }
0113
0114 void CaloSubdetectorGeometry::allocatePar(ParVec::size_type n, unsigned int m) {
0115   assert(nullptr == m_parMgr);
0116   m_parMgr = new ParMgr(n * m, m);
0117 }
0118
0119 void CaloSubdetectorGeometry::getSummary(CaloSubdetectorGeometry::TrVec& tVec,
0120                                          CaloSubdetectorGeometry::IVec& iVec,
0121                                          CaloSubdetectorGeometry::DimVec& dVec,
0122                                          CaloSubdetectorGeometry::IVec& /*dins*/) const {
0123   tVec.reserve(m_validIds.size() * numberOfTransformParms());
0124   iVec.reserve(numberOfShapes() == 1 ? 1 : m_validIds.size());
0125   dVec.reserve(numberOfShapes() * numberOfParametersPerShape());
0126
0127   for (const auto& pv : parVecVec()) {
0128     for (float iv : pv) {
0129       dVec.emplace_back(iv);
0130     }
0131   }
0132
0133   for (uint32_t i(0); i != m_validIds.size(); ++i) {
0134     Tr3D tr;
0135     auto ptr{cellGeomPtr(i)};
0136     assert(nullptr != ptr);
0137     ptr->getTransform(tr, (Pt3DVec*)nullptr);
0138
0139     if (Tr3D() == tr) {  // for preshower there is no rotation
0140       const GlobalPoint& gp(ptr->getPosition());
0141       tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
0142     }
0143
0144     const CLHEP::Hep3Vector tt(tr.getTranslation());
0145     tVec.emplace_back(tt.x());
0146     tVec.emplace_back(tt.y());
0147     tVec.emplace_back(tt.z());
0148     if (6 == numberOfTransformParms()) {
0149       const CLHEP::HepRotation rr(tr.getRotation());
0150       const ROOT::Math::Transform3D rtr(
0151           rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
0152       ROOT::Math::EulerAngles ea;
0153       rtr.GetRotation(ea);
0154       tVec.emplace_back(ea.Phi());
0155       tVec.emplace_back(ea.Theta());
0156       tVec.emplace_back(ea.Psi());
0157     }
0158
0159     const CCGFloat* par(ptr->param());
0160
0161     unsigned int ishape(9999);
0162     for (unsigned int ivv(0); ivv != parVecVec().size(); ++ivv) {
0163       bool ok(true);
0164       const CCGFloat* pv(&(*parVecVec()[ivv].begin()));
0165       for (unsigned int k(0); k != numberOfParametersPerShape(); ++k) {
0166         ok = ok && (fabs(par[k] - pv[k]) < 1.e-6);
0167       }
0168       if (ok) {
0169         ishape = ivv;
0170         break;
0171       }
0172     }
0173     assert(9999 != ishape);
0174
0175     const unsigned int nn((numberOfShapes() == 1) ? (unsigned int)1 : m_validIds.size());
0176     if (iVec.size() < nn)
0177       iVec.emplace_back(ishape);
0178   }
0179 }
0180
0181 CCGFloat CaloSubdetectorGeometry::deltaPhi(const DetId& detId) const {
0182   const CaloGenericDetId cgId(detId);
0183
0184   if (!m_deltaPhi.load(std::memory_order_acquire)) {
0185     const uint32_t kSize(sizeForDenseIndex(detId));
0186     auto ptr = new std::vector<CCGFloat>(kSize);
0187     for (uint32_t i(0); i != kSize; ++i) {
0188       auto cellPtr{cellGeomPtr(i)};
0189       if (nullptr != cellPtr) {
0190         CCGFloat dPhi1(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[1].x()) / 2.,
0191                                         (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[1].y()) / 2.,
0192                                         (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[1].z()) / 2.)
0193                                 .phi() -
0194                             GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[3].x()) / 2.,
0195                                         (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[3].y()) / 2.,
0196                                         (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[3].z()) / 2.)
0197                                 .phi()));
0198         CCGFloat dPhi2(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[3].x()) / 2.,
0199                                         (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[3].y()) / 2.,
0200                                         (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[3].z()) / 2.)
0201                                 .phi() -
0202                             GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[1].x()) / 2.,
0203                                         (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[1].y()) / 2.,
0204                                         (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[1].z()) / 2.)
0205                                 .phi()));
0206         if (M_PI < dPhi1)
0207           dPhi1 = fabs(dPhi1 - 2. * M_PI);
0208         if (M_PI < dPhi2)
0209           dPhi2 = fabs(dPhi2 - 2. * M_PI);
0210         (*ptr)[i] = dPhi1 > dPhi2 ? dPhi1 : dPhi2;
0211       }
0212     }
0213     std::vector<CCGFloat>* expect = nullptr;
0214     bool exchanged = m_deltaPhi.compare_exchange_strong(expect, ptr, std::memory_order_acq_rel);
0215     if (!exchanged)
0216       delete ptr;
0217   }
0218   return (*m_deltaPhi.load(std::memory_order_acquire))[indexFor(detId)];
0219 }
0220
0221 CCGFloat CaloSubdetectorGeometry::deltaEta(const DetId& detId) const {
0222   if (!m_deltaEta.load(std::memory_order_acquire)) {
0223     const uint32_t kSize(sizeForDenseIndex(detId));
0224     auto ptr = new std::vector<CCGFloat>(kSize);
0225     for (uint32_t i(0); i != kSize; ++i) {
0226       auto cellPtr{cellGeomPtr(i)};
0227       if (nullptr != cellPtr) {
0228         const CCGFloat dEta1(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[1].x()) / 2.,
0229                                               (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[1].y()) / 2.,
0230                                               (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[1].z()) / 2.)
0231                                       .eta() -
0232                                   GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[3].x()) / 2.,
0233                                               (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[3].y()) / 2.,
0234                                               (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[3].z()) / 2.)
0235                                       .eta()));
0236         const CCGFloat dEta2(fabs(GlobalPoint((cellPtr->getCorners()[0].x() + cellPtr->getCorners()[3].x()) / 2.,
0237                                               (cellPtr->getCorners()[0].y() + cellPtr->getCorners()[3].y()) / 2.,
0238                                               (cellPtr->getCorners()[0].z() + cellPtr->getCorners()[3].z()) / 2.)
0239                                       .eta() -
0240                                   GlobalPoint((cellPtr->getCorners()[2].x() + cellPtr->getCorners()[1].x()) / 2.,
0241                                               (cellPtr->getCorners()[2].y() + cellPtr->getCorners()[1].y()) / 2.,
0242                                               (cellPtr->getCorners()[2].z() + cellPtr->getCorners()[1].z()) / 2.)
0243                                       .eta()));
0244         (*ptr)[i] = dEta1 > dEta2 ? dEta1 : dEta2;
0245       }
0246     }
0247     std::vector<CCGFloat>* expect = nullptr;
0248     bool exchanged = m_deltaEta.compare_exchange_strong(expect, ptr, std::memory_order_acq_rel);
0249     if (!exchanged)
0250       delete ptr;
0251   }
0252   return (*m_deltaEta.load(std::memory_order_acquire))[indexFor(detId)];
0253 }
0254
0255 unsigned int CaloSubdetectorGeometry::indexFor(const DetId& id) const { return CaloGenericDetId(id).denseIndex(); }
0256
0257 unsigned int CaloSubdetectorGeometry::sizeForDenseIndex(const DetId& id) const {
0258   return CaloGenericDetId(id).sizeForDenseIndexing();
0259 }
0260
0261 CaloSubdetectorGeometry::CellPtr CaloSubdetectorGeometry::cellGeomPtr(uint32_t index) const {
0262   // Default version
0263   return getGeometryRawPtr(index);
0264 }