EcalSelectiveReadoutAlgos/src/EcalSelectiveReadout.cc

0001 //emacs settings:-*- mode: c++; c-basic-offset: 2; indent-tabs-mode: nil -*-"
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 //#include <iostream> //for debugging
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   //init superCrystalInterest (sets all elts to 'UNKNOWN'):
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   //printDccChMap(std::cout);
0036
0037   //classifies the trigger towers (single,center,neighbor,low interest)
0038   classifyTriggerTowers(ttFlags);
0039
0040   //count number of TT in each interest class for debugging display
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           //in endcaps
0050           ++nTriggerTowerE[towerInterest[iEta][iPhi]];
0051         } else {  //in barrel
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   //end TT interest class composition debugging display
0080
0081   //For the endcap the TT classification must be mapped to the SC:
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     // for each superCrystal, the interest is the highest interest
0089     // of any trigger tower associated with at least one crystal from this SC.
0090     // The forced bit must be set if the flag of one of these trigger towers has
0091     // the forced bit set.
0092     EcalTrigTowerDetId trigTower = theTriggerMap->towerOf(*eeDetIdItr);
0093     assert(trigTower.rawId() != 0);
0094     EEDetId eeDetId(*eeDetIdItr);
0095     int iz = (eeDetId.zside() > 0) ? 1 : 0;
0096     //Following statement will set properly the actual 2-bit flag value
0097     //and the forced bit: TTF forced bit is propagated to every RU that
0098     //overlaps with the corresponding TT.
0099     combineFlags(eeRuInterest(eeDetId), getTowerInterest(trigTower));
0100   }
0101 #else   //ECALSELECTIVEREADOUT_NOGEOM defined
0102   EEDetId xtal;
0103   for (int iZ0 = 0; iZ0 < 2; ++iZ0) {  //0->EE-, 1->EE+
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         //works around a EEDetId bug. To remove once the bug fixed.
0111         if (39 <= iX0 && iX0 <= 60 && 45 <= iY0 && iY0 <= 54) {
0112           continue;
0113         }
0114         // for each superCrystal, the interest is the highest interest
0115         // of any trigger tower associated with any crystal in this SC
0116         EcalTrigTowerDetId trigTower = theTriggerMap->towerOf(xtal);
0117         assert(trigTower.rawId() != 0);
0118         //Following statement will set properly the actual 2-bit flag value
0119         //and the forced bit: TTF forced bit is propagated to every RU that
0120         //overlaps with the corresponding TT.
0121         combineFlags(eeRuInterest(xtal), getTowerInterest(trigTower));
0122
0123         assert(0 <= eeRuInterest(xtal) && eeRuInterest(xtal) <= 0x7);
0124
0125       }  //next iY0
0126     }    //next iX0
0127   }      //next iZ0
0128 #endif  //ECALSELECTIVEREADOUT_NOGEOM not defined
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   //   int iz = (eeDetId.zside() > 0) ? 1 : 0;
0138   //   int superCrystalX = (eeDetId.ix()-1) / 5;
0139   //   int superCrystalY = (eeDetId.iy()-1) / 5;
0140   //   return supercrystalInterest[iz][superCrystalX][superCrystalY];
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   // remember, array indices start at zero
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   //starts with a all low interest map:
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       //copy forced bit from ttFlags to towerInterests:
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         //flags this tower as a center tower
0201         towerInterest[iEta][iPhi] = CENTER;
0202         //flags the neighbours of this tower
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             //beware, iPhiNeigh must be moved to [0,72] interval
0208             //=> %nTriggerTowersInPhi required
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   //dealing with pseudo-TT in the two innest eta-ring of the endcaps
0223   //=>choose the highest priority  SRF of the 2 pseudo-TT constituting
0224   //a TT. Note that for S and C, the 2 pseudo-TT must already have the
0225   //same mask.
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   //EE-
0268   printEndcap(0, os);
0269
0270   //EB
0271   printBarrel(os);
0272
0273   //EE+
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";  //one phi per line
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         //        srFlag = UNKNOWN;
0294         c = ' ';
0295       } else {
0296         srFlag = getSuperCrystalInterest(
0297             EcalScDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1));  //supercrystalInterest[endcap][iX][iY];
0298         c = srFlag == UNKNOWN ? '?' : srpFlagMarker[srFlag];
0299       }
0300       os << c;
0301     }
0302     os << "\n";  //one Y supercystal column per line
0303   }              //next supercrystal X-index
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         //if(iX) os << ",";
0330         if (!EcalScDetId::validDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1)) {
0331           //os << std::setw(2) << -1;
0332           os << (char)('0' - 1);
0333         } else {
0334           //os << std::setw(2) << theElecMap->getDCCandSC(EcalScDetId(iX+1, iY+1, endcap>0?1:-1)).second-1;
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 }