Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-02-20 03:45:40

0001 // This function fetches Layer1 ECAL and HCAL LUTs from CMSSW configuration
0002 // It is provided as a global helper function outside of class structure
0003 // so that it can be shared by L1CaloLayer1 and L1CaloLayer1Spy
0004 
0005 #include <vector>
0006 
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include "L1Trigger/L1TCalorimeter/interface/CaloParamsHelper.h"
0012 #include "CondFormats/L1TObjects/interface/CaloParams.h"
0013 #include "CondFormats/DataRecord/interface/L1TCaloParamsRcd.h"
0014 
0015 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0016 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0017 #include "DataFormats/HcalDetId/interface/HcalTrigTowerDetId.h"
0018 #include "DataFormats/HcalDigi/interface/HcalTriggerPrimitiveSample.h"
0019 #include "DataFormats/HcalDigi/interface/HcalTriggerPrimitiveDigi.h"
0020 
0021 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0022 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0023 
0024 #include "L1TCaloLayer1FetchLUTs.hh"
0025 #include "UCTLogging.hh"
0026 
0027 using namespace l1tcalo;
0028 
0029 bool L1TCaloLayer1FetchLUTs(
0030     const L1TCaloLayer1FetchLUTsTokens &iTokens,
0031     const edm::EventSetup &iSetup,
0032     std::vector<std::array<std::array<std::array<uint32_t, nEtBins>, nCalSideBins>, nCalEtaBins> > &eLUT,
0033     std::vector<std::array<std::array<std::array<uint32_t, nEtBins>, nCalSideBins>, nCalEtaBins> > &hLUT,
0034     std::vector<std::array<std::array<uint32_t, nEtBins>, nHfEtaBins> > &hfLUT,
0035     std::vector<unsigned long long int> &hcalFBLUT,
0036     std::vector<unsigned int> &ePhiMap,
0037     std::vector<unsigned int> &hPhiMap,
0038     std::vector<unsigned int> &hfPhiMap,
0039     bool useLSB,
0040     bool useCalib,
0041     bool useECALLUT,
0042     bool useHCALLUT,
0043     bool useHFLUT,
0044     bool useHCALFBLUT,
0045     int fwVersion) {
0046   int hfValid = 1;
0047   const HcalTrigTowerGeometry &pG = iSetup.getData(iTokens.geom_);
0048   if (!pG.use1x1()) {
0049     edm::LogError("L1TCaloLayer1FetchLUTs")
0050         << "Using Stage2-Layer1 but HCAL Geometry has use1x1 = 0! HF will be suppressed.  Check Global Tag, etc.";
0051     hfValid = 0;
0052   }
0053 
0054   // CaloParams contains all persisted parameters for Layer 1
0055   edm::ESHandle<l1t::CaloParams> paramsHandle = iSetup.getHandle(iTokens.params_);
0056   if (not paramsHandle.isValid()) {
0057     edm::LogError("L1TCaloLayer1FetchLUTs") << "Missing CaloParams object! Check Global Tag, etc.";
0058     return false;
0059   }
0060   l1t::CaloParamsHelper caloParams(*paramsHandle.product());
0061 
0062   // Calo Trigger Layer1 output LSB Real ET value
0063   double caloLSB = caloParams.towerLsbSum();
0064   if (caloLSB != 0.5) {
0065     // Lots of things expect this, better give fair warning if not
0066     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloLSB (caloParams.towerLsbSum()) != 0.5, actually = " << caloLSB;
0067   }
0068 
0069   // ECal/HCal scale factors will be a x*y*28 array:
0070   //   ieta = 28 eta scale factors (1 .. 28)
0071   //   etBin = size of Real ET Bins vector
0072   //   phiBin = max(Real Phi Bins vector)
0073   //   So, index = phiBin*etBin*28+etBin*28+ieta
0074   auto ecalScaleETBins = caloParams.layer1ECalScaleETBins();
0075   auto ecalScalePhiBins = caloParams.layer1ECalScalePhiBins();
0076   if (ecalScalePhiBins.empty()) {
0077     // Backwards-compatibility (no phi binning)
0078     ecalScalePhiBins.resize(36, 0);
0079   } else if (ecalScalePhiBins.size() % 36 != 0) {
0080     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1ECalScaleETBins().size() is not multiple of 36 !!";
0081     return false;
0082   }
0083   size_t numEcalPhiBins = (*std::max_element(ecalScalePhiBins.begin(), ecalScalePhiBins.end())) + 1;
0084   auto ecalSF = caloParams.layer1ECalScaleFactors();
0085   auto ecalZSF = caloParams.layer1ECalZSFactors();
0086   if (ecalSF.size() != ecalScaleETBins.size() * numEcalPhiBins * 28) {
0087     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1ECalScaleFactors().size() != "
0088                                                "caloParams.layer1ECalScaleETBins().size()*numEcalPhiBins*28 !!";
0089     return false;
0090   }
0091   auto hcalScaleETBins = caloParams.layer1HCalScaleETBins();
0092   auto hcalScalePhiBins = caloParams.layer1HCalScalePhiBins();
0093   if (hcalScalePhiBins.empty()) {
0094     hcalScalePhiBins.resize(36, 0);
0095   } else if (hcalScalePhiBins.size() % 36 != 0) {
0096     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1HCalScaleETBins().size() is not multiple of 36 !!";
0097     return false;
0098   }
0099   size_t numHcalPhiBins = (*std::max_element(hcalScalePhiBins.begin(), hcalScalePhiBins.end())) + 1;
0100   auto hcalSF = caloParams.layer1HCalScaleFactors();
0101   auto hcalZSF = caloParams.layer1HCalZSFactors();
0102   if (hcalSF.size() != hcalScaleETBins.size() * numHcalPhiBins * 28) {
0103     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1HCalScaleFactors().size() != "
0104                                                "caloParams.layer1HCalScaleETBins().size()*numHcalPhiBins*28 !!";
0105     return false;
0106   }
0107 
0108   // HF 1x1 scale factors will be a x*y*12 array:
0109   //   ieta = 12 eta scale factors (30 .. 41)
0110   //   etBin = size of Real ET Bins vector
0111   //   phiBin = max(Real Phi Bins vector)
0112   //   So, index = phiBin*etBin*12+etBin*12+ieta
0113   auto hfScaleETBins = caloParams.layer1HFScaleETBins();
0114   auto hfScalePhiBins = caloParams.layer1HFScalePhiBins();
0115   if (hfScalePhiBins.empty()) {
0116     hfScalePhiBins.resize(36, 0);
0117   } else if (hfScalePhiBins.size() % 36 != 0) {
0118     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1HFScaleETBins().size() is not multiple of 36 !!";
0119     return false;
0120   }
0121   size_t numHFPhiBins = (*std::max_element(hfScalePhiBins.begin(), hfScalePhiBins.end())) + 1;
0122   auto hfSF = caloParams.layer1HFScaleFactors();
0123   if (hfSF.size() != hfScaleETBins.size() * numHFPhiBins * 12) {
0124     edm::LogError("L1TCaloLayer1FetchLUTs")
0125         << "caloParams.layer1HFScaleFactors().size() != caloParams.layer1HFScaleETBins().size()*numHFPhiBins*12 !!";
0126     return false;
0127   }
0128 
0129   // Sanity check scale factors exist
0130   if (useCalib && (ecalSF.empty() || hcalSF.empty() || hfSF.empty())) {
0131     edm::LogError("L1TCaloLayer1FetchLUTs") << "Layer 1 calibrations requested (useCalib = True) but there are missing "
0132                                                "scale factors in CaloParams!  Please check conditions setup.";
0133     return false;
0134   }
0135 
0136   // HCAL FB LUT will be a 1*28 array:
0137   //   ieta = 28 eta scale factors (1 .. 28)
0138   //   So, index = ieta
0139   auto fbLUTUpper = caloParams.layer1HCalFBLUTUpper();
0140   auto fbLUTLower = caloParams.layer1HCalFBLUTLower();
0141   // Only check for HCAL FB LUT if useHCALFBLUT = true
0142   if (useHCALFBLUT) {
0143     if (fbLUTUpper.size() != nCalEtaBins) {
0144       edm::LogError("L1TCaloLayer1FetchLUTs")
0145           << "caloParams.layer1HCalFBLUTUpper().size() " << fbLUTUpper.size() << " != " << nCalEtaBins << " !!";
0146       return false;
0147     }
0148     if (fbLUTLower.size() != nCalEtaBins) {
0149       edm::LogError("L1TCaloLayer1FetchLUTs")
0150           << "caloParams.layer1HCalFBLUTLower().size() " << fbLUTLower.size() << " != " << nCalEtaBins << " !!";
0151       return false;
0152     }
0153   }
0154 
0155   // get energy scale to convert input from ECAL - this should be linear with LSB = 0.5 GeV
0156   const double ecalLSB = 0.5;
0157 
0158   // get energy scale to convert input from HCAL
0159   edm::ESHandle<CaloTPGTranscoder> decoder = iSetup.getHandle(iTokens.decoder_);
0160   if (not decoder.isValid()) {
0161     edm::LogError("L1TCaloLayer1FetchLUTs") << "Missing CaloTPGTranscoder object! Check Global Tag, etc.";
0162     return false;
0163   }
0164 
0165   // TP compression scale is always phi symmetric
0166   // We default to 3 since HF has no ieta=41 iphi=1,2
0167   auto decodeHcalEt = [&decoder](int iEta, uint32_t compressedEt, uint32_t iPhi = 3) -> double {
0168     HcalTriggerPrimitiveSample sample(compressedEt);
0169     HcalTrigTowerDetId id(iEta, iPhi);
0170     if (std::abs(iEta) >= 30) {
0171       id.setVersion(1);
0172     }
0173     return decoder->hcaletValue(id, sample);
0174   };
0175 
0176   // Make ECal LUT
0177   for (uint32_t phiBin = 0; phiBin < numEcalPhiBins; phiBin++) {
0178     std::array<std::array<std::array<uint32_t, nEtBins>, nCalSideBins>, nCalEtaBins> phiLUT;
0179 
0180     for (uint32_t etaBin = 0; etaBin < nCalEtaBins; etaBin++) {
0181       for (uint32_t fb = 0; fb < nCalSideBins; fb++) {
0182         for (uint32_t ecalInput = 0; ecalInput <= 0xFF; ecalInput++) {
0183           uint32_t value = ecalInput;
0184           if (useECALLUT) {
0185             double linearizedECalInput = ecalInput * ecalLSB;  // in GeV
0186 
0187             uint32_t etBin = 0;
0188             for (; etBin < ecalScaleETBins.size(); etBin++) {
0189               if (linearizedECalInput < ecalScaleETBins[etBin])
0190                 break;
0191             }
0192             if (etBin >= ecalScaleETBins.size())
0193               etBin = ecalScaleETBins.size() - 1;
0194 
0195             double calibratedECalInput = linearizedECalInput;
0196             if (useCalib)
0197               calibratedECalInput *= ecalSF.at(phiBin * ecalScaleETBins.size() * 28 + etBin * 28 + etaBin);
0198             if (useLSB)
0199               calibratedECalInput /= caloLSB;
0200             if (ecalZSF.size() >= 28 && ecalInput < ecalZSF.size() / 28 - 1)
0201               calibratedECalInput *= ecalZSF.at(ecalInput * 28 + etaBin);
0202 
0203             value = calibratedECalInput;
0204             if (fwVersion > 2) {
0205               // Saturate if either decompressed value is over 127.5 GeV or input saturated
0206               // (meaningless for ecal, since ecalLSB == caloLSB)
0207               if (value > 0xFF || ecalInput == 0xFF) {
0208                 value = 0xFF;
0209               }
0210             } else {
0211               if (value > 0xFF) {
0212                 value = 0xFF;
0213               }
0214             }
0215           }
0216           if (value == 0) {
0217             value = (1 << 11);
0218           } else {
0219             uint32_t et_log2 = ((uint32_t)log2(value)) & 0x7;
0220             value |= (et_log2 << 12);
0221           }
0222           value |= (fb << 10);
0223           phiLUT[etaBin][fb][ecalInput] = value;
0224         }
0225       }
0226     }
0227 
0228     eLUT.push_back(phiLUT);
0229   }
0230 
0231   // Make HCal LUT
0232   for (uint32_t phiBin = 0; phiBin < numHcalPhiBins; phiBin++) {
0233     std::array<std::array<std::array<uint32_t, nEtBins>, nCalSideBins>, nCalEtaBins> phiLUT;
0234     for (uint32_t etaBin = 0; etaBin < nCalEtaBins; etaBin++) {
0235       int caloEta = etaBin + 1;
0236       int iPhi = 3;
0237       auto pos = std::find(hcalScalePhiBins.begin(), hcalScalePhiBins.end(), phiBin);
0238       if (pos != hcalScalePhiBins.end()) {
0239         // grab an iPhi bin
0240         auto index = std::distance(hcalScalePhiBins.begin(), pos);
0241         if (index < 18) {
0242           caloEta *= -1;
0243           iPhi = index * 4 + 1;
0244         } else {
0245           iPhi = (index - 18) * 4 + 1;
0246         }
0247       }
0248       for (uint32_t fb = 0; fb < nCalSideBins; fb++) {
0249         for (uint32_t hcalInput = 0; hcalInput <= 0xFF; hcalInput++) {
0250           uint32_t value = hcalInput;
0251           if (useHCALLUT) {
0252             // hcaletValue defined in L137 of CalibCalorimetry/CaloTPG/src/CaloTPGTranscoderULUT.cc
0253             double linearizedHcalInput = decodeHcalEt(caloEta, hcalInput, iPhi);  // in GeV
0254 
0255             uint32_t etBin = 0;
0256             for (; etBin < hcalScaleETBins.size(); etBin++) {
0257               if (linearizedHcalInput < hcalScaleETBins[etBin])
0258                 break;
0259             }
0260             if (etBin >= hcalScaleETBins.size())
0261               etBin = hcalScaleETBins.size() - 1;
0262 
0263             double calibratedHcalInput = linearizedHcalInput;
0264             if (useCalib)
0265               calibratedHcalInput *= hcalSF.at(phiBin * hcalScaleETBins.size() * 28 + etBin * 28 + etaBin);
0266             if (useLSB)
0267               calibratedHcalInput /= caloLSB;
0268             if (hcalZSF.size() >= 28 && hcalInput < hcalZSF.size() / 28 - 1)
0269               calibratedHcalInput *= hcalZSF.at(hcalInput * 28 + etaBin);
0270             value = calibratedHcalInput;
0271             if (fwVersion > 2) {
0272               // Saturate if either decompressed value is over 127.5 GeV or input saturated
0273               if (value > 0xFF || hcalInput == 0xFF) {
0274                 value = 0xFF;
0275               }
0276             } else {
0277               if (value > 0xFF) {
0278                 value = 0xFF;
0279               }
0280             }
0281           }
0282           if (value == 0) {
0283             value = (1 << 11);
0284           } else {
0285             uint32_t et_log2 = ((uint32_t)log2(value)) & 0x7;
0286             value |= (et_log2 << 12);
0287           }
0288           value |= (fb << 10);
0289           phiLUT[etaBin][fb][hcalInput] = value;
0290         }
0291       }
0292     }
0293 
0294     hLUT.push_back(phiLUT);
0295   }
0296 
0297   // Make HF LUT
0298   for (uint32_t phiBin = 0; phiBin < numHFPhiBins; phiBin++) {
0299     std::array<std::array<uint32_t, nEtBins>, nHfEtaBins> phiLUT;
0300 
0301     for (uint32_t etaBin = 0; etaBin < nHfEtaBins; etaBin++) {
0302       int caloEta = etaBin + 30;
0303       int iPhi = 3;
0304       auto pos = std::find(hfScalePhiBins.begin(), hfScalePhiBins.end(), phiBin);
0305       if (pos != hfScalePhiBins.end()) {
0306         auto index = std::distance(hfScalePhiBins.begin(), pos);
0307         if (index < 18) {
0308           caloEta *= -1;
0309           iPhi = index * 4 - 1;
0310         } else {
0311           iPhi = (index - 18) * 4 - 1;
0312         }
0313         if (iPhi < 0)
0314           iPhi = 71;
0315       }
0316       for (uint32_t etCode = 0; etCode < nEtBins; etCode++) {
0317         uint32_t value = etCode;
0318         if (useHFLUT) {
0319           double linearizedHFInput = 0;
0320           if (hfValid) {
0321             linearizedHFInput = decodeHcalEt(caloEta, value, iPhi);  // in GeV
0322           }
0323 
0324           uint32_t etBin = 0;
0325           for (; etBin < hfScaleETBins.size(); etBin++) {
0326             if (linearizedHFInput < hfScaleETBins[etBin])
0327               break;
0328           }
0329           if (etBin >= hfScaleETBins.size())
0330             etBin = hfScaleETBins.size() - 1;
0331 
0332           double calibratedHFInput = linearizedHFInput;
0333           if (useCalib)
0334             calibratedHFInput *= hfSF.at(phiBin * hfScalePhiBins.size() * 12 + etBin * 12 + etaBin);
0335           if (useLSB)
0336             calibratedHFInput /= caloLSB;
0337 
0338           if (fwVersion > 2) {
0339             uint32_t absCaloEta = std::abs(caloEta);
0340             if (absCaloEta > 29 && absCaloEta < 40) {
0341               // Divide by two (since two duplicate towers are sent)
0342               calibratedHFInput *= 0.5;
0343             } else if (absCaloEta == 40 || absCaloEta == 41) {
0344               // Divide by four
0345               calibratedHFInput *= 0.25;
0346             }
0347             value = calibratedHFInput;
0348             // Saturate if either decompressed value is over 127.5 GeV or input saturated
0349             if (value >= 0xFF || etCode == 0xFF) {
0350               value = 0x1FD;
0351             }
0352           } else {
0353             value = calibratedHFInput;
0354             if (value > 0xFF) {
0355               value = 0xFF;
0356             }
0357           }
0358         }
0359         phiLUT[etaBin][etCode] = value;
0360       }
0361     }
0362     hfLUT.push_back(phiLUT);
0363   }
0364 
0365   // Make HCal FB LUT
0366   for (uint32_t etaBin = 0; etaBin < nCalEtaBins; etaBin++) {
0367     uint64_t value = 0xFFFFFFFFFFFFFFFF;
0368     if (useHCALFBLUT) {
0369       value = (((uint64_t)fbLUTUpper.at(etaBin)) << 32) | fbLUTLower.at(etaBin);
0370     }
0371     hcalFBLUT.push_back(value);
0372   }
0373 
0374   // plus/minus, 18 CTP7, 4 iPhi each
0375   for (uint32_t isPos = 0; isPos < 2; isPos++) {
0376     for (uint32_t iPhi = 1; iPhi <= 72; iPhi++) {
0377       uint32_t card = floor((iPhi + 1) / 4);
0378       if (card > 17)
0379         card -= 18;
0380       ePhiMap[isPos * 72 + iPhi - 1] = ecalScalePhiBins[isPos * 18 + card];
0381       hPhiMap[isPos * 72 + iPhi - 1] = hcalScalePhiBins[isPos * 18 + card];
0382       hfPhiMap[isPos * 72 + iPhi - 1] = hfScalePhiBins[isPos * 18 + card];
0383     }
0384   }
0385 
0386   return true;
0387 }
0388 /* vim: set ts=8 sw=2 tw=0 et :*/