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