Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-03 01:05:45 

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 bit12_energy =
0402         0;  // defaults for energy requirement for bits 12-15 are high / low to avoid FG bit 0-4 being set when not intended
0403     unsigned int bit13_energy = 0;
0404     unsigned int bit14_energy = 999;
0405     unsigned int bit15_energy = 999;
0406 
0407     bool is2018OrLater = topo_->triggerMode() >= HcalTopologyMode::TriggerMode_2018 or
0408                          topo_->triggerMode() == HcalTopologyMode::TriggerMode_2018legacy;
0409     if (is2018OrLater or topo_->dddConstants()->isPlan1(cell)) {
0410       bit12_energy = 16;  // depths 1,2 max energy
0411       bit13_energy = 80;  // depths 3+ min energy
0412       bit14_energy = 64;  // prompt min energy
0413       bit15_energy = 64;  // delayed min energy
0414     }
0415 
0416     int lutId = getLUTId(cell);
0417     Lut& lut = inputLUT_[lutId];
0418     float ped = 0;
0419     float gain = 0;
0420     uint32_t status = 0;
0421 
0422     if (LUTGenerationMode_) {
0423       const HcalCalibrations& calibrations = conditions.getHcalCalibrations(cell);
0424       for (auto capId : {0, 1, 2, 3}) {
0425         ped += calibrations.effpedestal(capId);
0426         gain += calibrations.LUTrespcorrgain(capId);
0427       }
0428       ped /= 4.0;
0429       gain /= 4.0;
0430 
0431       //Get Channel Quality
0432       const HcalChannelStatus* channelStatus = conditions.getHcalChannelStatus(cell);
0433       status = channelStatus->getValue();
0434 
0435     } else {
0436       const HcalL1TriggerObject* myL1TObj = conditions.getHcalL1TriggerObject(cell);
0437       ped = myL1TObj->getPedestal();
0438       gain = myL1TObj->getRespGain();
0439       status = myL1TObj->getFlag();
0440     }  // LUTGenerationMode_
0441 
0442     ped_[lutId] = ped;
0443     gain_[lutId] = gain;
0444     bool isMasked = ((status & bitToMask_) > 0);
0445     float rcalib = meta->getRCalib();
0446 
0447     auto adc2fC = [channelCoder, shape](unsigned int adc) {
0448       float fC = 0;
0449       for (auto capId : {0, 1, 2, 3})
0450         fC += channelCoder->charge(*shape, adc, capId);
0451       return fC / 4;
0452     };
0453 
0454     int qieType = conditions.getHcalQIEType(cell)->getValue();
0455 
0456     const size_t SIZE = qieType == QIE8 ? INPUT_LUT_SIZE : UPGRADE_LUT_SIZE;
0457     const int MASK = qieType == QIE8 ? QIE8_LUT_BITMASK : qieType == QIE10 ? QIE10_LUT_BITMASK : QIE11_LUT_BITMASK;
0458     double linearLSB = linearLSB_QIE8_;
0459     if (qieType == QIE11 and cell.ietaAbs() == topo_->lastHBRing())
0460       linearLSB = linearLSB_QIE11Overlap_;
0461     else if (qieType == QIE11)
0462       linearLSB = linearLSB_QIE11_;
0463 
0464     lut.resize(SIZE, 0);
0465 
0466     // Input LUT for HB/HE/HF
0467     if (subdet == HcalBarrel || subdet == HcalEndcap) {
0468       int granularity = meta->getLutGranularity();
0469 
0470       double correctionPhaseNS = conditions.getHcalRecoParam(cell)->correctionPhaseNS();
0471 
0472       if (qieType == QIE11) {
0473         if (overrideDBweightsAndFilterHB_ and cell.ietaAbs() <= lastHBRing)
0474           correctionPhaseNS = containPhaseNSHB_;
0475         else if (overrideDBweightsAndFilterHE_ and cell.ietaAbs() > lastHBRing)
0476           correctionPhaseNS = containPhaseNSHE_;
0477       }
0478       for (unsigned int adc = 0; adc < SIZE; ++adc) {
0479         if (isMasked)
0480           lut[adc] = 0;
0481         else {
0482           double nonlinearityCorrection = 1.0;
0483           double containmentCorrection = 1.0;
0484           // SiPM nonlinearity was not corrected in 2017
0485           // and containment corrections  were not
0486           // ET-dependent prior to 2018
0487           if (is2018OrLater) {
0488             double containmentCorrection1TS = pulseCorr_->correction(cell, 1, correctionPhaseNS, adc2fC(adc));
0489             // Use the 1-TS containment correction to estimate the charge of the pulse
0490             // from the individual samples
0491             double correctedCharge = containmentCorrection1TS * adc2fC(adc);
0492             double containmentCorrection2TSCorrected =
0493                 pulseCorr_->correction(cell, 2, correctionPhaseNS, correctedCharge);
0494             if (qieType == QIE11) {
0495               // When contain1TS_ is set, it should still only apply for QIE11-related things
0496               if ((((contain1TSHB_ and overrideDBweightsAndFilterHB_) or newHBtp) and cell.ietaAbs() <= lastHBRing) or
0497                   (((contain1TSHE_ and overrideDBweightsAndFilterHE_) or newHEtp) and cell.ietaAbs() > lastHBRing)) {
0498                 containmentCorrection = containmentCorrection1TS;
0499               } else {
0500                 containmentCorrection = containmentCorrection2TSCorrected;
0501               }
0502 
0503               const HcalSiPMParameter& siPMParameter(*conditions.getHcalSiPMParameter(cell));
0504               HcalSiPMnonlinearity corr(
0505                   conditions.getHcalSiPMCharacteristics()->getNonLinearities(siPMParameter.getType()));
0506               const double fcByPE = siPMParameter.getFCByPE();
0507               const double effectivePixelsFired = correctedCharge / fcByPE;
0508               nonlinearityCorrection = corr.getRecoCorrectionFactor(effectivePixelsFired);
0509             } else {
0510               containmentCorrection = containmentCorrection2TSCorrected;
0511             }
0512           }
0513           if (allLinear_)
0514             lut[adc] = (LutElement)std::min(
0515                 std::max(0,
0516                          int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection * containmentCorrection /
0517                              linearLSB / cosh_ieta(cell.ietaAbs(), cell.depth(), HcalEndcap))),
0518                 MASK);
0519           else
0520             lut[adc] = (LutElement)std::min(std::max(0,
0521                                                      int((adc2fC(adc) - ped) * gain * rcalib * nonlinearityCorrection *
0522                                                          containmentCorrection / nominalgain_ / granularity)),
0523                                             MASK);
0524 
0525           unsigned int linearizedADC =
0526               lut[adc];  // used for bits 12, 13, 14, 15 for Group 0 LUT for LLP time and depth bits that rely on linearized energies
0527 
0528           if (qieType == QIE11) {
0529             if (subdet == HcalBarrel) {  // edit since bits 12-15 not supported in HE yet
0530               if ((linearizedADC < bit12_energy and cell.depth() <= 2) or (cell.depth() >= 3))
0531                 lut[adc] |= 1 << 12;
0532               if (linearizedADC >= bit13_energy and cell.depth() >= 3)
0533                 lut[adc] |= 1 << 13;
0534               if (linearizedADC >= bit14_energy)
0535                 lut[adc] |= 1 << 14;
0536               if (linearizedADC >= bit15_energy)
0537                 lut[adc] |= 1 << 15;
0538             }
0539           }
0540 
0541           //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.
0542           if (abs(cell.ieta()) == 16 && cell.depth() == 4 &&
0543               topo_->triggerMode() >= HcalTopologyMode::TriggerMode_2021) {
0544             lut[adc] = 0;
0545           }
0546         }
0547       }
0548     } else if (subdet == HcalForward) {
0549       for (unsigned int adc = 0; adc < SIZE; ++adc) {
0550         if (isMasked)
0551           lut[adc] = 0;
0552         else {
0553           lut[adc] =
0554               std::min(std::max(0, int((adc2fC(adc) - ped) * gain * rcalib / lsb_ / cosh_ieta_[cell.ietaAbs()])), MASK);
0555           if (adc > FG_HF_thresholds_[0])
0556             lut[adc] |= QIE10_LUT_MSB0;
0557           if (adc > FG_HF_thresholds_[1])
0558             lut[adc] |= QIE10_LUT_MSB1;
0559         }
0560       }
0561     }
0562   }
0563 }
0564 
0565 void HcaluLUTTPGCoder::adc2Linear(const HBHEDataFrame& df, IntegerCaloSamples& ics) const {
0566   int lutId = getLUTId(df.id());
0567   const Lut& lut = inputLUT_.at(lutId);
0568   for (int i = 0; i < df.size(); i++) {
0569     ics[i] = (lut.at(df[i].adc()) & QIE8_LUT_BITMASK);
0570   }
0571 }
0572 
0573 void HcaluLUTTPGCoder::adc2Linear(const HFDataFrame& df, IntegerCaloSamples& ics) const {
0574   int lutId = getLUTId(df.id());
0575   const Lut& lut = inputLUT_.at(lutId);
0576   for (int i = 0; i < df.size(); i++) {
0577     ics[i] = (lut.at(df[i].adc()) & QIE8_LUT_BITMASK);
0578   }
0579 }
0580 
0581 void HcaluLUTTPGCoder::adc2Linear(const QIE10DataFrame& df, IntegerCaloSamples& ics) const {
0582   int lutId = getLUTId(HcalDetId(df.id()));
0583   const Lut& lut = inputLUT_.at(lutId);
0584   for (int i = 0; i < df.samples(); i++) {
0585     ics[i] = (lut.at(df[i].adc()) & QIE10_LUT_BITMASK);
0586   }
0587 }
0588 
0589 void HcaluLUTTPGCoder::adc2Linear(const QIE11DataFrame& df, IntegerCaloSamples& ics) const {
0590   int lutId = getLUTId(HcalDetId(df.id()));
0591   const Lut& lut = inputLUT_.at(lutId);
0592   for (int i = 0; i < df.samples(); i++) {
0593     ics[i] = (lut.at(df[i].adc()) & QIE11_LUT_BITMASK);
0594   }
0595 }
0596 
0597 unsigned short HcaluLUTTPGCoder::adc2Linear(HcalQIESample sample, HcalDetId id) const {
0598   int lutId = getLUTId(id);
0599   return ((inputLUT_.at(lutId)).at(sample.adc()) & QIE8_LUT_BITMASK);
0600 }
0601 
0602 std::vector<unsigned short> HcaluLUTTPGCoder::group0FGbits(const QIE11DataFrame& df) const {
0603   int lutId = getLUTId(HcalDetId(df.id()));
0604   const Lut& lut = inputLUT_.at(lutId);
0605   std::vector<unsigned short> group0LLPbits;
0606   group0LLPbits.reserve(df.samples());
0607   for (int i = 0; i < df.samples(); i++) {
0608     group0LLPbits.push_back((lut.at(df[i].adc()) >> 12) &
0609                             0xF);  // four bits (12-15) of LUT used to set 6 finegrain bits from uHTR
0610   }
0611   return group0LLPbits;
0612 }
0613 
0614 float HcaluLUTTPGCoder::getLUTPedestal(HcalDetId id) const {
0615   int lutId = getLUTId(id);
0616   return ped_.at(lutId);
0617 }
0618 
0619 float HcaluLUTTPGCoder::getLUTGain(HcalDetId id) const {
0620   int lutId = getLUTId(id);
0621   return gain_.at(lutId);
0622 }
0623 
0624 std::vector<unsigned short> HcaluLUTTPGCoder::getLinearizationLUT(HcalDetId id) const {
0625   int lutId = getLUTId(id);
0626   return inputLUT_.at(lutId);
0627 }
0628 
0629 void HcaluLUTTPGCoder::lookupMSB(const HBHEDataFrame& df, std::vector<bool>& msb) const {
0630   msb.resize(df.size());
0631   for (int i = 0; i < df.size(); ++i)
0632     msb[i] = getMSB(df.id(), df.sample(i).adc());
0633 }
0634 
0635 bool HcaluLUTTPGCoder::getMSB(const HcalDetId& id, int adc) const {
0636   int lutId = getLUTId(id);
0637   const Lut& lut = inputLUT_.at(lutId);
0638   return (lut.at(adc) & QIE8_LUT_MSB);
0639 }
0640 
0641 void HcaluLUTTPGCoder::lookupMSB(const QIE10DataFrame& df, std::vector<std::bitset<2>>& msb) const {
0642   msb.resize(df.samples());
0643   int lutId = getLUTId(HcalDetId(df.id()));
0644   const Lut& lut = inputLUT_.at(lutId);
0645   for (int i = 0; i < df.samples(); ++i) {
0646     msb[i][0] = lut.at(df[i].adc()) & QIE10_LUT_MSB0;
0647     msb[i][1] = lut.at(df[i].adc()) & QIE10_LUT_MSB1;
0648   }
0649 }
0650 
0651 void HcaluLUTTPGCoder::lookupMSB(const QIE11DataFrame& df, std::vector<std::bitset<2>>& msb) const {
0652   int lutId = getLUTId(HcalDetId(df.id()));
0653   const Lut& lut = inputLUT_.at(lutId);
0654   for (int i = 0; i < df.samples(); ++i) {
0655     msb[i][0] = lut.at(df[i].adc()) & QIE11_LUT_MSB0;
0656     msb[i][1] = lut.at(df[i].adc()) & QIE11_LUT_MSB1;
0657   }
0658 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!