Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:10:19 

0001 #include "EgammaAnalysis/ElectronTools/interface/EnergyScaleCorrection_class.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include <RooDataSet.h>
0004 #include <RooArgSet.h>
0005 
0006 #define NGENMAX 10000
0007 
0008 // for exit(0)
0009 #include <cstdlib>
0010 #include <cfloat>
0011 // for setw
0012 #include <iomanip>
0013 //for istreamstring
0014 #include <sstream>
0015 
0016 //#define DEBUG
0017 //#define PEDANTIC_OUTPUT
0018 
0019 EnergyScaleCorrection_class::EnergyScaleCorrection_class(std::string correctionFileName, unsigned int genSeed)
0020     : doScale(false), doSmearings(false), smearingType_(ECALELF) {
0021   if (!correctionFileName.empty()) {
0022     std::string filename = correctionFileName + "_scales.dat";
0023     ReadFromFile(filename);
0024     if (scales.empty()) {
0025       std::cerr << "[ERROR] scale correction map empty" << std::endl;
0026       exit(1);
0027     }
0028   }
0029 
0030   if (!correctionFileName.empty()) {
0031     std::string filename = correctionFileName + "_smearings.dat";
0032     ReadSmearingFromFile(filename);
0033     if (smearings.empty()) {
0034       std::cerr << "[ERROR] smearing correction map empty" << std::endl;
0035       exit(1);
0036     }
0037   }
0038 
0039   return;
0040 }
0041 
0042 EnergyScaleCorrection_class::~EnergyScaleCorrection_class(void) { return; }
0043 
0044 float EnergyScaleCorrection_class::ScaleCorrection(
0045     unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const {
0046   float correction = 1;
0047   if (doScale)
0048     correction *= getScaleOffset(runNumber, isEBEle, R9Ele, etaSCEle, EtEle);
0049 
0050   return correction;
0051 }
0052 
0053 float EnergyScaleCorrection_class::ScaleCorrectionUncertainty(
0054     unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const {
0055   float statUncert = getScaleStatUncertainty(runNumber, isEBEle, R9Ele, etaSCEle, EtEle);
0056   float systUncert = getScaleSystUncertainty(runNumber, isEBEle, R9Ele, etaSCEle, EtEle);
0057 
0058   return sqrt(statUncert * statUncert + systUncert * systUncert);
0059 }
0060 
0061 float EnergyScaleCorrection_class::getScaleStatUncertainty(
0062     unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const {
0063   return getScaleCorrection(runNumber, isEBEle, R9Ele, etaSCEle, EtEle).scale_err;
0064 }
0065 
0066 float EnergyScaleCorrection_class::getScaleSystUncertainty(
0067     unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const {
0068   return getScaleCorrection(runNumber, isEBEle, R9Ele, etaSCEle, EtEle).scale_err_syst;
0069 }
0070 
0071 correctionValue_class EnergyScaleCorrection_class::getScaleCorrection(
0072     unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const {
0073   // buld the category based on the values of the object
0074   correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
0075   correction_map_t::const_iterator corr_itr =
0076       scales.find(category);  // find the correction value in the map that associates the category to the correction
0077 
0078   if (corr_itr == scales.end()) {  // if not in the standard classes, add it in the list of not defined classes
0079 
0080     // this part is commented because it makes the method not constant
0081     // if(scales_not_defined.count(category) == 0) {
0082     //  correctionValue_class corr;
0083     //  scales_not_defined[category] = corr;
0084     // }
0085     // corr_itr = scales_not_defined.find(category);
0086     /// \todo this can be switched to an exeption
0087     std::cout << "[ERROR] Scale category not found: " << std::endl;
0088     std::cout << category << std::endl;
0089     std::cout << "Returning uncorrected value." << std::endl;
0090     //     exit(1);
0091     correctionValue_class nocorr;
0092     std::cout << nocorr << std::endl;
0093     return nocorr;
0094   }
0095 
0096 #ifdef DEBUG
0097   std::cout << "[DEBUG] Checking scale correction for category: " << category << std::endl;
0098   std::cout << "[DEBUG] Correction is: " << corr_itr->second << std::endl
0099             << "        given for category " << corr_itr->first << std::endl;
0100   ;
0101 #endif
0102   return corr_itr->second;
0103 }
0104 
0105 float EnergyScaleCorrection_class::getScaleOffset(
0106     unsigned int runNumber, bool isEBEle, double R9Ele, double etaSCEle, double EtEle) const {
0107   // buld the category based on the values of the object
0108   correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
0109   correction_map_t::const_iterator corr_itr =
0110       scales.find(category);  // find the correction value in the map that associates the category to the correction
0111 
0112   if (corr_itr == scales.end()) {  // if not in the standard classes, add it in the list of not defined classes
0113 
0114     // this part is commented because it makes the method not constant
0115     // if(scales_not_defined.count(category) == 0) {
0116     //  correctionValue_class corr;
0117     //  scales_not_defined[category] = corr;
0118     // }
0119     // corr_itr = scales_not_defined.find(category);
0120     /// \todo this can be switched to an exeption
0121     std::cout << "[ERROR] Scale offset category not found: " << std::endl;
0122     std::cout << category << std::endl;
0123     std::cout << "Returning uncorrected value." << std::endl;
0124     //     exit(1);
0125     correctionValue_class nocorr;
0126     return nocorr.scale;
0127   }
0128 
0129 #ifdef DEBUG
0130   std::cout << "[DEBUG] Checking scale offset correction for category: " << category << std::endl;
0131   std::cout << "[DEBUG] Correction is: " << corr_itr->second << std::endl
0132             << "        given for category " << corr_itr->first << std::endl;
0133   ;
0134 #endif
0135 
0136   return corr_itr->second.scale;
0137 }
0138 
0139 /**
0140  *  Input file structure:
0141  *  category  "runNumber"   runMin  runMax   deltaP  err_deltaP(stat on single bins)  err_deltaP_stat(to be used) err_deltaP_syst(to be used)
0142  *
0143  */
0144 void EnergyScaleCorrection_class::ReadFromFile(TString filename) {
0145 #ifdef PEDANTIC_OUTPUT
0146   std::cout << "[STATUS] Reading energy scale correction  values from file: " << filename << std::endl;
0147 #endif
0148   std::cout << "[STATUS] Reading energy scale correction  values from file: " << filename << std::endl;
0149 
0150   //std::ifstream Ccufile(edm::FileInPath(Ccufilename).fullPath().c_str(),std::ios::in);
0151   std::ifstream f_in(edm::FileInPath(filename).fullPath().c_str());
0152 
0153   if (!f_in.good()) {
0154     std::cerr << "[ERROR] file " << filename << " not readable" << std::endl;
0155     exit(1);
0156     return;
0157   }
0158 
0159   int runMin, runMax;
0160   TString category, region2;
0161   double deltaP, err_deltaP, err_deltaP_stat, err_deltaP_syst;
0162 
0163   for (f_in >> category; f_in.good(); f_in >> category) {
0164     f_in >> region2 >> runMin >> runMax >> deltaP >> err_deltaP >> err_deltaP_stat >> err_deltaP_syst;
0165 
0166     AddScale(category, runMin, runMax, deltaP, err_deltaP_stat, err_deltaP_syst);
0167   }
0168 
0169   f_in.close();
0170 
0171   return;
0172 }
0173 
0174 // this method adds the correction values read from the txt file to the map
0175 void EnergyScaleCorrection_class::AddScale(
0176     TString category_, int runMin_, int runMax_, double deltaP_, double err_deltaP_, double err_syst_deltaP) {
0177   correctionCategory_class cat(category_);  // build the category from the string
0178   cat.runmin = runMin_;
0179   cat.runmax = runMax_;
0180 
0181   // the following check is not needed, can be removed
0182   if (scales.count(cat) != 0) {
0183     std::cerr << "[ERROR] Category already defined!" << std::endl;
0184     std::cerr << "        Adding category:  " << cat << std::endl;
0185     std::cerr << "        Defined category: " << scales[cat] << std::endl;
0186     exit(1);
0187   }
0188 
0189   correctionValue_class corr;  // define the correction values
0190   corr.scale = deltaP_;
0191   corr.scale_err = err_deltaP_;
0192   corr.scale_err_syst = err_syst_deltaP;
0193   scales[cat] = corr;
0194 
0195 #ifdef PEDANTIC_OUTPUT
0196   std::cout << "[INFO:scale correction] " << cat << corr << std::endl;
0197 #endif
0198   return;
0199 }
0200 
0201 //============================== SMEARING
0202 void EnergyScaleCorrection_class::AddSmearing(TString category_,
0203                                               int runMin_,
0204                                               int runMax_,
0205                                               double rho,
0206                                               double err_rho,
0207                                               double phi,
0208                                               double err_phi,
0209                                               double Emean,
0210                                               double err_Emean) {
0211   correctionCategory_class cat(category_);
0212   cat.runmin = (runMin_ < 0) ? 0 : runMin_;
0213   cat.runmax = runMax_;
0214 
0215   if (smearings.count(cat) != 0) {
0216     std::cerr << "[ERROR] Smearing category already defined!" << std::endl;
0217     std::cerr << "        Adding category:  " << cat << std::endl;
0218     std::cerr << "        Defined category: " << smearings[cat] << std::endl;
0219     exit(1);
0220   }
0221 
0222   correctionValue_class corr;
0223   corr.rho = rho;
0224   corr.rho_err = err_rho;
0225   corr.phi = phi;
0226   corr.phi_err = err_phi;
0227   corr.Emean = Emean;
0228   corr.Emean_err = err_Emean;
0229   smearings[cat] = corr;
0230 
0231 #ifdef PEDANTIC_OUTPUT
0232 #ifndef CMSSW
0233   std::cout << "[INFO:smearings] " << cat << corr << std::endl;
0234 #else
0235   edm::LogInfo("[INFO:smearings] ") << cat << corr;
0236 #endif
0237 #endif
0238 
0239   return;
0240 }
0241 
0242 /**
0243  *  File structure:
0244  EBlowEtaBad8TeV    0 0.0 1.0 -999. 0.94 -999999 999999 6.73 0. 7.7e-3  6.32e-4 0.00 0.16
0245  EBlowEtaGold8TeV   0 0.0 1.0 0.94  999. -999999 999999 6.60 0. 7.4e-3  6.50e-4 0.00 0.16
0246  EBhighEtaBad8TeV   0 1.0 1.5 -999. 0.94 -999999 999999 6.73 0. 1.26e-2 1.03e-3 0.00 0.07
0247  EBhighEtaGold8TeV  0 1.0 1.5 0.94  999. -999999 999999 6.52 0. 1.12e-2 1.32e-3 0.00 0.22
0248  ##################################################################################################
0249  EElowEtaBad8TeV    0 1.5 2.0 -999. 0.94 -999999 999999 0.   0. 1.98e-2 3.03e-3 0.  0.
0250  EElowEtaGold8TeV   0 1.5 2.0 0.94  999. -999999 999999 0.   0. 1.63e-2 1.22e-3 0.  0.
0251  EEhighEtaBad8TeV   0 2.0 3.0 -999. 0.94 -999999 999999 0.   0. 1.92e-2 9.22e-4 0.  0.
0252  EEhighEtaGold8TeV  0 2.0 3.0 0.94  999. -999999 999999 0.   0. 1.86e-2 7.81e-4 0.  0.
0253  ##################################################################################################
0254  *
0255  */
0256 
0257 void EnergyScaleCorrection_class::ReadSmearingFromFile(TString filename) {
0258 #ifdef PEDANTIC_OUTPUT
0259   std::cout << "[STATUS] Reading smearing values from file: " << filename << std::endl;
0260 #endif
0261   //edm::FileInPath(Ccufilename).fullPath().c_str(),std::ios::in); .fullPath().c_str()
0262   std::ifstream f_in(edm::FileInPath(filename).fullPath().c_str());
0263   if (!f_in.good()) {
0264     std::cerr << "[ERROR] file " << filename << " not readable" << std::endl;
0265     exit(1);
0266   }
0267 
0268   int runMin = 0, runMax = 900000;
0269   int unused = 0;
0270   TString category, region2;
0271   //double smearing, err_smearing;
0272   double rho, phi, Emean, err_rho, err_phi, err_Emean;
0273   double etaMin, etaMax, r9Min, r9Max;
0274   std::string phi_string, err_phi_string;
0275 
0276   while (f_in.peek() != EOF && f_in.good()) {
0277     if (f_in.peek() == 10) {  // 10 = \n
0278       f_in.get();
0279       continue;
0280     }
0281 
0282     if (f_in.peek() == 35) {  // 35 = #
0283       f_in.ignore(1000, 10);  // ignore the rest of the line until \n
0284       continue;
0285     }
0286 
0287     if (smearingType_ == UNKNOWN) {  // trying to guess: not recommended
0288       std::cerr << "[ERROR] Not implemented" << std::endl;
0289       assert(false);
0290 
0291     } else if (smearingType_ == GLOBE) {
0292       f_in >> category >> unused >> etaMin >> etaMax >> r9Min >> r9Max >> runMin >> runMax >> Emean >> err_Emean >>
0293           rho >> err_rho >> phi >> err_phi;
0294 
0295       AddSmearing(category, runMin, runMax, rho, err_rho, phi, err_phi, Emean, err_Emean);
0296 
0297     } else if (smearingType_ == ECALELF) {
0298       f_in >> category >> Emean >> err_Emean >> rho >> err_rho >> phi_string >> err_phi_string;
0299 #ifdef DEBUG
0300       std::cout
0301           << category
0302           //<< "\t" << etaMin << "\t" << etaMax << "\t" << r9Min << "\t" << r9Max << "\t" << runMin << "\t" << runMax
0303           << "\tEmean=" << Emean << "\t" << rho << "\t" << err_rho << "\tphi_string=" << phi_string
0304           << "#\terr_phi_string=" << err_phi_string << std::endl;
0305 #endif
0306 
0307       if (phi_string == "M_PI_2")
0308         phi = M_PI_2;
0309       else
0310         phi = std::stod(phi_string);
0311 
0312       if (err_phi_string == "M_PI_2")
0313         err_phi = M_PI_2;
0314       else
0315         err_phi = std::stod(err_phi_string);
0316 
0317       AddSmearing(category, runMin, runMax, rho, err_rho, phi, err_phi, Emean, err_Emean);
0318 
0319     } else {
0320       f_in >> category >> rho >> phi;
0321       err_rho = err_phi = Emean = err_Emean = 0;
0322       AddSmearing(category, runMin, runMax, rho, err_rho, phi, err_phi, Emean, err_Emean);
0323     }
0324 #ifdef DEBUG
0325     std::cout << category << "\t" << etaMin << "\t" << etaMax << "\t" << r9Min << "\t" << r9Max << "\t" << runMin
0326               << "\t" << runMax << "\tEmean=" << Emean << "\t" << rho << "\t" << phi << std::endl;
0327 #endif
0328   }
0329 
0330   f_in.close();
0331   //  runCorrection_itr=runMin_map.begin();
0332 
0333   return;
0334 }
0335 
0336 float EnergyScaleCorrection_class::getSmearingSigma(
0337     int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle, paramSmear_t par, float nSigma) const {
0338   if (par == kRho)
0339     return getSmearingSigma(runNumber, isEBEle, R9Ele, etaSCEle, EtEle, nSigma, 0.);
0340   if (par == kPhi)
0341     return getSmearingSigma(runNumber, isEBEle, R9Ele, etaSCEle, EtEle, 0., nSigma);
0342   return getSmearingSigma(runNumber, isEBEle, R9Ele, etaSCEle, EtEle, 0., 0.);
0343 }
0344 
0345 float EnergyScaleCorrection_class::getSmearingSigma(
0346     int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle, float nSigma_rho, float nSigma_phi) const {
0347   correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
0348   correction_map_t::const_iterator corr_itr = smearings.find(category);
0349   if (corr_itr == smearings.end()) {  // if not in the standard classes, add it in the list of not defined classes
0350     // the following commented part makes the method non const
0351     // if(smearings_not_defined.count(category) == 0) {
0352     //  correctionValue_class corr;
0353     //  smearings_not_defined[category] = corr;
0354     // }
0355     corr_itr = smearings_not_defined.find(category);
0356     std::cerr << "[WARNING] Smearing category not found: " << std::endl;
0357     std::cerr << category << std::endl;
0358     //     exit(1);
0359   }
0360 
0361 #ifdef DEBUG
0362   std::cout << "[DEBUG] Checking smearing correction for category: " << category << std::endl;
0363   std::cout << "[DEBUG] Correction is: " << corr_itr->second << std::endl
0364             << "        given for category " << corr_itr->first;
0365 #endif
0366 
0367   double rho = corr_itr->second.rho + corr_itr->second.rho_err * nSigma_rho;
0368   double phi = corr_itr->second.phi + corr_itr->second.phi_err * nSigma_phi;
0369 
0370   double constTerm = rho * sin(phi);
0371   double alpha = rho * corr_itr->second.Emean * cos(phi);
0372 
0373   return sqrt(constTerm * constTerm + alpha * alpha / EtEle);
0374 }
0375 
0376 float EnergyScaleCorrection_class::getSmearingRho(
0377     int runNumber, bool isEBEle, float R9Ele, float etaSCEle, float EtEle) const {
0378   correctionCategory_class category(runNumber, etaSCEle, R9Ele, EtEle);
0379   correction_map_t::const_iterator corr_itr = smearings.find(category);
0380   if (corr_itr == smearings.end()) {  // if not in the standard classes, add it in the list of not defined classes
0381     // if(smearings_not_defined.count(category) == 0) {
0382     //  correctionValue_class corr;
0383     //  smearings_not_defined[category] = corr;
0384     // }
0385     corr_itr = smearings_not_defined.find(category);
0386   }
0387 
0388   return corr_itr->second.rho;
0389 }
0390 
0391 bool correctionCategory_class::operator<(const correctionCategory_class& b) const {
0392   if (runmin < b.runmin && runmax < b.runmax)
0393     return true;
0394   if (runmax > b.runmax && runmin > b.runmin)
0395     return false;
0396 
0397   if (etamin < b.etamin && etamax < b.etamax)
0398     return true;
0399   if (etamax > b.etamax && etamin > b.etamin)
0400     return false;
0401 
0402   if (r9min < b.r9min && r9max < b.r9max)
0403     return true;
0404   if (r9max > b.r9max && r9min > b.r9min)
0405     return false;
0406 
0407   if (etmin < b.etmin && etmax < b.etmax)
0408     return true;
0409   if (etmax > b.etmax && etmin > b.etmin)
0410     return false;
0411   return false;
0412 }
0413 
0414 correctionCategory_class::correctionCategory_class(TString category_) {
0415   std::string category(category_.Data());
0416 #ifdef DEBUG
0417   std::cout << "[DEBUG] correctionClass defined for category: " << category << std::endl;
0418 #endif
0419   // default values (corresponding to an existing category -- the worst one)
0420   runmin = 0;
0421   runmax = 999999;
0422   etamin = 2;
0423   etamax = 7;
0424   r9min = -1;
0425   r9max = 0.94;
0426   etmin = -1;
0427   etmax = 99999.;
0428 
0429   size_t p1, p2;  // boundary
0430   // eta region
0431   p1 = category.find("absEta_");
0432   if (category.find("absEta_0_1") != std::string::npos) {
0433     etamin = 0;
0434     etamax = 1;
0435   } else if (category.find("absEta_1_1.4442") != std::string::npos) {
0436     etamin = 1;
0437     etamax = 1.479;
0438   }  //etamax=1.4442; }
0439   else if (category.find("absEta_1.566_2") != std::string::npos) {
0440     etamin = 1.479;
0441     etamax = 2;
0442   }  //etamin=1.566; etamax=2;}
0443   else if (category.find("absEta_2_2.5") != std::string::npos) {
0444     etamin = 2;
0445     etamax = 3;
0446   } else {
0447     if (p1 != std::string::npos) {
0448       p1 = category.find('_', p1);
0449       p2 = category.find('_', p1 + 1);
0450       etamin = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
0451       p1 = p2;
0452       p2 = category.find('-', p1);
0453       etamax = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
0454     }
0455   }
0456 
0457   if (category.find("EBlowEta") != std::string::npos) {
0458     etamin = 0;
0459     etamax = 1;
0460   };
0461   if (category.find("EBhighEta") != std::string::npos) {
0462     etamin = 1;
0463     etamax = 1.479;
0464   };
0465   if (category.find("EElowEta") != std::string::npos) {
0466     etamin = 1.479;
0467     etamax = 2;
0468   };
0469   if (category.find("EEhighEta") != std::string::npos) {
0470     etamin = 2;
0471     etamax = 7;
0472   };
0473 
0474   // Et region
0475   p1 = category.find("-Et_");
0476   //  p2 = p1 + 1;//this is needed only for the printouts
0477   //    std::cout << "p1 = " << p1 << "\t" << std::string::npos << "\t" << category.size() << std::endl;
0478   //std::cout << etmin << "\t" << etmax << "\t" << category.substr(p1+1, p2-p1-1) << "\t" << p1 << "\t" << p2 << std::endl;
0479   if (p1 != std::string::npos) {
0480     p1 = category.find('_', p1);
0481     p2 = category.find('_', p1 + 1);
0482     etmin = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
0483     p1 = p2;
0484     p2 = category.find('-', p1);
0485     etmax = TString(category.substr(p1 + 1, p2 - p1 - 1)).Atof();
0486     //  std::cout << etmin << "\t" << etmax << "\t" << category.substr(p1 + 1, p2 - p1 - 1) << std::endl;
0487   }
0488 
0489   if (category.find("gold") != std::string::npos || category.find("Gold") != std::string::npos ||
0490       category.find("highR9") != std::string::npos) {
0491     r9min = 0.94;
0492     r9max = FLT_MAX;
0493   } else if (category.find("bad") != std::string::npos || category.find("Bad") != std::string::npos ||
0494              category.find("lowR9") != std::string::npos) {
0495     r9min = -1;
0496     r9max = 0.94;
0497   };
0498 }
[ 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!