HcalDDDGeometryLoader.cc

Go to the documentation of this file.

#include "Geometry/HcalTowerAlgo/interface/HcalDDDGeometryLoader.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "Geometry/HcalCommonData/interface/HcalParametersFromDD.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
#include "Geometry/CaloGeometry/interface/IdealZPrism.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <string>

typedef CaloCellGeometry::CCGFloat CCGFloat;

//#define EDM_ML_DEBUG

HcalDDDGeometryLoader::HcalDDDGeometryLoader(const HcalDDDRecConstants* hcons) : hcalConstants_(hcons) {
  isBH_ = hcalConstants_->isBH();
}

HcalDDDGeometryLoader::~HcalDDDGeometryLoader() {}

HcalDDDGeometryLoader::ReturnType HcalDDDGeometryLoader::load(const HcalTopology& topo,
                                                              DetId::Detector det,
                                                              int subdet) {
  HcalSubdetector hsub = static_cast<HcalSubdetector>(subdet);

  HcalDDDGeometry* geom(new HcalDDDGeometry(topo));

  if (geom->cornersMgr() == nullptr) {
    const unsigned int count(hcalConstants_->numberOfCells(HcalBarrel) + hcalConstants_->numberOfCells(HcalEndcap) +
                             2 * hcalConstants_->numberOfCells(HcalForward) + hcalConstants_->numberOfCells(HcalOuter));
    geom->allocateCorners(count);
  }

  //  if( geom->cornersMgr() == 0 )  geom->allocateCorners( 2592 ) ;

  if (geom->parMgr() == nullptr)
    geom->allocatePar(500, 3);

  fill(hsub, geom);
  //fast insertion of valid ids requires sort at end
  geom->sortValidIds();
  return geom;
}

HcalDDDGeometryLoader::ReturnType HcalDDDGeometryLoader::load(const HcalTopology& topo) {
  HcalDDDGeometry* geom(new HcalDDDGeometry(topo));

  if (geom->cornersMgr() == nullptr) {
    const unsigned int count(hcalConstants_->numberOfCells(HcalBarrel) + hcalConstants_->numberOfCells(HcalEndcap) +
                             2 * hcalConstants_->numberOfCells(HcalForward) + hcalConstants_->numberOfCells(HcalOuter));
    geom->allocateCorners(count);
  }
  if (geom->parMgr() == nullptr)
    geom->allocatePar(500, 3);

  fill(HcalBarrel, geom);
  fill(HcalEndcap, geom);
  fill(HcalForward, geom);
  fill(HcalOuter, geom);
  //fast insertion of valid ids requires sort at end
  geom->sortValidIds();
  return geom;
}

void HcalDDDGeometryLoader::fill(HcalSubdetector subdet, HcalDDDGeometry* geom) {
  // start by making the new HcalDetIds
  std::vector<HcalCellType> hcalCells = hcalConstants_->HcalCellTypes(subdet);
  geom->insertCell(hcalCells);
#ifdef EDM_ML_DEBUG
  std::cout << "HcalDDDGeometryLoader::fill gets " << hcalCells.size() << " cells for subdetector " << subdet
            << std::endl;
#endif
  // Make the new HcalDetIds and the cells

  std::vector<HcalDetId> hcalIds;
  for (auto& hcalCell : hcalCells) {
    int etaRing = hcalCell.etaBin();
    int iside = hcalCell.zside();
    int depthBin = hcalCell.depthSegment();
    double dphi = hcalCell.phiBinWidth();
    std::vector<std::pair<int, double> > phis = hcalCell.phis();
#ifdef EDM_ML_DEBUG
    std::cout << "HcalDDDGeometryLoader: Subdet " << subdet << " side " << iside << " eta " << etaRing << " depth "
              << depthBin << " with " << phis.size() << "modules:" << std::endl;
#endif
    geom->increaseReserve(phis.size());
    for (auto& phi : phis) {
#ifdef EDM_ML_DEBUG
      std::cout << "HcalDDDGeometryLoader::fill Cell " << i << " eta " << iside * etaRing << " phi " << phis[k].first
                << "(" << phis[k].second / CLHEP::deg << ", " << dphi / CLHEP::deg << ") depth " << depthBin
                << std::endl;
#endif
      HcalDetId id(subdet, iside * etaRing, phi.first, depthBin);
      hcalIds.emplace_back(id);
      makeCell(id, hcalCell, phi.second, dphi, geom);
    }
  }

  edm::LogInfo("HCalGeom") << "Number of HCAL DetIds made for " << subdet << " is " << hcalIds.size();
}

void HcalDDDGeometryLoader::makeCell(
    const HcalDetId& detId, const HcalCellType& hcalCell, double phi, double dphi, HcalDDDGeometry* geom) const {
  // the two eta boundaries of the cell
  double eta1 = hcalCell.etaMin();
  double eta2 = hcalCell.etaMax();
  HcalSubdetector subdet = detId.subdet();
  double eta = 0.5 * (eta1 + eta2) * detId.zside();
  double deta = (eta2 - eta1);
  double theta = 2.0 * atan(exp(-eta));

  // barrel vs forward
  bool rzType = hcalCell.depthType();
  bool isBarrel = (subdet == HcalBarrel || subdet == HcalOuter);

  double z, r, thickness;
#ifdef EDM_ML_DEBUG
  double r0, r1, r2;
#endif
  if (rzType) {
    r = hcalCell.depthMin();
    if (isBarrel) {
      z = r * sinh(eta);                                  // sinh(eta) == 1/tan(theta)
      thickness = (hcalCell.depthMax() - r) * cosh(eta);  // cosh(eta) == 1/sin(theta)
#ifdef EDM_ML_DEBUG
      r1 = r;
      r2 = hcalCell.depthMax();
      r0 = 0.5 * (r1 + r2);
#endif
    } else {
      z = r * sinh(eta2);
      thickness = 2. * hcalCell.halfSize();
      r = z / sinh(std::abs(eta));
#ifdef EDM_ML_DEBUG
      r0 = z / sinh(std::abs(eta));
      r1 = z / sinh(std::abs(eta) + 0.5 * deta);
      r2 = z / sinh(std::abs(eta) - 0.5 * deta);
#endif
    }
#ifdef EDM_ML_DEBUG
    std::cout << "HcalDDDGeometryLoader::makeCell SubDet " << subdet << " eta = " << eta << " theta = " << theta
              << " r = " << r << " thickness = " << thickness << " r0-r2 (" << r0 << ":" << r1 << ":" << r2 << ")"
              << std::endl;
#endif
  } else {
    z = hcalCell.depthMin();
    thickness = hcalCell.depthMax() - z;
    if (isBH_)
      z += (0.5 * thickness);
    z *= detId.zside();  // get the sign right.
    r = z * tan(theta);
    thickness /= std::abs(cos(theta));
#ifdef EDM_ML_DEBUG
    r0 = z / sinh(std::abs(eta));
    r1 = z / sinh(std::abs(eta) + 0.5 * deta);
    r2 = z / sinh(std::abs(eta) - 0.5 * deta);
    std::cout << "HcalDDDGeometryLoader::makeCell SubDet " << subdet << " eta = " << eta << " theta = " << theta
              << " z = " << z << " r = " << r << " thickness = " << thickness << " r0-r2 (" << r0 << ":" << r1 << ":"
              << r2 << ")" << std::endl;
#endif
  }

  double x = r * cos(phi);
  double y = r * sin(phi);
  GlobalPoint point(x, y, z);

#ifdef EDM_ML_DEBUG
  std::cout << "HcalDDDGeometryLoader::makeCell for " << detId << " Point " << point << " deta = " << deta
            << " dphi = " << dphi << " thickness = " << thickness << " isBarrel = " << isBarrel << " " << rzType
            << std::endl;
#endif

  std::vector<CCGFloat> hp;
  hp.reserve(3);

  if (subdet == HcalForward) {
    hp.emplace_back(deta / 2.);
    hp.emplace_back(dphi / 2.);
    hp.emplace_back(thickness / 2.);
  } else {
    const double sign(isBarrel ? 1 : -1);
    hp.emplace_back(deta / 2.);
    hp.emplace_back(dphi / 2.);
    hp.emplace_back(sign * thickness / 2.);
  }
  geom->newCellFast(point, point, point, CaloCellGeometry::getParmPtr(hp, geom->parMgr(), geom->parVecVec()), detId);
}