EcalSelectiveReadout.cc

Go to the documentation of this file.

//emacs settings:-*- mode: c++; c-basic-offset: 2; indent-tabs-mode: nil -*-"
/*
 */

#include "SimCalorimetry/EcalSelectiveReadoutAlgos/interface/EcalSelectiveReadout.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <string>
#include <iomanip>
#include <cassert>
#include <atomic>
//#include <iostream> //for debugging

using std::vector;

const char EcalSelectiveReadout::srpFlagMarker[] = {'.', 'S', 'N', 'C', '4', '5', '6', '7'};

EcalSelectiveReadout::EcalSelectiveReadout(int dEta_, int dPhi_)
    : theTriggerMap(nullptr), theElecMap(nullptr), dEta(dEta_), dPhi(dPhi_) {}

void EcalSelectiveReadout::resetEeRuInterest() {
  //init superCrystalInterest (sets all elts to 'UNKNOWN'):
  for (size_t iCap = 0; iCap < nEndcaps; ++iCap) {
    for (int iDccPhi = 0; iDccPhi < nDccPerEe; ++iDccPhi) {
      for (int iDccCh = 0; iDccCh < maxDccChs; ++iDccCh) {
        eeRuInterest_[iCap][iDccPhi][iDccCh] = UNKNOWN;
      }
    }
  }
}

void EcalSelectiveReadout::runSelectiveReadout0(const ttFlag_t ttFlags[nTriggerTowersInEta][nTriggerTowersInPhi]) {
  //printDccChMap(std::cout);

  //classifies the trigger towers (single,center,neighbor,low interest)
  classifyTriggerTowers(ttFlags);

  //count number of TT in each interest class for debugging display
  int nTriggerTowerE[] = {0, 0, 0, 0, 0, 0, 0, 0};
  int nTriggerTowerB[] = {0, 0, 0, 0, 0, 0, 0, 0};
  static std::atomic<int> ncall{0};
  if (ncall < 10) {
    ++ncall;
    for (size_t iPhi = 0; iPhi < nTriggerTowersInPhi; ++iPhi) {
      for (size_t iEta = 0; iEta < nTriggerTowersInEta; ++iEta) {
        if (iEta < nEndcapTriggerTowersInEta || iEta >= nBarrelTriggerTowersInEta + nEndcapTriggerTowersInEta) {
          //in endcaps
          ++nTriggerTowerE[towerInterest[iEta][iPhi]];
        } else {  //in barrel
          ++nTriggerTowerB[towerInterest[iEta][iPhi]];
        }

        assert(towerInterest[iEta][iPhi] >= 0 && towerInterest[iEta][iPhi] <= 0x7);
      }
    }
    edm::LogInfo("EcalSelectiveReadout") << "without forced bit + with forced bit set:\n"
                                         << nTriggerTowerB[LOWINTEREST] << " + "
                                         << nTriggerTowerB[LOWINTEREST | FORCED_MASK]
                                         << " low interest TT(s) in barrel\n"
                                         << nTriggerTowerB[SINGLE] << " + " << nTriggerTowerB[SINGLE | FORCED_MASK]
                                         << " single TT(s) in barrel\n"
                                         << nTriggerTowerB[NEIGHBOUR] << " + "
                                         << nTriggerTowerB[NEIGHBOUR | FORCED_MASK]
                                         << " neighbor interest TT(s) in barrel\n"
                                         << nTriggerTowerB[CENTER] << " + " << nTriggerTowerB[CENTER | FORCED_MASK]
                                         << " centre interest TT(s) in barrel\n"
                                         << nTriggerTowerE[LOWINTEREST] << " + "
                                         << nTriggerTowerE[LOWINTEREST | FORCED_MASK]
                                         << " low interest TT(s) in endcap\n"
                                         << nTriggerTowerE[SINGLE] << " + " << nTriggerTowerE[SINGLE | FORCED_MASK]
                                         << " single TT(s) in endcap\n"
                                         << nTriggerTowerE[NEIGHBOUR] << " + "
                                         << nTriggerTowerE[NEIGHBOUR | FORCED_MASK] << " neighbor TT(s) in endcap\n"
                                         << nTriggerTowerE[CENTER] << " + " << nTriggerTowerE[CENTER | FORCED_MASK]
                                         << " center TT(s) in endcap\n";
  }
  //end TT interest class composition debugging display

  //For the endcap the TT classification must be mapped to the SC:
  resetEeRuInterest();

#ifndef ECALSELECTIVEREADOUT_NOGEOM
  const std::vector<DetId>& endcapDetIds = theGeometry->getValidDetIds(DetId::Ecal, EcalEndcap);
  for (std::vector<DetId>::const_iterator eeDetIdItr = endcapDetIds.begin(); eeDetIdItr != endcapDetIds.end();
       ++eeDetIdItr) {
    // for each superCrystal, the interest is the highest interest
    // of any trigger tower associated with at least one crystal from this SC.
    // The forced bit must be set if the flag of one of these trigger towers has
    // the forced bit set.
    EcalTrigTowerDetId trigTower = theTriggerMap->towerOf(*eeDetIdItr);
    assert(trigTower.rawId() != 0);
    EEDetId eeDetId(*eeDetIdItr);
    int iz = (eeDetId.zside() > 0) ? 1 : 0;
    //Following statement will set properly the actual 2-bit flag value
    //and the forced bit: TTF forced bit is propagated to every RU that
    //overlaps with the corresponding TT.
    combineFlags(eeRuInterest(eeDetId), getTowerInterest(trigTower));
  }
#else   //ECALSELECTIVEREADOUT_NOGEOM defined
  EEDetId xtal;
  for (int iZ0 = 0; iZ0 < 2; ++iZ0) {  //0->EE-, 1->EE+
    for (unsigned iX0 = 0; iX0 < nEndcapXBins; ++iX0) {
      for (unsigned iY0 = 0; iY0 < nEndcapYBins; ++iY0) {
        if (!(xtal.validDetId(iX0 + 1, iY0 + 1, (iZ0 > 0 ? 1 : -1)))) {
          continue;
        }
        xtal = EEDetId(iX0 + 1, iY0 + 1, (iZ0 > 0 ? 1 : -1));
        //works around a EEDetId bug. To remove once the bug fixed.
        if (39 <= iX0 && iX0 <= 60 && 45 <= iY0 && iY0 <= 54) {
          continue;
        }
        // for each superCrystal, the interest is the highest interest
        // of any trigger tower associated with any crystal in this SC
        EcalTrigTowerDetId trigTower = theTriggerMap->towerOf(xtal);
        assert(trigTower.rawId() != 0);
        //Following statement will set properly the actual 2-bit flag value
        //and the forced bit: TTF forced bit is propagated to every RU that
        //overlaps with the corresponding TT.
        combineFlags(eeRuInterest(xtal), getTowerInterest(trigTower));

        assert(0 <= eeRuInterest(xtal) && eeRuInterest(xtal) <= 0x7);

      }  //next iY0
    }    //next iX0
  }      //next iZ0
#endif  //ECALSELECTIVEREADOUT_NOGEOM not defined
}

EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getCrystalInterest(const EBDetId& ebDetId) const {
  EcalTrigTowerDetId thisTower = theTriggerMap->towerOf(ebDetId);
  return getTowerInterest(thisTower);
}

EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getCrystalInterest(const EEDetId& eeDetId) const {
  //   int iz = (eeDetId.zside() > 0) ? 1 : 0;
  //   int superCrystalX = (eeDetId.ix()-1) / 5;
  //   int superCrystalY = (eeDetId.iy()-1) / 5;
  //   return supercrystalInterest[iz][superCrystalX][superCrystalY];
  return const_cast<EcalSelectiveReadout*>(this)->eeRuInterest(eeDetId);
}

EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getSuperCrystalInterest(const EcalScDetId& scDetId) const {
  return const_cast<EcalSelectiveReadout*>(this)->eeRuInterest(scDetId);
}

EcalSelectiveReadout::towerInterest_t& EcalSelectiveReadout::eeRuInterest(const EEDetId& eeDetId) {
  const EcalElectronicsId& id = theElecMap->getElectronicsId(eeDetId);
  const int iZ0 = id.zside() > 0 ? 1 : 0;
  const int iDcc0 = id.dccId() - 1;
  const int iDccPhi0 = (iDcc0 < 9) ? iDcc0 : (iDcc0 - 45);
  const int iDccCh0 = id.towerId() - 1;
  assert(0 <= iDccPhi0 && iDccPhi0 < nDccPerEe);
  assert(0 <= iDccCh0 && iDccCh0 < maxDccChs);

  assert(eeRuInterest_[iZ0][iDccPhi0][iDccCh0] == UNKNOWN ||
         (0 <= eeRuInterest_[iZ0][iDccPhi0][iDccCh0] && eeRuInterest_[iZ0][iDccPhi0][iDccCh0] <= 7));

  return eeRuInterest_[iZ0][iDccPhi0][iDccCh0];
}

EcalSelectiveReadout::towerInterest_t& EcalSelectiveReadout::eeRuInterest(const EcalScDetId& scDetId) {
  std::pair<int, int> dccAndDccCh = theElecMap->getDCCandSC(scDetId);
  const int iZ0 = (scDetId.zside() > 0) ? 1 : 0;
  const int iDcc0 = dccAndDccCh.first - 1;
  const int iDccPhi0 = (iDcc0 < 9) ? iDcc0 : (iDcc0 - 45);
  const int iDccCh0 = dccAndDccCh.second - 1;
  assert(0 <= iDccPhi0 && iDccPhi0 <= nDccPerEe);
  assert(0 <= iDccCh0 && iDccCh0 <= maxDccChs);

  assert(-1 <= eeRuInterest_[iZ0][iDccPhi0][iDccCh0] && eeRuInterest_[iZ0][iDccPhi0][iDccCh0] <= 7);

  return eeRuInterest_[iZ0][iDccPhi0][iDccCh0];
}

EcalSelectiveReadout::towerInterest_t EcalSelectiveReadout::getTowerInterest(const EcalTrigTowerDetId& tower) const {
  // remember, array indices start at zero
  int iEta = tower.ieta() < 0 ? tower.ieta() + nTriggerTowersInEta / 2 : tower.ieta() + nTriggerTowersInEta / 2 - 1;
  int iPhi = tower.iphi() - 1;

  assert(-1 <= towerInterest[iEta][iPhi] && int(towerInterest[iEta][iPhi]) < 8);

  return towerInterest[iEta][iPhi];
}

void EcalSelectiveReadout::classifyTriggerTowers(const ttFlag_t ttFlags[nTriggerTowersInEta][nTriggerTowersInPhi]) {
  //starts with a all low interest map:
  for (int iEta = 0; iEta < (int)nTriggerTowersInEta; ++iEta) {
    for (int iPhi = 0; iPhi < (int)nTriggerTowersInPhi; ++iPhi) {
      towerInterest[iEta][iPhi] = LOWINTEREST;
    }
  }

  for (int iEta = 0; iEta < (int)nTriggerTowersInEta; ++iEta) {
    for (int iPhi = 0; iPhi < (int)nTriggerTowersInPhi; ++iPhi) {
      //copy forced bit from ttFlags to towerInterests:
      towerInterest[iEta][iPhi] = (towerInterest_t)(towerInterest[iEta][iPhi] | (ttFlags[iEta][iPhi] & FORCED_MASK));
      if ((ttFlags[iEta][iPhi] & ~FORCED_MASK) == TTF_HIGH_INTEREST) {
        //flags this tower as a center tower
        towerInterest[iEta][iPhi] = CENTER;
        //flags the neighbours of this tower
        for (int iEtaNeigh = std::max<int>(0, iEta - dEta);
             iEtaNeigh <= std::min<int>(nTriggerTowersInEta - 1, iEta + dEta);
             ++iEtaNeigh) {
          for (int iPhiNeigh = iPhi - dPhi; iPhiNeigh <= iPhi + dPhi; ++iPhiNeigh) {
            //beware, iPhiNeigh must be moved to [0,72] interval
            //=> %nTriggerTowersInPhi required
            int iPhiNeigh_ = iPhiNeigh % (int)nTriggerTowersInPhi;
            if (iPhiNeigh_ < 0) {
              iPhiNeigh_ += nTriggerTowersInPhi;
            }
            combineFlags(towerInterest[iEtaNeigh][iPhiNeigh_], NEIGHBOUR);
          }
        }
      } else if ((ttFlags[iEta][iPhi] & ~FORCED_MASK) == TTF_MID_INTEREST) {
        combineFlags(towerInterest[iEta][iPhi], SINGLE);
      }
    }
  }

  //dealing with pseudo-TT in the two innest eta-ring of the endcaps
  //=>choose the highest priority  SRF of the 2 pseudo-TT constituting
  //a TT. Note that for S and C, the 2 pseudo-TT must already have the
  //same mask.
  const size_t innerEtas[] = {0, 1, nTriggerTowersInEta - 2, nTriggerTowersInEta - 1};
  for (size_t i = 0; i < 4; ++i) {
    size_t iEta = innerEtas[i];
    for (size_t iPhi = 0; iPhi < nTriggerTowersInPhi; iPhi += 2) {
      const towerInterest_t srf = std::max(towerInterest[iEta][iPhi], towerInterest[iEta][iPhi + 1]);
      towerInterest[iEta][iPhi] = srf;
      towerInterest[iEta][iPhi + 1] = srf;
    }
  }
}

void EcalSelectiveReadout::printHeader(std::ostream& os) const {
  os << "#SRP flag map\n#\n"
        "# +-->Phi/Y "
     << srpFlagMarker[0]
     << ": low interest\n"
        "# |         "
     << srpFlagMarker[1]
     << ": single\n"
        "# |         "
     << srpFlagMarker[2]
     << ": neighbour\n"
        "# V Eta/X   "
     << srpFlagMarker[3]
     << ": center\n"
        "#           "
     << srpFlagMarker[4]
     << ": forced low interest\n"
        "#           "
     << srpFlagMarker[5]
     << ": forced single\n"
        "#           "
     << srpFlagMarker[6]
     << ": forced neighbout\n"
        "#           "
     << srpFlagMarker[7]
     << ": forced center\n"
        "#\n";
}

void EcalSelectiveReadout::print(std::ostream& os) const {
  //EE-
  printEndcap(0, os);

  //EB
  printBarrel(os);

  //EE+
  printEndcap(1, os);
}

void EcalSelectiveReadout::printBarrel(std::ostream& os) const {
  for (size_t iEta = nEndcapTriggerTowersInEta; iEta < nEndcapTriggerTowersInEta + nBarrelTriggerTowersInEta; ++iEta) {
    for (size_t iPhi = 0; iPhi < nTriggerTowersInPhi; ++iPhi) {
      towerInterest_t srFlag = towerInterest[iEta][iPhi];
      os << srpFlagMarker[srFlag];
    }
    os << "\n";  //one phi per line
  }
}

void EcalSelectiveReadout::printEndcap(int endcap, std::ostream& os) const {
  for (size_t iX = 0; iX < nSupercrystalXBins; ++iX) {
    for (size_t iY = 0; iY < nSupercrystalYBins; ++iY) {
      towerInterest_t srFlag;
      char c;
      if (!EcalScDetId::validDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1)) {
        //        srFlag = UNKNOWN;
        c = ' ';
      } else {
        srFlag = getSuperCrystalInterest(
            EcalScDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1));  //supercrystalInterest[endcap][iX][iY];
        c = srFlag == UNKNOWN ? '?' : srpFlagMarker[srFlag];
      }
      os << c;
    }
    os << "\n";  //one Y supercystal column per line
  }              //next supercrystal X-index
}

std::ostream& operator<<(std::ostream& os, const EcalSelectiveReadout& selectiveReadout) {
  selectiveReadout.print(os);
  return os;
}

void EcalSelectiveReadout::printDccChMap(std::ostream& os) const {
  for (int i = -1; i <= 68; ++i) {
    if ((i + 1) % 10 == 0)
      os << "//";
    os << std::setw(2) << i << ": " << (char)('0' + i);
    if (i % 10 == 9)
      os << "\n";
    else
      os << " ";
  }

  os << "\n";

  for (int endcap = 0; endcap < 2; ++endcap) {
    os << "Sc2DCCch0: " << (endcap ? "EE+" : "EE-") << "\n";
    for (size_t iY = 0; iY < nSupercrystalYBins; ++iY) {
      os << "Sc2DCCch0: ";
      for (size_t iX = 0; iX < nSupercrystalXBins; ++iX) {
        //if(iX) os << ",";
        if (!EcalScDetId::validDetId(iX + 1, iY + 1, endcap >= 1 ? 1 : -1)) {
          //os << std::setw(2) << -1;
          os << (char)('0' - 1);
        } else {
          //os << std::setw(2) << theElecMap->getDCCandSC(EcalScDetId(iX+1, iY+1, endcap>0?1:-1)).second-1;
          os << (char)('0' + (theElecMap->getDCCandSC(EcalScDetId(iX + 1, iY + 1, endcap > 0 ? 1 : -1)).second - 1));
        }
      }
      os << "\n";
    }
    os << "\n";
  }
  os << "\n";
}