Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-16 05:03:57 

0001 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0002 #include <cmath>
0003 #include <iostream>
0004 #include <cassert>
0005 #include <algorithm>
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "DataFormats/Math/interface/GeantUnits.h"
0009 
0010 using namespace geant_units;
0011 using namespace geant_units::operators;
0012 
0013 static const int IPHI_MAX = 72;
0014 
0015 //#define EDM_ML_DEBUG
0016 
0017 HcalTopology::HcalTopology(const HcalDDDRecConstants* hcons, const bool mergePosition)
0018     : hcons_(hcons),
0019       mergePosition_(mergePosition),
0020       excludeHB_(false),
0021       excludeHE_(false),
0022       excludeHO_(false),
0023       excludeHF_(false),
0024       firstHBRing_(1),
0025       firstHERing_(999),
0026       lastHERing_(0),
0027       firstHFRing_(29),
0028       lastHFRing_(41),
0029       firstHORing_(1),
0030       lastHORing_(15),
0031       firstHEDoublePhiRing_(999),
0032       firstHEQuadPhiRing_(999),
0033       firstHFQuadPhiRing_(40),
0034       firstHETripleDepthRing_(999),
0035       singlePhiBins_(IPHI_MAX),
0036       doublePhiBins_(36),
0037       maxPhiHE_(IPHI_MAX) {
0038   mode_ = (HcalTopologyMode::Mode)(hcons_->getTopoMode());
0039   triggerMode_ = (HcalTopologyMode::TriggerMode)(hcons_->getTriggerMode());
0040   maxDepthHB_ = hcons_->getMaxDepth(0);
0041   maxDepthHE_ = hcons_->getMaxDepth(1);
0042   maxDepthHF_ = hcons_->getMaxDepth(2);
0043   etaBinsHB_ = hcons_->getEtaBins(0);
0044   etaBinsHE_ = hcons_->getEtaBins(1);
0045   nEtaHB_ = (hcons_->getEtaRange(0)).second - (hcons_->getEtaRange(0)).first + 1;
0046   lastHBRing_ = firstHBRing_ + nEtaHB_ - 1;
0047   if (hcons_->getNPhi(1) > maxPhiHE_)
0048     maxPhiHE_ = hcons_->getNPhi(1);
0049   for (auto& i : etaBinsHE_) {
0050     if (firstHERing_ > i.ieta)
0051       firstHERing_ = i.ieta;
0052     if (lastHERing_ < i.ieta)
0053       lastHERing_ = i.ieta;
0054     int unit = static_cast<int>((i.dphi / 5.0_deg) + 0.01);
0055     if (unit == 2 && firstHEDoublePhiRing_ > i.ieta)
0056       firstHEDoublePhiRing_ = i.ieta;
0057     if (unit == 4 && firstHEQuadPhiRing_ > i.ieta)
0058       firstHEQuadPhiRing_ = i.ieta;
0059     if (i.layer.size() > 2 && firstHETripleDepthRing_ > i.ieta)
0060       firstHETripleDepthRing_ = i.ieta;
0061   }
0062   if (firstHERing_ > lastHERing_) {
0063     firstHERing_ = lastHERing_ = firstHEDoublePhiRing_ = firstHEQuadPhiRing_ = firstHETripleDepthRing_ = nEtaHE_ = 0;
0064   } else {
0065     nEtaHE_ = (lastHERing_ - firstHERing_ + 1);
0066   }
0067   if (phase1()) {
0068     topoVersion_ = 0;         //DL
0069     HBSize_ = kHBSizePreLS1;  // qie-per-fiber * fiber/rm * rm/rbx * rbx/barrel * barrel/hcal
0070     HESize_ = kHESizePreLS1;  // qie-per-fiber * fiber/rm * rm/rbx * rbx/endcap * endcap/hcal
0071     HOSize_ = kHOSizePreLS1;  // ieta * iphi * 2
0072     HFSize_ = kHFSizePreLS1;  // ieta * iphi * depth * 2
0073     CALIBSize_ = kCALIBSizePreLS1;
0074     numberOfShapes_ = 87;
0075   } else {  // need to know more eventually
0076     topoVersion_ = 10;
0077     HBSize_ = nEtaHB_ * IPHI_MAX * maxDepthHB_ * 2;
0078     HESize_ = nEtaHE_ * maxPhiHE_ * maxDepthHE_ * 2;
0079     HOSize_ = (lastHORing_ - firstHORing_ + 1) * IPHI_MAX * 2;                // ieta * iphi * 2
0080     HFSize_ = (lastHFRing_ - firstHFRing_ + 1) * IPHI_MAX * maxDepthHF_ * 2;  // ieta * iphi * depth * 2
0081     CALIBSize_ = kOffCalibHFX_;
0082     numberOfShapes_ = (maxPhiHE_ > 72) ? 1200 : 500;
0083   }
0084   maxEta_ = (lastHERing_ > lastHFRing_) ? lastHERing_ : lastHFRing_;
0085   if (triggerMode_ == HcalTopologyMode::TriggerMode_2009) {
0086     HTSize_ = kHTSizePreLS1;
0087   } else {
0088     HTSize_ = kHTSizePhase1;
0089   }
0090 
0091 #ifdef EDM_ML_DEBUG
0092   edm::LogVerbatim("HCalGeom") << "Topo sizes " << HBSize_ << ":" << HESize_ << ":" << HOSize_ << ":" << HFSize_ << ":"
0093                                << HTSize_ << ":" << CALIBSize_ << " for mode " << mode_ << ":" << triggerMode_;
0094 #endif
0095 
0096   //The transition between HE/HF in eta
0097   etaTableHF = hcons_->getEtaTableHF();
0098   etaTable = hcons_->getEtaTable();
0099   dPhiTableHF = hcons_->getPhiTableHF();
0100   dPhiTable = hcons_->getPhiTable();
0101   phioff = hcons_->getPhiOffs();
0102   std::pair<int, int> ietaHF = hcons_->getEtaRange(2);
0103   etaHE2HF_ = firstHFRing_;
0104   etaHF2HE_ = lastHERing_;
0105   if (etaBinsHE_.size() > 1) {
0106     double eta = etaBinsHE_[etaBinsHE_.size() - 1].etaMax;
0107     for (unsigned int i = 1; i < etaTableHF.size(); ++i) {
0108       if (eta < etaTableHF[i]) {
0109         etaHE2HF_ = ietaHF.first + i - 1;
0110         break;
0111       }
0112     }
0113     eta = etaTableHF[0];
0114     for (auto& i : etaBinsHE_) {
0115       if (eta < i.etaMax) {
0116         etaHF2HE_ = i.ieta;
0117         break;
0118       }
0119     }
0120   }
0121   const double fiveDegInRad = 5.0_deg;
0122   for (double k : dPhiTable) {
0123     int units = (int)(k / fiveDegInRad + 0.5);
0124     unitPhi.emplace_back(units);
0125   }
0126   for (double k : dPhiTableHF) {
0127     int units = (int)(k / fiveDegInRad + 0.5);
0128     unitPhiHF.emplace_back(units);
0129   }
0130   int nEta = hcons_->getNEta();
0131   for (int ring = 1; ring <= nEta; ++ring) {
0132     std::vector<int> segmentation = hcons_->getDepth(ring - 1, false);
0133     setDepthSegmentation(ring, segmentation, false);
0134 #ifdef EDM_ML_DEBUG
0135     edm::LogVerbatim("HCalGeom") << "Set segmentation for ring " << ring << " with " << segmentation.size()
0136                                  << " elements:";
0137     for (unsigned int k = 0; k < segmentation.size(); ++k)
0138       edm::LogVerbatim("HCalGeom") << "[" << k << "] " << segmentation[k];
0139 #endif
0140     segmentation = hcons_->getDepth(ring - 1, true);
0141     setDepthSegmentation(ring, segmentation, true);
0142 #ifdef EDM_ML_DEBUG
0143     edm::LogVerbatim("HCalGeom") << "Set Plan-1 segmentation for ring " << ring << " with " << segmentation.size()
0144                                  << " elements:";
0145     for (unsigned int k = 0; k < segmentation.size(); ++k)
0146       edm::LogVerbatim("HCalGeom") << "[" << k << "] " << segmentation[k];
0147 #endif
0148   }
0149 
0150 #ifdef EDM_ML_DEBUG
0151   edm::LogVerbatim("HCalGeom") << "Constants in HcalTopology " << firstHBRing_ << ":" << lastHBRing_ << " "
0152                                << firstHERing_ << ":" << lastHERing_ << ":" << firstHEDoublePhiRing_ << ":"
0153                                << firstHEQuadPhiRing_ << ":" << firstHETripleDepthRing_ << " " << firstHFRing_ << ":"
0154                                << lastHFRing_ << ":" << firstHFQuadPhiRing_ << " " << firstHORing_ << ":" << lastHORing_
0155                                << " " << maxDepthHB_ << ":" << maxDepthHE_ << " " << nEtaHB_ << ":" << nEtaHE_ << " "
0156                                << etaHE2HF_ << ":" << etaHF2HE_ << " " << maxPhiHE_;
0157 #endif
0158 
0159   // Now redfine some of the sizes
0160 #ifdef EDM_ML_DEBUG
0161   edm::LogVerbatim("HCalGeom") << "Redefined sizes could be for HB:" << ncells(HcalBarrel)
0162                                << " HE:" << ncells(HcalEndcap) << " HF: " << ncells(HcalForward);
0163 #endif
0164 }
0165 
0166 bool HcalTopology::valid(const DetId& id) const {
0167   assert(id.det() == DetId::Hcal);
0168   return validHcal(id);
0169 }
0170 
0171 bool HcalTopology::validHcal(const HcalDetId& id) const {
0172   // check the raw rules
0173   return validRaw(id, false) && !isExcluded(id);
0174 }
0175 
0176 bool HcalTopology::validDetId(HcalSubdetector subdet, int ieta, int iphi, int depth) const {
0177   return validHcal(HcalDetId(subdet, ieta, iphi, depth));
0178 }
0179 
0180 bool HcalTopology::validHT(const HcalTrigTowerDetId& id) const {
0181   if (id.iphi() < 1 || id.iphi() > IPHI_MAX || id.ieta() == 0)
0182     return false;
0183   if (id.depth() != 0)
0184     return false;
0185   if (maxDepthHE_ == 0) {
0186     if (id.ietaAbs() > lastHBRing_ && id.ietaAbs() < firstHFRing_)
0187       return false;
0188   }
0189   // Version 2 TPs should be for HBHE when using 1TS filter scheme
0190   if (id.version() == 0 or id.version() == 2) {
0191     if (id.ietaAbs() > 28) {
0192       if (triggerMode_ >= HcalTopologyMode::TriggerMode_2017)
0193         return false;
0194       if (triggerMode_ == HcalTopologyMode::TriggerMode_2018legacy)
0195         return false;
0196       if ((id.iphi() % 4) != 1)
0197         return false;
0198       if (id.ietaAbs() > 32)
0199         return false;
0200     }
0201   } else if (id.version() == 1) {
0202     if (triggerMode_ == HcalTopologyMode::TriggerMode_2009)
0203       return false;
0204     if (id.ietaAbs() < 30 || id.ietaAbs() > 41)
0205       return false;
0206     if (id.ietaAbs() > 29 && ((id.iphi() % 2) == 0))
0207       return false;
0208     if (id.ietaAbs() > 39 && ((id.iphi() % 4) != 3))
0209       return false;
0210   } else if (id.version() > 2) {
0211     // only versions 0, 1, and 2 are supported
0212     return false;
0213   }
0214 
0215   return true;
0216 }
0217 
0218 bool HcalTopology::validCalib(const HcalCalibDetId& tid) const {
0219   bool ok(false);
0220   if (tid.calibFlavor() == HcalCalibDetId::CalibrationBox) {
0221     HcalSubdetector subdet = tid.hcalSubdet();
0222     int ieta = tid.ieta();
0223     int chan = tid.cboxChannel();
0224     unsigned int iphi = static_cast<unsigned int>(tid.iphi());
0225     if (subdet == HcalBarrel) {
0226       if ((std::find(etaCalibHB_, etaCalibHB_ + nEtaCalibHB_, ieta) != (etaCalibHB_ + nEtaCalibHB_)) &&
0227           (std::find(chanCalibHB_, chanCalibHB_ + nchanCalibHB_, chan) != (chanCalibHB_ + nchanCalibHB_)) &&
0228           (iphi >= minPhi_) && (iphi <= maxPhi_))
0229         ok = true;
0230     } else if (subdet == HcalEndcap) {
0231       if ((std::find(etaCalibHE_, etaCalibHE_ + nEtaCalibHE_, ieta) != (etaCalibHE_ + nEtaCalibHE_)) &&
0232           (std::find(chanCalibHE1_, chanCalibHE1_ + nchanCalibHE1_, chan) != (chanCalibHE1_ + nchanCalibHE1_) ||
0233            (chan == chanCalibHE2_)) &&
0234           (iphi >= minPhi_) && (iphi <= maxPhi_))
0235         ok = true;
0236     } else if (subdet == HcalForward) {
0237       if ((std::find(etaCalibHF_, etaCalibHF_ + nEtaCalibHF_, ieta) != (etaCalibHF_ + nEtaCalibHF_)) &&
0238           (std::find(chanCalibHF1_, chanCalibHF1_ + nchanCalibHF1_, chan) != (chanCalibHF1_ + nchanCalibHF1_) ||
0239            (chan == chanCalibHF2_)) &&
0240           (iphi >= minPhi_) && (iphi <= maxPhi_))
0241         ok = true;
0242     } else if (subdet == HcalOuter) {
0243       if ((std::find(etaCalibHO_, etaCalibHO_ + nEtaCalibHO_, ieta) != (etaCalibHO_ + nEtaCalibHO_)) &&
0244           (std::find(chanCalibHO_, chanCalibHO_ + nchanCalibHO_, chan) != (chanCalibHO_ + nchanCalibHO_) ||
0245            (chan == chanCalibHOs_)) &&
0246           (iphi >= minPhi_) && (iphi <= maxPhi_))
0247         ok = true;
0248     }
0249   } else if (tid.calibFlavor() == HcalCalibDetId::HOCrosstalk) {
0250     int ieta = std::abs(tid.ieta());
0251     unsigned int iphi = static_cast<unsigned int>(tid.iphi());
0252     if ((std::find(etaCalibHOX_, etaCalibHOX_ + nEtaCalibHOX_, ieta) != (etaCalibHOX_ + nEtaCalibHOX_)) &&
0253         (iphi >= minPhi_) && (iphi <= maxPhi_))
0254       ok = true;
0255   } else if (tid.calibFlavor() == HcalCalibDetId::HBX) {
0256     int ieta = std::abs(tid.ieta());
0257     unsigned int iphi = static_cast<unsigned int>(tid.iphi());
0258     if ((ieta == etaCalibHBX_) && (iphi >= minPhi_) && (iphi <= maxPhi_))
0259       ok = true;
0260   } else if (tid.calibFlavor() == HcalCalibDetId::HEX) {
0261     int ieta = std::abs(tid.ieta());
0262     unsigned int iphi = static_cast<unsigned int>(tid.iphi());
0263     if ((std::find(etaCalibHEX_, etaCalibHEX_ + nEtaCalibHEX_, ieta) != (etaCalibHEX_ + nEtaCalibHEX_)) &&
0264         (iphi >= minPhi_) && (iphi <= maxPhi_))
0265       ok = true;
0266   } else if ((tid.calibFlavor() == HcalCalibDetId::uMNqie) || (tid.calibFlavor() == HcalCalibDetId::LASERMON) ||
0267              (tid.calibFlavor() == HcalCalibDetId::CastorRadFacility)) {
0268     ok = true;
0269   }
0270   return ok;
0271 }
0272 
0273 bool HcalTopology::validHcal(const HcalDetId& id, const unsigned int flag) const {
0274   /* original logic show here because condensed form below is rather terse
0275   // check the raw rules
0276   bool ok = validHcal(id);
0277   if (flag == 0) { // This is all what is needed
0278   } else if (flag == 1) { // See if it is in the to be merged list and merged list
0279     if (hcons_->isPlan1MergedId(id))          ok = true;
0280     else if (hcons_->isPlan1ToBeMergedId(id)) ok = false;
0281   } else if (!ok) {
0282     ok = hcons_->isPlan1MergedId(id);
0283   }
0284   return ok;
0285   */
0286   return (flag > 0 and hcons_->isPlan1MergedId(id)) or
0287          ((flag != 1 or !hcons_->isPlan1ToBeMergedId(id)) and validHcal(id));
0288 }
0289 
0290 bool HcalTopology::isExcluded(const HcalDetId& id) const {
0291   bool exed = false;
0292   // first, check the full detector exclusions...  (fast)
0293   switch (id.subdet()) {
0294     case (HcalBarrel):
0295       exed = excludeHB_;
0296       break;
0297     case (HcalEndcap):
0298       exed = excludeHE_;
0299       break;
0300     case (HcalOuter):
0301       exed = excludeHO_;
0302       break;
0303     case (HcalForward):
0304       exed = excludeHF_;
0305       break;
0306     default:
0307       exed = false;
0308   }
0309   // next, check the list (slower)
0310   if (!exed && !exclusionList_.empty()) {
0311     std::vector<HcalDetId>::const_iterator i = std::lower_bound(exclusionList_.begin(), exclusionList_.end(), id);
0312     if (i != exclusionList_.end() && *i == id)
0313       exed = true;
0314   }
0315   return exed;
0316 }
0317 
0318 void HcalTopology::exclude(const HcalDetId& id) {
0319   std::vector<HcalDetId>::iterator i = std::lower_bound(exclusionList_.begin(), exclusionList_.end(), id);
0320   if (i == exclusionList_.end() || *i != id) {
0321     exclusionList_.insert(i, id);
0322   }
0323 }
0324 
0325 void HcalTopology::excludeSubdetector(HcalSubdetector subdet) {
0326   switch (subdet) {
0327     case (HcalBarrel):
0328       excludeHB_ = true;
0329       break;
0330     case (HcalEndcap):
0331       excludeHE_ = true;
0332       break;
0333     case (HcalOuter):
0334       excludeHO_ = true;
0335       break;
0336     case (HcalForward):
0337       excludeHF_ = true;
0338       break;
0339     default:
0340       break;
0341   }
0342 }
0343 
0344 std::vector<DetId> HcalTopology::east(const DetId& id) const {
0345   std::vector<DetId> vNeighborsDetId;
0346   HcalDetId neighbors[2];
0347   for (int i = 0; i < decIEta(HcalDetId(id), neighbors); i++) {
0348     if (neighbors[i].oldFormat())
0349       neighbors[i].changeForm();
0350     vNeighborsDetId.emplace_back(DetId(neighbors[i].rawId()));
0351   }
0352   return vNeighborsDetId;
0353 }
0354 
0355 std::vector<DetId> HcalTopology::west(const DetId& id) const {
0356   std::vector<DetId> vNeighborsDetId;
0357   HcalDetId neighbors[2];
0358   for (int i = 0; i < incIEta(HcalDetId(id), neighbors); i++) {
0359     if (neighbors[i].oldFormat())
0360       neighbors[i].changeForm();
0361     vNeighborsDetId.emplace_back(DetId(neighbors[i].rawId()));
0362   }
0363   return vNeighborsDetId;
0364 }
0365 
0366 std::vector<DetId> HcalTopology::north(const DetId& id) const {
0367   std::vector<DetId> vNeighborsDetId;
0368   HcalDetId neighbor;
0369   if (incIPhi(HcalDetId(id), neighbor)) {
0370     if (neighbor.oldFormat())
0371       neighbor.changeForm();
0372     vNeighborsDetId.emplace_back(DetId(neighbor.rawId()));
0373   }
0374   return vNeighborsDetId;
0375 }
0376 
0377 std::vector<DetId> HcalTopology::south(const DetId& id) const {
0378   std::vector<DetId> vNeighborsDetId;
0379   HcalDetId neighbor;
0380   if (decIPhi(HcalDetId(id), neighbor)) {
0381     if (neighbor.oldFormat())
0382       neighbor.changeForm();
0383     vNeighborsDetId.emplace_back(DetId(neighbor.rawId()));
0384   }
0385   return vNeighborsDetId;
0386 }
0387 
0388 std::vector<DetId> HcalTopology::up(const DetId& id) const {
0389   HcalDetId neighbor = id;
0390   std::vector<DetId> vNeighborsDetId;
0391   if (incrementDepth(neighbor)) {
0392     if (neighbor.oldFormat())
0393       neighbor.changeForm();
0394     vNeighborsDetId.emplace_back(neighbor);
0395   }
0396   return vNeighborsDetId;
0397 }
0398 
0399 std::vector<DetId> HcalTopology::down(const DetId& id) const {
0400   HcalDetId neighbor = id;
0401   std::vector<DetId> vNeighborsDetId;
0402   if (decrementDepth(neighbor)) {
0403     if (neighbor.oldFormat())
0404       neighbor.changeForm();
0405     vNeighborsDetId.emplace_back(neighbor);
0406   }
0407   return vNeighborsDetId;
0408 }
0409 
0410 int HcalTopology::exclude(HcalSubdetector subdet, int ieta1, int ieta2, int iphi1, int iphi2, int depth1, int depth2) {
0411   bool exed = false;
0412   // first, check the full detector exclusions...  (fast)
0413   switch (subdet) {
0414     case (HcalBarrel):
0415       exed = excludeHB_;
0416       break;
0417     case (HcalEndcap):
0418       exed = excludeHE_;
0419       break;
0420     case (HcalOuter):
0421       exed = excludeHO_;
0422       break;
0423     case (HcalForward):
0424       exed = excludeHF_;
0425       break;
0426     default:
0427       exed = false;
0428   }
0429   if (exed)
0430     return 0;  // if the whole detector is excluded...
0431 
0432   int ieta_l = std::min(ieta1, ieta2);
0433   int ieta_h = std::max(ieta1, ieta2);
0434   int iphi_l = std::min(iphi1, iphi2);
0435   int iphi_h = std::max(iphi1, iphi2);
0436   int depth_l = std::min(depth1, depth2);
0437   int depth_h = std::max(depth1, depth2);
0438 
0439   int n = 0;
0440   for (int ieta = ieta_l; ieta <= ieta_h; ieta++)
0441     for (int iphi = iphi_l; iphi <= iphi_h; iphi++)
0442       for (int depth = depth_l; depth <= depth_h; depth++) {
0443         HcalDetId id(subdet, ieta, iphi, depth);
0444         if (validRaw(id, false)) {  // use 'validRaw' to include check validity in "uncut" detector
0445           exclude(id);
0446           n++;
0447         }
0448       }
0449   return n;
0450 }
0451 
0452 /** Basic rules used to derive this code:
0453       
0454   HB has 72 towers in iphi.  Ieta 1-14 have depth=1, Ieta 15-16 have depth=1 or 2.
0455 
0456   HE ieta=16-20 have 72 towers in iphi
0457      ieta=21-29 have 36 towers in iphi
0458      ieta=16 is depth 3 only
0459      ieta=17 is depth 1 only
0460      ieta=18-26 & 29 have depth 1 and 2
0461      ieta=27-28 has depth 1-3
0462 
0463   HF ieta=29-39 have 36 in iphi
0464      ieta=40-41 have 18 in iphi (71,3,7,11...)
0465      all have two depths
0466 
0467 
0468   HO has 15 towers in ieta and 72 in iphi and depth = 4 (one value)
0469 
0470   At H2:
0471 
0472   HE ieta 17 is two depths
0473   HE ieta 22- have 36 towers in iphi (starts one higher)
0474   HE ieta 24- has three depths
0475 
0476   */
0477 
0478 bool HcalTopology::validDetIdPreLS1(const HcalDetId& id) const {
0479   const HcalSubdetector sd(id.subdet());
0480   const int ie(id.ietaAbs());
0481   const int ip(id.iphi());
0482   const int dp(id.depth());
0483 
0484   return ((ip >= 1) && (ip <= IPHI_MAX) && (dp >= 1) && (ie >= 1) &&
0485           (((sd == HcalBarrel) && (((ie <= 14) && (dp == 1)) || (((ie == 15) || (ie == 16)) && (dp <= 2)))) ||
0486            ((sd == HcalEndcap) &&
0487             (((ie == firstHERing()) && (dp == 3)) || ((ie == 17) && (dp == 1)) ||
0488              ((ie >= 18) && (ie <= 20) && (dp <= 2)) || ((ie >= 21) && (ie <= 26) && (dp <= 2) && (ip % 2 == 1)) ||
0489              ((ie >= 27) && (ie <= 28) && (dp <= 3) && (ip % 2 == 1)) || ((ie == 29) && (dp <= 2) && (ip % 2 == 1)))) ||
0490            ((sd == HcalOuter) && (ie <= 15) && (dp == 4)) ||
0491            ((sd == HcalForward) && (dp <= 2) &&
0492             (((ie >= firstHFRing()) && (ie < firstHFQuadPhiRing()) && (ip % 2 == 1)) ||
0493              ((ie >= firstHFQuadPhiRing()) && (ie <= lastHFRing()) && (ip % 4 == 3))))));
0494 }
0495 
0496 /** Is this a valid cell id? */
0497 bool HcalTopology::validRaw(const HcalDetId& id,
0498 #ifdef EDM_ML_DEBUG
0499                             const bool debug) const {
0500 #else
0501                             const bool) const {
0502 #endif
0503   bool ok = true;
0504   int ieta = id.ieta();
0505   int aieta = id.ietaAbs();
0506   int depth = id.depth();
0507   int iphi = id.iphi();
0508   int zside = id.zside();
0509   HcalSubdetector subdet = id.subdet();
0510   int maxPhi = (subdet == HcalEndcap) ? maxPhiHE_ : IPHI_MAX;
0511   if ((ieta == 0 || iphi <= 0 || iphi > maxPhi) || aieta > maxEta_)
0512     ok = false;  // outer limits
0513 
0514   if (ok) {
0515     if (subdet == HcalBarrel) {
0516       if (phase1()) {
0517         if (aieta > lastHBRing() || depth > 2 || (aieta <= 14 && depth > 1))
0518           ok = false;
0519       } else {
0520         if ((aieta > lastHBRing()) || (depth > hcons_->getMaxDepth(0, aieta, iphi, zside)) ||
0521             (depth < hcons_->getMinDepth(0, aieta, iphi, zside)))
0522           ok = false;
0523       }
0524     } else if (subdet == HcalEndcap) {
0525       if (phase1A()) {
0526         if (depth > hcons_->getMaxDepth(1, aieta, iphi, zside) || aieta < firstHERing() || aieta > lastHERing() ||
0527             (aieta == firstHERing() && depth != hcons_->getDepthEta16(2, iphi, zside)) ||
0528             (aieta == 17 && depth != 1 && mode_ != HcalTopologyMode::H2) ||  // special case at H2
0529             (((aieta >= 17 && aieta < firstHETripleDepthRing()) || aieta == lastHERing()) && depth > 2) ||
0530             (aieta >= firstHEDoublePhiRing() && (iphi % 2) == 0))
0531           ok = false;
0532       } else {
0533         if ((depth > hcons_->getMaxDepth(1, aieta, iphi, zside)) ||
0534             (depth < hcons_->getMinDepth(1, aieta, iphi, zside)) || (aieta < firstHERing()) || (aieta > lastHERing())) {
0535           ok = false;
0536         } else {
0537           for (const auto& i : etaBinsHE_) {
0538             if (aieta == i.ieta) {
0539               if (aieta >= firstHEDoublePhiRing() && (iphi % 2) == 0)
0540                 ok = false;
0541               if (aieta >= firstHEQuadPhiRing() && (iphi % 4) != 3)
0542                 ok = false;
0543               if (aieta + 1 == hcons_->getNoff(1)) {
0544                 if (depth < 1)
0545                   ok = false;
0546               } else {
0547                 if (depth < i.depthStart)
0548                   ok = false;
0549               }
0550               break;
0551             }
0552           }
0553         }
0554       }
0555     } else if (subdet == HcalOuter) {
0556       if (aieta > lastHORing() || iphi > IPHI_MAX || depth != 4)
0557         ok = false;
0558     } else if (subdet == HcalForward) {
0559       if (aieta < firstHFRing() || aieta > lastHFRing() || ((iphi % 2) == 0) ||
0560           (depth > hcons_->maxHFDepth(ieta, iphi)) || (aieta >= firstHFQuadPhiRing() && ((iphi + 1) % 4) != 0))
0561         ok = false;
0562     } else if (subdet == HcalTriggerTower) {
0563       ok = validHT(HcalTrigTowerDetId(id.rawId()));
0564     } else if (subdet == HcalOther) {
0565       ok = validCalib(HcalCalibDetId(id.rawId()));
0566     } else {
0567       ok = false;
0568     }
0569   }
0570 #ifdef EDM_ML_DEBUG
0571   if (debug)
0572     edm::LogVerbatim("HCalGeom") << "HcalValdRaw ID: " << subdet << ":" << aieta << ":" << depth << ":" << iphi << ":"
0573                                  << zside << " Mode " << mode_ << ":" << phase1() << ":" << phase1A() << ":"
0574                                  << phase1B() << ":" << phase2() << " Limits HB " << lastHBRing() << " HE "
0575                                  << firstHERing() << ":" << lastHERing() << ":" << firstHEDoublePhiRing() << ":"
0576                                  << firstHEQuadPhiRing() << " Depth " << hcons_->getMaxDepth(0, aieta, iphi, zside)
0577                                  << ":" << hcons_->getMinDepth(0, aieta, iphi, zside) << ":"
0578                                  << hcons_->getDepthEta16(2, iphi, zside) << " OK " << ok;
0579 #endif
0580   return ok;
0581 }
0582 
0583 bool HcalTopology::incIPhi(const HcalDetId& id, HcalDetId& neighbor) const {
0584   bool ok = valid(id);
0585   if (ok) {
0586     switch (id.subdet()) {
0587       case (HcalBarrel):
0588       case (HcalOuter):
0589         if (id.iphi() == IPHI_MAX)
0590           neighbor = HcalDetId(id.subdet(), id.ieta(), 1, id.depth());
0591         else
0592           neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() + 1, id.depth());
0593         break;
0594       case (HcalEndcap):
0595         if (id.ietaAbs() >= firstHEQuadPhiRing()) {
0596           if (id.iphi() == IPHI_MAX - 1)
0597             neighbor = HcalDetId(id.subdet(), id.ieta(), 3, id.depth());
0598           else
0599             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() + 4, id.depth());
0600         } else if (id.ietaAbs() >= firstHEDoublePhiRing()) {
0601           if (id.iphi() == IPHI_MAX - 1)
0602             neighbor = HcalDetId(id.subdet(), id.ieta(), 1, id.depth());
0603           else
0604             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() + 2, id.depth());
0605         } else {
0606           if (id.iphi() == maxPhiHE_)
0607             neighbor = HcalDetId(id.subdet(), id.ieta(), 1, id.depth());
0608           else
0609             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() + 1, id.depth());
0610         }
0611         break;
0612       case (HcalForward):
0613         if (id.ietaAbs() >= firstHFQuadPhiRing()) {
0614           if (id.iphi() == IPHI_MAX - 1)
0615             neighbor = HcalDetId(id.subdet(), id.ieta(), 3, id.depth());
0616           else
0617             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() + 4, id.depth());
0618         } else {
0619           if (id.iphi() == IPHI_MAX - 1)
0620             neighbor = HcalDetId(id.subdet(), id.ieta(), 1, id.depth());
0621           else
0622             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() + 2, id.depth());
0623         }
0624         if (!validRaw(neighbor))
0625           ok = false;
0626         break;
0627       default:
0628         ok = false;
0629     }
0630   }
0631   return ok;
0632 }
0633 
0634 /** Get the neighbor (if present) of the given cell with lower iphi */
0635 bool HcalTopology::decIPhi(const HcalDetId& id, HcalDetId& neighbor) const {
0636   bool ok = valid(id);
0637   if (ok) {
0638     switch (id.subdet()) {
0639       case (HcalBarrel):
0640       case (HcalOuter):
0641         if (id.iphi() == 1)
0642           neighbor = HcalDetId(id.subdet(), id.ieta(), IPHI_MAX, id.depth());
0643         else
0644           neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() - 1, id.depth());
0645         break;
0646       case (HcalEndcap):
0647         if (id.ietaAbs() >= firstHEQuadPhiRing()) {
0648           if (id.iphi() == 3)
0649             neighbor = HcalDetId(id.subdet(), id.ieta(), IPHI_MAX - 1, id.depth());
0650           else
0651             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() - 4, id.depth());
0652         } else if (id.ietaAbs() >= firstHEDoublePhiRing()) {
0653           if (id.iphi() == 1)
0654             neighbor = HcalDetId(id.subdet(), id.ieta(), IPHI_MAX - 1, id.depth());
0655           else
0656             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() - 2, id.depth());
0657         } else {
0658           if (id.iphi() == 1)
0659             neighbor = HcalDetId(id.subdet(), id.ieta(), maxPhiHE_, id.depth());
0660           else
0661             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() - 1, id.depth());
0662         }
0663         break;
0664       case (HcalForward):
0665         if (id.ietaAbs() >= firstHFQuadPhiRing()) {
0666           if (id.iphi() == 3)
0667             neighbor = HcalDetId(id.subdet(), id.ieta(), IPHI_MAX - 1, id.depth());
0668           else
0669             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() - 4, id.depth());
0670         } else {
0671           if (id.iphi() == 1)
0672             neighbor = HcalDetId(id.subdet(), id.ieta(), IPHI_MAX - 1, id.depth());
0673           else
0674             neighbor = HcalDetId(id.subdet(), id.ieta(), id.iphi() - 2, id.depth());
0675         }
0676         if (!validRaw(neighbor))
0677           ok = false;
0678         break;
0679       default:
0680         ok = false;
0681     }
0682   }
0683   return ok;
0684 }
0685 
0686 int HcalTopology::incIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
0687   if (id.zside() == 1)
0688     return incAIEta(id, neighbors);
0689   else
0690     return decAIEta(id, neighbors);
0691 }
0692 
0693 int HcalTopology::decIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
0694   if (id.zside() == 1)
0695     return decAIEta(id, neighbors);
0696   else
0697     return incAIEta(id, neighbors);
0698 }
0699 
0700 /** Increasing in |ieta|, there is always at most one neighbor */
0701 int HcalTopology::incAIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
0702   int n = 1;
0703   int aieta = id.ietaAbs();
0704 
0705   if (aieta == firstHEDoublePhiRing() - 1 && (id.iphi() % 2) == 0)
0706     neighbors[0] = HcalDetId(id.subdet(), (aieta + 1) * id.zside(), id.iphi() - 1, id.depth());
0707   else if (aieta == firstHFQuadPhiRing() - 1 && ((id.iphi() + 1) % 4) != 0)
0708     neighbors[0] =
0709         HcalDetId(id.subdet(), (aieta + 1) * id.zside(), ((id.iphi() == 1) ? (71) : (id.iphi() - 2)), id.depth());
0710   else if (aieta == firstHEQuadPhiRing() - 1 && ((id.iphi() + 1) % 4) != 0)
0711     neighbors[0] =
0712         HcalDetId(id.subdet(), (aieta + 1) * id.zside(), ((id.iphi() == 1) ? (71) : (id.iphi() - 2)), id.depth());
0713   else if (aieta == lastHBRing() && id.subdet() == HcalBarrel)
0714     neighbors[0] = HcalDetId(HcalEndcap, (aieta + 1) * id.zside(), id.iphi(), 1);
0715   else if (aieta == lastHERing() && id.subdet() == HcalEndcap)
0716     neighbors[0] = HcalDetId(HcalForward, etaHE2HF_ * id.zside(), id.iphi(), 1);
0717   else
0718     neighbors[0] = HcalDetId(id.subdet(), (aieta + 1) * id.zside(), id.iphi(), id.depth());
0719 
0720   if (!valid(neighbors[0]))
0721     n = 0;
0722   return n;
0723 }
0724 
0725 /** Decreasing in |ieta|, there are two neighbors of 40 and 21*/
0726 int HcalTopology::decAIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
0727   int n = 1;
0728   int aieta = id.ietaAbs();
0729 
0730   if (aieta == firstHEDoublePhiRing()) {
0731     n = 2;
0732     neighbors[0] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi(), id.depth());
0733     neighbors[1] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi() + 1, id.depth());
0734   } else if (aieta == firstHFQuadPhiRing()) {
0735     n = 2;
0736     neighbors[0] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi(), id.depth());
0737     if (id.iphi() == IPHI_MAX - 1)
0738       neighbors[1] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), 1, id.depth());
0739     else
0740       neighbors[1] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi() + 2, id.depth());
0741   } else if (aieta == firstHEQuadPhiRing()) {
0742     n = 2;
0743     neighbors[0] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi(), id.depth());
0744     if (id.iphi() == IPHI_MAX - 1)
0745       neighbors[1] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), 1, id.depth());
0746     else
0747       neighbors[1] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi() + 2, id.depth());
0748   } else if (aieta == 1) {
0749     neighbors[0] = HcalDetId(id.subdet(), -aieta * id.zside(), id.iphi(), id.depth());
0750   } else if (aieta == firstHERing() && id.subdet() == HcalEndcap) {
0751     neighbors[0] = HcalDetId(HcalBarrel, (aieta - 1) * id.zside(), id.iphi(), 1);
0752   } else if (aieta == firstHFRing() && id.subdet() == HcalForward) {
0753     neighbors[0] = HcalDetId(HcalEndcap, etaHF2HE_ * id.zside(), id.iphi(), 1);
0754   } else
0755     neighbors[0] = HcalDetId(id.subdet(), (aieta - 1) * id.zside(), id.iphi(), id.depth());
0756 
0757   if (!valid(neighbors[0]) && n == 2) {
0758     if (!valid(neighbors[1]))
0759       n = 0;
0760     else {
0761       n = 1;
0762       neighbors[0] = neighbors[1];
0763     }
0764   }
0765   if (n == 2 && !valid(neighbors[1]))
0766     n = 1;
0767   if (n == 1 && !valid(neighbors[0]))
0768     n = 0;
0769 
0770   return n;
0771 }
0772 
0773 void HcalTopology::depthBinInformation(
0774     HcalSubdetector subdet, int etaRing, int iphi, int zside, int& nDepthBins, int& startingBin) const {
0775   if (subdet == HcalBarrel) {
0776     if (phase1() || (mode_ == HcalTopologyMode::H2)) {
0777       if (etaRing <= 14) {
0778         nDepthBins = 1;
0779         startingBin = 1;
0780       } else {
0781         nDepthBins = 2;
0782         startingBin = 1;
0783       }
0784     } else {
0785       startingBin = hcons_->getMinDepth(0, etaRing, iphi, zside);
0786       if (etaRing == lastHBRing()) {
0787         nDepthBins = hcons_->getDepthEta16(1, iphi, zside) - startingBin + 1;
0788       } else {
0789         nDepthBins = hcons_->getMaxDepth(0, etaRing, iphi, zside) - startingBin + 1;
0790       }
0791     }
0792   } else if (subdet == HcalEndcap) {
0793     if (phase1A()) {
0794       if (etaRing == firstHERing()) {
0795         nDepthBins = 1;
0796         startingBin = 3;
0797       } else if (etaRing == 17) {
0798         nDepthBins = 1;
0799         startingBin = 1;
0800       } else if (etaRing == lastHERing()) {
0801         nDepthBins = 2;
0802         startingBin = 1;
0803       } else {
0804         nDepthBins = (etaRing >= firstHETripleDepthRing()) ? 3 : 2;
0805         startingBin = 1;
0806       }
0807     } else {
0808       if (etaRing == firstHERing()) {
0809         startingBin = hcons_->getDepthEta16(2, iphi, zside);
0810       } else {
0811         startingBin = hcons_->getMinDepth(1, etaRing, iphi, zside);
0812       }
0813       nDepthBins = hcons_->getMaxDepth(1, etaRing, iphi, zside) - startingBin + 1;
0814     }
0815   } else if (subdet == HcalForward) {
0816     nDepthBins = maxDepthHF_;
0817     startingBin = 1;
0818   } else if (subdet == HcalOuter) {
0819     nDepthBins = 1;
0820     startingBin = 4;
0821   } else {
0822     edm::LogWarning("HCalGeom") << "Bad HCAL subdetector " << subdet;
0823   }
0824 }
0825 
0826 bool HcalTopology::incrementDepth(HcalDetId& detId) const {
0827   HcalSubdetector subdet = detId.subdet();
0828   int ieta = detId.ieta();
0829   int etaRing = detId.ietaAbs();
0830   int depth = detId.depth();
0831   int iphi = detId.iphi();
0832   int zside = detId.zside();
0833   int nDepthBins(0), startingBin(0);
0834   depthBinInformation(subdet, etaRing, iphi, zside, nDepthBins, startingBin);
0835 
0836   // see if the new depth bin exists
0837   ++depth;
0838   if (depth >= (startingBin + nDepthBins)) {
0839     // handle on a case-by-case basis
0840     if (subdet == HcalBarrel && etaRing < lastHORing()) {
0841       // HO
0842       subdet = HcalOuter;
0843       depth = 4;
0844     } else if (subdet == HcalBarrel && etaRing == lastHBRing()) {
0845       // overlap
0846       if (phase2()) {  // No more chance
0847         detId = HcalDetId();
0848         return false;
0849       } else {
0850         subdet = HcalEndcap;
0851         if (phase1B())
0852           depth = hcons_->getDepthEta16(2, iphi, zside);
0853       }
0854     } else if ((subdet == HcalEndcap) && (etaRing == lastHERing() - 1) && !phase1B()) {
0855       // guard ring HF29 is behind HE 28
0856       subdet = HcalForward;
0857       (ieta > 0) ? ++ieta : --ieta;
0858       depth = 1;
0859     } else if ((subdet == HcalEndcap) && (etaRing == lastHERing()) && !phase1B()) {
0860       // split cells go to bigger granularity.  Ring 29 -> 28
0861       (ieta > 0) ? --ieta : ++ieta;
0862     } else {
0863       // no more chances
0864       detId = HcalDetId();
0865       return false;
0866     }
0867   }
0868   detId = HcalDetId(subdet, ieta, iphi, depth);
0869   return validRaw(detId);
0870 }
0871 
0872 bool HcalTopology::decrementDepth(HcalDetId& detId) const {
0873   HcalSubdetector subdet = detId.subdet();
0874   int ieta = detId.ieta();
0875   int etaRing = detId.ietaAbs();
0876   int depth = detId.depth();
0877   int iphi = detId.iphi();
0878   int zside = detId.zside();
0879   int nDepthBins, startingBin;
0880   depthBinInformation(subdet, etaRing, iphi, zside, nDepthBins, startingBin);
0881 
0882   // see if the new depth bin exists
0883   --depth;
0884   if ((subdet == HcalOuter) || (subdet == HcalEndcap && etaRing == firstHERing())) {
0885     subdet = HcalBarrel;
0886     for (int i = 0; i < nEtaHB_; ++i) {
0887       if (etaRing == etaBinsHB_[i].ieta) {
0888         depth = etaBinsHB_[i].depthStart + etaBinsHB_[i].layer.size() - 1;
0889         break;
0890       }
0891     }
0892   } else if ((subdet == HcalEndcap) && (etaRing == lastHERing()) && (depth == hcons_->getDepthEta29(iphi, zside, 0)) &&
0893              !phase1B()) {
0894     (ieta > 0) ? --ieta : ++ieta;
0895   } else if (depth <= 0) {
0896     if (subdet == HcalForward && etaRing == firstHFRing()) {
0897       // overlap
0898       subdet = HcalEndcap;
0899       etaRing = etaHF2HE_;
0900       ieta = (ieta > 0) ? etaRing : -etaRing;
0901       for (const auto& i : etaBinsHE_) {
0902         if (etaRing == i.ieta) {
0903           depth = i.depthStart + i.layer.size() - 1;
0904           break;
0905         }
0906       }
0907     } else {
0908       // no more chances
0909       detId = HcalDetId();
0910       return false;
0911     }
0912   }
0913   detId = HcalDetId(subdet, ieta, detId.iphi(), depth);
0914   return validRaw(detId);
0915 }
0916 
0917 int HcalTopology::nPhiBins(int etaRing) const {
0918   int lastPhiBin = singlePhiBins_;
0919   if (etaRing >= firstHFQuadPhiRing() || etaRing >= firstHEQuadPhiRing())
0920     lastPhiBin = doublePhiBins_ / 2;
0921   else if (etaRing >= firstHEDoublePhiRing())
0922     lastPhiBin = doublePhiBins_;
0923   if (hcons_ && etaRing >= hcons_->getEtaRange(1).first && etaRing <= hcons_->getEtaRange(1).second) {
0924     return nPhiBins(HcalBarrel, etaRing);
0925   }
0926   return lastPhiBin;
0927 }
0928 
0929 int HcalTopology::nPhiBins(HcalSubdetector bc, int etaRing) const {
0930   double phiTableVal;
0931   if (bc == HcalForward) {
0932     phiTableVal = dPhiTableHF[etaRing - firstHFRing_];
0933   } else {
0934     phiTableVal = dPhiTable[etaRing - firstHBRing_];
0935   }
0936   int lastPhiBin = 0;
0937   if (phiTableVal != 0.0)
0938     lastPhiBin = static_cast<int>((2._pi / phiTableVal) + 0.001);
0939   return lastPhiBin;
0940 }
0941 
0942 int HcalTopology::maxDepth() const {
0943   int maxd1 = std::max(maxDepthHB_, maxDepthHE_);
0944   int maxd2 = std::max(maxDepthHF_, minMaxDepth_);
0945   return std::max(maxd1, maxd2);
0946 }
0947 
0948 int HcalTopology::maxDepth(HcalSubdetector bc) const {
0949   if (bc == HcalBarrel)
0950     return maxDepthHB_;
0951   else if (bc == HcalEndcap)
0952     return maxDepthHE_;
0953   else if (bc == HcalForward)
0954     return maxDepthHF_;
0955   else
0956     return 4;
0957 }
0958 
0959 int HcalTopology::etaRing(HcalSubdetector bc, double abseta) const {
0960   int etaring = firstHBRing_;
0961   if (bc == HcalForward) {
0962     etaring = firstHFRing_;
0963     for (unsigned int k = 0; k < etaTableHF.size() - 1; ++k) {
0964       if (abseta < etaTableHF[k + 1]) {
0965         etaring += k;
0966         break;
0967       }
0968     }
0969   } else {
0970     for (unsigned int k = 0; k < etaTable.size() - 1; ++k) {
0971       if (abseta < etaTable[k + 1]) {
0972         etaring += k;
0973         break;
0974       }
0975     }
0976     if (abseta >= etaTable[etaTable.size() - 1])
0977       etaring = lastHERing_;
0978   }
0979   return etaring;
0980 }
0981 
0982 int HcalTopology::phiBin(HcalSubdetector bc, int etaring, double phi) const {
0983   //put phi in correct range (0->2pi)
0984   int index(0);
0985   if (bc == HcalBarrel) {
0986     index = (etaring - firstHBRing_);
0987     phi -= phioff[0];
0988   } else if (bc == HcalEndcap) {
0989     index = (etaring - firstHBRing_);
0990     phi -= phioff[1];
0991   } else if (bc == HcalForward) {
0992     index = (etaring - firstHFRing_);
0993     if (index < static_cast<int>(dPhiTableHF.size())) {
0994       if (index >= 0 && unitPhiHF[index] > 2)
0995         phi -= phioff[4];
0996       else
0997         phi -= phioff[2];
0998     }
0999   }
1000   if (index < 0)
1001     index = 0;
1002   if (phi < 0.0)
1003     phi += 2._pi;
1004   else if (phi > 2._pi)
1005     phi -= 2._pi;
1006   int phibin(1), unit(1);
1007   if (bc == HcalForward) {
1008     if (index < (int)(dPhiTableHF.size())) {
1009       unit = unitPhiHF[index];
1010       phibin = static_cast<int>(phi / dPhiTableHF[index]) + 1;
1011     }
1012   } else {
1013     if (index < (int)(dPhiTable.size())) {
1014       phibin = static_cast<int>(phi / dPhiTable[index]) + 1;
1015       unit = unitPhi[index];
1016     }
1017   }
1018   int iphi(phibin);
1019   if (unit == 2)
1020     iphi = (phibin - 1) * 2 + 1;
1021   else if (unit == 4)
1022     iphi = (phibin - 1) * 4 + 3;
1023   return iphi;
1024 }
1025 
1026 void HcalTopology::getDepthSegmentation(const unsigned ring, std::vector<int>& readoutDepths, const bool one) const {
1027   // if it doesn't exist, return the first entry with a lower index.  So if we only
1028   // have entries for 1 and 17, any input from 1-16 should return the entry for ring 1
1029   SegmentationMap::const_iterator pos;
1030   if (!one) {
1031     pos = depthSegmentation_.upper_bound(ring);
1032     if (pos == depthSegmentation_.begin()) {
1033       throw cms::Exception("HcalTopology") << "No depth segmentation found for ring" << ring;
1034     }
1035   } else {
1036     pos = depthSegmentationOne_.upper_bound(ring);
1037     if (pos == depthSegmentationOne_.begin()) {
1038       throw cms::Exception("HcalTopology") << "No depth segmentation found for ring" << ring;
1039     }
1040   }
1041   --pos;
1042   // pos now refers to the last element with key <= ring.
1043   readoutDepths = pos->second;
1044 }
1045 
1046 void HcalTopology::setDepthSegmentation(const unsigned ring, const std::vector<int>& readoutDepths, const bool one) {
1047   if (one) {
1048     depthSegmentationOne_[ring] = readoutDepths;
1049   } else {
1050     depthSegmentation_[ring] = readoutDepths;
1051   }
1052 }
1053 
1054 std::pair<int, int> HcalTopology::segmentBoundaries(const unsigned ring, const unsigned depth, const bool one) const {
1055   std::vector<int> readoutDepths;
1056   getDepthSegmentation(ring, readoutDepths, one);
1057   int d1 = std::lower_bound(readoutDepths.begin(), readoutDepths.end(), depth) - readoutDepths.begin();
1058   int d2 = std::upper_bound(readoutDepths.begin(), readoutDepths.end(), depth) - readoutDepths.begin();
1059   return std::pair<int, int>(d1, d2);
1060 }
1061 
1062 double HcalTopology::etaMax(HcalSubdetector subdet) const {
1063   double eta(0);
1064   switch (subdet) {
1065     case (HcalBarrel):
1066       if (lastHBRing_ < (int)(etaTable.size()))
1067         eta = etaTable[lastHBRing_];
1068       break;
1069     case (HcalEndcap):
1070       if (lastHERing_ < (int)(etaTable.size()) && nEtaHE_ > 0)
1071         eta = etaTable[lastHERing_];
1072       break;
1073     case (HcalOuter):
1074       if (lastHORing_ < (int)(etaTable.size()))
1075         eta = etaTable[lastHORing_];
1076       break;
1077     case (HcalForward):
1078       if (!etaTableHF.empty())
1079         eta = etaTableHF[etaTableHF.size() - 1];
1080       break;
1081     default:
1082       eta = 0;
1083   }
1084   return eta;
1085 }
1086 
1087 std::pair<double, double> HcalTopology::etaRange(HcalSubdetector subdet, int keta) const {
1088   int ieta = (keta > 0) ? keta : -keta;
1089   if (subdet == HcalForward) {
1090     if (ieta >= firstHFRing_) {
1091       unsigned int ii = static_cast<unsigned int>(ieta - firstHFRing_);
1092       if (ii + 1 < etaTableHF.size())
1093         return std::pair<double, double>(etaTableHF[ii], etaTableHF[ii + 1]);
1094     }
1095   } else {
1096     int ietal = (phase1() && ieta == lastHERing_ - 1) ? (ieta + 1) : ieta;
1097     if ((ietal < static_cast<int>(etaTable.size())) && (ieta > 0))
1098       return std::pair<double, double>(etaTable[ieta - 1], etaTable[ietal]);
1099   }
1100   return std::pair<double, double>(0, 0);
1101 }
1102 
1103 unsigned int HcalTopology::detId2denseIdPreLS1(const DetId& id) const {
1104   HcalDetId hid(id);
1105   const HcalSubdetector sd(hid.subdet());
1106   const int ip(hid.iphi());
1107   const int ie(hid.ietaAbs());
1108   const int dp(hid.depth());
1109   const int zn(hid.zside() < 0 ? 1 : 0);
1110   unsigned int retval =
1111       ((sd == HcalBarrel)
1112            ? (ip - 1) * 18 + dp - 1 + ie - (ie < 16 ? 1 : 0) + zn * kHBhalf
1113            : ((sd == HcalEndcap)
1114                   ? 2 * kHBhalf + (ip - 1) * 8 + (ip / 2) * 20 +
1115                         ((ie == 16 || ie == 17)
1116                              ? ie - 16
1117                              : ((ie >= 18 && ie <= 20)
1118                                     ? 2 + 2 * (ie - 18) + dp - 1
1119                                     : ((ie >= 21 && ie <= 26)
1120                                            ? 8 + 2 * (ie - 21) + dp - 1
1121                                            : ((ie >= 27 && ie <= 28) ? 20 + 3 * (ie - 27) + dp - 1
1122                                                                      : 26 + 2 * (ie - 29) + dp - 1)))) +
1123                         zn * kHEhalf
1124                   : ((sd == HcalOuter)
1125                          ? 2 * kHBhalf + 2 * kHEhalf + (ip - 1) * 15 + (ie - 1) + zn * kHOhalf
1126                          : ((sd == HcalForward) ? 2 * kHBhalf + 2 * kHEhalf + 2 * kHOhalf + ((ip - 1) / 4) * 4 +
1127                                                       ((ip - 1) / 2) * 22 + 2 * (ie - 29) + (dp - 1) + zn * kHFhalf
1128                                                 : 0xFFFFFFFFu))));
1129   return retval;
1130 }
1131 
1132 unsigned int HcalTopology::detId2denseIdHB(const DetId& id) const {
1133   HcalDetId hid(id);
1134   const int ip(hid.iphi());
1135   const int ie(hid.ietaAbs());
1136   const int dp(hid.depth());
1137   const int zn(hid.zside() < 0 ? 1 : 0);
1138   unsigned int retval = 0xFFFFFFFFu;
1139   if (topoVersion_ == 0) {
1140     retval = (ip - 1) * 18 + dp - 1 + ie - (ie < 16 ? 1 : 0) + zn * kHBhalf;
1141   } else if (topoVersion_ == 10) {
1142     retval = (dp - 1) + maxDepthHB_ * (ip - 1);
1143     if (hid.ieta() > 0)
1144       retval += maxDepthHB_ * IPHI_MAX * (hid.ieta() - firstHBRing());
1145     else
1146       retval += maxDepthHB_ * IPHI_MAX * (hid.ieta() + lastHBRing() + nEtaHB_);
1147   }
1148   return retval;
1149 }
1150 
1151 unsigned int HcalTopology::detId2denseIdHE(const DetId& id) const {
1152   HcalDetId hid(id);
1153   const int ip(hid.iphi());
1154   const int ie(hid.ietaAbs());
1155   const int dp(hid.depth());
1156   const int zn(hid.zside() < 0 ? 1 : 0);
1157   unsigned int retval = 0xFFFFFFFFu;
1158   if (topoVersion_ == 0) {
1159     retval = (ip - 1) * 8 + (ip / 2) * 20 +
1160              ((ie == 16 || ie == 17)
1161                   ? ie - 16
1162                   : ((ie >= 18 && ie <= 20)
1163                          ? 2 + 2 * (ie - 18) + dp - 1
1164                          : ((ie >= 21 && ie <= 26) ? 8 + 2 * (ie - 21) + dp - 1
1165                                                    : ((ie >= 27 && ie <= 28) ? 20 + 3 * (ie - 27) + dp - 1
1166                                                                              : 26 + 2 * (ie - 29) + dp - 1)))) +
1167              zn * kHEhalf;
1168   } else if (topoVersion_ == 10) {
1169     retval = (dp - 1) + maxDepthHE_ * (ip - 1);
1170     if (hid.ieta() > 0)
1171       retval += maxDepthHE_ * maxPhiHE_ * (hid.ieta() - firstHERing());
1172     else
1173       retval += maxDepthHE_ * maxPhiHE_ * (hid.ieta() + lastHERing() + nEtaHE_);
1174   }
1175   return retval;
1176 }
1177 
1178 unsigned int HcalTopology::detId2denseIdHO(const DetId& id) const {
1179   HcalDetId hid(id);
1180   const int ip(hid.iphi());
1181   const int ie(hid.ietaAbs());
1182   const int zn(hid.zside() < 0 ? 1 : 0);
1183 
1184   unsigned int retval = 0xFFFFFFFFu;
1185   if (topoVersion_ == 0) {
1186     retval = (ip - 1) * 15 + (ie - 1) + zn * kHOhalf;
1187   } else if (topoVersion_ == 10) {
1188     if (hid.ieta() > 0)
1189       retval = (ip - 1) + IPHI_MAX * (hid.ieta() - 1);
1190     else
1191       retval = (ip - 1) + IPHI_MAX * (30 + hid.ieta());
1192   }
1193   return retval;
1194 }
1195 
1196 unsigned int HcalTopology::detId2denseIdHF(const DetId& id) const {
1197   HcalDetId hid(id);
1198   const int ip(hid.iphi());
1199   const int ie(hid.ietaAbs());
1200   const int dp(hid.depth());
1201   const int zn(hid.zside() < 0 ? 1 : 0);
1202 
1203   unsigned int retval = 0xFFFFFFFFu;
1204   if (topoVersion_ == 0) {
1205     retval = ((ip - 1) / 4) * 4 + ((ip - 1) / 2) * 22 + 2 * (ie - 29) + (dp - 1) + zn * kHFhalf;
1206   } else if (topoVersion_ == 10) {
1207     retval = dp - 1 + 2 * (ip - 1);
1208     if (hid.ieta() > 0)
1209       retval += maxDepthHF_ * IPHI_MAX * (hid.ieta() - 29);
1210     else
1211       retval += maxDepthHF_ * IPHI_MAX * ((41 + 13) + hid.ieta());
1212   }
1213   return retval;
1214 }
1215 
1216 unsigned int HcalTopology::detId2denseIdHT(const DetId& id) const {
1217   HcalTrigTowerDetId tid(id);
1218   int zside = tid.zside();
1219   unsigned int ietaAbs = tid.ietaAbs();
1220   unsigned int iphi = tid.iphi();
1221   unsigned int ivers = tid.version();
1222 
1223   unsigned int index;
1224   if (ivers == 0) {
1225     if ((iphi - 1) % 4 == 0)
1226       index = (iphi - 1) * 32 + (ietaAbs - 1) - (12 * ((iphi - 1) / 4));
1227     else
1228       index = (iphi - 1) * 28 + (ietaAbs - 1) + (4 * (((iphi - 1) / 4) + 1));
1229     if (zside == -1)
1230       index += kHThalf;
1231   } else {
1232     index = kHTSizePreLS1;
1233     if (zside == -1)
1234       index += ((kHTSizePhase1 - kHTSizePreLS1) / 2);
1235     index += (36 * (ietaAbs - 30) + ((iphi - 1) / 2));
1236   }
1237 
1238   return index;
1239 }
1240 
1241 unsigned int HcalTopology::detId2denseIdCALIB(const DetId& id) const {
1242   HcalCalibDetId tid(id);
1243   int channel = tid.cboxChannel();
1244   int ieta = tid.ieta();
1245   int iphi = tid.iphi();
1246   int zside = tid.zside();
1247   unsigned int index = 0xFFFFFFFFu;
1248 
1249   if (tid.calibFlavor() == HcalCalibDetId::CalibrationBox) {
1250     HcalSubdetector subDet = tid.hcalSubdet();
1251 
1252     if (subDet == HcalBarrel) {
1253 #ifdef EDM_ML_DEBUG
1254       edm::LogVerbatim("HCalGeom") << "CALIB_HB:  dphi = 4 (18 phi values), 3 channel types (0,1,2), eta = -1 or 1\n   "
1255                                       "        total of 18*3*2=108 channels";
1256 #endif
1257       auto indx = std::find(etaCalibHB_, etaCalibHB_ + nEtaCalibHB_, ieta);
1258       auto kndx = std::find(chanCalibHB_, chanCalibHB_ + nchanCalibHB_, channel);
1259       if (indx != etaCalibHB_ + nEtaCalibHB_ && kndx != chanCalibHB_ + nchanCalibHB_) {
1260         int keta = static_cast<int>(indx - etaCalibHB_);
1261         int kchn = static_cast<int>(kndx - chanCalibHB_);
1262         index = ((iphi + 1) / mPhiCalibHB_ - 1) + kPhiCalibHB_ * kchn + keta * kchanCalibHB_ + kOffCalibHB_;
1263       }
1264     } else if (subDet == HcalEndcap) {
1265 #ifdef EDM_ML_DEBUG
1266       edm::LogVerbatim("HCalGeom") << "CALIB_HE:  dphi = 4 (18 phi values), 7 channel types (0,1,2,3,4,5,6), eta = "
1267                                       "-1/+1\n           total of 18*7*2=252 channels    if (channel > 2) channel -= 1";
1268 #endif
1269       auto indx = std::find(etaCalibHE_, etaCalibHE_ + nEtaCalibHE_, ieta);
1270       if (indx != etaCalibHE_ + nEtaCalibHE_) {
1271         int keta = static_cast<int>(indx - etaCalibHE_);
1272         auto kndx = std::find(chanCalibHE1_, chanCalibHE1_ + nchanCalibHE1_, channel);
1273         if (kndx != chanCalibHE1_ + nchanCalibHE1_) {
1274           int kchn = static_cast<int>(kndx - chanCalibHE1_);
1275           index = ((iphi + 1) / mPhiCalibHE_ - 1) + kPhiCalibHE_ * kchn + keta * kchanCalibHE1_ + kOffCalibHE1_;
1276         } else if (channel == chanCalibHE2_) {
1277           index = ((iphi + 1) / mPhiCalibHE_ - 1) + keta * kchanCalibHE2_ + kOffCalibHE2_;
1278         }
1279       }
1280     } else if (subDet == HcalForward) {
1281 #ifdef EDM_ML_DEBUG
1282       edm::LogVerbatim("HCalGeom")
1283           << "CALIB_HF:  dphi = 18 (4 phi values), 3 channel types (0,1,8), eta = -1 or +1\n           or iphi = 1; "
1284              "channel = 9, eta = -1 or +1      total channels 4*3*2=24 + 2";
1285 #endif
1286       auto indx = std::find(etaCalibHF_, etaCalibHF_ + nEtaCalibHF_, ieta);
1287       if (indx != etaCalibHF_ + nEtaCalibHF_) {
1288         int keta = static_cast<int>(indx - etaCalibHF_);
1289         auto kndx = std::find(chanCalibHF1_, chanCalibHF1_ + nchanCalibHF1_, channel);
1290         if (kndx != chanCalibHF1_ + nchanCalibHF1_) {
1291           int kchn = static_cast<int>(kndx - chanCalibHF1_);
1292           index = ((iphi - 1) / mPhiCalibHF1_) + kPhiCalibHF1_ * kchn + keta * kchanCalibHF1_ + kOffCalibHF1_;
1293         } else if (channel == chanCalibHF2_) {
1294           index = keta * kchanCalibHF2_ + kOffCalibHF2_;
1295         }
1296       }
1297     } else if (subDet == HcalOuter) {
1298 #ifdef EDM_ML_DEBUG
1299       edm::LogVerbatim("HCalGeom") << "CALIB_HO:  ";
1300 #endif
1301       auto indx = std::find(etaCalibHO_, etaCalibHO_ + nEtaCalibHO_, ieta);
1302       if (indx != etaCalibHO_ + nEtaCalibHO_) {
1303         int keta = static_cast<int>(indx - etaCalibHO_);
1304         auto kndx = std::find(chanCalibHO_, chanCalibHO_ + nchanCalibHO_, channel);
1305         //there are 5 special calib crosstalk channels, one in each ring
1306         if (channel == chanCalibHOs_) {
1307           index = keta + kOffCalibHO2_;
1308         } else if (kndx != chanCalibHO_ + nchanCalibHO_) {
1309           //for HOM/HOP dphi = 6 (12 phi values),  2 channel types (0,1), eta = -2,-1 or 1,2
1310           //for HO0/YB0 dphi = 12 (6 phi values),  2 channel types (0,1), eta = 0
1311           int kchn = static_cast<int>(kndx - chanCalibHO_);
1312           int kphi = (ieta == 0) ? ((iphi + 1) / mPhiCalibHO0_ - 1) : ((iphi + 1) / mPhiCalibHO1_ - 1);
1313           if (ieta < 0) {
1314             index = kphi + kPhiCalibHO2_ * kchn + kPhiCalibHO1_ * keta + kOffCalibHO1_;
1315           } else if (ieta > 0) {
1316             index = kphi + kPhiCalibHO2_ * kchn + kPhiCalibHO1_ * (keta - 1) + kPhiCalibHO0_ + kOffCalibHO1_;
1317           } else {
1318             index = kphi + kPhiCalibHO2_ * kchn + kPhiCalibHO1_ * keta + kOffCalibHO1_;
1319           }
1320         }
1321       }
1322     } else {
1323       edm::LogWarning("HGCalGeom") << "HCAL Det Id not valid!";
1324     }
1325   } else if (tid.calibFlavor() == HcalCalibDetId::HOCrosstalk) {
1326 #ifdef EDM_ML_DEBUG
1327     edm::LogVerbatim("HCalGeom") << "HX:  for YB0/HO0 phi is grouped in 6 groups of 6 with dphi=2 but the transitions "
1328                                     "are 1 or 3 in such a way that the %36 operation yeilds unique values for every "
1329                                     "iphi\n     ieta = 0 for HO2M/HO1M ieta=2 for HO1P/HO2P; /ieta = 1 YB0/HO0";
1330 #endif
1331     int kphi = ((std::abs(ieta) == etaCalibHOX_[0]) ? ((iphi % 2 == 0) ? (iphi / 2 - 1) : (iphi - 1) / 2) : (iphi - 1));
1332     if (std::abs(ieta) == etaCalibHOX_[0]) {
1333       index = kphi + (((zside + 1) * nPhiCalibHOX_[0]) / 2) + nPhiCalibHOX_[1] + kOffCalibHOX_;
1334     } else if (std::abs(ieta) == etaCalibHOX_[1]) {
1335       index = kphi + ((zside + 1) * nPhiCalibHOX_[1]) + kOffCalibHOX_;
1336     }
1337   } else if (tid.calibFlavor() == HcalCalibDetId::HBX) {
1338     if (std::abs(ieta) == etaCalibHBX_) {
1339       index = kOffCalibHBX_ + (iphi - 1) + (zside + 1) * kPhiCalibHBX_ / 2;
1340     }
1341   } else if (tid.calibFlavor() == HcalCalibDetId::HEX) {
1342     if (std::abs(ieta) == etaCalibHEX_[0]) {
1343       index = kOffCalibHEX_ + (iphi - 1) / mPhiCalibHEX_ + (zside + 1) * kPhiCalibHEX_ / 2;
1344     } else if (std::abs(ieta) == etaCalibHEX_[1]) {
1345       index = kOffCalibHEX_ + (iphi - 1) / mPhiCalibHEX_ + (zside + 1) * kPhiCalibHEX_ / 2 + 2 * kPhiCalibHEX_;
1346     }
1347   }
1348 #ifdef EDM_ML_DEBUG
1349   edm::LogVerbatim("HCalGeom") << "  " << ieta << "  " << zside << "  " << iphi << "  " << index;
1350 #endif
1351   return index;
1352 }
1353 
1354 HcalCalibDetId HcalTopology::denseId2detIdCALIB(const unsigned int& hid) const {
1355   HcalCalibDetId id;
1356   unsigned int hid0(hid);
1357   if (hid0 < kOffCalibHOX_) {
1358     HcalSubdetector subdet(HcalEmpty);
1359     int ieta(0), iphi(0), ichan(0), ctype(0);
1360     if (hid0 < kOffCalibHE1_) {
1361       int id0 = static_cast<int>(hid0);
1362       subdet = HcalBarrel;
1363       iphi = hid0 % kPhiCalibHB_;
1364       int keta = (hid0 < kchanCalibHB_) ? 0 : 1;
1365       ieta = etaCalibHB_[keta];
1366       ichan = (id0 - iphi - keta * kchanCalibHB_) / kPhiCalibHB_;
1367       iphi = mPhiCalibHB_ * (iphi + 1) - 1;
1368       ctype = chanCalibHB_[ichan];
1369     } else if (hid0 < kOffCalibHF1_) {
1370       hid0 -= kOffCalibHE1_;
1371       int id0 = static_cast<int>(hid0);
1372       subdet = HcalEndcap;
1373       iphi = hid0 % kPhiCalibHE_;
1374       int keta = (hid0 < kchanCalibHE1_) ? 0 : 1;
1375       ieta = etaCalibHE_[keta];
1376       ichan = (id0 - iphi - keta * kchanCalibHE1_) / kPhiCalibHE_;
1377       iphi = mPhiCalibHE_ * (iphi + 1) - 1;
1378       ctype = chanCalibHE1_[ichan];
1379     } else if (hid0 < kOffCalibHO1_) {
1380       hid0 -= kOffCalibHF1_;
1381       int id0 = static_cast<int>(hid0);
1382       subdet = HcalForward;
1383       iphi = hid0 % kPhiCalibHF1_;
1384       int keta = (hid0 < kchanCalibHF1_) ? 0 : 1;
1385       ieta = etaCalibHF_[keta];
1386       ichan = (id0 - iphi - keta * kchanCalibHF1_) / kPhiCalibHF1_;
1387       iphi = mPhiCalibHF1_ * iphi + 1;
1388       ctype = chanCalibHF1_[ichan];
1389     } else if (hid0 < kOffCalibHO2_) {
1390       hid0 -= kOffCalibHO1_;
1391       int id0 = static_cast<int>(hid0);
1392       subdet = HcalOuter;
1393       unsigned int kphi = hid0 % kPhiCalibHO2_;
1394       if (kphi < 2 * kPhiCalibHO1_) {
1395         ieta = (kphi >= kPhiCalibHO1_) ? etaCalibHO_[1] : etaCalibHO_[0];
1396         iphi = kphi % kPhiCalibHO1_;
1397         ichan = (id0 - iphi - (ieta + 2) * kPhiCalibHO1_) / kPhiCalibHO2_;
1398         iphi = (iphi + 1) * mPhiCalibHO1_ - 1;
1399       } else if (kphi < (2 * kPhiCalibHO1_ + kPhiCalibHO0_)) {
1400         ieta = etaCalibHO_[2];
1401         iphi = kphi % kPhiCalibHO0_;
1402         ichan = (id0 - iphi - (ieta + 2) * kPhiCalibHO1_) / kPhiCalibHO2_;
1403         iphi = (iphi + 1) * mPhiCalibHO0_ - 1;
1404       } else {
1405         ieta = (kphi >= 3 * kPhiCalibHO1_ + kPhiCalibHO0_) ? etaCalibHO_[4] : etaCalibHO_[3];
1406         iphi = kphi % kPhiCalibHO1_;
1407         ichan = (id0 - iphi - (ieta + 1) * kPhiCalibHO1_ - kPhiCalibHO0_) / kPhiCalibHO2_;
1408         iphi = (iphi + 1) * mPhiCalibHO1_ - 1;
1409       }
1410       ctype = chanCalibHO_[ichan];
1411     } else if (hid0 < kOffCalibHE2_) {
1412       hid0 -= kOffCalibHO2_;
1413       subdet = HcalOuter;
1414       iphi = phiCalibHO_[hid0];
1415       ctype = static_cast<int>(chanCalibHOs_);
1416       ieta = etaCalibHO_[hid0];
1417     } else if (hid0 < kOffCalibHF2_) {
1418       hid0 -= kOffCalibHE2_;
1419       subdet = HcalEndcap;
1420       iphi = hid0 % kPhiCalibHE_;
1421       int keta = (hid0 < kchanCalibHE2_) ? 0 : 1;
1422       ieta = etaCalibHE_[keta];
1423       iphi = mPhiCalibHE_ * (iphi + 1) - 1;
1424       ctype = chanCalibHE2_;
1425     } else {
1426       hid0 -= kOffCalibHF2_;
1427       subdet = HcalForward;
1428       int keta = (hid0 < kchanCalibHF2_) ? 0 : 1;
1429       ieta = etaCalibHF_[keta];
1430       iphi = phiCalibHF2_;
1431       ctype = chanCalibHF2_;
1432     }
1433     id = HcalCalibDetId(subdet, ieta, iphi, ctype);
1434 #ifdef EDM_ML_DEBUG
1435     edm::LogVerbatim("HCalGeom") << "CalibrationBox: " << hid << " o/p " << ieta << ":" << iphi << ":" << ichan << ":"
1436                                  << ctype << " " << id;
1437 #endif
1438   } else if (hid < kOffCalibHBX_) {
1439     hid0 -= kOffCalibHOX_;
1440     int iphi, ieta;
1441     if (hid0 < nPhiCalibHOX_[1]) {
1442       iphi = static_cast<int>(hid0) + 1;
1443       ieta = -etaCalibHOX_[1];
1444     } else if (hid0 < (nPhiCalibHOX_[1] + nPhiCalibHOX_[0])) {
1445       hid0 -= nPhiCalibHOX_[1];
1446       iphi = ((hid0 + phiCalibHOX1_) % phiCalibHOX2_ < phiCalibHOX3_) ? 2 * hid0 + 1 : 2 * hid0 + 2;
1447       ieta = -etaCalibHOX_[0];
1448     } else if (hid0 < (nPhiCalibHOX_[1] + 2 * nPhiCalibHOX_[0])) {
1449       hid0 -= (nPhiCalibHOX_[1] + nPhiCalibHOX_[0]);
1450       iphi = ((hid0 + phiCalibHOX1_) % phiCalibHOX2_ < phiCalibHOX3_) ? 2 * hid0 + 1 : 2 * hid0 + 2;
1451       ieta = etaCalibHOX_[0];
1452     } else {
1453       hid0 -= (nPhiCalibHOX_[1] + 2 * nPhiCalibHOX_[0]);
1454       iphi = static_cast<int>(hid0) + 1;
1455       ieta = etaCalibHOX_[1];
1456     }
1457     id = HcalCalibDetId(HcalCalibDetId::HOCrosstalk, ieta, iphi);
1458 #ifdef EDM_ML_DEBUG
1459     edm::LogVerbatim("HCalGeom") << "HOCrossTalk: " << hid << ":" << hid0 << " o/p " << ieta << ":" << iphi << " "
1460                                  << id;
1461 #endif
1462   } else if (hid < kOffCalibHEX_) {
1463     hid0 -= kOffCalibHBX_;
1464     int ieta = (hid0 >= kPhiCalibHBX_) ? etaCalibHBX_ : -etaCalibHBX_;
1465     int iphi = (hid0 % kPhiCalibHBX_) + 1;
1466     id = HcalCalibDetId(HcalCalibDetId::HBX, ieta, iphi);
1467 #ifdef EDM_ML_DEBUG
1468     edm::LogVerbatim("HCalGeom") << "HBX: " << hid << ":" << hid0 << " o/p " << ieta << ":" << iphi << " " << id;
1469 #endif
1470   } else if (hid < kOffCalibHFX_) {
1471     hid0 -= kOffCalibHEX_;
1472     int iphi = 2 * (hid0 % kPhiCalibHEX_) + 1;
1473     int ieta = ((hid0 < kPhiCalibHEX_)
1474                     ? -etaCalibHEX_[0]
1475                     : ((hid0 < 2 * kPhiCalibHEX_) ? etaCalibHEX_[0]
1476                                                   : ((hid0 < 3 * kPhiCalibHEX_) ? -etaCalibHEX_[1] : etaCalibHEX_[1])));
1477     id = HcalCalibDetId(HcalCalibDetId::HEX, ieta, iphi);
1478 #ifdef EDM_ML_DEBUG
1479     edm::LogVerbatim("HCalGeom") << "HEX: " << hid << ":" << hid0 << " o/p " << ieta << ":" << iphi << " " << id;
1480 #endif
1481   }
1482   return id;
1483 }
1484 
1485 unsigned int HcalTopology::detId2denseId(const DetId& id) const {
1486   unsigned int retval(0);
1487   if (topoVersion_ == 0) {  // pre-LS1
1488     retval = detId2denseIdPreLS1(id);
1489   } else if (topoVersion_ == 10) {
1490     HcalDetId hid(id);
1491     if (hid.subdet() == HcalBarrel) {
1492       retval = (hid.depth() - 1) + maxDepthHB_ * (hid.iphi() - 1);
1493       if (hid.ieta() > 0)
1494         retval += maxDepthHB_ * IPHI_MAX * (hid.ieta() - firstHBRing());
1495       else
1496         retval += maxDepthHB_ * IPHI_MAX * (hid.ieta() + lastHBRing() + nEtaHB_);
1497     } else if (hid.subdet() == HcalEndcap) {
1498       retval = HBSize_;
1499       retval += (hid.depth() - 1) + maxDepthHE_ * (hid.iphi() - 1);
1500       if (hid.ieta() > 0)
1501         retval += maxDepthHE_ * maxPhiHE_ * (hid.ieta() - firstHERing());
1502       else
1503         retval += maxDepthHE_ * maxPhiHE_ * (hid.ieta() + lastHERing() + nEtaHE_);
1504     } else if (hid.subdet() == HcalOuter) {
1505       retval = HBSize_ + HESize_;
1506       if (hid.ieta() > 0)
1507         retval += (hid.iphi() - 1) + IPHI_MAX * (hid.ieta() - 1);
1508       else
1509         retval += (hid.iphi() - 1) + IPHI_MAX * (30 + hid.ieta());
1510     } else if (hid.subdet() == HcalForward) {
1511       retval = HBSize_ + HESize_ + HOSize_;
1512       retval += (hid.depth() - 1) + maxDepthHF_ * (hid.iphi() - 1);
1513       if (hid.ieta() > 0)
1514         retval += maxDepthHF_ * IPHI_MAX * (hid.ieta() - 29);
1515       else
1516         retval += maxDepthHF_ * IPHI_MAX * ((41 + 13) + hid.ieta());
1517     } else {
1518       retval = 0xFFFFFFFu;
1519     }
1520   }
1521 #ifdef EDM_ML_DEBUG
1522   edm::LogVerbatim("HcalGeom") << "DetId2Dense " << topoVersion_ << " ID " << std::hex << id.rawId() << std::dec
1523                                << " | " << HcalDetId(id) << " : " << std::hex << retval << std::dec;
1524 #endif
1525   return retval;
1526 }
1527 
1528 DetId HcalTopology::denseId2detId(unsigned int denseid) const {
1529   HcalSubdetector sd(HcalBarrel);
1530   int ie(0);
1531   int ip(0);
1532   int dp(0);
1533   int in(denseid);
1534   int iz(1);
1535   if (topoVersion_ == 0) {  //DL// pre-LS1
1536     if (denseid < kSizeForDenseIndexingPreLS1) {
1537       if (in > 2 * (kHBhalf + kHEhalf + kHOhalf) - 1) {  // HF
1538         sd = HcalForward;
1539         in -= 2 * (kHBhalf + kHEhalf + kHOhalf);
1540         iz = (in < kHFhalf ? 1 : -1);
1541         in %= kHFhalf;
1542         ip = 4 * (in / 48);
1543         in %= 48;
1544         ip += 1 + (in > 21 ? 2 : 0);
1545         if (3 == ip % 4)
1546           in -= 22;
1547         ie = 29 + in / 2;
1548         dp = 1 + in % 2;
1549       } else if (in > 2 * (kHBhalf + kHEhalf) - 1) {  // HO
1550         sd = HcalOuter;
1551         in -= 2 * (kHBhalf + kHEhalf);
1552         iz = (in < kHOhalf ? 1 : -1);
1553         in %= kHOhalf;
1554         dp = 4;
1555         ip = 1 + in / 15;
1556         ie = 1 + (in - 15 * (ip - 1));
1557       } else if (in > 2 * kHBhalf - 1) {  // Endcap
1558         sd = HcalEndcap;
1559         in -= 2 * kHBhalf;
1560         iz = (in < kHEhalf ? 1 : -1);
1561         in %= kHEhalf;
1562         ip = 2 * (in / 36);
1563         in %= 36;
1564         ip += 1 + in / 28;
1565         if (0 == ip % 2)
1566           in %= 28;
1567         ie = 15 + (in < 2 ? 1 + in : 2 + (in < 20 ? 1 + (in - 2) / 2 : 9 + (in < 26 ? 1 + (in - 20) / 3 : 3)));
1568         dp = (in < 1
1569                   ? 3
1570                   : (in < 2 ? 1 : (in < 20 ? 1 + (in - 2) % 2 : (in < 26 ? 1 + (in - 20) % 3 : (1 + (in - 26) % 2)))));
1571       } else {  // barrel
1572         iz = (in < kHBhalf ? 1 : -1);
1573         in %= kHBhalf;
1574         ip = in / 18 + 1;
1575         in %= 18;
1576         if (in < 14) {
1577           dp = 1;
1578           ie = in + 1;
1579         } else {
1580           in %= 14;
1581           dp = 1 + in % 2;
1582           ie = 15 + in / 2;
1583         }
1584       }
1585     }
1586   } else if (topoVersion_ == 10) {
1587     if (denseid < ncells()) {
1588       if (denseid >= (HBSize_ + HESize_ + HOSize_)) {
1589         sd = HcalForward;
1590         in -= (HBSize_ + HESize_ + HOSize_);
1591         dp = (in % maxDepthHF_) + 1;
1592         ip = (in - dp + 1) % (maxDepthHF_ * IPHI_MAX);
1593         ip = (ip / maxDepthHF_) + 1;
1594         ie = (in - dp + 1 - maxDepthHF_ * (ip - 1)) / (IPHI_MAX * maxDepthHF_);
1595         if (ie > 12) {
1596           ie = 54 - ie;
1597           iz = -1;
1598         } else {
1599           ie += 29;
1600           iz = 1;
1601         }
1602       } else if (denseid >= (HBSize_ + HESize_)) {
1603         sd = HcalOuter;
1604         in -= (HBSize_ + HESize_);
1605         dp = 4;
1606         ip = (in % IPHI_MAX) + 1;
1607         ie = (in - ip + 1) / IPHI_MAX;
1608         if (ie > 14) {
1609           ie = 30 - ie;
1610           iz = -1;
1611         } else {
1612           ie += 1;
1613           iz = 1;
1614         }
1615       } else if (denseid >= (HBSize_)) {
1616         sd = HcalEndcap;
1617         in -= (HBSize_);
1618         dp = (in % maxDepthHE_) + 1;
1619         ip = (in - dp + 1) % (maxDepthHE_ * maxPhiHE_);
1620         ip = (ip / maxDepthHE_) + 1;
1621         ie = (in - dp + 1 - maxDepthHE_ * (ip - 1)) / (maxPhiHE_ * maxDepthHE_);
1622         if (ie >= nEtaHE_) {
1623           ie = lastHERing() + nEtaHE_ - ie;
1624           iz = -1;
1625         } else {
1626           ie = firstHERing() + ie;
1627           iz = 1;
1628         }
1629       } else {
1630         sd = HcalBarrel;
1631         dp = (in % maxDepthHB_) + 1;
1632         ip = (in - dp + 1) % (maxDepthHB_ * IPHI_MAX);
1633         ip = (ip / maxDepthHB_) + 1;
1634         ie = (in - dp + 1 - maxDepthHB_ * (ip - 1)) / (IPHI_MAX * maxDepthHB_);
1635         if (ie >= nEtaHB_) {
1636           ie = lastHBRing() + nEtaHB_ - ie;
1637           iz = -1;
1638         } else {
1639           ie = firstHBRing() + ie;
1640           iz = 1;
1641         }
1642       }
1643     }
1644   }
1645   HcalDetId hid(sd, iz * int(ie), ip, dp);
1646 #ifdef EDM_ML_DEBUG
1647   edm::LogVerbatim("HcalGeom") << "Dens2Det " << topoVersion_ << " i/p " << std::hex << denseid << " : " << hid.rawId()
1648                                << std::dec << " | " << hid;
1649 #endif
1650   return hid;
1651 }
1652 
1653 unsigned int HcalTopology::ncells() const { return HBSize_ + HESize_ + HOSize_ + HFSize_; }
1654 
1655 unsigned int HcalTopology::ncells(int subdet) const {
1656   unsigned int ncell(0);
1657   for (int eta = -maxEta_; eta <= maxEta_; eta++) {
1658     for (unsigned int phi = 1; phi <= maxPhi_; phi++) {
1659       for (int depth = 1; depth <= 7; depth++) {
1660         HcalDetId cell((HcalSubdetector)subdet, eta, phi, depth);
1661         if (validRaw(cell, true))
1662           ++ncell;
1663       }
1664     }
1665   }
1666   return ncell;
1667 }
1668 
1669 int HcalTopology::topoVersion() const { return topoVersion_; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!