L1GctHardwareJetFinder.cc

Go to the documentation of this file.

#include "L1Trigger/GlobalCaloTrigger/interface/L1GctHardwareJetFinder.h"

//DEFINE STATICS
const unsigned int L1GctHardwareJetFinder::N_COLS = 2;
const unsigned int L1GctHardwareJetFinder::CENTRAL_COL0 = 0;
const unsigned int L1GctHardwareJetFinder::MAX_REGIONS_IN =
    (((L1CaloRegionDetId::N_ETA) / 2) + N_EXTRA_REGIONS_ETA00) * L1GctHardwareJetFinder::N_COLS;

L1GctHardwareJetFinder::L1GctHardwareJetFinder(int id)
    : L1GctJetFinderBase(id),
      m_localMaxima(MAX_JETS_OUT),
      m_clusters(MAX_JETS_OUT),
      m_numberOfClusters(0),
      m_localMax00(2),
      m_cluster00(2) {
  this->reset();
  // Initialise parameters for Region input calculations in the
  // derived class so we get the right values of constants.
  static const unsigned NPHI = L1CaloRegionDetId::N_PHI;
  m_minColThisJf = (NPHI + m_id * 2 - CENTRAL_COL0) % NPHI;
}

L1GctHardwareJetFinder::~L1GctHardwareJetFinder() {}

std::ostream& operator<<(std::ostream& os, const L1GctHardwareJetFinder& algo) {
  os << "===L1GctHardwareJetFinder===" << std::endl;
  const L1GctJetFinderBase* temp = &algo;
  os << *temp;
  return os;
}

void L1GctHardwareJetFinder::reset() { L1GctJetFinderBase::reset(); }

void L1GctHardwareJetFinder::fetchInput() {
  if (setupOk()) {
    findProtoJets();
  }
}

void L1GctHardwareJetFinder::process() {
  if (setupOk()) {
    fetchProtoJetsFromNeighbour(TOPBOT);
    findJets();
    sortJets();
    doEnergySums();
  }
}


void L1GctHardwareJetFinder::findProtoJets() {
  findLocalMaxima();
  findProtoClusters();
  convertClustersToProtoJets();
}

void L1GctHardwareJetFinder::findJets() {
  findFinalClusters();
  convertClustersToOutputJets();
}

//  Find the local et maxima in the 2x11 array of regions
void L1GctHardwareJetFinder::findLocalMaxima() {
  m_localMaxima.clear();
  m_localMaxima.resize(MAX_JETS_OUT);
  m_localMax00.clear();
  m_localMax00.resize(2);

  UShort jetNum = 0;  //holds the number of jets currently found
  UShort centreIndex = COL_OFFSET * this->centralCol0();
  for (UShort column = 0; column < 2; ++column)  //Find jets in the central search region
  {
    // The input regions include two extra bins on the other side of eta=0. This allows "seamless"
    // jetfinding across the eta=0 boundary. We skip the first input region in each row. We perform
    // the full pre-clustering on the next region but store the resulting clusters separately
    // from the main list of output pre-clusters - they will be used in the final cluster stage to
    // make sure we do not produce jets in adjacent regions on opposite sides of eta=0.
    ++centreIndex;
    for (UShort row = 1; row < COL_OFFSET; ++row) {
      // Here's the array of greater-than and greater-or-equal tests
      // to ensure each localMaximum appears once and only once in the list
      // It is different for forward and backward eta.
      unsigned JET_THRESHOLD = ((row > m_EtaBoundry) ? m_FwdJetSeed : m_CenJetSeed);
      bool localMax = !m_inputRegions.at(centreIndex).empty() && (m_inputRegions.at(centreIndex).et() >= JET_THRESHOLD);
      if (m_positiveEtaWheel) {  // Forward eta
        localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - 1).et());
        if (row < (COL_OFFSET - 1)) {
          localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + 1).et());
        }
        if (column == 0) {
          localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + COL_OFFSET).et());
          localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + COL_OFFSET - 1).et());
          if (row < (COL_OFFSET - 1)) {
            localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex + COL_OFFSET + 1).et());
          }
        } else {
          localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - COL_OFFSET).et());
          localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - COL_OFFSET - 1).et());
          if (row < (COL_OFFSET - 1)) {
            localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex - COL_OFFSET + 1).et());
          }
        }
      } else {  // Backward eta
        localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - 1).et());
        if (row < (COL_OFFSET - 1)) {
          localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + 1).et());
        }
        if (column == 0) {
          localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + COL_OFFSET).et());
          localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + COL_OFFSET - 1).et());
          if (row < (COL_OFFSET - 1)) {
            localMax &= (m_inputRegions.at(centreIndex).et() >= m_inputRegions.at(centreIndex + COL_OFFSET + 1).et());
          }
        } else {
          localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - COL_OFFSET).et());
          localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - COL_OFFSET - 1).et());
          if (row < (COL_OFFSET - 1)) {
            localMax &= (m_inputRegions.at(centreIndex).et() > m_inputRegions.at(centreIndex - COL_OFFSET + 1).et());
          }
        }
      }
      if (localMax) {
        if (row > 1) {
          if (jetNum < MAX_JETS_OUT) {
            m_localMaxima.at(jetNum) = m_inputRegions.at(centreIndex);
            ++jetNum;
          }
        }
        // Treat row 1 as a separate case. It's not required for jetfinding but
        // is used for vetoing of jets double counted across the eta=0 boundary
        else {
          unsigned phi = m_inputRegions.at(centreIndex).rctPhi();
          m_localMax00.at(phi) = m_inputRegions.at(centreIndex);
        }
      }
      ++centreIndex;
    }
  }

  m_numberOfClusters = jetNum;
}

//  For each local maximum, find the cluster et in a 2x3 region.
//  The logic ensures that a given region et cannot be used in more than one cluster.
//  The sorting of the local maxima ensures the highest et maximum has priority.
void L1GctHardwareJetFinder::findProtoClusters() {
  m_clusters.clear();
  m_clusters.resize(MAX_JETS_OUT);
  m_cluster00.clear();
  m_cluster00.resize(2);

  RegionsVector topJets(MAX_JETS_OUT), botJets(MAX_JETS_OUT);
  std::vector<unsigned> topJetsPosition(MAX_JETS_OUT), botJetsPosition(MAX_JETS_OUT);
  unsigned numberOfTopJets = 0, numberOfBotJets = 0;

  // Loop over local maxima
  for (unsigned j = 0; j < m_numberOfClusters; ++j) {
    // Make a proto-jet cluster
    L1GctRegion temp = makeProtoJet(m_localMaxima.at(j));

    if (m_localMaxima.at(j).rctPhi() == 0) {
      // Store "top edge" jets
      topJets.at(numberOfTopJets) = temp;
      topJetsPosition.at(numberOfTopJets) = 0;
      for (unsigned k = 0; k < numberOfTopJets; ++k) {
        if (topJets.at(numberOfTopJets).et() >= topJets.at(k).et()) {
          ++topJetsPosition.at(k);
        }
        if (topJets.at(numberOfTopJets).et() <= topJets.at(k).et()) {
          ++topJetsPosition.at(numberOfTopJets);
        }
      }
      ++numberOfTopJets;
    } else {
      // Store "bottom edge" jets
      botJets.at(numberOfBotJets) = temp;
      botJetsPosition.at(numberOfBotJets) = 0;
      for (unsigned k = 0; k < numberOfBotJets; ++k) {
        if (botJets.at(numberOfBotJets).et() >= botJets.at(k).et()) {
          ++botJetsPosition.at(k);
        }
        if (botJets.at(numberOfBotJets).et() <= botJets.at(k).et()) {
          ++botJetsPosition.at(numberOfBotJets);
        }
      }
      ++numberOfBotJets;
    }
  }
  // Now we've found all the proto-jets, copy the best ones to the output array
  //
  // We fill the first half of the array with "bottom jets"
  // and the remainder with "top jets". For cases where
  // we have found too many jets in one phi column,
  // we keep those with the highest Et.
  static const unsigned int MAX_TOPBOT_JETS = MAX_JETS_OUT / 2;
  unsigned pos = 0;
  for (unsigned j = 0; j < numberOfBotJets; ++j) {
    if (botJetsPosition.at(j) < MAX_TOPBOT_JETS) {
      m_clusters.at(pos++) = botJets.at(j);
    }
  }
  pos = MAX_TOPBOT_JETS;
  for (unsigned j = 0; j < numberOfTopJets; ++j) {
    if (topJetsPosition.at(j) < MAX_TOPBOT_JETS) {
      m_clusters.at(pos++) = topJets.at(j);
    }
  }
  // Finally, deal with eta00 maxima
  if (!m_localMax00.at(0).empty())
    m_cluster00.at(0) = makeProtoJet(m_localMax00.at(0));
  if (!m_localMax00.at(1).empty())
    m_cluster00.at(1) = makeProtoJet(m_localMax00.at(1));
}

L1GctRegion L1GctHardwareJetFinder::makeProtoJet(L1GctRegion localMax) {
  unsigned eta = localMax.gctEta();
  unsigned phi = localMax.gctPhi();
  int16_t bx = localMax.bx();

  unsigned localEta = localMax.rctEta();
  unsigned localPhi = localMax.rctPhi();

  unsigned etCluster = 0;
  bool ovrFlowOr = false;
  bool tauVetoOr = false;
  unsigned rgnsAboveIsoThreshold = 0;

  // check for row00
  const unsigned midEta = (L1CaloRegionDetId::N_ETA) / 2;
  bool wrongEtaWheel = ((!m_positiveEtaWheel) && (eta >= midEta)) || ((m_positiveEtaWheel) && (eta < midEta));

  // Which rows are we looking over?
  unsigned rowStart, rowEnd, rowMid;
  static const unsigned row0 = N_EXTRA_REGIONS_ETA00 - 1;
  if (wrongEtaWheel) {
    if (localEta > row0 - 1) {
      rowStart = 0;
      rowMid = 0;
    } else {
      rowStart = row0 - 1 - localEta;
      rowMid = rowStart + 1;
    }
    if (localEta > row0 + 2) {  // Shouldn't happen, but big problems if it does
      rowEnd = 0;
    } else {
      rowEnd = row0 + 2 - localEta;
    }
  } else {
    rowStart = row0 + localEta;
    rowMid = rowStart + 1;
    if (localEta < COL_OFFSET - row0 - 2) {
      rowEnd = rowStart + 3;
    } else {
      rowEnd = COL_OFFSET;
    }
  }

  for (unsigned row = rowStart; row < rowEnd; ++row) {
    for (unsigned column = 0; column < 2; ++column) {
      unsigned index = column * COL_OFFSET + row;
      etCluster += m_inputRegions.at(index).et();
      ovrFlowOr |= m_inputRegions.at(index).overFlow();
      // Distinguish between central and tau-flagged jets. Two versions of the algorithm.
      if (m_useImprovedTauAlgo) {
        //===========================================================================================
        // "Old" version of improved tau algorithm tests the tau veto for the central region always
        //  if ((row==(localEta+N_EXTRA_REGIONS_ETA00)) && (column==localPhi)) {
        //    // central region - check the tau veto
        //    tauVetoOr |= m_inputRegions.at(index).tauVeto();
        //  } else {
        //    // other regions - check the tau veto if required
        //    if (!m_ignoreTauVetoBitsForIsolation) {
        //      tauVetoOr |= m_inputRegions.at(index).tauVeto();
        //    }
        //    // check the region energy against the isolation threshold
        //    if (m_inputRegions.at(index).et() >= m_tauIsolationThreshold) {
        //      rgnsAboveIsoThreshold++;
        //    }
        //  }
        //===========================================================================================

        // In the hardware, the ignoreTauVetoBitsForIsolation switch ignores all the veto bits,
        // including the one for the central region.
        if (!((row == rowMid) && (column == localPhi))) {
          // non-central region - check the region energy against the isolation threshold
          if (m_inputRegions.at(index).et() >= m_tauIsolationThreshold) {
            rgnsAboveIsoThreshold++;
          }
        }
        // all regions - check the tau veto if required
        if (!m_ignoreTauVetoBitsForIsolation) {
          tauVetoOr |= m_inputRegions.at(index).tauVeto();
        }
        // End of improved tau algorithm
      } else {
        // Original tau algorithm
        tauVetoOr |= m_inputRegions.at(index).tauVeto();
      }
    }
  }
  // Encode the number of towers over threshold for the isolated tau algorithm
  bool tauFeatureBit = false;
  if (m_useImprovedTauAlgo) {
    tauVetoOr |= (rgnsAboveIsoThreshold > 1);
    tauFeatureBit |= (rgnsAboveIsoThreshold == 1);
  }

  L1GctRegion temp(L1GctRegion::makeProtoJetRegion(etCluster, ovrFlowOr, tauVetoOr, tauFeatureBit, eta, phi, bx));
  return temp;
}

void L1GctHardwareJetFinder::findFinalClusters() {
  m_clusters.clear();
  m_clusters.resize(MAX_JETS_OUT);

  // Loop over proto-jets received from neighbours.
  // Form a jet to send to the output if there is no proto-jet nearby in the
  // list of jets found locally. If local jets are found nearby, form a jet
  // if the received jet has higher Et than any one of the local ones.
  for (unsigned j = 0; j < MAX_JETS_OUT; ++j) {
    unsigned et0 = m_rcvdProtoJets.at(j).et();
    unsigned localEta0 = m_rcvdProtoJets.at(j).rctEta();
    unsigned localPhi0 = m_rcvdProtoJets.at(j).rctPhi();
    unsigned JET_THRESHOLD = ((localEta0 >= m_EtaBoundry) ? m_FwdJetSeed : m_CenJetSeed);
    if (et0 >= JET_THRESHOLD) {
      bool storeJet = false;
      bool isolated = true;
      // eta00 boundary check/veto
      if (localEta0 == 0) {
        unsigned neighbourEt = m_cluster00.at(1 - localPhi0).et();
        isolated &= et0 >= neighbourEt;
      }
      // If the jet is NOT vetoed, look at the jets found locally (m_keptProtoJets).
      // We accept the jet if there are no local jets nearby, or if the local jet
      // (there should be no more than one) has lower Et.
      if (isolated) {
        for (unsigned k = 0; k < MAX_JETS_OUT; ++k) {
          unsigned et1 = m_keptProtoJets.at(k).et();
          unsigned localEta1 = m_keptProtoJets.at(k).rctEta();
          unsigned localPhi1 = m_keptProtoJets.at(k).rctPhi();
          if (et1 > 0) {
            bool distantJet = ((localPhi0 == localPhi1) || (localEta1 > localEta0 + 1) || (localEta0 > localEta1 + 1));

            isolated &= distantJet;
            storeJet |= !distantJet && ((et0 > et1) || ((et0 == et1) && localPhi0 == 1));
          }
        }
      }

      storeJet |= isolated;

      if (storeJet) {
        // Start with the et sum, tau veto and overflow flags of the protoJet (2x3 regions)
        unsigned etCluster = et0;
        bool ovrFlowOr = m_rcvdProtoJets.at(j).overFlow();
        bool tauVetoOr = m_rcvdProtoJets.at(j).tauVeto();
        unsigned rgnsAboveIsoThreshold = (m_rcvdProtoJets.at(j).featureBit0() ? 1 : 0);

        // Combine with the corresponding regions from
        // the local array to make a 3x3 jet cluster
        unsigned column = 1 - localPhi0;
        // Which rows are we looking over?
        unsigned rowStart, rowEnd;
        static const unsigned row0 = N_EXTRA_REGIONS_ETA00 - 1;
        rowStart = row0 + localEta0;
        if (localEta0 < COL_OFFSET - row0 - 2) {
          rowEnd = rowStart + 3;
        } else {
          rowEnd = COL_OFFSET;
        }
        unsigned index = COL_OFFSET * (this->centralCol0() + column) + rowStart;
        for (unsigned row = rowStart; row < rowEnd; ++row) {
          etCluster += m_inputRegions.at(index).et();
          ovrFlowOr |= m_inputRegions.at(index).overFlow();
          if (m_useImprovedTauAlgo) {
            if (!m_ignoreTauVetoBitsForIsolation) {
              tauVetoOr |= m_inputRegions.at(index).tauVeto();
            }
            // check the region energy against the isolation threshold
            if (m_inputRegions.at(index).et() >= m_tauIsolationThreshold) {
              rgnsAboveIsoThreshold++;
            }
          } else {
            tauVetoOr |= m_inputRegions.at(index).tauVeto();
          }

          ++index;
        }

        // Store the new jet
        unsigned eta = m_rcvdProtoJets.at(j).gctEta();
        unsigned phi = m_rcvdProtoJets.at(j).gctPhi();
        int16_t bx = m_rcvdProtoJets.at(j).bx();

        // Use the number of towers over threshold for the isolated tau algorithm
        if (m_useImprovedTauAlgo) {
          tauVetoOr |= (rgnsAboveIsoThreshold > 1);
        }

        L1GctRegion temp(L1GctRegion::makeFinalJetRegion(etCluster, ovrFlowOr, tauVetoOr, eta, phi, bx));
        m_clusters.at(j) = temp;
      }
    }
  }
}

void L1GctHardwareJetFinder::convertClustersToProtoJets() {
  for (unsigned j = 0; j < MAX_JETS_OUT; ++j) {
    bool isForward = (m_clusters.at(j).rctEta() >= m_EtaBoundry);
    unsigned JET_THRESHOLD = (isForward ? m_FwdJetSeed : m_CenJetSeed);
    if (m_clusters.at(j).et() >= JET_THRESHOLD) {
      m_keptProtoJets.at(j) = m_clusters.at(j);
      m_sentProtoJets.at(j) = m_clusters.at(j);
    }
  }
}

void L1GctHardwareJetFinder::convertClustersToOutputJets() {
  for (unsigned j = 0; j < MAX_JETS_OUT; ++j) {
    bool isForward = (m_clusters.at(j).rctEta() >= m_EtaBoundry);
    unsigned JET_THRESHOLD = (isForward ? m_FwdJetSeed : m_CenJetSeed);
    if (m_clusters.at(j).et() >= JET_THRESHOLD) {
      L1GctJet temp(m_clusters.at(j).et(),
                    m_clusters.at(j).gctEta(),
                    m_clusters.at(j).gctPhi(),
                    m_clusters.at(j).overFlow(),
                    isForward,
                    m_clusters.at(j).tauVeto(),
                    m_clusters.at(j).bx());
      m_outputJets.at(j) = temp;
    }
  }
}