L1TCaloLayer1/src/L1TCaloLayer1FetchLUTs.cc

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   if (ecalSF.size() != ecalScaleETBins.size() * numEcalPhiBins * 28) {
0086     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1ECalScaleFactors().size() != "
0087                                                "caloParams.layer1ECalScaleETBins().size()*numEcalPhiBins*28 !!";
0088     return false;
0089   }
0090   auto hcalScaleETBins = caloParams.layer1HCalScaleETBins();
0091   auto hcalScalePhiBins = caloParams.layer1HCalScalePhiBins();
0092   if (hcalScalePhiBins.empty()) {
0093     hcalScalePhiBins.resize(36, 0);
0094   } else if (hcalScalePhiBins.size() % 36 != 0) {
0095     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1HCalScaleETBins().size() is not multiple of 36 !!";
0096     return false;
0097   }
0098   size_t numHcalPhiBins = (*std::max_element(hcalScalePhiBins.begin(), hcalScalePhiBins.end())) + 1;
0099   auto hcalSF = caloParams.layer1HCalScaleFactors();
0100   if (hcalSF.size() != hcalScaleETBins.size() * numHcalPhiBins * 28) {
0101     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1HCalScaleFactors().size() != "
0102                                                "caloParams.layer1HCalScaleETBins().size()*numHcalPhiBins*28 !!";
0103     return false;
0104   }
0105
0106   // HF 1x1 scale factors will be a x*y*12 array:
0107   //   ieta = 12 eta scale factors (30 .. 41)
0108   //   etBin = size of Real ET Bins vector
0109   //   phiBin = max(Real Phi Bins vector)
0110   //   So, index = phiBin*etBin*12+etBin*12+ieta
0111   auto hfScaleETBins = caloParams.layer1HFScaleETBins();
0112   auto hfScalePhiBins = caloParams.layer1HFScalePhiBins();
0113   if (hfScalePhiBins.empty()) {
0114     hfScalePhiBins.resize(36, 0);
0115   } else if (hfScalePhiBins.size() % 36 != 0) {
0116     edm::LogError("L1TCaloLayer1FetchLUTs") << "caloParams.layer1HFScaleETBins().size() is not multiple of 36 !!";
0117     return false;
0118   }
0119   size_t numHFPhiBins = (*std::max_element(hfScalePhiBins.begin(), hfScalePhiBins.end())) + 1;
0120   auto hfSF = caloParams.layer1HFScaleFactors();
0121   if (hfSF.size() != hfScaleETBins.size() * numHFPhiBins * 12) {
0122     edm::LogError("L1TCaloLayer1FetchLUTs")
0123         << "caloParams.layer1HFScaleFactors().size() != caloParams.layer1HFScaleETBins().size()*numHFPhiBins*12 !!";
0124     return false;
0125   }
0126
0127   // Sanity check scale factors exist
0128   if (useCalib && (ecalSF.empty() || hcalSF.empty() || hfSF.empty())) {
0129     edm::LogError("L1TCaloLayer1FetchLUTs") << "Layer 1 calibrations requested (useCalib = True) but there are missing "
0130                                                "scale factors in CaloParams!  Please check conditions setup.";
0131     return false;
0132   }
0133
0134   // HCAL FB LUT will be a 1*28 array:
0135   //   ieta = 28 eta scale factors (1 .. 28)
0136   //   So, index = ieta
0137   auto fbLUTUpper = caloParams.layer1HCalFBLUTUpper();
0138   auto fbLUTLower = caloParams.layer1HCalFBLUTLower();
0139   // Only check for HCAL FB LUT if useHCALFBLUT = true
0140   if (useHCALFBLUT) {
0141     if (fbLUTUpper.size() != nCalEtaBins) {
0142       edm::LogError("L1TCaloLayer1FetchLUTs")
0143           << "caloParams.layer1HCalFBLUTUpper().size() " << fbLUTUpper.size() << " != " << nCalEtaBins << " !!";
0144       return false;
0145     }
0146     if (fbLUTLower.size() != nCalEtaBins) {
0147       edm::LogError("L1TCaloLayer1FetchLUTs")
0148           << "caloParams.layer1HCalFBLUTLower().size() " << fbLUTLower.size() << " != " << nCalEtaBins << " !!";
0149       return false;
0150     }
0151   }
0152
0153   // get energy scale to convert input from ECAL - this should be linear with LSB = 0.5 GeV
0154   const double ecalLSB = 0.5;
0155
0156   // get energy scale to convert input from HCAL
0157   edm::ESHandle<CaloTPGTranscoder> decoder = iSetup.getHandle(iTokens.decoder_);
0158   if (not decoder.isValid()) {
0159     edm::LogError("L1TCaloLayer1FetchLUTs") << "Missing CaloTPGTranscoder object! Check Global Tag, etc.";
0160     return false;
0161   }
0162
0163   // TP compression scale is always phi symmetric
0164   // We default to 3 since HF has no ieta=41 iphi=1,2
0165   auto decodeHcalEt = [&decoder](int iEta, uint32_t compressedEt, uint32_t iPhi = 3) -> double {
0166     HcalTriggerPrimitiveSample sample(compressedEt);
0167     HcalTrigTowerDetId id(iEta, iPhi);
0168     if (std::abs(iEta) >= 30) {
0169       id.setVersion(1);
0170     }
0171     return decoder->hcaletValue(id, sample);
0172   };
0173
0174   // Make ECal LUT
0175   for (uint32_t phiBin = 0; phiBin < numEcalPhiBins; phiBin++) {
0176     std::array<std::array<std::array<uint32_t, nEtBins>, nCalSideBins>, nCalEtaBins> phiLUT;
0177
0178     for (uint32_t etaBin = 0; etaBin < nCalEtaBins; etaBin++) {
0179       for (uint32_t fb = 0; fb < nCalSideBins; fb++) {
0180         for (uint32_t ecalInput = 0; ecalInput <= 0xFF; ecalInput++) {
0181           uint32_t value = ecalInput;
0182           if (useECALLUT) {
0183             double linearizedECalInput = ecalInput * ecalLSB;  // in GeV
0184
0185             uint32_t etBin = 0;
0186             for (; etBin < ecalScaleETBins.size(); etBin++) {
0187               if (linearizedECalInput < ecalScaleETBins[etBin])
0188                 break;
0189             }
0190             if (etBin >= ecalScaleETBins.size())
0191               etBin = ecalScaleETBins.size() - 1;
0192
0193             double calibratedECalInput = linearizedECalInput;
0194             if (useCalib)
0195               calibratedECalInput *= ecalSF.at(phiBin * ecalScaleETBins.size() * 28 + etBin * 28 + etaBin);
0196             if (useLSB)
0197               calibratedECalInput /= caloLSB;
0198
0199             value = calibratedECalInput;
0200             if (fwVersion > 2) {
0201               // Saturate if either decompressed value is over 127.5 GeV or input saturated
0202               // (meaningless for ecal, since ecalLSB == caloLSB)
0203               if (value > 0xFF || ecalInput == 0xFF) {
0204                 value = 0xFF;
0205               }
0206             } else {
0207               if (value > 0xFF) {
0208                 value = 0xFF;
0209               }
0210             }
0211           }
0212           if (value == 0) {
0213             value = (1 << 11);
0214           } else {
0215             uint32_t et_log2 = ((uint32_t)log2(value)) & 0x7;
0216             value |= (et_log2 << 12);
0217           }
0218           value |= (fb << 10);
0219           phiLUT[etaBin][fb][ecalInput] = value;
0220         }
0221       }
0222     }
0223
0224     eLUT.push_back(phiLUT);
0225   }
0226
0227   // Make HCal LUT
0228   for (uint32_t phiBin = 0; phiBin < numHcalPhiBins; phiBin++) {
0229     std::array<std::array<std::array<uint32_t, nEtBins>, nCalSideBins>, nCalEtaBins> phiLUT;
0230     for (uint32_t etaBin = 0; etaBin < nCalEtaBins; etaBin++) {
0231       int caloEta = etaBin + 1;
0232       int iPhi = 3;
0233       auto pos = std::find(hcalScalePhiBins.begin(), hcalScalePhiBins.end(), phiBin);
0234       if (pos != hcalScalePhiBins.end()) {
0235         // grab an iPhi bin
0236         auto index = std::distance(hcalScalePhiBins.begin(), pos);
0237         if (index < 18) {
0238           caloEta *= -1;
0239           iPhi = index * 4 + 1;
0240         } else {
0241           iPhi = (index - 18) * 4 + 1;
0242         }
0243       }
0244       for (uint32_t fb = 0; fb < nCalSideBins; fb++) {
0245         for (uint32_t hcalInput = 0; hcalInput <= 0xFF; hcalInput++) {
0246           uint32_t value = hcalInput;
0247           if (useHCALLUT) {
0248             // hcaletValue defined in L137 of CalibCalorimetry/CaloTPG/src/CaloTPGTranscoderULUT.cc
0249             double linearizedHcalInput = decodeHcalEt(caloEta, hcalInput, iPhi);  // in GeV
0250
0251             uint32_t etBin = 0;
0252             for (; etBin < hcalScaleETBins.size(); etBin++) {
0253               if (linearizedHcalInput < hcalScaleETBins[etBin])
0254                 break;
0255             }
0256             if (etBin >= hcalScaleETBins.size())
0257               etBin = hcalScaleETBins.size() - 1;
0258
0259             double calibratedHcalInput = linearizedHcalInput;
0260             if (useCalib)
0261               calibratedHcalInput *= hcalSF.at(phiBin * hcalScaleETBins.size() * 28 + etBin * 28 + etaBin);
0262             if (useLSB)
0263               calibratedHcalInput /= caloLSB;
0264
0265             value = calibratedHcalInput;
0266             if (fwVersion > 2) {
0267               // Saturate if either decompressed value is over 127.5 GeV or input saturated
0268               if (value > 0xFF || hcalInput == 0xFF) {
0269                 value = 0xFF;
0270               }
0271             } else {
0272               if (value > 0xFF) {
0273                 value = 0xFF;
0274               }
0275             }
0276           }
0277           if (value == 0) {
0278             value = (1 << 11);
0279           } else {
0280             uint32_t et_log2 = ((uint32_t)log2(value)) & 0x7;
0281             value |= (et_log2 << 12);
0282           }
0283           value |= (fb << 10);
0284           phiLUT[etaBin][fb][hcalInput] = value;
0285         }
0286       }
0287     }
0288
0289     hLUT.push_back(phiLUT);
0290   }
0291
0292   // Make HF LUT
0293   for (uint32_t phiBin = 0; phiBin < numHFPhiBins; phiBin++) {
0294     std::array<std::array<uint32_t, nEtBins>, nHfEtaBins> phiLUT;
0295
0296     for (uint32_t etaBin = 0; etaBin < nHfEtaBins; etaBin++) {
0297       int caloEta = etaBin + 30;
0298       int iPhi = 3;
0299       auto pos = std::find(hfScalePhiBins.begin(), hfScalePhiBins.end(), phiBin);
0300       if (pos != hfScalePhiBins.end()) {
0301         auto index = std::distance(hfScalePhiBins.begin(), pos);
0302         if (index < 18) {
0303           caloEta *= -1;
0304           iPhi = index * 4 - 1;
0305         } else {
0306           iPhi = (index - 18) * 4 - 1;
0307         }
0308         if (iPhi < 0)
0309           iPhi = 71;
0310       }
0311       for (uint32_t etCode = 0; etCode < nEtBins; etCode++) {
0312         uint32_t value = etCode;
0313         if (useHFLUT) {
0314           double linearizedHFInput = 0;
0315           if (hfValid) {
0316             linearizedHFInput = decodeHcalEt(caloEta, value, iPhi);  // in GeV
0317           }
0318
0319           uint32_t etBin = 0;
0320           for (; etBin < hfScaleETBins.size(); etBin++) {
0321             if (linearizedHFInput < hfScaleETBins[etBin])
0322               break;
0323           }
0324           if (etBin >= hfScaleETBins.size())
0325             etBin = hfScaleETBins.size() - 1;
0326
0327           double calibratedHFInput = linearizedHFInput;
0328           if (useCalib)
0329             calibratedHFInput *= hfSF.at(phiBin * hfScalePhiBins.size() * 12 + etBin * 12 + etaBin);
0330           if (useLSB)
0331             calibratedHFInput /= caloLSB;
0332
0333           if (fwVersion > 2) {
0334             uint32_t absCaloEta = std::abs(caloEta);
0335             if (absCaloEta > 29 && absCaloEta < 40) {
0336               // Divide by two (since two duplicate towers are sent)
0337               calibratedHFInput *= 0.5;
0338             } else if (absCaloEta == 40 || absCaloEta == 41) {
0339               // Divide by four
0340               calibratedHFInput *= 0.25;
0341             }
0342             value = calibratedHFInput;
0343             // Saturate if either decompressed value is over 127.5 GeV or input saturated
0344             if (value >= 0xFF || etCode == 0xFF) {
0345               value = 0x1FD;
0346             }
0347           } else {
0348             value = calibratedHFInput;
0349             if (value > 0xFF) {
0350               value = 0xFF;
0351             }
0352           }
0353         }
0354         phiLUT[etaBin][etCode] = value;
0355       }
0356     }
0357     hfLUT.push_back(phiLUT);
0358   }
0359
0360   // Make HCal FB LUT
0361   for (uint32_t etaBin = 0; etaBin < nCalEtaBins; etaBin++) {
0362     uint64_t value = 0xFFFFFFFFFFFFFFFF;
0363     if (useHCALFBLUT) {
0364       value = (((uint64_t)fbLUTUpper.at(etaBin)) << 32) | fbLUTLower.at(etaBin);
0365     }
0366     hcalFBLUT.push_back(value);
0367   }
0368
0369   // plus/minus, 18 CTP7, 4 iPhi each
0370   for (uint32_t isPos = 0; isPos < 2; isPos++) {
0371     for (uint32_t iPhi = 1; iPhi <= 72; iPhi++) {
0372       uint32_t card = floor((iPhi + 1) / 4);
0373       if (card > 17)
0374         card -= 18;
0375       ePhiMap[isPos * 72 + iPhi - 1] = ecalScalePhiBins[isPos * 18 + card];
0376       hPhiMap[isPos * 72 + iPhi - 1] = hcalScalePhiBins[isPos * 18 + card];
0377       hfPhiMap[isPos * 72 + iPhi - 1] = hfScalePhiBins[isPos * 18 + card];
0378     }
0379   }
0380
0381   return true;
0382 }
0383 /* vim: set ts=8 sw=2 tw=0 et :*/