Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 //
0002 /// \class l1t::Stage2Layer2JetAlgorithmFirmwareImp1
0003 ///
0004 /// \author: Adam Elwood and Matthew Citron
0005 ///
0006 /// Description: Implementation of Jad's asymmetric map overlap algorithm with donut subtraction
0007 
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "L1Trigger/L1TCalorimeter/interface/Stage2Layer2JetAlgorithmFirmware.h"
0010 #include "DataFormats/Math/interface/LorentzVector.h"
0011 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0012 #include "L1Trigger/L1TCalorimeter/interface/AccumulatingSort.h"
0013 #include "L1Trigger/L1TCalorimeter/interface/BitonicSort.h"
0014 #include "CondFormats/L1TObjects/interface/CaloParams.h"
0015 #include "TMath.h"
0016 
0017 #include <vector>
0018 #include <algorithm>
0019 #include <cmath>
0020 
0021 namespace l1t {
0022   bool operator>(const l1t::Jet& a, const l1t::Jet& b) { return a.hwPt() > b.hwPt(); }
0023 }  // namespace l1t
0024 
0025 // jet mask, needs to be configurable at some point
0026 // just a square for now
0027 // for 1 do greater than, for 2 do greater than equal to
0028 
0029 namespace {
0030   constexpr int mask_[9][9] = {
0031       {1, 2, 2, 2, 2, 2, 2, 2, 2},
0032       {1, 1, 2, 2, 2, 2, 2, 2, 2},
0033       {1, 1, 1, 2, 2, 2, 2, 2, 2},
0034       {1, 1, 1, 1, 2, 2, 2, 2, 2},
0035       {1, 1, 1, 1, 0, 2, 2, 2, 2},
0036       {1, 1, 1, 1, 1, 2, 2, 2, 2},
0037       {1, 1, 1, 1, 1, 1, 2, 2, 2},
0038       {1, 1, 1, 1, 1, 1, 1, 2, 2},
0039       {1, 1, 1, 1, 1, 1, 1, 1, 2},
0040   };
0041 
0042 }
0043 
0044 l1t::Stage2Layer2JetAlgorithmFirmwareImp1::Stage2Layer2JetAlgorithmFirmwareImp1(CaloParamsHelper const* params)
0045     : params_(params) {}
0046 
0047 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::processEvent(const std::vector<l1t::CaloTower>& towers,
0048                                                              std::vector<l1t::Jet>& jets,
0049                                                              std::vector<l1t::Jet>& alljets) {
0050   // find jets
0051   create(towers, jets, alljets, params_->jetPUSType());
0052 
0053   // calibrate all jets
0054   calibrate(alljets, 0, true);  // pass all jets and the hw threshold above which to calibrate
0055 
0056   // jets accumulated sort
0057   accuSort(jets);
0058 }
0059 
0060 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::create(const std::vector<l1t::CaloTower>& towers,
0061                                                        std::vector<l1t::Jet>& jets,
0062                                                        std::vector<l1t::Jet>& alljets,
0063                                                        std::string PUSubMethod) {
0064   // etaSide=1 is positive eta, etaSide=-1 is negative eta
0065   for (int etaSide = 1; etaSide >= -1; etaSide -= 2) {
0066     // the 4 groups of rings
0067     std::vector<int> ringGroup1, ringGroup2, ringGroup3, ringGroup4;
0068     for (int i = 1; i <= CaloTools::mpEta(CaloTools::kHFEnd); i++) {
0069       if (!((i - 1) % 4))
0070         ringGroup1.push_back(i * etaSide);
0071       else if (!((i - 2) % 4))
0072         ringGroup2.push_back(i * etaSide);
0073       else if (!((i - 3) % 4))
0074         ringGroup3.push_back(i * etaSide);
0075       else if (!((i - 4) % 4))
0076         ringGroup4.push_back(i * etaSide);
0077     }
0078     std::vector<std::vector<int> > theRings = {ringGroup1, ringGroup2, ringGroup3, ringGroup4};
0079 
0080     // the 24 jets in this eta side
0081     std::vector<l1t::Jet> jetsHalf;
0082 
0083     // loop over the 4 groups of rings
0084     for (unsigned ringGroupIt = 1; ringGroupIt <= theRings.size(); ringGroupIt++) {
0085       // the 6 accumulated jets
0086       std::vector<l1t::Jet> jetsAccu;
0087 
0088       // loop over the 10 rings in this group
0089       for (unsigned ringIt = 0; ringIt < theRings.at(ringGroupIt - 1).size(); ringIt++) {
0090         int ieta = theRings.at(ringGroupIt - 1).at(ringIt);
0091 
0092         // the jets in this ring
0093         std::vector<l1t::Jet> jetsRing;
0094 
0095         // loop over phi in the ring
0096         for (int iphi = 1; iphi <= CaloTools::kHBHENrPhi; ++iphi) {
0097           // no more than 18 jets per ring
0098           if (jetsRing.size() == 18)
0099             break;
0100 
0101           // seed tower
0102           const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(ieta), iphi);
0103 
0104           int seedEt = tow.hwPt();
0105           int iDelay = (tow.hwQual() & 0b0100) >> 2;
0106           int iEt = seedEt;
0107           bool satSeed = false;
0108           bool vetoCandidate = false;
0109 
0110           // check it passes the seed threshold
0111           if (iEt < floor(params_->jetSeedThreshold() / params_->towerLsbSum()))
0112             continue;
0113 
0114           if (seedEt == CaloTools::kSatHcal || seedEt == CaloTools::kSatEcal || seedEt == CaloTools::kSatTower)
0115             satSeed = true;
0116 
0117           // loop over towers in this jet
0118           for (int deta = -4; deta < 5; ++deta) {
0119             for (int dphi = -4; dphi < 5; ++dphi) {
0120               int towEt = 0;
0121               int towDelay = 0;
0122               int ietaTest = ieta + deta;
0123               int iphiTest = iphi + dphi;
0124 
0125               // wrap around phi
0126               while (iphiTest > CaloTools::kHBHENrPhi)
0127                 iphiTest -= CaloTools::kHBHENrPhi;
0128               while (iphiTest < 1)
0129                 iphiTest += CaloTools::kHBHENrPhi;
0130 
0131               // wrap over eta=0
0132               if (ieta > 0 && ietaTest <= 0)
0133                 ietaTest -= 1;
0134               if (ieta < 0 && ietaTest >= 0)
0135                 ietaTest += 1;
0136 
0137               // check jet mask and sum tower et
0138               const CaloTower& towTest = CaloTools::getTower(towers, CaloTools::caloEta(ietaTest), iphiTest);
0139               towEt = towTest.hwPt();
0140               towDelay = (towTest.hwQual() & 0b0100) >> 2;
0141 
0142               if (mask_[8 - (dphi + 4)][deta + 4] == 0)
0143                 continue;
0144               else if (mask_[8 - (dphi + 4)][deta + 4] == 1)
0145                 vetoCandidate = (seedEt < towEt);
0146               else if (mask_[8 - (dphi + 4)][deta + 4] == 2)
0147                 vetoCandidate = (seedEt <= towEt);
0148 
0149               if (vetoCandidate)
0150                 break;
0151               else {
0152                 iEt += towEt;
0153                 if (abs(ieta) < 29 && abs(ietaTest) < 29)
0154                   iDelay += towDelay;  // don't include HF feature bits in HBHE flagged jets
0155               }
0156             }
0157             if (vetoCandidate)
0158               break;
0159           }
0160 
0161           // add the jet to the list
0162           if (!vetoCandidate) {
0163             int rawEt = iEt;
0164             int puEt(0);
0165 
0166             math::XYZTLorentzVector p4;
0167             int caloEta = CaloTools::caloEta(ieta);
0168             l1t::Jet jet(p4, -999, caloEta, iphi, 0);
0169 
0170             if (!params_->jetBypassPUS() && !satSeed) {
0171               // comment out for emulator studies
0172               if (params_->jetPUSUsePhiRing())
0173                 PUSubMethod = "PhiRing1";
0174               else
0175                 PUSubMethod = "ChunkyDonut";
0176 
0177               if (PUSubMethod == "Donut") {
0178                 puEt = donutPUEstimate(ieta, iphi, 5, towers);
0179                 iEt -= puEt;
0180               }
0181 
0182               if (PUSubMethod == "ChunkyDonut") {
0183                 puEt = chunkyDonutPUEstimate(jet, 5, towers);
0184                 iEt -= puEt;
0185               }
0186 
0187               if (PUSubMethod.find("ChunkySandwich") != std::string::npos) {
0188                 puEt = chunkySandwichPUEstimate(jet, 5, towers, PUSubMethod);
0189                 iEt -= puEt;
0190               }
0191 
0192               if (PUSubMethod == "PhiRing1") {
0193                 int size = 5;
0194                 std::map<int, int> SumEtEtaMap = getSumEtEtaMap(towers);
0195                 // int jetPhi=jet.hwPhi();
0196                 int jetEta = CaloTools::mpEta(jet.hwEta());
0197                 int jetEtaLow = jetEta - size + 1;
0198                 int jetEtaHigh = jetEta + size - 1;
0199                 if (jetEta > 0 && jetEtaLow <= 0)
0200                   jetEtaLow -= 1;
0201                 if (jetEta < 0 && jetEtaHigh >= 0)
0202                   jetEtaHigh += 1;
0203                 puEt = 0;
0204                 for (int ieta = jetEtaLow; ieta <= jetEtaHigh; ++ieta) {
0205                   if (ieta == 0)
0206                     continue;  // eta start from +- 1
0207                   auto iter = SumEtEtaMap.find(ieta);
0208                   puEt += iter->second;
0209                 }
0210                 puEt = (puEt * (size * 2 - 1)) / CaloTools::kHBHENrPhi;
0211                 iEt -= puEt;
0212               }
0213 
0214               if (PUSubMethod == "PhiRing2") {
0215                 int size = 5;
0216                 std::map<int, int> SumEtEtaMap = getSumEtEtaMap(towers);
0217                 // int jetPhi=jet.hwPhi();
0218                 int jetEta = CaloTools::mpEta(jet.hwEta());
0219                 int jetEtaLow = jetEta - size + 1;
0220                 int jetEtaHigh = jetEta + size - 1;
0221 
0222                 if (jetEta > 0 && jetEtaLow <= 0)
0223                   jetEtaLow -= 1;
0224                 if (jetEta < 0 && jetEtaHigh >= 0)
0225                   jetEtaHigh += 1;
0226                 puEt = 0;
0227 
0228                 for (int ieta = jetEtaLow; ieta <= jetEtaHigh; ++ieta) {
0229                   if (ieta == 0)
0230                     continue;
0231                   auto iter = SumEtEtaMap.find(ieta);
0232                   puEt += iter->second;
0233                 }
0234                 puEt -= iEt;
0235                 puEt = (puEt * (size * 2 - 1)) / (CaloTools::kHBHENrPhi - (size * 2 - 1));
0236                 iEt -= puEt;
0237               }
0238             }
0239 
0240             if (iEt <= 0)
0241               continue;
0242 
0243             // if tower Et is saturated, saturate jet Et
0244             if (seedEt == CaloTools::kSatHcal || seedEt == CaloTools::kSatEcal || seedEt == CaloTools::kSatTower)
0245               iEt = CaloTools::kSatJet;
0246 
0247             jet.setHwPt(iEt);
0248             if (iDelay >= 2)
0249               jet.setHwQual(1);
0250             jet.setRawEt((short int)rawEt);
0251             jet.setSeedEt((short int)seedEt);
0252             jet.setTowerIEta((short int)caloEta);
0253             jet.setTowerIPhi((short int)iphi);
0254             jet.setPUEt((short int)puEt);
0255 
0256             jetsRing.push_back(jet);
0257             alljets.push_back(jet);
0258           }
0259         }
0260 
0261         // jet energy corrections
0262         calibrate(jetsRing, 0, false);  // pass the jet collection and the hw threshold above which to calibrate
0263 
0264         // sort these jets and keep top 6
0265         auto start_ = jetsRing.begin();
0266         auto end_ = jetsRing.end();
0267         BitonicSort<l1t::Jet>(down, start_, end_);
0268         if (jetsRing.size() > 6)
0269           jetsRing.resize(6);
0270 
0271         // update jets
0272         jets.insert(jets.end(), jetsRing.begin(), jetsRing.end());
0273       }
0274     }
0275   }
0276 }
0277 
0278 //Accumulating sort
0279 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::accuSort(std::vector<l1t::Jet>& jets) {
0280   math::PtEtaPhiMLorentzVector emptyP4;
0281   l1t::Jet tempJet(emptyP4, 0, 0, 0, 0);
0282   std::vector<std::vector<l1t::Jet> > jetEtaPos(41, std::vector<l1t::Jet>(18, tempJet));
0283   std::vector<std::vector<l1t::Jet> > jetEtaNeg(41, std::vector<l1t::Jet>(18, tempJet));
0284 
0285   for (unsigned int iJet = 0; iJet < jets.size(); iJet++) {
0286     if (jets.at(iJet).hwEta() > 0)
0287       jetEtaPos.at(jets.at(iJet).hwEta() - 1).at((72 - jets.at(iJet).hwPhi()) / 4) = jets.at(iJet);
0288     else
0289       jetEtaNeg.at(-(jets.at(iJet).hwEta() + 1)).at((72 - jets.at(iJet).hwPhi()) / 4) = jets.at(iJet);
0290   }
0291 
0292   AccumulatingSort<l1t::Jet> etaPosSorter(7);
0293   AccumulatingSort<l1t::Jet> etaNegSorter(7);
0294   std::vector<l1t::Jet> accumEtaPos;
0295   std::vector<l1t::Jet> accumEtaNeg;
0296 
0297   for (int ieta = 0; ieta < 41; ++ieta) {
0298     // eta +
0299     std::vector<l1t::Jet>::iterator start_, end_;
0300     start_ = jetEtaPos.at(ieta).begin();
0301     end_ = jetEtaPos.at(ieta).end();
0302     BitonicSort<l1t::Jet>(down, start_, end_);
0303     etaPosSorter.Merge(jetEtaPos.at(ieta), accumEtaPos);
0304 
0305     // eta -
0306     start_ = jetEtaNeg.at(ieta).begin();
0307     end_ = jetEtaNeg.at(ieta).end();
0308     BitonicSort<l1t::Jet>(down, start_, end_);
0309     etaNegSorter.Merge(jetEtaNeg.at(ieta), accumEtaNeg);
0310   }
0311 
0312   //check for 6 & 7th jets with same et and eta. Keep jet with larger phi
0313   if (accumEtaPos.at(6).hwPt() == accumEtaPos.at(5).hwPt() && accumEtaPos.at(6).hwEta() == accumEtaPos.at(5).hwEta() &&
0314       accumEtaPos.at(6).hwPhi() > accumEtaPos.at(5).hwPhi()) {
0315     accumEtaPos.at(5) = accumEtaPos.at(6);
0316   }
0317   if (accumEtaNeg.at(6).hwPt() == accumEtaNeg.at(5).hwPt() && accumEtaNeg.at(6).hwEta() == accumEtaNeg.at(5).hwEta() &&
0318       accumEtaNeg.at(6).hwPhi() > accumEtaNeg.at(5).hwPhi()) {
0319     accumEtaNeg.at(5) = accumEtaNeg.at(6);
0320   }
0321 
0322   //truncate
0323   accumEtaPos.resize(6);
0324   accumEtaNeg.resize(6);
0325 
0326   // put all 12 candidates in the original jet vector, removing zero energy ones
0327   jets.clear();
0328   for (const l1t::Jet& accjet : accumEtaPos) {
0329     if (accjet.hwPt() > 0)
0330       jets.push_back(accjet);
0331   }
0332   for (const l1t::Jet& accjet : accumEtaNeg) {
0333     if (accjet.hwPt() > 0)
0334       jets.push_back(accjet);
0335   }
0336 }
0337 
0338 //A function to return the value for donut subtraction around an ieta and iphi position for donut subtraction
0339 //Also pass it a vector to store the individual values of the strip for later testing
0340 //The size is the number of ieta/iphi units out the ring is (ie for 9x9 jets, we want the 11x11 for PUS therefore we want to go 5 out, so size is 5)
0341 int l1t::Stage2Layer2JetAlgorithmFirmwareImp1::donutPUEstimate(int jetEta,
0342                                                                int jetPhi,
0343                                                                int size,
0344                                                                const std::vector<l1t::CaloTower>& towers) {
0345   //ring is a vector with 4 ring strips, one for each side of the ring
0346   std::vector<int> ring(4, 0);
0347 
0348   int iphiUp = jetPhi + size;
0349   while (iphiUp > CaloTools::kHBHENrPhi)
0350     iphiUp -= CaloTools::kHBHENrPhi;
0351   int iphiDown = jetPhi - size;
0352   while (iphiDown < 1)
0353     iphiDown += CaloTools::kHBHENrPhi;
0354 
0355   int ietaUp = jetEta + size;    //(jetEta + size > CaloTools::mpEta(CaloTools::kHFEnd)) ? 999 : jetEta+size;
0356   int ietaDown = jetEta - size;  //(abs(jetEta - size) > CaloTools::mpEta(CaloTools::kHFEnd)) ? 999 : jetEta-size;
0357 
0358   for (int ieta = jetEta - size + 1; ieta < jetEta + size; ++ieta) {
0359     if (abs(ieta) > CaloTools::mpEta(CaloTools::kHFEnd) || abs(ieta) < 1)
0360       continue;
0361     int towerEta;
0362 
0363     if (jetEta > 0 && ieta <= 0) {
0364       towerEta = ieta - 1;
0365     } else if (jetEta < 0 && ieta >= 0) {
0366       towerEta = ieta + 1;
0367     } else {
0368       towerEta = ieta;
0369     }
0370 
0371     const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(towerEta), iphiUp);
0372     int towEt = tow.hwPt();
0373     ring[0] += towEt;
0374 
0375     const CaloTower& tow2 = CaloTools::getTower(towers, CaloTools::caloEta(towerEta), iphiDown);
0376     towEt = tow2.hwPt();
0377     ring[1] += towEt;
0378   }
0379 
0380   for (int iphi = jetPhi - size + 1; iphi < jetPhi + size; ++iphi) {
0381     int towerPhi = iphi;
0382     while (towerPhi > CaloTools::kHBHENrPhi)
0383       towerPhi -= CaloTools::kHBHENrPhi;
0384     while (towerPhi < 1)
0385       towerPhi += CaloTools::kHBHENrPhi;
0386 
0387     const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(ietaUp), towerPhi);
0388     int towEt = tow.hwPt();
0389     ring[2] += towEt;
0390 
0391     const CaloTower& tow2 = CaloTools::getTower(towers, CaloTools::caloEta(ietaDown), towerPhi);
0392     towEt = tow2.hwPt();
0393     ring[3] += towEt;
0394   }
0395 
0396   //for the Donut Subtraction we only use the middle 2 (in energy) ring strips
0397   std::sort(ring.begin(), ring.end(), std::greater<int>());
0398 
0399   return 4 * (ring[1] + ring[2]);  // This should really be multiplied by 4.5 not 4.
0400 }
0401 
0402 std::map<int, int> l1t::Stage2Layer2JetAlgorithmFirmwareImp1::getSumEtEtaMap(const std::vector<l1t::CaloTower>& towers) {
0403   std::map<int, int> SumEtEtaMap;
0404   int EtaMin = -41;
0405   int EtaMax = 41;
0406   int phiMin = 1;
0407   int phiMax = 72;
0408 
0409   for (int iEta = EtaMin; iEta <= EtaMax; iEta++) {
0410     int SumEt = 0;
0411     for (int iPhi = phiMin; iPhi <= phiMax; iPhi++) {
0412       const CaloTower& tow = CaloTools::getTower(towers, CaloTools::caloEta(iEta), iPhi);
0413       int towEt = tow.hwPt();
0414       SumEt += towEt;
0415 
0416     }  // end for iPih<=phiMax
0417     SumEtEtaMap[iEta] = SumEt;
0418 
0419   }  // end for iEta<=EtaMax
0420 
0421   return SumEtEtaMap;
0422 }
0423 
0424 std::vector<int> l1t::Stage2Layer2JetAlgorithmFirmwareImp1::getChunkyRing(l1t::Jet& jet,
0425                                                                           int size,
0426                                                                           const std::vector<l1t::CaloTower>& towers,
0427                                                                           const std::string chunkyString) {
0428   int jetPhi = jet.hwPhi();
0429   int jetEta = CaloTools::mpEta(jet.hwEta());
0430 
0431   // ring is a vector with 4 or 8 eta or strips, one for each side of the ring in ChunkyDonut, else pure phi
0432   int ringSize = 4;
0433   if (chunkyString == "ChunkySandwich8")
0434     ringSize = 8;
0435 
0436   std::vector<int> ring(ringSize, 0);
0437 
0438   // number of strips in donut - should make this configurable
0439   int nStrips = 3;
0440 
0441   // loop over strips
0442   for (int stripIt = 0; stripIt < nStrips; stripIt++) {
0443     int iphiUp = jetPhi + size + stripIt;
0444     int iphiDown = jetPhi - size - stripIt;
0445     while (iphiUp > CaloTools::kHBHENrPhi)
0446       iphiUp -= CaloTools::kHBHENrPhi;
0447     while (iphiDown < 1)
0448       iphiDown += CaloTools::kHBHENrPhi;
0449 
0450     // we will always do phiUp and PhiDown
0451     // do PhiUp and PhiDown
0452     for (int ieta = jetEta - size + 1; ieta < jetEta + size; ++ieta) {
0453       if (abs(ieta) > CaloTools::mpEta(CaloTools::kHFEnd))
0454         continue;
0455 
0456       int towEta = ieta;
0457       if (jetEta > 0 && towEta <= 0)
0458         towEta -= 1;
0459       if (jetEta < 0 && towEta >= 0)
0460         towEta += 1;
0461 
0462       const CaloTower& towPhiUp = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiUp);
0463       int towEt = towPhiUp.hwPt();
0464       ring[0] += towEt;
0465 
0466       const CaloTower& towPhiDown = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiDown);
0467       towEt = towPhiDown.hwPt();
0468       ring[1] += towEt;
0469     }
0470 
0471     //Only consider eta for chunkydonut
0472     if (chunkyString == "ChunkyDonut") {
0473       int ietaUp = jetEta + size + stripIt;
0474       int ietaDown = jetEta - size - stripIt;
0475       if (jetEta < 0 && ietaUp >= 0)
0476         ietaUp += 1;
0477       if (jetEta > 0 && ietaDown <= 0)
0478         ietaDown -= 1;
0479 
0480       // do EtaUp
0481       for (int iphi = jetPhi - size + 1; iphi < jetPhi + size; ++iphi) {
0482         if (abs(ietaUp) <= CaloTools::mpEta(CaloTools::kHFEnd)) {
0483           int towPhi = iphi;
0484           while (towPhi > CaloTools::kHBHENrPhi)
0485             towPhi -= CaloTools::kHBHENrPhi;
0486           while (towPhi < 1)
0487             towPhi += CaloTools::kHBHENrPhi;
0488 
0489           const CaloTower& towEtaUp = CaloTools::getTower(towers, CaloTools::caloEta(ietaUp), towPhi);
0490           int towEt = towEtaUp.hwPt();
0491           ring[2] += towEt;
0492         }
0493       }
0494 
0495       // do EtaDown
0496       for (int iphi = jetPhi - size + 1; iphi < jetPhi + size; ++iphi) {
0497         if (abs(ietaDown) <= CaloTools::mpEta(CaloTools::kHFEnd)) {
0498           int towPhi = iphi;
0499           while (towPhi > CaloTools::kHBHENrPhi)
0500             towPhi -= CaloTools::kHBHENrPhi;
0501           while (towPhi < 1)
0502             towPhi += CaloTools::kHBHENrPhi;
0503 
0504           const CaloTower& towEtaDown = CaloTools::getTower(towers, CaloTools::caloEta(ietaDown), towPhi);
0505           int towEt = towEtaDown.hwPt();
0506           ring[3] += towEt;
0507         }
0508       }
0509     } else if (chunkyString == "ChunkySandwich2") {  //simplest of the sandwiches - we will just reuse phi flaps
0510       for (int i = 0; i < 2; ++i) {
0511         ring[i + 2] = ring[i];
0512       }
0513     } else if (chunkyString == "ChunkySandwich2Ave" ||
0514                chunkyString == "ChunkySandwich") {  //simplest of the sandwiches - we will just reuse phi flaps
0515       ring[2] = ((ring[0] + ring[1]) >> 1);         //bitwise, effective div by 2
0516       ring[3] = (ring[0] + ring[1]);                // just make sure it is >= max val
0517     } else if (chunkyString == "ChunkySandwich4" ||
0518                chunkyString ==
0519                    "ChunkySandwich8") {  //Need to calculate two or siz additional phiflaps. Note that ChunkySandwich defaults to ChunkySandwich4 variant
0520 
0521       for (int rI = 0; rI < (ringSize - 2) / 2; ++rI) {
0522         int iphiUp2 = jetPhi + size + (stripIt + nStrips * (rI + 1));
0523         int iphiDown2 = jetPhi - size - (stripIt + nStrips * (rI + 1));
0524         while (iphiUp2 > CaloTools::kHBHENrPhi)
0525           iphiUp2 -= CaloTools::kHBHENrPhi;
0526         while (iphiDown2 < 1)
0527           iphiDown2 += CaloTools::kHBHENrPhi;
0528 
0529         for (int ieta = jetEta - size + 1; ieta < jetEta + size; ++ieta) {
0530           if (abs(ieta) > CaloTools::mpEta(CaloTools::kHFEnd))
0531             continue;
0532 
0533           int towEta = ieta;
0534           if (jetEta > 0 && towEta <= 0)
0535             towEta -= 1;
0536           if (jetEta < 0 && towEta >= 0)
0537             towEta += 1;
0538 
0539           const CaloTower& towPhiUp = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiUp2);
0540           int towEt = towPhiUp.hwPt();
0541           ring[2 + rI * 2] += towEt;
0542 
0543           const CaloTower& towPhiDown = CaloTools::getTower(towers, CaloTools::caloEta(towEta), iphiDown2);
0544           towEt = towPhiDown.hwPt();
0545           ring[3 + rI * 2] += towEt;
0546         }
0547       }
0548     }
0549   }
0550 
0551   //sort our final ring and return;
0552   std::sort(ring.begin(), ring.end());
0553 
0554   return ring;
0555 }
0556 
0557 int l1t::Stage2Layer2JetAlgorithmFirmwareImp1::chunkyDonutPUEstimate(l1t::Jet& jet,
0558                                                                      int size,
0559                                                                      const std::vector<l1t::CaloTower>& towers) {
0560   // ring is a vector with 4 ring strips, one for each side of the ring
0561   // PhiUp, PhiDown, EtaUp, EtaDown, sorted w/ call to chunkyring
0562   std::vector<int> ring = getChunkyRing(jet, size, towers, "ChunkyDonut");
0563   for (unsigned int i = 0; i < 4; ++i)
0564     jet.setPUDonutEt(i, (short int)ring[i]);
0565   return (ring[0] + ring[1] + ring[2]);
0566 }
0567 
0568 //Following example of chunky donut, but considering only phi flaps
0569 //Useful for simple handling of zeroed ieta strips and for high PU environments like PbPb
0570 int l1t::Stage2Layer2JetAlgorithmFirmwareImp1::chunkySandwichPUEstimate(l1t::Jet& jet,
0571                                                                         int size,
0572                                                                         const std::vector<l1t::CaloTower>& towers,
0573                                                                         const std::string chunkySandwichStr) {
0574   //ring will be handled identically by all variants for now - could consider a different picking w/ ChunkySandwich8, say to exclude a downward fluctuation w/ 1,2,3 from ring, etc.
0575   std::vector<int> ring = getChunkyRing(jet, size, towers, chunkySandwichStr);
0576   for (unsigned int i = 0; i < 4; ++i)
0577     jet.setPUDonutEt(i, (short int)ring[i]);
0578   return (ring[0] + ring[1] + ring[2]);
0579 }
0580 
0581 void l1t::Stage2Layer2JetAlgorithmFirmwareImp1::calibrate(std::vector<l1t::Jet>& jets,
0582                                                           int calibThreshold,
0583                                                           bool isAllJets) {
0584   if (params_->jetCalibrationType() == "function6PtParams22EtaBins") {  //One eta bin per region
0585 
0586     //Check the vector of calibration parameters is the correct size
0587     //Order the vector in terms of the parameters per eta bin, starting in -ve eta
0588     //So first 6 entries are very negative eta, next 6 are the next bin etc.
0589 
0590     if (params_->jetCalibrationParams().size() != 6 * 22) {
0591       edm::LogError("l1t|stage 2") << "Invalid input vector to calo params. Input vector of size: "
0592                                    << params_->jetCalibrationParams().size()
0593                                    << "  Require size: 132  Not calibrating Stage 2 Jets" << std::endl;
0594       return;
0595     }
0596 
0597     //Loop over jets and apply corrections
0598     for (std::vector<l1t::Jet>::iterator jet = jets.begin(); jet != jets.end(); jet++) {
0599       //Check jet is above the calibration threshold, if not do nothing
0600       if (jet->hwPt() < calibThreshold)
0601         continue;
0602       if (jet->hwPt() >= 0xFFFF)
0603         continue;
0604 
0605       int etaBin = CaloTools::regionEta(jet->hwEta());
0606 
0607       //Get the parameters from the vector
0608       //Each 6 values are the parameters for an eta bin
0609       double params[6];
0610       for (int i = 0; i < 6; i++) {
0611         params[i] = params_->jetCalibrationParams()[etaBin * 6 + i];
0612       }
0613 
0614       //Perform the correction based on the calibration function defined
0615       //in calibFit
0616       //This is derived from the actual physical pt of the jets, not the hwEt
0617       //This needs to be addressed in the future
0618       double ptPhys = jet->hwPt() * params_->jetLsb();
0619       double correction = calibFit(ptPhys, params);
0620 
0621       jet->setHwPt(correction * jet->hwPt());
0622     }
0623 
0624   } else if (params_->jetCalibrationType() == "function8PtParams22EtaBins") {
0625     // as above but with cap on max correction at low pT
0626 
0627     if (params_->jetCalibrationParams().size() != 8 * 22) {
0628       edm::LogError("l1t|stage 2") << "Invalid input vector to calo params. Input vector of size: "
0629                                    << params_->jetCalibrationParams().size()
0630                                    << "  Require size: 176  Not calibrating Stage 2 Jets" << std::endl;
0631       return;
0632     }
0633 
0634     for (std::vector<l1t::Jet>::iterator jet = jets.begin(); jet != jets.end(); jet++) {
0635       if (jet->hwPt() < calibThreshold)
0636         continue;
0637       if (jet->hwPt() >= 0xFFFF)
0638         continue;
0639 
0640       int etaBin = CaloTools::regionEta(jet->hwEta());
0641 
0642       double params[8];
0643       for (int i = 0; i < 8; i++) {
0644         params[i] = params_->jetCalibrationParams()[etaBin * 8 + i];
0645       }
0646 
0647       double ptPhys = jet->hwPt() * params_->jetLsb();
0648       double correction = params[6];
0649       if (ptPhys > params[7])
0650         correction = calibFit(ptPhys, params);
0651 
0652       jet->setHwPt(correction * jet->hwPt());
0653     }
0654 
0655   } else if (params_->jetCalibrationType() == "functionErf11PtParams16EtaBins") {
0656     // as above but with cap on max correction at low pT
0657 
0658     if (params_->jetCalibrationParams().size() != 11 * 16) {
0659       edm::LogError("l1t|stage 2") << "Invalid input vector to calo params. Input vector of size: "
0660                                    << params_->jetCalibrationParams().size()
0661                                    << "  Require size: 176  Not calibrating Stage 2 Jets" << std::endl;
0662       return;
0663     }
0664 
0665     for (std::vector<l1t::Jet>::iterator jet = jets.begin(); jet != jets.end(); jet++) {
0666       if (jet->hwPt() < calibThreshold)
0667         continue;
0668       if (jet->hwPt() >= 0xFFFF)
0669         continue;
0670 
0671       int etaBin = CaloTools::bin16Eta(jet->hwEta());
0672 
0673       double params[11];
0674       for (int i = 0; i < 11; i++) {
0675         params[i] = params_->jetCalibrationParams()[etaBin * 11 + i];
0676       }
0677 
0678       double ptPhys = jet->hwPt() * params_->jetLsb();
0679       double correction;
0680 
0681       if (ptPhys < params[8])
0682         correction = params[7];
0683       else if (ptPhys > params[10])
0684         correction = params[9];
0685       else
0686         correction = calibFitErr(ptPhys, params);
0687 
0688       jet->setHwPt(correction * jet->hwPt());
0689     }
0690 
0691   } else if (params_->jetCalibrationType() == "LUT") {
0692     // Calibrate using 3 LUTs: pt and eta compression LUTs, and a multiplicand/addend LUT.
0693     // The pt and eta are each converted to a compressed scale using individual LUTs
0694     // pt : 8 -> 4 bits, eta 6 -> 4 bits
0695     // **updated** pt : 8 -> 6 bits, eta 6 -> 6 bits (uncompressed)
0696     // This then forms an address. Using the third LUT, we get a
0697     // multiplicand & addend, so we can do y = m*x + c on the original
0698     // (i.e. non-compressed) jet pt.
0699     // The multiplicand is 10-bit unsigned, addend is 8-bit signed.
0700 
0701     //Loop over jets and apply corrections
0702     for (std::vector<l1t::Jet>::iterator jet = jets.begin(); jet != jets.end(); jet++) {
0703       //Check jet is above the calibration threshold, if not do nothing
0704       if (jet->hwPt() < calibThreshold)
0705         continue;
0706 
0707       //don't calibrate saturated jets for HT
0708       if (isAllJets && (jet->hwPt() == CaloTools::kSatJet))
0709         continue;
0710 
0711       // In the firmware, we take bits 1 to 8 of the hwPt.
0712       // To avoid getting nonsense by only taking smaller bits,
0713       // any values larger than 511 are automatically set to 511.
0714       int jetHwPt = jet->hwPt();
0715       if (jetHwPt >= 0x200) {
0716         jetHwPt = 0x1FF;
0717       }
0718 
0719       unsigned int ptCompNrBits = params_->jetCompressPtLUT()->nrBitsData();
0720       unsigned int ptBin = params_->jetCompressPtLUT()->data(jetHwPt >> 1);
0721       unsigned int etaBin = params_->jetCompressEtaLUT()->data(abs(CaloTools::mpEta(jet->hwEta())));
0722       unsigned int compBin = (etaBin << ptCompNrBits) | ptBin;
0723 
0724       unsigned int addPlusMult = params_->jetCalibrationLUT()->data(compBin);
0725       unsigned int multiplier = addPlusMult & 0x3ff;
0726       // handles -ve numbers correctly
0727       int8_t addend = (addPlusMult >> 10);
0728       unsigned int jetPtCorr = ((jet->hwPt() * multiplier) >> 9) + addend;
0729 
0730       if (jetPtCorr < 0xFFFF) {
0731         jet->setHwPt(jetPtCorr);
0732       } else {
0733         jet->setHwPt(0xFFFF);
0734       }
0735     }
0736 
0737   } else {
0738     if (params_->jetCalibrationType() != "None" && params_->jetCalibrationType() != "none")
0739       edm::LogError("l1t|stage 2") << "Invalid calibration type in calo params. Not calibrating Stage 2 Jets"
0740                                    << std::endl;
0741     return;
0742   }
0743 }
0744 
0745 //Function for the calibration, correct as a function of pT in bins of eta
0746 double l1t::Stage2Layer2JetAlgorithmFirmwareImp1::calibFit(double pt, double* par) {
0747   double logX = log10(pt);
0748 
0749   double term1 = par[1] / (logX * logX + par[2]);
0750   double term2 = par[3] * exp(-par[4] * ((logX - par[5]) * (logX - par[5])));
0751 
0752   // Final fitting function, with sanity check
0753   double f = par[0] + term1 + term2;
0754   if (f < 0)
0755     f = 0;
0756   if (f != f)  // stop NaN
0757     f = 1;
0758   return f;
0759 }
0760 
0761 //NEW Function for the calibration, correct as a function of pT in bins of eta
0762 double l1t::Stage2Layer2JetAlgorithmFirmwareImp1::calibFitErr(double pt, double* par) {
0763   double f = par[0] + par[1] * TMath::Erf(par[2] * (log10(pt) - par[3]) +
0764                                           par[4] * exp(par[5] * (log10(pt) - par[6]) * (log10(pt) - par[6])));
0765   // sanity check
0766   if (f < 0)
0767     f = 0;
0768   if (f != f)  // stop NaN
0769     f = 1;
0770   return f;
0771 }
[ 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!