File indexing completed on 2024-04-06 11:57:33
0001 #include "CalibCalorimetry/CaloTPG/interface/CaloTPGTranscoderULUT.h"
0002 #include "DataFormats/HcalDetId/interface/HcalGenericDetId.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/Framework/interface/ESTransientHandle.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008 #include <iostream>
0009 #include <fstream>
0010 #include <cmath>
0011
0012
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "CondFormats/DataRecord/interface/HcalLutMetadataRcd.h"
0015 #include "CalibCalorimetry/HcalTPGAlgos/interface/HcaluLUTTPGCoder.h"
0016
0017 using namespace std;
0018
0019 CaloTPGTranscoderULUT::CaloTPGTranscoderULUT(const std::string& compressionFile, const std::string& decompressionFile)
0020 : theTopology(nullptr),
0021 nominal_gain_(0.),
0022 lsb_factor_(0.),
0023 rct_factor_(1.),
0024 nct_factor_(1.),
0025 lin8_factor_(1.),
0026 lin11_factor_(1.),
0027 compressionFile_(compressionFile),
0028 decompressionFile_(decompressionFile) {}
0029
0030 CaloTPGTranscoderULUT::~CaloTPGTranscoderULUT() {}
0031
0032 void CaloTPGTranscoderULUT::loadHCALCompress(HcalLutMetadata const& lutMetadata,
0033 HcalTrigTowerGeometry const& theTrigTowerGeometry) {
0034 if (!theTopology) {
0035 throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
0036 }
0037
0038 std::array<unsigned int, OUTPUT_LUT_SIZE> analyticalLUT;
0039 std::array<unsigned int, OUTPUT_LUT_SIZE> linearQIE8LUT;
0040 std::array<unsigned int, OUTPUT_LUT_SIZE> linearQIE11LUT;
0041 std::array<unsigned int, OUTPUT_LUT_SIZE> linearRctLUT;
0042 std::array<unsigned int, OUTPUT_LUT_SIZE> linearNctLUT;
0043
0044
0045 for (unsigned int i = 0; i < OUTPUT_LUT_SIZE; i++) {
0046 analyticalLUT[i] = min(static_cast<unsigned int>(sqrt(14.94 * log(1. + i / 14.94) * i) + 0.5), TPGMAX - 1);
0047 linearQIE8LUT[i] = min(static_cast<unsigned int>((i + .5) / lin8_factor_), TPGMAX - 1);
0048 linearQIE11LUT[i] = min(static_cast<unsigned int>((i + .5) / lin11_factor_), TPGMAX - 1);
0049 linearRctLUT[i] = min(static_cast<unsigned int>((i + .5) / rct_factor_), TPGMAX - 1);
0050 linearNctLUT[i] = min(static_cast<unsigned int>((i + .5) / nct_factor_), TPGMAX - 1);
0051 }
0052
0053 std::vector<DetId> allChannels = lutMetadata.getAllChannels();
0054
0055 for (std::vector<DetId>::iterator i = allChannels.begin(); i != allChannels.end(); ++i) {
0056 if (not HcalGenericDetId(*i).isHcalTrigTowerDetId()) {
0057 if ((not HcalGenericDetId(*i).isHcalDetId()) and (not HcalGenericDetId(*i).isHcalZDCDetId()) and
0058 (not HcalGenericDetId(*i).isHcalCastorDetId()))
0059 edm::LogWarning("CaloTPGTranscoderULUT") << "Encountered invalid HcalDetId " << HcalGenericDetId(*i);
0060 continue;
0061 }
0062
0063 HcalTrigTowerDetId id(*i);
0064 if (!theTopology->validHT(id))
0065 continue;
0066
0067 unsigned int index = getOutputLUTId(id);
0068
0069 const HcalLutMetadatum* meta = lutMetadata.getValues(id);
0070 unsigned int threshold = meta->getOutputLutThreshold();
0071
0072 int ieta = id.ieta();
0073 int version = id.version();
0074 bool isHBHE = (abs(ieta) < theTrigTowerGeometry.firstHFTower(version));
0075
0076 unsigned int lutsize = getOutputLUTSize(id);
0077 outputLUT_[index].resize(lutsize);
0078
0079 for (unsigned int i = 0; i < threshold; ++i)
0080 outputLUT_[index][i] = 0;
0081
0082 if (isHBHE) {
0083 for (unsigned int i = threshold; i < lutsize; ++i)
0084 if (allLinear_) {
0085 outputLUT_[index][i] = isOnlyQIE11(id) ? linearQIE11LUT[i] : linearQIE8LUT[i];
0086
0087 if (abs(ieta) > 20 && isOnlyQIE11(id) && linearQIE11LUT[i] >= (TPGMAX - 2) / 2.) {
0088 outputLUT_[index][i] = TPGMAX - 1;
0089 }
0090 } else {
0091 outputLUT_[index][i] = analyticalLUT[i];
0092 }
0093 } else {
0094 for (unsigned int i = threshold; i < lutsize; ++i)
0095 outputLUT_[index][i] = version == 1 ? linearNctLUT[i] : linearRctLUT[i];
0096 }
0097
0098 double eta_low = 0., eta_high = 0.;
0099 theTrigTowerGeometry.towerEtaBounds(ieta, version, eta_low, eta_high);
0100 double cosh_ieta = fabs(cosh((eta_low + eta_high) / 2.));
0101 double granularity = meta->getLutGranularity();
0102
0103 if (isHBHE) {
0104 if (allLinear_) {
0105 LUT tpg = outputLUT_[index][0];
0106 hcaluncomp_[index][tpg] = 0;
0107 for (unsigned int i = 0; i < lutsize; ++i) {
0108 if (outputLUT_[index][i] != tpg) {
0109 tpg = outputLUT_[index][i];
0110 hcaluncomp_[index][tpg] = lsb_factor_ * i / (isOnlyQIE11(id) ? lin11_factor_ : lin8_factor_);
0111
0112 if (abs(ieta) > 20 && isOnlyQIE11(id) && linearQIE11LUT[i] >= (TPGMAX - 2) / 2.) {
0113 hcaluncomp_[index][tpg] = (TPGMAX - 1) / 2.;
0114 }
0115 }
0116 }
0117 } else {
0118 double factor = nominal_gain_ / cosh_ieta * granularity;
0119 LUT tpg = outputLUT_[index][0];
0120 int low = 0;
0121 for (unsigned int i = 0; i < lutsize; ++i) {
0122 if (outputLUT_[index][i] != tpg) {
0123 unsigned int mid = (low + i) / 2;
0124 hcaluncomp_[index][tpg] = (tpg == 0 ? low : factor * mid);
0125 low = i;
0126 tpg = outputLUT_[index][i];
0127 }
0128 }
0129 hcaluncomp_[index][tpg] = factor * low;
0130 }
0131 } else {
0132 LUT tpg = outputLUT_[index][0];
0133 hcaluncomp_[index][tpg] = 0;
0134 for (unsigned int i = 0; i < lutsize; ++i) {
0135 if (outputLUT_[index][i] != tpg) {
0136 tpg = outputLUT_[index][i];
0137 hcaluncomp_[index][tpg] = lsb_factor_ * i / (version == 1 ? nct_factor_ : rct_factor_);
0138 }
0139 }
0140 }
0141 }
0142 }
0143
0144 HcalTriggerPrimitiveSample CaloTPGTranscoderULUT::hcalCompress(const HcalTrigTowerDetId& id,
0145 unsigned int sample,
0146 int fineGrain) const {
0147 unsigned int itower = getOutputLUTId(id);
0148
0149 if (sample >= getOutputLUTSize(id))
0150 throw cms::Exception("Out of Range") << "LUT has " << getOutputLUTSize(id) << " entries for " << id << " but "
0151 << sample << " was requested.";
0152
0153 if (itower >= outputLUT_.size())
0154 throw cms::Exception("Out of Range") << "No decompression LUT found for " << id;
0155
0156 return HcalTriggerPrimitiveSample(outputLUT_[itower][sample], fineGrain);
0157 }
0158
0159 double CaloTPGTranscoderULUT::hcaletValue(const int& ieta,
0160 const int& iphi,
0161 const int& version,
0162 const int& compET) const {
0163 double etvalue = 0.;
0164 int itower = getOutputLUTId(ieta, iphi, version);
0165 if (itower < 0) {
0166 edm::LogError("CaloTPGTranscoderULUT") << "No decompression LUT found for ieta, iphi = " << ieta << ", " << iphi;
0167 } else if (compET < 0 || compET >= (int)TPGMAX) {
0168 edm::LogError("CaloTPGTranscoderULUT")
0169 << "Compressed value out of range: eta, phi, cET = " << ieta << ", " << iphi << ", " << compET;
0170 } else {
0171 etvalue = hcaluncomp_[itower][compET];
0172 }
0173 return (etvalue);
0174 }
0175
0176 double CaloTPGTranscoderULUT::hcaletValue(const HcalTrigTowerDetId& hid, const HcalTriggerPrimitiveSample& hc) const {
0177 int compET = hc.compressedEt();
0178 int itower = getOutputLUTId(hid);
0179 double etvalue = hcaluncomp_[itower][compET];
0180 return etvalue;
0181 }
0182
0183 EcalTriggerPrimitiveSample CaloTPGTranscoderULUT::ecalCompress(const EcalTrigTowerDetId& id,
0184 unsigned int sample,
0185 bool fineGrain) const {
0186 throw cms::Exception("Not Implemented") << "CaloTPGTranscoderULUT::ecalCompress";
0187 }
0188
0189 void CaloTPGTranscoderULUT::rctEGammaUncompress(const HcalTrigTowerDetId& hid,
0190 const HcalTriggerPrimitiveSample& hc,
0191 const EcalTrigTowerDetId& eid,
0192 const EcalTriggerPrimitiveSample& ec,
0193 unsigned int& et,
0194 bool& egVecto,
0195 bool& activity) const {
0196 throw cms::Exception("Not Implemented") << "CaloTPGTranscoderULUT::rctEGammaUncompress";
0197 }
0198 void CaloTPGTranscoderULUT::rctJetUncompress(const HcalTrigTowerDetId& hid,
0199 const HcalTriggerPrimitiveSample& hc,
0200 const EcalTrigTowerDetId& eid,
0201 const EcalTriggerPrimitiveSample& ec,
0202 unsigned int& et) const {
0203 throw cms::Exception("Not Implemented") << "CaloTPGTranscoderULUT::rctJetUncompress";
0204 }
0205
0206 bool CaloTPGTranscoderULUT::HTvalid(const int ieta, const int iphiin, const int version) const {
0207 HcalTrigTowerDetId id(ieta, iphiin);
0208 id.setVersion(version);
0209 if (!theTopology) {
0210 throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
0211 }
0212 return theTopology->validHT(id);
0213 }
0214
0215 int CaloTPGTranscoderULUT::getOutputLUTId(const HcalTrigTowerDetId& id) const {
0216 if (!theTopology) {
0217 throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
0218 }
0219 return theTopology->detId2denseIdHT(id);
0220 }
0221
0222 int CaloTPGTranscoderULUT::getOutputLUTId(const int ieta, const int iphiin, const int version) const {
0223 if (!theTopology) {
0224 throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
0225 }
0226 HcalTrigTowerDetId id(ieta, iphiin);
0227 id.setVersion(version);
0228 return theTopology->detId2denseIdHT(id);
0229 }
0230
0231 unsigned int CaloTPGTranscoderULUT::getOutputLUTSize(const HcalTrigTowerDetId& id) const {
0232 if (!theTopology)
0233 throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
0234
0235 switch (theTopology->triggerMode()) {
0236 case HcalTopologyMode::TriggerMode_2009:
0237 case HcalTopologyMode::TriggerMode_2016:
0238 return QIE8_OUTPUT_LUT_SIZE;
0239 case HcalTopologyMode::TriggerMode_2017:
0240 if (id.ietaAbs() <= theTopology->lastHERing())
0241 return QIE8_OUTPUT_LUT_SIZE;
0242 else
0243 return QIE10_OUTPUT_LUT_SIZE;
0244 case HcalTopologyMode::TriggerMode_2017plan1:
0245 if (plan1_towers_.find(id) != plan1_towers_.end())
0246 return QIE11_OUTPUT_LUT_SIZE;
0247 else if (id.ietaAbs() <= theTopology->lastHERing())
0248 return QIE8_OUTPUT_LUT_SIZE;
0249 else
0250 return QIE10_OUTPUT_LUT_SIZE;
0251 case HcalTopologyMode::TriggerMode_2018legacy:
0252 case HcalTopologyMode::TriggerMode_2018:
0253 if (id.ietaAbs() <= theTopology->lastHBRing())
0254 return QIE8_OUTPUT_LUT_SIZE;
0255 else if (id.ietaAbs() <= theTopology->lastHERing())
0256 return QIE11_OUTPUT_LUT_SIZE;
0257 else
0258 return QIE10_OUTPUT_LUT_SIZE;
0259 case HcalTopologyMode::TriggerMode_2021:
0260 if (id.ietaAbs() <= theTopology->lastHBHERing())
0261 return QIE11_OUTPUT_LUT_SIZE;
0262 else
0263 return QIE10_OUTPUT_LUT_SIZE;
0264 default:
0265 throw cms::Exception("CaloTPGTranscoderULUT") << "Unknown trigger mode used by the topology!";
0266 }
0267 }
0268
0269 bool CaloTPGTranscoderULUT::isOnlyQIE11(const HcalTrigTowerDetId& id) const {
0270 if (!theTopology)
0271 throw cms::Exception("CaloTPGTranscoderULUT") << "Topology not set! Use CaloTPGTranscoderULUT::setup(...) first!";
0272
0273 switch (theTopology->triggerMode()) {
0274 case HcalTopologyMode::TriggerMode_2017plan1:
0275 if (plan1_towers_.find(id) != plan1_towers_.end() and id.ietaAbs() != theTopology->lastHBRing())
0276 return true;
0277 return false;
0278 case HcalTopologyMode::TriggerMode_2018legacy:
0279 case HcalTopologyMode::TriggerMode_2018:
0280 if (id.ietaAbs() <= theTopology->lastHBRing())
0281 return false;
0282 return true;
0283 case HcalTopologyMode::TriggerMode_2021:
0284 return true;
0285 default:
0286 throw cms::Exception("CaloTPGTranscoderULUT") << "Unknown trigger mode used by the topology!";
0287 }
0288 }
0289
0290 const std::vector<unsigned int> CaloTPGTranscoderULUT::getCompressionLUT(const HcalTrigTowerDetId& id) const {
0291 int itower = getOutputLUTId(id);
0292 auto lut = outputLUT_[itower];
0293 std::vector<unsigned int> result(lut.begin(), lut.end());
0294 return result;
0295 }
0296
0297 void CaloTPGTranscoderULUT::setup(HcalLutMetadata const& lutMetadata,
0298 HcalTrigTowerGeometry const& theTrigTowerGeometry,
0299 int nctScaleShift,
0300 int rctScaleShift,
0301 double lsbQIE8,
0302 double lsbQIE11,
0303 bool allLinear) {
0304 theTopology = lutMetadata.topo();
0305 nominal_gain_ = lutMetadata.getNominalGain();
0306 lsb_factor_ = lutMetadata.getRctLsb();
0307
0308 allLinear_ = allLinear;
0309
0310 rct_factor_ = lsb_factor_ / (HcaluLUTTPGCoder::lsb_ * (1 << rctScaleShift));
0311 nct_factor_ = lsb_factor_ / (HcaluLUTTPGCoder::lsb_ * (1 << nctScaleShift));
0312 lin8_factor_ = lsb_factor_ / lsbQIE8;
0313 lin11_factor_ = lsb_factor_ / lsbQIE11;
0314
0315 outputLUT_.resize(theTopology->getHTSize());
0316 hcaluncomp_.resize(theTopology->getHTSize());
0317
0318 plan1_towers_.clear();
0319 for (const auto& id : lutMetadata.getAllChannels()) {
0320 if (not(id.det() == DetId::Hcal and theTopology->valid(id)))
0321 continue;
0322 HcalDetId cell(id);
0323 if (not theTopology->dddConstants()->isPlan1(cell))
0324 continue;
0325 for (const auto& tower : theTrigTowerGeometry.towerIds(cell))
0326 plan1_towers_.emplace(tower);
0327 }
0328
0329 if (compressionFile_.empty() && decompressionFile_.empty()) {
0330 loadHCALCompress(lutMetadata, theTrigTowerGeometry);
0331 } else {
0332 throw cms::Exception("Not Implemented") << "setup of CaloTPGTranscoderULUT from text files";
0333 }
0334 }
