Commit 4ff95d48 authored by Nils Erik Krumnack's avatar Nils Erik Krumnack
Browse files

Merge branch '21.2-JetCalibTool-fixes' into '21.2'

Implementations for precision JES calibration

See merge request atlas/athena!47912
parents e5f4f411 2421e20c
/*
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#ifndef JETCALIBTOOLS_ETAJESCORRECTION_H
......@@ -16,6 +16,7 @@
#include <TAxis.h>
#include "TH1.h"
#include "TString.h"
#include "TList.h"
#include "JetCalibTools/IJetCalibrationTool.h"
#include "JetCalibTools/JetCalibrationToolBase.h"
......@@ -89,7 +90,7 @@ class EtaJESCorrection
// When using p-splines, the calibrations are stored as a finely binned TH1 for simplicity.
// This avoids importing new packages into Athena.
std::vector<std::unique_ptr<TH1D> > m_etajesFactors;
std::vector<std::unique_ptr<TH1> > m_etajesFactors;
......
......@@ -13,7 +13,7 @@
#include <map>
#include <tuple>
#include "TH1.h"
#include "TList.h"
//xAOD EDM classes
#include "xAODEventInfo/EventInfo.h"
......@@ -154,7 +154,7 @@ class GlobalNNCalibration : public ::JetCalibrationToolBase {
std::vector<std::unique_ptr<lwt::LightweightGraph> > m_lwnns;
std::vector<std::unique_ptr<TH1D> > m_ptCorrFactors;
std::vector<std::unique_ptr<TH1> > m_ptCorrFactors;
std::vector<double> m_nnEtaBins;
std::vector<double> m_closureEtaBins;
std::vector<TString> m_NNInputs;
......
......@@ -100,6 +100,7 @@ class GlobalSequentialCorrection
bool m_caloBased;
bool m_pTResponseRequirementOff;
bool m_nTrkwTrk_4PFlow;
bool m_TileGap3_4PFlow;
double m_turnOffStartingpT, m_turnOffEndpT;
};
......
/*
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
#include "JetCalibTools/CalibrationMethods/EtaJESCorrection.h"
......@@ -92,7 +92,7 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) {
for (uint ieta=0;ieta<etaBins.size()-1;++ieta) {
if(!m_isSpline){
// Read in absolute JES calibration factors
// Read in absolute JES calibration factors
TString key=Form("JES.%s_Bin%d",m_jetAlgo.Data(),ieta);
std::vector<double> params = JetCalibUtils::VectorizeD(m_config->GetValue(key,""));
m_nPar = params.size();
......@@ -105,38 +105,41 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) {
//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]);
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] );
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;
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 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
TString key=Form("EtaCorr.%s_Bin%d",m_jetAlgo.Data(),ieta);
std::vector<double> 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; }
if (!m_isSpline && params.size()!=m_nPar) { ATH_MSG_FATAL( "Cannot read jet eta calib constants " << key ); return StatusCode::FAILURE; }
if (m_isSpline){
if(params.size() == 0) { ATH_MSG_FATAL( "Cannot read jet eta calib constants " << key ); return StatusCode::FAILURE; }
m_nPar = params.size();
}
for (uint ipar=0;ipar<m_nPar;++ipar) m_etaCorrFactors[ieta][ipar] = params[ipar];
if(m_freezeJESatHighE){ // Read starting energy values to freeze JES correction
......@@ -202,7 +205,7 @@ StatusCode EtaJESCorrection::initializeTool(const std::string&) {
void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::string &etajes_name)
{
std::unique_ptr<TFile> tmpF(TFile::Open( fileName ));
TList *etajes_l = static_cast<TList*>( tmpF->Get(etajes_name.c_str()));
TList *etajes_l = dynamic_cast<TList*>( tmpF->Get(etajes_name.c_str()));
m_etajesFactors.resize( etajes_l->GetSize() );
if(etajes_l->GetSize() != m_etaBinAxis->GetNbins()+1){
......@@ -210,8 +213,7 @@ void EtaJESCorrection::loadSplineHists(const TString & fileName, const std::stri
}
for(int i=0 ; i<m_etaBinAxis->GetNbins(); i++){
//m_etajesFactors[i].reset((TH1D*)etajes_l->At(i));
m_etajesFactors[i].reset(static_cast<TH1D*>(etajes_l->At(i)));
m_etajesFactors[i].reset(dynamic_cast<TH1*>(etajes_l->At(i)));
m_etajesFactors[i]->SetDirectory(nullptr);
}
tmpF->Close();
......
/*
Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
/* ***********************************************************************************\
......@@ -165,7 +165,7 @@ void GlobalNNCalibration::loadSplineHists(const TString & fileName, const std::s
}
for(unsigned int i=0 ; i<m_closureEtaBins.size()-1; i++){
m_ptCorrFactors[i].reset(dynamic_cast<TH1D*>(ptCorr_l->At(i)));
m_ptCorrFactors[i].reset(dynamic_cast<TH1*>(ptCorr_l->At(i)));
m_ptCorrFactors[i]->SetDirectory(nullptr);
}
tmpF->Close();
......
......@@ -81,6 +81,10 @@ StatusCode GlobalSequentialCorrection::initializeTool(const std::string&) {
// For AFII calibrations, EM3 correction should be applied up to |eta|=3.2
m_EM3MaxEtaBin = m_config->GetValue("EM3MaxEtaBin", 35);
// In release 21 (precision), TileGap3 is used for PFlow jets for the precision recommendations
// The default is set to false to maintain backwards compatibility
m_TileGap3_4PFlow = m_config->GetValue("TileGap3_4PFlow", false);
if ( !m_config ) { ATH_MSG_FATAL("Config file not specified. Aborting."); return StatusCode::FAILURE; }
if ( m_jetAlgo.EqualTo("") ) { ATH_MSG_FATAL("No jet algorithm specified. Aborting."); return StatusCode::FAILURE; }
......@@ -152,13 +156,19 @@ StatusCode GlobalSequentialCorrection::initializeTool(const std::string&) {
else if ( depthString.Contains("ChargedFraction") ) m_depth = ApplyChargedFraction;
else { ATH_MSG_FATAL("depthString flag not properly set, please check your config file."); return StatusCode::FAILURE; }
} else {
if ( depthString.Contains("PunchThrough") || depthString.Contains("Full") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk | ApplytrackWIDTH | ApplyPunchThrough;
else if ( depthString.Contains("trackWIDTH") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk | ApplytrackWIDTH;
else if ( depthString.Contains("nTrk") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk;
else if ( depthString.Contains("EM3") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3;
else if ( depthString.Contains("Tile0") ) m_depth = ApplyChargedFraction | ApplyTile0;
else if ( depthString.Contains("ChargedFraction") ) m_depth = ApplyChargedFraction;
else { ATH_MSG_FATAL("depthString flag not properly set, please check your config file."); return StatusCode::FAILURE; }
if( m_TileGap3_4PFlow ){
if (depthString.Contains("TileGap3") || depthString.Contains("Full")) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk | ApplytrackWIDTH | ApplyPunchThrough | ApplyTileGap3;
else { ATH_MSG_FATAL("depthString flag not properly set, please check your config file."); return StatusCode::FAILURE; }
}
else{
if ( depthString.Contains("PunchThrough") || depthString.Contains("Full") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk | ApplytrackWIDTH | ApplyPunchThrough;
else if ( depthString.Contains("trackWIDTH") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk | ApplytrackWIDTH;
else if ( depthString.Contains("nTrk") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3 | ApplynTrk;
else if ( depthString.Contains("EM3") ) m_depth = ApplyChargedFraction | ApplyTile0 | ApplyEM3;
else if ( depthString.Contains("Tile0") ) m_depth = ApplyChargedFraction | ApplyTile0;
else if ( depthString.Contains("ChargedFraction") ) m_depth = ApplyChargedFraction;
else { ATH_MSG_FATAL("depthString flag not properly set, please check your config file."); return StatusCode::FAILURE; }
}
}
}
......@@ -267,6 +277,10 @@ StatusCode GlobalSequentialCorrection::initializeTool(const std::string&) {
ATH_MSG_FATAL("Vector of PunchThrough histograms may be empty. Please check your GSCFactors file: " << GSCFile);
return StatusCode::FAILURE;
}
else if ( (m_depth & ApplyTileGap3) && m_respFactorsTileGap3.size() < 3 ) {
ATH_MSG_FATAL("Vector of TileGap3 histograms may be empty. Please check your GSCFactors file: " << GSCFile);
return StatusCode::FAILURE;
}
else ATH_MSG_INFO("GSC Tool has been initialized with binning and eta fit factors from: " << fileName);
}
return StatusCode::SUCCESS;
......@@ -418,12 +432,20 @@ double GlobalSequentialCorrection::getGSCCorrection(xAOD::JetFourMom_t jetP4, do
if (m_depth & ApplyEM3) Corr*=1./getEM3Response(jetP4.pt()/m_GeV*Corr, etabin, EM3);
if ( m_nTrkwTrk_4PFlow && (m_depth & ApplynTrk) ) Corr*=1./getNTrkResponse(jetP4.pt()/m_GeV*Corr, etabin, nTrk);
if ( m_nTrkwTrk_4PFlow && (m_depth & ApplytrackWIDTH) ) Corr*=1./getTrackWIDTHResponse(jetP4.pt()/m_GeV*Corr,etabin,trackWIDTH);
if ( jetP4.pt() < m_punchThroughMinPt ) return Corr; //Applying punch through correction to low pT jets introduces a bias, default threshold is 50 GeV
if ( jetP4.pt() < m_punchThroughMinPt ) {
if (m_TileGap3_4PFlow && (m_depth & ApplyTileGap3) ) {
Corr*=1/getTileGap3Response(jetP4.pt()/m_GeV*Corr,etabin, TileGap3);
}
return Corr; //Applying punch through correction to low pT jets introduces a bias, default threshold is 50 GeV
}
//eta binning for the punch through correction differs from the rest of the GSC, so the eta bin is determined in the GetPunchThroughResponse method
else if (m_depth & ApplyPunchThrough) {
jetP4*=Corr; //The punch through correction is binned in E instead of pT, so we determine E from the corrected jet here
Corr*=1/getPunchThroughResponse(jetP4.e()/m_GeV,eta,Nsegments);
}
if (m_TileGap3_4PFlow && (m_depth & ApplyTileGap3) ) {
Corr*=1/getTileGap3Response(jetP4.pt()/m_GeV*Corr,etabin, TileGap3);
}
}
return Corr;
}
......@@ -540,6 +562,12 @@ StatusCode GlobalSequentialCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInf
}
}
if(m_TileGap3_4PFlow){
if (m_depth & ApplyTileGap3) {
TG3 = (samplingFrac[17])/jetE_constitscale;
}
}
xAOD::JetFourMom_t calibP4 = jetStartP4*getGSCCorrection( jetStartP4, fabs(detectorEta), trackWIDTHPV0, nTrkPV0, Tile0, EM3, Nsegments, ChargedFraction, caloWIDTH, N90Constituents, TG3);
//Transfer calibrated jet properties to the Jet object
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment