HcalTowerAlgo/src/CaloTowerHardcodeGeometryLoader.cc

0001 #include "Geometry/HcalTowerAlgo/interface/CaloTowerHardcodeGeometryLoader.h"
0002 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0003 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0004 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include <algorithm>
0007 #include <iostream>
0008 #include <iterator>
0009 #include <sstream>
0010
0011 //#define EDM_ML_DEUG
0012
0013 typedef CaloCellGeometry::CCGFloat CCGFloat;
0014
0015 std::unique_ptr<CaloSubdetectorGeometry> CaloTowerHardcodeGeometryLoader::load(const CaloTowerTopology* limits,
0016                                                                                const HcalTopology* hcaltopo,
0017                                                                                const HcalDDDRecConstants* hcons) {
0018   m_limits = limits;
0019   m_hcaltopo = hcaltopo;
0020   m_hcons = hcons;
0021
0022   //get eta limits from hcal rec constants
0023   theHFEtaBounds = m_hcons->getEtaTableHF();
0024   theHBHEEtaBounds = m_hcons->getEtaTable();
0025
0026 #ifdef EDM_ML_DEUG
0027   std::ostringstream st1;
0028   st1 << "CaloTowerHardcodeGeometryLoader: theHBHEEtaBounds = ";
0029   std::copy(theHBHEEtaBounds.begin(), theHBHEEtaBounds.end(), std::ostream_iterator<double>(st1, ","));
0030   edm::LogVerbatim("HCalGeom") << st1.str();
0031
0032   edm::LogVerbatim("HCalGeom") << "CaloTowerHardcodeGeometryLoader: lastHBRing = " << m_limits->lastHBRing()
0033                                << ", lastHERing = " << m_limits->lastHERing();
0034   edm::LogVerbatim("HCalGeom") << "CaloTowerHardcodeGeometryLoader: HcalTopology: firstHBRing = "
0035                                << hcaltopo->firstHBRing() << ", lastHBRing = " << hcaltopo->lastHBRing()
0036                                << ", firstHERing = " << hcaltopo->firstHERing()
0037                                << ", lastHERing = " << hcaltopo->lastHERing()
0038                                << ", firstHFRing = " << hcaltopo->firstHFRing()
0039                                << ", lastHFRing = " << hcaltopo->lastHFRing();
0040 #endif
0041
0042   CaloTowerGeometry* geom = new CaloTowerGeometry(m_limits);
0043
0044   if (nullptr == geom->cornersMgr())
0045     geom->allocateCorners(geom->numberOfCellsForCorners());
0046   if (nullptr == geom->parMgr())
0047     geom->allocatePar(geom->numberOfParametersPerShape() * geom->numberOfShapes(), geom->numberOfParametersPerShape());
0048
0049   // simple loop
0050   for (uint32_t din = 0; din < m_limits->sizeForDenseIndexing(); ++din) {
0051     makeCell(din, geom);
0052   }
0053   edm::LogInfo("Geometry") << "CaloTowersHardcodeGeometry made " << m_limits->sizeForDenseIndexing() << " towers.";
0054   return std::unique_ptr<CaloSubdetectorGeometry>(geom);
0055 }
0056
0057 void CaloTowerHardcodeGeometryLoader::makeCell(uint32_t din, CaloSubdetectorGeometry* geom) const {
0058   const double EBradius = 143.0;  // cm
0059   const double HOradius = 406.0 + 1.0;
0060   const double EEz = 320.0;  // rough (cm)
0061   const double HEz = 568.0;  // back (cm)
0062   const double HFz = 1100.0;
0063   const double HFthick = 165;
0064   // Tower 17 is the last EB tower
0065
0066   //use CT topology to get the DetId for this dense index
0067   CaloTowerDetId id = m_limits->detIdFromDenseIndex(din);
0068   int ieta = id.ieta();
0069   int iphi = id.iphi();
0070
0071   //use CT topology to get proper ieta for hcal
0072   int etaRing = m_limits->convertCTtoHcal(abs(ieta));
0073   int sign = (ieta > 0) ? (1) : (-1);
0074   double eta1, eta2;
0075
0076 #ifdef EDM_ML_DEUG
0077   edm::LogVerbatim("HCalGeom") << "CaloTowerHardcodeGeometryLoader: ieta = " << ieta << ", iphi = " << iphi
0078                                << ", etaRing = " << etaRing;
0079 #endif
0080
0081   if (abs(ieta) > m_limits->lastHERing()) {
0082     eta1 = theHFEtaBounds.at(etaRing - m_hcaltopo->firstHFRing());
0083     eta2 = theHFEtaBounds.at(etaRing - m_hcaltopo->firstHFRing() + 1);
0084   } else {
0085     eta1 = theHBHEEtaBounds.at(etaRing - 1);
0086     eta2 = theHBHEEtaBounds.at(etaRing);
0087   }
0088   double eta = 0.5 * (eta1 + eta2);
0089   double deta = (eta2 - eta1);
0090
0091   // in radians
0092   double dphi_nominal = 2.0 * M_PI / m_hcaltopo->nPhiBins(1);  // always the same
0093   double dphi_half = M_PI / m_hcaltopo->nPhiBins(etaRing);     // half-width
0094
0095   double phi_low = dphi_nominal * (iphi - 1);  // low-edge boundaries are constant...
0096   double phi = phi_low + dphi_half;
0097
0098 #ifdef EDM_ML_DEUG
0099   edm::LogVerbatim("HCalGeom") << "CaloTowerHardcodeGeometryLoader: eta1 = " << eta1 << ", eta2 = " << eta2
0100                                << ", eta = " << eta << ", phi = " << phi;
0101 #endif
0102
0103   double x, y, z, thickness;
0104   bool alongZ = true;
0105   if (abs(ieta) > m_limits->lastHERing()) {  // forward
0106     z = HFz;
0107     double r = z / sinh(eta);
0108     x = r * cos(phi);
0109     y = r * sin(phi);
0110     thickness = HFthick / tanh(eta);
0111   } else if (abs(ieta) > m_limits->firstHERing() + 1) {  // EE-containing
0112     z = EEz;
0113     double r = z / sinh(eta);
0114     x = r * cos(phi);
0115     y = r * sin(phi);
0116     thickness = (HEz - EEz) / tanh(eta);
0117   } else {  // EB-containing
0118     x = EBradius * cos(phi);
0119     y = EBradius * sin(phi);
0120     alongZ = false;
0121     z = EBradius * sinh(eta);
0122     thickness = (HOradius - EBradius) * cosh(eta);
0123   }
0124
0125   z *= sign;
0126   GlobalPoint point(x, y, z);
0127
0128   const double mysign(!alongZ ? 1 : -1);
0129   std::vector<CCGFloat> hh;
0130   hh.reserve(5);
0131   hh.emplace_back(deta / 2);
0132   hh.emplace_back(dphi_half);
0133   hh.emplace_back(mysign * thickness / 2.);
0134
0135   hh.emplace_back(fabs(eta));
0136   hh.emplace_back(fabs(z));
0137
0138 #ifdef EDM_ML_DEUG
0139   edm::LogVerbatim("HCalGeom") << "CaloTowerHardcodeGeometryLoader: x = " << x << ", y = " << y << ", z = " << z
0140                                << ", thickness = " << thickness;
0141 #endif
0142
0143   geom->newCell(point, point, point, CaloCellGeometry::getParmPtr(hh, geom->parMgr(), geom->parVecVec()), id);
0144 }