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