HcalTowerAlgo/src/HcalTrigTowerGeometry.cc

0001 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0003 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0004 #include "DataFormats/HcalDetId/interface/HcalTrigTowerDetId.h"
0005
0006 #include <iostream>
0007 #include <cassert>
0008
0009 HcalTrigTowerGeometry::HcalTrigTowerGeometry(const HcalTopology* topology) : theTopology(topology) {
0010   auto tmode = theTopology->triggerMode();
0011   useRCT_ = tmode <= HcalTopologyMode::TriggerMode_2016;
0012   use1x1_ = tmode >= HcalTopologyMode::TriggerMode_2016;
0013   use2017_ = tmode >= HcalTopologyMode::TriggerMode_2017 or tmode == HcalTopologyMode::TriggerMode_2018legacy;
0014 }
0015
0016 std::vector<HcalTrigTowerDetId> HcalTrigTowerGeometry::towerIds(const HcalDetId& cellId) const {
0017   std::vector<HcalTrigTowerDetId> results;
0018
0019   if (cellId.subdet() == HcalForward) {
0020     if (useRCT_) {
0021       // first do eta

0022       int hfRing = cellId.ietaAbs();
0023       int ieta = firstHFTower(0);
0024       // find the tower that contains this ring

0025       while (hfRing >= firstHFRingInTower(ieta + 1)) {
0026         ++ieta;
0027       }
0028
0029       ieta *= cellId.zside();
0030
0031       // now for phi

0032       // HF towers are quad, 18 in phi.

0033       // go two cells per trigger tower.

0034
0035       int iphi = (((cellId.iphi() + 1) / 4) * 4 + 1) % 72;  // 71+1 --> 1, 3+5 --> 5

0036       results.emplace_back(HcalTrigTowerDetId(ieta, iphi));
0037     }
0038     if (use1x1_) {
0039       int hfRing = cellId.ietaAbs();
0040       if (hfRing == 29)
0041         hfRing = 30;  // sum 29 into 30.

0042
0043       int ieta = hfRing * cellId.zside();
0044       int iphi = cellId.iphi();
0045
0046       HcalTrigTowerDetId id(ieta, iphi);
0047       id.setVersion(1);  // version 1 for 1x1 HF granularity

0048       results.emplace_back(id);
0049     }
0050
0051   } else {
0052     // the first twenty rings are one-to-one

0053     if (cellId.ietaAbs() <= theTopology->lastHBRing()) {
0054       results.emplace_back(HcalTrigTowerDetId(cellId.ieta(), cellId.iphi()));
0055     } else if (theTopology->maxDepthHE() == 0) {
0056       // Ignore these

0057     } else if (cellId.ietaAbs() < theTopology->firstHEDoublePhiRing()) {
0058       results.emplace_back(HcalTrigTowerDetId(cellId.ieta(), cellId.iphi()));
0059     } else {
0060       // the remaining rings are two-to-one in phi

0061       int iphi1 = cellId.iphi();
0062       int ieta = cellId.ieta();
0063       int depth = cellId.depth();
0064       // the last eta ring in HE is split.  Recombine.

0065       if (ieta == theTopology->lastHERing())
0066         --ieta;
0067       if (ieta == -theTopology->lastHERing())
0068         ++ieta;
0069
0070       if (use2017_) {
0071         if (ieta == 26 and depth == 7)
0072           ++ieta;
0073         if (ieta == -26 and depth == 7)
0074           --ieta;
0075       }
0076
0077       results.emplace_back(HcalTrigTowerDetId(ieta, iphi1));
0078       results.emplace_back(HcalTrigTowerDetId(ieta, iphi1 + 1));
0079     }
0080   }
0081
0082   return results;
0083 }
0084
0085 std::vector<HcalTrigTowerDetId> HcalTrigTowerGeometry::towerIds_ZDC(const HcalZDCDetId& cellId) const {
0086   std::vector<HcalTrigTowerDetId> results;
0087
0088   int ieta = cellId.zside() > 0 ? 42 : -42;
0089   int iphi = 99;
0090   int channelId = cellId.channel();
0091   if (cellId.depth() != 1)
0092     channelId += 5;
0093
0094   results.emplace_back(HcalTrigTowerDetId(ieta, channelId));
0095   results.emplace_back(HcalTrigTowerDetId(ieta, iphi));
0096
0097   return results;
0098 }
0099
0100 std::vector<HcalDetId> HcalTrigTowerGeometry::detIds(const HcalTrigTowerDetId& hcalTrigTowerDetId) const {
0101   // Written, tested by E. Berry (Princeton)

0102   std::vector<HcalDetId> results;
0103
0104   int tower_ieta = hcalTrigTowerDetId.ieta();
0105   int tower_iphi = hcalTrigTowerDetId.iphi();
0106
0107   int cell_ieta = tower_ieta;
0108   int cell_iphi = tower_iphi;
0109
0110   int min_depth, n_depths;
0111
0112   // HB

0113
0114   if (abs(cell_ieta) <= theTopology->lastHBRing()) {
0115     theTopology->depthBinInformation(
0116         HcalBarrel, abs(tower_ieta), tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0117     for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++)
0118       results.emplace_back(HcalDetId(HcalBarrel, cell_ieta, cell_iphi, cell_depth));
0119   }
0120
0121   // HO

0122
0123   if (abs(cell_ieta) <= theTopology->lastHORing()) {
0124     theTopology->depthBinInformation(
0125         HcalOuter, abs(tower_ieta), tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0126     for (int ho_depth = min_depth; ho_depth <= min_depth + n_depths - 1; ho_depth++)
0127       results.emplace_back(HcalDetId(HcalOuter, cell_ieta, cell_iphi, ho_depth));
0128   }
0129
0130   // HE

0131
0132   if (abs(cell_ieta) >= theTopology->firstHERing() && abs(cell_ieta) < theTopology->lastHERing()) {
0133     theTopology->depthBinInformation(
0134         HcalEndcap, abs(tower_ieta), tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0135
0136     // Special for double-phi cells

0137     if (abs(cell_ieta) >= theTopology->firstHEDoublePhiRing())
0138       if (tower_iphi % 2 == 0)
0139         cell_iphi = tower_iphi - 1;
0140
0141     if (use2017_) {
0142       if (abs(tower_ieta) == 26)
0143         --n_depths;
0144       if (tower_ieta == 27)
0145         results.emplace_back(HcalDetId(HcalEndcap, cell_ieta - 1, cell_iphi, 7));
0146       if (tower_ieta == -27)
0147         results.emplace_back(HcalDetId(HcalEndcap, cell_ieta + 1, cell_iphi, 7));
0148     }
0149
0150     for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++)
0151       results.emplace_back(HcalDetId(HcalEndcap, cell_ieta, cell_iphi, cell_depth));
0152
0153     // Special for split-eta cells

0154     if (abs(tower_ieta) == 28) {
0155       theTopology->depthBinInformation(
0156           HcalEndcap, abs(tower_ieta) + 1, tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0157       for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++) {
0158         if (tower_ieta < 0)
0159           results.emplace_back(HcalDetId(HcalEndcap, tower_ieta - 1, cell_iphi, cell_depth));
0160         if (tower_ieta > 0)
0161           results.emplace_back(HcalDetId(HcalEndcap, tower_ieta + 1, cell_iphi, cell_depth));
0162       }
0163     }
0164   }
0165
0166   // HF

0167
0168   if (abs(cell_ieta) >= theTopology->firstHFRing()) {
0169     if (hcalTrigTowerDetId.version() == 0) {
0170       int HfTowerPhiSize = 72 / nPhiBins(tower_ieta, 0);
0171
0172       int HfTowerEtaSize = hfTowerEtaSize(tower_ieta);
0173       int FirstHFRingInTower = firstHFRingInTower(abs(tower_ieta));
0174
0175       for (int iHFTowerPhiSegment = 0; iHFTowerPhiSegment < HfTowerPhiSize; iHFTowerPhiSegment++) {
0176         cell_iphi = (tower_iphi / HfTowerPhiSize) * HfTowerPhiSize;  // Find the minimum phi segment

0177         cell_iphi -= 2;                   // The first trigger tower starts at HCAL iphi = 71, not HCAL iphi = 1

0178         cell_iphi += iHFTowerPhiSegment;  // Get all of the HCAL iphi values in this trigger tower

0179         cell_iphi += 72;                  // Don't want to take the mod of a negative number

0180         cell_iphi = cell_iphi % 72;       // There are, at most, 72 cells.

0181         cell_iphi += 1;                   // There is no cell at iphi = 0

0182
0183         if (cell_iphi % 2 == 0)
0184           continue;  // These cells don't exist.

0185
0186         for (int iHFTowerEtaSegment = 0; iHFTowerEtaSegment < HfTowerEtaSize; iHFTowerEtaSegment++) {
0187           cell_ieta = FirstHFRingInTower + iHFTowerEtaSegment;
0188
0189           if (cell_ieta >= 40 && cell_iphi % 4 == 1)
0190             continue;  // These cells don't exist.

0191
0192           theTopology->depthBinInformation(
0193               HcalForward, cell_ieta, tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0194
0195           // Negative tower_ieta -> negative cell_ieta

0196           int zside = 1;
0197           if (tower_ieta < 0)
0198             zside = -1;
0199
0200           cell_ieta *= zside;
0201
0202           for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++)
0203             results.emplace_back(HcalDetId(HcalForward, cell_ieta, cell_iphi, cell_depth));
0204
0205           if (zside * cell_ieta == 30) {
0206             theTopology->depthBinInformation(
0207                 HcalForward, 29 * zside, tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0208             for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++)
0209               results.emplace_back(HcalDetId(HcalForward, 29 * zside, cell_iphi, cell_depth));
0210           }
0211         }
0212       }
0213     } else if (hcalTrigTowerDetId.version() == 1) {
0214       theTopology->depthBinInformation(
0215           HcalForward, tower_ieta, tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0216       for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++)
0217         results.emplace_back(HcalDetId(HcalForward, tower_ieta, tower_iphi, cell_depth));
0218       if (abs(tower_ieta) == 30) {
0219         int i29 = 29;
0220         if (tower_ieta < 0)
0221           i29 = -29;
0222         theTopology->depthBinInformation(HcalForward, i29, tower_iphi, hcalTrigTowerDetId.zside(), n_depths, min_depth);
0223         for (int cell_depth = min_depth; cell_depth <= min_depth + n_depths - 1; cell_depth++)
0224           results.emplace_back(HcalDetId(HcalForward, i29, tower_iphi, cell_depth));
0225       }
0226     }
0227   }
0228
0229   return results;
0230 }
0231
0232 int HcalTrigTowerGeometry::hfTowerEtaSize(int ieta) const {
0233   int ietaAbs = abs(ieta);
0234   assert(ietaAbs >= firstHFTower(0) && ietaAbs <= nTowers(0));
0235   // the first three come from rings 29-31, 32-34, 35-37. The last has 4 rings: 38-41

0236   return (ietaAbs == nTowers(0)) ? 4 : 3;
0237 }
0238
0239 int HcalTrigTowerGeometry::firstHFRingInTower(int ietaTower) const {
0240   // count up to the correct HF ring

0241   int inputTower = abs(ietaTower);
0242   int result = theTopology->firstHFRing();
0243   for (int iTower = firstHFTower(0); iTower != inputTower; ++iTower) {
0244     result += hfTowerEtaSize(iTower);
0245   }
0246
0247   // negative in, negative out.

0248   if (ietaTower < 0)
0249     result *= -1;
0250   return result;
0251 }
0252
0253 void HcalTrigTowerGeometry::towerEtaBounds(int ieta, int version, double& eta1, double& eta2) const {
0254   int ietaAbs = abs(ieta);
0255   std::pair<double, double> etas = (ietaAbs < firstHFTower(version)) ? theTopology->etaRange(HcalBarrel, ietaAbs)
0256                                                                      : theTopology->etaRange(HcalForward, ietaAbs);
0257   eta1 = etas.first;
0258   eta2 = etas.second;
0259
0260   // get the signs and order right

0261   if (ieta < 0) {
0262     double tmp = eta1;
0263     eta1 = -eta2;
0264     eta2 = -tmp;
0265   }
0266 }