HcalFlexiHardcodeGeometryLoader.cc

Go to the documentation of this file.

#include "Geometry/HcalTowerAlgo/interface/HcalFlexiHardcodeGeometryLoader.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
#include "Geometry/CaloGeometry/interface/IdealZPrism.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
#include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <algorithm>
#include <sstream>
#include <vector>

using CCGFloat = CaloCellGeometry::CCGFloat;

//#define EDM_ML_DEBUG
// ==============> Loader Itself <==========================

HcalFlexiHardcodeGeometryLoader::HcalFlexiHardcodeGeometryLoader() {
  MAX_HCAL_PHI = 72;
  DEGREE2RAD = M_PI / 180.;
}

CaloSubdetectorGeometry* HcalFlexiHardcodeGeometryLoader::load(const HcalTopology& fTopology,
                                                               const HcalDDDRecConstants& hcons) {
  HcalGeometry* hcalGeometry = new HcalGeometry(fTopology);
  if (nullptr == hcalGeometry->cornersMgr())
    hcalGeometry->allocateCorners(fTopology.ncells() + fTopology.getHFSize());
  if (nullptr == hcalGeometry->parMgr())
    hcalGeometry->allocatePar(hcalGeometry->numberOfShapes(), HcalGeometry::k_NumberOfParametersPerShape);
  isBH_ = hcons.isBH();
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader initialize with ncells " << fTopology.ncells() << " and shapes "
                               << hcalGeometry->numberOfShapes() << ":" << HcalGeometry::k_NumberOfParametersPerShape
                               << " with BH Flag " << isBH_;
#endif
  if (fTopology.mode() == HcalTopologyMode::H2) {  // TB geometry
    fillHBHO(hcalGeometry, makeHBCells(hcons), true);
    fillHBHO(hcalGeometry, makeHOCells(), false);
    fillHE(hcalGeometry, makeHECells_H2());
  } else {  // regular geometry
    fillHBHO(hcalGeometry, makeHBCells(hcons), true);
    fillHBHO(hcalGeometry, makeHOCells(), false);
    fillHF(hcalGeometry, makeHFCells(hcons));
    fillHE(hcalGeometry, makeHECells(hcons));
  }
  //fast insertion of valid ids requires sort at end
  hcalGeometry->sortValidIds();
  return hcalGeometry;
}

// ----------> HB <-----------
std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> HcalFlexiHardcodeGeometryLoader::makeHBCells(
    const HcalDDDRecConstants& hcons) {
  std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> result;
  std::vector<std::pair<double, double> > gconsHB = hcons.getConstHBHE(0);
  std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = hcons.getEtaBins(0);

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader called for " << etabins.size() << " Eta Bins";
  for (unsigned int k = 0; k < gconsHB.size(); ++k) {
    edm::LogVerbatim("HCalGeom") << "gconsHB[" << k << "] = " << gconsHB[k].first << " +- " << gconsHB[k].second;
  }
#endif
  for (auto& etabin : etabins) {
    int iring = (etabin.zside >= 0) ? etabin.ieta : -etabin.ieta;
    int depth = etabin.depthStart;
    double dphi =
        (etabin.phis.size() > 1) ? (etabin.phis[1].second - etabin.phis[0].second) : ((2.0 * M_PI) / MAX_HCAL_PHI);
    for (unsigned int k = 0; k < etabin.layer.size(); ++k) {
      int layf = etabin.layer[k].first - 1;
      int layl = etabin.layer[k].second - 1;
      double rmin = gconsHB[layf].first - gconsHB[layf].second;
      double rmax = gconsHB[layl].first + gconsHB[layl].second;
      for (unsigned int j = 0; j < etabin.phis.size(); ++j) {
#ifdef EDM_ML_DEBUG
        edm::LogVerbatim("HCalGeom") << "HBRing " << iring << " eta " << etabin.etaMin << ":" << etabin.etaMax
                                     << " depth " << depth << " R " << rmin << ":" << rmax << " Phi "
                                     << etabin.phis[j].first << ":" << etabin.phis[j].second << ":" << dphi << " layer["
                                     << k << "]: " << etabin.layer[k].first - 1 << ":" << etabin.layer[k].second;
#endif
        result.emplace_back(HcalFlexiHardcodeGeometryLoader::HBHOCellParameters(
            iring, depth, etabin.phis[j].first, etabin.phis[j].second, dphi, rmin, rmax, etabin.etaMin, etabin.etaMax));
      }
      depth++;
    }
  }
  return result;
}

// ----------> HO <-----------
std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> HcalFlexiHardcodeGeometryLoader::makeHOCells() {
  const double HORMIN0 = 390.0;
  const double HORMIN1 = 412.6;
  const double HORMAX = 413.6;
  const int nCells = 15;
  const int nPhi = 72;
  const double etamin[nCells] = {
      0.000, 0.087, 0.174, 0.261, 0.3395, 0.435, 0.522, 0.609, 0.696, 0.783, 0.873, 0.957, 1.044, 1.131, 1.218};
  const double etamax[nCells] = {
      0.087, 0.174, 0.261, 0.3075, 0.435, 0.522, 0.609, 0.696, 0.783, 0.8494, 0.957, 1.044, 1.131, 1.218, 1.305};
  std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> result;
  result.reserve(nCells * nPhi);
  double dphi = ((2.0 * M_PI) / nPhi);
  for (int i = 0; i < nCells; ++i) {
    double rmin = ((i < 4) ? HORMIN0 : HORMIN1);
    for (int iside = -1; iside <= 1; iside += 2) {
      for (int j = 0; j < nPhi; ++j) {
        double phi = (j + 0.5) * dphi;
        // eta, depth, phi, phi0, deltaPhi, rMin, rMax, etaMin, etaMax
        result.emplace_back(HcalFlexiHardcodeGeometryLoader::HBHOCellParameters(
            iside * (i + 1), 4, j + 1, phi, dphi, rmin, HORMAX, etamin[i], etamax[i]));
      }
    }
  }
  return result;
}

//
// Convert constants to appropriate cells
//
void HcalFlexiHardcodeGeometryLoader::fillHBHO(
    HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters>& fCells, bool fHB) {
  fGeometry->increaseReserve(fCells.size());
  for (const auto& param : fCells) {
    HcalDetId hid(fHB ? HcalBarrel : HcalOuter, param.ieta, param.iphi, param.depth);
    float phiCenter = param.phi;  // middle of the cell
    float etaCenter = 0.5 * (param.etaMin + param.etaMax);
    float x = param.rMin * cos(phiCenter);
    float y = param.rMin * sin(phiCenter);
    float z = (param.ieta < 0) ? -(param.rMin * sinh(etaCenter)) : (param.rMin * sinh(etaCenter));
    // make cell geometry
    GlobalPoint refPoint(x, y, z);  // center of the cell's face
    std::vector<CCGFloat> cellParams;
    cellParams.reserve(5);
    cellParams.emplace_back(0.5 * (param.etaMax - param.etaMin));                // deta_half
    cellParams.emplace_back(0.5 * param.dphi);                                   // dphi_half
    cellParams.emplace_back(0.5 * (param.rMax - param.rMin) * cosh(etaCenter));  // dr_half
    cellParams.emplace_back(fabs(refPoint.eta()));
    cellParams.emplace_back(fabs(refPoint.z()));
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHBHO-> " << hid << " " << hid.rawId() << " "
                                 << std::hex << hid.rawId() << std::dec << " " << hid << " " << refPoint << '/'
                                 << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
#endif
    fGeometry->newCellFast(refPoint,
                           refPoint,
                           refPoint,
                           CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
                           hid);
  }
}

// ----------> HE <-----------
std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> HcalFlexiHardcodeGeometryLoader::makeHECells(
    const HcalDDDRecConstants& hcons) {
  std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> result;
  std::vector<std::pair<double, double> > gconsHE = hcons.getConstHBHE(1);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader:HE with " << gconsHE.size() << " cells";
#endif
  if (!gconsHE.empty()) {
    std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = hcons.getEtaBins(1);

#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader called for HE with " << etabins.size() << " Eta Bins and "
                                 << gconsHE.size() << " depths";
    std::ostringstream st1;
    for (unsigned int i = 0; i < gconsHE.size(); ++i)
      st1 << " Depth[" << i << "] = " << gconsHE[i].first << " +- " << gconsHE[i].second;
    edm::LogVerbatim("HCalGeom") << st1.str();
#endif
    for (auto& etabin : etabins) {
      int iring = (etabin.zside >= 0) ? etabin.ieta : -etabin.ieta;
      int depth = etabin.depthStart;
      double dphi =
          (etabin.phis.size() > 1) ? (etabin.phis[1].second - etabin.phis[0].second) : ((2.0 * M_PI) / MAX_HCAL_PHI);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader::Ring " << iring << " nphi " << etabin.phis.size()
                                   << " dstart " << depth << " dphi " << dphi << " layers " << etabin.layer.size();
      std::ostringstream st2;
      for (unsigned int j = 0; j < etabin.phis.size(); ++j)
        st2 << " [" << j << "] " << etabin.phis[j].first << ":" << etabin.phis[j].second;
      edm::LogVerbatim("HCalGeom") << st2.str();
#endif
      for (unsigned int k = 0; k < etabin.layer.size(); ++k) {
        int layf = etabin.layer[k].first - 1;
        int layl = etabin.layer[k].second - 1;
        double zmin = gconsHE[layf].first - gconsHE[layf].second;
        double zmax = gconsHE[layl].first + gconsHE[layl].second;
        if (zmin < 1.0) {
          for (int k2 = layf; k2 <= layl; ++k2) {
            if (gconsHE[k2].first > 10) {
              zmin = gconsHE[k2].first - gconsHE[k2].second;
              break;
            }
          }
        }
        if (zmin >= zmax)
          zmax = zmin + 10.;
        for (unsigned int j = 0; j < etabin.phis.size(); ++j) {
#ifdef EDM_ML_DEBUG
          edm::LogVerbatim("HCalGeom") << "HERing " << iring << " eta " << etabin.etaMin << ":" << etabin.etaMax
                                       << " depth " << depth << " Z " << zmin << ":" << zmax
                                       << " Phi :" << etabin.phis[j].first << ":" << etabin.phis[j].second << ":"
                                       << dphi << " layer[" << k << "]: " << etabin.layer[k].first - 1 << ":"
                                       << etabin.layer[k].second - 1;
#endif
          result.emplace_back(HcalFlexiHardcodeGeometryLoader::HECellParameters(
              iring, depth, etabin.phis[j].first, etabin.phis[j].second, dphi, zmin, zmax, etabin.etaMin, etabin.etaMax));
        }
        depth++;
      }
    }
  }
  return result;
}

// ----------> HE @ H2 <-----------
std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> HcalFlexiHardcodeGeometryLoader::makeHECells_H2() {
  const double HEZMIN_H2 = 400.715;
  const double HEZMID_H2 = 436.285;
  const double HEZMAX_H2 = 541.885;
  const int nEtas = 10;
  const int nDepth[nEtas] = {1, 2, 2, 2, 2, 2, 2, 2, 3, 3};
  const int dStart[nEtas] = {3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
  const int nPhis[nEtas] = {8, 8, 8, 8, 8, 8, 4, 4, 4, 4};
  const double etas[nEtas + 1] = {1.305, 1.373, 1.444, 1.521, 1.603, 1.693, 1.790, 1.880, 1.980, 2.090, 2.210};
  const double zval[4 * nEtas] = {
      409.885,   462.685,   0.,        0.,        HEZMIN_H2, 427.485,   506.685,   0.0,       HEZMIN_H2, HEZMID_H2,
      524.285,   0.,        HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,        HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,
      HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,        HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,        HEZMIN_H2, HEZMID_H2,
      HEZMAX_H2, 0.,        HEZMIN_H2, 418.685,   HEZMID_H2, HEZMAX_H2, HEZMIN_H2, 418.685,   HEZMID_H2, HEZMAX_H2};
  std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> result;

  for (int i = 0; i < nEtas; ++i) {
    int ieta = 16 + i;
    for (int k = 0; k < nDepth[i]; ++k) {
      int depth = dStart[i] + k;
      for (int j = 0; j < nPhis[i]; ++j) {
        int iphi = (nPhis[i] == 8) ? (j + 1) : (2 * j + 1);
        double dphi = (40.0 * DEGREE2RAD) / nPhis[i];
        double phi0 = (j + 0.5) * dphi;
        // ieta, depth, iphi, phi0, deltaPhi, zMin, zMax, etaMin, etaMax
        result.emplace_back(HcalFlexiHardcodeGeometryLoader::HECellParameters(
            ieta, depth, iphi, phi0, dphi, zval[4 * i + k + 1], zval[4 * i + k + 2], etas[i], etas[i + 1]));
      }
    }
  }
  return result;
}

// ----------> HF <-----------
std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters> HcalFlexiHardcodeGeometryLoader::makeHFCells(
    const HcalDDDRecConstants& hcons) {
  const float HFZMIN1 = 1115.;
  const float HFZMIN2 = 1137.;
  const float HFZMAX = 1280.1;
  std::vector<HcalDDDRecConstants::HFCellParameters> cells = hcons.getHFCellParameters();
  unsigned int nCells = cells.size();
  std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters> result;
  result.reserve(nCells);
  for (unsigned int i = 0; i < nCells; ++i) {
    HcalFlexiHardcodeGeometryLoader::HFCellParameters cell1(cells[i].ieta,
                                                            cells[i].depth,
                                                            cells[i].firstPhi,
                                                            cells[i].stepPhi,
                                                            cells[i].nPhi,
                                                            5 * cells[i].stepPhi,
                                                            HFZMIN1,
                                                            HFZMAX,
                                                            cells[i].rMin,
                                                            cells[i].rMax);
    result.emplace_back(cell1);
    HcalFlexiHardcodeGeometryLoader::HFCellParameters cell2(cells[i].ieta,
                                                            1 + cells[i].depth,
                                                            cells[i].firstPhi,
                                                            cells[i].stepPhi,
                                                            cells[i].nPhi,
                                                            5 * cells[i].stepPhi,
                                                            HFZMIN2,
                                                            HFZMAX,
                                                            cells[i].rMin,
                                                            cells[i].rMax);
    result.emplace_back(cell2);
  }
  return result;
}

void HcalFlexiHardcodeGeometryLoader::fillHE(
    HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters>& fCells) {
  fGeometry->increaseReserve(fCells.size());
  for (const auto& param : fCells) {
    HcalDetId hid(HcalEndcap, param.ieta, param.iphi, param.depth);
    float phiCenter = param.phi;  // middle of the cell
    float etaCenter = 0.5 * (param.etaMin + param.etaMax);
    int iside = (param.ieta >= 0) ? 1 : -1;
    float z = (isBH_) ? (iside * 0.5 * (param.zMin + param.zMax)) : (iside * param.zMin);
    float perp = std::abs(z) / sinh(etaCenter);
    float x = perp * cos(phiCenter);
    float y = perp * sin(phiCenter);
    // make cell geometry
    GlobalPoint refPoint(x, y, z);  // center of the cell's face
    std::vector<CCGFloat> cellParams;
    cellParams.reserve(5);
    cellParams.emplace_back(0.5 * (param.etaMax - param.etaMin));                 //deta_half
    cellParams.emplace_back(0.5 * param.dphi);                                    // dphi_half
    cellParams.emplace_back(-0.5 * (param.zMax - param.zMin) / tanh(etaCenter));  // dz_half, "-" means edges in Z
    cellParams.emplace_back(fabs(refPoint.eta()));
    cellParams.emplace_back(fabs(refPoint.z()));
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHE-> " << hid << " " << hid.rawId() << " "
                                 << std::hex << hid.rawId() << std::dec << " " << hid << refPoint << '/'
                                 << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
#endif
    fGeometry->newCellFast(refPoint,
                           refPoint,
                           refPoint,
                           CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
                           hid);
  }
}

void HcalFlexiHardcodeGeometryLoader::fillHF(
    HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters>& fCells) {
  fGeometry->increaseReserve(fCells.size());
  for (const auto& param : fCells) {
    for (int kPhi = 0; kPhi < param.nPhi; ++kPhi) {
      int iPhi = param.phiFirst + kPhi * param.phiStep;
      HcalDetId hid(HcalForward, param.eta, iPhi, param.depth);
      // middle of the cell
      float phiCenter = ((iPhi - 1) * 360. / MAX_HCAL_PHI + 0.5 * param.dphi) * DEGREE2RAD;
      GlobalPoint inner(param.rMin, 0, param.zMin);
      GlobalPoint outer(param.rMax, 0, param.zMin);
      float iEta = inner.eta();
      float oEta = outer.eta();
      float etaCenter = 0.5 * (iEta + oEta);

      float perp = param.zMin / sinh(etaCenter);
      float x = perp * cos(phiCenter);
      float y = perp * sin(phiCenter);
      float z = (param.eta > 0) ? param.zMin : -param.zMin;
      // make cell geometry
      GlobalPoint refPoint(x, y, z);  // center of the cell's face
      std::vector<CCGFloat> cellParams;
      cellParams.reserve(5);
      cellParams.emplace_back(0.5 * (iEta - oEta));              // deta_half
      cellParams.emplace_back(0.5 * param.dphi * DEGREE2RAD);    // dphi_half
      cellParams.emplace_back(0.5 * (param.zMax - param.zMin));  // dz_half
      cellParams.emplace_back(fabs(refPoint.eta()));
      cellParams.emplace_back(fabs(refPoint.z()));
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHF-> " << hid << " " << hid.rawId() << " "
                                   << std::hex << hid.rawId() << std::dec << " " << hid << " " << refPoint << '/'
                                   << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
#endif
      fGeometry->newCellFast(refPoint,
                             refPoint,
                             refPoint,
                             CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
                             hid);
    }
  }
}