Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:02:22 

0001 #include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
0002 
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <cmath>
0006 
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "DataFormats/L1GlobalCaloTrigger/interface/L1GctStaticParameters.h"
0011 
0012 using std::ios;
0013 
0014 const unsigned L1GctJetFinderParams::NUMBER_ETA_VALUES = 11;
0015 const unsigned L1GctJetFinderParams::N_CENTRAL_ETA_VALUES = 7;
0016 
0017 L1GctJetFinderParams::L1GctJetFinderParams()
0018     : rgnEtLsb_(0.),
0019       htLsb_(0.),
0020       cenJetEtSeed_(0.),
0021       forJetEtSeed_(0.),
0022       tauJetEtSeed_(0.),
0023       tauIsoEtThreshold_(0.),
0024       htJetEtThreshold_(0.),
0025       mhtJetEtThreshold_(0.),
0026       cenForJetEtaBoundary_(0),
0027       corrType_(0),
0028       jetCorrCoeffs_(),
0029       tauCorrCoeffs_(),
0030       convertToEnergy_(false),
0031       energyConversionCoeffs_() {}
0032 
0033 L1GctJetFinderParams::L1GctJetFinderParams(double rgnEtLsb,
0034                                            double htLsb,
0035                                            double cJetSeed,
0036                                            double fJetSeed,
0037                                            double tJetSeed,
0038                                            double tauIsoEtThresh,
0039                                            double htJetEtThresh,
0040                                            double mhtJetEtThresh,
0041                                            unsigned etaBoundary,
0042                                            unsigned corrType,
0043                                            const std::vector<std::vector<double> >& jetCorrCoeffs,
0044                                            const std::vector<std::vector<double> >& tauCorrCoeffs,
0045                                            bool convertToEnergy,
0046                                            const std::vector<double>& energyConvCoeffs)
0047     : rgnEtLsb_(rgnEtLsb),
0048       htLsb_(htLsb),
0049       cenJetEtSeed_(cJetSeed),
0050       forJetEtSeed_(fJetSeed),
0051       tauJetEtSeed_(tJetSeed),
0052       tauIsoEtThreshold_(tauIsoEtThresh),
0053       htJetEtThreshold_(htJetEtThresh),
0054       mhtJetEtThreshold_(mhtJetEtThresh),
0055       cenForJetEtaBoundary_(etaBoundary),
0056       corrType_(corrType),
0057       jetCorrCoeffs_(jetCorrCoeffs),
0058       tauCorrCoeffs_(tauCorrCoeffs),
0059       convertToEnergy_(convertToEnergy),
0060       energyConversionCoeffs_(energyConvCoeffs) {
0061   // check number of eta bins
0062   if (jetCorrCoeffs_.size() != NUMBER_ETA_VALUES || tauCorrCoeffs_.size() != N_CENTRAL_ETA_VALUES ||
0063       energyConversionCoeffs_.size() != NUMBER_ETA_VALUES) {
0064     LogDebug("L1-O2O") << "GCT jet corrections constructed with " << jetCorrCoeffs_.size() << " bins, expected "
0065                        << NUMBER_ETA_VALUES << std::endl;
0066     LogDebug("L1-O2O") << "GCT tau corrections constructed with " << tauCorrCoeffs_.size() << " bins, expected "
0067                        << N_CENTRAL_ETA_VALUES << std::endl;
0068     LogDebug("L1-O2O") << "GCT energy corrections constructed with " << energyConversionCoeffs_.size()
0069                        << " bins, expected " << NUMBER_ETA_VALUES << std::endl;
0070 
0071     throw cms::Exception("InconsistentConfig")
0072         << "L1GctJetFinderParams constructed with wrong number of eta bins : " << jetCorrCoeffs_.size() << " jets, "
0073         << tauCorrCoeffs_.size() << " taus, " << energyConversionCoeffs_.size() << " energy conversion bins"
0074         << std::endl;
0075   }
0076 
0077   // check number of coefficients against expectation
0078   unsigned expCoeffs = 0;
0079   if (corrType_ == 2)
0080     expCoeffs = 8;  // ORCA style
0081   if (corrType_ == 3)
0082     expCoeffs = 4;  // Simple
0083   if (corrType_ == 4)
0084     expCoeffs = 15;  // piecewise-cubic
0085   if (corrType_ == 5)
0086     expCoeffs = 6;  // PF
0087 
0088   // correction types 1 and 4 can have any number of parameters
0089   if (corrType_ != 1 && corrType_ != 4) {
0090     std::vector<std::vector<double> >::const_iterator itr;
0091     for (itr = jetCorrCoeffs_.begin(); itr != jetCorrCoeffs_.end(); ++itr) {
0092       if (itr->size() != expCoeffs) {
0093         throw cms::Exception("InconsistentConfig")
0094             << "L1GctJetFinderParams constructed with " << itr->size() << " jet correction coefficients, when "
0095             << expCoeffs << " expected" << std::endl;
0096       }
0097     }
0098     for (itr = tauCorrCoeffs_.begin(); itr != tauCorrCoeffs_.end(); ++itr) {
0099       if (itr->size() != expCoeffs) {
0100         throw cms::Exception("InconsistentConfig")
0101             << "L1GctJetFinderParams constructed with " << itr->size() << " tau correction coefficients, when "
0102             << expCoeffs << " expected" << std::endl;
0103       }
0104     }
0105   }
0106 }
0107 
0108 L1GctJetFinderParams::~L1GctJetFinderParams() {}
0109 
0110 //---------------------------------------------------------------------------------------------
0111 //
0112 // set methods
0113 //
0114 
0115 void L1GctJetFinderParams::setRegionEtLsb(const double rgnEtLsb) { rgnEtLsb_ = rgnEtLsb; }
0116 
0117 void L1GctJetFinderParams::setSlidingWindowParams(const double cJetSeed,
0118                                                   const double fJetSeed,
0119                                                   const double tJetSeed,
0120                                                   const unsigned etaBoundary) {
0121   cenJetEtSeed_ = cJetSeed;
0122   forJetEtSeed_ = fJetSeed;
0123   tauJetEtSeed_ = tJetSeed;
0124   cenForJetEtaBoundary_ = etaBoundary;
0125 }
0126 
0127 void L1GctJetFinderParams::setJetEtCalibrationParams(const unsigned corrType,
0128                                                      const std::vector<std::vector<double> >& jetCorrCoeffs,
0129                                                      const std::vector<std::vector<double> >& tauCorrCoeffs) {
0130   corrType_ = corrType;
0131   jetCorrCoeffs_ = jetCorrCoeffs;
0132   tauCorrCoeffs_ = tauCorrCoeffs;
0133 }
0134 
0135 void L1GctJetFinderParams::setJetEtConvertToEnergyOn(const std::vector<double>& energyConvCoeffs) {
0136   convertToEnergy_ = true;
0137   energyConversionCoeffs_ = energyConvCoeffs;
0138 }
0139 
0140 void L1GctJetFinderParams::setJetEtConvertToEnergyOff() {
0141   convertToEnergy_ = false;
0142   energyConversionCoeffs_.clear();
0143 }
0144 
0145 void L1GctJetFinderParams::setHtSumParams(const double htLsb, const double htJetEtThresh, const double mhtJetEtThresh) {
0146   htLsb_ = htLsb;
0147   htJetEtThreshold_ = htJetEtThresh;
0148   mhtJetEtThreshold_ = mhtJetEtThresh;
0149 }
0150 
0151 void L1GctJetFinderParams::setTauAlgorithmParams(const double tauIsoEtThresh) { tauIsoEtThreshold_ = tauIsoEtThresh; }
0152 
0153 void L1GctJetFinderParams::setParams(const double rgnEtLsb,
0154                                      const double htLsb,
0155                                      const double cJetSeed,
0156                                      const double fJetSeed,
0157                                      const double tJetSeed,
0158                                      const double tauIsoEtThresh,
0159                                      const double htJetEtThresh,
0160                                      const double mhtJetEtThresh,
0161                                      const unsigned etaBoundary,
0162                                      const unsigned corrType,
0163                                      const std::vector<std::vector<double> >& jetCorrCoeffs,
0164                                      const std::vector<std::vector<double> >& tauCorrCoeffs) {
0165   setRegionEtLsb(rgnEtLsb);
0166   setSlidingWindowParams(cJetSeed, fJetSeed, tJetSeed, etaBoundary);
0167   setJetEtCalibrationParams(corrType, jetCorrCoeffs, tauCorrCoeffs);
0168   setHtSumParams(htLsb, htJetEtThresh, mhtJetEtThresh);
0169   setTauAlgorithmParams(tauIsoEtThresh);
0170 }
0171 
0172 //---------------------------------------------------------------------------------------------
0173 //
0174 // Jet Et correction methods
0175 //
0176 
0177 double L1GctJetFinderParams::correctedEtGeV(const double et, const unsigned eta, const bool tauVeto) const {
0178   if (eta >= NUMBER_ETA_VALUES) {
0179     return 0;
0180   } else {
0181     double result = 0;
0182     if ((eta >= cenForJetEtaBoundary_) || tauVeto) {
0183       // Use jetCorrCoeffs for central and forward jets.
0184       // In forward eta bins we ignore the tau flag (as in the firmware)
0185       result = correctionFunction(et, jetCorrCoeffs_.at(eta));
0186     } else {
0187       // Use tauCorrCoeffs for tau jets (in central eta bins)
0188       result = correctionFunction(et, tauCorrCoeffs_.at(eta));
0189     }
0190     if (convertToEnergy_) {
0191       result *= energyConversionCoeffs_.at(eta);
0192     }
0193     return result;
0194   }
0195 }
0196 
0197 /// Convert the corrected Et value to an integer Et for Ht summing
0198 uint16_t L1GctJetFinderParams::correctedEtGct(const double correctedEt) const {
0199   double scaledEt = correctedEt / htLsb_;
0200 
0201   uint16_t jetEtOut = static_cast<uint16_t>(scaledEt);
0202 
0203   if (jetEtOut > L1GctStaticParameters::jetCalibratedEtMax) {
0204     return L1GctStaticParameters::jetCalibratedEtMax;
0205   } else {
0206     return jetEtOut;
0207   }
0208 }
0209 
0210 // private methods
0211 
0212 double L1GctJetFinderParams::correctionFunction(const double Et, const std::vector<double>& coeffs) const {
0213   double result = 0;
0214   switch (corrType_) {
0215     case 0:  // no correction
0216       result = Et;
0217       break;
0218     case 1:  // power series correction
0219       result = powerSeriesCorrect(Et, coeffs);
0220       break;
0221     case 2:  // ORCA style correction
0222       result = orcaStyleCorrect(Et, coeffs);
0223       break;
0224     case 3:  // simple correction (JetMEt style)
0225       result = simpleCorrect(Et, coeffs);
0226       break;
0227     case 4:  // piecwise cubic correction a la Greg Landsberg et al
0228       result = piecewiseCubicCorrect(Et, coeffs);
0229       break;
0230     case 5:  // PF correction
0231       result = pfCorrect(Et, coeffs);
0232       break;
0233     default:
0234       result = Et;
0235   }
0236   return result;
0237 }
0238 
0239 double L1GctJetFinderParams::powerSeriesCorrect(const double Et, const std::vector<double>& coeffs) const {
0240   double corrEt = Et;
0241   for (unsigned i = 0; i < coeffs.size(); i++) {
0242     corrEt += coeffs.at(i) * pow(Et, (int)i);
0243   }
0244   return corrEt;
0245 }
0246 
0247 double L1GctJetFinderParams::orcaStyleCorrect(const double Et, const std::vector<double>& coeffs) const {
0248   // The coefficients are arranged in groups of four
0249   // The first in each group is a threshold value of Et.
0250 
0251   std::vector<double>::const_iterator next_coeff = coeffs.begin();
0252   while (next_coeff != coeffs.end()) {
0253     double threshold = *next_coeff++;
0254     double A = *next_coeff++;
0255     double B = *next_coeff++;
0256     double C = *next_coeff++;
0257     if (Et > threshold) {
0258       // This function is an inverse quadratic:
0259       //   (input Et) = A + B*(output Et) + C*(output Et)^2
0260       return 2 * (Et - A) / (B + sqrt(B * B - 4 * A * C + 4 * Et * C));
0261     }
0262     // If we are below all specified thresholds (or the vector is empty), return output=input.
0263   }
0264   return Et;
0265 }
0266 
0267 double L1GctJetFinderParams::simpleCorrect(const double Et, const std::vector<double>& coeffs) const {
0268   // function is :
0269   //    Et_out = A + B/[ (log10(Et_in))^C + D ] + E/Et_in
0270   //
0271   //    fitcor_as_str = "[0]+[1]/(pow(log10(x),[2])+[3]) + [4]/x"
0272   //
0273 
0274   return coeffs.at(0) + coeffs.at(1) / (pow(log10(Et), coeffs.at(2)) + coeffs.at(3)) + (coeffs.at(4) / Et);
0275 }
0276 
0277 double L1GctJetFinderParams::piecewiseCubicCorrect(const double Et, const std::vector<double>& coeffs) const {
0278   // The correction fuction is a set of 3rd order polynomials
0279   //    Et_out = Et_in + (p0 + p1*Et_in + p2*Et_in^2 + p3*Et_in^3)
0280   // with different coefficients for different energy ranges.
0281   // The parameters are arranged in groups of five.
0282   // The first in each group is a threshold value of input Et,
0283   // followed by the four coefficients for the cubic function.
0284   double etOut = Et;
0285   std::vector<double>::const_iterator next_coeff = coeffs.begin();
0286   while (next_coeff != coeffs.end()) {
0287     // Read the coefficients from the vector
0288     double threshold = *next_coeff++;
0289     double A = *next_coeff++;  //p0
0290     double B = *next_coeff++;  //p1
0291     double C = *next_coeff++;  //p2
0292     double D = *next_coeff++;  //p3
0293 
0294     // Check we are in the right energy range and make correction
0295     if (Et > threshold) {
0296       etOut += (A + etOut * (B + etOut * (C + etOut * D)));
0297       break;
0298     }
0299   }
0300   return etOut;
0301 }
0302 
0303 double L1GctJetFinderParams::pfCorrect(const double et, const std::vector<double>& coeffs) const {
0304   //
0305   // corr_factor = [0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))
0306   // Et_out = Et_in * corr_factor
0307   //
0308 
0309   return et * (coeffs.at(0) + coeffs.at(1) / (pow(log10(et), 2) + coeffs.at(2)) +
0310                coeffs.at(3) * exp(-coeffs.at(4) * (log10(et) - coeffs.at(5)) * (log10(et) - coeffs.at(5))));
0311 }
0312 
0313 std::ostream& operator<<(std::ostream& os, const L1GctJetFinderParams& fn) {
0314   //  os << std::setprecision(2);
0315 
0316   os << "=== Level-1 GCT : Jet Finder Parameters  ===" << std::endl;
0317   os << "RCT region LSB               : " << std::fixed << fn.getRgnEtLsbGeV() << " GeV" << std::endl;
0318   os << "Central jet seed threshold   : " << std::fixed << fn.getCenJetEtSeedGeV() << " GeV" << std::endl;
0319   os << "Tau jet seed threshold       : " << std::fixed << fn.getTauJetEtSeedGeV() << " GeV" << std::endl;
0320   os << "Forward jet seed threshold   : " << std::fixed << fn.getForJetEtSeedGeV() << " GeV" << std::endl;
0321   os << "Tau isolation threshold      : " << std::fixed << fn.getTauIsoEtThresholdGeV() << " GeV" << std::endl;
0322   os << "Ht jet Et threshold          : " << std::fixed << fn.getHtJetEtThresholdGeV() << " GeV" << std::endl;
0323   os << "MHt jet Et threshold         : " << std::fixed << fn.getMHtJetEtThresholdGeV() << " GeV" << std::endl;
0324   os << "Ht LSB                       : " << std::fixed << fn.getHtLsbGeV() << " GeV" << std::endl;
0325   os << "Central/Forward boundary     : " << std::fixed << fn.getCenForJetEtaBoundary() << std::endl;
0326 
0327   os << std::endl;
0328 
0329   os << std::setprecision(6);
0330   os << ios::scientific;
0331 
0332   os << "=== Level-1 GCT : Jet Et Calibration Function  ===" << std::endl;
0333   if (fn.getCorrType() == 0) {
0334     os << "No jet energy corrections applied" << std::endl;
0335   } else {
0336     switch (fn.getCorrType()) {
0337       case 1:
0338         os << "Function = Power series" << std::endl;
0339         break;
0340       case 2:
0341         os << "Function = ORCA" << std::endl;
0342         break;
0343       case 3:
0344         os << "Function = Simple" << std::endl;
0345         break;
0346       case 4:
0347         os << "Function = PiecewiseCubic" << std::endl;
0348         break;
0349       case 5:
0350         os << "Function = PF" << std::endl;
0351         break;
0352       default:
0353         os << "Unrecognised" << std::endl;
0354         break;
0355     }
0356     std::vector<std::vector<double> > jetCoeffs = fn.getJetCorrCoeffs();
0357     std::vector<std::vector<double> > tauCoeffs = fn.getTauCorrCoeffs();
0358 
0359     os << "Non-tau jet correction coefficients" << std::endl;
0360     for (unsigned i = 0; i < jetCoeffs.size(); i++) {
0361       os << "Eta =" << std::setw(2) << i;
0362       if (jetCoeffs.at(i).empty()) {
0363         os << ", no coefficients";
0364       } else {
0365         os << " Coefficients = ";
0366         for (unsigned j = 0; j < jetCoeffs.at(i).size(); j++) {
0367           os << jetCoeffs.at(i).at(j) << ", ";
0368         }
0369       }
0370       os << std::endl;
0371     }
0372     os << "Tau jet correction coefficients" << std::endl;
0373     for (unsigned i = 0; i < tauCoeffs.size(); i++) {
0374       os << "Eta =" << std::setw(2) << i;
0375       if (tauCoeffs.at(i).empty()) {
0376         os << ", no coefficients";
0377       } else {
0378         os << " Coefficients = ";
0379         for (unsigned j = 0; j < tauCoeffs.at(i).size(); j++) {
0380           os << tauCoeffs.at(i).at(j) << ", ";
0381         }
0382       }
0383       os << std::endl;
0384     }
0385   }
0386 
0387   os.unsetf(ios::fixed | ios::scientific);
0388 
0389   return os;
0390 }
[ 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!