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 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   // stupid implementation not to be really used...
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& /*dins*/) 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) {  // for preshower there is no rotation
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   // Default version
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 }