Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-04 02:52:34

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 //#include "FWCore/Framework/interface/Frameworkfwd.h"
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   // Compute compression LUT
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           //Modifying the saturation (127 -> 255) for the 'split cells'.
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             //Modifying the saturation for the 'split cells'
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();  // to be within the range by the class
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 }