HcalTowerAlgo/src/HcalDDDGeometry.cc

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 //#define EDM_ML_DEBUG
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     //another thread already did the work
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   // Now find the closest eta_bin, eta value of a bin i is average
0073   // of eta[i] and eta[i-1]
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           // rings 40 and 41 are offset wrt the other phi numbering
0095           //  1        1         1         2
0096           //  ------------------------------
0097           //  72       36        36        1
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           // convert to the convention of numbering 1,3,5, in 36 phi bins
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 {  // assuming HcalForward here!
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   // Modify the RawPtr class
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()); }