From 80b1b92b53e80d1e6f2907fd9efb2884c2be6263 Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Tue, 20 Apr 2021 14:28:58 +0200 Subject: [PATCH 1/9] Fixing high-pileup config for residual and adding ability to use p-splines using hists --- .../CalibrationMethods/EtaJESCorrection.h | 12 +- .../JetCalibTools/Root/EtaJESCorrection.cxx | 157 +++++++++++------- .../Root/PUResidual3DCorrection.h | 19 +++ 3 files changed, 130 insertions(+), 58 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h index 879c1368ed4..0b6deb95646 100644 --- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h +++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h @@ -14,6 +14,7 @@ #include #include +#include "TH1.h" #include "JetCalibTools/IJetCalibrationTool.h" #include "JetCalibTools/JetCalibrationToolBase.h" @@ -44,6 +45,7 @@ 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 getSpline(const int etaBin, double E) const; int getEtaBin(double eta_det) const; private: @@ -53,6 +55,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; @@ -66,7 +69,7 @@ class EtaJESCorrection TAxis * m_etaBinAxis; // 90 eta bins, and up to 9 parameter for the pol-fit - const static unsigned int s_nEtaBins=90; + const static unsigned int s_nEtaBins=2; const static unsigned int s_nParMin=7; const static unsigned int s_nParMax=9; unsigned int m_nPar; // number of parameters in config file @@ -81,6 +84,13 @@ 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; + + void loadParameters(const std::string & fileName, const std::string &etajes_name = "etaJes"); + + }; #endif diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index f008dc6910f..15449a33f63 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","")); @@ -88,53 +89,23 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { else { ATH_MSG_FATAL( "You can't apply the mass correction unless you specify ApplyMassCorrection: true in the configuration file!"); return StatusCode::FAILURE; } } + m_config->Print(); 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 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; - } - } - - // Read in jet eta calibration factors - key=Form("EtaCorr.%s_Bin%d",m_jetAlgo.Data(),ieta); - params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); if (params.size()!=m_nPar) { ATH_MSG_FATAL( "Cannot read jet eta calib constants " << key ); return StatusCode::FAILURE; } for (uint ipar=0;iparGetValue(key,"")); + if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} + for (uint ipar=0;iparGetValue(key,"")); @@ -142,19 +113,77 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { for (uint ipar=0;ipar<1;++ipar) m_energyFreezeJES[ieta] = params[ipar]; } - if(m_applyMassCorrection) { - // Read in absolute JMS calibration factors - key=Form("MassCorr.%s_Bin%d",m_jetAlgo.Data(),ieta); - params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); - if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS 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; + } + } } + } + // For splines, only use the spline for the JES calibration + if(m_isSpline){ + std::string fileName = PathResolverFindCalibFile("JetCalibTools/calibrationFactors_response.root"); + loadParameters(fileName, "etaJes"); } return StatusCode::SUCCESS; } +/// Loads the calib constants from histograms in TFile named fileName. +// TODO: Should figure out how to deal with low-pT extrapolation +void EtaJESCorrection::loadParameters(const std::string & fileName, const std::string &etajes_name) +{ + std::unique_ptr tmpF(TFile::Open( fileName.c_str() )); + TList *etajes_l = (TList*) tmpF->Get(etajes_name.c_str()); + + m_etajesFactors.resize( etajes_l->GetSize() ); + if(etajes_l->GetSize() != m_etaBinAxis->GetNbins()){ + ATH_MSG_FATAL("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() << "\t" << s_nEtaBins ); + } + + for(size_t i=0 ; iGetNbins(); i++){ + m_etajesFactors[i].reset((TH1D*)etajes_l->At(i)); + m_etajesFactors[i]->SetDirectory(nullptr); + } +} + + StatusCode EtaJESCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo&) const { xAOD::JetFourMom_t jetStartP4; @@ -190,8 +219,20 @@ 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 { + // Get the factors + int ieta = getEtaBin(eta_det); double E = E_uncorr/m_GeV; // E in GeV + + // 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 = getSpline(ieta, E); + return 1.0/R; + } + //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,16 +250,13 @@ 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]; // Calculate the jet response and then the JES as 1/R - double R = getLogPolN(factors,E); - return 1.0/R; + 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 { @@ -288,6 +326,11 @@ double EtaJESCorrection::getLogPolNSlope(const double *factors, double x) const return y; } +double EtaJESCorrection::getSpline(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 2ced214f398..46d018cc921 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); -- GitLab From f42a746690d64757437f792b94321850d4562cc5 Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Wed, 21 Apr 2021 21:52:04 +0200 Subject: [PATCH 2/9] Fixing a couple of remaining issues --- .../CalibrationMethods/EtaJESCorrection.h | 7 +++-- .../JetCalibTools/Root/EtaJESCorrection.cxx | 31 +++++++++++++------ 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h index 0b6deb95646..6544af764f7 100644 --- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h +++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h @@ -15,6 +15,7 @@ #include #include #include "TH1.h" +#include "TString.h" #include "JetCalibTools/IJetCalibrationTool.h" #include "JetCalibTools/JetCalibrationToolBase.h" @@ -45,7 +46,8 @@ 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 getSpline(const int etaBin, double E) const; + double getSplineCorr(const int etaBin, double E) const; + void loadSplineHists(const TString & fileName, const std::string &etajes_name = "etaJes"); int getEtaBin(double eta_det) const; private: @@ -69,7 +71,7 @@ class EtaJESCorrection TAxis * m_etaBinAxis; // 90 eta bins, and up to 9 parameter for the pol-fit - const static unsigned int s_nEtaBins=2; + const static unsigned int s_nEtaBins=90; const static unsigned int s_nParMin=7; const static unsigned int s_nParMax=9; unsigned int m_nPar; // number of parameters in config file @@ -88,7 +90,6 @@ class EtaJESCorrection // This avoids importing new packages into Athena. std::vector > m_etajesFactors; - void loadParameters(const std::string & fileName, const std::string &etajes_name = "etaJes"); }; diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index 15449a33f63..aee4ff60225 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -156,10 +156,20 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { } } - // For splines, only use the spline for the JES calibration if(m_isSpline){ - std::string fileName = PathResolverFindCalibFile("JetCalibTools/calibrationFactors_response.root"); - loadParameters(fileName, "etaJes"); + TString absoluteJESCalibHists = m_config->GetValue("AbsoluteJES.CalibHists",""); + if(m_dev){ + absoluteJESCalibHists.Remove(0,33); + absoluteJESCalibHists.Insert(0,"JetCalibTools/"); + } + else{ + absoluteJESCalibHists.Insert(14,m_calibAreaTag); + } + TString calibHistFile = PathResolverFindCalibFile(absoluteJESCalibHists.Data()); + + + //std::string fileName = PathResolverFindCalibFile("JetCalibTools/calibrationFactors_response.root"); + loadSplineHists(calibHistFile, "etaJes"); } return StatusCode::SUCCESS; @@ -167,9 +177,9 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { /// Loads the calib constants from histograms in TFile named fileName. // TODO: Should figure out how to deal with low-pT extrapolation -void EtaJESCorrection::loadParameters(const std::string & fileName, const std::string &etajes_name) +void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::string &etajes_name) { - std::unique_ptr tmpF(TFile::Open( fileName.c_str() )); + std::unique_ptr tmpF(TFile::Open( fileName )); TList *etajes_l = (TList*) tmpF->Get(etajes_name.c_str()); m_etajesFactors.resize( etajes_l->GetSize() ); @@ -229,7 +239,7 @@ double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { // The low pT extrapolation doesn't work for the spline yet, so putting this code here. if(m_isSpline){ - double R = getSpline(ieta, E); + double R = getSplineCorr(ieta, E); return 1.0/R; } @@ -253,9 +263,9 @@ double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { // Calculate the jet response and then the JES as 1/R - const double *factors = m_JESFactors[ieta]; - double R = getLogPolN(factors,E); - return 1.0/R; + const double *factors = m_JESFactors[ieta]; + double R = getLogPolN(factors,E); + return 1.0/R; } @@ -284,6 +294,7 @@ double EtaJESCorrection::getLowPtJES(double E_uncorr, double eta_det) const { } double EtaJESCorrection::getEtaCorr(double E_corr, double eta_det) const { + //return 0; int ieta = getEtaBin(eta_det); const double *factors = m_etaCorrFactors[ieta]; @@ -326,7 +337,7 @@ double EtaJESCorrection::getLogPolNSlope(const double *factors, double x) const return y; } -double EtaJESCorrection::getSpline(const int etaBin, double E) const { +double EtaJESCorrection::getSplineCorr(const int etaBin, double E) const { double R = m_etajesFactors[ etaBin ]->Interpolate(E); return R; } -- GitLab From be9bf8df9e9bf7b43a7973c493a5efac12ab5341 Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Wed, 21 Apr 2021 22:01:16 +0200 Subject: [PATCH 3/9] Moving some code around --- .../JetCalibTools/Root/EtaJESCorrection.cxx | 54 +++++++++---------- 1 file changed, 25 insertions(+), 29 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index aee4ff60225..b2a563d0bd0 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -89,30 +89,7 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { else { ATH_MSG_FATAL( "You can't apply the mass correction unless you specify ApplyMassCorrection: true in the configuration file!"); return StatusCode::FAILURE; } } - m_config->Print(); for (uint ieta=0;ieta 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(key,"")); - if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} - for (uint ipar=0;iparGetValue(key,"")); - if (params.size()!=1) { ATH_MSG_FATAL( "Cannot read starting energy for the freezing of JES correction " << key ); return StatusCode::FAILURE; } - for (uint ipar=0;ipar<1;++ipar) m_energyFreezeJES[ieta] = params[ipar]; - } - if(!m_isSpline){ TString key=Form("JES.%s_Bin%d",m_jetAlgo.Data(),ieta); @@ -154,6 +131,30 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { } } } + + // Read in jet eta calibration factors + 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(key,"")); + if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} + for (uint ipar=0;iparGetValue(key,"")); + if (params.size()!=1) { ATH_MSG_FATAL( "Cannot read starting energy for the freezing of JES correction " << key ); return StatusCode::FAILURE; } + for (uint ipar=0;ipar<1;++ipar) m_energyFreezeJES[ieta] = params[ipar]; + } + + } if(m_isSpline){ @@ -166,9 +167,6 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { absoluteJESCalibHists.Insert(14,m_calibAreaTag); } TString calibHistFile = PathResolverFindCalibFile(absoluteJESCalibHists.Data()); - - - //std::string fileName = PathResolverFindCalibFile("JetCalibTools/calibrationFactors_response.root"); loadSplineHists(calibHistFile, "etaJes"); } @@ -229,7 +227,6 @@ 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 { - // Get the factors int ieta = getEtaBin(eta_det); double E = E_uncorr/m_GeV; // E in GeV @@ -263,10 +260,10 @@ double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { // 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 { @@ -294,7 +291,6 @@ double EtaJESCorrection::getLowPtJES(double E_uncorr, double eta_det) const { } double EtaJESCorrection::getEtaCorr(double E_corr, double eta_det) const { - //return 0; int ieta = getEtaBin(eta_det); const double *factors = m_etaCorrFactors[ieta]; -- GitLab From f29a4f0b4bc8516187a98b6c44d29ed47233fda7 Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Wed, 21 Apr 2021 22:03:41 +0200 Subject: [PATCH 4/9] Whitespace changes --- .../Jet/JetCalibTools/Root/EtaJESCorrection.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index b2a563d0bd0..3b0f2f33e12 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -140,11 +140,11 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { for (uint ipar=0;iparGetValue(key,"")); - if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} - for (uint ipar=0;iparGetValue(key,"")); + if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} + for (uint ipar=0;ipar Date: Wed, 21 Apr 2021 22:05:28 +0200 Subject: [PATCH 5/9] Whitespace changes --- .../Jet/JetCalibTools/Root/EtaJESCorrection.cxx | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index 3b0f2f33e12..866b067023e 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -139,14 +139,6 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { if (params.size()!=m_nPar) { ATH_MSG_FATAL( "Cannot read jet eta calib constants " << key ); return StatusCode::FAILURE; } for (uint ipar=0;iparGetValue(key,"")); - if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} - for (uint ipar=0;iparGetValue(key,"")); @@ -154,6 +146,13 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { for (uint ipar=0;ipar<1;++ipar) m_energyFreezeJES[ieta] = params[ipar]; } + if(m_applyMassCorrection) { + // Read in absolute JMS calibration factors + key=Form("MassCorr.%s_Bin%d",m_jetAlgo.Data(),ieta); + params = JetCalibUtils::VectorizeD(m_config->GetValue(key,"")); + if (params.size()!=m_nPar) {ATH_MSG_FATAL( "Cannot read JMS calib constants " << key ); return StatusCode::FAILURE;} + for (uint ipar=0;ipar Date: Wed, 21 Apr 2021 23:33:48 +0200 Subject: [PATCH 6/9] Some basic implementation for low-pT extrapolation --- .../JetCalibTools/Root/EtaJESCorrection.cxx | 66 +++++++++++++++---- 1 file changed, 52 insertions(+), 14 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index 866b067023e..b4dc8ccddb7 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -92,8 +92,8 @@ 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; } @@ -167,13 +167,38 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { } 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_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; + } + } } return StatusCode::SUCCESS; } /// Loads the calib constants from histograms in TFile named fileName. -// TODO: Should figure out how to deal with low-pT extrapolation void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::string &etajes_name) { std::unique_ptr tmpF(TFile::Open( fileName )); @@ -181,13 +206,14 @@ void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::stri m_etajesFactors.resize( etajes_l->GetSize() ); if(etajes_l->GetSize() != m_etaBinAxis->GetNbins()){ - ATH_MSG_FATAL("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() << "\t" << s_nEtaBins ); + ATH_MSG_FATAL("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() ); } for(size_t i=0 ; iGetNbins(); i++){ m_etajesFactors[i].reset((TH1D*)etajes_l->At(i)); m_etajesFactors[i]->SetDirectory(nullptr); } + tmpF->Close(); } @@ -226,19 +252,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 { - int ieta = getEtaBin(eta_det); - double E = E_uncorr/m_GeV; // E in GeV - // 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; - } - //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){ @@ -256,6 +271,15 @@ double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { } } + int ieta = getEtaBin(eta_det); + // 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 @@ -267,6 +291,7 @@ double EtaJESCorrection::getJES(double E_uncorr, double eta_det) const { 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); @@ -332,6 +357,19 @@ 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; -- GitLab From 2623f89f27a46017f305963a67dc851e1ae916e6 Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Thu, 22 Apr 2021 09:51:36 +0200 Subject: [PATCH 7/9] Some minor fixes --- .../JetCalibTools/CalibrationMethods/EtaJESCorrection.h | 1 + Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h index 6544af764f7..ed5f97e8dc1 100644 --- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h +++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/EtaJESCorrection.h @@ -47,6 +47,7 @@ class EtaJESCorrection 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; diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index b4dc8ccddb7..e55c381496c 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -181,7 +181,7 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { 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]); + if(std::abs(etaBins[ieta]) < m_etaSecondaryminPt_JES) Ecutoff = m_minPt_JES*cosh(etaBins[ieta]); else{ Ecutoff = m_secondaryminPt_JES*cosh(etaBins[ieta]);} } const double Rcutoff = getSplineCorr(ieta, Ecutoff); @@ -209,7 +209,7 @@ void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::stri ATH_MSG_FATAL("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() ); } - for(size_t i=0 ; iGetNbins(); i++){ + for(int i=0 ; iGetNbins(); i++){ m_etajesFactors[i].reset((TH1D*)etajes_l->At(i)); m_etajesFactors[i]->SetDirectory(nullptr); } -- GitLab From 6cab0d52c9a332f800e12003fc55341839b6d1ea Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Thu, 22 Apr 2021 17:18:34 +0200 Subject: [PATCH 8/9] Small fix for eta binning check --- Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index e55c381496c..1a4dfb09ef1 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -205,7 +205,7 @@ void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::stri TList *etajes_l = (TList*) tmpF->Get(etajes_name.c_str()); m_etajesFactors.resize( etajes_l->GetSize() ); - if(etajes_l->GetSize() != m_etaBinAxis->GetNbins()){ + if(etajes_l->GetSize() != m_etaBinAxis->GetNbins()+1){ ATH_MSG_FATAL("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() ); } -- GitLab From 6c2c90537899100da2aa623165f78edd8a8e4e03 Mon Sep 17 00:00:00 2001 From: Jennifer Kathryn Roloff Date: Thu, 22 Apr 2021 17:31:26 +0200 Subject: [PATCH 9/9] Adding a couple extra style fixes --- .../Jet/JetCalibTools/Root/EtaJESCorrection.cxx | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx index 1a4dfb09ef1..4d9e7cb3a71 100644 --- a/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx +++ b/Reconstruction/Jet/JetCalibTools/Root/EtaJESCorrection.cxx @@ -186,7 +186,7 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { } const double Rcutoff = getSplineCorr(ieta, Ecutoff); const double Slope = getSplineSlope(ieta, 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] ); + if(Slope > 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; @@ -202,15 +202,16 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) { void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::string &etajes_name) { std::unique_ptr tmpF(TFile::Open( fileName )); - TList *etajes_l = (TList*) tmpF->Get(etajes_name.c_str()); + 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_FATAL("Do not have the correct number of eta bins for " << fileName << "\t" << etajes_name << "\t" << etajes_l->GetSize() ); + 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((TH1D*)etajes_l->At(i)); + m_etajesFactors[i].reset(static_cast(etajes_l->At(i))); m_etajesFactors[i]->SetDirectory(nullptr); } tmpF->Close(); -- GitLab