diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h index 879c1368ed4ff1960216f9b5490ab21ccce609cf..ed5f97e8dc1d4bf1cd381a27d822f08eccfdae47 100644 --- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h +++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h @@ -14,6 +14,8 @@ #include #include +#include "TH1.h" +#include "TString.h" #include "JetCalibTools/IJetCalibrationTool.h" #include "JetCalibTools/JetCalibrationToolBase.h" @@ -44,6 +46,9 @@ class EtaJESCorrection double getMassCorr(double E_corr, double eta_det) const; double getLogPolN(const double *factors, double x) const; double getLogPolNSlope(const double *factors, double x) const; + double getSplineCorr(const int etaBin, double E) const; + double getSplineSlope(const int ieta, const double minE) const; + void loadSplineHists(const TString & fileName, const std::string &etajes_name = "etaJes"); int getEtaBin(double eta_det) const; private: @@ -53,6 +58,7 @@ class EtaJESCorrection bool m_mass; bool m_dev; bool m_freezeJESatHighE; + bool m_isSpline; TString m_jesDesc; double m_minPt_JES, m_minPt_EtaCorr, m_maxE_EtaCorr; @@ -81,6 +87,12 @@ class EtaJESCorrection double m_JMSFactors[s_nEtaBins][s_nParMax]; double m_energyFreezeJES[s_nEtaBins]; + // When using p-splines, the calibrations are stored as a finely binned TH1 for simplicity. + // This avoids importing new packages into Athena. + std::vector > m_etajesFactors; + + + }; #endif diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index f008dc6910fc70692b547b61b89f1ef71f8c51d3..4d9e7cb3a71145f8a161db3ddd900deb60872a2e 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -7,7 +7,7 @@ EtaJESCorrection::EtaJESCorrection() : JetCalibrationToolBase::JetCalibrationToolBase("EtaJESCorrection::EtaJESCorrection"), - m_config(NULL), m_jetAlgo(""), m_calibAreaTag(""), m_mass(false), m_dev(false), + m_config(NULL), m_jetAlgo(""), m_calibAreaTag(""), m_mass(false), m_dev(false), m_isSpline(false), m_minPt_JES(10), m_minPt_EtaCorr(8), m_maxE_EtaCorr(2500), m_lowPtExtrap(0), m_lowPtMinR(0.25), m_etaBinAxis(NULL) @@ -15,7 +15,7 @@ EtaJESCorrection::EtaJESCorrection() EtaJESCorrection::EtaJESCorrection(const std::string& name) : JetCalibrationToolBase::JetCalibrationToolBase( name ), - m_config(NULL), m_jetAlgo(""), m_calibAreaTag(""), m_mass(false), m_dev(false), + m_config(NULL), m_jetAlgo(""), m_calibAreaTag(""), m_mass(false), m_dev(false), m_isSpline(false), m_minPt_JES(10), m_minPt_EtaCorr(8), m_maxE_EtaCorr(2500), m_lowPtExtrap(0), m_lowPtMinR(0.25), m_etaBinAxis(NULL) @@ -23,7 +23,7 @@ EtaJESCorrection::EtaJESCorrection(const std::string& name) EtaJESCorrection::EtaJESCorrection(const std::string& name, TEnv * config, TString jetAlgo, TString calibAreaTag, bool mass, bool dev) : JetCalibrationToolBase::JetCalibrationToolBase( name ), - m_config(config), m_jetAlgo(jetAlgo), m_calibAreaTag(calibAreaTag), m_mass(mass), m_dev(dev), + m_config(config), m_jetAlgo(jetAlgo), m_calibAreaTag(calibAreaTag), m_mass(mass), m_dev(dev), m_isSpline(false), m_minPt_JES(10), m_minPt_EtaCorr(8), m_maxE_EtaCorr(2500), m_lowPtExtrap(0), m_lowPtMinR(0.25), m_etaBinAxis(NULL) @@ -67,6 +67,7 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { m_secondaryminPt_JES = m_config->GetValue(m_jetAlgo+".SecondaryMinPtForETAJES",7); // Freeze JES correction at maximum values of energy for each eta bin m_freezeJESatHighE = m_config->GetValue(m_jetAlgo+".FreezeJEScorrectionatHighE", false); + m_isSpline = m_config->GetValue(m_jetAlgo+".isSpline", false); // From mswiatlo, help from dag: variable eta binning std::vector etaBins = JetCalibUtils::VectorizeD(m_config->GetValue("JES.EtaBins","")); @@ -89,49 +90,52 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { } for (uint ieta=0;ieta params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); - m_nPar = params.size(); - if (m_nPars_nParMax) { ATH_MSG_FATAL( "Cannot read JES calib constants " << key ); return StatusCode::FAILURE; } - for (uint ipar=0;ipar params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); + m_nPar = params.size(); + if (m_nPars_nParMax) { ATH_MSG_FATAL( "Cannot read JES calib constants " << key ); return StatusCode::FAILURE; } + for (uint ipar=0;ipar 0) { - //Calculate the slope of the response curve at the minPt for each eta bin - //Used in the GetLowPtJES method when Pt < minPt - const double *factors = m_JESFactors[ieta]; - double Ecutoff; - if(!m_useSecondaryminPt_JES) Ecutoff = m_minPt_JES*cosh(etaBins[ieta]); - else { - if(fabs(etaBins[ieta]) < m_etaSecondaryminPt_JES) Ecutoff = m_minPt_JES*cosh(etaBins[ieta]); - else{ Ecutoff = m_secondaryminPt_JES*cosh(etaBins[ieta]);} - } - const double Rcutoff = getLogPolN(factors,Ecutoff); - const double Slope = getLogPolNSlope(factors,Ecutoff); - if(Slope > Rcutoff/Ecutoff) ATH_MSG_FATAL("Slope of calibration curve at minimum ET is too steep for the JES factors of etabin " << ieta << ", eta = " << etaBins[ieta] ); - - m_JES_MinPt_E[ieta] = Ecutoff; - m_JES_MinPt_R[ieta] = Rcutoff; - m_JES_MinPt_Slopes[ieta] = Slope; - - //Calculate the parameters for a 2nd order polynomial extension to the calibration curve below minimum ET - //Used in the GetLowPtJES method when Pt < minPt - if(m_lowPtExtrap == 2) { - const double h = m_lowPtMinR; - const double Param1 = (2/Ecutoff)*(Rcutoff-h)-Slope; - const double Param2 = (0.5/Ecutoff)*(Slope-Param1); - //Slope of the calibration curve should always be positive - if( Param1 < 0 || Param1 + 2*Param2*Ecutoff < 0) ATH_MSG_FATAL("Polynomial extension to calibration curve below minimum ET is not monotonically increasing for etabin " << ieta << ", eta = " << etaBins[ieta] ); - m_JES_MinPt_Param1[ieta] = Param1; - m_JES_MinPt_Param2[ieta] = Param2; - } + //Calculate the slope of the response curve at the minPt for each eta bin + //Used in the GetLowPtJES method when Pt < minPt + const double *factors = m_JESFactors[ieta]; + double Ecutoff; + if(!m_useSecondaryminPt_JES) Ecutoff = m_minPt_JES*cosh(etaBins[ieta]); + else { + if(fabs(etaBins[ieta]) < m_etaSecondaryminPt_JES) Ecutoff = m_minPt_JES*cosh(etaBins[ieta]); + else{ Ecutoff = m_secondaryminPt_JES*cosh(etaBins[ieta]);} + } + const double Rcutoff = getLogPolN(factors,Ecutoff); + const double Slope = getLogPolNSlope(factors,Ecutoff); + if(Slope > Rcutoff/Ecutoff) ATH_MSG_FATAL("Slope of calibration curve at minimum ET is too steep for the JES factors of etabin " << ieta << ", eta = " << etaBins[ieta] ); + + m_JES_MinPt_E[ieta] = Ecutoff; + m_JES_MinPt_R[ieta] = Rcutoff; + m_JES_MinPt_Slopes[ieta] = Slope; + + //Calculate the parameters for a 2nd order polynomial extension to the calibration curve below minimum ET + //Used in the GetLowPtJES method when Pt < minPt + if(m_lowPtExtrap == 2) { + const double h = m_lowPtMinR; + const double Param1 = (2/Ecutoff)*(Rcutoff-h)-Slope; + const double Param2 = (0.5/Ecutoff)*(Slope-Param1); + //Slope of the calibration curve should always be positive + if( Param1 < 0 || Param1 + 2*Param2*Ecutoff < 0) ATH_MSG_FATAL("Polynomial extension to calibration curve below minimum ET is not monotonically increasing for etabin " << ieta << ", eta = " << etaBins[ieta] ); + m_JES_MinPt_Param1[ieta] = Param1; + m_JES_MinPt_Param2[ieta] = Param2; + } } + } // Read in jet eta calibration factors - key=Form("EtaCorr.%s_Bin%d",m_jetAlgo.Data(),ieta); - params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); + TString key=Form("EtaCorr.%s_Bin%d",m_jetAlgo.Data(),ieta); + std::vector params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); + m_nPar = params.size(); if (params.size()!=m_nPar) { ATH_MSG_FATAL( "Cannot read jet eta calib constants " << key ); return StatusCode::FAILURE; } for (uint ipar=0;iparGetValue("AbsoluteJES.CalibHists",""); + if(m_dev){ + absoluteJESCalibHists.Remove(0,33); + absoluteJESCalibHists.Insert(0,"JetCalibTools/"); + } + else{ + absoluteJESCalibHists.Insert(14,m_calibAreaTag); + } + TString calibHistFile = PathResolverFindCalibFile(absoluteJESCalibHists.Data()); + loadSplineHists(calibHistFile, "etaJes"); + + //Protections for high order extrapolation methods at low Et (Et < _minPt_JES) + if(m_lowPtExtrap != 1) { + ATH_MSG_ERROR("Only linear extrapolations are supported for p-splines currently. Please change the config file to reflect this"); + return StatusCode::FAILURE; + } + + if(m_lowPtExtrap > 0) { + for (uint ieta=0;ieta Rcutoff/Ecutoff) ATH_MSG_WARNING("Slope of calibration curve at minimum ET is too steep for the JES factors of etabin " << ieta << ", eta = " << etaBins[ieta] ); + + m_JES_MinPt_E[ieta] = Ecutoff; + m_JES_MinPt_R[ieta] = Rcutoff; + m_JES_MinPt_Slopes[ieta] = Slope; + } + } + } + return StatusCode::SUCCESS; } +/// Loads the calib constants from histograms in TFile named fileName. +void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::string &etajes_name) +{ + std::unique_ptr tmpF(TFile::Open( fileName )); + TList *etajes_l = static_cast( tmpF->Get(etajes_name.c_str())); + + m_etajesFactors.resize( etajes_l->GetSize() ); + if(etajes_l->GetSize() != m_etaBinAxis->GetNbins()+1){ + ATH_MSG_WARNING("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() ); + } + + for(int i=0 ; iGetNbins(); i++){ + //m_etajesFactors[i].reset((TH1D*)etajes_l->At(i)); + m_etajesFactors[i].reset(static_cast(etajes_l->At(i))); + m_etajesFactors[i]->SetDirectory(nullptr); + } + tmpF->Close(); +} + + StatusCode EtaJESCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const { xAOD::JetFourMom_t jetStartP4; @@ -190,8 +253,8 @@ StatusCode EtaJESCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const // E_uncorr is the EM-scale, or LC-scale jet energy // eta_det is the eta of the raw, constituent-scale jet (constscale_eta) double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { - double E = E_uncorr/m_GeV; // E in GeV + //Check if the Pt goes below the minimum value, if so use the special GetLowPtJES method if(m_useSecondaryminPt_JES){ if(fabs(eta_det) < m_etaSecondaryminPt_JES && E/cosh(eta_det) < m_minPt_JES){ @@ -209,20 +272,27 @@ double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { } } - // Get the factors int ieta = getEtaBin(eta_det); - const double *factors = m_JESFactors[ieta]; - // Freeze correction if(m_freezeJESatHighE && E>m_energyFreezeJES[ieta] && m_energyFreezeJES[ieta]!=-1) E = m_energyFreezeJES[ieta]; + + // The low pT extrapolation doesn't work for the spline yet, so putting this code here. + if(m_isSpline){ + double R = getSplineCorr(ieta, E); + return 1.0/R; + } + // Calculate the jet response and then the JES as 1/R + // Get the factors + const double *factors = m_JESFactors[ieta]; double R = getLogPolN(factors,E); return 1.0/R; } double EtaJESCorrection::getLowPtJES(double E_uncorr, double eta_det) const { int ieta = getEtaBin(eta_det); + if (m_lowPtExtrap == 0) { const double *factors = m_JESFactors[ieta]; double E = m_minPt_JES*cosh(eta_det); @@ -288,6 +358,24 @@ double EtaJESCorrection::getLogPolNSlope(const double *factors, double x) const return y; } +double EtaJESCorrection::getSplineSlope(const int ieta, const double minE) const { + // Don't want to use interpolation here, so instead just use the values at the bin centers near the cutoff + int minBin = m_etajesFactors[ieta]->FindBin(minE); + + double rFirst = m_etajesFactors[ ieta ]->GetBinContent(minBin); + double rSecond = m_etajesFactors[ ieta ]->GetBinContent(minBin+1); + double binWidth = m_etajesFactors[ ieta ]->GetBinCenter(minBin+1) - m_etajesFactors[ ieta ]->GetBinCenter(minBin); + double slope = (rSecond - rFirst) / binWidth; + + return slope; +} + + +double EtaJESCorrection::getSplineCorr(const int etaBin, double E) const { + double R = m_etajesFactors[ etaBin ]->Interpolate(E); + return R; +} + int EtaJESCorrection::getEtaBin(double eta_det) const { int bin = m_etaBinAxis->FindBin(eta_det); if (bin<=0) return 0; diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h index 2ced214f398173ca66fc36830edb591bcdffbcdd..46d018cc9212eba1d4ff0614214054c01369a947 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h +++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h @@ -130,6 +130,14 @@ namespace PUCorrection { m_Dptp1_vs_eta.reset( (TH1F*) paramDelta_l->At(1) ) ; m_Dptp1_vs_eta->SetDirectory(nullptr); setupClosestNonEmptyBins(); + + // Finding the maximum values of mu and NPV for each eta bin, + // which is used to avoid using overflow bins in the residual calculation. + // The bin center is used to avoid any effects on the boundaries. + for(size_t i=0 ; i<(etaBins_v->size()-1); i++){ + m_maxMu.push_back(m_3Dp0_vs_muNPV[i]->GetXaxis()->GetBinCenter(m_3Dp0_vs_muNPV[i]->GetNbinsX())); + m_maxNPV.push_back(m_3Dp0_vs_muNPV[i]->GetYaxis()->GetBinCenter(m_3Dp0_vs_muNPV[i]->GetNbinsY())); + } } @@ -193,6 +201,8 @@ namespace PUCorrection { std::unique_ptr m_Dptp0_vs_eta=nullptr; std::unique_ptr m_Dptp1_vs_eta=nullptr; + std::vector m_maxMu; + std::vector m_maxNPV; // float m_maxPt=170.0 ; // GeV !! @@ -226,6 +236,15 @@ namespace PUCorrection { float correction3D_noextrap(float pt, float eta , float mu, int NPV) const { int muNPVbin = m_ref3DHisto->FindBin(mu, NPV); int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1; + + // Don't want to use the overflow bins, so using this instead + if(mu > m_maxMu[etaBin]){ + mu=m_maxMu[etaBin]; + } + if(NPV > m_maxNPV[etaBin]){ + NPV=m_maxNPV[etaBin]; + } + float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin); float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);