src/firmware/Stage1Layer2TauAlgorithmImpHW.cc

0001 /// \class l1t::Stage1Layer2TauAlgorithm
0002 ///
0003 /// Description: interface for MP firmware
0004 ///
0005 /// Implementation:
0006 ///
0007 /// \author: Kalanand Mishra - Fermilab
0008 ///
0009 /// Tau definition: 4x8 towers.
0010
0011 #include "L1Trigger/L1TCalorimeter/interface/Stage1Layer2TauAlgorithmImp.h"
0012 #include "DataFormats/L1TCalorimeter/interface/CaloRegion.h"
0013 #include "DataFormats/L1CaloTrigger/interface/L1CaloRegionDetId.h"
0014 #include "L1Trigger/L1TCalorimeter/interface/PUSubtractionMethods.h"
0015 #include "L1Trigger/L1TCalorimeter/interface/JetFinderMethods.h"
0016 #include "L1Trigger/L1TCalorimeter/interface/legacyGtHelper.h"
0017 #include "L1Trigger/L1TCalorimeter/interface/HardwareSortingMethods.h"
0018
0019 using namespace std;
0020 using namespace l1t;
0021
0022 Stage1Layer2TauAlgorithmImpHW::Stage1Layer2TauAlgorithmImpHW(CaloParamsHelper const* params)
0023     : params_(params), isoTauLut{make_unique<Stage1TauIsolationLUT>(params_)} {}
0024
0025 void l1t::Stage1Layer2TauAlgorithmImpHW::processEvent(const std::vector<l1t::CaloEmCand>& EMCands,
0026                                                       const std::vector<l1t::CaloRegion>& regions,
0027                                                       std::vector<l1t::Tau>* isoTaus,
0028                                                       std::vector<l1t::Tau>* taus) {
0029   double towerLsb = params_->towerLsbSum();
0030
0031   int tauSeedThreshold = floor(params_->tauSeedThreshold() / towerLsb + 0.5);            // convert GeV to HW units
0032   int tauNeighbourThreshold = floor(params_->tauNeighbourThreshold() / towerLsb + 0.5);  // convert GeV to HW units
0033   int tauMaxPtTauVeto = floor(params_->tauMaxPtTauVeto() / towerLsb + 0.5);
0034   int isoTauEtaMin = params_->isoTauEtaMin();
0035   int isoTauEtaMax = params_->isoTauEtaMax();
0036
0037   std::vector<l1t::CaloRegion> subRegions;
0038
0039   //Region Correction will return uncorrected subregions if
0040   //regionPUSType is set to None in the config
0041   RegionCorrection(regions, &subRegions, params_);
0042
0043   // ----- need to cluster jets in order to compute jet isolation ----
0044   std::vector<l1t::Jet> unCorrJets;
0045   TwelveByTwelveFinder(0, &subRegions, &unCorrJets);
0046
0047   std::vector<l1t::Tau> preGtTaus;
0048   std::vector<l1t::Tau> preSortTaus;
0049   std::vector<l1t::Tau> sortedTaus;
0050   std::vector<l1t::Tau> preGtIsoTaus;
0051   std::vector<l1t::Tau> preSortIsoTaus;
0052   std::vector<l1t::Tau> sortedIsoTaus;
0053
0054   for (CaloRegionBxCollection::const_iterator region = subRegions.begin(); region != subRegions.end(); region++) {
0055     if (region->hwEta() < 4 || region->hwEta() > 17)
0056       continue;
0057
0058     int regionEt = region->hwPt();
0059     if (regionEt < tauSeedThreshold)
0060       continue;
0061
0062     int regionEta = region->hwEta();
0063     int regionPhi = region->hwPhi();
0064
0065     int tauEt = regionEt;
0066     int isoFlag = 0;  // is 1 if it passes the relative jet iso requirement
0067     int quality = 0;  //doesn't really mean anything and isn't used
0068
0069     int highestNeighborEt = -1;
0070     int highestNeighborEta = -1;
0071     int highestNeighborPhi = -1;
0072     int highestNeighborTauVeto = -1;
0073
0074     //Find neighbor with highest Et
0075     for (CaloRegionBxCollection::const_iterator neighbor = subRegions.begin(); neighbor != subRegions.end();
0076          neighbor++) {
0077       int neighborPhi = neighbor->hwPhi();
0078       int neighborEta = neighbor->hwEta();
0079       if (neighborEta < 4 || neighborEta > 17)
0080         continue;
0081
0082       int deltaPhi = regionPhi - neighborPhi;
0083       if (std::abs(deltaPhi) == L1CaloRegionDetId::N_PHI - 1)
0084         deltaPhi = -deltaPhi / std::abs(deltaPhi);  //18 regions in phi
0085
0086       deltaPhi = std::abs(deltaPhi);
0087       int deltaEta = std::abs(regionEta - neighborEta);
0088
0089       if (deltaPhi + deltaEta > 0 && deltaPhi + deltaEta < 2) {  //nondiagonal neighbors
0090         if (neighbor->hwPt() > highestNeighborEt) {
0091           highestNeighborEt = neighbor->hwPt();
0092           highestNeighborEta = neighbor->hwEta();
0093           highestNeighborPhi = neighbor->hwPhi();
0094           highestNeighborTauVeto = neighbor->hwQual() & 0x1;  // tauVeto should be the first bit of quality integer
0095         }
0096       }
0097     }
0098
0099     string NESW = findNESW(regionEta, regionPhi, highestNeighborEta, highestNeighborPhi);
0100
0101     if ((tauEt > highestNeighborEt && (NESW == "isEast" || NESW == "isNorth")) ||
0102         (tauEt >= highestNeighborEt && (NESW == "isSouth" || NESW == "isWest")) || highestNeighborEt == 0) {
0103       if (highestNeighborEt >= tauNeighbourThreshold)
0104         tauEt += highestNeighborEt;
0105
0106       int regionTauVeto = region->hwQual() & 0x1;  // tauVeto should be the first bit of quality integer
0107
0108       // compute relative jet isolation
0109       if (region->hwEta() >= isoTauEtaMin && region->hwEta() <= isoTauEtaMax) {
0110         if ((highestNeighborTauVeto == 0 && regionTauVeto == 0) || tauEt > tauMaxPtTauVeto) {
0111           int jetEt = AssociatedJetPt(region->hwEta(), region->hwPhi(), &unCorrJets);
0112           if (jetEt > 0) {
0113             unsigned int MAX_LUT_ADDRESS = params_->tauIsolationLUT()->maxSize() - 1;
0114             // unsigned int lutAddress = isoLutIndex(tauEt,jetEt);
0115             unsigned lutAddress = isoTauLut->lutAddress(tauEt, jetEt);
0116             if (tauEt > 0) {
0117               if (lutAddress > MAX_LUT_ADDRESS)
0118                 lutAddress = MAX_LUT_ADDRESS;
0119               isoFlag = params_->tauIsolationLUT()->data(lutAddress);
0120             }
0121
0122           } else {  // no associated jet
0123             isoFlag = 1;
0124           }
0125         }
0126       }
0127
0128       ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > tauLorentz(0, 0, 0, 0);
0129
0130       l1t::Tau theTau(*&tauLorentz, tauEt, region->hwEta(), region->hwPhi(), quality, isoFlag);
0131
0132       preGtTaus.push_back(theTau);
0133       if (isoFlag)
0134         preGtIsoTaus.push_back(theTau);
0135     }
0136   }
0137
0138   // TauToGtPtScales(params_, preGtTaus, preRankTaus);
0139   // TauToGtPtScales(params_, preGtIsoTaus, preRankIsoTaus);
0140   calibrateAndRankTaus(params_, &preGtTaus, &preSortTaus);
0141   calibrateAndRankTaus(params_, &preGtIsoTaus, &preSortIsoTaus);
0142
0143   SortTaus(&preSortTaus, &sortedTaus);
0144   SortTaus(&preSortIsoTaus, &sortedIsoTaus);
0145
0146   TauToGtEtaScales(params_, &sortedTaus, taus);
0147   TauToGtEtaScales(params_, &sortedIsoTaus, isoTaus);
0148
0149   // set all filler taus to have isolation bit set
0150   for (std::vector<l1t::Tau>::iterator iTau = isoTaus->begin(); iTau != isoTaus->end(); ++iTau)
0151     iTau->setHwIso(1);
0152 }
0153
0154 //  Compute jet isolation.
0155 double l1t::Stage1Layer2TauAlgorithmImpHW::JetIsolation(int et,
0156                                                         int ieta,
0157                                                         int iphi,
0158                                                         const std::vector<l1t::Jet>& jets) const {
0159   for (JetBxCollection::const_iterator jet = jets.begin(); jet != jets.end(); jet++) {
0160     if (ieta == jet->hwEta() && iphi == jet->hwPhi()) {
0161       double isolation = (double)(jet->hwPt() - et);
0162       return isolation / et;
0163     }
0164   }
0165
0166   // set output
0167   return 999.;
0168 }
0169
0170 //  Find if the neighbor with the highest Et is N, E, S, or W
0171 string l1t::Stage1Layer2TauAlgorithmImpHW::findNESW(int ieta, int iphi, int neta, int nphi) const {
0172   int deltaPhi = iphi - nphi;
0173   if (std::abs(deltaPhi) == L1CaloRegionDetId::N_PHI - 1)
0174     deltaPhi = -deltaPhi / std::abs(deltaPhi);  //18 regions in phi
0175
0176   int deltaEta = ieta - neta;
0177
0178   if ((std::abs(deltaPhi) + std::abs(deltaEta)) < 2) {
0179     if (deltaEta == -1) {
0180       return "isEast";
0181     } else if (deltaEta == 0) {
0182       if (deltaPhi == -1) {
0183         return "isNorth";
0184       }
0185       if (deltaPhi == 1) {
0186         return "isSouth";
0187       }
0188     } else {
0189       return "isWest";
0190     }
0191   }
0192
0193   return "999";
0194 }
0195
0196 int l1t::Stage1Layer2TauAlgorithmImpHW::AssociatedJetPt(int ieta, int iphi, const std::vector<l1t::Jet>* jets) const {
0197   int pt = -1;
0198
0199   for (JetBxCollection::const_iterator itJet = jets->begin(); itJet != jets->end(); ++itJet) {
0200     int jetEta = itJet->hwEta();
0201     int jetPhi = itJet->hwPhi();
0202     if ((jetEta == ieta) && (jetPhi == iphi)) {
0203       pt = itJet->hwPt();
0204       break;
0205     }
0206   }
0207
0208   // set output
0209   return pt;
0210 }
0211
0212 unsigned l1t::Stage1Layer2TauAlgorithmImpHW::isoLutIndex(unsigned int tauPt, unsigned int jetPt) const {
0213   //const unsigned int nbitsTau=9;  // number of bits used for et in LUT file (needed for left shift operation)
0214   //const unsigned int nbitsJet=9;
0215
0216   const unsigned int nbitsTau = 8;  // number of bits used for et in LUT file (needed for left shift operation)
0217   const unsigned int nbitsJet = 8;
0218
0219   const unsigned int maxJet = pow(2, nbitsJet) - 1;
0220   const unsigned int maxTau = pow(2, nbitsTau) - 1;
0221
0222   if (nbitsTau < 9) {
0223     if (nbitsTau == 6) {
0224       tauPt = tauPt >> 3;
0225     } else if (nbitsTau == 7) {
0226       tauPt = tauPt >> 2;
0227     } else if (nbitsTau == 8) {
0228       tauPt = tauPt >> 1;
0229     }
0230   }
0231
0232   if (nbitsJet < 9)  // no need to do shift if nbits>=9
0233   {
0234     if (nbitsJet == 6) {
0235       jetPt = jetPt >> 3;
0236     } else if (nbitsJet == 7) {
0237       jetPt = jetPt >> 2;
0238     } else if (nbitsJet == 8) {
0239       jetPt = jetPt >> 1;
0240     }
0241   }
0242
0243   if (jetPt > maxJet)
0244     jetPt = maxJet;
0245   if (tauPt > maxTau)
0246     tauPt = maxTau;
0247
0248   unsigned int address = (jetPt << nbitsTau) + tauPt;
0249
0250   return address;
0251 }