L1GctTdrJetFinder.cc

Go to the documentation of this file.

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

using namespace std;

//DEFINE STATICS
// *** Note the following definition in terms of COL_OFFSET appears not to work ***
// ***          for some deep C++ reason that I don't understand - GPH          ***
// const unsigned int L1GctTdrJetFinder::MAX_REGIONS_IN = L1GctJetFinderBase::COL_OFFSET*L1GctTdrJetFinder::N_COLS;
// ***                        So - use the following instead                    ***
const unsigned int L1GctTdrJetFinder::MAX_REGIONS_IN = (((L1CaloRegionDetId::N_ETA)/2)+1)*L1GctTdrJetFinder::N_COLS;

const unsigned int L1GctTdrJetFinder::N_COLS = 4;
const unsigned int L1GctTdrJetFinder::CENTRAL_COL0 = 1;

L1GctTdrJetFinder::L1GctTdrJetFinder(int id):
  L1GctJetFinderBase(id)
{
  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;
}

L1GctTdrJetFinder::~L1GctTdrJetFinder()
{
}

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

void L1GctTdrJetFinder::fetchInput()
{
}

void L1GctTdrJetFinder::process() 
{
  if (setupOk()) {
    findJets();
    sortJets();
    doEnergySums();
  }
}


void L1GctTdrJetFinder::findJets() 
{
  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
  {
    //don't include row zero as it is not in the search region
    ++centreIndex;
    for (UShort row = 1; row < COL_OFFSET; ++row)  
    {
      //Determine if we are at end of the HF or not (so need 3*2 window)
      bool hfBoundary = (row == COL_OFFSET-1);
      //Determine if we are at the end of the endcap HCAL regions, so need boundary condition tauveto
      bool heBoundary = (row == COL_OFFSET-5);

      //debug checks for improper input indices
      if ((centreIndex % COL_OFFSET != 0)  //Don't want the 4 regions from other half of detector
      && (centreIndex >= COL_OFFSET)  //Don't want the shared column to left of jet finding area
      && (centreIndex < (MAX_REGIONS_IN - COL_OFFSET))) { //Don't want column to the right either
                        
    if(detectJet(centreIndex, hfBoundary))
      {
        if (jetNum < MAX_JETS_OUT) {
            
          m_outputJets.at(jetNum).setRawsum(calcJetEnergy(centreIndex, hfBoundary));
          m_outputJets.at(jetNum).setDetId(calcJetPosition(centreIndex));
          m_outputJets.at(jetNum).setBx(m_inputRegions.at(centreIndex).bx());
          if(row < COL_OFFSET-4)  //if we are not in the HF, perform tauVeto analysis
        {
          m_outputJets.at(jetNum).setForward(false);
          m_outputJets.at(jetNum).setTauVeto(calcJetTauVeto(centreIndex,heBoundary));
        }
          else //can't be a tau jet because we are in the HF
        {
          m_outputJets.at(jetNum).setForward(true);
          m_outputJets.at(jetNum).setTauVeto(true);
        }
          ++jetNum;
        }
      }
    ++centreIndex;
      }
    }
  }
}

// Returns true if region index is the centre of a jet. Set boundary = true if at edge of HCAL.
bool L1GctTdrJetFinder::detectJet(const UShort centreIndex, const bool boundary) const
{
  if(!boundary)  //Not at boundary, so use 3*3 window of regions to determine if a jet
  {
    // Get the energy of the central region
    ULong testEt = m_inputRegions.at(centreIndex).et();
        
    //Test if our region qualifies as a jet by comparing its energy with the energies of the
    //surrounding eight regions.  In the event of neighbouring regions with identical energy,
    //this will locate the jet in the lower-most (furthest away from eta=0), left-most (least phi) region.
    if(testEt >  m_inputRegions.at(centreIndex-1-COL_OFFSET).et() &&
       testEt >  m_inputRegions.at(centreIndex - COL_OFFSET).et() &&
       testEt >  m_inputRegions.at(centreIndex+1-COL_OFFSET).et() &&
           
       testEt >= m_inputRegions.at(centreIndex - 1).et() &&
       testEt >  m_inputRegions.at(centreIndex + 1).et() &&
           
       testEt >= m_inputRegions.at(centreIndex-1+COL_OFFSET).et() &&
       testEt >= m_inputRegions.at(centreIndex + COL_OFFSET).et() &&
       testEt >= m_inputRegions.at(centreIndex+1+COL_OFFSET).et())
    {
      return true;
    }
//USE THIS BLOCK INSTEAD IF YOU WANT OVERFLOW BIT FUNCTIONALITY        
//*** BUT IT WILL NEED MODIFICATION SINCE L1GctRegion IS OBSOLETE ***
/*    // Get the energy of the central region & OR the overflow bit to become the MSB
    ULong testEt = (m_inputRegions.at(centreIndex).et() | (m_inputRegions.at(centreIndex).getOverFlow() << L1GctRegion::ET_BITWIDTH));
        
    //Test if our region qualifies as a jet by comparing its energy with the energies of the
    //surrounding eight regions.  In the event of neighbouring regions with identical energy,
    //this will locate the jet in the lower-most (furthest away from eta=0), left-most (least phi) region.
    if(testEt >  (m_inputRegions.at(centreIndex-1-COL_OFFSET).et() | (m_inputRegions.at(centreIndex-1-COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
       testEt >  (m_inputRegions.at(centreIndex - COL_OFFSET).et() | (m_inputRegions.at(centreIndex - COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
       testEt >  (m_inputRegions.at(centreIndex+1-COL_OFFSET).et() | (m_inputRegions.at(centreIndex+1-COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
           
       testEt >= (m_inputRegions.at(centreIndex - 1).et() | (m_inputRegions.at(centreIndex - 1).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
       testEt >  (m_inputRegions.at(centreIndex + 1).et() | (m_inputRegions.at(centreIndex + 1).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
           
       testEt >= (m_inputRegions.at(centreIndex-1+COL_OFFSET).et() | (m_inputRegions.at(centreIndex-1+COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
       testEt >= (m_inputRegions.at(centreIndex + COL_OFFSET).et() | (m_inputRegions.at(centreIndex + COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)) &&
       testEt >= (m_inputRegions.at(centreIndex+1+COL_OFFSET).et() | (m_inputRegions.at(centreIndex+1+COL_OFFSET).getOverFlow() << L1GctRegion::ET_BITWIDTH)))
    {
      return true;
    }
*/  //END OVERFLOW FUNCTIONALITY       
  }
  else    //...so only test surround 5 regions in our jet testing.
  {    
    // Get the energy of the central region
    // Don't need all the overflow bit adjustments as above, since we are in the HF here
    ULong testEt = m_inputRegions.at(centreIndex).et();        
        
    if(testEt >  m_inputRegions.at(centreIndex-1-COL_OFFSET).et() &&
       testEt >  m_inputRegions.at(centreIndex - COL_OFFSET).et() &&
       
       testEt >= m_inputRegions.at(centreIndex - 1).et() &&
           
       testEt >= m_inputRegions.at(centreIndex-1+COL_OFFSET).et() &&
       testEt >= m_inputRegions.at(centreIndex + COL_OFFSET).et())
    {
      return true;
    }
  }
  return false;           
}

//returns the energy sum of the nine regions centred (physically) about centreIndex
L1GctJetFinderBase::ULong L1GctTdrJetFinder::calcJetEnergy(const UShort centreIndex, const bool boundary) const
{
  ULong energy = 0;
    
  if(!boundary)
  {
    for(int column = -1; column <= +1; ++column)
    {
      energy += m_inputRegions.at(centreIndex-1 + (column*COL_OFFSET)).et() +
                m_inputRegions.at( centreIndex  + (column*COL_OFFSET)).et() +
                m_inputRegions.at(centreIndex+1 + (column*COL_OFFSET)).et();
    }
  }
  else
  {
    for(int column = -1; column <= +1; ++column)
    {
      energy += m_inputRegions.at(centreIndex-1 + (column*COL_OFFSET)).et() +
                m_inputRegions.at( centreIndex  + (column*COL_OFFSET)).et();
    }
  }

  return energy;                                   
}

// returns the encoded (eta, phi) position of the centre region
L1CaloRegionDetId L1GctTdrJetFinder::calcJetPosition(const UShort centreIndex) const
{
  return m_inputRegions.at(centreIndex).id();
}

// returns the combined tauveto of the nine regions centred (physically) about centreIndex. Set boundary = true if at edge of Endcap.
bool L1GctTdrJetFinder::calcJetTauVeto(const UShort centreIndex, const bool boundary) const
{
  bool partial[3] = {false, false, false};
    
  if(!boundary)
  {
    for(int column = -1; column <= +1; ++column)
    {
      partial[column+1] = m_inputRegions.at(centreIndex-1 + (column*COL_OFFSET)).tauVeto() ||
                          m_inputRegions.at( centreIndex  + (column*COL_OFFSET)).tauVeto() ||
                          m_inputRegions.at(centreIndex+1 + (column*COL_OFFSET)).tauVeto();
    }
  }
  else
  {
    for(int column = -1; column <= +1; ++column)
    {
      partial[column+1] = m_inputRegions.at(centreIndex-1 + (column*COL_OFFSET)).tauVeto() ||
                          m_inputRegions.at( centreIndex  + (column*COL_OFFSET)).tauVeto();
    }
  }
  return partial[0] || partial[1] || partial[2];
}