L1GctMet.cc

Go to the documentation of this file.

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

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

#include <math.h>

L1GctMet::L1GctMet(const unsigned ex, const unsigned ey, const L1GctMet::metAlgoType algo) :
  m_exComponent(ex),
  m_eyComponent(ey),
  m_algoType(algo),
  m_bitShift(0),
  m_htMissLut(new L1GctHtMissLut())
{ }

L1GctMet::L1GctMet(const etComponentType& ex, const etComponentType& ey, const metAlgoType algo) :
  m_exComponent(ex),
  m_eyComponent(ey),
  m_algoType(algo),
  m_bitShift(0),
  m_htMissLut(new L1GctHtMissLut())
{ }

L1GctMet::~L1GctMet() {}

// the Etmiss algorithm - external entry point
L1GctMet::etmiss_vec
L1GctMet::metVector () const
{
  etmiss_vec      result;
  etmiss_internal algoResult;
  switch (m_algoType)
    {
    case cordicTranslate:
      algoResult = cordicTranslateAlgo (m_exComponent.value(), m_eyComponent.value(), (m_exComponent.overFlow() || m_eyComponent.overFlow()) );
      break;

    case useHtMissLut:
      algoResult = useHtMissLutAlgo    (m_exComponent.value(), m_eyComponent.value(), (m_exComponent.overFlow() || m_eyComponent.overFlow()) );
      break;

    case oldGct:
      algoResult = oldGctAlgo          (m_exComponent.value(), m_eyComponent.value());
      break;

    case floatingPoint:
      algoResult = floatingPointAlgo   (m_exComponent.value(), m_eyComponent.value());
      break;

    default:
      algoResult.mag = 0;
      algoResult.phi = 0;
      break;
    }

  // The parameter m_bitShift allows us to discard additional LSB
  // in order to change the output scale. 
  result.mag.setValue(algoResult.mag>>(m_bitShift));
  result.phi.setValue(algoResult.phi);

  result.mag.setOverFlow( result.mag.overFlow() || inputOverFlow() );

  return result;

}

void L1GctMet::setExComponent(const unsigned ex) {
  etComponentType temp(ex);
  setExComponent(temp);
}

void L1GctMet::setEyComponent(const unsigned ey) {
  etComponentType temp(ey);
  setEyComponent(temp);
}

// private member functions - the different algorithms:

L1GctMet::etmiss_internal
L1GctMet::cordicTranslateAlgo (const int ex, const int ey, const bool of) const
{
  //---------------------------------------------------------------------------------
  //
  // This is an implementation of the CORDIC algorithm (COordinate Rotation for DIgital Computers)
  //
  // Starting from an initial two-component vector ex, ey, we perform successively smaller rotations
  // to transform the y component to zero. At the end of the procedure, the x component is the magnitude
  // of the original vector, scaled by a known constant factor. The azimuth angle phi is the sum of the
  // rotations applied.
  //
  // The algorithm can be used in a number of different variants for calculation of trigonometric
  // and hyperbolic functions as well as exponentials, logarithms and square roots. This variant
  // is called the "vector translation" mode in the Xilinx documentation.
  //
  // Original references:
  // Volder, J., "The CORDIC Trigonometric Computing Technique" IRE Trans. Electronic Computing, Vol.
  // EC-8, Sept. 1959, pp330-334
  // Walther, J.S., "A Unified Algorithm for Elementary Functions," Spring Joint computer conf., 1971,
  // proc., pp379-385
  //
  // Other information sources: http://www.xilinx.com/support/documentation/ip_documentation/cordic.pdf;
  // http://www.fpga-guru.com/files/crdcsrvy.pdf; and http://en.wikipedia.org/wiki/CORDIC
  //
  //---------------------------------------------------------------------------------

  etmiss_internal result;

  static const int of_val = 0x1FFF; // set components to 8191 (decimal) if there's an overflow on the input

  static const int n_iterations = 6;
  // The angle units here are 1/32 of a 5 degree bin.
  // So a 90 degree rotation is 32*18=576 or 240 hex.
  const int cordic_angles[n_iterations] = { 0x120, 0x0AA, 0x05A, 0x02E, 0x017, 0x00B };
  const int cordic_starting_angle_090 = 0x240;
  const int cordic_starting_angle_270 = 0x6C0;
  const int cordic_angle_360 = 0x900;

  const int cordic_scale_factor = 0x26E; // decimal 622

  int  x, y;
  int dx,dy;
  int  z;

  if (of) {
    x =  of_val;
    y = -of_val;
    z = cordic_starting_angle_090;
  } else {
    if (ey>=0) {
      x =  ey;
      y = -ex;
      z = cordic_starting_angle_090;
    } else {
      x = -ey;
      y =  ex;
      z = cordic_starting_angle_270;
    }
  }

  for (int i=0; i<n_iterations; i++) {
    dx = cordicShiftAndRoundBits (y,i);
    dy = cordicShiftAndRoundBits (x,i);
    if (y>=0) {
      x = x + dx;
      y = y - dy;
      z = z + cordic_angles[i];
    } else {
      x = x - dx;
      y = y + dy;
      z = z - cordic_angles[i];
    }
  }

  int scaled_magnitude = x * cordic_scale_factor;
  int adjusted_angle   = ( (z < 0) ? (z + cordic_angle_360) : z ) % cordic_angle_360;
  result.mag = scaled_magnitude >> 10;
  result.phi = adjusted_angle   >> 5;
  if (result.mag > (unsigned) of_val) result.mag = (unsigned) of_val;
  return result;
}

int L1GctMet::cordicShiftAndRoundBits (const int e, const unsigned nBits) const
{
  int r;
  if (nBits==0) {
    r = e;
  } else {
    r = (((e>>(nBits-1)) + 1)>>1);
  }
  return r;
}


L1GctMet::etmiss_internal
L1GctMet::useHtMissLutAlgo (const int ex, const int ey, const bool of) const
{
  // The firmware discards the LSB of the input values, before forming
  // the address for the LUT. We do the same here.
  static const int maxComponent    = 1<<L1GctHtMissLut::kHxOrHyMissComponentNBits;
  static const int componentMask   = maxComponent-1;
  static const int maxPosComponent = componentMask >> 1;

  static const int maxInput = 1<<(L1GctHtMissLut::kHxOrHyMissComponentNBits+kExOrEyMissComponentShift-1);

  static const unsigned resultMagMask = (1<<L1GctHtMissLut::kHtMissMagnitudeNBits) - 1;
  static const unsigned resultPhiMask = (1<<L1GctHtMissLut::kHtMissAngleNBits)     - 1;

  etmiss_internal result;

  if (m_htMissLut == 0) {

    result.mag = 0;
    result.phi = 0;

  } else {

    // Extract the bit fields of the input components to be used for the LUT address
    int hxCompBits = (ex >> kExOrEyMissComponentShift) & componentMask;
    int hyCompBits = (ey >> kExOrEyMissComponentShift) & componentMask;

    if (of || (abs(ex) >= maxInput) || (abs(ey) >= maxInput)) {
      hxCompBits = maxPosComponent;
      hyCompBits = maxPosComponent;
    }

    // Perform the table lookup to get the missing Ht magnitude and phi
    uint16_t lutAddress = static_cast<uint16_t> ( (hxCompBits << L1GctHtMissLut::kHxOrHyMissComponentNBits) | hyCompBits );

    uint16_t lutData = m_htMissLut->lutValue(lutAddress);

    result.mag = static_cast<unsigned>(lutData>>L1GctHtMissLut::kHtMissAngleNBits) & resultMagMask;
    result.phi = static_cast<unsigned>(lutData) & resultPhiMask;

  }

  return result;
}

L1GctMet::etmiss_internal
L1GctMet::oldGctAlgo (const int ex, const int ey) const
{
  //---------------------------------------------------------------------------------
  //
  // Calculates magnitude and direction of missing Et, given measured Ex and Ey.
  //
  // The algorithm used is suitable for implementation in hardware, using integer
  // multiplication, addition and comparison and bit shifting operations.
  //
  // Proceed in two stages. The first stage gives a result that lies between
  // 92% and 100% of the true Et, with the direction measured in 45 degree bins.
  // The final precision depends on the number of factors used in corrFact.
  // The present version with eleven factors gives a precision of 1% on Et, and
  // finds the direction to the nearest 5 degrees.
  //
  //---------------------------------------------------------------------------------
  etmiss_internal result;

  unsigned eneCoarse, phiCoarse;
  unsigned eneCorect, phiCorect;

  const unsigned root2fact = 181;
  const unsigned corrFact[11] = {24, 39, 51, 60, 69, 77, 83, 89, 95, 101, 106};
  const unsigned corrDphi[11] = { 0,  1,  1,  2,  2,  3,  3,  3,  3,   4,   4};

  std::vector<bool> s(3);
  unsigned Mx, My, Mw;

  unsigned Dx, Dy;
  unsigned eFact;

  unsigned b,phibin;
  bool midphi=false;

  // Here's the coarse calculation, with just one multiply operation
  //
  My = static_cast<unsigned>(abs(ey));
  Mx = static_cast<unsigned>(abs(ex));
  Mw = (((Mx+My)*root2fact)+0x80)>>8;

  s.at(0) = (ey<0);
  s.at(1) = (ex<0);
  s.at(2) = (My>Mx);

  phibin = 0; b = 0;
  for (int i=0; i<3; i++) {
    if (s.at(i)) { b=1-b;} phibin = 2*phibin + b;
  }

  eneCoarse = std::max(std::max(Mx, My), Mw);
  phiCoarse = phibin*9;

  // For the fine calculation we multiply both input components
  // by all the factors in the corrFact list in order to find
  // the required corrections to the energy and angle
  //
  for (eFact=0; eFact<10; eFact++) {
    Dx = (Mx*corrFact[eFact])>>8;
    Dy = (My*corrFact[eFact])>>8;
    if         ((Dx>=My) || (Dy>=Mx))         {midphi=false; break;}
    if ((Mx+Dx)>=(My-Dy) && (My+Dy)>=(Mx-Dx)) {midphi=true;  break;}
  }
  eneCorect = (eneCoarse*(128+eFact))>>7;
  if (midphi ^ (b==1)) {
    phiCorect = phiCoarse + 8 - corrDphi[eFact];
  } else {
    phiCorect = phiCoarse + corrDphi[eFact];
  }

  // Store the result of the calculation
  //
  result.mag = eneCorect;
  result.phi = phiCorect;

  return result;
}

L1GctMet::etmiss_internal
L1GctMet::floatingPointAlgo (const int ex, const int ey) const
{

  etmiss_internal result;

  double fx = static_cast<double>(ex);
  double fy = static_cast<double>(ey);
  double fmag = sqrt(fx*fx + fy*fy);
  double fphi = 36.*atan2(fy, fx)/M_PI;

  result.mag = static_cast<unsigned>(fmag);
  if (fphi>=0) {
    result.phi = static_cast<unsigned>(fphi);
  } else {
    result.phi = static_cast<unsigned>(fphi+72.);
  }

  return result;

}

void L1GctMet::setEtScale(const L1CaloEtScale* const fn)
{
  m_htMissLut->setEtScale(fn);
}

void L1GctMet::setEtComponentLsb(const double lsb)
{
  m_htMissLut->setExEyLsb( lsb * static_cast<double>( 1<<kExOrEyMissComponentShift ) );
}

const L1CaloEtScale* L1GctMet::etScale() const
{ 
  return m_htMissLut->etScale();
}

const double L1GctMet::componentLsb() const
{ 
  return m_htMissLut->componentLsb();
}

const bool L1GctMet::inputOverFlow() const
{
  bool result = m_exComponent.overFlow() || m_eyComponent.overFlow();

  if (m_algoType == useHtMissLut) {
    static const int maxComponentInput = ( 1 << (L1GctHtMissLut::kHxOrHyMissComponentNBits+kExOrEyMissComponentShift-1) ) - 1;

    // Emulate the (symmetric) overflow condition used in the firmware
    result |= (m_exComponent.value() > maxComponentInput) || (m_exComponent.value() < -maxComponentInput) ||
      (m_eyComponent.value() > maxComponentInput) || (m_eyComponent.value() < -maxComponentInput);
  }

  return result;
}