HcalTowerAlgo/src/HcalDDDGeometryLoader.cc

0001 #include "Geometry/HcalTowerAlgo/interface/HcalDDDGeometryLoader.h"
0002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0003 #include "Geometry/HcalCommonData/interface/HcalParametersFromDD.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
0006 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include <string>
0010
0011 typedef CaloCellGeometry::CCGFloat CCGFloat;
0012
0013 //#define EDM_ML_DEBUG
0014
0015 HcalDDDGeometryLoader::HcalDDDGeometryLoader(const HcalDDDRecConstants* hcons) : hcalConstants_(hcons) {
0016   isBH_ = hcalConstants_->isBH();
0017 }
0018
0019 HcalDDDGeometryLoader::~HcalDDDGeometryLoader() {}
0020
0021 HcalDDDGeometryLoader::ReturnType HcalDDDGeometryLoader::load(const HcalTopology& topo,
0022                                                               DetId::Detector det,
0023                                                               int subdet) {
0024   HcalSubdetector hsub = static_cast<HcalSubdetector>(subdet);
0025
0026   HcalDDDGeometry* geom(new HcalDDDGeometry(topo));
0027
0028   if (geom->cornersMgr() == nullptr) {
0029     const unsigned int count(hcalConstants_->numberOfCells(HcalBarrel) + hcalConstants_->numberOfCells(HcalEndcap) +
0030                              2 * hcalConstants_->numberOfCells(HcalForward) + hcalConstants_->numberOfCells(HcalOuter));
0031     geom->allocateCorners(count);
0032   }
0033
0034   //  if( geom->cornersMgr() == 0 )  geom->allocateCorners( 2592 ) ;
0035
0036   if (geom->parMgr() == nullptr)
0037     geom->allocatePar(500, 3);
0038
0039   fill(hsub, geom);
0040   //fast insertion of valid ids requires sort at end
0041   geom->sortValidIds();
0042   return geom;
0043 }
0044
0045 HcalDDDGeometryLoader::ReturnType HcalDDDGeometryLoader::load(const HcalTopology& topo) {
0046   HcalDDDGeometry* geom(new HcalDDDGeometry(topo));
0047
0048   if (geom->cornersMgr() == nullptr) {
0049     const unsigned int count(hcalConstants_->numberOfCells(HcalBarrel) + hcalConstants_->numberOfCells(HcalEndcap) +
0050                              2 * hcalConstants_->numberOfCells(HcalForward) + hcalConstants_->numberOfCells(HcalOuter));
0051     geom->allocateCorners(count);
0052   }
0053   if (geom->parMgr() == nullptr)
0054     geom->allocatePar(500, 3);
0055
0056   fill(HcalBarrel, geom);
0057   fill(HcalEndcap, geom);
0058   fill(HcalForward, geom);
0059   fill(HcalOuter, geom);
0060   //fast insertion of valid ids requires sort at end
0061   geom->sortValidIds();
0062   return geom;
0063 }
0064
0065 void HcalDDDGeometryLoader::fill(HcalSubdetector subdet, HcalDDDGeometry* geom) {
0066   // start by making the new HcalDetIds
0067   std::vector<HcalCellType> hcalCells = hcalConstants_->HcalCellTypes(subdet);
0068   geom->insertCell(hcalCells);
0069 #ifdef EDM_ML_DEBUG
0070   edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::fill gets " << hcalCells.size() << " cells for subdetector "
0071                                << subdet;
0072 #endif
0073   // Make the new HcalDetIds and the cells
0074
0075   std::vector<HcalDetId> hcalIds;
0076   for (auto& hcalCell : hcalCells) {
0077     int etaRing = hcalCell.etaBin();
0078     int iside = hcalCell.zside();
0079     int depthBin = hcalCell.depthSegment();
0080     double dphi = hcalCell.phiBinWidth();
0081     std::vector<std::pair<int, double> > phis = hcalCell.phis();
0082 #ifdef EDM_ML_DEBUG
0083     edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader: Subdet " << subdet << " side " << iside << " eta "
0084                                  << etaRing << " depth " << depthBin << " with " << phis.size() << "modules:";
0085     size_t i(0);
0086 #endif
0087     geom->increaseReserve(phis.size());
0088     for (auto& phi : phis) {
0089 #ifdef EDM_ML_DEBUG
0090       edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::fill Cell " << i << " eta " << iside * etaRing << " phi "
0091                                    << phi.first << "(" << phi.second / CLHEP::deg << ", " << dphi / CLHEP::deg
0092                                    << ") depth " << depthBin;
0093       i++;
0094 #endif
0095       HcalDetId id(subdet, iside * etaRing, phi.first, depthBin);
0096       hcalIds.emplace_back(id);
0097       makeCell(id, hcalCell, phi.second, dphi, geom);
0098     }
0099   }
0100
0101   edm::LogInfo("HCalGeom") << "Number of HCAL DetIds made for " << subdet << " is " << hcalIds.size();
0102 }
0103
0104 void HcalDDDGeometryLoader::makeCell(
0105     const HcalDetId& detId, const HcalCellType& hcalCell, double phi, double dphi, HcalDDDGeometry* geom) const {
0106   // the two eta boundaries of the cell
0107   double eta1 = hcalCell.etaMin();
0108   double eta2 = hcalCell.etaMax();
0109   HcalSubdetector subdet = detId.subdet();
0110   double eta = 0.5 * (eta1 + eta2) * detId.zside();
0111   double deta = (eta2 - eta1);
0112   double theta = 2.0 * atan(exp(-eta));
0113
0114   // barrel vs forward
0115   bool rzType = hcalCell.depthType();
0116   bool isBarrel = (subdet == HcalBarrel || subdet == HcalOuter);
0117
0118   double z, r, thickness;
0119 #ifdef EDM_ML_DEBUG
0120   double r0, r1, r2;
0121 #endif
0122   if (rzType) {
0123     r = hcalCell.depthMin();
0124     if (isBarrel) {
0125       z = r * sinh(eta);                                  // sinh(eta) == 1/tan(theta)
0126       thickness = (hcalCell.depthMax() - r) * cosh(eta);  // cosh(eta) == 1/sin(theta)
0127 #ifdef EDM_ML_DEBUG
0128       r1 = r;
0129       r2 = hcalCell.depthMax();
0130       r0 = 0.5 * (r1 + r2);
0131 #endif
0132     } else {
0133       z = r * sinh(eta2);
0134       thickness = 2. * hcalCell.halfSize();
0135       r = z / sinh(std::abs(eta));
0136 #ifdef EDM_ML_DEBUG
0137       r0 = z / sinh(std::abs(eta));
0138       r1 = z / sinh(std::abs(eta) + 0.5 * deta);
0139       r2 = z / sinh(std::abs(eta) - 0.5 * deta);
0140 #endif
0141     }
0142 #ifdef EDM_ML_DEBUG
0143     edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::makeCell SubDet " << subdet << " eta = " << eta
0144                                  << " theta = " << theta << " r = " << r << " thickness = " << thickness << " r0-r2 ("
0145                                  << r0 << ":" << r1 << ":" << r2 << ")";
0146 #endif
0147   } else {
0148     z = hcalCell.depthMin();
0149     thickness = hcalCell.depthMax() - z;
0150     if (isBH_)
0151       z += (0.5 * thickness);
0152     z *= detId.zside();  // get the sign right.
0153     r = z * tan(theta);
0154     thickness /= std::abs(cos(theta));
0155 #ifdef EDM_ML_DEBUG
0156     r0 = z / sinh(std::abs(eta));
0157     r1 = z / sinh(std::abs(eta) + 0.5 * deta);
0158     r2 = z / sinh(std::abs(eta) - 0.5 * deta);
0159     edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::makeCell SubDet " << subdet << " eta = " << eta
0160                                  << " theta = " << theta << " z = " << z << " r = " << r << " thickness = " << thickness
0161                                  << " r0-r2 (" << r0 << ":" << r1 << ":" << r2 << ")";
0162 #endif
0163   }
0164
0165   double x = r * cos(phi);
0166   double y = r * sin(phi);
0167   GlobalPoint point(x, y, z);
0168
0169 #ifdef EDM_ML_DEBUG
0170   edm::LogVerbatim("HCalGeom") << "HcalDDDGeometryLoader::makeCell for " << detId << " Point " << point
0171                                << " deta = " << deta << " dphi = " << dphi << " thickness = " << thickness
0172                                << " isBarrel = " << isBarrel << " " << rzType;
0173 #endif
0174
0175   std::vector<CCGFloat> hp;
0176   hp.reserve(3);
0177
0178   if (subdet == HcalForward) {
0179     hp.emplace_back(deta / 2.);
0180     hp.emplace_back(dphi / 2.);
0181     hp.emplace_back(thickness / 2.);
0182   } else {
0183     const double sign(isBarrel ? 1 : -1);
0184     hp.emplace_back(deta / 2.);
0185     hp.emplace_back(dphi / 2.);
0186     hp.emplace_back(sign * thickness / 2.);
0187   }
0188   geom->newCellFast(point, point, point, CaloCellGeometry::getParmPtr(hp, geom->parMgr(), geom->parVecVec()), detId);
0189 }