CaloTowerHardcodeGeometryLoader.cc

Go to the documentation of this file.

#include "Geometry/HcalTowerAlgo/interface/CaloTowerHardcodeGeometryLoader.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
#include <iterator>
#include <algorithm>
 
//#define DebugLog
 
typedef CaloCellGeometry::CCGFloat CCGFloat;
 
std::unique_ptr<CaloSubdetectorGeometry> CaloTowerHardcodeGeometryLoader::load(const CaloTowerTopology* limits,
                                                                               const HcalTopology* hcaltopo,
                                                                               const HcalDDDRecConstants* hcons) {
  m_limits = limits;
  m_hcaltopo = hcaltopo;
  m_hcons = hcons;
 
  //get eta limits from hcal rec constants
  theHFEtaBounds = m_hcons->getEtaTableHF();
  theHBHEEtaBounds = m_hcons->getEtaTable();
 
#ifdef DebugLog
  std::cout << "CaloTowerHardcodeGeometryLoader: theHBHEEtaBounds = ";
  std::copy(theHBHEEtaBounds.begin(), theHBHEEtaBounds.end(), std::ostream_iterator<double>(std::cout, ","));
  std::cout << std::endl;
 
  std::cout << "CaloTowerHardcodeGeometryLoader: lastHBRing = " << m_limits->lastHBRing()
            << ", lastHERing = " << m_limits->lastHERing() << std::endl;
  std::cout << "CaloTowerHardcodeGeometryLoader: HcalTopology: firstHBRing = " << hcaltopo->firstHBRing()
            << ", lastHBRing = " << hcaltopo->lastHBRing() << ", firstHERing = " << hcaltopo->firstHERing()
            << ", lastHERing = " << hcaltopo->lastHERing() << ", firstHFRing = " << hcaltopo->firstHFRing()
            << ", lastHFRing = " << hcaltopo->lastHFRing() << std::endl;
#endif
 
  CaloTowerGeometry* geom = new CaloTowerGeometry(m_limits);
 
  if (nullptr == geom->cornersMgr())
    geom->allocateCorners(geom->numberOfCellsForCorners());
  if (nullptr == geom->parMgr())
    geom->allocatePar(geom->numberOfParametersPerShape() * geom->numberOfShapes(), geom->numberOfParametersPerShape());
 
  // simple loop
  for (uint32_t din = 0; din < m_limits->sizeForDenseIndexing(); ++din) {
    makeCell(din, geom);
  }
  edm::LogInfo("Geometry") << "CaloTowersHardcodeGeometry made " << m_limits->sizeForDenseIndexing() << " towers.";
  return std::unique_ptr<CaloSubdetectorGeometry>(geom);
}
 
void CaloTowerHardcodeGeometryLoader::makeCell(uint32_t din, CaloSubdetectorGeometry* geom) const {
  const double EBradius = 143.0;  // cm
  const double HOradius = 406.0 + 1.0;
  const double EEz = 320.0;  // rough (cm)
  const double HEz = 568.0;  // back (cm)
  const double HFz = 1100.0;
  const double HFthick = 165;
  // Tower 17 is the last EB tower
 
  //use CT topology to get the DetId for this dense index
  CaloTowerDetId id = m_limits->detIdFromDenseIndex(din);
  int ieta = id.ieta();
  int iphi = id.iphi();
 
  //use CT topology to get proper ieta for hcal
  int etaRing = m_limits->convertCTtoHcal(abs(ieta));
  int sign = (ieta > 0) ? (1) : (-1);
  double eta1, eta2;
 
#ifdef DebugLog
  std::cout << "CaloTowerHardcodeGeometryLoader: ieta = " << ieta << ", iphi = " << iphi << ", etaRing = " << etaRing
            << std::endl;
#endif
 
  if (abs(ieta) > m_limits->lastHERing()) {
    eta1 = theHFEtaBounds.at(etaRing - m_hcaltopo->firstHFRing());
    eta2 = theHFEtaBounds.at(etaRing - m_hcaltopo->firstHFRing() + 1);
  } else {
    eta1 = theHBHEEtaBounds.at(etaRing - 1);
    eta2 = theHBHEEtaBounds.at(etaRing);
  }
  double eta = 0.5 * (eta1 + eta2);
  double deta = (eta2 - eta1);
 
  // in radians
  double dphi_nominal = 2.0 * M_PI / m_hcaltopo->nPhiBins(1);  // always the same
  double dphi_half = M_PI / m_hcaltopo->nPhiBins(etaRing);     // half-width
 
  double phi_low = dphi_nominal * (iphi - 1);  // low-edge boundaries are constant...
  double phi = phi_low + dphi_half;
 
#ifdef DebugLog
  std::cout << "CaloTowerHardcodeGeometryLoader: eta1 = " << eta1 << ", eta2 = " << eta2 << ", eta = " << eta
            << ", phi = " << phi << std::endl;
#endif
 
  double x, y, z, thickness;
  bool alongZ = true;
  if (abs(ieta) > m_limits->lastHERing()) {  // forward
    z = HFz;
    double r = z / sinh(eta);
    x = r * cos(phi);
    y = r * sin(phi);
    thickness = HFthick / tanh(eta);
  } else if (abs(ieta) > m_limits->firstHERing() + 1) {  // EE-containing
    z = EEz;
    double r = z / sinh(eta);
    x = r * cos(phi);
    y = r * sin(phi);
    thickness = (HEz - EEz) / tanh(eta);
  } else {  // EB-containing
    x = EBradius * cos(phi);
    y = EBradius * sin(phi);
    alongZ = false;
    z = EBradius * sinh(eta);
    thickness = (HOradius - EBradius) * cosh(eta);
  }
 
  z *= sign;
  GlobalPoint point(x, y, z);
 
  const double mysign(!alongZ ? 1 : -1);
  std::vector<CCGFloat> hh;
  hh.reserve(5);
  hh.emplace_back(deta / 2);
  hh.emplace_back(dphi_half);
  hh.emplace_back(mysign * thickness / 2.);
 
  hh.emplace_back(fabs(eta));
  hh.emplace_back(fabs(z));
 
#ifdef DebugLog
  std::cout << "CaloTowerHardcodeGeometryLoader: x = " << x << ", y = " << y << ", z = " << z
            << ", thickness = " << thickness << std::endl;
#endif
 
  geom->newCell(point, point, point, CaloCellGeometry::getParmPtr(hh, geom->parMgr(), geom->parVecVec()), id);
}