src/firmware/Stage1Layer2EtSumAlgorithmImpPP.cc

0001 ///
0002 /// \class l1t::Stage1Layer2EtSumAlgorithmImpPP
0003 ///
0004 /// \author: L. Apanasevich
0005 ///
0006 /// Description: first iteration of stage 1 jet sums algo
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
0018 l1t::Stage1Layer2EtSumAlgorithmImpPP::Stage1Layer2EtSumAlgorithmImpPP(CaloParamsHelper const* params)
0019     : params_(params) {
0020   //now do what ever initialization is needed
0021   for (unsigned int i = 0; i < L1CaloRegionDetId::N_PHI; i++) {
0022     sinPhi.push_back(sin(2. * 3.1415927 * i * 1.0 / L1CaloRegionDetId::N_PHI));
0023     cosPhi.push_back(cos(2. * 3.1415927 * i * 1.0 / L1CaloRegionDetId::N_PHI));
0024   }
0025 }
0026
0027 //double l1t::Stage1Layer2EtSumAlgorithmImpPP::regionPhysicalEt(const l1t::CaloRegion& cand) const {
0028 //  return jetLsb*cand.hwPt();
0029 //}
0030
0031 void l1t::Stage1Layer2EtSumAlgorithmImpPP::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   unsigned int sumET = 0;
0036   double sumEx = 0;
0037   double sumEy = 0;
0038   unsigned int sumHT = 0;
0039   double sumHx = 0;
0040   double sumHy = 0;
0041
0042   std::vector<l1t::CaloRegion> subRegions;
0043
0044   //Region Correction will return uncorrected subregions if
0045   //regionPUSType is set to None in the config
0046   double jetLsb = params_->jetLsb();
0047
0048   int etSumEtaMinEt = params_->etSumEtaMin(0);
0049   int etSumEtaMaxEt = params_->etSumEtaMax(0);
0050   //double etSumEtThresholdEt = params_->etSumEtThreshold(0);
0051   int etSumEtThresholdEt = (int)(params_->etSumEtThreshold(0) / jetLsb);
0052
0053   int etSumEtaMinHt = params_->etSumEtaMin(1);
0054   int etSumEtaMaxHt = params_->etSumEtaMax(1);
0055   //double etSumEtThresholdHt = params_->etSumEtThreshold(1);
0056   int etSumEtThresholdHt = (int)(params_->etSumEtThreshold(1) / jetLsb);
0057
0058   RegionCorrection(regions, &subRegions, params_);
0059
0060   double towerLsb = params_->towerLsbSum();
0061   int jetSeedThreshold = floor(params_->jetSeedThreshold() / towerLsb + 0.5);
0062   // ----- cluster jets for repurposing of MHT phi (use if for angle between leading jet)
0063   std::vector<l1t::Jet> unCorrJets;
0064   std::vector<l1t::Jet> unSortedJets;
0065   std::vector<l1t::Jet> SortedJets;
0066   slidingWindowJetFinder(jetSeedThreshold, &subRegions, &unCorrJets);
0067
0068   //if jetCalibrationType is set to None in the config
0069   std::string jetCalibrationType = params_->jetCalibrationType();
0070   std::vector<double> jetCalibrationParams = params_->jetCalibrationParams();
0071   JetCalibration(&unCorrJets, jetCalibrationParams, &unSortedJets, jetCalibrationType, towerLsb);
0072
0073   SortJets(&unSortedJets, &SortedJets);
0074   int dijet_phi = DiJetPhi(&SortedJets);
0075
0076   for (std::vector<CaloRegion>::const_iterator region = subRegions.begin(); region != subRegions.end(); region++) {
0077     if (region->hwEta() < etSumEtaMinEt || region->hwEta() > etSumEtaMaxEt) {
0078       continue;
0079     }
0080
0081     //double regionET= regionPhysicalEt(*region);
0082     int regionET = region->hwPt();
0083
0084     if (regionET >= etSumEtThresholdEt) {
0085       sumET += regionET;
0086       sumEx += (((double)regionET) * cosPhi[region->hwPhi()]);
0087       sumEy += (((double)regionET) * sinPhi[region->hwPhi()]);
0088     }
0089   }
0090
0091   for (std::vector<CaloRegion>::const_iterator region = subRegions.begin(); region != subRegions.end(); region++) {
0092     if (region->hwEta() < etSumEtaMinHt || region->hwEta() > etSumEtaMaxHt) {
0093       continue;
0094     }
0095
0096     //double regionET= regionPhysicalEt(*region);
0097     int regionET = region->hwPt();
0098
0099     if (regionET >= etSumEtThresholdHt) {
0100       sumHT += regionET;
0101       sumHx += (((double)regionET) * cosPhi[region->hwPhi()]);
0102       sumHy += (((double)regionET) * sinPhi[region->hwPhi()]);
0103     }
0104   }
0105
0106   unsigned int MET = ((unsigned int)sqrt(sumEx * sumEx + sumEy * sumEy));
0107   unsigned int MHT = ((unsigned int)sqrt(sumHx * sumHx + sumHy * sumHy));
0108
0109   double physicalPhi = atan2(sumEy, sumEx) + 3.1415927;
0110   // Global Trigger expects MET phi to be 0-71 (e.g. tower granularity)
0111   // Although we calculated it with regions, there is some benefit to interpolation.
0112   unsigned int iPhiET = 4 * L1CaloRegionDetId::N_PHI * physicalPhi / (2 * 3.1415927);
0113
0114   const ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > etLorentz(0, 0, 0, 0);
0115
0116   // Set quality (i.e. overflow) bits appropriately
0117   int METqual = 0;
0118   int MHTqual = 0;
0119   int ETTqual = 0;
0120   int HTTqual = 0;
0121   if (MET >= 0xfff)  // MET 12 bits
0122     METqual = 1;
0123   if (MHT >= 0x7f)  // MHT 7 bits
0124     MHTqual = 1;
0125   if (sumET >= 0xfff)
0126     ETTqual = 1;
0127   if (sumHT >= 0xfff)
0128     HTTqual = 1;
0129
0130   // scale MHT by sumHT
0131   // int mtmp = floor (((double) MHT / (double) sumHT)*100 + 0.5);
0132   // double mtmp = ((double) MHT / (double) sumHT);
0133   // int rank=params_->HtMissScale().rank(mtmp);
0134   // MHT=rank;
0135
0136   uint16_t MHToHT = MHToverHT(MHT, sumHT);
0137   unsigned int iPhiHT = (unsigned int)dijet_phi;
0138
0139   l1t::EtSum etMiss(*&etLorentz, EtSum::EtSumType::kMissingEt, MET, 0, iPhiET, METqual);
0140   l1t::EtSum htMiss(*&etLorentz, EtSum::EtSumType::kMissingHt, MHToHT & 0x7f, 0, iPhiHT, MHTqual);
0141   l1t::EtSum etTot(*&etLorentz, EtSum::EtSumType::kTotalEt, sumET & 0xfff, 0, 0, ETTqual);
0142   l1t::EtSum htTot(*&etLorentz, EtSum::EtSumType::kTotalHt, sumHT & 0xfff, 0, 0, HTTqual);
0143
0144   std::vector<l1t::EtSum> preGtEtSums = {etMiss, htMiss, etTot, htTot};
0145
0146   EtSumToGtScales(params_, &preGtEtSums, etsums);
0147 }
0148
0149 int l1t::Stage1Layer2EtSumAlgorithmImpPP::DiJetPhi(const std::vector<l1t::Jet>* jets) const {
0150   int dphi = 10;  // initialize to negative physical dphi value
0151   if (jets->size() < 2)
0152     return dphi;  // size() not really reliable as we pad the size to 8 (4cen+4for) in the sorter
0153   if ((*jets).at(0).hwPt() == 0)
0154     return dphi;
0155   if ((*jets).at(1).hwPt() == 0)
0156     return dphi;
0157
0158   int iphi1 = (*jets).at(0).hwPhi();
0159   int iphi2 = (*jets).at(1).hwPhi();
0160
0161   int difference = abs(iphi1 - iphi2);
0162
0163   if (difference > 9)
0164     difference = L1CaloRegionDetId::N_PHI - difference;  // make Physical dphi always positive
0165   return difference;
0166 }
0167
0168 uint16_t l1t::Stage1Layer2EtSumAlgorithmImpPP::MHToverHT(uint16_t num, uint16_t den) const {
0169   uint16_t result;
0170   uint32_t numerator(num), denominator(den);
0171   assert(den != 0);
0172
0173   if (numerator == denominator)
0174     result = 0x7f;
0175   else {
0176     numerator = numerator << 7;
0177     result = numerator / denominator;
0178     result = result & 0x7f;
0179   }
0180   return result;
0181 }