Go to the documentation of this file.00001 #include "EgammaAnalysis/ElectronTools/interface/ElectronEnergyCalibrator.h"
00002
00003 #include <CLHEP/Random/RandGaussQ.h>
00004 #include <CLHEP/Random/RandFlat.h>
00005 #include <CLHEP/Random/Random.h>
00006 #include "FWCore/ServiceRegistry/interface/Service.h"
00007 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00008 #include "FWCore/Utilities/interface/Exception.h"
00009
00010
00011
00012
00013
00014
00015
00016 #include <string>
00017 #include <vector>
00018 #include <fstream>
00019 #include <sstream>
00020 #include <iostream>
00021
00022 using std::string;
00023 using std::vector;
00024 using std::ifstream;
00025 using std::istringstream;
00026 using std::cout;
00027 using namespace edm;
00028
00029 void ElectronEnergyCalibrator::init()
00030 {
00031 if ( !isMC_ )
00032 {
00033 if ( verbose_ )
00034 {
00035 std::cout << "[ElectronEnergyCalibrator] Initialization in DATA mode" << std::endl;
00036 }
00037
00038 ifstream fin(pathData_.c_str());
00039
00040 if (!fin){
00041 throw cms::Exception("Configuration")
00042 << "[ElectronEnergyCalibrator] Cannot open the file "
00043 << pathData_ << "\n It is not found, missed or corrupted" ;
00044 } else
00045 {
00046 if ( verbose_ )
00047 {
00048 std::cout << "[ElectronEnergyCalibrator] File "
00049 << pathData_ << " succesfully opened" << std::endl;
00050 }
00051
00052 string s;
00053 vector<string> selements;
00054 string delimiter = ",";
00055 nCorrValRaw = 0;
00056
00057 while ( !fin.eof() )
00058 {
00059 getline(fin, s);
00060 if ( !s.empty() )
00061 {
00062 splitString(s, selements, delimiter);
00063 corrValArray[nCorrValRaw].nRunMin = stringToDouble(selements[0]);
00064 corrValArray[nCorrValRaw].nRunMax = stringToDouble(selements[1]);
00065 corrValArray[nCorrValRaw].corrCat0 = stringToDouble(selements[2]);
00066 corrValArray[nCorrValRaw].corrCat1 = stringToDouble(selements[3]);
00067 corrValArray[nCorrValRaw].corrCat2 = stringToDouble(selements[4]);
00068 corrValArray[nCorrValRaw].corrCat3 = stringToDouble(selements[5]);
00069 corrValArray[nCorrValRaw].corrCat4 = stringToDouble(selements[6]);
00070 corrValArray[nCorrValRaw].corrCat5 = stringToDouble(selements[7]);
00071 corrValArray[nCorrValRaw].corrCat6 = stringToDouble(selements[8]);
00072 corrValArray[nCorrValRaw].corrCat7 = stringToDouble(selements[9]);
00073
00074 nCorrValRaw++;
00075
00076 selements.clear();
00077 }
00078 }
00079
00080 fin.close();
00081
00082 if ( verbose_ )
00083 {
00084 std::cout << "[ElectronEnergyCalibrator] File closed" << std::endl;
00085 }
00086
00087 }
00088
00089 if(applyLinearityCorrection_)
00090 {
00091 ifstream finlin(pathLinData_.c_str());
00092
00093 if (!finlin)
00094 {
00095 throw cms::Exception("Configuration")
00096 << "[ElectronEnergyCalibrator] Cannot open the file "<< pathLinData_ << "\n It is not found, missed or corrupted" ;
00097 } else
00098 {
00099 if (verbose_)
00100 {
00101 std::cout<<"[ElectronEnergyCalibrator] File with Linearity Corrections "<<pathLinData_<<" succesfully opened"<<std::endl;
00102 }
00103
00104 string s;
00105 vector<string> selements;
00106 string delimiter = ",";
00107 nLinCorrValRaw = 0;
00108
00109 while ( !finlin.eof() )
00110 {
00111 getline(finlin, s);
00112 if ( !s.empty() )
00113 {
00114 splitString(s, selements, delimiter);
00115
00116 linCorrValArray[nLinCorrValRaw].ptMin = stringToDouble(selements[0]);
00117 linCorrValArray[nLinCorrValRaw].ptMax = stringToDouble(selements[1]);
00118 linCorrValArray[nLinCorrValRaw].corrCat0 = stringToDouble(selements[2]);
00119 linCorrValArray[nLinCorrValRaw].corrCat1 = stringToDouble(selements[3]);
00120 linCorrValArray[nLinCorrValRaw].corrCat2 = stringToDouble(selements[4]);
00121 linCorrValArray[nLinCorrValRaw].corrCat3 = stringToDouble(selements[5]);
00122 linCorrValArray[nLinCorrValRaw].corrCat4 = stringToDouble(selements[6]);
00123 linCorrValArray[nLinCorrValRaw].corrCat5 = stringToDouble(selements[7]);
00124
00125 nLinCorrValRaw++;
00126
00127 selements.clear();
00128 }
00129 }
00130
00131 finlin.close();
00132
00133 if (verbose_)
00134 {
00135 std::cout<<"[ElectronEnergyCalibrator] File closed"<<std::endl;
00136 }
00137 }
00138 }
00139 } else
00140 {
00141 if ( verbose_ )
00142 {
00143 std::cout << "[ElectronEnergyCalibrator] Initialization in MC mode" << std::endl;
00144 }
00145 }
00146 }
00147
00148 void ElectronEnergyCalibrator::splitString(const string &fullstr, vector<string> &elements, const string &delimiter)
00149 {
00150 string::size_type lastpos = fullstr.find_first_not_of(delimiter, 0);
00151 string::size_type pos = fullstr.find_first_of(delimiter, lastpos);
00152
00153 while ( ( string::npos != pos ) || ( string::npos != lastpos ) )
00154 {
00155 elements.push_back(fullstr.substr(lastpos, pos-lastpos));
00156 lastpos = fullstr.find_first_not_of(delimiter, pos);
00157 pos = fullstr.find_first_of(delimiter, lastpos);
00158 }
00159 }
00160
00161 double ElectronEnergyCalibrator::stringToDouble(const string &str)
00162 {
00163 istringstream stm;
00164 double val = 0;
00165 stm.str(str);
00166 stm >> val;
00167 return val;
00168 }
00169
00170 void ElectronEnergyCalibrator::calibrate(SimpleElectron &electron)
00171 {
00172 double scale = 1.0;
00173 double dsigMC=0.;
00174 double corrMC=0.;
00175 double run_ = electron.getRunNumber();
00176 bool isEB = electron.isEB();
00177 double eta = electron.getEta();
00178 double r9 = electron.getR9();
00179
00180 switch ( correctionsType_ )
00181 {
00182 case 1:
00183 if ( verbose_ )
00184 {
00185 std::cout << "[ElectronEnergyCalibrator] Using regression energy for calibration" << std::endl;
00186 }
00187 newEnergy_ = electron.getRegEnergy();
00188 newEnergyError_ = electron.getRegEnergyError();
00189 break;
00190 case 2:
00191 if( verbose_ )
00192 {
00193 std::cout << "[ElectronEnergyCalibrator] Using scale corrections for new regression" << std::endl;
00194 }
00195 newEnergy_ = electron.getRegEnergy();
00196 newEnergyError_ = electron.getRegEnergyError();
00197 break;
00198 case 3:
00199 if ( verbose_ )
00200 {
00201 std::cout << "[ElectronEnergyCalibrator] Using standard ecal energy for calibration" << std::endl;
00202 }
00203 newEnergy_ = electron.getSCEnergy();
00204 newEnergyError_ = electron.getSCEnergyError();
00205 break;
00206 }
00207
00208 edm::Service<edm::RandomNumberGenerator> rng;
00209 if ( !rng.isAvailable() )
00210 {
00211 throw cms::Exception("Configuration")
00212 << "XXXXXXX requires the RandomNumberGeneratorService\n"
00213 "which is not present in the configuration file. You must add the service\n"
00214 "in the configuration file or remove the modules that require it.";
00215 }
00216
00217 if (!isMC_ )
00218 {
00219 for ( int i=0; i < nCorrValRaw; i++ )
00220 {
00221 if ( (run_ >= corrValArray[i].nRunMin ) && ( run_ <= corrValArray[i].nRunMax ) )
00222 {
00223 if ( isEB )
00224 {
00225 if ( fabs(eta) < 1 )
00226 {
00227 if ( r9<0.94 )
00228 {
00229 scale = corrValArray[i].corrCat0;
00230 } else
00231 {
00232 scale = corrValArray[i].corrCat1;
00233 }
00234 } else
00235 {
00236 if ( r9<0.94 )
00237 {
00238 scale = corrValArray[i].corrCat2;
00239 } else
00240 {
00241 scale = corrValArray[i].corrCat3;
00242 }
00243 }
00244 } else
00245 {
00246 if ( fabs(eta) < 2 )
00247 {
00248 if ( r9<0.94 )
00249 {
00250 scale = corrValArray[i].corrCat4;
00251 } else
00252 {
00253 scale = corrValArray[i].corrCat5;
00254 }
00255 } else
00256 {
00257 if ( r9<0.94 )
00258 {
00259 scale = corrValArray[i].corrCat6;
00260 } else
00261 {
00262 scale = corrValArray[i].corrCat7;
00263 }
00264 }
00265 }
00266 }
00267 }
00268 newEnergy_ = newEnergy_*scale;
00269 }
00270
00271 switch ( correctionsType_ )
00272 {
00273 case 1:
00274
00275 if ( dataset_ == "Summer12_DR53X_HCP2012" || dataset_ == "Moriond2013" )
00276 {
00277 if ( !isMC_ )
00278 {
00279 if ( run_ <= 203002 )
00280 {
00281 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
00282 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
00283 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
00284 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
00285 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
00286 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
00287 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
00288 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
00289 } else
00290 {
00291 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
00292 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
00293 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
00294 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
00295 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
00296 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
00297 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
00298 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
00299 }
00300 } else
00301 {
00302 CLHEP::RandFlat flatRandom(rng->getEngine());
00303 double rn = flatRandom.fire();
00304 if ( rn > lumiRatio_ )
00305 {
00306 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
00307 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
00308 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
00309 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
00310 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
00311 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
00312 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
00313 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
00314 } else
00315 {
00316 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
00317 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
00318 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
00319 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
00320 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
00321 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
00322 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
00323 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
00324 }
00325 if ( lumiRatio_ == 0.0 )
00326 {
00327 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
00328 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
00329 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
00330 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
00331 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
00332 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
00333 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
00334 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
00335 }
00336 if ( lumiRatio_ == 1.0 )
00337 {
00338 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
00339 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
00340 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
00341 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
00342 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
00343 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
00344 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
00345 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
00346 }
00347 }
00348 }
00349 break;
00350
00351 case 2:
00352
00353 if ( dataset_ == "Summer12_LegacyPaper" || dataset_ == "22Jan2013ReReco" )
00354 {
00355 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0094;
00356 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0092;
00357 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
00358 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0139;
00359 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0220;
00360 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0229;
00361 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0290;
00362 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0234;
00363 }
00364 break;
00365
00366 case 3:
00367
00368 if ( dataset_ == "Summer11" || dataset_ == "ReReco" )
00369 {
00370 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.01;
00371 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
00372 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0217;
00373 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0157;
00374 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0326;
00375 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0330;
00376 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0331;
00377 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0378;
00378 } else if ( dataset_ == "Fall11" || dataset_ == "Jan16ReReco" )
00379 {
00380 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0096;
00381 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0074;
00382 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0196;
00383 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0141;
00384 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0279;
00385 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0268;
00386 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0301;
00387 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0293;
00388 } else if ( dataset_ == "Summer12" || dataset_ == "ICHEP2012" )
00389 {
00390 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0119;
00391 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0107;
00392 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0240;
00393 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0149;
00394 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0330;
00395 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0375;
00396 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0602;
00397 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0607;
00398 } else if ( dataset_ == "Summer12_DR53X_HCP2012" || dataset_ == "Moriond2013" )
00399 {
00400 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0099;
00401 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0103;
00402 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0219;
00403 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0158;
00404 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0222;
00405 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0298;
00406 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0318;
00407 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0302;
00408 }
00409 break;
00410 }
00411
00412 if ( isMC_ )
00413 {
00414 CLHEP::RandGaussQ gaussDistribution(rng->getEngine(), 1.,dsigMC);
00415 corrMC = gaussDistribution.fire();
00416 if ( verbose_ )
00417 {
00418 std::cout << "[ElectronEnergyCalibrator] unsmeared energy " << newEnergy_ << std::endl;
00419 }
00420 if ( synchronization_ )
00421 {
00422 std::cout << "[ElectronEnergyCalibrator] "
00423 << "======================= SYNCRONIZATION MODE! ======================="
00424 << std::endl;
00425 newEnergy_ = newEnergy_*(1+dsigMC);
00426 } else
00427 {
00428 newEnergy_ = newEnergy_*corrMC;
00429 }
00430 if ( verbose_ )
00431 {
00432 std::cout << "[ElectronEnergyCalibrator] smeared energy " << newEnergy_ << std::endl;
00433 }
00434 }
00435
00436
00437 if ( updateEnergyErrors_ )
00438 {
00439 newEnergyError_ = sqrt(newEnergyError_*newEnergyError_ + dsigMC*dsigMC*newEnergy_*newEnergy_);
00440 }
00441 if ( verbose_ )
00442 {
00443 std::cout << "[ElectronEnergyCalibrator] initial energy "
00444 << electron.getRegEnergy() << " recalibrated energy " << newEnergy_ << std::endl;
00445 }
00446 if ( verbose_ )
00447 {
00448 std::cout << "[ElectronEnergyCalibrator] initial energy error "
00449 << electron.getRegEnergyError() << " recalibrated energy error "
00450 << newEnergyError_ << std::endl;
00451 }
00452
00453 electron.setNewEnergy(newEnergy_);
00454 electron.setNewEnergyError(newEnergyError_);
00455 }
00456
00457 void ElectronEnergyCalibrator::correctLinearity(SimpleElectron &electron)
00458 {
00459 if(!isMC_ && applyLinearityCorrection_)
00460 {
00461 bool isEB = electron.isEB();
00462 double eta = electron.getEta();
00463 double theta = 2*atan(exp(-eta));
00464 double p = electron.getCombinedMomentum();
00465 double pt = p * fabs(sin(theta));
00466 int classification = electron.getElClass();
00467 double linscale = 0.;
00468
00469 for (int i=0; i < nLinCorrValRaw; i++)
00470 {
00471 if ((pt >= linCorrValArray[i].ptMin) && (pt <= linCorrValArray[i].ptMax))
00472 {
00473 if (isEB)
00474 {
00475 if (fabs(eta) < 1)
00476 {
00477 if (classification<2)
00478 {
00479 linscale = linCorrValArray[i].corrCat0;
00480 } else
00481 {
00482 linscale = linCorrValArray[i].corrCat3;
00483 }
00484 } else
00485 {
00486 if (classification<2)
00487 {
00488 linscale = linCorrValArray[i].corrCat1;
00489 } else
00490 {
00491 linscale = linCorrValArray[i].corrCat4;
00492 }
00493 }
00494 } else
00495 {
00496 if (classification<2)
00497 {
00498 linscale = linCorrValArray[i].corrCat2;
00499 } else
00500 {
00501 linscale = linCorrValArray[i].corrCat5;
00502 }
00503 }
00504 }
00505 }
00506 double newP = p/(1.+linscale);
00507 if (verbose_)
00508 {
00509 std::cout << "[ElectronEnergyCalibrator] Applying a linearity correction of " << 1./(1.+linscale) << " to " << pt << " GeV in pt" << std::endl;
00510 }
00511 electron.setCombinedMomentum(newP);
00512 if (verbose_)
00513 {
00514 std::cout << "[ElectronEnergyCalibrator] calibrated transverse momentum " << pt << " GeV recalibrated for linearity to momentum " << electron.getCombinedMomentum()*fabs(sin(theta)) << " GeV" << std::endl;
00515 }
00516 }
00517 }