src/firmware/Stage1Layer2EGammaAlgorithmImpPP.cc

0001 ///step03
0002 /// \class l1t::Stage1Layer2EGammaAlgorithm
0003 ///
0004 /// Description: interface for MP firmware
0005 ///
0006 /// Implementation:
0007 ///
0008 /// \author: Kalanand Mishra - Fermilab
0009 ///
0010
0011 #include "L1Trigger/L1TCalorimeter/interface/Stage1Layer2EGammaAlgorithmImp.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/HardwareSortingMethods.h"
0016 #include "L1Trigger/L1TCalorimeter/interface/JetFinderMethods.h"
0017 #include "L1Trigger/L1TCalorimeter/interface/legacyGtHelper.h"
0018
0019 #include <bitset>
0020
0021 using namespace std;
0022 using namespace l1t;
0023
0024 Stage1Layer2EGammaAlgorithmImpPP::Stage1Layer2EGammaAlgorithmImpPP(CaloParamsHelper const* params) : params_(params){};
0025
0026 void l1t::Stage1Layer2EGammaAlgorithmImpPP::processEvent(const std::vector<l1t::CaloEmCand>& EMCands,
0027                                                          const std::vector<l1t::CaloRegion>& regions,
0028                                                          const std::vector<l1t::Jet>* jets,
0029                                                          std::vector<l1t::EGamma>* egammas) {
0030   double egLsb = params_->egLsb();
0031   double jetLsb = params_->jetLsb();
0032   int egSeedThreshold = floor(params_->egSeedThreshold() / egLsb + 0.5);
0033   int jetSeedThreshold = floor(params_->jetSeedThreshold() / jetLsb + 0.5);
0034   int egMinPtJetIsolation = params_->egMinPtJetIsolation();
0035   int egMaxPtJetIsolation = params_->egMaxPtJetIsolation();
0036   int egMinPtHOverEIsolation = params_->egMinPtHOverEIsolation();
0037   int egMaxPtHOverEIsolation = params_->egMaxPtHOverEIsolation();
0038
0039   std::vector<l1t::CaloRegion> subRegions;
0040   std::vector<l1t::EGamma> preSortEGammas;
0041   std::vector<l1t::EGamma> preGtEGammas;
0042
0043   //Region Correction will return uncorrected subregions if
0044   //regionPUSType is set to None in the config
0045   RegionCorrection(regions, &subRegions, params_);
0046
0047   // ----- need to cluster jets in order to compute jet isolation ----
0048   std::vector<l1t::Jet>* unCorrJets = new std::vector<l1t::Jet>();
0049   // slidingWindowJetFinder(jetSeedThreshold, subRegions, unCorrJets);
0050   TwelveByTwelveFinder(jetSeedThreshold, &subRegions, unCorrJets);
0051
0052   for (CaloEmCandBxCollection::const_iterator egCand = EMCands.begin(); egCand != EMCands.end(); egCand++) {
0053     int eg_et = egCand->hwPt();
0054     int eg_eta = egCand->hwEta();
0055     int eg_phi = egCand->hwPhi();
0056     int index = (egCand->hwIso() * 4 + egCand->hwQual());
0057
0058     if (eg_et <= egSeedThreshold)
0059       continue;
0060
0061     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > egLorentz(0, 0, 0, 0);
0062
0063     //int quality = 1;
0064     int isoFlag = 0;
0065     int isoFlagRct = 0;
0066
0067     // 3x3 HoE, computed in 3x3
0068     if (eg_et >= egMinPtHOverEIsolation && eg_et < egMaxPtHOverEIsolation) {
0069       if (egCand->hwIso())
0070         isoFlagRct = 1;
0071     } else {
0072       isoFlagRct = 1;
0073     }
0074
0075     int ijet_pt = AssociatedJetPt(eg_eta, eg_phi, unCorrJets);
0076     // double jet_pt=ijet_pt*jetLsb;
0077     bool isinBarrel = (eg_eta >= 7 && eg_eta <= 14);
0078     if (ijet_pt > 0 && eg_et >= egMinPtJetIsolation && eg_et < egMaxPtJetIsolation) {
0079       // Combined Barrel/Endcap LUT uses upper bit to indicate Barrel / Endcap:
0080       enum { MAX_LUT_ADDRESS = 0x7fff };
0081       enum { LUT_BARREL_OFFSET = 0x0, LUT_ENDCAP_OFFSET = 0x8000 };
0082
0083       unsigned int lutAddress = isoLutIndex(eg_et, ijet_pt);
0084
0085       if (eg_et > 0) {
0086         if (lutAddress > MAX_LUT_ADDRESS)
0087           lutAddress = MAX_LUT_ADDRESS;
0088         if (isinBarrel) {
0089           isoFlag = params_->egIsolationLUT()->data(LUT_BARREL_OFFSET + lutAddress);
0090         } else {
0091           isoFlag = params_->egIsolationLUT()->data(LUT_ENDCAP_OFFSET + lutAddress);
0092         }
0093       }
0094
0095     } else {  // no associated jet; assume it's an isolated eg
0096       isoFlag = 1;
0097     }
0098
0099     int fullIsoFlag = isoFlag * isoFlagRct;
0100
0101     // ------- fill the EG candidate vector ---------
0102     l1t::EGamma theEG(*&egLorentz, eg_et, eg_eta, eg_phi, index, fullIsoFlag);
0103     //?? if( hoe < HoverECut) egammas->push_back(theEG);
0104     preSortEGammas.push_back(theEG);
0105   }
0106
0107   SortEGammas(&preSortEGammas, &preGtEGammas);
0108
0109   EGammaToGtScales(params_, &preGtEGammas, egammas);
0110 }
0111
0112 /// -----  Compute isolation sum --------------------
0113 double l1t::Stage1Layer2EGammaAlgorithmImpPP::Isolation(int ieta,
0114                                                         int iphi,
0115                                                         const std::vector<l1t::CaloRegion>& regions) const {
0116   double isolation = 0;
0117
0118   for (CaloRegionBxCollection::const_iterator region = regions.begin(); region != regions.end(); region++) {
0119     int regionPhi = region->hwPhi();
0120     int regionEta = region->hwEta();
0121     ///    unsigned int deltaPhi = iphi - regionPhi;
0122     int deltaPhi = iphi - regionPhi;
0123     if (std::abs(deltaPhi) == L1CaloRegionDetId::N_PHI - 1)
0124       deltaPhi = -deltaPhi / std::abs(deltaPhi);  //18 regions in phi
0125
0126     unsigned int deltaEta = std::abs(ieta - regionEta);
0127
0128     if ((deltaPhi + deltaEta) > 0 && deltaPhi < 2 && deltaEta < 2)
0129       isolation += region->hwPt();
0130   }
0131
0132   // set output
0133   return isolation;
0134 }
0135
0136 //ieta =-28, nrTowers 0 is 0, increases to ieta28, nrTowers=kNrTowersInSum
0137 unsigned l1t::Stage1Layer2EGammaAlgorithmImpPP::isoLutIndex(unsigned int egPt, unsigned int jetPt) const {
0138   const unsigned int nbitsEG = 6;  // number of bits used for EG bins in LUT file (needed for left shift operation)
0139   //  const unsigned int nbitsJet=9; // not used but here for info  number of bits used for Jet bins in LUT file
0140
0141   unsigned int address = (jetPt << nbitsEG) + egPt;
0142   return address;
0143 }
0144
0145 int l1t::Stage1Layer2EGammaAlgorithmImpPP::AssociatedJetPt(int ieta,
0146                                                            int iphi,
0147                                                            const std::vector<l1t::Jet>* jets) const {
0148   int pt = -1;
0149
0150   for (JetBxCollection::const_iterator itJet = jets->begin(); itJet != jets->end(); ++itJet) {
0151     int jetEta = itJet->hwEta();
0152     int jetPhi = itJet->hwPhi();
0153     if ((jetEta == ieta) && (jetPhi == iphi)) {
0154       pt = itJet->hwPt();
0155       break;
0156     }
0157   }
0158
0159   // set output
0160   return pt;
0161 }
0162
0163 /// -----  Compute H/E --------------------
0164 double l1t::Stage1Layer2EGammaAlgorithmImpPP::HoverE(int et,
0165                                                      int ieta,
0166                                                      int iphi,
0167                                                      const std::vector<l1t::CaloRegion>& regions) const {
0168   int hadronicET = 0;
0169
0170   for (CaloRegionBxCollection::const_iterator region = regions.begin(); region != regions.end(); region++) {
0171     int regionET = region->hwPt();
0172     int regionPhi = region->hwPhi();
0173     int regionEta = region->hwEta();
0174
0175     if (iphi == regionPhi && ieta == regionEta) {
0176       hadronicET = regionET;
0177       break;
0178     }
0179   }
0180
0181   hadronicET -= et;
0182
0183   double hoe = 0.0;
0184
0185   if (hadronicET > 0 && et > 0)
0186     hoe = (double)hadronicET / (double)et;
0187
0188   // set output
0189   return hoe;
0190 }