src/firmware/Stage1Layer2EtSumAlgorithmImpHI.cc

0001 ///
0002 /// \class l1t::Stage1Layer2EtSumAlgorithmImpHI
0003 ///
0004 /// \author: Nick Smith (nick.smith@cern.ch)
0005 ///
0006 /// Description: hardware emulation of et sum algorithm
0007
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "L1Trigger/L1TCalorimeter/interface/Stage1Layer2EtSumAlgorithmImp.h"
0010 #include "L1Trigger/L1TCalorimeter/interface/PUSubtractionMethods.h"
0011 #include "L1Trigger/L1TCalorimeter/interface/legacyGtHelper.h"
0012 #include "DataFormats/L1Trigger/interface/EtSum.h"
0013 #include "DataFormats/L1CaloTrigger/interface/L1CaloRegionDetId.h"
0014 #include "L1Trigger/L1TCalorimeter/interface/JetFinderMethods.h"
0015 #include "L1Trigger/L1TCalorimeter/interface/JetCalibrationMethods.h"
0016 #include "L1Trigger/L1TCalorimeter/interface/HardwareSortingMethods.h"
0017 #include <cassert>
0018
0019 l1t::Stage1Layer2EtSumAlgorithmImpHI::Stage1Layer2EtSumAlgorithmImpHI(CaloParamsHelper const* params)
0020     : params_(params) {
0021   //now do what ever initialization is needed
0022   for (size_t i = 0; i < cordicPhiValues.size(); ++i) {
0023     cordicPhiValues[i] = static_cast<int>(pow(2., 16) * (((float)i) - 36) * M_PI / 36);
0024   }
0025   for (size_t i = 0; i < sines.size(); ++i) {
0026     sines[i] = static_cast<long>(pow(2, 30) * sin(i * 20 * M_PI / 180));
0027     cosines[i] = static_cast<long>(pow(2, 30) * cos(i * 20 * M_PI / 180));
0028   }
0029 }
0030
0031 void l1t::Stage1Layer2EtSumAlgorithmImpHI::processEvent(const std::vector<l1t::CaloRegion>& regions,
0032                                                         const std::vector<l1t::CaloEmCand>& EMCands,
0033                                                         const std::vector<l1t::Jet>* jets,
0034                                                         std::vector<l1t::EtSum>* etsums) {
0035   //Region Correction will return uncorrected subregions if
0036   //regionPUSType is set to None in the config
0037   //double jetLsb=params_->jetLsb();
0038   double jetLsb = 0.5;  // HI O2O does not set this, and it will never change.
0039
0040   //int etSumEtaMinEt = params_->etSumEtaMin(0);
0041   //int etSumEtaMaxEt = params_->etSumEtaMax(0);
0042   //double etSumEtThresholdEt = params_->etSumEtThreshold(0);
0043   int etSumEtThresholdEt = (int)(params_->etSumEtThreshold(0) / jetLsb);
0044
0045   //int etSumEtaMinHt = params_->etSumEtaMin(1);
0046   //int etSumEtaMaxHt = params_->etSumEtaMax(1);
0047   //double etSumEtThresholdHt = params_->etSumEtThreshold(1);
0048   int etSumEtThresholdHt = (int)(params_->etSumEtThreshold(1) / jetLsb);
0049
0050   // These values are not changeable online. O2O code for HI does not set a default,
0051   // so previous results were garbage.
0052   // Boundaries of 4 and 17 correspond to all non-HF regions.
0053   int etSumEtaMinEt = 4;
0054   int etSumEtaMaxEt = 17;
0055   int etSumEtaMinHt = 4;
0056   int etSumEtaMaxHt = 17;
0057
0058   //RegionCorrection(regions, subRegions, params_);
0059
0060   std::vector<SimpleRegion> regionEtVect;
0061   std::vector<SimpleRegion> regionHtVect;
0062
0063   // check the un-subtracted regions for overflow
0064   bool regionOverflowEt(false);
0065   bool regionOverflowHt(false);
0066   for (auto& region : regions) {
0067     if (region.hwEta() >= etSumEtaMinEt && region.hwEta() <= etSumEtaMaxEt) {
0068       if (region.hwPt() >= 1023) {
0069         regionOverflowEt = true;
0070       }
0071     }
0072     if (region.hwEta() >= etSumEtaMinHt && region.hwEta() <= etSumEtaMaxHt) {
0073       if (region.hwPt() >= 1023) {
0074         regionOverflowHt = true;
0075       }
0076     }
0077   }
0078
0079   // hwPt() is the sum ET+HT in region, for stage 1 this will be
0080   // the region sum input to MET algorithm
0081   // In stage 2, we would move to hwEtEm() and hwEtHad() for separate MET/MHT
0082   // Thresholds will be hardware values not physical
0083   //for (auto& region : *subRegions) {
0084   for (auto& region : regions) {
0085     if (region.hwEta() >= etSumEtaMinEt && region.hwEta() <= etSumEtaMaxEt) {
0086       if (region.hwPt() >= etSumEtThresholdEt) {
0087         SimpleRegion r;
0088         r.ieta = region.hwEta();
0089         r.iphi = region.hwPhi();
0090         r.et = region.hwPt();
0091         regionEtVect.push_back(r);
0092       }
0093     }
0094     if (region.hwEta() >= etSumEtaMinHt && region.hwEta() <= etSumEtaMaxHt) {
0095       if (region.hwPt() >= etSumEtThresholdHt) {
0096         SimpleRegion r;
0097         r.ieta = region.hwEta();
0098         r.iphi = region.hwPhi();
0099         r.et = region.hwPt();
0100         regionHtVect.push_back(r);
0101       }
0102     }
0103   }
0104
0105   int sumET, MET, iPhiET;
0106   std::tie(sumET, MET, iPhiET) = doSumAndMET(regionEtVect, ETSumType::kEmSum);
0107
0108   int sumHT, MHT, iPhiHT;
0109   std::tie(sumHT, MHT, iPhiHT) = doSumAndMET(regionHtVect, ETSumType::kHadronicSum);
0110
0111   // Set quality (i.e. overflow) bits appropriately
0112   int METqual = 0;
0113   int MHTqual = 0;
0114   int ETTqual = 0;
0115   int HTTqual = 0;
0116   if (MET >= 0xfff || regionOverflowEt)  // MET 12 bits
0117     METqual = 1;
0118   if (MHT >= 0x7f || regionOverflowHt)  // MHT 7 bits
0119     MHTqual = 1;
0120   if (sumET >= 0xfff || regionOverflowEt)
0121     ETTqual = 1;
0122   if (sumHT >= 0xfff || regionOverflowHt)
0123     HTTqual = 1;
0124
0125   MHT &= 127;  // limit MHT to 7 bits as the firmware does, but only after checking for overflow.
0126   //MHT is replaced with MHT/HT
0127   uint16_t MHToHT = MHToverHT(MHT, sumHT);
0128   //iPhiHt is replaced by the dPhi between two most energetic jets
0129   iPhiHT = DiJetPhi(jets);
0130
0131   const ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > etLorentz(0, 0, 0, 0);
0132   l1t::EtSum etMiss(*&etLorentz, EtSum::EtSumType::kMissingEt, MET & 0xfff, 0, iPhiET, METqual);
0133   l1t::EtSum htMiss(*&etLorentz, EtSum::EtSumType::kMissingHt, MHToHT & 0x7f, 0, iPhiHT, MHTqual);
0134   l1t::EtSum etTot(*&etLorentz, EtSum::EtSumType::kTotalEt, sumET & 0xfff, 0, 0, ETTqual);
0135   l1t::EtSum htTot(*&etLorentz, EtSum::EtSumType::kTotalHt, sumHT & 0xfff, 0, 0, HTTqual);
0136
0137   std::vector<l1t::EtSum> preGtEtSums;
0138   preGtEtSums.reserve(4);
0139
0140   preGtEtSums.push_back(etMiss);
0141   preGtEtSums.push_back(htMiss);
0142   preGtEtSums.push_back(etTot);
0143   preGtEtSums.push_back(htTot);
0144
0145   EtSumToGtScales(params_, &preGtEtSums, etsums);
0146 }
0147
0148 std::tuple<int, int, int> l1t::Stage1Layer2EtSumAlgorithmImpHI::doSumAndMET(const std::vector<SimpleRegion>& regionEt,
0149                                                                             ETSumType sumType) {
0150   std::array<int, 18> sumEtaPos{};
0151   std::array<int, 18> sumEtaNeg{};
0152   for (const auto& r : regionEt) {
0153     if (r.ieta < 11)
0154       sumEtaNeg[r.iphi] += r.et;
0155     else
0156       sumEtaPos[r.iphi] += r.et;
0157   }
0158
0159   std::array<int, 18> sumEta{};
0160   int sumEt(0);
0161   for (size_t i = 0; i < sumEta.size(); ++i) {
0162     assert(sumEtaPos[i] >= 0 && sumEtaNeg[i] >= 0);
0163     sumEta[i] = sumEtaPos[i] + sumEtaNeg[i];
0164     sumEt += sumEta[i];
0165   }
0166
0167   // 0, 20, 40, 60, 80 degrees
0168   std::array<int, 5> sumsForCos{};
0169   std::array<int, 5> sumsForSin{};
0170   for (size_t iphi = 0; iphi < sumEta.size(); ++iphi) {
0171     if (iphi < 5) {
0172       sumsForCos[iphi] += sumEta[iphi];
0173       sumsForSin[iphi] += sumEta[iphi];
0174     } else if (iphi < 9) {
0175       sumsForCos[9 - iphi] -= sumEta[iphi];
0176       sumsForSin[9 - iphi] += sumEta[iphi];
0177     } else if (iphi < 14) {
0178       sumsForCos[iphi - 9] -= sumEta[iphi];
0179       sumsForSin[iphi - 9] -= sumEta[iphi];
0180     } else {
0181       sumsForCos[18 - iphi] += sumEta[iphi];
0182       sumsForSin[18 - iphi] -= sumEta[iphi];
0183     }
0184   }
0185
0186   long sumX(0l);
0187   long sumY(0l);
0188   for (int i = 0; i < 5; ++i) {
0189     sumX += sumsForCos[i] * cosines[i];
0190     sumY += sumsForSin[i] * sines[i];
0191   }
0192   assert(abs(sumX) < (1l << 48) && abs(sumY) < (1l << 48));
0193   int cordicX = sumX >> 25;
0194   int cordicY = sumY >> 25;
0195
0196   uint32_t cordicMag(0);
0197   int cordicPhase(0);
0198   cordic(cordicX, cordicY, cordicPhase, cordicMag);
0199
0200   int met(0);
0201   int metPhi(0);
0202   if (sumType == ETSumType::kHadronicSum) {
0203     met = (cordicMag % (1 << 7)) | ((cordicMag >= (1 << 7)) ? (1 << 7) : 0);
0204     metPhi = cordicToMETPhi(cordicPhase) >> 2;
0205     assert(metPhi >= 0 && metPhi < 18);
0206   } else {
0207     met = (cordicMag % (1 << 12)) | ((cordicMag >= (1 << 12)) ? (1 << 12) : 0);
0208     metPhi = cordicToMETPhi(cordicPhase);
0209     assert(metPhi >= 0 && metPhi < 72);
0210   }
0211
0212   return std::make_tuple(sumEt, met, metPhi);
0213 }
0214
0215 // converts phase from 3Q16 to 0-71
0216 // Expects abs(phase) <= 205887 (pi*2^16)
0217 int l1t::Stage1Layer2EtSumAlgorithmImpHI::cordicToMETPhi(int phase) {
0218   assert(abs(phase) <= 205887);
0219   for (size_t i = 0; i < cordicPhiValues.size() - 1; ++i)
0220     if (phase >= cordicPhiValues[i] && phase < cordicPhiValues[i + 1])
0221       return i;
0222   // if phase == +205887 (+pi), return zero
0223   return 0;
0224 }
0225
0226 int l1t::Stage1Layer2EtSumAlgorithmImpHI::DiJetPhi(const std::vector<l1t::Jet>* jets) const {
0227   int dphi = 10;  // initialize to negative physical dphi value
0228   if (jets->size() < 2)
0229     return dphi;  // size() not really reliable as we pad the size to 8 (4cen+4for) in the sorter
0230   if ((*jets).at(0).hwPt() == 0)
0231     return dphi;
0232   if ((*jets).at(1).hwPt() == 0)
0233     return dphi;
0234
0235   int iphi1 = (*jets).at(0).hwPhi();
0236   int iphi2 = (*jets).at(1).hwPhi();
0237
0238   int difference = abs(iphi1 - iphi2);
0239
0240   if (difference > 9)
0241     difference = L1CaloRegionDetId::N_PHI - difference;  // make Physical dphi always positive
0242   return difference;
0243 }
0244
0245 uint16_t l1t::Stage1Layer2EtSumAlgorithmImpHI::MHToverHT(uint16_t num, uint16_t den) const {
0246   uint16_t result;
0247   uint32_t numerator(num), denominator(den);
0248
0249   if (numerator == denominator)
0250     result = 0x7f;
0251   else {
0252     numerator = numerator << 7;
0253     result = numerator / denominator;
0254     result = result & 0x7f;
0255   }
0256
0257   return result;
0258 }