L1THGCal/src/HGCalTriggerTowerGeometryHelper.cc

0001
0002 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTowerGeometryHelper.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/Utilities/interface/EDMException.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include <cmath>
0009 #include <iostream>
0010 #include <fstream>
0011 #include <algorithm>
0012
0013 HGCalTriggerTowerGeometryHelper::HGCalTriggerTowerGeometryHelper(const edm::ParameterSet& conf)
0014     : doNose_(conf.getParameter<bool>("doNose")),
0015       minEta_(conf.getParameter<double>("minEta")),
0016       maxEta_(conf.getParameter<double>("maxEta")),
0017       minPhi_(conf.getParameter<double>("minPhi")),
0018       maxPhi_(conf.getParameter<double>("maxPhi")),
0019       nBinsEta_(conf.getParameter<int>("nBinsEta")),
0020       nBinsPhi_(conf.getParameter<int>("nBinsPhi")),
0021       binsEta_(conf.getParameter<std::vector<double> >("binsEta")),
0022       binsPhi_(conf.getParameter<std::vector<double> >("binsPhi")),
0023       splitModuleSum_(conf.getParameter<bool>("splitModuleSum")) {
0024   if (!binsEta_.empty() && ((unsigned int)(binsEta_.size()) != nBinsEta_ + 1)) {
0025     throw edm::Exception(edm::errors::Configuration, "Configuration")
0026         << "HGCalTriggerTowerGeometryHelper nBinsEta for the tower map not consistent with binsEta size" << std::endl;
0027   }
0028
0029   if (!binsPhi_.empty() && ((unsigned int)(binsPhi_.size()) != nBinsPhi_ + 1)) {
0030     throw edm::Exception(edm::errors::Configuration, "Configuration")
0031         << "HGCalTriggerTowerGeometryHelper nBinsPhi for the tower map not consistent with binsPhi size" << std::endl;
0032   }
0033
0034   // if the bin vecctor is empty we assume the bins to be regularly spaced
0035   if (binsEta_.empty()) {
0036     for (unsigned int bin1 = 0; bin1 != nBinsEta_ + 1; bin1++) {
0037       binsEta_.push_back(minEta_ + bin1 * ((maxEta_ - minEta_) / nBinsEta_));
0038     }
0039   }
0040
0041   // if the bin vecctor is empty we assume the bins to be regularly spaced
0042   if (binsPhi_.empty()) {
0043     for (unsigned int bin2 = 0; bin2 != nBinsPhi_ + 1; bin2++) {
0044       binsPhi_.push_back(minPhi_ + bin2 * ((maxPhi_ - minPhi_) / nBinsPhi_));
0045     }
0046   }
0047
0048   for (int zside = -1; zside <= 1; zside += 2) {
0049     for (unsigned int bin1 = 0; bin1 != nBinsEta_; bin1++) {
0050       for (unsigned int bin2 = 0; bin2 != nBinsPhi_; bin2++) {
0051         l1t::HGCalTowerID towerId(doNose_, zside, bin1, bin2);
0052         tower_coords_.emplace_back(towerId.rawId(),
0053                                    zside * ((binsEta_[bin1 + 1] + binsEta_[bin1]) / 2),
0054                                    (binsPhi_[bin2 + 1] + binsPhi_[bin2]) / 2);
0055       }
0056     }
0057   }
0058
0059   if (conf.getParameter<bool>("readMappingFile")) {
0060     // We read the TC to TT mapping from file,
0061     // otherwise we derive the TC to TT mapping on the fly from eta-phi coord. of the TCs
0062     std::ifstream l1tTriggerTowerMappingStream(conf.getParameter<edm::FileInPath>("L1TTriggerTowerMapping").fullPath());
0063     if (!l1tTriggerTowerMappingStream.is_open()) {
0064       throw cms::Exception("MissingDataFile") << "Cannot open HGCalTriggerGeometry L1TTriggerTowerMapping file\n";
0065     }
0066
0067     unsigned trigger_cell_id = 0;
0068     unsigned short iEta = 0;
0069     unsigned short iPhi = 0;
0070
0071     for (; l1tTriggerTowerMappingStream >> trigger_cell_id >> iEta >> iPhi;) {
0072       if (iEta >= nBinsEta_ || iPhi >= nBinsPhi_) {
0073         throw edm::Exception(edm::errors::Configuration, "Configuration")
0074             << "HGCalTriggerTowerGeometryHelper warning inconsistent mapping TC : " << trigger_cell_id
0075             << " to TT iEta: " << iEta << " iPhi: " << iPhi << " when max #bins eta: " << nBinsEta_
0076             << " phi: " << nBinsPhi_ << std::endl;
0077       }
0078       l1t::HGCalTowerID towerId(doNose_, triggerTools_.zside(DetId(trigger_cell_id)), iEta, iPhi);
0079       cells_to_trigger_towers_[trigger_cell_id] = towerId.rawId();
0080     }
0081     l1tTriggerTowerMappingStream.close();
0082   }
0083
0084   if (splitModuleSum_) {
0085     //variables for transforming towers
0086     rotate180Deg_ = int(nBinsPhi_) / 2;
0087     rotate120Deg_ = int(nBinsPhi_) / 3;
0088     reverseX_ = int(nBinsPhi_) / 2 - 1;
0089
0090     std::ifstream moduleTowerMappingStream(conf.getParameter<edm::FileInPath>("moduleTowerMapping").fullPath());
0091     if (!moduleTowerMappingStream.is_open()) {
0092       throw cms::Exception("MissingDataFile") << "Cannot open HGCalTowerMapProducer moduleTowerMapping file\n";
0093     }
0094     //get split divisors
0095     std::string line;
0096     std::getline(moduleTowerMappingStream, line);  //Skip row
0097     std::getline(moduleTowerMappingStream, line);
0098     std::stringstream ss(line);
0099     ss >> splitDivisorSilic_ >> splitDivisorScint_;
0100
0101     //get towers and module sum shares
0102     std::getline(moduleTowerMappingStream, line);  //Skip row
0103     std::getline(moduleTowerMappingStream, line);  //Skip row
0104     const int minNumOfWordsPerRow = 5;
0105     const int numOfWordsPerTower = 3;
0106     for (std::string line; std::getline(moduleTowerMappingStream, line);) {
0107       int numOfWordsInThisRow = 0;
0108       for (std::string::size_type i = 0; i < line.size(); i++) {
0109         if (line[i] != ' ' && line[i + 1] == ' ') {
0110           numOfWordsInThisRow++;
0111         }
0112       }
0113       if (numOfWordsInThisRow < minNumOfWordsPerRow) {
0114         throw edm::Exception(edm::errors::Configuration, "Configuration")
0115             << "HGCalTriggerTowerGeometryHelper warning: Incorrect/incomplete values for module ID in the mapping "
0116                "file.\n"
0117             << "The incorrect line is:" << line << std::endl;
0118       }
0119       int subdet = 0;
0120       int layer = 0;
0121       int moduleU = 0;
0122       int moduleV = 0;
0123       int numTowers = 0;
0124       std::stringstream ss(line);
0125       ss >> subdet >> layer >> moduleU >> moduleV >> numTowers;
0126       if (numOfWordsInThisRow != (numTowers * numOfWordsPerTower + minNumOfWordsPerRow)) {
0127         throw edm::Exception(edm::errors::Configuration, "Configuration")
0128             << "HGCalTriggerTowerGeometryHelper warning: Incorrect/incomplete values for module ID or tower "
0129                "share/eta/phi in the mapping file.\n"
0130             << "The incorrect line is:" << line << std::endl;
0131       }
0132       unsigned packed_modID = packLayerSubdetWaferId(subdet, layer, moduleU, moduleV);
0133       std::vector<unsigned> towers;
0134       for (int itr_tower = 0; itr_tower < numTowers; itr_tower++) {
0135         int iEta_raw = 0;
0136         int iPhi_raw = 0;
0137         int towerShare = 0;
0138         ss >> iEta_raw >> iPhi_raw >> towerShare;
0139         int splitDivisor = (subdet == 2) ? splitDivisorScint_ : splitDivisorSilic_;
0140         if ((towerShare > splitDivisor) || (towerShare < 1)) {
0141           throw edm::Exception(edm::errors::Configuration, "Configuration")
0142               << "HGCalTriggerTowerGeometryHelper warning: invalid tower share in the mapping file.\n"
0143               << "Tower share must be a positive integer and less than splitDivisor. The incorrect values found for "
0144                  "module ID:"
0145               << std::endl
0146               << "subdet=" << subdet << ", l=" << layer << ", u=" << moduleU << ", v=" << moduleV << std::endl;
0147         }
0148         towers.push_back(packTowerIDandShare(iEta_raw, iPhi_raw, towerShare));
0149       }
0150       modules_to_trigger_towers_[packed_modID] = towers;
0151     }
0152     moduleTowerMappingStream.close();
0153   }
0154 }
0155
0156 unsigned HGCalTriggerTowerGeometryHelper::packLayerSubdetWaferId(int subdet, int layer, int moduleU, int moduleV) const {
0157   unsigned packed_modID = 0;
0158   packed_modID |= ((subdet & HGCalTriggerModuleDetId::kHGCalTriggerSubdetMask)
0159                    << HGCalTriggerModuleDetId::kHGCalTriggerSubdetOffset);
0160   packed_modID |= ((layer & HGCalTriggerModuleDetId::kHGCalLayerMask) << HGCalTriggerModuleDetId::kHGCalLayerOffset);
0161   packed_modID |=
0162       ((moduleU & HGCalTriggerModuleDetId::kHGCalModuleUMask) << HGCalTriggerModuleDetId::kHGCalModuleUOffset);
0163   packed_modID |=
0164       ((moduleV & HGCalTriggerModuleDetId::kHGCalModuleVMask) << HGCalTriggerModuleDetId::kHGCalModuleVOffset);
0165   return packed_modID;
0166 }
0167
0168 unsigned HGCalTriggerTowerGeometryHelper::packTowerIDandShare(int iEta_raw, int iPhi_raw, int towerShare) const {
0169   unsigned packed_towerIDandShare = 0;
0170   unsigned iEtaAbs = std::abs(iEta_raw);
0171   unsigned iEtaSign = std::signbit(iEta_raw);
0172   unsigned iPhiAbs = std::abs(iPhi_raw);
0173   unsigned iPhiSign = std::signbit(iPhi_raw);
0174   packed_towerIDandShare |= ((iEtaAbs & l1t::HGCalTowerID::coordMask) << l1t::HGCalTowerID::coord1Shift);
0175   packed_towerIDandShare |= ((iEtaSign & signMask) << sign1Shift);
0176   packed_towerIDandShare |= ((iPhiAbs & l1t::HGCalTowerID::coordMask) << l1t::HGCalTowerID::coord2Shift);
0177   packed_towerIDandShare |= ((iPhiSign & signMask) << sign2Shift);
0178   packed_towerIDandShare |= ((towerShare & towerShareMask) << towerShareShift);
0179   return packed_towerIDandShare;
0180 }
0181
0182 void HGCalTriggerTowerGeometryHelper::unpackTowerIDandShare(unsigned towerIDandShare,
0183                                                             int& iEta_raw,
0184                                                             int& iPhi_raw,
0185                                                             int& towerShare) const {
0186   //eta
0187   iEta_raw = (towerIDandShare >> l1t::HGCalTowerID::coord1Shift) & l1t::HGCalTowerID::coordMask;
0188   unsigned iEtaSign = (towerIDandShare >> sign1Shift) & signMask;
0189   iEta_raw = (iEtaSign) ? -1 * iEta_raw : iEta_raw;
0190   //phi
0191   iPhi_raw = (towerIDandShare >> l1t::HGCalTowerID::coord2Shift) & l1t::HGCalTowerID::coordMask;
0192   unsigned iPhiSign = (towerIDandShare >> sign2Shift) & signMask;
0193   iPhi_raw = (iPhiSign) ? -1 * iPhi_raw : iPhi_raw;
0194   //tower share
0195   towerShare = (towerIDandShare >> towerShareShift) & towerShareMask;
0196 }
0197
0198 int HGCalTriggerTowerGeometryHelper::moveToCorrectSector(int iPhi_raw, int sector) const {
0199   int iPhi = (iPhi_raw + sector * rotate120Deg_ + rotate180Deg_) % int(nBinsPhi_);
0200   return iPhi;
0201 }
0202
0203 void HGCalTriggerTowerGeometryHelper::reverseXaxis(int& iPhi) const {
0204   iPhi = reverseX_ - iPhi;                          //correct x -> -x in z>0
0205   iPhi = (int(nBinsPhi_) + iPhi) % int(nBinsPhi_);  // make all phi between 0 to nBinsPhi_-1
0206 }
0207
0208 const std::vector<l1t::HGCalTowerCoord>& HGCalTriggerTowerGeometryHelper::getTowerCoordinates() const {
0209   return tower_coords_;
0210 }
0211
0212 unsigned short HGCalTriggerTowerGeometryHelper::getTriggerTowerFromEtaPhi(const float& eta, const float& phi) const {
0213   auto bin_eta_l = std::lower_bound(binsEta_.begin(), binsEta_.end(), fabs(eta));
0214   unsigned int bin_eta = 0;
0215   // we add a protection for TCs in Hadron part which are outside the boundaries and possible rounding effects
0216   if (bin_eta_l == binsEta_.end()) {
0217     if (fabs(eta) < minEta_) {
0218       bin_eta = 0;
0219     } else if (fabs(eta) >= maxEta_) {
0220       bin_eta = nBinsEta_;
0221     } else {
0222       edm::LogError("HGCalTriggerTowerGeometryHelper")
0223           << " did not manage to map eta " << eta << " to any Trigger Tower\n";
0224     }
0225   } else {
0226     bin_eta = bin_eta_l - binsEta_.begin() - 1;
0227   }
0228
0229   auto bin_phi_l = std::lower_bound(binsPhi_.begin(), binsPhi_.end(), phi);
0230   unsigned int bin_phi = 0;
0231   if (bin_phi_l == binsPhi_.end()) {
0232     if (phi < minPhi_) {
0233       bin_phi = nBinsPhi_;
0234     } else if (phi >= maxPhi_) {
0235       bin_phi = 0;
0236     } else {
0237       edm::LogError("HGCalTriggerTowerGeometryHelper")
0238           << " did not manage to map phi " << phi << " to any Trigger Tower\n";
0239     }
0240   } else {
0241     bin_phi = bin_phi_l - binsPhi_.begin() - 1;
0242   }
0243   int zside = eta < 0 ? -1 : 1;
0244   return l1t::HGCalTowerID(doNose_, zside, bin_eta, bin_phi).rawId();
0245 }
0246
0247 std::unordered_map<unsigned short, float> HGCalTriggerTowerGeometryHelper::getTriggerTower(
0248     const l1t::HGCalTriggerCell& thecell) const {
0249   std::unordered_map<unsigned short, float> towerIDandShares = {};
0250   unsigned int trigger_cell_id = thecell.detId();
0251   // NOTE: if the TC is not found in the map than it is mapped via eta-phi coords.
0252   // this can be considered dangerous (silent failure of the map) but it actually allows to save
0253   // memory mapping explicitly only what is actually needed
0254   auto tower_id_itr = cells_to_trigger_towers_.find(trigger_cell_id);
0255   if (tower_id_itr != cells_to_trigger_towers_.end()) {
0256     towerIDandShares.insert({tower_id_itr->second, 1.0});
0257     return towerIDandShares;
0258   }
0259   towerIDandShares.insert({getTriggerTowerFromEtaPhi(thecell.position().eta(), thecell.position().phi()), 1.0});
0260   return towerIDandShares;
0261 }
0262
0263 std::unordered_map<unsigned short, float> HGCalTriggerTowerGeometryHelper::getTriggerTower(
0264     const l1t::HGCalTriggerSums& thesum) const {
0265   std::unordered_map<unsigned short, float> towerIDandShares = {};
0266   if (!splitModuleSum_) {
0267     towerIDandShares.insert({getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi()), 1.0});
0268     return towerIDandShares;
0269   } else {
0270     HGCalTriggerModuleDetId detid(thesum.detId());
0271     int moduleU = detid.moduleU();
0272     int moduleV = detid.moduleV();
0273     int layer = detid.layer();
0274     int sector = detid.sector();
0275     int zside = detid.zside();
0276     int subdet = 0;
0277     int splitDivisor;
0278     if (detid.isHScintillator()) {
0279       subdet = 2;
0280       splitDivisor = splitDivisorScint_;
0281     } else if (detid.isEE()) {
0282       subdet = 0;
0283       splitDivisor = splitDivisorSilic_;
0284     } else if (detid.isHSilicon()) {
0285       subdet = 1;
0286       splitDivisor = splitDivisorSilic_;
0287     } else {  //HFNose
0288       towerIDandShares.insert({getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi()), 1.0});
0289       return towerIDandShares;
0290     }
0291
0292     unsigned packed_modID = packLayerSubdetWaferId(subdet, layer, moduleU, moduleV);
0293     auto module_id_itr = modules_to_trigger_towers_.find(packed_modID);
0294     if (module_id_itr != modules_to_trigger_towers_.end()) {
0295       //eta variables
0296       int iEta = -999;
0297       int iEta_raw = -999;
0298       int offsetEta = 2;
0299       //phi variables
0300       int iPhi = -999;
0301       int iPhi_raw = -999;
0302       int towerShare = -999;  //the share each tower gets from module sum
0303       for (auto towerIDandShare : module_id_itr->second) {
0304         unpackTowerIDandShare(towerIDandShare, iEta_raw, iPhi_raw, towerShare);
0305         iEta = offsetEta + iEta_raw;
0306         iPhi = moveToCorrectSector(iPhi_raw, sector);
0307         if (zside == 1) {
0308           reverseXaxis(iPhi);
0309         }
0310         towerIDandShares.insert(
0311             {l1t::HGCalTowerID(doNose_, zside, iEta, iPhi).rawId(), double(towerShare) / splitDivisor});
0312       }
0313       return towerIDandShares;
0314     } else {  // for modules not found in the mapping file (currently a few partial modules) use the traditional method.
0315       towerIDandShares.insert({getTriggerTowerFromEtaPhi(thesum.position().eta(), thesum.position().phi()), 1.0});
0316       return towerIDandShares;
0317     }
0318   }
0319 }