File indexing completed on 2024-09-10 02:59:04
0001
0002
0003
0004
0005 #include "SimCalorimetry/EcalSelectiveReadoutAlgos/interface/EcalSelectiveReadout.h"
0006 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0007 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0008 #include "FWCore/Utilities/interface/EDMException.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include <string>
0011 #include <iomanip>
0012 #include <cassert>
0013 #include <atomic>
0014
0015
0016 using std::vector;
0017
0018 const char EcalSelectiveReadout::srpFlagMarker[] = {'.', 'S', 'N', 'C', '4', '5', '6', '7'};
0019
0020 EcalSelectiveReadout::EcalSelectiveReadout(int dEta_, int dPhi_)
0021 : theTriggerMap(nullptr), theElecMap(nullptr), dEta(dEta_), dPhi(dPhi_) {}
0022
0023 void EcalSelectiveReadout::resetEeRuInterest() {
0024
0025 for (size_t iCap = 0; iCap < nEndcaps; ++iCap) {
0026 for (int iDccPhi = 0; iDccPhi < nDccPerEe; ++iDccPhi) {
0027 for (int iDccCh = 0; iDccCh < maxDccChs; ++iDccCh) {
0028 eeRuInterest_[iCap][iDccPhi][iDccCh] = UNKNOWN;
0029 }
0030 }
0031 }
0032 }
0033
0034 void EcalSelectiveReadout::runSelectiveReadout0(const ttFlag_t ttFlags[nTriggerTowersInEta][nTriggerTowersInPhi]) {
0035
0036
0037
0038 classifyTriggerTowers(ttFlags);
0039
0040
0041 int nTriggerTowerE[] = {0, 0, 0, 0, 0, 0, 0, 0};
0042 int nTriggerTowerB[] = {0, 0, 0, 0, 0, 0, 0, 0};
0043 static std::atomic<int> ncall{0};
0044 if (ncall < 10) {
0045 ++ncall;
0046 for (size_t iPhi = 0; iPhi < nTriggerTowersInPhi; ++iPhi) {
0047 for (size_t iEta = 0; iEta < nTriggerTowersInEta; ++iEta) {
0048 if (iEta < nEndcapTriggerTowersInEta || iEta >= nBarrelTriggerTowersInEta + nEndcapTriggerTowersInEta) {
0049
0050 ++nTriggerTowerE[towerInterest[iEta][iPhi]];
0051 } else {
0052 ++nTriggerTowerB[towerInterest[iEta][iPhi]];
0053 }
0054
0055 assert(towerInterest[iEta][iPhi] >= 0 && towerInterest[iEta][iPhi] <= 0x7);
0056 }
0057 }
0058 edm::LogInfo("EcalSelectiveReadout") << "without forced bit + with forced bit set:\n"
0059 << nTriggerTowerB[LOWINTEREST] << " + "
0060 << nTriggerTowerB[LOWINTEREST | FORCED_MASK]
0061 << " low interest TT(s) in barrel\n"
0062 << nTriggerTowerB[SINGLE] << " + " << nTriggerTowerB[SINGLE | FORCED_MASK]
0063 << " single TT(s) in barrel\n"
0064 << nTriggerTowerB[NEIGHBOUR] << " + "
0065 << nTriggerTowerB[NEIGHBOUR | FORCED_MASK]
0066 << " neighbor interest TT(s) in barrel\n"
0067 << nTriggerTowerB[CENTER] << " + " << nTriggerTowerB[CENTER | FORCED_MASK]
0068 << " centre interest TT(s) in barrel\n"
0069 << nTriggerTowerE[LOWINTEREST] << " + "
0070 << nTriggerTowerE[LOWINTEREST | FORCED_MASK]
0071 << " low interest TT(s) in endcap\n"
0072 << nTriggerTowerE[SINGLE] << " + " << nTriggerTowerE[SINGLE | FORCED_MASK]
0073 << " single TT(s) in endcap\n"
0074 << nTriggerTowerE[NEIGHBOUR] << " + "
0075 << nTriggerTowerE[NEIGHBOUR | FORCED_MASK] << " neighbor TT(s) in endcap\n"
0076 << nTriggerTowerE[CENTER] << " + " << nTriggerTowerE[CENTER | FORCED_MASK]
0077 << " center TT(s) in endcap\n";
0078 }
0079
0080
0081
0082 resetEeRuInterest();
0083
0084 #ifndef ECALSELECTIVEREADOUT_NOGEOM
0085 const std::vector<DetId>& endcapDetIds = theGeometry->getValidDetIds(DetId::Ecal, EcalEndcap);
0086 for (std::vector<DetId>::const_iterator eeDetIdItr = endcapDetIds.begin(); eeDetIdItr != endcapDetIds.end();
0087 ++eeDetIdItr) {
0088
0089
0090
0091
0092 EcalTrigTowerDetId trigTower = theTriggerMap->towerOf(*eeDetIdItr);
0093 assert(trigTower.rawId() != 0);
0094 EEDetId eeDetId(*eeDetIdItr);
0095 int iz = (eeDetId.zside() > 0) ? 1 : 0;
0096
0097
0098
0099 combineFlags(eeRuInterest(eeDetId), getTowerInterest(trigTower));
0100 }
0101 #else
0102 EEDetId xtal;
0103 for (int iZ0 = 0; iZ0 < 2; ++iZ0) {
0104 for (unsigned iX0 = 0; iX0 < nEndcapXBins; ++iX0) {
0105 for (unsigned iY0 = 0; iY0 < nEndcapYBins; ++iY0) {
0106 if (!(xtal.validDetId(iX0 + 1, iY0 + 1, (iZ0 > 0 ? 1 : -1)))) {
0107 continue;
0108 }
0109 xtal = EEDetId(iX0 + 1, iY0 + 1, (iZ0 > 0 ? 1 : -1));
0110
0111 if (39 <= iX0 && iX0 <= 60 && 45 <= iY0 && iY0 <= 54) {
0112 continue;
0113 }
0114
0115
0116 EcalTrigTowerDetId trigTower = theTriggerMap->towerOf(xtal);
0117 assert(trigTower.rawId() != 0);
0118
0119
0120
0121 combineFlags(eeRuInterest(xtal), getTowerInterest(trigTower));
0122
0123 assert(0 <= eeRuInterest(xtal) && eeRuInterest(xtal) <= 0x7);
0124
0125 }
0126 }
0127 }
0128 #endif
0129 }
0130
0131 EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getCrystalInterest(const EBDetId& ebDetId) const {
0132 EcalTrigTowerDetId thisTower = theTriggerMap->towerOf(ebDetId);
0133 return getTowerInterest(thisTower);
0134 }
0135
0136 EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getCrystalInterest(const EEDetId& eeDetId) const {
0137
0138
0139
0140
0141 return const_cast<EcalSelectiveReadout*>(this)->eeRuInterest(eeDetId);
0142 }
0143
0144 EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getSuperCrystalInterest(const EcalScDetId& scDetId) const {
0145 return const_cast<EcalSelectiveReadout*>(this)->eeRuInterest(scDetId);
0146 }
0147
0148 EcalSelectiveReadout::towerInterest_t& EcalSelectiveReadout::eeRuInterest(const EEDetId& eeDetId) {
0149 const EcalElectronicsId& id = theElecMap->getElectronicsId(eeDetId);
0150 const int iZ0 = id.zside() > 0 ? 1 : 0;
0151 const int iDcc0 = id.dccId() - 1;
0152 const int iDccPhi0 = (iDcc0 < 9) ? iDcc0 : (iDcc0 - 45);
0153 const int iDccCh0 = id.towerId() - 1;
0154 assert(0 <= iDccPhi0 && iDccPhi0 < nDccPerEe);
0155 assert(0 <= iDccCh0 && iDccCh0 < maxDccChs);
0156
0157 assert(eeRuInterest_[iZ0][iDccPhi0][iDccCh0] == UNKNOWN ||
0158 (0 <= eeRuInterest_[iZ0][iDccPhi0][iDccCh0] && eeRuInterest_[iZ0][iDccPhi0][iDccCh0] <= 7));
0159
0160 return eeRuInterest_[iZ0][iDccPhi0][iDccCh0];
0161 }
0162
0163 EcalSelectiveReadout::towerInterest_t& EcalSelectiveReadout::eeRuInterest(const EcalScDetId& scDetId) {
0164 std::pair<int, int> dccAndDccCh = theElecMap->getDCCandSC(scDetId);
0165 const int iZ0 = (scDetId.zside() > 0) ? 1 : 0;
0166 const int iDcc0 = dccAndDccCh.first - 1;
0167 const int iDccPhi0 = (iDcc0 < 9) ? iDcc0 : (iDcc0 - 45);
0168 const int iDccCh0 = dccAndDccCh.second - 1;
0169 assert(0 <= iDccPhi0 && iDccPhi0 <= nDccPerEe);
0170 assert(0 <= iDccCh0 && iDccCh0 <= maxDccChs);
0171
0172 assert(-1 <= eeRuInterest_[iZ0][iDccPhi0][iDccCh0] && eeRuInterest_[iZ0][iDccPhi0][iDccCh0] <= 7);
0173
0174 return eeRuInterest_[iZ0][iDccPhi0][iDccCh0];
0175 }
0176
0177 EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getTowerInterest(const EcalTrigTowerDetId& tower) const {
0178
0179 int iEta = tower.ieta() < 0 ? tower.ieta() + nTriggerTowersInEta / 2 : tower.ieta() + nTriggerTowersInEta / 2 - 1;
0180 int iPhi = tower.iphi() - 1;
0181
0182 assert(-1 <= towerInterest[iEta][iPhi] && int(towerInterest[iEta][iPhi]) < 8);
0183
0184 return towerInterest[iEta][iPhi];
0185 }
0186
0187 void EcalSelectiveReadout::classifyTriggerTowers(const ttFlag_t ttFlags[nTriggerTowersInEta][nTriggerTowersInPhi]) {
0188
0189 for (int iEta = 0; iEta < (int)nTriggerTowersInEta; ++iEta) {
0190 for (int iPhi = 0; iPhi < (int)nTriggerTowersInPhi; ++iPhi) {
0191 towerInterest[iEta][iPhi] = LOWINTEREST;
0192 }
0193 }
0194
0195 for (int iEta = 0; iEta < (int)nTriggerTowersInEta; ++iEta) {
0196 for (int iPhi = 0; iPhi < (int)nTriggerTowersInPhi; ++iPhi) {
0197
0198 towerInterest[iEta][iPhi] = (towerInterest_t)(towerInterest[iEta][iPhi] | (ttFlags[iEta][iPhi] & FORCED_MASK));
0199 if ((ttFlags[iEta][iPhi] & ~FORCED_MASK) == TTF_HIGH_INTEREST) {
0200
0201 towerInterest[iEta][iPhi] = CENTER;
0202
0203 for (int iEtaNeigh = std::max<int>(0, iEta - dEta);
0204 iEtaNeigh <= std::min<int>(nTriggerTowersInEta - 1, iEta + dEta);
0205 ++iEtaNeigh) {
0206 for (int iPhiNeigh = iPhi - dPhi; iPhiNeigh <= iPhi + dPhi; ++iPhiNeigh) {
0207
0208
0209 int iPhiNeigh_ = iPhiNeigh % (int)nTriggerTowersInPhi;
0210 if (iPhiNeigh_ < 0) {
0211 iPhiNeigh_ += nTriggerTowersInPhi;
0212 }
0213 combineFlags(towerInterest[iEtaNeigh][iPhiNeigh_], NEIGHBOUR);
0214 }
0215 }
0216 } else if ((ttFlags[iEta][iPhi] & ~FORCED_MASK) == TTF_MID_INTEREST) {
0217 combineFlags(towerInterest[iEta][iPhi], SINGLE);
0218 }
0219 }
0220 }
0221
0222
0223
0224
0225
0226 const size_t innerEtas[] = {0, 1, nTriggerTowersInEta - 2, nTriggerTowersInEta - 1};
0227 for (size_t i = 0; i < 4; ++i) {
0228 size_t iEta = innerEtas[i];
0229 for (size_t iPhi = 0; iPhi < nTriggerTowersInPhi; iPhi += 2) {
0230 const towerInterest_t srf = std::max(towerInterest[iEta][iPhi], towerInterest[iEta][iPhi + 1]);
0231 towerInterest[iEta][iPhi] = srf;
0232 towerInterest[iEta][iPhi + 1] = srf;
0233 }
0234 }
0235 }
0236
0237 void EcalSelectiveReadout::printHeader(std::ostream& os) const {
0238 os << "#SRP flag map\n#\n"
0239 "# +-->Phi/Y "
0240 << srpFlagMarker[0]
0241 << ": low interest\n"
0242 "# | "
0243 << srpFlagMarker[1]
0244 << ": single\n"
0245 "# | "
0246 << srpFlagMarker[2]
0247 << ": neighbour\n"
0248 "# V Eta/X "
0249 << srpFlagMarker[3]
0250 << ": center\n"
0251 "# "
0252 << srpFlagMarker[4]
0253 << ": forced low interest\n"
0254 "# "
0255 << srpFlagMarker[5]
0256 << ": forced single\n"
0257 "# "
0258 << srpFlagMarker[6]
0259 << ": forced neighbout\n"
0260 "# "
0261 << srpFlagMarker[7]
0262 << ": forced center\n"
0263 "#\n";
0264 }
0265
0266 void EcalSelectiveReadout::print(std::ostream& os) const {
0267
0268 printEndcap(0, os);
0269
0270
0271 printBarrel(os);
0272
0273
0274 printEndcap(1, os);
0275 }
0276
0277 void EcalSelectiveReadout::printBarrel(std::ostream& os) const {
0278 for (size_t iEta = nEndcapTriggerTowersInEta; iEta < nEndcapTriggerTowersInEta + nBarrelTriggerTowersInEta; ++iEta) {
0279 for (size_t iPhi = 0; iPhi < nTriggerTowersInPhi; ++iPhi) {
0280 towerInterest_t srFlag = towerInterest[iEta][iPhi];
0281 os << srpFlagMarker[srFlag];
0282 }
0283 os << "\n";
0284 }
0285 }
0286
0287 void EcalSelectiveReadout::printEndcap(int endcap, std::ostream& os) const {
0288 for (size_t iX = 0; iX < nSupercrystalXBins; ++iX) {
0289 for (size_t iY = 0; iY < nSupercrystalYBins; ++iY) {
0290 towerInterest_t srFlag;
0291 char c;
0292 if (!EcalScDetId::validDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1)) {
0293
0294 c = ' ';
0295 } else {
0296 srFlag = getSuperCrystalInterest(
0297 EcalScDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1));
0298 c = srFlag == UNKNOWN ? '?' : srpFlagMarker[srFlag];
0299 }
0300 os << c;
0301 }
0302 os << "\n";
0303 }
0304 }
0305
0306 std::ostream& operator<<(std::ostream& os, const EcalSelectiveReadout& selectiveReadout) {
0307 selectiveReadout.print(os);
0308 return os;
0309 }
0310
0311 void EcalSelectiveReadout::printDccChMap(std::ostream& os) const {
0312 for (int i = -1; i <= 68; ++i) {
0313 if ((i + 1) % 10 == 0)
0314 os << "//";
0315 os << std::setw(2) << i << ": " << (char)('0' + i);
0316 if (i % 10 == 9)
0317 os << "\n";
0318 else
0319 os << " ";
0320 }
0321
0322 os << "\n";
0323
0324 for (int endcap = 0; endcap < 2; ++endcap) {
0325 os << "Sc2DCCch0: " << (endcap ? "EE+" : "EE-") << "\n";
0326 for (size_t iY = 0; iY < nSupercrystalYBins; ++iY) {
0327 os << "Sc2DCCch0: ";
0328 for (size_t iX = 0; iX < nSupercrystalXBins; ++iX) {
0329
0330 if (!EcalScDetId::validDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1)) {
0331
0332 os << (char)('0' - 1);
0333 } else {
0334
0335 os << (char)('0' + (theElecMap->getDCCandSC(EcalScDetId(iX + 1, iY + 1, endcap > 0 ? 1 : -1)).second - 1));
0336 }
0337 }
0338 os << "\n";
0339 }
0340 os << "\n";
0341 }
0342 os << "\n";
0343 }