/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/L1Trigger/GlobalCaloTrigger/src/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;
    }
  }
}