GlobalCaloTrigger/src/L1GctMet.cc

0001 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctMet.h"
0002
0003 #include "L1Trigger/GlobalCaloTrigger/interface/L1GctHtMissLut.h"
0004
0005 #include <cmath>
0006
0007 L1GctMet::L1GctMet(const unsigned ex, const unsigned ey, const L1GctMet::metAlgoType algo)
0008     : m_exComponent(ex), m_eyComponent(ey), m_algoType(algo), m_bitShift(0), m_htMissLut(new L1GctHtMissLut()) {}
0009
0010 L1GctMet::L1GctMet(const etComponentType& ex, const etComponentType& ey, const metAlgoType algo)
0011     : m_exComponent(ex), m_eyComponent(ey), m_algoType(algo), m_bitShift(0), m_htMissLut(new L1GctHtMissLut()) {}
0012
0013 L1GctMet::~L1GctMet() {}
0014
0015 // the Etmiss algorithm - external entry point
0016 L1GctMet::etmiss_vec L1GctMet::metVector() const {
0017   etmiss_vec result;
0018   etmiss_internal algoResult;
0019   switch (m_algoType) {
0020     case cordicTranslate:
0021       algoResult = cordicTranslateAlgo(
0022           m_exComponent.value(), m_eyComponent.value(), (m_exComponent.overFlow() || m_eyComponent.overFlow()));
0023       break;
0024
0025     case useHtMissLut:
0026       algoResult = useHtMissLutAlgo(
0027           m_exComponent.value(), m_eyComponent.value(), (m_exComponent.overFlow() || m_eyComponent.overFlow()));
0028       break;
0029
0030     case oldGct:
0031       algoResult = oldGctAlgo(m_exComponent.value(), m_eyComponent.value());
0032       break;
0033
0034     case floatingPoint:
0035       algoResult = floatingPointAlgo(m_exComponent.value(), m_eyComponent.value());
0036       break;
0037
0038     default:
0039       algoResult.mag = 0;
0040       algoResult.phi = 0;
0041       break;
0042   }
0043
0044   // The parameter m_bitShift allows us to discard additional LSB
0045   // in order to change the output scale.
0046   result.mag.setValue(algoResult.mag >> (m_bitShift));
0047   result.phi.setValue(algoResult.phi);
0048
0049   result.mag.setOverFlow(result.mag.overFlow() || inputOverFlow());
0050
0051   return result;
0052 }
0053
0054 void L1GctMet::setExComponent(const unsigned ex) {
0055   etComponentType temp(ex);
0056   setExComponent(temp);
0057 }
0058
0059 void L1GctMet::setEyComponent(const unsigned ey) {
0060   etComponentType temp(ey);
0061   setEyComponent(temp);
0062 }
0063
0064 // private member functions - the different algorithms:
0065
0066 L1GctMet::etmiss_internal L1GctMet::cordicTranslateAlgo(const int ex, const int ey, const bool of) const {
0067   //---------------------------------------------------------------------------------
0068   //
0069   // This is an implementation of the CORDIC algorithm (COordinate Rotation for DIgital Computers)
0070   //
0071   // Starting from an initial two-component vector ex, ey, we perform successively smaller rotations
0072   // to transform the y component to zero. At the end of the procedure, the x component is the magnitude
0073   // of the original vector, scaled by a known constant factor. The azimuth angle phi is the sum of the
0074   // rotations applied.
0075   //
0076   // The algorithm can be used in a number of different variants for calculation of trigonometric
0077   // and hyperbolic functions as well as exponentials, logarithms and square roots. This variant
0078   // is called the "vector translation" mode in the Xilinx documentation.
0079   //
0080   // Original references:
0081   // Volder, J., "The CORDIC Trigonometric Computing Technique" IRE Trans. Electronic Computing, Vol.
0082   // EC-8, Sept. 1959, pp330-334
0083   // Walther, J.S., "A Unified Algorithm for Elementary Functions," Spring Joint computer conf., 1971,
0084   // proc., pp379-385
0085   //
0086   // Other information sources: http://www.xilinx.com/support/documentation/ip_documentation/cordic.pdf;
0087   // http://www.fpga-guru.com/files/crdcsrvy.pdf; and http://en.wikipedia.org/wiki/CORDIC
0088   //
0089   //---------------------------------------------------------------------------------
0090
0091   etmiss_internal result;
0092
0093   static const int of_val = 0x1FFF;  // set components to 8191 (decimal) if there's an overflow on the input
0094
0095   static const int n_iterations = 6;
0096   // The angle units here are 1/32 of a 5 degree bin.
0097   // So a 90 degree rotation is 32*18=576 or 240 hex.
0098   const int cordic_angles[n_iterations] = {0x120, 0x0AA, 0x05A, 0x02E, 0x017, 0x00B};
0099   const int cordic_starting_angle_090 = 0x240;
0100   const int cordic_starting_angle_270 = 0x6C0;
0101   const int cordic_angle_360 = 0x900;
0102
0103   const int cordic_scale_factor = 0x26E;  // decimal 622
0104
0105   int x, y;
0106   int dx, dy;
0107   int z;
0108
0109   if (of) {
0110     x = of_val;
0111     y = -of_val;
0112     z = cordic_starting_angle_090;
0113   } else {
0114     if (ey >= 0) {
0115       x = ey;
0116       y = -ex;
0117       z = cordic_starting_angle_090;
0118     } else {
0119       x = -ey;
0120       y = ex;
0121       z = cordic_starting_angle_270;
0122     }
0123   }
0124
0125   for (int i = 0; i < n_iterations; i++) {
0126     dx = cordicShiftAndRoundBits(y, i);
0127     dy = cordicShiftAndRoundBits(x, i);
0128     if (y >= 0) {
0129       x = x + dx;
0130       y = y - dy;
0131       z = z + cordic_angles[i];
0132     } else {
0133       x = x - dx;
0134       y = y + dy;
0135       z = z - cordic_angles[i];
0136     }
0137   }
0138
0139   int scaled_magnitude = x * cordic_scale_factor;
0140   int adjusted_angle = ((z < 0) ? (z + cordic_angle_360) : z) % cordic_angle_360;
0141   result.mag = scaled_magnitude >> 10;
0142   result.phi = adjusted_angle >> 5;
0143   if (result.mag > (unsigned)of_val)
0144     result.mag = (unsigned)of_val;
0145   return result;
0146 }
0147
0148 int L1GctMet::cordicShiftAndRoundBits(const int e, const unsigned nBits) const {
0149   int r;
0150   if (nBits == 0) {
0151     r = e;
0152   } else {
0153     r = (((e >> (nBits - 1)) + 1) >> 1);
0154   }
0155   return r;
0156 }
0157
0158 L1GctMet::etmiss_internal L1GctMet::useHtMissLutAlgo(const int ex, const int ey, const bool of) const {
0159   // The firmware discards the LSB of the input values, before forming
0160   // the address for the LUT. We do the same here.
0161   static const int maxComponent = 1 << L1GctHtMissLut::kHxOrHyMissComponentNBits;
0162   static const int componentMask = maxComponent - 1;
0163   static const int maxPosComponent = componentMask >> 1;
0164
0165   static const int maxInput = 1 << (L1GctHtMissLut::kHxOrHyMissComponentNBits + kExOrEyMissComponentShift - 1);
0166
0167   static const unsigned resultMagMask = (1 << L1GctHtMissLut::kHtMissMagnitudeNBits) - 1;
0168   static const unsigned resultPhiMask = (1 << L1GctHtMissLut::kHtMissAngleNBits) - 1;
0169
0170   etmiss_internal result;
0171
0172   if (m_htMissLut == nullptr) {
0173     result.mag = 0;
0174     result.phi = 0;
0175
0176   } else {
0177     // Extract the bit fields of the input components to be used for the LUT address
0178     int hxCompBits = (ex >> kExOrEyMissComponentShift) & componentMask;
0179     int hyCompBits = (ey >> kExOrEyMissComponentShift) & componentMask;
0180
0181     if (of || (abs(ex) >= maxInput) || (abs(ey) >= maxInput)) {
0182       hxCompBits = maxPosComponent;
0183       hyCompBits = maxPosComponent;
0184     }
0185
0186     // Perform the table lookup to get the missing Ht magnitude and phi
0187     uint16_t lutAddress = static_cast<uint16_t>((hxCompBits << L1GctHtMissLut::kHxOrHyMissComponentNBits) | hyCompBits);
0188
0189     uint16_t lutData = m_htMissLut->lutValue(lutAddress);
0190
0191     result.mag = static_cast<unsigned>(lutData >> L1GctHtMissLut::kHtMissAngleNBits) & resultMagMask;
0192     result.phi = static_cast<unsigned>(lutData) & resultPhiMask;
0193   }
0194
0195   return result;
0196 }
0197
0198 L1GctMet::etmiss_internal L1GctMet::oldGctAlgo(const int ex, const int ey) const {
0199   //---------------------------------------------------------------------------------
0200   //
0201   // Calculates magnitude and direction of missing Et, given measured Ex and Ey.
0202   //
0203   // The algorithm used is suitable for implementation in hardware, using integer
0204   // multiplication, addition and comparison and bit shifting operations.
0205   //
0206   // Proceed in two stages. The first stage gives a result that lies between
0207   // 92% and 100% of the true Et, with the direction measured in 45 degree bins.
0208   // The final precision depends on the number of factors used in corrFact.
0209   // The present version with eleven factors gives a precision of 1% on Et, and
0210   // finds the direction to the nearest 5 degrees.
0211   //
0212   //---------------------------------------------------------------------------------
0213   etmiss_internal result;
0214
0215   unsigned eneCoarse, phiCoarse;
0216   unsigned eneCorect, phiCorect;
0217
0218   const unsigned root2fact = 181;
0219   const unsigned corrFact[11] = {24, 39, 51, 60, 69, 77, 83, 89, 95, 101, 106};
0220   const unsigned corrDphi[11] = {0, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4};
0221
0222   std::vector<bool> s(3);
0223   unsigned Mx, My, Mw;
0224
0225   unsigned Dx, Dy;
0226   unsigned eFact;
0227
0228   unsigned b, phibin;
0229   bool midphi = false;
0230
0231   // Here's the coarse calculation, with just one multiply operation
0232   //
0233   My = static_cast<unsigned>(abs(ey));
0234   Mx = static_cast<unsigned>(abs(ex));
0235   Mw = (((Mx + My) * root2fact) + 0x80) >> 8;
0236
0237   s.at(0) = (ey < 0);
0238   s.at(1) = (ex < 0);
0239   s.at(2) = (My > Mx);
0240
0241   phibin = 0;
0242   b = 0;
0243   for (int i = 0; i < 3; i++) {
0244     if (s.at(i)) {
0245       b = 1 - b;
0246     }
0247     phibin = 2 * phibin + b;
0248   }
0249
0250   eneCoarse = std::max(std::max(Mx, My), Mw);
0251   phiCoarse = phibin * 9;
0252
0253   // For the fine calculation we multiply both input components
0254   // by all the factors in the corrFact list in order to find
0255   // the required corrections to the energy and angle
0256   //
0257   for (eFact = 0; eFact < 10; eFact++) {
0258     Dx = (Mx * corrFact[eFact]) >> 8;
0259     Dy = (My * corrFact[eFact]) >> 8;
0260     if ((Dx >= My) || (Dy >= Mx)) {
0261       midphi = false;
0262       break;
0263     }
0264     if ((Mx + Dx) >= (My - Dy) && (My + Dy) >= (Mx - Dx)) {
0265       midphi = true;
0266       break;
0267     }
0268   }
0269   eneCorect = (eneCoarse * (128 + eFact)) >> 7;
0270   if (midphi ^ (b == 1)) {
0271     phiCorect = phiCoarse + 8 - corrDphi[eFact];
0272   } else {
0273     phiCorect = phiCoarse + corrDphi[eFact];
0274   }
0275
0276   // Store the result of the calculation
0277   //
0278   result.mag = eneCorect;
0279   result.phi = phiCorect;
0280
0281   return result;
0282 }
0283
0284 L1GctMet::etmiss_internal L1GctMet::floatingPointAlgo(const int ex, const int ey) const {
0285   etmiss_internal result;
0286
0287   double fx = static_cast<double>(ex);
0288   double fy = static_cast<double>(ey);
0289   double fmag = sqrt(fx * fx + fy * fy);
0290   double fphi = 36. * atan2(fy, fx) / M_PI;
0291
0292   result.mag = static_cast<unsigned>(fmag);
0293   if (fphi >= 0) {
0294     result.phi = static_cast<unsigned>(fphi);
0295   } else {
0296     result.phi = static_cast<unsigned>(fphi + 72.);
0297   }
0298
0299   return result;
0300 }
0301
0302 void L1GctMet::setEtScale(const L1CaloEtScale* const fn) { m_htMissLut->setEtScale(fn); }
0303
0304 void L1GctMet::setEtComponentLsb(const double lsb) {
0305   m_htMissLut->setExEyLsb(lsb * static_cast<double>(1 << kExOrEyMissComponentShift));
0306 }
0307
0308 const L1CaloEtScale* L1GctMet::etScale() const { return m_htMissLut->etScale(); }
0309
0310 const double L1GctMet::componentLsb() const { return m_htMissLut->componentLsb(); }
0311
0312 /// Private method to check for an overflow condition on the input components
0313 /// Allows the check to depend on the algorithm type
0314 const bool L1GctMet::inputOverFlow() const {
0315   bool result = m_exComponent.overFlow() || m_eyComponent.overFlow();
0316
0317   if (m_algoType == useHtMissLut) {
0318     static const int maxComponentInput =
0319         (1 << (L1GctHtMissLut::kHxOrHyMissComponentNBits + kExOrEyMissComponentShift - 1)) - 1;
0320
0321     // Emulate the (symmetric) overflow condition used in the firmware
0322     result |= (m_exComponent.value() > maxComponentInput) || (m_exComponent.value() < -maxComponentInput) ||
0323               (m_eyComponent.value() > maxComponentInput) || (m_eyComponent.value() < -maxComponentInput);
0324   }
0325
0326   return result;
0327 }