3 #include <CLHEP/Random/RandGaussQ.h>
4 #include <CLHEP/Random/RandFlat.h>
5 #include <CLHEP/Random/Random.h>
25 using std::istringstream;
35 std::cout <<
"[ElectronEnergyCalibrator] Initialization in DATA mode" << std::endl;
38 ifstream
fin(pathData_.c_str());
42 <<
"[ElectronEnergyCalibrator] Cannot open the file "
43 << pathData_ <<
"\n It is not found, missed or corrupted" ;
48 std::cout <<
"[ElectronEnergyCalibrator] File "
49 << pathData_ <<
" succesfully opened" << std::endl;
53 vector<string> selements;
54 string delimiter =
",";
63 corrValArray[nCorrValRaw].nRunMin = stringToDouble(selements[0]);
64 corrValArray[nCorrValRaw].nRunMax = stringToDouble(selements[1]);
65 corrValArray[nCorrValRaw].corrCat0 = stringToDouble(selements[2]);
66 corrValArray[nCorrValRaw].corrCat1 = stringToDouble(selements[3]);
67 corrValArray[nCorrValRaw].corrCat2 = stringToDouble(selements[4]);
68 corrValArray[nCorrValRaw].corrCat3 = stringToDouble(selements[5]);
69 corrValArray[nCorrValRaw].corrCat4 = stringToDouble(selements[6]);
70 corrValArray[nCorrValRaw].corrCat5 = stringToDouble(selements[7]);
71 corrValArray[nCorrValRaw].corrCat6 = stringToDouble(selements[8]);
72 corrValArray[nCorrValRaw].corrCat7 = stringToDouble(selements[9]);
84 std::cout <<
"[ElectronEnergyCalibrator] File closed" << std::endl;
89 if(applyLinearityCorrection_)
91 ifstream finlin(pathLinData_.c_str());
96 <<
"[ElectronEnergyCalibrator] Cannot open the file "<< pathLinData_ <<
"\n It is not found, missed or corrupted" ;
101 std::cout<<
"[ElectronEnergyCalibrator] File with Linearity Corrections "<<pathLinData_<<
" succesfully opened"<<std::endl;
105 vector<string> selements;
106 string delimiter =
",";
109 while ( !finlin.eof() )
116 linCorrValArray[nLinCorrValRaw].ptMin = stringToDouble(selements[0]);
117 linCorrValArray[nLinCorrValRaw].ptMax = stringToDouble(selements[1]);
118 linCorrValArray[nLinCorrValRaw].corrCat0 = stringToDouble(selements[2]);
119 linCorrValArray[nLinCorrValRaw].corrCat1 = stringToDouble(selements[3]);
120 linCorrValArray[nLinCorrValRaw].corrCat2 = stringToDouble(selements[4]);
121 linCorrValArray[nLinCorrValRaw].corrCat3 = stringToDouble(selements[5]);
122 linCorrValArray[nLinCorrValRaw].corrCat4 = stringToDouble(selements[6]);
123 linCorrValArray[nLinCorrValRaw].corrCat5 = stringToDouble(selements[7]);
135 std::cout<<
"[ElectronEnergyCalibrator] File closed"<<std::endl;
143 std::cout <<
"[ElectronEnergyCalibrator] Initialization in MC mode" << std::endl;
153 while ( ( string::npos != pos ) || ( string::npos != lastpos ) )
155 elements.push_back(fullstr.substr(lastpos, pos-lastpos));
156 lastpos = fullstr.find_first_not_of(delimiter, pos);
157 pos = fullstr.find_first_of(delimiter, lastpos);
176 bool isEB = electron.
isEB();
178 double r9 = electron.
getR9();
180 switch ( correctionsType_ )
185 std::cout <<
"[ElectronEnergyCalibrator] Using regression energy for calibration" << std::endl;
193 std::cout <<
"[ElectronEnergyCalibrator] Using scale corrections for new regression" << std::endl;
201 std::cout <<
"[ElectronEnergyCalibrator] Using standard ecal energy for calibration" << std::endl;
212 <<
"XXXXXXX requires the RandomNumberGeneratorService\n"
213 "which is not present in the configuration file. You must add the service\n"
214 "in the configuration file or remove the modules that require it.";
219 for (
int i=0;
i < nCorrValRaw;
i++ )
221 if ( (run_ >= corrValArray[
i].nRunMin ) && ( run_ <= corrValArray[
i].nRunMax ) )
229 scale = corrValArray[
i].corrCat0;
232 scale = corrValArray[
i].corrCat1;
238 scale = corrValArray[
i].corrCat2;
241 scale = corrValArray[
i].corrCat3;
250 scale = corrValArray[
i].corrCat4;
253 scale = corrValArray[
i].corrCat5;
259 scale = corrValArray[
i].corrCat6;
262 scale = corrValArray[
i].corrCat7;
268 newEnergy_ = newEnergy_*
scale;
271 switch ( correctionsType_ )
275 if ( dataset_ ==
"Summer12_DR53X_HCP2012" || dataset_ ==
"Moriond2013" )
279 if ( run_ <= 203002 )
281 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
282 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
283 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
284 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
285 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
286 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
287 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
288 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
291 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
292 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
293 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
294 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
295 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
296 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
297 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
298 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
302 CLHEP::RandFlat flatRandom(rng->
getEngine());
303 double rn = flatRandom.fire();
304 if ( rn > lumiRatio_ )
306 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
307 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
308 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
309 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
310 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
311 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
312 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
313 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
316 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
317 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
318 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
319 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
320 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
321 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
322 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
323 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
325 if ( lumiRatio_ == 0.0 )
327 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0103;
328 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0090;
329 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0190;
330 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0156;
331 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0269;
332 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0287;
333 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0364;
334 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0321;
336 if ( lumiRatio_ == 1.0 )
338 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0109;
339 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
340 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
341 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0200;
342 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0282;
343 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0309;
344 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0386;
345 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0359;
353 if ( dataset_ ==
"Summer12_LegacyPaper" || dataset_ ==
"22Jan2013ReReco" )
355 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0094;
356 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0092;
357 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0182;
358 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0139;
359 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0220;
360 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0229;
361 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0290;
362 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0234;
368 if ( dataset_ ==
"Summer11" || dataset_ ==
"ReReco" )
370 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.01;
371 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0099;
372 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0217;
373 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0157;
374 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0326;
375 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0330;
376 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0331;
377 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0378;
378 }
else if ( dataset_ ==
"Fall11" || dataset_ ==
"Jan16ReReco" )
380 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0096;
381 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0074;
382 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0196;
383 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0141;
384 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0279;
385 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0268;
386 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0301;
387 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0293;
388 }
else if ( dataset_ ==
"Summer12" || dataset_ ==
"ICHEP2012" )
390 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0119;
391 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0107;
392 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0240;
393 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0149;
394 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0330;
395 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0375;
396 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0602;
397 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0607;
398 }
else if ( dataset_ ==
"Summer12_DR53X_HCP2012" || dataset_ ==
"Moriond2013" )
400 if ( isEB && fabs(eta) < 1 && r9 < 0.94 ) dsigMC = 0.0099;
401 if ( isEB && fabs(eta) < 1 && r9 >= 0.94 ) dsigMC = 0.0103;
402 if ( isEB && fabs(eta) >= 1 && r9 < 0.94 ) dsigMC = 0.0219;
403 if ( isEB && fabs(eta) >= 1 && r9 >= 0.94 ) dsigMC = 0.0158;
404 if ( !isEB && fabs(eta) < 2 && r9 < 0.94 ) dsigMC = 0.0222;
405 if ( !isEB && fabs(eta) < 2 && r9 >= 0.94 ) dsigMC = 0.0298;
406 if ( !isEB && fabs(eta) >= 2 && r9 < 0.94 ) dsigMC = 0.0318;
407 if ( !isEB && fabs(eta) >= 2 && r9 >= 0.94 ) dsigMC = 0.0302;
414 CLHEP::RandGaussQ gaussDistribution(rng->
getEngine(), 1.,dsigMC);
415 corrMC = gaussDistribution.fire();
418 std::cout <<
"[ElectronEnergyCalibrator] unsmeared energy " << newEnergy_ << std::endl;
420 if ( synchronization_ )
422 std::cout <<
"[ElectronEnergyCalibrator] "
423 <<
"======================= SYNCRONIZATION MODE! ======================="
425 newEnergy_ = newEnergy_*(1+dsigMC);
428 newEnergy_ = newEnergy_*corrMC;
432 std::cout <<
"[ElectronEnergyCalibrator] smeared energy " << newEnergy_ << std::endl;
437 if ( updateEnergyErrors_ )
439 newEnergyError_ =
sqrt(newEnergyError_*newEnergyError_ + dsigMC*dsigMC*newEnergy_*newEnergy_);
443 std::cout <<
"[ElectronEnergyCalibrator] initial energy "
444 << electron.
getRegEnergy() <<
" recalibrated energy " << newEnergy_ << std::endl;
448 std::cout <<
"[ElectronEnergyCalibrator] initial energy error "
450 << newEnergyError_ << std::endl;
459 if(!isMC_ && applyLinearityCorrection_)
461 bool isEB = electron.
isEB();
465 double pt = p * fabs(
sin(theta));
467 double linscale = 0.;
469 for (
int i=0;
i < nLinCorrValRaw;
i++)
471 if ((pt >= linCorrValArray[
i].
ptMin) && (pt <= linCorrValArray[
i].
ptMax))
477 if (classification<2)
479 linscale = linCorrValArray[
i].corrCat0;
482 linscale = linCorrValArray[
i].corrCat3;
486 if (classification<2)
488 linscale = linCorrValArray[
i].corrCat1;
491 linscale = linCorrValArray[
i].corrCat4;
496 if (classification<2)
498 linscale = linCorrValArray[
i].corrCat2;
501 linscale = linCorrValArray[
i].corrCat5;
506 double newP = p/(1.+linscale);
509 std::cout <<
"[ElectronEnergyCalibrator] Applying a linearity correction of " << 1./(1.+linscale) <<
" to " << pt <<
" GeV in pt" << std::endl;
514 std::cout <<
"[ElectronEnergyCalibrator] calibrated transverse momentum " << pt <<
" GeV recalibrated for linearity to momentum " << electron.
getCombinedMomentum()*fabs(
sin(theta)) <<
" GeV" << std::endl;
void calibrate(SimpleElectron &electron)
double getRegEnergyError()
void setNewEnergyError(double newEnergyError)
void setCombinedMomentum(double combinedMomentum)
void setNewEnergy(double newEnergy)
Sin< T >::type sin(const T &t)
Geom::Theta< T > theta() const
std::vector< std::string > splitString(const std::string &fLine)
double getCombinedMomentum()
void splitString(const string &fullstr, vector< string > &elements, const string &delimiter)
double stringToDouble(const string &str)
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
void correctLinearity(SimpleElectron &electron)
double getSCEnergyError()