HcalDDDGeometry.cc

Go to the documentation of this file.

#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/HcalTowerAlgo/interface/HcalDDDGeometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <algorithm>
#include <mutex>

//#define EDM_ML_DEBUG

static std::mutex s_fillLock;

HcalDDDGeometry::HcalDDDGeometry(const HcalTopology& topo)
    : topo_(topo),
      etaMax_(0),
      m_hbCellVec(topo.getHBSize()),
      m_heCellVec(topo.getHESize()),
      m_hoCellVec(topo.getHOSize()),
      m_hfCellVec(topo.getHFSize()),
      m_filledDetIds(false) {}

HcalDDDGeometry::~HcalDDDGeometry() {}

void HcalDDDGeometry::fillDetIds() const {
  std::lock_guard<std::mutex> guard(s_fillLock);
  if (m_filledDetIds) {
    //another thread already did the work
    return;
  }
  const std::vector<DetId>& baseIds(CaloSubdetectorGeometry::getValidDetIds());
  for (unsigned int i(0); i != baseIds.size(); ++i) {
    const DetId id(baseIds[i]);
    if (id.subdetId() == HcalBarrel) {
      m_hbIds.emplace_back(id);
    } else {
      if (id.subdetId() == HcalEndcap) {
        m_heIds.emplace_back(id);
      } else {
        if (id.subdetId() == HcalOuter) {
          m_hoIds.emplace_back(id);
        } else {
          if (id.subdetId() == HcalForward) {
            m_hfIds.emplace_back(id);
          }
        }
      }
    }
  }
  std::sort(m_hbIds.begin(), m_hbIds.end());
  std::sort(m_heIds.begin(), m_heIds.end());
  std::sort(m_hoIds.begin(), m_hoIds.end());
  std::sort(m_hfIds.begin(), m_hfIds.end());

  m_emptyIds.resize(0);
  m_filledDetIds = true;
}

std::vector<DetId> const& HcalDDDGeometry::getValidDetIds(DetId::Detector det, int subdet) const {
  if (0 != subdet && not m_filledDetIds)
    fillDetIds();
  return (0 == subdet
              ? CaloSubdetectorGeometry::getValidDetIds()
              : (HcalBarrel == subdet
                     ? m_hbIds
                     : (HcalEndcap == subdet
                            ? m_heIds
                            : (HcalOuter == subdet ? m_hoIds : (HcalForward == subdet ? m_hfIds : m_emptyIds)))));
}

DetId HcalDDDGeometry::getClosestCell(const GlobalPoint& r) const {
  constexpr double twopi = M_PI + M_PI;

  // Now find the closest eta_bin, eta value of a bin i is average
  // of eta[i] and eta[i-1]
  double abseta = fabs(r.eta());
  double phi = r.phi();
  if (phi < 0)
    phi += twopi;
  double radius = r.mag();
  double z = fabs(r.z());

#ifdef EDM_ML_DEBUG
  constexpr double deg = M_PI / 180.;
  edm::LogVerbatim("HCalGeom") << "HcalDDDGeometry::getClosestCell for eta " << r.eta() << " phi " << phi / deg << " z "
                               << r.z() << " radius " << radius;
#endif
  HcalDetId bestId;
  if (abseta <= etaMax_) {
    for (const auto& hcalCell : hcalCells_) {
      if (abseta >= hcalCell.etaMin() && abseta <= hcalCell.etaMax()) {
        HcalSubdetector bc = hcalCell.detType();
        int etaring = hcalCell.etaBin();
        int phibin = 0;
        if (hcalCell.unitPhi() == 4) {
          // rings 40 and 41 are offset wrt the other phi numbering
          //  1        1         1         2
          //  ------------------------------
          //  72       36        36        1
          phibin =
              static_cast<int>((phi + hcalCell.phiOffset() + 0.5 * hcalCell.phiBinWidth()) / hcalCell.phiBinWidth());
          if (phibin == 0)
            phibin = hcalCell.nPhiBins();
          phibin = phibin * 4 - 1;
        } else {
          phibin = static_cast<int>((phi + hcalCell.phiOffset()) / hcalCell.phiBinWidth()) + 1;
          // convert to the convention of numbering 1,3,5, in 36 phi bins
          phibin = (phibin - 1) * (hcalCell.unitPhi()) + 1;
        }

        int dbin = 1;
        int etabin = (r.z() > 0) ? etaring : -etaring;
        if (bc == HcalForward) {
          bestId = HcalDetId(bc, etabin, phibin, dbin);
          break;
        } else {
          double rz = z;
          if (hcalCell.depthType())
            rz = radius;
          if (rz < hcalCell.depthMax()) {
            dbin = hcalCell.depthSegment();
            bestId = HcalDetId(bc, etabin, phibin, dbin);
            break;
          }
        }
      }
    }
  }

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalDDDGeometry::getClosestCell " << bestId;
#endif
  return bestId;
}

int HcalDDDGeometry::insertCell(std::vector<HcalCellType> const& cells) {
  hcalCells_.insert(hcalCells_.end(), cells.begin(), cells.end());
  int num = static_cast<int>(hcalCells_.size());
  for (const auto& cell : cells) {
    if (cell.etaMax() > etaMax_)
      etaMax_ = cell.etaMax();
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalDDDGeometry::insertCell " << cells.size() << " cells inserted == Total " << num
                               << " EtaMax = " << etaMax_;
#endif
  return num;
}

void HcalDDDGeometry::newCellImpl(
    const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
  assert(detId.det() == DetId::Hcal);

  const unsigned int din(topo_.detId2denseId(detId));

  HcalDetId hId(detId);

  if (hId.subdet() == HcalBarrel) {
    m_hbCellVec[din] = IdealObliquePrism(f1, cornersMgr(), parm);
  } else {
    if (hId.subdet() == HcalEndcap) {
      const unsigned int index(din - m_hbCellVec.size());
      m_heCellVec[index] = IdealObliquePrism(f1, cornersMgr(), parm);
    } else {
      if (hId.subdet() == HcalOuter) {
        const unsigned int index(din - m_hbCellVec.size() - m_heCellVec.size());
        m_hoCellVec[index] = IdealObliquePrism(f1, cornersMgr(), parm);
      } else {  // assuming HcalForward here!
        const unsigned int index(din - m_hbCellVec.size() - m_heCellVec.size() - m_hoCellVec.size());
        m_hfCellVec[index] =
            IdealZPrism(f1, cornersMgr(), parm, hId.depth() == 1 ? IdealZPrism::EM : IdealZPrism::HADR);
      }
    }
  }
}

void HcalDDDGeometry::newCell(
    const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
  newCellImpl(f1, f2, f3, parm, detId);
  addValidID(detId);
}

void HcalDDDGeometry::newCellFast(
    const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
  newCellImpl(f1, f2, f3, parm, detId);
  m_validIds.emplace_back(detId);
}

const CaloCellGeometry* HcalDDDGeometry::getGeometryRawPtr(uint32_t din) const {
  // Modify the RawPtr class
  const CaloCellGeometry* cell(nullptr);
  if (m_hbCellVec.size() > din) {
    cell = (&m_hbCellVec[din]);
  } else if (m_hbCellVec.size() + m_heCellVec.size() > din) {
    const unsigned int ind(din - m_hbCellVec.size());
    cell = (&m_heCellVec[ind]);
  } else if (m_hbCellVec.size() + m_heCellVec.size() + m_hoCellVec.size() > din) {
    const unsigned int ind(din - m_hbCellVec.size() - m_heCellVec.size());
    cell = (&m_hoCellVec[ind]);
  } else if (m_hbCellVec.size() + m_heCellVec.size() + m_hoCellVec.size() + m_hfCellVec.size() > din) {
    const unsigned int ind(din - m_hbCellVec.size() - m_heCellVec.size() - m_hoCellVec.size());
    cell = (&m_hfCellVec[ind]);
  }

  return ((nullptr == cell || nullptr == cell->param()) ? nullptr : cell);
}

void HcalDDDGeometry::increaseReserve(unsigned int extra) { m_validIds.reserve(m_validIds.size() + extra); }

void HcalDDDGeometry::sortValidIds() { std::sort(m_validIds.begin(), m_validIds.end()); }