Reference Manual

Go to the documentation of this file.
00001 #include "CondFormats/L1TObjects/interface/L1GctJetFinderParams.h"
00002 
00003 #include <iostream>
00004 #include <iomanip>
00005 #include <math.h>
00006 
00007 #include "FWCore/Utilities/interface/Exception.h"
00008 
00009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00010 #include "DataFormats/L1GlobalCaloTrigger/interface/L1GctStaticParameters.h"
00011 
00012 using std::ios;
00013 
00014 const unsigned L1GctJetFinderParams::NUMBER_ETA_VALUES = 11;
00015 const unsigned L1GctJetFinderParams::N_CENTRAL_ETA_VALUES = 7;
00016 
00017 L1GctJetFinderParams::L1GctJetFinderParams() :
00018   rgnEtLsb_(0.),
00019   htLsb_(0.),
00020   cenJetEtSeed_(0.),
00021   forJetEtSeed_(0.),
00022   tauJetEtSeed_(0.),
00023   tauIsoEtThreshold_(0.),
00024   htJetEtThreshold_(0.),
00025   mhtJetEtThreshold_(0.),
00026   cenForJetEtaBoundary_(0),
00027   corrType_(0),
00028   jetCorrCoeffs_(),
00029   tauCorrCoeffs_(),
00030   convertToEnergy_(false),
00031   energyConversionCoeffs_()
00032 { }
00033 
00034 L1GctJetFinderParams::L1GctJetFinderParams(double rgnEtLsb,
00035                                            double htLsb,
00036                                            double cJetSeed,
00037                                            double fJetSeed,
00038                                            double tJetSeed,
00039                                            double tauIsoEtThresh,
00040                                            double htJetEtThresh,
00041                                            double mhtJetEtThresh,
00042                                            unsigned etaBoundary,
00043                                            unsigned corrType,
00044                                            const std::vector< std::vector<double> >& jetCorrCoeffs,
00045                                            const std::vector< std::vector<double> >& tauCorrCoeffs,
00046                                            bool convertToEnergy,
00047                                            const std::vector<double>& energyConvCoeffs) :
00048   rgnEtLsb_(rgnEtLsb),
00049   htLsb_(htLsb),
00050   cenJetEtSeed_(cJetSeed),
00051   forJetEtSeed_(fJetSeed),
00052   tauJetEtSeed_(tJetSeed),
00053   tauIsoEtThreshold_(tauIsoEtThresh),
00054   htJetEtThreshold_(htJetEtThresh),
00055   mhtJetEtThreshold_(mhtJetEtThresh),
00056   cenForJetEtaBoundary_(etaBoundary),
00057   corrType_(corrType),
00058   jetCorrCoeffs_(jetCorrCoeffs),
00059   tauCorrCoeffs_(tauCorrCoeffs),
00060   convertToEnergy_(convertToEnergy),
00061   energyConversionCoeffs_(energyConvCoeffs)
00062 { 
00063   // check number of eta bins
00064   if (jetCorrCoeffs_.size() != NUMBER_ETA_VALUES ||
00065       tauCorrCoeffs_.size() != N_CENTRAL_ETA_VALUES ||
00066       energyConversionCoeffs_.size() != NUMBER_ETA_VALUES) {
00067 
00068     LogDebug("L1-O2O") << "GCT jet corrections constructed with " << jetCorrCoeffs_.size() << " bins, expected " << NUMBER_ETA_VALUES << std::endl;
00069     LogDebug("L1-O2O") << "GCT tau corrections constructed with " << tauCorrCoeffs_.size() << " bins, expected " << N_CENTRAL_ETA_VALUES << std::endl;
00070     LogDebug("L1-O2O") << "GCT energy corrections constructed with " << energyConversionCoeffs_.size() << " bins, expected " << NUMBER_ETA_VALUES << std::endl;
00071 
00072     throw cms::Exception("InconsistentConfig") << "L1GctJetFinderParams constructed with wrong number of eta bins : " << jetCorrCoeffs_.size() << " jets, " << tauCorrCoeffs_.size() << " taus, " << energyConversionCoeffs_.size() << " energy conversion bins" << std::endl;
00073 
00074   }
00075 
00076   // check number of coefficients against expectation
00077   unsigned expCoeffs = 0;
00078   if (corrType_ == 2) expCoeffs=8;  // ORCA style
00079   if (corrType_ == 3) expCoeffs=4;  // Simple
00080   if (corrType_ == 4) expCoeffs=15;  // piecewise-cubic
00081   if (corrType_ == 5) expCoeffs=6;  // PF
00082 
00083   // correction types 1 and 4 can have any number of parameters
00084   if (corrType_ != 1 &&
00085       corrType_ != 4) {
00086     std::vector< std::vector<double> >::const_iterator itr;      
00087     for (itr=jetCorrCoeffs_.begin(); itr!=jetCorrCoeffs_.end(); ++itr) {
00088       if (itr->size() != expCoeffs) {
00089         throw cms::Exception("InconsistentConfig") << "L1GctJetFinderParams constructed with " << itr->size() << " jet correction coefficients, when " << expCoeffs << " expected" << std::endl;
00090       }
00091     }
00092     for (itr=tauCorrCoeffs_.begin(); itr!=tauCorrCoeffs_.end(); ++itr) {
00093       if (itr->size() != expCoeffs) {
00094         throw cms::Exception("InconsistentConfig") << "L1GctJetFinderParams constructed with " << itr->size() << " tau correction coefficients, when " << expCoeffs << " expected"<< std::endl;
00095       }
00096     }
00097   }
00098   
00099 }
00100 
00101 
00102 L1GctJetFinderParams::~L1GctJetFinderParams() {}
00103 
00104 //---------------------------------------------------------------------------------------------
00105 //
00106 // set methods
00107 //
00108 
00109 void L1GctJetFinderParams::setRegionEtLsb (const double rgnEtLsb)
00110 {
00111   rgnEtLsb_ = rgnEtLsb;
00112 }
00113 
00114 void L1GctJetFinderParams::setSlidingWindowParams(const double cJetSeed,
00115                                                   const double fJetSeed,
00116                                                   const double tJetSeed,
00117                                                   const unsigned etaBoundary)
00118 {
00119   cenJetEtSeed_ = cJetSeed;
00120   forJetEtSeed_ = fJetSeed;
00121   tauJetEtSeed_ = tJetSeed;
00122   cenForJetEtaBoundary_ = etaBoundary;
00123 }
00124 
00125 void L1GctJetFinderParams::setJetEtCalibrationParams(const unsigned corrType,
00126                                                      const std::vector< std::vector<double> >& jetCorrCoeffs,
00127                                                      const std::vector< std::vector<double> >& tauCorrCoeffs)
00128 {
00129   corrType_ = corrType;
00130   jetCorrCoeffs_ = jetCorrCoeffs;
00131   tauCorrCoeffs_ = tauCorrCoeffs;
00132 }
00133 
00134 void L1GctJetFinderParams::setJetEtConvertToEnergyOn(const std::vector<double>& energyConvCoeffs)
00135 {
00136   convertToEnergy_ = true;
00137   energyConversionCoeffs_ = energyConvCoeffs;
00138 }
00139 
00140 void L1GctJetFinderParams::setJetEtConvertToEnergyOff()
00141 {
00142   convertToEnergy_ = false;
00143   energyConversionCoeffs_.clear();
00144 }
00145 
00146 void L1GctJetFinderParams::setHtSumParams(const double htLsb,
00147                                           const double htJetEtThresh,
00148                                           const double mhtJetEtThresh)
00149 {
00150   htLsb_ = htLsb;
00151   htJetEtThreshold_ = htJetEtThresh;
00152   mhtJetEtThreshold_ = mhtJetEtThresh;
00153 }
00154 
00155 void L1GctJetFinderParams::setTauAlgorithmParams(const double tauIsoEtThresh)
00156 {
00157   tauIsoEtThreshold_ = tauIsoEtThresh;
00158 }
00159 
00160 void L1GctJetFinderParams::setParams(const double rgnEtLsb,
00161                                      const double htLsb,
00162                                      const double cJetSeed,
00163                                      const double fJetSeed,
00164                                      const double tJetSeed,
00165                                      const double tauIsoEtThresh,
00166                                      const double htJetEtThresh,
00167                                      const double mhtJetEtThresh,
00168                                      const unsigned etaBoundary,
00169                                      const unsigned corrType,
00170                                      const std::vector< std::vector<double> >& jetCorrCoeffs,
00171                                      const std::vector< std::vector<double> >& tauCorrCoeffs)
00172 {
00173   setRegionEtLsb (rgnEtLsb);
00174   setSlidingWindowParams(cJetSeed, fJetSeed, tJetSeed, etaBoundary);
00175   setJetEtCalibrationParams(corrType, jetCorrCoeffs, tauCorrCoeffs);
00176   setHtSumParams(htLsb, htJetEtThresh, mhtJetEtThresh);
00177   setTauAlgorithmParams(tauIsoEtThresh);
00178 }
00179 
00180 //---------------------------------------------------------------------------------------------
00181 //
00182 // Jet Et correction methods
00183 //
00184 
00185 double L1GctJetFinderParams::correctedEtGeV(const double et, 
00186                                          const unsigned eta, 
00187                                          const bool tauVeto) const {
00188 
00189   if (eta>=NUMBER_ETA_VALUES) {
00190     return 0;
00191   } else {
00192     double result=0;
00193     if ((eta>=cenForJetEtaBoundary_) || tauVeto) {
00194       // Use jetCorrCoeffs for central and forward jets.
00195       // In forward eta bins we ignore the tau flag (as in the firmware)
00196       result=correctionFunction(et, jetCorrCoeffs_.at(eta));
00197     } else {
00198       // Use tauCorrCoeffs for tau jets (in central eta bins)
00199       result=correctionFunction(et, tauCorrCoeffs_.at(eta));
00200     }
00201     if (convertToEnergy_)  { result *= energyConversionCoeffs_.at(eta); }
00202     return result;
00203   }
00204 
00205 }
00206 
00207 
00209 uint16_t L1GctJetFinderParams::correctedEtGct(const double correctedEt) const
00210 {
00211   double scaledEt = correctedEt / htLsb_;
00212 
00213   uint16_t jetEtOut = static_cast<uint16_t>(scaledEt);
00214   
00215   if(jetEtOut > L1GctStaticParameters::jetCalibratedEtMax) {
00216     return L1GctStaticParameters::jetCalibratedEtMax;
00217   } else {
00218     return jetEtOut;
00219   }
00220 }
00221 
00222 
00223 
00224 // private methods
00225  
00226 double L1GctJetFinderParams::correctionFunction(const double Et, const std::vector<double>& coeffs) const
00227 {
00228   double result=0;
00229   switch (corrType_)
00230   {
00231   case 0:  // no correction
00232     result = Et;
00233     break;
00234   case 1:   // power series correction
00235     result = powerSeriesCorrect(Et, coeffs);
00236     break;
00237   case 2:  // ORCA style correction
00238     result = orcaStyleCorrect(Et, coeffs);
00239     break;
00240   case 3:  // simple correction (JetMEt style)
00241     result = simpleCorrect(Et, coeffs);
00242     break;
00243   case 4:  // piecwise cubic correction a la Greg Landsberg et al
00244     result = piecewiseCubicCorrect(Et, coeffs);
00245     break;
00246   case 5:  // PF correction
00247     result = pfCorrect(Et, coeffs);
00248     break;
00249   default:
00250     result = Et;      
00251   }
00252   return result;
00253 }
00254 
00255 double L1GctJetFinderParams::powerSeriesCorrect(const double Et, const std::vector<double>& coeffs) const
00256 {
00257   double corrEt = Et;
00258   for (unsigned i=0; i<coeffs.size();i++) {
00259     corrEt += coeffs.at(i)*pow(Et,(int)i); 
00260   }
00261   return corrEt;
00262 }
00263 
00264 double L1GctJetFinderParams::orcaStyleCorrect(const double Et, const std::vector<double>& coeffs) const
00265 {
00266 
00267   // The coefficients are arranged in groups of four
00268   // The first in each group is a threshold value of Et.
00269 
00270   std::vector<double>::const_iterator next_coeff=coeffs.begin();
00271   while (next_coeff != coeffs.end()) {
00272     double threshold = *next_coeff++;
00273     double A = *next_coeff++;
00274     double B = *next_coeff++;
00275     double C = *next_coeff++;
00276     if (Et>threshold) {
00277       // This function is an inverse quadratic:
00278       //   (input Et) = A + B*(output Et) + C*(output Et)^2
00279       return 2*(Et-A)/(B+sqrt(B*B-4*A*C+4*Et*C));
00280     }
00281     // If we are below all specified thresholds (or the vector is empty), return output=input.
00282   }
00283   return Et;
00284 }
00285 
00286 double L1GctJetFinderParams::simpleCorrect(const double Et, const std::vector<double>& coeffs) const
00287 {
00288   // function is :
00289   //    Et_out = A + B/[ (log10(Et_in))^C + D ] + E/Et_in
00290   // 
00291   //    fitcor_as_str = "[0]+[1]/(pow(log10(x),[2])+[3]) + [4]/x"
00292   // 
00293 
00294   return coeffs.at(0) + coeffs.at(1)/(pow(log10(Et),coeffs.at(2))+coeffs.at(3)) + (coeffs.at(4)/Et);
00295 
00296 }
00297 
00298 double L1GctJetFinderParams::piecewiseCubicCorrect(const double Et, const std::vector<double>& coeffs) const
00299 {
00300 
00301   // The correction fuction is a set of 3rd order polynomials
00302   //    Et_out = Et_in + (p0 + p1*Et_in + p2*Et_in^2 + p3*Et_in^3)
00303   // with different coefficients for different energy ranges.
00304   // The parameters are arranged in groups of five.
00305   // The first in each group is a threshold value of input Et,
00306   // followed by the four coefficients for the cubic function.
00307   double etOut = Et;
00308   std::vector<double>::const_iterator next_coeff=coeffs.begin();
00309   while (next_coeff != coeffs.end()) {
00310 
00311     // Read the coefficients from the vector
00312     double threshold = *next_coeff++;
00313     double A = *next_coeff++; //p0
00314     double B = *next_coeff++; //p1
00315     double C = *next_coeff++; //p2
00316     double D = *next_coeff++; //p3
00317 
00318     // Check we are in the right energy range and make correction
00319     if (Et>threshold) {
00320       etOut += (A + etOut*(B + etOut*(C + etOut*D))) ;
00321       break;
00322     }
00323 
00324   }
00325   return etOut;
00326 
00327 }
00328 
00329 double  L1GctJetFinderParams::pfCorrect(const double et, const std::vector<double>& coeffs) const
00330 {  
00331   //
00332   // corr_factor = [0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5])) 
00333   // Et_out = Et_in * corr_factor
00334   //
00335 
00336   return et * (coeffs.at(0)+coeffs.at(1)/(pow(log10(et),2)+coeffs.at(2))+coeffs.at(3)*exp(-coeffs.at(4)*(log10(et)-coeffs.at(5))*(log10(et)-coeffs.at(5))));
00337 
00338 }
00339 
00340 
00341 std::ostream& operator << (std::ostream& os, const L1GctJetFinderParams& fn)
00342 {
00343   //  os << std::setprecision(2);
00344 
00345   os << "=== Level-1 GCT : Jet Finder Parameters  ===" << std::endl;
00346   os << "RCT region LSB               : " << std::fixed << fn.getRgnEtLsbGeV() << " GeV" << std::endl;
00347   os << "Central jet seed threshold   : " << std::fixed << fn.getCenJetEtSeedGeV() << " GeV" << std::endl;
00348   os << "Tau jet seed threshold       : " << std::fixed << fn.getTauJetEtSeedGeV() << " GeV" << std::endl;
00349   os << "Forward jet seed threshold   : " << std::fixed << fn.getForJetEtSeedGeV() << " GeV" << std::endl;
00350   os << "Tau isolation threshold      : " << std::fixed << fn.getTauIsoEtThresholdGeV() << " GeV" << std::endl;
00351   os << "Ht jet Et threshold          : " << std::fixed << fn.getHtJetEtThresholdGeV() << " GeV" << std::endl;
00352   os << "MHt jet Et threshold         : " << std::fixed << fn.getMHtJetEtThresholdGeV() << " GeV" << std::endl;
00353   os << "Ht LSB                       : " << std::fixed << fn.getHtLsbGeV() << " GeV" << std::endl;
00354   os << "Central/Forward boundary     : " << std::fixed << fn.getCenForJetEtaBoundary() << std::endl;
00355 
00356   os << std::endl;
00357 
00358   os << std::setprecision(6);
00359   os << ios::scientific;
00360 
00361   os << "=== Level-1 GCT : Jet Et Calibration Function  ===" << std::endl;
00362   if (fn.getCorrType() == 0) {
00363     os << "No jet energy corrections applied" << std::endl;
00364   } else { 
00365     switch (fn.getCorrType())
00366     {
00367       case 1:
00368         os << "Function = Power series" << std::endl;
00369         break;
00370       case 2:
00371         os << "Function = ORCA" << std::endl;
00372         break;
00373       case 3:
00374         os << "Function = Simple" << std::endl;
00375         break;
00376       case 4:
00377         os << "Function = PiecewiseCubic" << std::endl;
00378         break;
00379       case 5:
00380         os << "Function = PF" << std::endl;
00381         break;
00382       default:
00383         os << "Unrecognised" << std::endl;
00384         break; 
00385     }
00386     std::vector< std::vector<double> > jetCoeffs = fn.getJetCorrCoeffs();
00387     std::vector< std::vector<double> > tauCoeffs = fn.getTauCorrCoeffs();
00388 
00389     os << "Non-tau jet correction coefficients" << std::endl;
00390     for (unsigned i=0; i<jetCoeffs.size(); i++){
00391       os << "Eta =" << std::setw(2) << i;
00392       if (jetCoeffs.at(i).empty()) {
00393         os << ", no coefficients";
00394       } else {
00395         os << " Coefficients = ";
00396         for (unsigned j=0; j<jetCoeffs.at(i).size();j++){
00397           os << jetCoeffs.at(i).at(j) << ", "; 
00398         }
00399       }
00400       os << std::endl;
00401     }
00402     os << "Tau jet correction coefficients" << std::endl;
00403     for (unsigned i=0; i<tauCoeffs.size(); i++){
00404       os << "Eta =" << std::setw(2) << i;
00405       if (tauCoeffs.at(i).empty()) {
00406         os << ", no coefficients";
00407       } else {
00408         os << " Coefficients = ";
00409         for (unsigned j=0; j<tauCoeffs.at(i).size();j++){
00410           os << tauCoeffs.at(i).at(j) << ", "; 
00411         }
00412       }
00413       os << std::endl;
00414     }
00415   }
00416 
00417   os.unsetf(ios::fixed | ios::scientific);
00418 
00419   return os;
00420 }
00421 
00422
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/CondFormats/L1TObjects/src/L1GctJetFinderParams.cc
