CPMTobAlgorithm.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
/***************************************************************************
                          CPMTobAlgorithm.cxx  -  description
                             -------------------
    begin                : Thu 6 Mar 2014
                           Derived from CPAlgorithm (Run 1 version)
    email                : Alan.Watson@cern.ch
 ***************************************************************************/
 
#include "TrigT1CaloUtils/CPMTobAlgorithm.h"
#include "TrigConfL1Data/L1DataDef.h"
#include "TrigT1Interfaces/TrigT1CaloDefs.h"
#include "TrigT1Interfaces/CoordinateRange.h"
#include "TrigT1CaloUtils/CoordToHardware.h"
#include "TrigConfData/L1Menu.h"
 
#include <math.h>
 
namespace LVL1 {
using namespace TrigConf;
 
const unsigned int CPMTobAlgorithm::m_maxClus = 0xFF;
const double CPMTobAlgorithm::m_maxEta = 2.5;
 
const unsigned int CPMTobAlgorithm::m_emLUT_ClusterFirstBit = 1;
const unsigned int CPMTobAlgorithm::m_emLUT_ClusterNBits = 6;
const unsigned int CPMTobAlgorithm::m_emLUT_EMIsolFirstBit = 0;
const unsigned int CPMTobAlgorithm::m_emLUT_EMIsolNBits = 4;
const unsigned int CPMTobAlgorithm::m_emLUT_HadVetoFirstBit = 0;
const unsigned int CPMTobAlgorithm::m_emLUT_HadVetoNBits = 3;
  
const unsigned int CPMTobAlgorithm::m_tauLUT_ClusterFirstBit = 1;
const unsigned int CPMTobAlgorithm::m_tauLUT_ClusterNBits = 7;
const unsigned int CPMTobAlgorithm::m_tauLUT_EMIsolFirstBit = 0;
const unsigned int CPMTobAlgorithm::m_tauLUT_EMIsolNBits = 6;
 
const unsigned int CPMTobAlgorithm::m_noIsol = 999;
 
LVL1::CPMTobAlgorithm::CPMTobAlgorithm( double eta, double phi, const xAOD::CPMTowerMap_t* ttContainer,
                                        const TrigConf::L1Menu * l1menu, int slice ):
  m_l1menu(l1menu),
  m_Core(0),
  m_EMClus(0),
  m_TauClus(0),
  m_EMIsolWord(0),
  m_TauIsolWord(0),
  m_EMCore(0),
  m_EMIsol(0),
  m_HadCore(0),
  m_HadIsol(0),
  m_EtMax(false),
  m_EMThresh(false),
  m_TauThresh(false)
  //m_debug(false)
{
  TriggerTowerKey get(phi-M_PI/128.,eta-0.025);
  Coordinate refCoord = get.coord();
  m_refEta = refCoord.eta();
  m_refPhi = refCoord.phi();
 
  // Set coordinate of centre of RoI, starting from reference tower coordinate
  setRoICoord(m_refEta, m_refPhi);
  
  std::vector<unsigned int> emEt(4);
  std::vector<unsigned int> cores(9);
  for (int etaOffset = -1; etaOffset <= 2; etaOffset++) {
    double tempEta = m_refEta + etaOffset*0.1;
    for (int phiOffset = -1; phiOffset <= 2; phiOffset++) {
      double tempPhi = m_refPhi + phiOffset*M_PI/32;
      int key = get.ttKey(tempPhi, tempEta);
      xAOD::CPMTowerMap_t::const_iterator tt = ttContainer->find(key);
      if (tt != ttContainer->end() && fabs(tempEta) < m_maxEta) {
        // Get the TT ET values once here, rather than repeat function calls
        int emTT = 0;
        int hadTT = 0;
        if (slice < 0) {  // Default to using peak slice
          emTT  = (tt->second)->emEnergy();
          hadTT = (tt->second)->hadEnergy();
        }
        else {            // Use user-specified slice
          emTT  = (tt->second)->emSliceEnergy(slice);
          hadTT = (tt->second)->hadSliceEnergy(slice);
        }
        // EM Trigger Clusters and Core clusters
        if (etaOffset >= 0 && etaOffset <= 1 && phiOffset >= 0 && phiOffset <= 1) {
          emEt[phiOffset+2*etaOffset] = emTT;
          m_EMCore += emTT;
          m_HadCore += hadTT;
        }
        // Isolation sums
        else {
          m_EMIsol  += emTT;
          m_HadIsol += hadTT;
        }
        // Each tower is part of up to 4 RoI core clusters
        if (etaOffset >= 0) {
          if (phiOffset >= 0) cores[phiOffset+3*etaOffset] += emTT + hadTT;
          if (phiOffset < 2)  cores[phiOffset+3*etaOffset+1] += emTT + hadTT;
        }
        if (etaOffset < 2) {
          if (phiOffset >= 0) cores[phiOffset+3*etaOffset+3] += emTT + hadTT;
          if (phiOffset < 2)  cores[phiOffset+3*etaOffset+4] += emTT + hadTT;
        }
 
      }
    } // phi offset loop
  }   // eta offset loop
 
  // Core of this window is cores[4]
  m_Core = cores[4];
 
  // Form most energetic 2-tower EM cluster
  m_EMClus = ( (emEt[0] > emEt[3]) ? emEt[0] : emEt[3] )
           + ( (emEt[1] > emEt[2]) ? emEt[1] : emEt[2] );
 
  // Tau cluster
  m_TauClus = m_EMClus + m_HadCore;
  
  // Check whether RoI condition is met.
  testEtMax(cores);
 
  // test cluster and isolation conditions
  emAlgorithm();
  tauAlgorithm();
  
}
 
LVL1::CPMTobAlgorithm::~CPMTobAlgorithm(){
}
 
void LVL1::CPMTobAlgorithm::setRoICoord(double eta, double phi) {
 
  // TriggerTowerKey::coord() returns centre of TriggerTower
  TriggerTowerKey keyLL(phi,eta);
  Coordinate lowerLeft = keyLL.coord();
 
  // Get centre of tower at eta+1, phi+1, i.e. opposite corner of RoI core
  double dEta = keyLL.dEta(lowerLeft);
  double dPhi = keyLL.dPhi(lowerLeft);  
 
  TriggerTowerKey keyUR(phi+dPhi, eta + dEta); 
  Coordinate upperRight = keyUR.coord();
 
  // CoordinateRange object will compute centre, correcting for wrap-around
  CoordinateRange roi(lowerLeft.phi(),upperRight.phi(),lowerLeft.eta(),upperRight.eta());
  m_eta = roi.eta();
  m_phi = roi.phi();
 
  // Range check (shouldn't be needed, but floating point precision errors possible)
  if       (m_eta > 2.5) m_eta = 2.5;
  else if (m_eta < -2.5) m_eta = -2.5;
}
 
void LVL1::CPMTobAlgorithm::testEtMax(const std::vector<unsigned int>& cores) {
 
  // RoI condition test
  m_EtMax = true;
  for (int i = 0; i < 4; i++) if (m_Core < cores[i])  m_EtMax = false;
  for (int i = 5; i < 9; i++) if (m_Core <= cores[i]) m_EtMax = false;
  
}
 
 
void LVL1::CPMTobAlgorithm::emAlgorithm() {
  
  // Initialise to correct defaults, to be sure
  m_EMThresh = false;
  m_EMIsolWord = 0;
  
  // Don't waste time if it isn't an RoI candidate
  if (!m_EtMax) return;
  
  // Does this pass min TOB pT cut?
  // Need to worry about digit scale not being GeV
  float emEnergyScale{1}; 
  unsigned int tobPtMinEM{0};
  emEnergyScale = m_l1menu->thrExtraInfo().EM().emScale();
  tobPtMinEM = m_l1menu->thrExtraInfo().EM().ptMinToTopo();
 
  unsigned int thresh = tobPtMinEM * emEnergyScale;
 
  if (m_EMClus <= thresh) return;
  
  // Passed, so there is a TOB here
  m_EMThresh = true;
 
  // Now calculate isolation word
  // Define cluster and isolation values with appropriate scales and range for isolation LUT  
  unsigned int clusMask    = ( (1<<m_emLUT_ClusterNBits)-1 )<<m_emLUT_ClusterFirstBit;
  unsigned int emIsolMask  = ( (1<<m_emLUT_EMIsolNBits) -1 )<<m_emLUT_EMIsolFirstBit;
  unsigned int hadIsolMask = ( (1<<m_emLUT_HadVetoNBits)-1 )<<m_emLUT_HadVetoFirstBit;
 
  // If cluster overflows, set all threshold bits
  if (m_EMClus >= clusMask) {
    m_EMIsolWord = 0x1F;
    return;
  }
  
  // For consistency with LUT filling and menu parameters, convert ET sums to 100 MeV units
  int clus    = (( (m_EMClus < clusMask)     ? (m_EMClus  & clusMask)    : clusMask ) * 10)/emEnergyScale;
  int emIsol  = (( (m_EMIsol < emIsolMask)   ? (m_EMIsol  & emIsolMask)  : emIsolMask ) * 10)/emEnergyScale;
  int hadVeto = (( (m_HadCore < hadIsolMask) ? (m_HadCore & hadIsolMask) : hadIsolMask ) * 10)/emEnergyScale;
 
  std::vector<int> emisolcuts(5,m_noIsol);
  std::vector<int> hadisolcuts(5,m_noIsol);
  for(size_t bit = 1; bit<=5; ++bit) {
    const TrigConf::IsolationLegacy & iso = m_l1menu->thrExtraInfo().EM().isolation("EMIsoForEMthr", bit);
    if(! iso.isDefined() ) {
      continue;
    }
    int offset = iso.offset();
    int slope  = iso.slope();
    int mincut = iso.mincut();
    int upperlimit = iso.upperlimit()*10;
    if (clus < upperlimit && bit >= 1 && bit <= 5) {
      // As "slope" = 10* slope, rescale clus too to get fraction right.
      int cut = offset + (slope != 0 ? 10*clus/slope : 0);
      if (cut < mincut) cut = mincut;
      emisolcuts[bit-1] = cut;
    }
  }
  for(size_t bit = 1; bit<=5; ++bit) {
    const TrigConf::IsolationLegacy & iso = m_l1menu->thrExtraInfo().EM().isolation("HAIsoForEMthr", bit);
    if(! iso.isDefined() ) {
      continue;
    }
    int offset = iso.offset();
    int slope  = iso.slope();
    int mincut = iso.mincut();
    int upperlimit = iso.upperlimit()*10;
    if (clus < upperlimit && bit >= 1 && bit <= 5) {
      // As "slope" = 10* slope, rescale clus too to get fraction right.
      int cut = offset + (slope != 0 ? 10*clus/slope : 0);
      if (cut < mincut) cut = mincut;
      hadisolcuts[bit-1] = cut;
    }
  }
 
  // Set isolation bits
  for (unsigned int isol = 0; isol < 5; ++isol) {
    if (emIsol <= emisolcuts[isol] && hadVeto <= hadisolcuts[isol]) m_EMIsolWord += (1 << isol);
  }
    
}
 
 
void LVL1::CPMTobAlgorithm::tauAlgorithm() {
  
  // Initialise to correct defaults, to be sure
  m_TauThresh = false;
  m_TauIsolWord = 0;
  
  // Don't waste time if it isn't an RoI candidate
  if (!m_EtMax) return;
  
  // Does this pass min TOB pT cut?
  // Need to worry about digit scale not being GeV
  float emEnergyScale{1}; 
  unsigned int tobPtMinTau{0};
  emEnergyScale = m_l1menu->thrExtraInfo().EM().emScale();
  tobPtMinTau = m_l1menu->thrExtraInfo().TAU().ptMinToTopo();
  unsigned int thresh = tobPtMinTau * emEnergyScale;
 
  if (m_TauClus <= thresh) return;
  
  // Passed, so there is a TOB here
  m_TauThresh = true;
  
  // Now calculate isolation word
  // Define cluster and isolation values with appropriate scales and range for isolation LUT
  unsigned int clusMask    = ( (1<<m_tauLUT_ClusterNBits)-1 )<<m_tauLUT_ClusterFirstBit;
  unsigned int emIsolMask  = ( (1<<m_tauLUT_EMIsolNBits) -1 )<<m_tauLUT_EMIsolFirstBit;
 
  // If cluster overflows, set all threshold bits
  if (m_TauClus >= clusMask) {
     m_TauIsolWord = 0x1F;
     return;
  }
                 
  // Convert to 100 MeV units to match parameters
  int clus    = (( (m_TauClus < clusMask)    ? (m_TauClus & clusMask)    : clusMask ) * 10)/emEnergyScale;
  int emIsol  = (( (m_EMIsol < emIsolMask)   ? (m_EMIsol  & emIsolMask)  : emIsolMask ) * 10)/emEnergyScale;
  
  // Get isolation values from menu (placeholder code - example logic)
  std::vector<int> emisolcuts(5,m_noIsol);
  for(size_t bit = 1; bit<=5; ++bit) {
    const TrigConf::IsolationLegacy & iso = m_l1menu->thrExtraInfo().TAU().isolation("EMIsoForTAUthr", bit);
    if(! iso.isDefined() ) {
      continue;
    }
    int offset = iso.offset();
    int slope  = iso.slope();
    int mincut = iso.mincut();
    int upperlimit = iso.upperlimit()*10;
    if (clus < upperlimit && bit >= 1 && bit <= 5) {
      // As "slope" = 10* slope, rescale clus too to get fraction right.
      int cut = offset + (slope != 0 ? 10*clus/slope : 0);
      if (cut < mincut) cut = mincut;
      emisolcuts[bit-1] = cut;
    }
  }
 
  // Set isolation bits
  for (unsigned int isol = 0; isol < 5; ++isol) 
    if (emIsol <= emisolcuts[isol]) m_TauIsolWord += (1 << isol);
 
}
 
bool LVL1::CPMTobAlgorithm::isEMRoI(){
  
  bool passed = ( m_EMThresh && m_EtMax );
 
  return passed;
}
 
bool LVL1::CPMTobAlgorithm::isTauRoI(){
 
  bool passed = ( m_TauThresh && m_EtMax );
 
  return passed;
}
 
// Public accessor methods follow
 
int LVL1::CPMTobAlgorithm::CoreET() {
  return m_Core;
}
 
int LVL1::CPMTobAlgorithm::EMClusET() {
  return ( (m_EMClus < m_maxClus) ? m_EMClus : m_maxClus );
}
 
int LVL1::CPMTobAlgorithm::TauClusET() {
  return ( (m_TauClus < m_maxClus) ? m_TauClus : m_maxClus );
}
 
int LVL1::CPMTobAlgorithm::EMIsolWord() {
  return m_EMIsolWord;
}
 
int LVL1::CPMTobAlgorithm::TauIsolWord() {
  return m_TauIsolWord;
}
 
int LVL1::CPMTobAlgorithm::EMIsolET() {
  return m_EMIsol;
}
 
int LVL1::CPMTobAlgorithm::HadIsolET() {
  return m_HadIsol;
}
 
int LVL1::CPMTobAlgorithm::EMCoreET() {
  return m_EMCore;
}
 
int LVL1::CPMTobAlgorithm::HadCoreET() {
  return m_HadCore;
}
 
bool LVL1::CPMTobAlgorithm::isEtMax() {
  return m_EtMax;
}
 
int LVL1::CPMTobAlgorithm::EMLUTClus() {
  int ET = (m_EMClus >> m_emLUT_ClusterFirstBit);
  if (ET > (1<< m_emLUT_ClusterNBits) - 1) ET = (1<<m_emLUT_ClusterNBits) - 1;
  return (ET << m_emLUT_ClusterFirstBit);
}
 
int LVL1::CPMTobAlgorithm::TauLUTClus() {
  int ET = (m_TauClus >> m_tauLUT_ClusterFirstBit);
  if (ET > (1<<m_tauLUT_ClusterNBits) - 1) ET = (1<<m_tauLUT_ClusterNBits) - 1;
  return (ET << m_tauLUT_ClusterFirstBit);
}
 
int LVL1::CPMTobAlgorithm::EMLUTEMIsol() {
  int ET = (m_EMIsol >> m_emLUT_EMIsolFirstBit);
  if (ET > (1<<m_emLUT_EMIsolNBits) - 1) ET = (1<<m_emLUT_EMIsolNBits) - 1;
  return (ET << m_emLUT_EMIsolFirstBit);
}
 
int LVL1::CPMTobAlgorithm::EMLUTHadVeto() {
  int ET = (m_HadCore >> m_emLUT_HadVetoFirstBit);
  if (ET > (1<<m_emLUT_HadVetoNBits) - 1) ET = (1<<m_emLUT_HadVetoNBits) - 1;
  return (ET << m_emLUT_HadVetoFirstBit);
}
 
int LVL1::CPMTobAlgorithm::TauLUTEMIsol() {
  int ET = (m_EMIsol >> m_tauLUT_EMIsolFirstBit);
  if (ET > (1<<m_tauLUT_EMIsolNBits) - 1) ET = (1<<m_tauLUT_EMIsolNBits) - 1;
  return (ET << m_tauLUT_EMIsolFirstBit);
}
 
double LVL1::CPMTobAlgorithm::eta() {
  return m_eta;
}
 
double LVL1::CPMTobAlgorithm::phi() {
  return ( (m_phi <= M_PI) ? m_phi : m_phi - 2.*M_PI);
}
 
xAOD::CPMTobRoI* LVL1::CPMTobAlgorithm::EMCPMTobRoI() {
  if (isEMRoI()) {
    xAOD::CPMTobRoI* pointer = createTobRoI(TrigT1CaloDefs::emTobType); 
    return pointer;
  }
  
  return 0;
 
}
 
xAOD::CPMTobRoI* LVL1::CPMTobAlgorithm::TauCPMTobRoI() {
  
  xAOD::CPMTobRoI* pointer = 0;
  
  if (isTauRoI()) pointer = createTobRoI(TrigT1CaloDefs::tauTobType);
  
  return pointer;
 
}
 
xAOD::CPMTobRoI* LVL1::CPMTobAlgorithm::createTobRoI(int type) {
  
  xAOD::CPMTobRoI* pointer = 0;
  
  int et = 0;
  int isol = 0;
  
  if (type == TrigT1CaloDefs::emTobType) {
    et   = EMClusET();
    isol = EMIsolWord();
  }
  else if (type == TrigT1CaloDefs::tauTobType) {
    et   = TauClusET();
    isol = TauIsolWord();
  }
  else return pointer;
  
  // Need to calculate hardware coordinate
  // Object coordinate is RoI centre, but this works better if use reference tower coord
  float eta = m_eta - 0.05;
  float phi = m_phi - M_PI/64;
  if (phi < 0) phi += 2*M_PI; 
  Coordinate coord(phi, eta);
  CoordToHardware convertor;
  
  int crate = convertor.cpCrate(coord);
  int module = convertor.cpModule(coord);
  int chip = convertor.cpModuleFPGA(coord);
  int location = convertor.cpModuleLocalRoI(coord);
  
  /*
  float phi = ( (m_phi >= 0) ?  m_phi : 2.*M_PI + m_phi );
  int crate = 2*phi/M_PI;
  int module = (m_eta + 2.8)/0.4 + 1;
  int iPhi = (phi - crate*M_PI/2)*32/M_PI;
  int chip = iPhi/2;
  int iEta = (m_eta + 2.8 - module*0.4)/0.1;
  int iSide = iEta/2;
  int location = (iSide << 2) + ((iPhi&1) << 1) + (iEta&1);
  */
  
  pointer = new xAOD::CPMTobRoI;
  pointer->makePrivateStore();
  pointer->initialize(crate, module, chip, location, type, et, isol);
  
  return pointer;
  
}
 
 
} // end of namespace bracket