HcalTowerAlgo/src/HcalGeometry.cc

0001 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0004 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0005 #include <algorithm>
0006 #include <iostream>
0007
0008 #include <Math/Transform3D.h>
0009 #include <Math/EulerAngles.h>
0010
0011 typedef CaloCellGeometry::CCGFloat CCGFloat;
0012 typedef CaloCellGeometry::Pt3D Pt3D;
0013 typedef CaloCellGeometry::Pt3DVec Pt3DVec;
0014 typedef CaloCellGeometry::Tr3D Tr3D;
0015
0016 //#define EDM_ML_DEBUG
0017
0018 HcalGeometry::HcalGeometry(const HcalTopology& topology)
0019     : m_topology(topology), m_mergePosition(topology.getMergePositionFlag()) {
0020   init();
0021 }
0022
0023 HcalGeometry::~HcalGeometry() {}
0024
0025 void HcalGeometry::init() {
0026   if (!m_topology.getMergePositionFlag())
0027     m_mergePosition = false;
0028 #ifdef EDM_ML_DEBUG
0029   edm::LogVerbatim("HCalGeom") << "HcalGeometry_init(): HBSize " << m_topology.getHBSize() << " HESize "
0030                                << m_topology.getHESize() << " HOSize " << m_topology.getHOSize() << " HFSize "
0031                                << m_topology.getHFSize();
0032 #endif
0033
0034   m_hbCellVec = HBCellVec(m_topology.getHBSize());
0035   m_heCellVec = HECellVec(m_topology.getHESize());
0036   m_hoCellVec = HOCellVec(m_topology.getHOSize());
0037   m_hfCellVec = HFCellVec(m_topology.getHFSize());
0038 }
0039
0040 void HcalGeometry::fillDetIds() const {
0041   // this must test the last record filled to avoid a race condition
0042   if (!m_emptyIds.isSet()) {
0043     std::unique_ptr<std::vector<DetId>> p_hbIds{new std::vector<DetId>};
0044     std::unique_ptr<std::vector<DetId>> p_heIds{new std::vector<DetId>};
0045     std::unique_ptr<std::vector<DetId>> p_hoIds{new std::vector<DetId>};
0046     std::unique_ptr<std::vector<DetId>> p_hfIds{new std::vector<DetId>};
0047     std::unique_ptr<std::vector<DetId>> p_emptyIds{new std::vector<DetId>};
0048
0049     const std::vector<DetId>& baseIds(CaloSubdetectorGeometry::getValidDetIds());
0050     for (unsigned int i(0); i != baseIds.size(); ++i) {
0051       const DetId id(baseIds[i]);
0052       if (id.subdetId() == HcalBarrel) {
0053         p_hbIds->emplace_back(id);
0054       } else if (id.subdetId() == HcalEndcap) {
0055         p_heIds->emplace_back(id);
0056       } else if (id.subdetId() == HcalOuter) {
0057         p_hoIds->emplace_back(id);
0058       } else if (id.subdetId() == HcalForward) {
0059         p_hfIds->emplace_back(id);
0060       }
0061     }
0062     std::sort(p_hbIds->begin(), p_hbIds->end());
0063     std::sort(p_heIds->begin(), p_heIds->end());
0064     std::sort(p_hoIds->begin(), p_hoIds->end());
0065     std::sort(p_hfIds->begin(), p_hfIds->end());
0066     p_emptyIds->resize(0);
0067
0068     m_hbIds.set(std::move(p_hbIds));
0069     m_heIds.set(std::move(p_heIds));
0070     m_hoIds.set(std::move(p_hoIds));
0071     m_hfIds.set(std::move(p_hfIds));
0072     m_emptyIds.set(std::move(p_emptyIds));
0073   }
0074 }
0075
0076 const std::vector<DetId>& HcalGeometry::getValidDetIds(DetId::Detector det, int subdet) const {
0077   if (0 != subdet)
0078     fillDetIds();
0079   return (0 == subdet
0080               ? CaloSubdetectorGeometry::getValidDetIds()
0081               : (HcalBarrel == subdet
0082                      ? *m_hbIds.load()
0083                      : (HcalEndcap == subdet
0084                             ? *m_heIds.load()
0085                             : (HcalOuter == subdet ? *m_hoIds.load()
0086                                                    : (HcalForward == subdet ? *m_hfIds.load() : *m_emptyIds.load())))));
0087 }
0088
0089 std::shared_ptr<const CaloCellGeometry> HcalGeometry::getGeometry(const DetId& id) const {
0090 #ifdef EDM_ML_DEBUG
0091   std::ostringstream st1;
0092   st1 << "HcalGeometry::getGeometry for " << HcalDetId(id) << "  " << m_mergePosition << " ";
0093 #endif
0094   if (!m_mergePosition) {
0095 #ifdef EDM_ML_DEBUG
0096     st1 << m_topology.detId2denseId(id) << " " << getGeometryBase(id) << "\n";
0097     edm::LogVerbatim("HCalGeom") << st1.str();
0098 #endif
0099     return getGeometryBase(id);
0100   } else {
0101 #ifdef EDM_ML_DEBUG
0102     st1 << m_topology.detId2denseId(m_topology.idFront(id)) << " " << getGeometryBase(m_topology.idFront(id));
0103     edm::LogVerbatim("HCalGeom") << st1.str();
0104 #endif
0105     return getGeometryBase(m_topology.idFront(id));
0106   }
0107 }
0108
0109 DetId HcalGeometry::getClosestCell(const GlobalPoint& r) const { return getClosestCell(r, false); }
0110
0111 DetId HcalGeometry::getClosestCell(const GlobalPoint& r, bool ignoreCorrect) const {
0112   // Now find the closest eta_bin, eta value of a bin i is average
0113   // of eta[i] and eta[i-1]
0114   static const double z_long = 1100.0;
0115   double abseta = fabs(r.eta());
0116   double absz = fabs(r.z());
0117
0118   // figure out subdetector, giving preference to HE in HE/HF overlap region
0119   HcalSubdetector bc = HcalEmpty;
0120   if (abseta <= m_topology.etaMax(HcalBarrel)) {
0121     bc = HcalBarrel;
0122   } else if (absz >= z_long) {
0123     bc = HcalForward;
0124   } else if (abseta <= m_topology.etaMax(HcalEndcap)) {
0125     bc = HcalEndcap;
0126   } else {
0127     bc = HcalForward;
0128   }
0129
0130   // find eta bin
0131   int etaring = etaRing(bc, abseta);
0132
0133   int phibin = phiBin(bc, etaring, r.phi());
0134
0135   // add a sign to the etaring
0136   int etabin = (r.z() > 0) ? etaring : -etaring;
0137
0138   if (bc == HcalForward) {
0139     static const double z_short = 1137.0;
0140     // Next line is premature depth 1 and 2 can coexist for large z-extent
0141     //    HcalDetId bestId(bc,etabin,phibin,((fabs(r.z())>=z_short)?(2):(1)));
0142     // above line is no good with finite precision
0143     HcalDetId bestId(bc, etabin, phibin, ((fabs(r.z()) - z_short > -0.1) ? (2) : (1)));
0144     return correctId(bestId);
0145   } else {
0146     //Now do depth if required
0147     int zside = (r.z() > 0) ? 1 : -1;
0148     int dbin = m_topology.dddConstants()->getMinDepth(((int)(bc)-1), etaring, phibin, zside);
0149     double pointrz(0), drz(99999.);
0150     HcalDetId currentId(bc, etabin, phibin, dbin);
0151     if (bc == HcalBarrel)
0152       pointrz = r.mag();
0153     else
0154       pointrz = std::abs(r.z());
0155     HcalDetId bestId;
0156     for (; currentId != HcalDetId(); m_topology.incrementDepth(currentId)) {
0157       std::shared_ptr<const CaloCellGeometry> cell = getGeometry(currentId);
0158       if (cell == nullptr) {
0159         assert(bestId != HcalDetId());
0160         break;
0161       } else {
0162         double rz;
0163         if (bc == HcalEndcap)
0164           rz = std::abs(cell->getPosition().z());
0165         else
0166           rz = cell->getPosition().mag();
0167         if (std::abs(pointrz - rz) < drz) {
0168           bestId = currentId;
0169           drz = std::abs(pointrz - rz);
0170         }
0171       }
0172     }
0173 #ifdef EDM_ML_DEBUG
0174     edm::LogVerbatim("HCalGeom") << bestId << " Corrected to " << HcalDetId(correctId(bestId));
0175 #endif
0176
0177     return (ignoreCorrect ? bestId : correctId(bestId));
0178   }
0179 }
0180
0181 GlobalPoint HcalGeometry::getPosition(const DetId& id) const {
0182   if (!m_mergePosition) {
0183     return (getGeometryBase(id)->getPosition());
0184   } else {
0185     return (getGeometryBase(m_topology.idFront(id))->getPosition());
0186   }
0187 }
0188
0189 GlobalPoint HcalGeometry::getPosition(uint32_t idx, bool test) const {
0190   if (test) {
0191     int subdet, z, depth, eta, phi, lay;
0192     HcalTestNumbering::unpackHcalIndex(idx, subdet, z, depth, eta, phi, lay);
0193     int sign = (z == 0) ? (-1) : (1);
0194     auto id = m_topology.dddConstants()->getHCID(subdet, (sign * eta), phi, lay, depth);
0195     auto hcId = ((id.subdet == static_cast<int>(HcalBarrel))
0196                      ? HcalDetId(HcalBarrel, sign * id.eta, id.phi, id.depth)
0197                      : ((id.subdet == static_cast<int>(HcalEndcap))
0198                             ? HcalDetId(HcalEndcap, sign * id.eta, id.phi, id.depth)
0199                             : ((id.subdet == static_cast<int>(HcalOuter))
0200                                    ? HcalDetId(HcalOuter, sign * id.eta, id.phi, id.depth)
0201                                    : ((id.subdet == static_cast<int>(HcalForward))
0202                                           ? HcalDetId(HcalForward, sign * id.eta, id.phi, id.depth)
0203                                           : HcalDetId()))));
0204     return getPosition(DetId(hcId));
0205   } else {
0206     return getPosition(DetId(idx));
0207   }
0208 }
0209
0210 GlobalPoint HcalGeometry::getBackPosition(const DetId& id) const {
0211   if (!m_mergePosition) {
0212     return (getGeometryBase(id)->getBackPoint());
0213   } else {
0214     std::vector<HcalDetId> ids;
0215     m_topology.unmergeDepthDetId(HcalDetId(id), ids);
0216     return (getGeometryBase(ids.back())->getBackPoint());
0217   }
0218 }
0219
0220 GlobalPoint HcalGeometry::getBackPosition(uint32_t idx, bool test) const {
0221   if (test) {
0222     int subdet, z, depth, eta, phi, lay;
0223     HcalTestNumbering::unpackHcalIndex(idx, subdet, z, depth, eta, phi, lay);
0224     int sign = (z == 0) ? (-1) : (1);
0225     auto id = m_topology.dddConstants()->getHCID(subdet, (sign * eta), phi, lay, depth);
0226     auto hcId = ((id.subdet == static_cast<int>(HcalBarrel))
0227                      ? HcalDetId(HcalBarrel, sign * id.eta, id.phi, id.depth)
0228                      : ((id.subdet == static_cast<int>(HcalEndcap))
0229                             ? HcalDetId(HcalEndcap, sign * id.eta, id.phi, id.depth)
0230                             : ((id.subdet == static_cast<int>(HcalOuter))
0231                                    ? HcalDetId(HcalOuter, sign * id.eta, id.phi, id.depth)
0232                                    : ((id.subdet == static_cast<int>(HcalForward))
0233                                           ? HcalDetId(HcalForward, sign * id.eta, id.phi, id.depth)
0234                                           : HcalDetId()))));
0235     return getBackPosition(DetId(hcId));
0236   } else {
0237     return getBackPosition(DetId(idx));
0238   }
0239 }
0240
0241 CaloCellGeometry::CornersVec HcalGeometry::getCorners(const DetId& id) const {
0242   if (!m_mergePosition) {
0243     return (getGeometryBase(id)->getCorners());
0244   } else {
0245     std::vector<HcalDetId> ids;
0246     m_topology.unmergeDepthDetId(HcalDetId(id), ids);
0247     CaloCellGeometry::CornersVec mcorners;
0248     CaloCellGeometry::CornersVec mcf = getGeometryBase(ids.front())->getCorners();
0249     CaloCellGeometry::CornersVec mcb = getGeometryBase(ids.back())->getCorners();
0250     for (unsigned int k = 0; k < 4; ++k) {
0251       mcorners[k] = mcf[k];
0252       mcorners[k + 4] = mcb[k + 4];
0253     }
0254     return mcorners;
0255   }
0256 }
0257
0258 int HcalGeometry::etaRing(HcalSubdetector bc, double abseta) const { return m_topology.etaRing(bc, abseta); }
0259
0260 int HcalGeometry::phiBin(HcalSubdetector bc, int etaring, double phi) const {
0261   return m_topology.phiBin(bc, etaring, phi);
0262 }
0263
0264 CaloSubdetectorGeometry::DetIdSet HcalGeometry::getCells(const GlobalPoint& r, double dR) const {
0265   CaloSubdetectorGeometry::DetIdSet dis;  // this is the return object
0266
0267   if (0.000001 < dR) {
0268     if (dR > M_PI / 2.) {                              // this version needs "small" dR
0269       dis = CaloSubdetectorGeometry::getCells(r, dR);  // base class version
0270     } else {
0271       const double dR2(dR * dR);
0272       const double reta(r.eta());
0273       const double rphi(r.phi());
0274       const double lowEta(reta - dR);
0275       const double highEta(reta + dR);
0276       const double lowPhi(rphi - dR);
0277       const double highPhi(rphi + dR);
0278
0279       const double hfEtaHi(m_topology.etaMax(HcalForward));
0280
0281       if (highEta > -hfEtaHi && lowEta < hfEtaHi) {  // in hcal
0282         const HcalSubdetector hs[] = {HcalBarrel, HcalOuter, HcalEndcap, HcalForward};
0283
0284         for (unsigned int is(0); is != 4; ++is) {
0285           const int sign(reta > 0 ? 1 : -1);
0286           const int ieta_center(sign * etaRing(hs[is], fabs(reta)));
0287           const int ieta_lo((0 < lowEta * sign ? sign : -sign) * etaRing(hs[is], fabs(lowEta)));
0288           const int ieta_hi((0 < highEta * sign ? sign : -sign) * etaRing(hs[is], fabs(highEta)));
0289           const int iphi_lo(phiBin(hs[is], ieta_center, lowPhi));
0290           const int iphi_hi(phiBin(hs[is], ieta_center, highPhi));
0291           const int jphi_lo(iphi_lo > iphi_hi ? iphi_lo - 72 : iphi_lo);
0292           const int jphi_hi(iphi_hi);
0293
0294           const int idep_lo(1 == is ? 4 : 1);
0295           const int idep_hi(m_topology.maxDepth(hs[is]));
0296           for (int ieta(ieta_lo); ieta <= ieta_hi; ++ieta) {  // over eta limits
0297             if (ieta != 0) {
0298               for (int jphi(jphi_lo); jphi <= jphi_hi; ++jphi) {  // over phi limits
0299                 const int iphi(1 > jphi ? jphi + 72 : jphi);
0300                 for (int idep(idep_lo); idep <= idep_hi; ++idep) {
0301                   const HcalDetId did(hs[is], ieta, iphi, idep);
0302                   if (m_topology.valid(did)) {
0303                     std::shared_ptr<const CaloCellGeometry> cell(getGeometryBase(did));
0304                     if (nullptr != cell) {
0305                       const GlobalPoint& p(cell->getPosition());
0306                       const double eta(p.eta());
0307                       const double phi(p.phi());
0308                       if (reco::deltaR2(eta, phi, reta, rphi) < dR2)
0309                         dis.insert(did);
0310                     }
0311                   }
0312                 }
0313               }
0314             }
0315           }
0316         }
0317       }
0318     }
0319   }
0320   return dis;
0321 }
0322
0323 unsigned int HcalGeometry::getHxSize(const int type) const {
0324   unsigned int hxsize(0);
0325   if (!m_hbIds.isSet())
0326     fillDetIds();
0327   if (type == 1)
0328     hxsize = (m_hbIds.isSet()) ? m_hbIds->size() : 0;
0329   else if (type == 2)
0330     hxsize = (m_heIds.isSet()) ? m_heIds->size() : 0;
0331   else if (type == 3)
0332     hxsize = (m_hoIds.isSet()) ? m_hoIds->size() : 0;
0333   else if (type == 4)
0334     hxsize = (m_hfIds.isSet()) ? m_hfIds->size() : 0;
0335   else if (type == 0)
0336     hxsize = (m_emptyIds.isSet()) ? m_emptyIds->size() : 0;
0337   return hxsize;
0338 }
0339
0340 DetId HcalGeometry::detIdFromBarrelAlignmentIndex(unsigned int i) {
0341   assert(i < numberOfBarrelAlignments());
0342   const int ieta(i < numberOfBarrelAlignments() / 2 ? -1 : 1);
0343   const int iphi(1 + (4 * i) % 72);
0344   return HcalDetId(HcalBarrel, ieta, iphi, 1);
0345 }
0346
0347 DetId HcalGeometry::detIdFromEndcapAlignmentIndex(unsigned int i) {
0348   assert(i < numberOfEndcapAlignments());
0349   const int ieta(i < numberOfEndcapAlignments() / 2 ? -16 : 16);
0350   const int iphi(1 + (4 * i) % 72);
0351   return HcalDetId(HcalEndcap, ieta, iphi, 1);
0352 }
0353
0354 DetId HcalGeometry::detIdFromForwardAlignmentIndex(unsigned int i) {
0355   assert(i < numberOfForwardAlignments());
0356   const int ieta(i < numberOfForwardAlignments() / 2 ? -29 : 29);
0357   const int iphi(1 + (4 * i) % 72);
0358   return HcalDetId(HcalForward, ieta, iphi, 1);
0359 }
0360
0361 DetId HcalGeometry::detIdFromOuterAlignmentIndex(unsigned int i) {
0362   assert(i < numberOfOuterAlignments());
0363   const int ring(i / 12);
0364   const int ieta(0 == ring ? -11 : 1 == ring ? -5 : 2 == ring ? 1 : 3 == ring ? 5 : 11);
0365   const int iphi(1 + (i - ring * 12) * 6);
0366   return HcalDetId(HcalOuter, ieta, iphi, 4);
0367 }
0368
0369 DetId HcalGeometry::detIdFromLocalAlignmentIndex(unsigned int i) {
0370   assert(i < numberOfAlignments());
0371
0372   const unsigned int nB(numberOfBarrelAlignments());
0373   const unsigned int nE(numberOfEndcapAlignments());
0374   const unsigned int nF(numberOfForwardAlignments());
0375   //   const unsigned int nO ( numberOfOuterAlignments()   ) ;
0376
0377   return (i < nB             ? detIdFromBarrelAlignmentIndex(i)
0378           : i < nB + nE      ? detIdFromEndcapAlignmentIndex(i - nB)
0379           : i < nB + nE + nF ? detIdFromForwardAlignmentIndex(i - nB - nE)
0380                              : detIdFromOuterAlignmentIndex(i - nB - nE - nF));
0381 }
0382
0383 unsigned int HcalGeometry::alignmentBarEndForIndexLocal(const DetId& id, unsigned int nD) {
0384   const HcalDetId hid(id);
0385   const unsigned int iphi(hid.iphi());
0386   const int ieta(hid.ieta());
0387   const unsigned int index((0 < ieta ? nD / 2 : 0) + (iphi + 1) % 72 / 4);
0388   assert(index < nD);
0389   return index;
0390 }
0391
0392 unsigned int HcalGeometry::alignmentBarrelIndexLocal(const DetId& id) {
0393   return alignmentBarEndForIndexLocal(id, numberOfBarrelAlignments());
0394 }
0395
0396 unsigned int HcalGeometry::alignmentEndcapIndexLocal(const DetId& id) {
0397   return alignmentBarEndForIndexLocal(id, numberOfEndcapAlignments());
0398 }
0399
0400 unsigned int HcalGeometry::alignmentForwardIndexLocal(const DetId& id) {
0401   return alignmentBarEndForIndexLocal(id, numberOfForwardAlignments());
0402 }
0403
0404 unsigned int HcalGeometry::alignmentOuterIndexLocal(const DetId& id) {
0405   const HcalDetId hid(id);
0406   const int ieta(hid.ieta());
0407   const int iphi(hid.iphi());
0408   const int ring(ieta < -10 ? 0 : (ieta < -4 ? 1 : (ieta < 5 ? 2 : (ieta < 11 ? 3 : 4))));
0409
0410   const unsigned int index(12 * ring + (iphi - 1) / 6);
0411   assert(index < numberOfOuterAlignments());
0412   return index;
0413 }
0414
0415 unsigned int HcalGeometry::alignmentTransformIndexLocal(const DetId& id) {
0416   assert(id.det() == DetId::Hcal);
0417
0418   const HcalDetId hid(id);
0419   bool isHB = (hid.subdet() == HcalBarrel);
0420   bool isHE = (hid.subdet() == HcalEndcap);
0421   bool isHF = (hid.subdet() == HcalForward);
0422   // bool isHO = (hid.subdet() == HcalOuter);
0423
0424   const unsigned int nB(numberOfBarrelAlignments());
0425   const unsigned int nE(numberOfEndcapAlignments());
0426   const unsigned int nF(numberOfForwardAlignments());
0427   // const unsigned int nO ( numberOfOuterAlignments()   ) ;
0428
0429   const unsigned int index(isHB   ? alignmentBarrelIndexLocal(id)
0430                            : isHE ? alignmentEndcapIndexLocal(id) + nB
0431                            : isHF ? alignmentForwardIndexLocal(id) + nB + nE
0432                                   : alignmentOuterIndexLocal(id) + nB + nE + nF);
0433
0434   assert(index < numberOfAlignments());
0435   return index;
0436 }
0437
0438 unsigned int HcalGeometry::alignmentTransformIndexGlobal(const DetId& id) { return (unsigned int)DetId::Hcal - 1; }
0439
0440 void HcalGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref) {
0441   HcalDetId hid = HcalDetId(m_topology.denseId2detId(i));
0442
0443   if (hid.subdet() == HcalForward) {
0444     IdealZPrism::localCorners(lc, pv, ref);
0445   } else {
0446     IdealObliquePrism::localCorners(lc, pv, ref);
0447   }
0448 }
0449
0450 unsigned int HcalGeometry::newCellImpl(
0451     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0452   assert(detId.det() == DetId::Hcal);
0453
0454   const HcalDetId hid(detId);
0455   unsigned int din = m_topology.detId2denseId(detId);
0456
0457 #ifdef EDM_ML_DEBUG
0458   edm::LogVerbatim("HCalGeom") << "HcalGeometry: newCell subdet " << detId.subdetId() << ", raw ID " << detId.rawId()
0459                                << ", hid " << hid << ", din " << din << ", index ";
0460 #endif
0461
0462   if (hid.subdet() == HcalBarrel) {
0463     m_hbCellVec.at(din) = IdealObliquePrism(f1, cornersMgr(), parm);
0464   } else if (hid.subdet() == HcalEndcap) {
0465     const unsigned int index(din - m_hbCellVec.size());
0466     m_heCellVec.at(index) = IdealObliquePrism(f1, cornersMgr(), parm);
0467   } else if (hid.subdet() == HcalOuter) {
0468     const unsigned int index(din - m_hbCellVec.size() - m_heCellVec.size());
0469     m_hoCellVec.at(index) = IdealObliquePrism(f1, cornersMgr(), parm);
0470   } else {
0471     const unsigned int index(din - m_hbCellVec.size() - m_heCellVec.size() - m_hoCellVec.size());
0472     m_hfCellVec.at(index) = IdealZPrism(f1, cornersMgr(), parm, hid.depth() == 1 ? IdealZPrism::EM : IdealZPrism::HADR);
0473   }
0474
0475   return din;
0476 }
0477
0478 void HcalGeometry::newCell(
0479     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0480   unsigned int din = newCellImpl(f1, f2, f3, parm, detId);
0481
0482   addValidID(detId);
0483   m_dins.emplace_back(din);
0484 }
0485
0486 void HcalGeometry::newCellFast(
0487     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0488   unsigned int din = newCellImpl(f1, f2, f3, parm, detId);
0489
0490   m_validIds.emplace_back(detId);
0491   m_dins.emplace_back(din);
0492 }
0493
0494 const CaloCellGeometry* HcalGeometry::getGeometryRawPtr(uint32_t din) const {
0495   // Modify the RawPtr class
0496   const CaloCellGeometry* cell(nullptr);
0497   if (m_hbCellVec.size() > din) {
0498     cell = (&m_hbCellVec[din]);
0499   } else if (m_hbCellVec.size() + m_heCellVec.size() > din) {
0500     const unsigned int ind(din - m_hbCellVec.size());
0501     cell = (&m_heCellVec[ind]);
0502   } else if (m_hbCellVec.size() + m_heCellVec.size() + m_hoCellVec.size() > din) {
0503     const unsigned int ind(din - m_hbCellVec.size() - m_heCellVec.size());
0504     cell = (&m_hoCellVec[ind]);
0505   } else if (m_hbCellVec.size() + m_heCellVec.size() + m_hoCellVec.size() + m_hfCellVec.size() > din) {
0506     const unsigned int ind(din - m_hbCellVec.size() - m_heCellVec.size() - m_hoCellVec.size());
0507     cell = (&m_hfCellVec[ind]);
0508   }
0509
0510   return ((nullptr == cell || nullptr == cell->param()) ? nullptr : cell);
0511 }
0512
0513 void HcalGeometry::getSummary(CaloSubdetectorGeometry::TrVec& tVec,
0514                               CaloSubdetectorGeometry::IVec& iVec,
0515                               CaloSubdetectorGeometry::DimVec& dVec,
0516                               CaloSubdetectorGeometry::IVec& dinsVec) const {
0517   tVec.reserve(m_topology.ncells() * numberOfTransformParms());
0518   iVec.reserve(numberOfShapes() == 1 ? 1 : m_topology.ncells());
0519   dVec.reserve(numberOfShapes() * numberOfParametersPerShape());
0520   dinsVec.reserve(m_topology.ncells());
0521
0522   for (const auto& pv : parVecVec()) {
0523     for (float iv : pv) {
0524       dVec.emplace_back(iv);
0525     }
0526   }
0527
0528   for (auto i : m_dins) {
0529     Tr3D tr;
0530     auto ptr = cellGeomPtr(i);
0531
0532     if (nullptr != ptr) {
0533       dinsVec.emplace_back(i);
0534
0535       const CCGFloat* par(ptr->param());
0536
0537       unsigned int ishape(9999);
0538
0539       for (unsigned int ivv(0); ivv != parVecVec().size(); ++ivv) {
0540         bool ok(true);
0541         const CCGFloat* pv(&(*parVecVec()[ivv].begin()));
0542         for (unsigned int k(0); k != numberOfParametersPerShape(); ++k) {
0543           ok = ok && (fabs(par[k] - pv[k]) < 1.e-6);
0544         }
0545         if (ok) {
0546           ishape = ivv;
0547           break;
0548         }
0549       }
0550       assert(9999 != ishape);
0551
0552       const unsigned int nn((numberOfShapes() == 1) ? (unsigned int)1 : m_dins.size());
0553       if (iVec.size() < nn)
0554         iVec.emplace_back(ishape);
0555
0556       ptr->getTransform(tr, (Pt3DVec*)nullptr);
0557
0558       if (Tr3D() == tr) {  // for preshower there is no rotation
0559         const GlobalPoint& gp(ptr->getPosition());
0560         tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
0561       }
0562
0563       const CLHEP::Hep3Vector tt(tr.getTranslation());
0564       tVec.emplace_back(tt.x());
0565       tVec.emplace_back(tt.y());
0566       tVec.emplace_back(tt.z());
0567       if (6 == numberOfTransformParms()) {
0568         const CLHEP::HepRotation rr(tr.getRotation());
0569         const ROOT::Math::Transform3D rtr(
0570             rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
0571         ROOT::Math::EulerAngles ea;
0572         rtr.GetRotation(ea);
0573         tVec.emplace_back(ea.Phi());
0574         tVec.emplace_back(ea.Theta());
0575         tVec.emplace_back(ea.Psi());
0576       }
0577     }
0578   }
0579 }
0580
0581 DetId HcalGeometry::correctId(const DetId& id) const {
0582   if (m_mergePosition) {
0583     HcalDetId hid(id);
0584     return ((DetId)(m_topology.mergedDepthDetId(hid)));
0585   } else {
0586     return id;
0587   }
0588 }
0589
0590 void HcalGeometry::increaseReserve(unsigned int extra) { m_validIds.reserve(m_validIds.size() + extra); }
0591
0592 void HcalGeometry::sortValidIds() { std::sort(m_validIds.begin(), m_validIds.end()); }