ForIA: /Users/jmonk/Physics/ForIA/ForIA/AnalysisTools/OffsetEtaJES/OffsetEtaJES.hh Source File

Go to the documentation of this file.
00001 /*
00002  *  OffsetEtaJES (formerly EtaJES)
00003  *
00004  *  - Derived from dijet MC by matching truth and reco jets
00005  *    Reconstructed energy and eta are corrected such that
00006  *    the corrected reco quantities match the true ones,
00007  *    when binned in truth E, and detctor eta of reco jet
00008  *    (="impact point" of reco jet)
00009  *    
00010  *  - The correction is only derived for isolated jets
00011  *    meaning no jet with pT > 7 GeV in cone of 2.5*R (= 1.5 here)
00012  *
00013  *  - For these constants, that were derived from the very first
00014  *    rel16 MC10 dijet MC (200k evts per JX). OTX correction is applied
00015  *    to rel 16 reco, but discrepancies of calorimeter response in the dead 
00016  *    regions are still observed. 
00017  *    Jets that fall in bad OTX regions were _NOT_ ignored.
00018  *
00019  *  - No weighting was appied to the sample
00020  *  
00021  *  - In the code, the quantities with name of the form
00022  *       *_uncorr
00023  *    is the quantities before correction, i.e. EM scale
00024  *    for EM+JES calibration, GCW scale for GCW+JES calibration
00025  *       eta_det
00026  *    meand detector eta, i.e. the (raw) EM scale eta for _all_
00027  *    calibrations.
00028  *       eta/phi
00029  *    without syntax can be eta/phi after origin correction, or
00030  *    or the "uncorr" quantities.
00031  *
00032  *  Dag Gillberg, Oct 2010
00033  *  David W. Miller, Jan 2011
00034  *
00035  *  Modification, Jan 20 2011, D. Miller
00036  *
00037  *  - Code extended to also be able to apply pile-up correction
00038  *
00039  */
00040 
00041 #ifndef OFFSETETAJES_H_
00042 #define OFFSETETAJES_H_
00043 
00044 #include <iostream>
00045 #include "TLorentzVector.h"
00046 #include "TString.h"
00047 
00048 // The calibration factors
00049 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt4TopoEM_EtaMassJES.hh"
00050 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt6TopoEM_EtaMassJES.hh"
00051 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt4TowerEM_EtaMassJES.hh"
00052 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt6TowerEM_EtaMassJES.hh"
00053 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt4TopoGCW_EtaMassJES.hh"
00054 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt6TopoGCW_EtaMassJES.hh"
00055 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt4TowerGCW_EtaMassJES.hh"
00056 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt6TowerGCW_EtaMassJES.hh"
00057 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt4LCTopo_EtaMassJES.hh"
00058 #include "ForIA/AnalysisTools/OffsetEtaJES/AntiKt6LCTopo_EtaMassJES.hh"
00059 
00060 #include "ForIA/AnalysisTools/OffsetEtaJES/OffsetConstants.hh"
00061 
00062 class OffsetEtaJES {
00063   
00064  public:
00065   OffsetEtaJES(TString jetAlgo="AntiKt6TopoEM")
00066     : useTowers(false) {
00067     init(jetAlgo);
00068   };
00069   
00070  private:
00071   void error(TString msg) { 
00072     printf("\nERROR - OffsetEtaJES:\n\n  %s\n\n",msg.Data()); 
00073     abort();
00074   }
00075 
00076   // 90 eta bins, and up to 7 parameter for the pol-fit
00077   const static int NetaBins=90;
00078   const static int Npar=7;
00079   
00080   double JES_factors[NetaBins][Npar];
00081   double etaCorr_factors[NetaBins][Npar];
00082   double JMS_factors[NetaBins][Npar];
00083   
00084   const static int NOffsetEtaBins=100;
00085   const static int NPVmax=5;
00086   double offset_factors[NOffsetEtaBins][NPVmax];
00087   bool useTowers;
00088 
00089   // Some constraints that apply to the Oct MC10 factors
00090   double _minPt_JES, _minPt_EtaCorr, _maxE_EtaCorr;
00091 
00092  public:
00093   
00094   double GetJES(double E_uncorr, double eta_det)
00095   {
00096     double GeV=1000;
00097     double E = E_uncorr/GeV;
00098     if ( E/cosh(eta_det) <_minPt_JES ) E= _minPt_JES*cosh(eta_det);
00099     
00100     // Get the factors
00101     int ieta = GetEtaBin(eta_det);
00102     const double *factors = JES_factors[ieta];
00103     
00104     // Calculate the jet response and then the JES as 1/R
00105     double R = 
00106       GetLogPolN(factors,E);
00107     //factors[0]*(1.0-factors[1]*TMath::Power(E/0.75,factors[2]-1.0));
00108     return 1.0/R;
00109   }
00110   
00111   int GetOffsetEtaBin(double eta_det)
00112   {
00113     int ieta=int(eta_det*10)+NOffsetEtaBins/2;
00114     if (eta_det<0) ieta-=1;
00115     if ( ieta <  0 ) return 0;
00116     if ( ieta > NOffsetEtaBins-1 ) return NOffsetEtaBins-1;
00117     return ieta;
00118   }
00119 
00120   double GetRel16EMscaleFactor(double eta_det, double emf) {
00121     if (fabs(eta_det)<2.5) return 1.0;
00122 
00123     // scale em-fraction of jet acc. to Z->ee meas.
00124     if (fabs(eta_det)<2.8) return (1.0-emf+emf/(1.0-0.034));
00125     if (fabs(eta_det)<3.2) return (1.0-emf+emf/(1.0-0.019));
00126 
00127     // For the FCal, hadronic layers calibrated wrt EM layer:
00128     return 1.0/(1+0.05);
00129   }
00130   
00131   double GetOffsetCorrection(double eta_det, 
00132                              int Npv,
00133                              double nTowers = 1.0) { 
00134   
00135     double jet_et_offset = 0.0;
00136   
00137     // Take care of the Npv = 1 case right away
00138     if (Npv < 1) {
00139       return jet_et_offset;
00140     }
00141   
00142     // Eta bin
00143     int ieta = GetOffsetEtaBin(eta_det);
00144   
00145     // Set the maximum number of vertices                              
00146     if (Npv > NPVmax) Npv = NPVmax;
00147     
00148     // Get the correction
00149     if (useTowers) jet_et_offset = nTowers * offset_factors[ieta][Npv-1];
00150     else jet_et_offset = offset_factors[ieta][Npv-1];
00151   
00152     // Done
00153     return jet_et_offset;
00154   }
00155   
00156   double GetLogPolN(const double *factors, double x) {
00157     double y=0;
00158     for (int i=0;i<Npar;++i)
00159       y += factors[i]*pow(log(x),i);
00160     return y;
00161   }
00162   
00163   int GetEtaBin(double eta_det)
00164   {
00165     int ieta=int(eta_det*10)+45;
00166     if (eta_det<0) ieta-=1;
00167     if ( ieta <  0 ) return 0;
00168     if ( ieta > 89 ) return 89;
00169     return ieta;
00170   }
00171 
00172   TLorentzVector ApplyJES(double E_uncorr, double eta_det, double eta, 
00173                           double phi, double mass_uncorr){
00174     TLorentzVector jet;
00175     double JES=GetJES(E_uncorr,eta_det);
00176     
00177     double p=E_uncorr>mass_uncorr?sqrt(E_uncorr*E_uncorr-mass_uncorr*mass_uncorr):0;
00178     jet.SetPtEtaPhiM(p/cosh(eta),eta,phi,mass_uncorr);
00179     
00180     jet*=JES; // Scale the jet by the JES
00181     return jet;
00182   }
00183   
00184   TLorentzVector ApplyEtaJES(double E_uncorr, double eta_det, 
00185                              double eta, double phi, double mass_uncorr) {
00186     // Get corrections
00187     double JES=GetJES(E_uncorr,eta_det);
00188     double etaCorr=GetEtaCorr(E_uncorr*JES,eta_det);
00189     
00190     // Apply the corection
00191     TLorentzVector jet;
00192     double E=E_uncorr*JES, m=mass_uncorr*JES, corrEta=eta+etaCorr;
00193     double p= ( E>m ? sqrt(E*E-m*m) : 0 );
00194 
00195     jet.SetPtEtaPhiM(p/cosh(corrEta),corrEta,phi,m);
00196     return jet;
00197   }
00198 
00199   TLorentzVector ApplyEtaMassJES(double E_uncorr, double eta_det, double eta, double phi, double mass_uncorr){
00200     // Get all corrections
00201     double JES=GetJES(E_uncorr,eta_det); // Jet Energy Scale
00202     double etaCorr=GetEtaCorr(E_uncorr*JES,eta_det); // Jet Eta correction
00203     double JMS=GetMassCorr(E_uncorr*JES,eta_det); // Jet Mass Scale
00204     
00205     // Apply them
00206     TLorentzVector jet;
00207     double E=E_uncorr*JES, m=mass_uncorr*JMS, corrEta=eta+etaCorr;
00208     double p= ( E>m ? sqrt(E*E-m*m) : 0 );
00209     //  printf("E: %.2f/m: %.2f %.2f\n",E,m,JMS);
00210     jet.SetPtEtaPhiM(p/cosh(corrEta),corrEta,phi,m);
00211     return jet;
00212   }
00213   
00214   double GetEtaCorr(double Ecorr, double eta_det) {
00215     int ieta = GetEtaBin(eta_det);
00216     const double *factors = etaCorr_factors[ieta];
00217     
00218     double GeV=1000;
00219     double E = Ecorr/GeV;
00220     if ( E < _minPt_EtaCorr*cosh(eta_det) ) 
00221       E = _minPt_EtaCorr*cosh(eta_det);
00222     if ( E>_maxE_EtaCorr ) E=_maxE_EtaCorr;
00223     
00224     double etaCorr = GetLogPolN(factors,E);
00225     
00226     // This is ( reco_eta - truth_eta )
00227     // to make it an additative correction return the negative value
00228     return -etaCorr;
00229   }
00230   
00231   double GetMassCorr(double Ecorr, double eta_det) {
00232     error("Mass correction is not supported");
00233     
00234     int ieta = GetEtaBin(eta_det);
00235     const double *factors = JMS_factors[ieta];
00236 
00237     double GeV=1000;
00238     double E = ( Ecorr/cosh(eta_det)<5.0*GeV ? 5.0*cosh(eta_det) : Ecorr/GeV );
00239     //double logE = log(E);
00240     
00241     double massR = GetLogPolN(factors,E);
00242     return 1.0/massR;
00243   }
00244 
00245   TLorentzVector ApplyOffsetEtaJES(double E_uncorr, double eta_det, 
00246                                    double eta, double phi, double mass_uncorr,
00247                                    int NPV, double nTowers = 1.0) {
00248     double jetET_offset = GetOffsetCorrection(eta_det, NPV, nTowers);
00249     double jetET = E_uncorr/cosh(eta_det);
00250     
00251     double offset_r = (jetET-jetET_offset)/jetET;
00252     if (offset_r<0) offset_r=0; // protect against neg energy!
00253 
00254     return ApplyEtaJES(E_uncorr*offset_r,eta_det,eta,phi,mass_uncorr*offset_r);
00255   }
00256 
00257   TLorentzVector ApplyOffsetEtaJES_Data(double E_uncorr, double eta_det, 
00258                                         double eta, double phi, double mass_uncorr,
00259                                         int NPV, double emf, double nTowers = 1.0) {
00260     double EMsf = GetRel16EMscaleFactor(eta_det,emf);
00261 
00262     double jetET_offset = GetOffsetCorrection(eta_det, NPV, nTowers);
00263     double jetET = E_uncorr/cosh(eta_det);
00264     
00265     double offset_r = (jetET-jetET_offset)/jetET;
00266     if (offset_r<0) offset_r=0; // protect against neg energy!
00267 
00268     return ApplyEtaJES(EMsf*E_uncorr*offset_r,eta_det,eta,phi,mass_uncorr*offset_r*EMsf);
00269   }
00270 
00271   
00272   void init(TString jetAlgo) {
00273     
00274     printf("===================================\n\n");
00275     printf("  Initializing the OffsetEtaJES jet calibration tool\n\n");
00276     //printf("  WARNING - This version of the code has not yet been validated (Jan 2011).\n");
00277     //printf("            Please upgrade to a later version (JetCalibTool tag) when available.\n\n");
00278     printf("  Calibration of %s jets\n",jetAlgo.Data());
00279     printf("  For data produced with rel 16.0 (the full 2010 dataset)\n");
00280     printf("\n===================================\n");
00281 
00282     // Used when deriving first rel 16 JES 
00283     // Oct, 2010
00284     _minPt_JES     = 10; // GeV
00285     _minPt_EtaCorr = 8.0; // GeV
00286     _maxE_EtaCorr = 2500; // GeV
00287     
00288     for (int ieta=0; ieta<NetaBins; ++ieta) {
00289       for (int par=0; par<Npar; ++par) {
00290         if (jetAlgo=="AntiKt6TopoEM") {
00291           JES_factors[ieta][par]     = AntiKt6TopoEM_JES[ieta][par];
00292           etaCorr_factors[ieta][par] = AntiKt6TopoEM_EtaCorr[ieta][par];
00293           //JMS_factors[ieta][par]     = AntiKt6EMTopo_MassCorr[ieta][par];
00294         } else if (jetAlgo=="AntiKt4TopoEM") {
00295           JES_factors[ieta][par]     = AntiKt4TopoEM_JES[ieta][par];
00296           etaCorr_factors[ieta][par] = AntiKt4TopoEM_EtaCorr[ieta][par];
00297         } else if (jetAlgo=="AntiKt6TowerEM") {
00298           JES_factors[ieta][par]     = AntiKt6TowerEM_JES[ieta][par];
00299           etaCorr_factors[ieta][par] = AntiKt6TowerEM_EtaCorr[ieta][par];
00300           //JMS_factors[ieta][par]     = AntiKt6EMTower_MassCorr[ieta][par];
00301         } else if (jetAlgo=="AntiKt4TowerEM") {
00302           JES_factors[ieta][par]     = AntiKt4TowerEM_JES[ieta][par];
00303           etaCorr_factors[ieta][par] = AntiKt4TowerEM_EtaCorr[ieta][par];
00304         } else if (jetAlgo=="AntiKt4TopoGCW") {
00305           JES_factors[ieta][par]     = AntiKt4TopoGCW_JES[ieta][par];
00306           etaCorr_factors[ieta][par] = AntiKt4TopoGCW_EtaCorr[ieta][par];
00307         } else if (jetAlgo=="AntiKt6TopoGCW") {
00308           JES_factors[ieta][par]     = AntiKt6TopoGCW_JES[ieta][par];
00309           etaCorr_factors[ieta][par] = AntiKt6TopoGCW_EtaCorr[ieta][par];
00310         } else if (jetAlgo=="AntiKt4TowerGCW") {
00311           JES_factors[ieta][par]     = AntiKt4TowerGCW_JES[ieta][par];
00312           etaCorr_factors[ieta][par] = AntiKt4TowerGCW_EtaCorr[ieta][par];
00313         } else if (jetAlgo=="AntiKt6TowerGCW") {
00314           JES_factors[ieta][par]     = AntiKt6TowerGCW_JES[ieta][par];
00315           etaCorr_factors[ieta][par] = AntiKt6TowerGCW_EtaCorr[ieta][par];
00316         } else if (jetAlgo=="AntiKt4LCTopo") {
00317           JES_factors[ieta][par]     = AntiKt4LCTopo_JES[ieta][par];
00318           etaCorr_factors[ieta][par] = AntiKt4LCTopo_EtaCorr[ieta][par];
00319         } else if (jetAlgo=="AntiKt6LCTopo") {
00320           JES_factors[ieta][par]     = AntiKt6LCTopo_JES[ieta][par];
00321           etaCorr_factors[ieta][par] = AntiKt6LCTopo_EtaCorr[ieta][par];
00322         } else 
00323           error("No calibration factors for jetalgo "+jetAlgo);
00324       }
00325     }
00326     
00327     for (int ieta=0; ieta<NOffsetEtaBins; ++ieta) {
00328       for (int par=0; par<NPVmax; ++par) {
00329         if (jetAlgo=="AntiKt6TopoEM") {
00330           useTowers = false;
00331           offset_factors[ieta][par] = jet_offset_new_periodD[ieta][par];
00332         } else if (jetAlgo=="AntiKt4TopoEM") {
00333           useTowers = false;
00334           offset_factors[ieta][par] = jet_offset_R4_new_periodD[ieta][par];
00335         } else if (jetAlgo=="AntiKt6TowerEM") {
00336           useTowers = true;
00337           offset_factors[ieta][par] = tower_offset_new_periodD[ieta][par];
00338         } else if (jetAlgo=="AntiKt4TowerEM") {
00339           useTowers = true;
00340           offset_factors[ieta][par] = tower_offset_new_periodD[ieta][par];
00341         } else if (jetAlgo=="AntiKt4TopoGCW") {
00342           useTowers = false;
00343           offset_factors[ieta][par] = 0.0;
00344         } else if (jetAlgo=="AntiKt6TopoGCW") {
00345           useTowers = false;
00346           offset_factors[ieta][par] = 0.0;
00347         } else if (jetAlgo=="AntiKt4TowerGCW") {
00348           useTowers = false;
00349           offset_factors[ieta][par] = 0.0;
00350         } else if (jetAlgo=="AntiKt6TowerGCW") {
00351           useTowers = false;
00352           offset_factors[ieta][par] = 0.0;
00353         } else if (jetAlgo=="AntiKt4LCTopo") {
00354           useTowers = false;
00355           offset_factors[ieta][par] = 0.0;
00356         } else if (jetAlgo=="AntiKt6LCTopo") {
00357           useTowers = false;
00358           offset_factors[ieta][par] = 0.0;
00359         } else 
00360           error("No offset factors for jetalgo "+jetAlgo);
00361       }
00362     }
00363     
00364   }
00365   
00366 };
00367 
00368 #endif