HcalTPGAlgos/src/HcaluLUTTPGCoder.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include <cmath>
0004 #include <string>
0005 #include "CalibCalorimetry/HcalTPGAlgos/interface/HcaluLUTTPGCoder.h"
0006 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0007 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0008 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0009 #include "DataFormats/HcalDigi/interface/HcalQIENum.h"
0010 #include "DataFormats/HcalDigi/interface/QIE10DataFrame.h"
0011 #include "DataFormats/HcalDigi/interface/QIE11DataFrame.h"
0012 #include "DataFormats/HcalDetId/interface/HcalTrigTowerDetId.h"
0013 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0016 #include "FWCore/Utilities/interface/Exception.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "CondFormats/HcalObjects/interface/HcalL1TriggerObjects.h"
0024 #include "CondFormats/HcalObjects/interface/HcalL1TriggerObject.h"
0025 #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
0026 #include "CalibCalorimetry/HcalAlgos/interface/HcalSiPMnonlinearity.h"
0027 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentCorrection.h"
0028 #include "CalibCalorimetry/HcalTPGAlgos/interface/XMLProcessor.h"
0029 #include "CalibCalorimetry/HcalTPGAlgos/interface/LutXml.h"
0030
0031 const float HcaluLUTTPGCoder::lsb_ = 1. / 16;
0032
0033 const int HcaluLUTTPGCoder::QIE8_LUT_BITMASK;
0034 const int HcaluLUTTPGCoder::QIE10_LUT_BITMASK;
0035 const int HcaluLUTTPGCoder::QIE11_LUT_BITMASK;
0036
0037 constexpr double MaximumFractionalError = 0.002;  // 0.2% error allowed from this source
0038
0039 HcaluLUTTPGCoder::HcaluLUTTPGCoder()
0040     : topo_{},
0041       delay_{},
0042       LUTGenerationMode_{},
0043       FG_HF_thresholds_{},
0044       bitToMask_{},
0045       firstHBEta_{},
0046       lastHBEta_{},
0047       nHBEta_{},
0048       maxDepthHB_{},
0049       sizeHB_{},
0050       firstHEEta_{},
0051       lastHEEta_{},
0052       nHEEta_{},
0053       maxDepthHE_{},
0054       sizeHE_{},
0055       firstHFEta_{},
0056       lastHFEta_{},
0057       nHFEta_{},
0058       maxDepthHF_{},
0059       sizeHF_{},
0060       cosh_ieta_28_HE_low_depths_{},
0061       cosh_ieta_28_HE_high_depths_{},
0062       cosh_ieta_29_HE_{},
0063       allLinear_{},
0064       contain1TSHB_{},
0065       contain1TSHE_{},
0066       applyFixPCC_{},
0067       linearLSB_QIE8_{},
0068       linearLSB_QIE11_{},
0069       linearLSB_QIE11Overlap_{} {}
0070
0071 HcaluLUTTPGCoder::HcaluLUTTPGCoder(const HcalTopology* top, const HcalTimeSlew* delay) { init(top, delay); }
0072
0073 void HcaluLUTTPGCoder::init(const HcalTopology* top, const HcalTimeSlew* delay) {
0074   topo_ = top;
0075   delay_ = delay;
0076   LUTGenerationMode_ = true;
0077   FG_HF_thresholds_ = {0, 0};
0078   bitToMask_ = 0;
0079   allLinear_ = false;
0080   contain1TSHB_ = false;
0081   contain1TSHE_ = false;
0082   applyFixPCC_ = false;
0083   linearLSB_QIE8_ = 1.;
0084   linearLSB_QIE11_ = 1.;
0085   linearLSB_QIE11Overlap_ = 1.;
0086   pulseCorr_ = std::make_unique<HcalPulseContainmentManager>(MaximumFractionalError, false);
0087   firstHBEta_ = topo_->firstHBRing();
0088   lastHBEta_ = topo_->lastHBRing();
0089   nHBEta_ = (lastHBEta_ - firstHBEta_ + 1);
0090   maxDepthHB_ = topo_->maxDepth(HcalBarrel);
0091   sizeHB_ = 2 * nHBEta_ * nFi_ * maxDepthHB_;
0092   firstHEEta_ = topo_->firstHERing();
0093   lastHEEta_ = topo_->lastHERing();
0094   nHEEta_ = (lastHEEta_ - firstHEEta_ + 1);
0095   maxDepthHE_ = topo_->maxDepth(HcalEndcap);
0096   sizeHE_ = 2 * nHEEta_ * nFi_ * maxDepthHE_;
0097   firstHFEta_ = topo_->firstHFRing();
0098   lastHFEta_ = topo_->lastHFRing();
0099   nHFEta_ = (lastHFEta_ - firstHFEta_ + 1);
0100   maxDepthHF_ = topo_->maxDepth(HcalForward);
0101   sizeHF_ = 2 * nHFEta_ * nFi_ * maxDepthHF_;
0102   size_t nluts = (size_t)(sizeHB_ + sizeHE_ + sizeHF_ + 1);
0103   inputLUT_ = std::vector<HcaluLUTTPGCoder::Lut>(nluts);
0104   gain_ = std::vector<float>(nluts, 0.);
0105   ped_ = std::vector<float>(nluts, 0.);
0106   make_cosh_ieta_map();
0107 }
0108
0109 void HcaluLUTTPGCoder::compress(const IntegerCaloSamples& ics,
0110                                 const std::vector<bool>& featureBits,
0111                                 HcalTriggerPrimitiveDigi& tp) const {
0112   throw cms::Exception("PROBLEM: This method should never be invoked!");
0113 }
0114
0115 HcaluLUTTPGCoder::~HcaluLUTTPGCoder() {}
0116
0117 int HcaluLUTTPGCoder::getLUTId(HcalSubdetector id, int ieta, int iphi, int depth) const {
0118   int retval(0);
0119   if (id == HcalBarrel) {
0120     retval = (depth - 1) + maxDepthHB_ * (iphi - 1);
0121     if (ieta > 0)
0122       retval += maxDepthHB_ * nFi_ * (ieta - firstHBEta_);
0123     else
0124       retval += maxDepthHB_ * nFi_ * (ieta + lastHBEta_ + nHBEta_);
0125   } else if (id == HcalEndcap) {
0126     retval = sizeHB_;
0127     retval += (depth - 1) + maxDepthHE_ * (iphi - 1);
0128     if (ieta > 0)
0129       retval += maxDepthHE_ * nFi_ * (ieta - firstHEEta_);
0130     else
0131       retval += maxDepthHE_ * nFi_ * (ieta + lastHEEta_ + nHEEta_);
0132   } else if (id == HcalForward) {
0133     retval = sizeHB_ + sizeHE_;
0134     retval += (depth - 1) + maxDepthHF_ * (iphi - 1);
0135     if (ieta > 0)
0136       retval += maxDepthHF_ * nFi_ * (ieta - firstHFEta_);
0137     else
0138       retval += maxDepthHF_ * nFi_ * (ieta + lastHFEta_ + nHFEta_);
0139   }
0140   return retval;
0141 }
0142
0143 int HcaluLUTTPGCoder::getLUTId(uint32_t rawid) const {
0144   HcalDetId detid(rawid);
0145   return getLUTId(detid.subdet(), detid.ieta(), detid.iphi(), detid.depth());
0146 }
0147
0148 int HcaluLUTTPGCoder::getLUTId(const HcalDetId& detid) const {
0149   return getLUTId(detid.subdet(), detid.ieta(), detid.iphi(), detid.depth());
0150 }
0151
0152 void HcaluLUTTPGCoder::update(const char* filename, bool appendMSB) {
0153   std::ifstream file(filename, std::ios::in);
0154   assert(file.is_open());
0155
0156   std::vector<HcalSubdetector> subdet;
0157   std::string buffer;
0158
0159   // Drop first (comment) line
0160   std::getline(file, buffer);
0161   std::getline(file, buffer);
0162
0163   unsigned int index = buffer.find('H', 0);
0164   while (index < buffer.length()) {
0165     std::string subdetStr = buffer.substr(index, 2);
0166     if (subdetStr == "HB")
0167       subdet.push_back(HcalBarrel);
0168     else if (subdetStr == "HE")
0169       subdet.push_back(HcalEndcap);
0170     else if (subdetStr == "HF")
0171       subdet.push_back(HcalForward);
0172     //TODO Check subdet
0173     //else exception
0174     index += 2;
0175     index = buffer.find('H', index);
0176   }
0177
0178   // Get upper/lower ranges for ieta/iphi/depth
0179   size_t nCol = subdet.size();
0180   assert(nCol > 0);
0181
0182   std::vector<int> ietaU;
0183   std::vector<int> ietaL;
0184   std::vector<int> iphiU;
0185   std::vector<int> iphiL;
0186   std::vector<int> depU;
0187   std::vector<int> depL;
0188   std::vector<Lut> lutFromFile(nCol);
0189   LutElement lutValue;
0190
0191   for (size_t i = 0; i < nCol; ++i) {
0192     int ieta;
0193     file >> ieta;
0194     ietaL.push_back(ieta);
0195   }
0196
0197   for (size_t i = 0; i < nCol; ++i) {
0198     int ieta;
0199     file >> ieta;
0200     ietaU.push_back(ieta);
0201   }
0202
0203   for (size_t i = 0; i < nCol; ++i) {
0204     int iphi;
0205     file >> iphi;
0206     iphiL.push_back(iphi);
0207   }
0208
0209   for (size_t i = 0; i < nCol; ++i) {
0210     int iphi;
0211     file >> iphi;
0212     iphiU.push_back(iphi);
0213   }
0214
0215   for (size_t i = 0; i < nCol; ++i) {
0216     int dep;
0217     file >> dep;
0218     depL.push_back(dep);
0219   }
0220
0221   for (size_t i = 0; i < nCol; ++i) {
0222     int dep;
0223     file >> dep;
0224     depU.push_back(dep);
0225   }
0226
0227   // Read Lut Entry
0228   for (size_t i = 0; file >> lutValue; i = (i + 1) % nCol) {
0229     lutFromFile[i].push_back(lutValue);
0230   }
0231
0232   // Check lut size
0233   for (size_t i = 0; i < nCol; ++i)
0234     assert(lutFromFile[i].size() == INPUT_LUT_SIZE);
0235
0236   for (size_t i = 0; i < nCol; ++i) {
0237     for (int ieta = ietaL[i]; ieta <= ietaU[i]; ++ieta) {
0238       for (int iphi = iphiL[i]; iphi <= iphiU[i]; ++iphi) {
0239         for (int depth = depL[i]; depth <= depU[i]; ++depth) {
0240           HcalDetId id(subdet[i], ieta, iphi, depth);
0241           if (!topo_->valid(id))
0242             continue;
0243
0244           int lutId = getLUTId(id);
0245           for (size_t adc = 0; adc < INPUT_LUT_SIZE; ++adc) {
0246             if (appendMSB) {
0247               // Append FG bit LUT to MSB
0248               // MSB = Most Significant Bit = bit 10
0249               // Overwrite bit 10
0250               LutElement msb = (lutFromFile[i][adc] != 0 ? QIE8_LUT_MSB : 0);
0251               inputLUT_[lutId][adc] = (msb | (inputLUT_[lutId][adc] & QIE8_LUT_BITMASK));
0252             } else
0253               inputLUT_[lutId][adc] = lutFromFile[i][adc];
0254           }  // for adc
0255         }    // for depth
0256       }      // for iphi
0257     }        // for ieta
0258   }          // for nCol
0259 }
0260
0261 void HcaluLUTTPGCoder::updateXML(const char* filename) {
0262   LutXml* _xml = new LutXml(filename);
0263   _xml->create_lut_map();
0264   HcalSubdetector subdet[3] = {HcalBarrel, HcalEndcap, HcalForward};
0265   for (int ieta = -HcalDetId::kHcalEtaMask2; ieta <= (int)(HcalDetId::kHcalEtaMask2); ++ieta) {
0266     for (unsigned int iphi = 0; iphi <= HcalDetId::kHcalPhiMask2; ++iphi) {
0267       for (unsigned int depth = 1; depth < HcalDetId::kHcalDepthMask2; ++depth) {
0268         for (int isub = 0; isub < 3; ++isub) {
0269           HcalDetId detid(subdet[isub], ieta, iphi, depth);
0270           if (!topo_->valid(detid))
0271             continue;
0272           int id = getLUTId(subdet[isub], ieta, iphi, depth);
0273           std::vector<unsigned int>* lut = _xml->getLutFast(detid);
0274           if (lut == nullptr)
0275             throw cms::Exception("PROBLEM: No inputLUT_ in xml file for ") << detid << std::endl;
0276           if (lut->size() != INPUT_LUT_SIZE)
0277             throw cms::Exception("PROBLEM: Wrong inputLUT_ size in xml file for ") << detid << std::endl;
0278           for (unsigned int i = 0; i < INPUT_LUT_SIZE; ++i)
0279             inputLUT_[id][i] = (LutElement)lut->at(i);
0280         }
0281       }
0282     }
0283   }
0284   delete _xml;
0285   XMLProcessor::getInstance()->terminate();
0286 }
0287
0288 double HcaluLUTTPGCoder::cosh_ieta(int ieta, int depth, HcalSubdetector subdet) {
0289   // ieta = 28 and 29 are both associated with trigger tower 28
0290   // so special handling is required. HF ieta=29 channels included in TT30
0291   // are already handled correctly in cosh_ieta_
0292   if (abs(ieta) >= 28 && subdet == HcalEndcap && allLinear_) {
0293     if (abs(ieta) == 29)
0294       return cosh_ieta_29_HE_;
0295     if (abs(ieta) == 28) {
0296       if (depth <= 3)
0297         return cosh_ieta_28_HE_low_depths_;
0298       else
0299         return cosh_ieta_28_HE_high_depths_;
0300     }
0301   }
0302
0303   return cosh_ieta_[ieta];
0304 }
0305
0306 void HcaluLUTTPGCoder::make_cosh_ieta_map(void) {
0307   cosh_ieta_ = std::vector<double>(lastHFEta_ + 1, -1.0);
0308
0309   HcalTrigTowerGeometry triggeo(topo_);
0310
0311   for (int i = 1; i <= firstHFEta_; ++i) {
0312     double eta_low = 0., eta_high = 0.;
0313     triggeo.towerEtaBounds(i, 0, eta_low, eta_high);
0314     cosh_ieta_[i] = cosh((eta_low + eta_high) / 2.);
0315   }
0316   for (int i = firstHFEta_; i <= lastHFEta_; ++i) {
0317     std::pair<double, double> etas = topo_->etaRange(HcalForward, i);
0318     double eta1 = etas.first;
0319     double eta2 = etas.second;
0320     cosh_ieta_[i] = cosh((eta1 + eta2) / 2.);
0321   }
0322
0323   // trigger tower 28 in HE has a more complicated geometry
0324   std::pair<double, double> eta28 = topo_->etaRange(HcalEndcap, 28);
0325   std::pair<double, double> eta29 = topo_->etaRange(HcalEndcap, 29);
0326   cosh_ieta_29_HE_ = cosh((eta29.first + eta29.second) / 2.);
0327   cosh_ieta_28_HE_low_depths_ = cosh((eta28.first + eta28.second) / 2.);
0328   // for higher depths in ieta = 28, the trigger tower extends past
0329   // the ieta = 29 channels
0330   cosh_ieta_28_HE_high_depths_ = cosh((eta28.first + eta29.second) / 2.);
0331 }
0332
0333 void HcaluLUTTPGCoder::update(const HcalDbService& conditions) {
0334   const HcalLutMetadata* metadata = conditions.getHcalLutMetadata();
0335   assert(metadata != nullptr);
0336   float nominalgain_ = metadata->getNominalGain();
0337
0338   pulseCorr_ = std::make_unique<HcalPulseContainmentManager>(MaximumFractionalError, applyFixPCC_);
0339   pulseCorr_->beginRun(&conditions, delay_);
0340
0341   make_cosh_ieta_map();
0342
0343   // Here we will determine if we are using new version of TPs (1TS)
0344   // i.e. are we using a new pulse filter scheme.
0345   const HcalElectronicsMap* emap = conditions.getHcalMapping();
0346
0347   int lastHBRing = topo_->lastHBRing();
0348   int lastHERing = topo_->lastHERing();
0349
0350   // First, determine if we should configure for the filter scheme
0351   // Check the tp version to make this determination
0352   bool foundHB = false;
0353   bool foundHE = false;
0354   bool newHBtp = false;
0355   bool newHEtp = false;
0356   std::vector<HcalElectronicsId> vIds = emap->allElectronicsIdTrigger();
0357   for (std::vector<HcalElectronicsId>::const_iterator eId = vIds.begin(); eId != vIds.end(); eId++) {
0358     // The first HB or HE id is enough to tell whether to use new scheme in HB or HE
0359     if (foundHB and foundHE)
0360       break;
0361
0362     HcalTrigTowerDetId hcalTTDetId(emap->lookupTrigger(*eId));
0363     if (hcalTTDetId.null())
0364       continue;
0365
0366     int aieta = abs(hcalTTDetId.ieta());
0367
0368     // The absence of TT channels in the HcalTPChannelParameters
0369     // is intepreted as to not use the new filter
0370     int weight = -1.0;
0371     auto tpParam = conditions.getHcalTPChannelParameter(hcalTTDetId, false);
0372     if (tpParam)
0373       weight = tpParam->getauxi1();
0374
0375     if (aieta <= lastHBRing) {
0376       foundHB = true;
0377       if (weight != -1.0)
0378         newHBtp = true;
0379     } else if (aieta > lastHBRing and aieta < lastHERing) {
0380       foundHE = true;
0381       if (weight != -1.0)
0382         newHEtp = true;
0383     }
0384   }
0385
0386   for (const auto& id : metadata->getAllChannels()) {
0387     if (not(id.det() == DetId::Hcal and topo_->valid(id)))
0388       continue;
0389
0390     HcalDetId cell(id);
0391     HcalSubdetector subdet = cell.subdet();
0392
0393     if (subdet != HcalBarrel and subdet != HcalEndcap and subdet != HcalForward)
0394       continue;
0395
0396     const HcalQIECoder* channelCoder = conditions.getHcalCoder(cell);
0397     const HcalQIEShape* shape = conditions.getHcalShape(cell);
0398     HcalCoderDb coder(*channelCoder, *shape);
0399     const HcalLutMetadatum* meta = metadata->getValues(cell);
0400
0401     unsigned int mipMax = 0;
0402     unsigned int mipMin = 0;
0403     unsigned int bit12_energy =
0404         0;  // defaults for energy requirement for bits 12-15 are high / low to avoid FG bit 0-4 being set when not intended
0405     unsigned int bit13_energy = 0;
0406     unsigned int bit14_energy = 999;
0407     unsigned int bit15_energy = 999;
0408
0409     bool is2018OrLater = topo_->triggerMode() >= HcalTopologyMode::TriggerMode_2018 or
0410                          topo_->triggerMode() == HcalTopologyMode::TriggerMode_2018legacy;
0411     if (is2018OrLater or topo_->dddConstants()->isPlan1(cell)) {
0412       const HcalTPChannelParameter* channelParameters = conditions.getHcalTPChannelParameter(cell);
0413       mipMax = channelParameters->getFGBitInfo() >> 16;
0414       mipMin = channelParameters->getFGBitInfo() & 0xFFFF;
0415       bit12_energy = 16;  // depths 1,2 max energy
0416       bit13_energy = 80;  // depths 3+ min energy
0417       bit14_energy = 64;  // prompt min energy
0418       bit15_energy = 64;  // delayed min energy
0419     }
0420
0421     int lutId = getLUTId(cell);
0422     Lut& lut = inputLUT_[lutId];
0423     float ped = 0;
0424     float gain = 0;
0425     uint32_t status = 0;
0426
0427     if (LUTGenerationMode_) {
0428       const HcalCalibrations& calibrations = conditions.getHcalCalibrations(cell);
0429       for (auto capId : {0, 1, 2, 3}) {
0430         ped += calibrations.effpedestal(capId);
0431         gain += calibrations.LUTrespcorrgain(capId);
0432       }
0433       ped /= 4.0;
0434       gain /= 4.0;
0435
0436       //Get Channel Quality
0437       const HcalChannelStatus* channelStatus = conditions.getHcalChannelStatus(cell);
0438       status = channelStatus->getValue();
0439
0440     } else {
0441       const HcalL1TriggerObject* myL1TObj = conditions.getHcalL1TriggerObject(cell);
0442       ped = myL1TObj->getPedestal();
0443       gain = myL1TObj->getRespGain();
0444       status = myL1TObj->getFlag();
0445     }  // LUTGenerationMode_
0446
0447     ped_[lutId] = ped;
0448     gain_[lutId] = gain;
0449     bool isMasked = ((status & bitToMask_) > 0);
0450     float rcalib = meta->getRCalib();
0451
0452     auto adc2fC = [channelCoder, shape](unsigned int adc) {
0453       float fC = 0;
0454       for (auto capId : {0, 1, 2, 3})
0455         fC += channelCoder->charge(*shape, adc, capId);
0456       return fC / 4;
0457     };
0458
0459     int qieType = conditions.getHcalQIEType(cell)->getValue();
0460
0461     const size_t SIZE = qieType == QIE8 ? INPUT_LUT_SIZE : UPGRADE_LUT_SIZE;
0462     const int MASK = qieType == QIE8 ? QIE8_LUT_BITMASK : qieType == QIE10 ? QIE10_LUT_BITMASK : QIE11_LUT_BITMASK;
0463     double linearLSB = linearLSB_QIE8_;
0464     if (qieType == QIE11 and cell.ietaAbs() == topo_->lastHBRing())
0465       linearLSB = linearLSB_QIE11Overlap_;
0466     else if (qieType == QIE11)
0467       linearLSB = linearLSB_QIE11_;
0468
0469     lut.resize(SIZE, 0);
0470
0471     // Input LUT for HB/HE/HF
0472     if (subdet == HcalBarrel || subdet == HcalEndcap) {
0473       int granularity = meta->getLutGranularity();
0474
0475       double correctionPhaseNS = conditions.getHcalRecoParam(cell)->correctionPhaseNS();
0476
0477       if (qieType == QIE11) {
0478         if (overrideDBweightsAndFilterHB_ and cell.ietaAbs() <= lastHBRing)
0479           correctionPhaseNS = containPhaseNSHB_;
0480         else if (overrideDBweightsAndFilterHE_ and cell.ietaAbs() > lastHBRing)
0481           correctionPhaseNS = containPhaseNSHE_;
0482       }
0483       for (unsigned int adc = 0; adc < SIZE; ++adc) {
0484         if (isMasked)
0485           lut[adc] = 0;
0486         else {
0487           double nonlinearityCorrection = 1.0;
0488           double containmentCorrection = 1.0;
0489           // SiPM nonlinearity was not corrected in 2017
0490           // and containment corrections  were not
0491           // ET-dependent prior to 2018
0492           if (is2018OrLater) {
0493             double containmentCorrection1TS = pulseCorr_->correction(cell, 1, correctionPhaseNS, adc2fC(adc));
0494             // Use the 1-TS containment correction to estimate the charge of the pulse
0495             // from the individual samples
0496             double correctedCharge = containmentCorrection1TS * adc2fC(adc);
0497             double containmentCorrection2TSCorrected =
0498                 pulseCorr_->correction(cell, 2, correctionPhaseNS, correctedCharge);
0499             if (qieType == QIE11) {
0500               // When contain1TS_ is set, it should still only apply for QIE11-related things
0501               if ((((contain1TSHB_ and overrideDBweightsAndFilterHB_) or newHBtp) and cell.ietaAbs() <= lastHBRing) or
0502                   (((contain1TSHE_ and overrideDBweightsAndFilterHE_) or newHEtp) and cell.ietaAbs() > lastHBRing)) {
0503                 containmentCorrection = containmentCorrection1TS;
0504               } else {
0505                 containmentCorrection = containmentCorrection2TSCorrected;
0506               }
0507
0508               const HcalSiPMParameter& siPMParameter(*conditions.getHcalSiPMParameter(cell));
0509               HcalSiPMnonlinearity corr(
0510                   conditions.getHcalSiPMCharacteristics()->getNonLinearities(siPMParameter.getType()));
0511               const double fcByPE = siPMParameter.getFCByPE();
0512               const double effectivePixelsFired = correctedCharge / fcByPE;
0513               nonlinearityCorrection = corr.getRecoCorrectionFactor(effectivePixelsFired);
0514             } else {
0515               containmentCorrection = containmentCorrection2TSCorrected;
0516             }
0517           }
0518           if (allLinear_)
0519             lut[adc] = (LutElement)std::min(
0520                 std::max(0,
0521                          int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection * containmentCorrection /
0522                              linearLSB / cosh_ieta(cell.ietaAbs(), cell.depth(), HcalEndcap))),
0523                 MASK);
0524           else
0525             lut[adc] = (LutElement)std::min(std::max(0,
0526                                                      int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection *
0527                                                          containmentCorrection / nominalgain_ / granularity)),
0528                                             MASK);
0529
0530           unsigned int linearizedADC =
0531               lut[adc];  // used for bits 12, 13, 14, 15 for Group 0 LUT for LLP time and depth bits that rely on linearized energies
0532
0533           if (qieType == QIE11) {
0534             if ((linearizedADC < bit12_energy and cell.depth() <= 2) or (cell.depth() >= 3))
0535               lut[adc] |= 1 << 12;
0536             if (linearizedADC >= bit13_energy and cell.depth() >= 3)
0537               lut[adc] |= 1 << 13;
0538             if (linearizedADC >= bit14_energy)
0539               lut[adc] |= 1 << 14;
0540             if (linearizedADC >= bit15_energy)
0541               lut[adc] |= 1 << 15;
0542             if (adc >= mipMin and adc < mipMax)
0543               lut[adc] |= QIE11_LUT_MSB0;
0544             else if (adc >= mipMax)
0545               lut[adc] |= QIE11_LUT_MSB1;
0546           }
0547
0548           //Zeroing the 4th depth in the trigger towers where |ieta| = 16 to match the behavior in the uHTR firmware in Run3, where the 4th depth is not included in the sum over depths when constructing the TP energy for this tower.
0549           if (abs(cell.ieta()) == 16 && cell.depth() == 4 &&
0550               topo_->triggerMode() >= HcalTopologyMode::TriggerMode_2021) {
0551             lut[adc] = 0;
0552           }
0553         }
0554       }
0555     } else if (subdet == HcalForward) {
0556       for (unsigned int adc = 0; adc < SIZE; ++adc) {
0557         if (isMasked)
0558           lut[adc] = 0;
0559         else {
0560           lut[adc] =
0561               std::min(std::max(0, int((adc2fC(adc) - ped) * gain * rcalib / lsb_ / cosh_ieta_[cell.ietaAbs()])), MASK);
0562           if (adc > FG_HF_thresholds_[0])
0563             lut[adc] |= QIE10_LUT_MSB0;
0564           if (adc > FG_HF_thresholds_[1])
0565             lut[adc] |= QIE10_LUT_MSB1;
0566         }
0567       }
0568     }
0569   }
0570 }
0571
0572 void HcaluLUTTPGCoder::adc2Linear(const HBHEDataFrame& df, IntegerCaloSamples& ics) const {
0573   int lutId = getLUTId(df.id());
0574   const Lut& lut = inputLUT_.at(lutId);
0575   for (int i = 0; i < df.size(); i++) {
0576     ics[i] = (lut.at(df[i].adc()) & QIE8_LUT_BITMASK);
0577   }
0578 }
0579
0580 void HcaluLUTTPGCoder::adc2Linear(const HFDataFrame& df, IntegerCaloSamples& ics) const {
0581   int lutId = getLUTId(df.id());
0582   const Lut& lut = inputLUT_.at(lutId);
0583   for (int i = 0; i < df.size(); i++) {
0584     ics[i] = (lut.at(df[i].adc()) & QIE8_LUT_BITMASK);
0585   }
0586 }
0587
0588 void HcaluLUTTPGCoder::adc2Linear(const QIE10DataFrame& df, IntegerCaloSamples& ics) const {
0589   int lutId = getLUTId(HcalDetId(df.id()));
0590   const Lut& lut = inputLUT_.at(lutId);
0591   for (int i = 0; i < df.samples(); i++) {
0592     ics[i] = (lut.at(df[i].adc()) & QIE10_LUT_BITMASK);
0593   }
0594 }
0595
0596 void HcaluLUTTPGCoder::adc2Linear(const QIE11DataFrame& df, IntegerCaloSamples& ics) const {
0597   int lutId = getLUTId(HcalDetId(df.id()));
0598   const Lut& lut = inputLUT_.at(lutId);
0599   for (int i = 0; i < df.samples(); i++) {
0600     ics[i] = (lut.at(df[i].adc()) & QIE11_LUT_BITMASK);
0601   }
0602 }
0603
0604 unsigned short HcaluLUTTPGCoder::adc2Linear(HcalQIESample sample, HcalDetId id) const {
0605   int lutId = getLUTId(id);
0606   return ((inputLUT_.at(lutId)).at(sample.adc()) & QIE8_LUT_BITMASK);
0607 }
0608
0609 std::vector<unsigned short> HcaluLUTTPGCoder::group0FGbits(const QIE11DataFrame& df) const {
0610   int lutId = getLUTId(HcalDetId(df.id()));
0611   const Lut& lut = inputLUT_.at(lutId);
0612   std::vector<unsigned short> group0LLPbits;
0613   group0LLPbits.reserve(df.samples());
0614   for (int i = 0; i < df.samples(); i++) {
0615     group0LLPbits.push_back((lut.at(df[i].adc()) >> 12) &
0616                             0xF);  // four bits (12-15) of LUT used to set 6 finegrain bits from uHTR
0617   }
0618   return group0LLPbits;
0619 }
0620
0621 float HcaluLUTTPGCoder::getLUTPedestal(HcalDetId id) const {
0622   int lutId = getLUTId(id);
0623   return ped_.at(lutId);
0624 }
0625
0626 float HcaluLUTTPGCoder::getLUTGain(HcalDetId id) const {
0627   int lutId = getLUTId(id);
0628   return gain_.at(lutId);
0629 }
0630
0631 std::vector<unsigned short> HcaluLUTTPGCoder::getLinearizationLUT(HcalDetId id) const {
0632   int lutId = getLUTId(id);
0633   return inputLUT_.at(lutId);
0634 }
0635
0636 void HcaluLUTTPGCoder::lookupMSB(const HBHEDataFrame& df, std::vector<bool>& msb) const {
0637   msb.resize(df.size());
0638   for (int i = 0; i < df.size(); ++i)
0639     msb[i] = getMSB(df.id(), df.sample(i).adc());
0640 }
0641
0642 bool HcaluLUTTPGCoder::getMSB(const HcalDetId& id, int adc) const {
0643   int lutId = getLUTId(id);
0644   const Lut& lut = inputLUT_.at(lutId);
0645   return (lut.at(adc) & QIE8_LUT_MSB);
0646 }
0647
0648 void HcaluLUTTPGCoder::lookupMSB(const QIE10DataFrame& df, std::vector<std::bitset<2>>& msb) const {
0649   msb.resize(df.samples());
0650   int lutId = getLUTId(HcalDetId(df.id()));
0651   const Lut& lut = inputLUT_.at(lutId);
0652   for (int i = 0; i < df.samples(); ++i) {
0653     msb[i][0] = lut.at(df[i].adc()) & QIE10_LUT_MSB0;
0654     msb[i][1] = lut.at(df[i].adc()) & QIE10_LUT_MSB1;
0655   }
0656 }
0657
0658 void HcaluLUTTPGCoder::lookupMSB(const QIE11DataFrame& df, std::vector<std::bitset<2>>& msb) const {
0659   int lutId = getLUTId(HcalDetId(df.id()));
0660   const Lut& lut = inputLUT_.at(lutId);
0661   for (int i = 0; i < df.samples(); ++i) {
0662     msb[i][0] = lut.at(df[i].adc()) & QIE11_LUT_MSB0;
0663     msb[i][1] = lut.at(df[i].adc()) & QIE11_LUT_MSB1;
0664   }
0665 }