CaloTopology/src/CaloTowerTopology.cc

0001 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
0004 #include <cassert>
0005 #include <algorithm>
0006
0007 //#define EDM_ML_DEBUG
0008
0009 CaloTowerTopology::CaloTowerTopology(const HcalTopology* topology) : hcaltopo(topology) {
0010   //get number of towers in each hcal subdet from hcaltopo
0011   int nEtaHB, nEtaHO, nEtaHF;
0012   nEtaHB = hcaltopo->lastHBRing() - hcaltopo->firstHBRing() + 1;
0013   nEtaHE_ = hcaltopo->lastHERing() - hcaltopo->firstHERing() + 1;
0014   nEtaHO = hcaltopo->lastHORing() - hcaltopo->firstHORing() + 1;
0015   nEtaHF = hcaltopo->lastHFRing() - hcaltopo->firstHFRing() + 1;
0016 #ifdef EDM_ML_DEBUG
0017   edm::LogVerbatim("CaloTower") << "CaloTowerTopology:(1) " << nEtaHB << ":" << nEtaHE_ << ":" << nEtaHO << ":"
0018                                 << nEtaHF << ":" << hcaltopo->isBH();
0019 #endif
0020   if (hcaltopo->isBH())
0021     nEtaHE_ = 0;
0022
0023   //setup continuous ieta
0024   firstHBRing_ = 1;
0025   lastHBRing_ = firstHBRing_ + nEtaHB - 1;
0026   firstHERing_ = lastHBRing_;                             //crossover
0027   lastHERing_ = firstHERing_ + std::max(nEtaHE_ - 1, 0);  //max = protection for case with no HE
0028   //no crossover for CaloTowers; HF crossover cells go in the subsequent non-crossover HF tower
0029   firstHFRing_ = lastHERing_ + 1;
0030   lastHFRing_ = firstHFRing_ + (nEtaHF - 1) - 1;  //nEtaHF - 1 to account for no crossover
0031   firstHORing_ = 1;
0032   lastHORing_ = firstHORing_ + nEtaHO - 1;
0033 #ifdef EDM_ML_DEBUG
0034   edm::LogVerbatim("CaloTower") << "CaloTowerTopology: (2) " << firstHBRing_ << ":" << lastHBRing_ << ":"
0035                                 << firstHERing_ << ":" << lastHERing_ << ":" << firstHFRing_ << ":" << lastHFRing_
0036                                 << ":" << firstHORing_ << ":" << lastHORing_;
0037 #endif
0038
0039   //translate phi segmentation boundaries into continuous ieta
0040   if (hcaltopo->firstHEDoublePhiRing() == 999 || nEtaHE_ == 0)
0041     firstHEDoublePhiRing_ = firstHFRing_;
0042   else
0043     firstHEDoublePhiRing_ = firstHERing_ + (hcaltopo->firstHEDoublePhiRing() - hcaltopo->firstHERing());
0044   firstHEQuadPhiRing_ = firstHERing_ + (hcaltopo->firstHEQuadPhiRing() - hcaltopo->firstHERing());
0045   firstHFQuadPhiRing_ = firstHFRing_ - 1 + (hcaltopo->firstHFQuadPhiRing() - hcaltopo->firstHFRing());
0046
0047   //number of etas per phi segmentation type
0048   int nEtaSinglePhi_, nEtaDoublePhi_, nEtaQuadPhi_;
0049   nEtaSinglePhi_ = firstHEDoublePhiRing_ - firstHBRing_;
0050   nEtaDoublePhi_ = firstHFQuadPhiRing_ - firstHEDoublePhiRing_;
0051   nEtaQuadPhi_ = lastHFRing_ - firstHFQuadPhiRing_ + 1;  //include lastHFRing
0052
0053   //total number of towers per phi segmentation type
0054   nSinglePhi_ = nEtaSinglePhi_ * 72;
0055   nDoublePhi_ = nEtaDoublePhi_ * 36;
0056   nQuadPhi_ = nEtaQuadPhi_ * 18;
0057
0058 #ifdef EDM_ML_DEBUG
0059   edm::LogVerbatim("CaloTower") << "CaloTowerTopology: (3) " << nEtaSinglePhi_ << ":" << nEtaDoublePhi_ << ":"
0060                                 << nEtaQuadPhi_ << ":" << nSinglePhi_ << ":" << nDoublePhi_ << ":" << nQuadPhi_;
0061 #endif
0062
0063   //calculate maximum dense index size
0064   kSizeForDenseIndexing = 2 * (nSinglePhi_ + nDoublePhi_ + nQuadPhi_);
0065 }
0066
0067 //convert CaloTowerTopology ieta to HcalTopology ieta
0068 int CaloTowerTopology::convertCTtoHcal(int ct_ieta) const {
0069   if (ct_ieta <= lastHBRing_)
0070     return ct_ieta - firstHBRing_ + hcaltopo->firstHBRing();
0071   else if (ct_ieta <= lastHERing_)
0072     return ct_ieta - firstHERing_ + hcaltopo->firstHERing();
0073   else if (ct_ieta <= lastHFRing_)
0074     return ct_ieta - firstHFRing_ + hcaltopo->firstHFRing() + 1;  //account for no HF crossover
0075   else
0076     return 0;  //if ct_ieta outside range
0077 }
0078
0079 //convert HcalTopology ieta to CaloTowerTopology ieta
0080 int CaloTowerTopology::convertHcaltoCT(int hcal_ieta, HcalSubdetector subdet) const {
0081   if (subdet == HcalBarrel && hcal_ieta >= hcaltopo->firstHBRing() && hcal_ieta <= hcaltopo->lastHBRing()) {
0082     return hcal_ieta - hcaltopo->firstHBRing() + firstHBRing_;
0083   } else if (subdet == HcalEndcap && hcal_ieta >= hcaltopo->firstHERing() && hcal_ieta <= hcaltopo->lastHERing()) {
0084     return ((nEtaHE_ == 0) ? 0 : (hcal_ieta - hcaltopo->firstHERing() + firstHERing_));
0085   } else if (subdet == HcalForward && hcal_ieta >= hcaltopo->firstHFRing() && hcal_ieta <= hcaltopo->lastHFRing()) {
0086     if (hcal_ieta == hcaltopo->firstHFRing())
0087       hcal_ieta++;  //account for no HF crossover
0088     return hcal_ieta - hcaltopo->firstHFRing() + firstHFRing_ - 1;
0089   } else if (subdet == HcalOuter && hcal_ieta >= hcaltopo->firstHORing() && hcal_ieta <= hcaltopo->lastHORing()) {
0090     return hcal_ieta - hcaltopo->firstHORing() + firstHORing_;
0091   } else
0092     return 0;  //if hcal_ieta outside range, or unknown subdet
0093 }
0094
0095 bool CaloTowerTopology::valid(const DetId& id) const {
0096   assert(id.det() == DetId::Calo && id.subdetId() == CaloTowerDetId::SubdetId);
0097   return validDetId(id);
0098 }
0099
0100 bool CaloTowerTopology::validDetId(const CaloTowerDetId& id) const {
0101   int ia = id.ietaAbs();
0102   int ip = id.iphi();
0103
0104   return ((ia >= firstHBRing_) && (ia <= lastHFRing_)               //eta range
0105           && (ip >= 1) && (ip <= 72)                                //phi range
0106           && ((ia < firstHEDoublePhiRing_)                          //72 phi segments
0107               || (ia < firstHFQuadPhiRing_ && (ip - 1) % 2 == 0)    //36 phi segments, numbered 1,3,...,33,35
0108               || (ia >= firstHFQuadPhiRing_ && (ip - 3) % 4 == 0))  //18 phi segments, numbered 71,3,7,11,...
0109   );
0110 }
0111
0112 //decreasing ieta
0113 std::vector<DetId> CaloTowerTopology::east(const DetId& id) const {
0114   std::vector<DetId> dd;
0115   CaloTowerDetId tid(id);
0116   int ieta = tid.ieta();
0117   int iphi = tid.iphi();
0118
0119   if (ieta == 1) {  //no ieta=0
0120     ieta = -1;
0121   } else if (ieta == firstHEDoublePhiRing_) {
0122     //currently double phi, going to single phi (positive eta) -> extra neighbor
0123     ieta--;
0124     dd.emplace_back(CaloTowerDetId(ieta, iphi + 1));
0125   } else if (ieta - 1 == -firstHEDoublePhiRing_) {
0126     //currently single phi, going to double phi (negative eta) -> change numbering
0127     if ((iphi % 2) == 0)
0128       iphi--;
0129     ieta--;
0130   } else if (ieta == firstHFQuadPhiRing_) {  //currently quad phi, going to double phi (positive eta) -> extra neighbor
0131     ieta--;
0132     dd.emplace_back(CaloTowerDetId(ieta, ((iphi + 1) % 72) + 1));
0133   } else if (ieta - 1 == -firstHFQuadPhiRing_) {
0134     //currently double phi, going to quad phi (negative eta) -> change numbering
0135     if (((iphi - 1) % 4) == 0) {
0136       if (iphi == 1)
0137         iphi = 71;
0138       else
0139         iphi -= 2;
0140     }
0141     ieta--;
0142   } else {  //general case
0143     ieta--;
0144   }
0145
0146   if (ieta >= -lastHFRing_)
0147     dd.emplace_back(CaloTowerDetId(ieta, iphi));
0148   return dd;
0149 }
0150
0151 //increasing ieta
0152 std::vector<DetId> CaloTowerTopology::west(const DetId& id) const {
0153   std::vector<DetId> dd;
0154   CaloTowerDetId tid(id);
0155
0156   int ieta = tid.ieta();
0157   int iphi = tid.iphi();
0158
0159   if (ieta == -1) {  //no ieta=0
0160     ieta = 1;
0161   } else if (ieta == -firstHEDoublePhiRing_) {
0162     //currently double phi, going to single phi (negative eta) -> extra neighbor
0163     ieta++;
0164     dd.emplace_back(CaloTowerDetId(ieta, iphi + 1));
0165   } else if (ieta + 1 == firstHEDoublePhiRing_) {
0166     //currently single phi, going to double phi (positive eta) -> change numbering
0167     if ((iphi % 2) == 0)
0168       iphi--;
0169     ieta++;
0170   } else if (ieta == -firstHFQuadPhiRing_) {  //currently quad phi, going to double phi (negative eta) -> extra neighbor
0171     ieta++;
0172     dd.emplace_back(CaloTowerDetId(ieta, ((iphi + 1) % 72) + 1));
0173   } else if (ieta + 1 == firstHFQuadPhiRing_) {
0174     //currently double phi, going to quad phi (positive eta) -> change numbering
0175     if (((iphi - 1) % 4) == 0) {
0176       if (iphi == 1)
0177         iphi = 71;
0178       else
0179         iphi -= 2;
0180     }
0181     ieta++;
0182   } else {
0183     ieta++;
0184   }
0185
0186   if (ieta <= lastHFRing_)
0187     dd.emplace_back(CaloTowerDetId(ieta, iphi));
0188
0189   return dd;
0190 }
0191
0192 //increasing iphi
0193 std::vector<DetId> CaloTowerTopology::north(const DetId& id) const {
0194   CaloTowerDetId tid(id);
0195   int iphi_n = tid.iphi() + 1;
0196   if (iphi_n > 72)
0197     iphi_n = 1;
0198   if (tid.ietaAbs() >= firstHFQuadPhiRing_) {  //18 phi segments, numbered 71,3,7,11,...
0199     iphi_n += 3;
0200     if (iphi_n > 72)
0201       iphi_n -= 72;
0202   } else if (tid.ietaAbs() >= firstHEDoublePhiRing_ && (iphi_n % 2) == 0) {  //36 phi segments, numbered 1,3,...,33,35
0203     iphi_n++;
0204     if (iphi_n > 72)
0205       iphi_n -= 72;
0206   }
0207
0208   std::vector<DetId> dd;
0209   dd.emplace_back(CaloTowerDetId(tid.ieta(), iphi_n));
0210   return dd;
0211 }
0212
0213 //decreasing iphi
0214 std::vector<DetId> CaloTowerTopology::south(const DetId& id) const {
0215   CaloTowerDetId tid(id);
0216   int iphi_s = tid.iphi() - 1;
0217   if (iphi_s == 0)
0218     iphi_s = 72;
0219   if (tid.ietaAbs() >= firstHFQuadPhiRing_) {  //18 phi segments, numbered 71,3,7,11,...
0220     iphi_s -= 3;
0221     if (iphi_s <= 0)
0222       iphi_s += 72;
0223   } else if (tid.ietaAbs() >= firstHEDoublePhiRing_ && (iphi_s % 2) == 0) {  //36 phi segments, numbered 1,3,...,33,35
0224     iphi_s--;
0225   }
0226
0227   std::vector<DetId> dd;
0228   dd.emplace_back(CaloTowerDetId(tid.ieta(), iphi_s));
0229   return dd;
0230 }
0231
0232 std::vector<DetId> CaloTowerTopology::up(const DetId& /*id*/) const { return std::vector<DetId>(); }
0233
0234 std::vector<DetId> CaloTowerTopology::down(const DetId& /*id*/) const { return std::vector<DetId>(); }
0235
0236 uint32_t CaloTowerTopology::denseIndex(const DetId& id) const {
0237   CaloTowerDetId tid(id);
0238   const int ie(tid.ietaAbs());
0239   const int ip(tid.iphi() - 1);
0240
0241   return ((0 > tid.zside() ? 0 : kSizeForDenseIndexing / 2) +
0242           ((firstHEDoublePhiRing_ > ie
0243                 ? (ie - 1) * 72 + ip
0244                 : (firstHFQuadPhiRing_ > ie ? nSinglePhi_ + (ie - firstHEDoublePhiRing_) * 36 + ip / 2
0245                                             : nSinglePhi_ + nDoublePhi_ + (ie - firstHFQuadPhiRing_) * 18 + ip / 4))));
0246 }
0247
0248 CaloTowerDetId CaloTowerTopology::detIdFromDenseIndex(uint32_t din) const {
0249   const int iz(din < kSizeForDenseIndexing / 2 ? -1 : 1);
0250   din %= kSizeForDenseIndexing / 2;
0251   const int ie(nSinglePhi_ + nDoublePhi_ - 1 < (int)(din)
0252                    ? firstHFQuadPhiRing_ + (din - nSinglePhi_ - nDoublePhi_) / 18
0253                    : (nSinglePhi_ - 1 < (int)din ? firstHEDoublePhiRing_ + (din - nSinglePhi_) / 36 : din / 72 + 1));
0254
0255   const int ip(nSinglePhi_ + nDoublePhi_ - 1 < (int)(din)
0256                    ? ((din - nSinglePhi_ - nDoublePhi_) % 18) * 4 + 3
0257                    : (nSinglePhi_ - 1 < (int)(din) ? ((din - nSinglePhi_) % 36) * 2 + 1 : din % 72 + 1));
0258
0259   return (validDenseIndex(din) ? CaloTowerDetId(iz * ie, ip) : CaloTowerDetId());
0260 }