Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:20:24 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!