HcalTPGAlgos/src/HcaluLUTTPGCoder.cc

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