CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/EgammaAnalysis/ElectronTools/src/ElectronEnergyCalibrator.cc

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  * Propagate SC calibration from Zee fit to the electrons
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_ ) // DATA
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         // linearity corrections data
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 // MC
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                     // Implementation of the MC smearing for regression energy type 1
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             // Implementation of the MC smearing for regression energy type 2
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             // standard SC energy scale corrections implementation
00368             if ( dataset_ == "Summer11" || dataset_ == "ReReco" ) 
00369             { // values from https://indico.cern.ch/conferenceDisplay.py?confId=146386
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             { // values from https://hypernews.cern.ch/HyperNews/CMS/get/higgs2g/634.html, consistant with Jan16ReReco corrections
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             { // new values from https://twiki.cern.ch/twiki/pub/CMS/EcalEnergyResolutionWithZee/oriented-ICHEP-scales_resolution.pdf
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     // correct energy error for MC and for data as error is obtained from (ideal) MC parametrisation
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 // !isEB
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 }