File indexing completed on 2024-04-06 12:10:18
0001 #include "EgammaAnalysis/ElectronTools/interface/ElectronEnergyCalibrator.h"
0002
0003 #include <CLHEP/Random/RandGaussQ.h>
0004 #include <CLHEP/Random/RandFlat.h>
0005 #include <CLHEP/Random/Random.h>
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009
0010
0011
0012
0013
0014
0015
0016 #include <string>
0017 #include <vector>
0018 #include <fstream>
0019 #include <sstream>
0020 #include <iostream>
0021
0022 void ElectronEnergyCalibrator::init() {
0023 if (!isMC_)
0024 {
0025 if (verbose_) {
0026 std::cout << "[ElectronEnergyCalibrator] Initialization in DATA mode" << std::endl;
0027 }
0028
0029 std::ifstream fin(pathData_.c_str());
0030
0031 if (!fin) {
0032 throw cms::Exception("Configuration") << "[ElectronEnergyCalibrator] Cannot open the file " << pathData_
0033 << "\n It is not found, missed or corrupted";
0034 } else {
0035 if (verbose_) {
0036 std::cout << "[ElectronEnergyCalibrator] File " << pathData_ << " succesfully opened" << std::endl;
0037 }
0038
0039 std::string s;
0040 std::vector<std::string> selements;
0041 std::string delimiter = ",";
0042 nCorrValRaw = 0;
0043
0044 while (!fin.eof()) {
0045 getline(fin, s);
0046 if (!s.empty()) {
0047 splitString(s, selements, delimiter);
0048 corrValArray[nCorrValRaw].nRunMin = stringToDouble(selements[0]);
0049 corrValArray[nCorrValRaw].nRunMax = stringToDouble(selements[1]);
0050 corrValArray[nCorrValRaw].corrCat0 = stringToDouble(selements[2]);
0051 corrValArray[nCorrValRaw].corrCat1 = stringToDouble(selements[3]);
0052 corrValArray[nCorrValRaw].corrCat2 = stringToDouble(selements[4]);
0053 corrValArray[nCorrValRaw].corrCat3 = stringToDouble(selements[5]);
0054 corrValArray[nCorrValRaw].corrCat4 = stringToDouble(selements[6]);
0055 corrValArray[nCorrValRaw].corrCat5 = stringToDouble(selements[7]);
0056 corrValArray[nCorrValRaw].corrCat6 = stringToDouble(selements[8]);
0057 corrValArray[nCorrValRaw].corrCat7 = stringToDouble(selements[9]);
0058
0059 nCorrValRaw++;
0060
0061 selements.clear();
0062 }
0063 }
0064
0065 fin.close();
0066
0067 if (verbose_) {
0068 std::cout << "[ElectronEnergyCalibrator] File closed" << std::endl;
0069 }
0070 }
0071
0072 if (applyLinearityCorrection_) {
0073 std::ifstream finlin(pathLinData_.c_str());
0074
0075 if (!finlin) {
0076 throw cms::Exception("Configuration") << "[ElectronEnergyCalibrator] Cannot open the file " << pathLinData_
0077 << "\n It is not found, missed or corrupted";
0078 } else {
0079 if (verbose_) {
0080 std::cout << "[ElectronEnergyCalibrator] File with Linearity Corrections " << pathLinData_
0081 << " succesfully opened" << std::endl;
0082 }
0083
0084 std::string s;
0085 std::vector<std::string> selements;
0086 std::string delimiter = ",";
0087 nLinCorrValRaw = 0;
0088
0089 while (!finlin.eof()) {
0090 getline(finlin, s);
0091 if (!s.empty()) {
0092 splitString(s, selements, delimiter);
0093
0094 linCorrValArray[nLinCorrValRaw].ptMin = stringToDouble(selements[0]);
0095 linCorrValArray[nLinCorrValRaw].ptMax = stringToDouble(selements[1]);
0096 linCorrValArray[nLinCorrValRaw].corrCat0 = stringToDouble(selements[2]);
0097 linCorrValArray[nLinCorrValRaw].corrCat1 = stringToDouble(selements[3]);
0098 linCorrValArray[nLinCorrValRaw].corrCat2 = stringToDouble(selements[4]);
0099 linCorrValArray[nLinCorrValRaw].corrCat3 = stringToDouble(selements[5]);
0100 linCorrValArray[nLinCorrValRaw].corrCat4 = stringToDouble(selements[6]);
0101 linCorrValArray[nLinCorrValRaw].corrCat5 = stringToDouble(selements[7]);
0102
0103 nLinCorrValRaw++;
0104
0105 selements.clear();
0106 }
0107 }
0108
0109 finlin.close();
0110
0111 if (verbose_) {
0112 std::cout << "[ElectronEnergyCalibrator] File closed" << std::endl;
0113 }
0114 }
0115 }
0116 } else
0117 {
0118 if (verbose_) {
0119 std::cout << "[ElectronEnergyCalibrator] Initialization in MC mode" << std::endl;
0120 }
0121 }
0122 }
0123
0124 void ElectronEnergyCalibrator::splitString(const std::string &fullstr,
0125 std::vector<std::string> &elements,
0126 const std::string &delimiter) {
0127 std::string::size_type lastpos = fullstr.find_first_not_of(delimiter, 0);
0128 std::string::size_type pos = fullstr.find_first_of(delimiter, lastpos);
0129
0130 while ((std::string::npos != pos) || (std::string::npos != lastpos)) {
0131 elements.push_back(fullstr.substr(lastpos, pos - lastpos));
0132 lastpos = fullstr.find_first_not_of(delimiter, pos);
0133 pos = fullstr.find_first_of(delimiter, lastpos);
0134 }
0135 }
0136
0137 double ElectronEnergyCalibrator::stringToDouble(const std::string &str) {
0138 std::istringstream stm;
0139 double val = 0;
0140 stm.str(str);
0141 stm >> val;
0142 return val;
0143 }
0144
0145 void ElectronEnergyCalibrator::calibrate(SimpleElectron &electron, edm::StreamID const &streamID) {
0146 double scale = 1.0;
0147 double dsigMC = 0.;
0148 double corrMC = 0.;
0149 double run_ = electron.getRunNumber();
0150 bool isEB = electron.isEB();
0151 double eta = electron.getEta();
0152 double r9 = electron.getR9();
0153
0154 switch (correctionsType_) {
0155 case 1:
0156 if (verbose_) {
0157 std::cout << "[ElectronEnergyCalibrator] Using regression energy for calibration" << std::endl;
0158 }
0159 newEnergy_ = electron.getRegEnergy();
0160 newEnergyError_ = electron.getRegEnergyError();
0161 break;
0162 case 2:
0163 if (verbose_) {
0164 std::cout << "[ElectronEnergyCalibrator] Using scale corrections for new regression" << std::endl;
0165 }
0166 newEnergy_ = electron.getRegEnergy();
0167 newEnergyError_ = electron.getRegEnergyError();
0168 break;
0169 case 3:
0170 if (verbose_) {
0171 std::cout << "[ElectronEnergyCalibrator] Using standard ecal energy for calibration" << std::endl;
0172 }
0173 newEnergy_ = electron.getSCEnergy();
0174 newEnergyError_ = electron.getSCEnergyError();
0175 break;
0176 }
0177
0178 edm::Service<edm::RandomNumberGenerator> rng;
0179 if (!rng.isAvailable()) {
0180 throw cms::Exception("Configuration")
0181 << "XXXXXXX requires the RandomNumberGeneratorService\n"
0182 "which is not present in the configuration file. You must add the service\n"
0183 "in the configuration file or remove the modules that require it.";
0184 }
0185
0186 if (!isMC_) {
0187 for (int i = 0; i < nCorrValRaw; i++) {
0188 if ((run_ >= corrValArray[i].nRunMin) && (run_ <= corrValArray[i].nRunMax)) {
0189 if (isEB) {
0190 if (fabs(eta) < 1) {
0191 if (r9 < 0.94) {
0192 scale = corrValArray[i].corrCat0;
0193 } else {
0194 scale = corrValArray[i].corrCat1;
0195 }
0196 } else {
0197 if (r9 < 0.94) {
0198 scale = corrValArray[i].corrCat2;
0199 } else {
0200 scale = corrValArray[i].corrCat3;
0201 }
0202 }
0203 } else {
0204 if (fabs(eta) < 2) {
0205 if (r9 < 0.94) {
0206 scale = corrValArray[i].corrCat4;
0207 } else {
0208 scale = corrValArray[i].corrCat5;
0209 }
0210 } else {
0211 if (r9 < 0.94) {
0212 scale = corrValArray[i].corrCat6;
0213 } else {
0214 scale = corrValArray[i].corrCat7;
0215 }
0216 }
0217 }
0218 }
0219 }
0220 newEnergy_ = newEnergy_ * scale;
0221 }
0222
0223 switch (correctionsType_) {
0224 case 1:
0225
0226 if (dataset_ == "Summer12_DR53X_HCP2012" || dataset_ == "Moriond2013") {
0227 if (!isMC_) {
0228 if (run_ <= 203002) {
0229 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0230 dsigMC = 0.0103;
0231 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0232 dsigMC = 0.0090;
0233 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0234 dsigMC = 0.0190;
0235 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0236 dsigMC = 0.0156;
0237 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0238 dsigMC = 0.0269;
0239 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0240 dsigMC = 0.0287;
0241 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0242 dsigMC = 0.0364;
0243 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0244 dsigMC = 0.0321;
0245 } else {
0246 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0247 dsigMC = 0.0109;
0248 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0249 dsigMC = 0.0099;
0250 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0251 dsigMC = 0.0182;
0252 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0253 dsigMC = 0.0200;
0254 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0255 dsigMC = 0.0282;
0256 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0257 dsigMC = 0.0309;
0258 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0259 dsigMC = 0.0386;
0260 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0261 dsigMC = 0.0359;
0262 }
0263 } else {
0264 CLHEP::RandFlat flatRandom(rng->getEngine(streamID));
0265 double rn = flatRandom.fire();
0266 if (rn > lumiRatio_) {
0267 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0268 dsigMC = 0.0109;
0269 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0270 dsigMC = 0.0099;
0271 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0272 dsigMC = 0.0182;
0273 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0274 dsigMC = 0.0200;
0275 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0276 dsigMC = 0.0282;
0277 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0278 dsigMC = 0.0309;
0279 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0280 dsigMC = 0.0386;
0281 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0282 dsigMC = 0.0359;
0283 } else {
0284 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0285 dsigMC = 0.0103;
0286 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0287 dsigMC = 0.0090;
0288 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0289 dsigMC = 0.0190;
0290 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0291 dsigMC = 0.0156;
0292 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0293 dsigMC = 0.0269;
0294 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0295 dsigMC = 0.0287;
0296 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0297 dsigMC = 0.0364;
0298 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0299 dsigMC = 0.0321;
0300 }
0301 if (lumiRatio_ == 0.0) {
0302 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0303 dsigMC = 0.0103;
0304 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0305 dsigMC = 0.0090;
0306 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0307 dsigMC = 0.0190;
0308 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0309 dsigMC = 0.0156;
0310 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0311 dsigMC = 0.0269;
0312 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0313 dsigMC = 0.0287;
0314 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0315 dsigMC = 0.0364;
0316 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0317 dsigMC = 0.0321;
0318 }
0319 if (lumiRatio_ == 1.0) {
0320 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0321 dsigMC = 0.0109;
0322 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0323 dsigMC = 0.0099;
0324 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0325 dsigMC = 0.0182;
0326 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0327 dsigMC = 0.0200;
0328 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0329 dsigMC = 0.0282;
0330 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0331 dsigMC = 0.0309;
0332 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0333 dsigMC = 0.0386;
0334 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0335 dsigMC = 0.0359;
0336 }
0337 }
0338 }
0339 break;
0340
0341 case 2:
0342
0343 if (dataset_ == "Summer12_LegacyPaper" || dataset_ == "22Jan2013ReReco") {
0344 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0345 dsigMC = 0.0094;
0346 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0347 dsigMC = 0.0092;
0348 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0349 dsigMC = 0.0182;
0350 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0351 dsigMC = 0.0139;
0352 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0353 dsigMC = 0.0220;
0354 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0355 dsigMC = 0.0229;
0356 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0357 dsigMC = 0.0290;
0358 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0359 dsigMC = 0.0234;
0360 }
0361 break;
0362
0363 case 3:
0364
0365 if (dataset_ == "Summer11" ||
0366 dataset_ == "ReReco") {
0367 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0368 dsigMC = 0.01;
0369 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0370 dsigMC = 0.0099;
0371 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0372 dsigMC = 0.0217;
0373 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0374 dsigMC = 0.0157;
0375 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0376 dsigMC = 0.0326;
0377 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0378 dsigMC = 0.0330;
0379 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0380 dsigMC = 0.0331;
0381 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0382 dsigMC = 0.0378;
0383 } else if (dataset_ == "Fall11" ||
0384 dataset_ ==
0385 "Jan16ReReco") {
0386 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0387 dsigMC = 0.0096;
0388 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0389 dsigMC = 0.0074;
0390 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0391 dsigMC = 0.0196;
0392 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0393 dsigMC = 0.0141;
0394 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0395 dsigMC = 0.0279;
0396 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0397 dsigMC = 0.0268;
0398 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0399 dsigMC = 0.0301;
0400 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0401 dsigMC = 0.0293;
0402 } else if (dataset_ == "Summer12" ||
0403 dataset_ ==
0404 "ICHEP2012") {
0405 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0406 dsigMC = 0.0119;
0407 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0408 dsigMC = 0.0107;
0409 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0410 dsigMC = 0.0240;
0411 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0412 dsigMC = 0.0149;
0413 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0414 dsigMC = 0.0330;
0415 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0416 dsigMC = 0.0375;
0417 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0418 dsigMC = 0.0602;
0419 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0420 dsigMC = 0.0607;
0421 } else if (dataset_ == "Summer12_DR53X_HCP2012" || dataset_ == "Moriond2013") {
0422 if (isEB && fabs(eta) < 1 && r9 < 0.94)
0423 dsigMC = 0.0099;
0424 if (isEB && fabs(eta) < 1 && r9 >= 0.94)
0425 dsigMC = 0.0103;
0426 if (isEB && fabs(eta) >= 1 && r9 < 0.94)
0427 dsigMC = 0.0219;
0428 if (isEB && fabs(eta) >= 1 && r9 >= 0.94)
0429 dsigMC = 0.0158;
0430 if (!isEB && fabs(eta) < 2 && r9 < 0.94)
0431 dsigMC = 0.0222;
0432 if (!isEB && fabs(eta) < 2 && r9 >= 0.94)
0433 dsigMC = 0.0298;
0434 if (!isEB && fabs(eta) >= 2 && r9 < 0.94)
0435 dsigMC = 0.0318;
0436 if (!isEB && fabs(eta) >= 2 && r9 >= 0.94)
0437 dsigMC = 0.0302;
0438 }
0439 break;
0440 }
0441
0442 if (isMC_) {
0443 CLHEP::RandGaussQ gaussDistribution(rng->getEngine(streamID), 1., dsigMC);
0444 corrMC = gaussDistribution.fire();
0445 if (verbose_) {
0446 std::cout << "[ElectronEnergyCalibrator] unsmeared energy " << newEnergy_ << std::endl;
0447 }
0448 if (synchronization_) {
0449 std::cout << "[ElectronEnergyCalibrator] "
0450 << "======================= SYNCRONIZATION MODE! =======================" << std::endl;
0451 newEnergy_ = newEnergy_ * (1 + dsigMC);
0452 } else {
0453 newEnergy_ = newEnergy_ * corrMC;
0454 }
0455 if (verbose_) {
0456 std::cout << "[ElectronEnergyCalibrator] smeared energy " << newEnergy_ << std::endl;
0457 }
0458 }
0459
0460
0461 if (updateEnergyErrors_) {
0462 newEnergyError_ = sqrt(newEnergyError_ * newEnergyError_ + dsigMC * dsigMC * newEnergy_ * newEnergy_);
0463 }
0464 if (verbose_) {
0465 std::cout << "[ElectronEnergyCalibrator] initial energy " << electron.getRegEnergy() << " recalibrated energy "
0466 << newEnergy_ << std::endl;
0467 }
0468 if (verbose_) {
0469 std::cout << "[ElectronEnergyCalibrator] initial energy error " << electron.getRegEnergyError()
0470 << " recalibrated energy error " << newEnergyError_ << std::endl;
0471 }
0472
0473 electron.setNewEnergy(newEnergy_);
0474 electron.setNewEnergyError(newEnergyError_);
0475 }
0476
0477 void ElectronEnergyCalibrator::correctLinearity(SimpleElectron &electron) {
0478 if (!isMC_ && applyLinearityCorrection_) {
0479 bool isEB = electron.isEB();
0480 double eta = electron.getEta();
0481 double theta = 2 * atan(exp(-eta));
0482 double p = electron.getCombinedMomentum();
0483 double pt = p * fabs(sin(theta));
0484 int classification = electron.getElClass();
0485 double linscale = 0.;
0486
0487 for (int i = 0; i < nLinCorrValRaw; i++) {
0488 if ((pt >= linCorrValArray[i].ptMin) && (pt <= linCorrValArray[i].ptMax)) {
0489 if (isEB) {
0490 if (fabs(eta) < 1) {
0491 if (classification < 2) {
0492 linscale = linCorrValArray[i].corrCat0;
0493 } else {
0494 linscale = linCorrValArray[i].corrCat3;
0495 }
0496 } else {
0497 if (classification < 2) {
0498 linscale = linCorrValArray[i].corrCat1;
0499 } else {
0500 linscale = linCorrValArray[i].corrCat4;
0501 }
0502 }
0503 } else
0504 {
0505 if (classification < 2) {
0506 linscale = linCorrValArray[i].corrCat2;
0507 } else {
0508 linscale = linCorrValArray[i].corrCat5;
0509 }
0510 }
0511 }
0512 }
0513 double newP = p / (1. + linscale);
0514 if (verbose_) {
0515 std::cout << "[ElectronEnergyCalibrator] Applying a linearity correction of " << 1. / (1. + linscale)
0516 << " to " << pt << " GeV in pt" << std::endl;
0517 }
0518 electron.setCombinedMomentum(newP);
0519 if (verbose_) {
0520 std::cout << "[ElectronEnergyCalibrator] calibrated transverse momentum " << pt
0521 << " GeV recalibrated for linearity to momentum " << electron.getCombinedMomentum() * fabs(sin(theta))
0522 << " GeV" << std::endl;
0523 }
0524 }
0525 }