src/firmware/Stage1Layer2EtSumAlgorithmImpHW.cc

0001 ///
0002 /// \class l1t::Stage1Layer2EtSumAlgorithmImpHW
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::Stage1Layer2EtSumAlgorithmImpHW::Stage1Layer2EtSumAlgorithmImpHW(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::Stage1Layer2EtSumAlgorithmImpHW::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   std::vector<l1t::CaloRegion> subRegions;
0036
0037   //Region Correction will return uncorrected subregions if
0038   //regionPUSType is set to None in the config
0039   double jetLsb = params_->jetLsb();
0040
0041   int etSumEtaMinEt = params_->etSumEtaMin(0);
0042   int etSumEtaMaxEt = params_->etSumEtaMax(0);
0043   //double etSumEtThresholdEt = params_->etSumEtThreshold(0);
0044   int etSumEtThresholdEt = (int)(params_->etSumEtThreshold(0) / jetLsb);
0045
0046   int etSumEtaMinHt = params_->etSumEtaMin(1);
0047   int etSumEtaMaxHt = params_->etSumEtaMax(1);
0048   //double etSumEtThresholdHt = params_->etSumEtThreshold(1);
0049   int etSumEtThresholdHt = (int)(params_->etSumEtThreshold(1) / jetLsb);
0050
0051   RegionCorrection(regions, &subRegions, params_);
0052
0053   std::vector<SimpleRegion> regionEtVect;
0054   std::vector<SimpleRegion> regionHtVect;
0055
0056   // check the un-subtracted regions for overflow
0057   bool regionOverflowEt(false);
0058   bool regionOverflowHt(false);
0059   for (auto& region : regions) {
0060     if (region.hwEta() >= etSumEtaMinEt && region.hwEta() <= etSumEtaMaxEt) {
0061       if (region.hwPt() >= 1023) {
0062         regionOverflowEt = true;
0063       }
0064     }
0065     if (region.hwEta() >= etSumEtaMinHt && region.hwEta() <= etSumEtaMaxHt) {
0066       if (region.hwPt() >= 1023) {
0067         regionOverflowHt = true;
0068       }
0069     }
0070   }
0071
0072   // hwPt() is the sum ET+HT in region, for stage 1 this will be
0073   // the region sum input to MET algorithm
0074   // In stage 2, we would move to hwEtEm() and hwEtHad() for separate MET/MHT
0075   // Thresholds will be hardware values not physical
0076   for (auto& region : subRegions) {
0077     if (region.hwEta() >= etSumEtaMinEt && region.hwEta() <= etSumEtaMaxEt) {
0078       if (region.hwPt() >= etSumEtThresholdEt) {
0079         SimpleRegion r;
0080         r.ieta = region.hwEta();
0081         r.iphi = region.hwPhi();
0082         r.et = region.hwPt();
0083         regionEtVect.push_back(r);
0084       }
0085     }
0086     if (region.hwEta() >= etSumEtaMinHt && region.hwEta() <= etSumEtaMaxHt) {
0087       if (region.hwPt() >= etSumEtThresholdHt) {
0088         SimpleRegion r;
0089         r.ieta = region.hwEta();
0090         r.iphi = region.hwPhi();
0091         r.et = region.hwPt();
0092         regionHtVect.push_back(r);
0093       }
0094     }
0095   }
0096
0097   int sumET, MET, iPhiET;
0098   std::tie(sumET, MET, iPhiET) = doSumAndMET(regionEtVect, ETSumType::kEmSum);
0099
0100   int sumHT, MHT, iPhiHT;
0101   std::tie(sumHT, MHT, iPhiHT) = doSumAndMET(regionHtVect, ETSumType::kHadronicSum);
0102
0103   // Set quality (i.e. overflow) bits appropriately
0104   int METqual = 0;
0105   int MHTqual = 0;
0106   int ETTqual = 0;
0107   int HTTqual = 0;
0108   if (MET > 0xfff || regionOverflowEt)  // MET 12 bits
0109     METqual = 1;
0110   if (MHT > 0x7f || regionOverflowHt)  // MHT 7 bits
0111     MHTqual = 1;
0112   if (sumET > 0xfff || regionOverflowEt)
0113     ETTqual = 1;
0114   if (sumHT > 0xfff || regionOverflowHt)
0115     HTTqual = 1;
0116
0117   MHT &= 127;  // limit MHT to 7 bits as the firmware does, but only after checking for overflow.
0118   //MHT is replaced with MHT/HT
0119   uint16_t MHToHT = MHToverHT(MHT, sumHT);
0120   //iPhiHt is replaced by the dPhi between two most energetic jets
0121   iPhiHT = DiJetPhi(jets);
0122
0123   const ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > etLorentz(0, 0, 0, 0);
0124   l1t::EtSum etMiss(*&etLorentz, EtSum::EtSumType::kMissingEt, MET & 0xfff, 0, iPhiET, METqual);
0125   l1t::EtSum htMiss(*&etLorentz, EtSum::EtSumType::kMissingHt, MHToHT & 0x7f, 0, iPhiHT, MHTqual);
0126   l1t::EtSum etTot(*&etLorentz, EtSum::EtSumType::kTotalEt, sumET & 0xfff, 0, 0, ETTqual);
0127   l1t::EtSum htTot(*&etLorentz, EtSum::EtSumType::kTotalHt, sumHT & 0xfff, 0, 0, HTTqual);
0128
0129   std::vector<l1t::EtSum> preGtEtSums = {etMiss, htMiss, etTot, htTot};
0130
0131   EtSumToGtScales(params_, &preGtEtSums, etsums);
0132 }
0133
0134 std::tuple<int, int, int> l1t::Stage1Layer2EtSumAlgorithmImpHW::doSumAndMET(const std::vector<SimpleRegion>& regionEt,
0135                                                                             ETSumType sumType) {
0136   std::array<int, 18> sumEtaPos{};
0137   std::array<int, 18> sumEtaNeg{};
0138   for (const auto& r : regionEt) {
0139     if (r.ieta < 11)
0140       sumEtaNeg[r.iphi] += r.et;
0141     else
0142       sumEtaPos[r.iphi] += r.et;
0143   }
0144
0145   std::array<int, 18> sumEta{};
0146   int sumEt(0);
0147   for (size_t i = 0; i < sumEta.size(); ++i) {
0148     assert(sumEtaPos[i] >= 0 && sumEtaNeg[i] >= 0);
0149     sumEta[i] = sumEtaPos[i] + sumEtaNeg[i];
0150     sumEt += sumEta[i];
0151   }
0152
0153   // 0, 20, 40, 60, 80 degrees
0154   std::array<int, 5> sumsForCos{};
0155   std::array<int, 5> sumsForSin{};
0156   for (size_t iphi = 0; iphi < sumEta.size(); ++iphi) {
0157     if (iphi < 5) {
0158       sumsForCos[iphi] += sumEta[iphi];
0159       sumsForSin[iphi] += sumEta[iphi];
0160     } else if (iphi < 9) {
0161       sumsForCos[9 - iphi] -= sumEta[iphi];
0162       sumsForSin[9 - iphi] += sumEta[iphi];
0163     } else if (iphi < 14) {
0164       sumsForCos[iphi - 9] -= sumEta[iphi];
0165       sumsForSin[iphi - 9] -= sumEta[iphi];
0166     } else {
0167       sumsForCos[18 - iphi] += sumEta[iphi];
0168       sumsForSin[18 - iphi] -= sumEta[iphi];
0169     }
0170   }
0171
0172   long sumX(0l);
0173   long sumY(0l);
0174   for (int i = 0; i < 5; ++i) {
0175     sumX += sumsForCos[i] * cosines[i];
0176     sumY += sumsForSin[i] * sines[i];
0177   }
0178   assert(abs(sumX) < (1l << 48) && abs(sumY) < (1l << 48));
0179   int cordicX = sumX >> 25;
0180   int cordicY = sumY >> 25;
0181
0182   uint32_t cordicMag(0);
0183   int cordicPhase(0);
0184   cordic(cordicX, cordicY, cordicPhase, cordicMag);
0185
0186   int met(0);
0187   int metPhi(0);
0188   if (sumType == ETSumType::kHadronicSum) {
0189     met = (cordicMag % (1 << 7)) | ((cordicMag >= (1 << 7)) ? (1 << 7) : 0);
0190     metPhi = cordicToMETPhi(cordicPhase) >> 2;
0191     assert(metPhi >= 0 && metPhi < 18);
0192   } else {
0193     met = (cordicMag % (1 << 12)) | ((cordicMag >= (1 << 12)) ? (1 << 12) : 0);
0194     metPhi = cordicToMETPhi(cordicPhase);
0195     assert(metPhi >= 0 && metPhi < 72);
0196   }
0197
0198   return std::make_tuple(sumEt, met, metPhi);
0199 }
0200
0201 // converts phase from 3Q16 to 0-71
0202 // Expects abs(phase) <= 205887 (pi*2^16)
0203 int l1t::Stage1Layer2EtSumAlgorithmImpHW::cordicToMETPhi(int phase) {
0204   assert(abs(phase) <= 205887);
0205   for (size_t i = 0; i < cordicPhiValues.size() - 1; ++i)
0206     if (phase >= cordicPhiValues[i] && phase < cordicPhiValues[i + 1])
0207       return i;
0208   // if phase == +205887 (+pi), return zero
0209   return 0;
0210 }
0211
0212 int l1t::Stage1Layer2EtSumAlgorithmImpHW::DiJetPhi(const std::vector<l1t::Jet>* jets) const {
0213   int dphi = 10;  // initialize to negative physical dphi value
0214   if (jets->size() < 2)
0215     return dphi;  // size() not really reliable as we pad the size to 8 (4cen+4for) in the sorter
0216   if ((*jets).at(0).hwPt() == 0)
0217     return dphi;
0218   if ((*jets).at(1).hwPt() == 0)
0219     return dphi;
0220
0221   int iphi1 = (*jets).at(0).hwPhi();
0222   int iphi2 = (*jets).at(1).hwPhi();
0223
0224   int difference = abs(iphi1 - iphi2);
0225
0226   if (difference > 9)
0227     difference = L1CaloRegionDetId::N_PHI - difference;  // make Physical dphi always positive
0228   return difference;
0229 }
0230
0231 uint16_t l1t::Stage1Layer2EtSumAlgorithmImpHW::MHToverHT(uint16_t num, uint16_t den) const {
0232   uint16_t result;
0233   uint32_t numerator(num), denominator(den);
0234
0235   if (numerator == denominator)
0236     result = 0x7f;
0237   else {
0238     numerator = numerator << 7;
0239     result = numerator / denominator;
0240     result = result & 0x7f;
0241   }
0242   return result;
0243 }