Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-12 04:16:16 

0001 /**
0002  * Author: Sridhara Dasu
0003  * Created: 04 July 2007
0004  * $Id: L1RCTParameters.cc,v 1.24 2010/05/24 11:54:35 efron Exp $
0005  **/
0006 
0007 #include <iostream>
0008 #include <fstream>
0009 #include <cmath>
0010 
0011 #include "CondFormats/L1TObjects/interface/L1RCTParameters.h"
0012 #include <iomanip>
0013 
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 using namespace std;
0016 
0017 L1RCTParameters::L1RCTParameters(double eGammaLSB,
0018                                  double jetMETLSB,
0019                                  double eMinForFGCut,
0020                                  double eMaxForFGCut,
0021                                  double hOeCut,
0022                                  double eMinForHoECut,
0023                                  double eMaxForHoECut,
0024                                  double hMinForHoECut,
0025                                  double eActivityCut,
0026                                  double hActivityCut,
0027                                  unsigned eicIsolationThreshold,
0028                                  unsigned jscQuietThresholdBarrel,
0029                                  unsigned jscQuietThresholdEndcap,
0030                                  bool noiseVetoHB,
0031                                  bool noiseVetoHEplus,
0032                                  bool noiseVetoHEminus,
0033                                  bool useCorrections,
0034                                  const std::vector<double>& eGammaECalScaleFactors,
0035                                  const std::vector<double>& eGammaHCalScaleFactors,
0036                                  const std::vector<double>& jetMETECalScaleFactors,
0037                                  const std::vector<double>& jetMETHCalScaleFactors,
0038                                  const std::vector<double>& ecal_calib,
0039                                  const std::vector<double>& hcal_calib,
0040                                  const std::vector<double>& hcal_high_calib,
0041                                  const std::vector<double>& cross_terms,
0042                                  const std::vector<double>& lowHoverE_smear,
0043                                  const std::vector<double>& highHoverE_smear)
0044     : eGammaLSB_(eGammaLSB),
0045       jetMETLSB_(jetMETLSB),
0046       eMinForFGCut_(eMinForFGCut),
0047       eMaxForFGCut_(eMaxForFGCut),
0048       hOeCut_(hOeCut),
0049       eMinForHoECut_(eMinForHoECut),
0050       eMaxForHoECut_(eMaxForHoECut),
0051       hMinForHoECut_(hMinForHoECut),
0052       eActivityCut_(eActivityCut),
0053       hActivityCut_(hActivityCut),
0054       eicIsolationThreshold_(eicIsolationThreshold),
0055       jscQuietThresholdBarrel_(jscQuietThresholdBarrel),
0056       jscQuietThresholdEndcap_(jscQuietThresholdEndcap),
0057       noiseVetoHB_(noiseVetoHB),
0058       noiseVetoHEplus_(noiseVetoHEplus),
0059       noiseVetoHEminus_(noiseVetoHEminus),
0060       useCorrections_(useCorrections),
0061       eGammaECalScaleFactors_(eGammaECalScaleFactors),
0062       eGammaHCalScaleFactors_(eGammaHCalScaleFactors),
0063       jetMETECalScaleFactors_(jetMETECalScaleFactors),
0064       jetMETHCalScaleFactors_(jetMETHCalScaleFactors),
0065       HoverE_smear_low_(lowHoverE_smear),
0066       HoverE_smear_high_(highHoverE_smear) {
0067   ecal_calib_.resize(28);
0068   hcal_calib_.resize(28);
0069   hcal_high_calib_.resize(28);
0070   cross_terms_.resize(28);
0071 
0072   for (unsigned i = 0; i < ecal_calib.size(); ++i)
0073     ecal_calib_[i / 3].push_back(ecal_calib[i]);
0074   for (unsigned i = 0; i < hcal_calib.size(); ++i)
0075     hcal_calib_[i / 3].push_back(hcal_calib[i]);
0076   for (unsigned i = 0; i < hcal_high_calib.size(); ++i)
0077     hcal_high_calib_[i / 3].push_back(hcal_high_calib[i]);
0078   for (unsigned i = 0; i < cross_terms.size(); ++i)
0079     cross_terms_[i / 6].push_back(cross_terms[i]);
0080 }
0081 
0082 // maps rct iphi, ieta of tower to crate
0083 unsigned short L1RCTParameters::calcCrate(unsigned short rct_iphi, short ieta) const {
0084   unsigned short crate = rct_iphi / 8;
0085   if (abs(ieta) > 28)
0086     crate = rct_iphi / 2;
0087   if (ieta > 0) {
0088     crate = crate + 9;
0089   }
0090   return crate;
0091 }
0092 
0093 //map digi rct iphi, ieta to card
0094 unsigned short L1RCTParameters::calcCard(unsigned short rct_iphi, unsigned short absIeta) const {
0095   unsigned short card = 999;
0096   // Note absIeta counts from 1-32 (not 0-31)
0097   if (absIeta <= 24) {
0098     card = ((absIeta - 1) / 8) * 2 + (rct_iphi % 8) / 4;
0099   }
0100   // 25 <= absIeta <= 28 (card 6)
0101   else if ((absIeta >= 25) && (absIeta <= 28)) {
0102     card = 6;
0103   } else {
0104   }
0105   return card;
0106 }
0107 
0108 //map digi rct iphi, ieta to tower
0109 unsigned short L1RCTParameters::calcTower(unsigned short rct_iphi, unsigned short absIeta) const {
0110   unsigned short tower = 999;
0111   unsigned short iphi = rct_iphi;
0112   unsigned short regionPhi = (iphi % 8) / 4;
0113 
0114   // Note absIeta counts from 1-32 (not 0-31)
0115   if (absIeta <= 24) {
0116     // assume iphi between 0 and 71; makes towers from 0-31, mod. 7Nov07
0117     tower = ((absIeta - 1) % 8) * 4 + (iphi % 4);  // REMOVED +1
0118   }
0119   // 25 <= absIeta <= 28 (card 6)
0120   else if ((absIeta >= 25) && (absIeta <= 28)) {
0121     if (regionPhi == 0) {
0122       // towers from 0-31, modified 7Nov07 Jessica Leonard
0123       tower = (absIeta - 25) * 4 + (iphi % 4);  // REMOVED +1
0124     } else {
0125       tower = 28 + iphi % 4 + (25 - absIeta) * 4;  // mod. 7Nov07 JLL
0126     }
0127   }
0128   // absIeta >= 29 (HF regions)
0129   else if ((absIeta >= 29) && (absIeta <= 32)) {
0130     // SPECIAL DEFINITION OF REGIONPHI FOR HF SINCE HF IPHI IS 0-17
0131     // Sept. 19 J. Leonard
0132     regionPhi = iphi % 2;
0133     // HF MAPPING, just regions now, don't need to worry about towers
0134     // just calling it "tower" for convenience
0135     tower = (regionPhi) * 4 + absIeta - 29;
0136   }
0137   return tower;
0138 }
0139 
0140 // iCrate 0-17, iCard 0-6, NEW iTower 0-31
0141 short L1RCTParameters::calcIEta(unsigned short iCrate, unsigned short iCard, unsigned short iTower) const {
0142   unsigned short absIEta = calcIAbsEta(iCrate, iCard, iTower);
0143   short iEta;
0144   if (iCrate < 9)
0145     iEta = -absIEta;
0146   else
0147     iEta = absIEta;
0148   return iEta;
0149 }
0150 
0151 // iCrate 0-17, iCard 0-6, NEW iTower 0-31
0152 unsigned short L1RCTParameters::calcIPhi(unsigned short iCrate, unsigned short iCard, unsigned short iTower) const {
0153   short iPhi;
0154   if (iCard < 6)
0155     iPhi = (iCrate % 9) * 8 + (iCard % 2) * 4 + (iTower % 4);  // rm -1 7Nov07
0156   else if (iCard == 6) {
0157     // region 0
0158     if (iTower < 16)                           // 17->16
0159       iPhi = (iCrate % 9) * 8 + (iTower % 4);  // rm -1
0160     // region 1
0161     else
0162       iPhi = (iCrate % 9) * 8 + ((iTower - 16) % 4) + 4;  // 17 -> 16
0163   }
0164   // HF regions
0165   else
0166     iPhi = (iCrate % 9) * 2 + iTower / 4;
0167   return iPhi;
0168 }
0169 
0170 // iCrate 0-17, iCard 0-6, NEW iTower 0-31
0171 unsigned short L1RCTParameters::calcIAbsEta(unsigned short iCrate, unsigned short iCard, unsigned short iTower) const {
0172   unsigned short absIEta;
0173   if (iCard < 6)
0174     absIEta = (iCard / 2) * 8 + (iTower / 4) + 1;  // rm -1 JLL 7Nov07
0175   else if (iCard == 6) {
0176     // card 6, region 0
0177     if (iTower < 16)              // 17->16
0178       absIEta = 25 + iTower / 4;  // rm -1
0179     // card 6, region 1
0180     else
0181       absIEta = 28 - ((iTower - 16) / 4);  // 17->16
0182   }
0183   // HF regions
0184   else
0185     absIEta = 29 + iTower % 4;
0186   return absIEta;
0187 }
0188 
0189 float L1RCTParameters::JetMETTPGSum(const float& ecal, const float& hcal, const unsigned& iAbsEta) const {
0190   // We never deal with HF in this function (see note below)
0191   if (iAbsEta < 1 || iAbsEta > 28)
0192     throw cms::Exception("L1RCTParameters invalid function call") << "Eta out of range in MET TPGSum: " << iAbsEta;
0193   float ecal_c = ecal * jetMETECalScaleFactors_.at(iAbsEta - 1);
0194   float hcal_c = hcal * jetMETHCalScaleFactors_.at(iAbsEta - 1);
0195 
0196   // scale factors will either be length 28 for legacy, or 28*(# et bins + 1) where the first set is an average over the et bins
0197   // The first set provides a legacy fallthrough option
0198   // Currently, # et bins is 9
0199   if (jetMETECalScaleFactors_.size() == 28 * 10) {
0200     int et_bin = ((int)floor(ecal) / 5);
0201     // lowest bin (1) is 0-10GeV
0202     if (et_bin < 1)
0203       et_bin = 1;
0204     // highest bin (9) is 45GeV and up
0205     if (et_bin > 9)
0206       et_bin = 9;
0207     ecal_c = ecal * jetMETECalScaleFactors_.at(et_bin * 28 + iAbsEta - 1);
0208   }
0209 
0210   // We may be interested in HF jets, in which case, there are four more scale factors
0211   // for the 4 HF regions.  HOWEVER, we will never expect to see them accessed from
0212   // this function.  HF scaling is done in L1Trigger/RegionalCaloTrigger/src/L1RCTLookupTables
0213   if (jetMETHCalScaleFactors_.size() == 32 * 10) {
0214     int ht_bin = ((int)floor(hcal) / 5);
0215     // lowest bin (1) is 0-10GeV
0216     if (ht_bin < 1)
0217       ht_bin = 1;
0218     // highest bin (9) is 45GeV and up
0219     if (ht_bin > 9)
0220       ht_bin = 9;
0221     hcal_c = hcal * jetMETHCalScaleFactors_.at(ht_bin * 32 + iAbsEta - 1);
0222   }
0223 
0224   float result = ecal_c + hcal_c;
0225 
0226   // defunct section (polynomial-parameterized corrections)
0227   if (useCorrections_) {
0228     if (jetMETHCalScaleFactors_.at(iAbsEta - 1) != 0)
0229       hcal_c = hcal;
0230 
0231     if (jetMETECalScaleFactors_.at(iAbsEta - 1) != 0)
0232       ecal_c = ecal * eGammaECalScaleFactors_.at(iAbsEta - 1);  // Use eGamma Corrections
0233 
0234     result = correctedTPGSum(ecal_c, hcal_c, iAbsEta - 1);
0235   }
0236 
0237   return result;
0238 }
0239 
0240 float L1RCTParameters::EGammaTPGSum(const float& ecal, const float& hcal, const unsigned& iAbsEta) const {
0241   // We never deal with HF in this function (EG objects don't use hcal at all, really)
0242   if (iAbsEta < 1 || iAbsEta > 28)
0243     throw cms::Exception("L1RCTParameters invalid function call") << "Eta out of range in MET TPGSum: " << iAbsEta;
0244   float ecal_c = ecal * eGammaECalScaleFactors_.at(iAbsEta - 1);
0245   float hcal_c = hcal * eGammaHCalScaleFactors_.at(iAbsEta - 1);
0246 
0247   // scale factors will either be length 28 for legacy, or 28*(# et bins + 1) where the first set of 28 is an average over the et bins
0248   // The first set of 28 provides a legacy fallthrough option
0249   // Currently, # et bins is 9
0250   if (eGammaECalScaleFactors_.size() == 28 * 10) {
0251     int et_bin = ((int)floor(ecal) / 5);
0252     // lowest bin (1) is 0-10GeV
0253     if (et_bin < 1)
0254       et_bin = 1;
0255     // highest bin (9) is 45GeV and up
0256     if (et_bin > 9)
0257       et_bin = 9;
0258     ecal_c = ecal * eGammaECalScaleFactors_.at(et_bin * 28 + iAbsEta - 1);
0259   }
0260   if (eGammaHCalScaleFactors_.size() == 28 * 10) {
0261     int ht_bin = ((int)floor(hcal) / 5);
0262     // lowest bin (1) is 0-10GeV
0263     if (ht_bin < 1)
0264       ht_bin = 1;
0265     // highest bin (9) is 45GeV and up
0266     if (ht_bin > 9)
0267       ht_bin = 9;
0268     hcal_c = hcal * eGammaHCalScaleFactors_.at(ht_bin * 28 + iAbsEta - 1);
0269   }
0270 
0271   float result = ecal_c + hcal_c;
0272 
0273   // defunct section (polynomial-parameterized corrections)
0274   if (useCorrections_) {
0275     if (eGammaHCalScaleFactors_.at(iAbsEta - 1) != 0)
0276       hcal_c = hcal;
0277 
0278     result = correctedTPGSum(ecal_c, hcal_c, iAbsEta - 1);
0279   }
0280 
0281   return result;
0282 }
0283 
0284 // index = iAbsEta - 1... make sure you call the function like so: "correctedTPGSum(ecal,hcal, iAbsEta - 1)"
0285 float L1RCTParameters::correctedTPGSum(const float& ecal, const float& hcal, const unsigned& index) const {
0286   if (index >= 28 && ecal > 120 && hcal > 120)
0287     return (ecal + hcal);  // return plain sum if outside of calibration range or index is too high
0288 
0289   // let's make sure we're asking for items that are there.
0290   if (ecal_calib_.at(index).size() != 3 || hcal_calib_.at(index).size() != 3 ||
0291       hcal_high_calib_.at(index).size() != 3 || cross_terms_.at(index).size() != 6 ||
0292       HoverE_smear_high_.size() <= index || HoverE_smear_low_.size() <= index)
0293     return (ecal + hcal);
0294 
0295   double e = ecal, h = hcal;
0296   double ec = 0.0, hc = 0.0, c = 0.0;
0297 
0298   ec = (ecal_calib_.at(index).at(0) * std::pow(e, 3.) + ecal_calib_.at(index).at(1) * std::pow(e, 2.) +
0299         ecal_calib_.at(index).at(2) * e);
0300 
0301   if (e + h < 23) {
0302     hc = (hcal_calib_.at(index).at(0) * std::pow(h, 3.) + hcal_calib_.at(index).at(1) * std::pow(h, 2.) +
0303           hcal_calib_.at(index).at(2) * h);
0304 
0305     c = (cross_terms_.at(index).at(0) * std::pow(e, 2.) * h + cross_terms_.at(index).at(1) * std::pow(h, 2.) * e +
0306          cross_terms_.at(index).at(2) * e * h + cross_terms_.at(index).at(3) * std::pow(e, 3.) * h +
0307          cross_terms_.at(index).at(4) * std::pow(h, 3.) * e +
0308          cross_terms_.at(index).at(5) * std::pow(h, 2.) * std::pow(e, 2.));
0309   } else {
0310     hc = (hcal_high_calib_.at(index).at(0) * std::pow(h, 3.) + hcal_high_calib_.at(index).at(1) * std::pow(h, 2.) +
0311           hcal_high_calib_.at(index).at(2) * h);
0312   }
0313 
0314   if (h / (e + h) >= 0.05) {
0315     ec *= HoverE_smear_high_.at(index);
0316     hc *= HoverE_smear_high_.at(index);
0317     c *= HoverE_smear_high_.at(index);
0318   } else {
0319     ec *= HoverE_smear_low_.at(index);
0320   }
0321   return ec + hc + c;
0322 }
0323 
0324 void L1RCTParameters::print(std::ostream& s) const {
0325   s << "\nPrinting record L1RCTParametersRcd" << endl;
0326   s << "\n\"Parameter description\" \n  \"Parameter name\"  \"Value\" " << endl;
0327   s << "\ne/gamma least significant bit energy transmitted from receiver cards to EIC cards. \n  "
0328     << "eGammaLSB = " << eGammaLSB_ << endl;
0329   s << "\nLSB of region Et scale from RCT to GCT (GeV) \n  "
0330     << "jetMETLSB = " << jetMETLSB_ << endl;
0331   s << "\nminimum ECAL Et for which fine-grain veto is applied (GeV) \n  "
0332     << " eMinForFGCut = " << eMinForFGCut_ << endl;
0333   s << "\nmaximum ECAL Et for which fine-grain veto is applied (GeV) \n  "
0334     << "eMaxForFGCut = " << eMaxForFGCut_ << endl;
0335   s << "\nmaximum value of (HCAL Et / ECAL Et) \n  "
0336     << "hOeCut = " << hOeCut_ << endl;
0337   s << "\nminimum ECAL Et for which H/E veto is applied (GeV) \n  "
0338     << "eMinForHoECut = " << eMinForHoECut_ << endl;
0339   s << "\nmaximum ECAL Et for which H/E veto is applied (GeV) \n  "
0340     << "eMaxForHoECut = " << eMaxForHoECut_ << endl;
0341   s << "\nminimum HCAL Et for which H/E veto is applied (GeV)  \n  "
0342     << "hMinForHoECut = " << hMinForHoECut_ << endl;
0343   s << "\nECAL Et threshold above which tau activity bit is set (GeV)  \n  "
0344     << "eActivityCut = " << eActivityCut_ << endl;
0345   s << "\nHCAL Et threshold above which tau activity bit is set (GeV)  \n  "
0346     << "hActivityCut = " << hActivityCut_ << endl;
0347   s << "\nNeighboring trigger tower energy minimum threshold that marks candidate as non-isolated. (LSB bits) \n  "
0348     << "eicIsolationThreshold = " << eicIsolationThreshold_ << endl;
0349   s << "\nIf jetMet energy in RCT Barrel Region is below this value, a quiet bit is set. (LSB bits)\n  "
0350     << "jscQuietThreshBarrel = " << jscQuietThresholdBarrel_ << endl;
0351   s << "\nIf jetMet energy in RCT Endcap Region is below this value, a quiet bit is set. (LSB bits) \n  "
0352     << "jscQuietThreshEndcap = " << jscQuietThresholdEndcap_ << endl;
0353   s << "\nWhen set to TRUE, HCAL energy is ignored if no ECAL energy is present in corresponding trigger tower for RCT "
0354        "Barrel \n  "
0355     << "noiseVetoHB = " << noiseVetoHB_ << endl;
0356   s << "\nWhen set to TRUE, HCAL energy is ignored if no ECAL energy is present in corresponding trigger tower for RCT "
0357        "Encap+ \n  "
0358     << "noiseVetoHEplus = " << noiseVetoHEplus_ << endl;
0359   s << "\nWhen set to TRUE, HCAL energy is ignored if no ECAL energy is present in corresponding trigger tower for RCT "
0360        "Endcap- \n  "
0361     << "noiseVetoHEminus = " << noiseVetoHEminus_ << endl;
0362 
0363   auto printScalefactors = [&s](const std::vector<double>& sf) {
0364     if (sf.size() == 10 * 28) {
0365       s << "et bin  ieta  ScaleFactor" << endl;
0366       for (unsigned i = 0; i < sf.size(); i++)
0367         s << setw(6) << i / 28 << "  " << setw(4) << i % 28 + 1 << "  " << sf.at(i) << endl;
0368     } else if (sf.size() == 10 * 32)  // jet HCAL (HF regions are 29-32)
0369     {
0370       s << "et bin  ieta  ScaleFactor" << endl;
0371       for (unsigned i = 0; i < sf.size(); i++)
0372         s << setw(6) << i / 32 << "  " << setw(4) << i % 32 + 1 << "  " << sf.at(i) << endl;
0373     } else {
0374       s << "ieta  ScaleFactor" << endl;
0375       for (unsigned i = 0; i < sf.size(); i++)
0376         s << setw(4) << i + 1 << "  " << sf.at(i) << endl;
0377     }
0378   };
0379 
0380   s << "\n\neta-dependent multiplicative factors for ECAL Et before summation \n  "
0381     << "eGammaECal Scale Factors " << endl;
0382   printScalefactors(eGammaECalScaleFactors_);
0383 
0384   s << "\n\neta-dependent multiplicative factors for HCAL Et before summation \n  "
0385     << "eGammaHCal Scale Factors " << endl;
0386   printScalefactors(eGammaHCalScaleFactors_);
0387 
0388   s << "\n\neta-dependent multiplicative factors for ECAL Et before summation \n  "
0389     << "jetMETECal Scale Factors " << endl;
0390   printScalefactors(jetMETECalScaleFactors_);
0391 
0392   s << "\n\neta-dependent multiplicative factors for HCAL Et before summation \n"
0393     << "jetMETHCal Scale Factors " << endl;
0394   printScalefactors(jetMETHCalScaleFactors_);
0395 
0396   if (useCorrections_) {
0397     s << "\n\nUSING calibration variables " << endl;
0398 
0399     s << "\n\nH over E smear low Correction Factors " << endl;
0400     s << "ieta  Correction Factor" << endl;
0401     for (int i = 0; i < 28; i++)
0402       s << setw(4) << i + 1 << "  " << HoverE_smear_low_.at(i) << endl;
0403 
0404     s << "\n\nH over E smear high Correction Factors " << endl;
0405     s << "ieta  Correction Factor" << endl;
0406     for (int i = 0; i < 28; i++)
0407       s << setw(4) << i + 1 << "  " << HoverE_smear_high_.at(i) << endl;
0408 
0409     s << "\n\necal calibrations " << endl;
0410     s << "ieta  CorrFactor1  CorrFactor2  CorrFactor3" << endl;
0411     int end = ecal_calib_[0].size();
0412     for (int i = 0; i < 28; i++) {
0413       s << setw(4) << i;
0414       for (int j = 0; j < end; j++)
0415         s << setw(11) << setprecision(8) << ecal_calib_[i][j];
0416 
0417       s << endl;
0418     }
0419 
0420     s << "\n\nhcal calibrations " << endl;
0421     s << "ieta  CorrFactor1  CorrFactor2  CorrFactor3" << endl;
0422     end = hcal_calib_[0].size();
0423     for (int i = 0; i < 28; i++) {
0424       s << setw(4) << i;
0425       for (int j = 0; j < end; j++)
0426         s << setw(11) << setprecision(8) << hcal_calib_[i][j];
0427 
0428       s << endl;
0429     }
0430     s << "\n\nhcal_high calibrations " << endl;
0431     s << "ieta  CorrFactor1  CorrFactor2  CorrFactor3" << endl;
0432     end = hcal_high_calib_[0].size();
0433     for (int i = 0; i < 28; i++) {
0434       s << setw(4) << i;
0435       for (int j = 0; j < end; j++)
0436         s << setw(11) << setprecision(8) << hcal_high_calib_[i][j];
0437 
0438       s << endl;
0439     }
0440     end = cross_terms_[0].size();
0441     s << "\n\ncross terms calibrations " << endl;
0442     s << "ieta  CorrFactor1  CorrFactor2  CorrFactor3  CorrFactor4  CorrFactor5  CorrFactor6" << endl;
0443     for (int i = 0; i < 28; i++) {
0444       s << setw(4) << i;
0445       for (int j = 0; j < end; j++)
0446         s << setw(11) << setprecision(8) << cross_terms_[i][j];
0447 
0448       s << endl;
0449     }
0450 
0451   } else
0452     s << "\n\nNOT USING calibration variables " << endl;
0453 
0454   s << "\n\n" << endl;
0455 }
[ 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!