RegionalCaloTrigger/src/L1RCTLookupTables.cc

0001 #include <iostream>
0002 using std::cout;
0003 using std::endl;
0004
0005 #include <cmath>
0006 #include <fstream>
0007 #include <string>
0008 #include <vector>
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012
0013 #include "L1Trigger/RegionalCaloTrigger/interface/L1RCTLookupTables.h"
0014
0015 #include "CondFormats/L1TObjects/interface/L1CaloEcalScale.h"
0016 #include "CondFormats/L1TObjects/interface/L1CaloEtScale.h"
0017 #include "CondFormats/L1TObjects/interface/L1CaloHcalScale.h"
0018 #include "CondFormats/L1TObjects/interface/L1RCTChannelMask.h"
0019 #include "CondFormats/L1TObjects/interface/L1RCTNoisyChannelMask.h"
0020 #include "CondFormats/L1TObjects/interface/L1RCTParameters.h"
0021
0022 unsigned int L1RCTLookupTables::lookup(unsigned short ecalInput,
0023                                        unsigned short hcalInput,
0024                                        unsigned short fgbit,
0025                                        unsigned short crtNo,
0026                                        unsigned short crdNo,
0027                                        unsigned short twrNo) const {
0028   if (rctParameters_ == nullptr)
0029     throw cms::Exception("L1RCTParameters Invalid") << "L1RCTParameters should be set every event" << rctParameters_;
0030   if (channelMask_ == nullptr)
0031     throw cms::Exception("L1RCTChannelMask Invalid") << "L1RCTChannelMask should be set every event" << channelMask_;
0032   if (noisyChannelMask_ == nullptr)
0033     throw cms::Exception("L1RCTNoisyChannelMask Invalid")
0034         << "L1RCTNoisyChannelMask should be set every event" << noisyChannelMask_;
0035   if (ecalInput > 0xFF)
0036     throw cms::Exception("Invalid Data") << "ECAL compressedET should be less than 0xFF, is " << ecalInput;
0037   if (hcalInput > 0xFF)
0038     throw cms::Exception("Invalid Data") << "HCAL compressedET should be less than 0xFF, is " << hcalInput;
0039   if (fgbit > 1)
0040     throw cms::Exception("Invalid Data") << "ECAL finegrain should be a single bit, is " << fgbit;
0041   short iEta = (short)rctParameters_->calcIEta(crtNo, crdNo, twrNo);
0042   unsigned short iAbsEta = (unsigned short)abs(iEta);
0043   short sign = iEta / iAbsEta;
0044   unsigned short iPhi = rctParameters_->calcIPhi(crtNo, crdNo, twrNo);
0045   unsigned short phiSide = (iPhi / 4) % 2;
0046   if (iAbsEta < 1 || iAbsEta > 28)
0047     throw cms::Exception("Invalid Data") << "1 <= |IEta| <= 28, is " << iAbsEta;
0048
0049   // Pre Input bits
0050   unsigned short ecalAfterMask = 0;
0051   unsigned short hcalAfterMask = 0;
0052
0053   // using channel mask to mask off ecal channels
0054   // Mike: Introducing the hot channel mask
0055   // If the Et is above the threshold then mask it as well
0056
0057   float ecalBeforeMask = convertEcal(ecalInput, iAbsEta, sign);
0058
0059   bool resetECAL = (channelMask_->ecalMask[crtNo][phiSide][iAbsEta - 1]) ||  // channel mask
0060                    (noisyChannelMask_->ecalMask[crtNo][phiSide][iAbsEta - 1] &&
0061                     ecalBeforeMask < noisyChannelMask_->ecalThreshold) ||  // hot mask
0062                    (rctParameters_->eGammaECalScaleFactors()[iAbsEta - 1] == 0. &&
0063                     rctParameters_->jetMETECalScaleFactors()[iAbsEta - 1] == 0.);
0064
0065   if (resetECAL) {
0066     ecalAfterMask = 0;
0067   } else {
0068     ecalAfterMask = ecalInput;
0069   }
0070
0071   float ecal = convertEcal(ecalAfterMask, iAbsEta, sign);
0072
0073   // masking off hcal for channels in channel mask
0074   float hcalBeforeMask = convertHcal(hcalInput, iAbsEta, sign);
0075
0076   bool resetHCAL = channelMask_->hcalMask[crtNo][phiSide][iAbsEta - 1] ||
0077                    (noisyChannelMask_->hcalMask[crtNo][phiSide][iAbsEta - 1] &&
0078                     hcalBeforeMask < noisyChannelMask_->hcalThreshold) ||  // hot mask
0079                    (rctParameters_->eGammaHCalScaleFactors()[iAbsEta - 1] == 0. &&
0080                     rctParameters_->jetMETHCalScaleFactors()[iAbsEta - 1] == 0.);
0081
0082   if (resetHCAL) {
0083     hcalAfterMask = 0;
0084   } else {
0085     hcalAfterMask = hcalInput;
0086   }
0087
0088   float hcal = convertHcal(hcalAfterMask, iAbsEta, sign);
0089
0090   unsigned long etIn7Bits;
0091   unsigned long etIn9Bits;
0092
0093   if ((ecalAfterMask == 0 && hcalAfterMask > 0) &&
0094       ((rctParameters_->noiseVetoHB() && iAbsEta > 0 && iAbsEta < 18) ||
0095        (rctParameters_->noiseVetoHEplus() && iAbsEta > 17 && crtNo > 8) ||
0096        (rctParameters_->noiseVetoHEminus() && iAbsEta > 17 && crtNo < 9))) {
0097     etIn7Bits = 0;
0098     etIn9Bits = 0;
0099   } else {
0100     etIn7Bits = eGammaETCode(ecal, hcal, iAbsEta);
0101     etIn9Bits = jetMETETCode(ecal, hcal, iAbsEta);
0102   }
0103   // Saturated input towers cause tower ET pegging at the highest value
0104   if ((ecalAfterMask == 0xFF && rctParameters_->eGammaECalScaleFactors()[iAbsEta - 1] != 0.) ||
0105       (hcalAfterMask == 0xFF && rctParameters_->eGammaHCalScaleFactors()[iAbsEta - 1] != 0.)) {
0106     etIn7Bits = 0x7F;  // egamma path
0107   }
0108   if ((ecalAfterMask == 0xFF && rctParameters_->jetMETECalScaleFactors()[iAbsEta - 1] != 0.) ||
0109       (hcalAfterMask == 0xFF && rctParameters_->jetMETHCalScaleFactors()[iAbsEta - 1] != 0.)) {
0110     etIn9Bits = 0x1FF;  // sums path
0111   }
0112
0113   unsigned long shiftEtIn9Bits = etIn9Bits << 8;
0114   unsigned long shiftHE_FGBit = hOeFGVetoBit(ecal, hcal, fgbit) << 7;
0115   unsigned long shiftActivityBit = 0;
0116   if (rctParameters_->jetMETECalScaleFactors()[iAbsEta - 1] == 0. &&
0117       rctParameters_->jetMETHCalScaleFactors()[iAbsEta - 1] == 0.) {
0118     // do nothing, it's already zero
0119   } else if (rctParameters_->jetMETECalScaleFactors()[iAbsEta - 1] == 0.) {
0120     shiftActivityBit = activityBit(0., hcal) << 17;
0121   } else if (rctParameters_->jetMETHCalScaleFactors()[iAbsEta - 1] == 0.) {
0122     shiftActivityBit = activityBit(ecal, 0.) << 17;
0123   } else {
0124     shiftActivityBit = activityBit(ecal, hcal) << 17;
0125   }
0126   unsigned long output = etIn7Bits + shiftHE_FGBit + shiftEtIn9Bits + shiftActivityBit;
0127   return output;
0128 }
0129
0130 unsigned int L1RCTLookupTables::lookup(unsigned short hfInput,
0131                                        unsigned short crtNo,
0132                                        unsigned short crdNo,
0133                                        unsigned short twrNo) const {
0134   if (rctParameters_ == nullptr)
0135     throw cms::Exception("L1RCTParameters Invalid") << "L1RCTParameters should be set every event" << rctParameters_;
0136   if (channelMask_ == nullptr)
0137     throw cms::Exception("L1RCTChannelMask Invalid") << "L1RCTChannelMask should be set every event" << channelMask_;
0138   if (hfInput > 0xFF)
0139     throw cms::Exception("Invalid Data") << "HF compressedET should be less than 0xFF, is " << hfInput;
0140   short iEta = rctParameters_->calcIEta(crtNo, crdNo, twrNo);
0141   unsigned short iAbsEta = abs(iEta);
0142   short sign = (iEta / iAbsEta);
0143   unsigned short phiSide = twrNo / 4;
0144   if (iAbsEta < 29 || iAbsEta > 32)
0145     throw cms::Exception("Invalid Data") << "29 <= |iEta| <= 32, is " << iAbsEta;
0146
0147   float et = convertHcal(hfInput, iAbsEta, sign);
0148   ;
0149
0150   if (channelMask_->hfMask[crtNo][phiSide][iAbsEta - 29] ||
0151       (noisyChannelMask_->hfMask[crtNo][phiSide][iAbsEta - 29] && et < noisyChannelMask_->hfThreshold)) {
0152     et = 0;
0153   }
0154
0155   // cout << (int) rctParameters_->jetMETHCalScaleFactors()[iAbsEta-1] << endl;
0156
0157   float scalehf = 1.;
0158   if (rctParameters_->jetMETHCalScaleFactors().size() == 32) {
0159     scalehf = (float)rctParameters_->jetMETHCalScaleFactors()[iAbsEta - 1];
0160   }  // The max eta for the various scale factors is 32, check to see if used.
0161   else if (rctParameters_->jetMETHCalScaleFactors().size() == 32 * 10) {
0162     int ht_bin = ((int)floor(et) / 5);
0163     // lowest bin (1) is 0-10GeV
0164     if (ht_bin < 1)
0165       ht_bin = 1;
0166     // highest bin (9) is 45GeV and up
0167     if (ht_bin > 9)
0168       ht_bin = 9;
0169     scalehf = (float)rctParameters_->jetMETHCalScaleFactors()[32 * ht_bin + iAbsEta - 1];
0170   }  // et-dependent scale factors (optional, of course, if set to 1)
0171
0172   et = scalehf * et;  // Allow for scaling the HF as well e.g. zero out if needed
0173
0174   unsigned int result = convertToInteger(et, rctParameters_->jetMETLSB(), 8);
0175   return result;
0176 }
0177
0178 bool L1RCTLookupTables::hOeFGVetoBit(float ecal, float hcal, bool fgbit) const {
0179   if (rctParameters_ == nullptr)
0180     throw cms::Exception("L1RCTParameters Invalid") << "L1RCTParameters should be set every event" << rctParameters_;
0181   bool veto = false;
0182   if (ecal > rctParameters_->eMinForFGCut() && ecal < rctParameters_->eMaxForFGCut()) {
0183     if (fgbit)
0184       veto = true;
0185   }
0186   if (ecal >= rctParameters_->eMinForHoECut() && ecal < rctParameters_->eMaxForHoECut()) {
0187     if ((hcal / ecal) > rctParameters_->hOeCut())
0188       veto = true;
0189   }
0190   //  else
0191   if (ecal < rctParameters_->eMinForHoECut()) {
0192     if (hcal >= rctParameters_->hMinForHoECut())
0193       veto = true;  // Changed from eMinForHoECut() - JLL 2008-Feb-13
0194   }
0195   return veto;
0196 }
0197
0198 bool L1RCTLookupTables::activityBit(float ecal, float hcal) const {
0199   // Redefined for upgrade as EM activity only
0200   if (rctParameters_ == nullptr)
0201     throw cms::Exception("L1RCTParameters Invalid") << "L1RCTParameters should be set every event" << rctParameters_;
0202   bool aBit = false;
0203   if (rctParameters_->eMinForHoECut() < rctParameters_->eMaxForHoECut()) {
0204     // For RCT operations HoE cut and tauVeto are used
0205     aBit = ((ecal > rctParameters_->eActivityCut()) || (hcal > rctParameters_->hActivityCut()));
0206   } else {
0207     // We redefine tauVeto() for upgrade as EM activity only  --
0208     // both EG and Tau make it through the EIC and JSC to CTP cards
0209     // In the CTP card we want to rechannel EG/Tau candidates to EG and Tau
0210     if (ecal > rctParameters_->eActivityCut()) {
0211       if ((hcal / ecal) < rctParameters_->hOeCut()) {
0212         aBit = true;
0213       }
0214     }
0215   }
0216   return aBit;
0217 }
0218
0219 // uses etScale
0220 unsigned int L1RCTLookupTables::emRank(unsigned short energy) const {
0221   if (etScale_) {
0222     return etScale_->rank(energy);
0223   } else
0224     //    edm::LogInfo("L1RegionalCaloTrigger")
0225     //      << "CaloEtScale was not used - energy instead of rank" << std::endl;
0226     return energy;
0227 }
0228
0229 // converts compressed ecal energy to linear (real) scale
0230 float L1RCTLookupTables::convertEcal(unsigned short ecal, unsigned short iAbsEta, short sign) const {
0231   if (ecalScale_) {
0232     // std::cout << "[luts] energy " << ecal << " sign " << sign
0233     //<< " iAbsEta " << iAbsEta << " iPhi " << iPhi << std::endl;
0234     float dummy = 0;
0235     dummy = float(ecalScale_->et(ecal, iAbsEta, sign));
0236     /*
0237     if (ecal > 0)
0238       {
0239         std::cout << "[luts] ecal converted from " << ecal << " to "
0240                   << dummy << " with iAbsEta " << iAbsEta << std::endl;
0241       }
0242     */
0243     return dummy;
0244   }
0245   // else if(rctParameters_ == 0)
0246   //  {
0247   //    throw cms::Exception("L1RCTParameters Invalid")
0248   //    << "L1RCTParameters should be set every event" << rctParameters_;
0249   //  }
0250   else {
0251     return ((float)ecal) * rctParameters_->eGammaLSB();
0252   }
0253 }
0254
0255 // converts compressed hcal energy to linear (real) scale
0256 float L1RCTLookupTables::convertHcal(unsigned short hcal, unsigned short iAbsEta, short sign) const {
0257   if (hcalScale_ != nullptr) {
0258     return (hcalScale_->et(hcal, iAbsEta, sign));
0259   } else {
0260     //      edm::LogInfo("L1RegionalCaloTrigger")
0261     //  << "CaloTPGTranscoder was not used" << std::endl;
0262     return ((float)hcal) * rctParameters_->jetMETLSB();
0263   }
0264 }
0265
0266 // integerize given an LSB and set maximum value of 2^precision-1
0267 unsigned long L1RCTLookupTables::convertToInteger(float et, float lsb, int precision) const {
0268   unsigned long etBits = (unsigned long)(et / lsb);
0269   unsigned long maxValue = (1 << precision) - 1;
0270   if (etBits > maxValue)
0271     return maxValue;
0272   else
0273     return etBits;
0274 }
0275
0276 unsigned int L1RCTLookupTables::eGammaETCode(float ecal, float hcal, int iAbsEta) const {
0277   if (rctParameters_ == nullptr)
0278     throw cms::Exception("L1RCTParameters Invalid") << "L1RCTParameters should be set every event" << rctParameters_;
0279   float etLinear = rctParameters_->EGammaTPGSum(ecal, hcal, iAbsEta);
0280   return convertToInteger(etLinear, rctParameters_->eGammaLSB(), 7);
0281 }
0282
0283 unsigned int L1RCTLookupTables::jetMETETCode(float ecal, float hcal, int iAbsEta) const {
0284   if (rctParameters_ == nullptr)
0285     throw cms::Exception("L1RCTParameters Invalid") << "L1RCTParameters should be set every event" << rctParameters_;
0286   float etLinear = rctParameters_->JetMETTPGSum(ecal, hcal, iAbsEta);
0287   return convertToInteger(etLinear, rctParameters_->jetMETLSB(), 9);
0288 }