From c6bdd1fd7fe6c9056766f2d1dc0f23219b8f4c90 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Tue, 7 Apr 2020 11:10:08 +0200
Subject: [PATCH 1/9] Adding implementation of 3D residual PU correction

---
 .../CalibrationMethods/JetPileupCorrection.h  |  11 +-
 .../Root/JetPileupCorrection.cxx              |  36 ++-
 .../Root/PUResidual3DCorrection.h             | 237 ++++++++++++++++++
 3 files changed, 276 insertions(+), 8 deletions(-)
 create mode 100644 Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h

diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h
index 90de31ac3a77..a91f3bdade9d 100644
--- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef JETCALIBTOOLS_JETPILEUPCORRECTION_H
@@ -18,6 +18,10 @@
 #include "JetCalibTools/JetCalibrationToolBase.h"
 #include "JetCalibTools/CalibrationMethods/ResidualOffsetCorrection.h"
 
+namespace PUCorrection {
+  struct PU3DCorrectionHelper;
+}
+
 class JetPileupCorrection 
   : virtual public ::IJetCalibrationTool,
     virtual public ::JetCalibrationToolBase
@@ -51,9 +55,14 @@ class JetPileupCorrection
   bool m_doNJetOnly;
   bool m_doSequentialResidual;
 
+  bool m_do3Dcorrection;
+
+  
   bool m_useFull4vectorArea;
   ResidualOffsetCorrection * m_residualOffsetCorr;
 
+  PUCorrection::PU3DCorrectionHelper *m_residual3DCorr=nullptr;
+  
   bool m_doOnlyResidual;
   
   std::string m_originScale;
diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
index 4e3ef3fd6921..aa49d5defbf5 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
@@ -1,8 +1,9 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "JetCalibTools/CalibrationMethods/JetPileupCorrection.h"
+#include "PUResidual3DCorrection.h"
 
 JetPileupCorrection::JetPileupCorrection()
   : JetCalibrationToolBase::JetCalibrationToolBase("JetPileupCorrection::JetPileupCorrection"),
@@ -37,7 +38,7 @@ JetPileupCorrection::JetPileupCorrection(const std::string& name, TEnv * config,
 JetPileupCorrection::~JetPileupCorrection() {
 
   if(m_residualOffsetCorr) delete m_residualOffsetCorr;
-
+  if(m_residual3DCorr) delete m_residual3DCorr;
 }
 
 
@@ -58,6 +59,8 @@ StatusCode JetPileupCorrection::initializeTool(const std::string& name) {
   m_doSequentialResidual = m_config->GetValue("DoSequentialResidual",false); // first mu and then NPV/NJet corrections
   bool useNjet           = m_config->GetValue("OffsetCorrection.UseNjet", false);
 
+  m_do3Dcorrection       = m_config->GetValue("Do3DCorrection", false);
+  
   if(m_doSequentialResidual) ATH_MSG_DEBUG("The pileup residual calibrations will be applied sequentially.");
   else                       ATH_MSG_DEBUG("The pileup residual calibrations will be applied simultaneously (default).");
   if(m_doMuOnly)             ATH_MSG_INFO("Only the pileup mu-based calibration will be applied.");
@@ -91,6 +94,12 @@ StatusCode JetPileupCorrection::initializeTool(const std::string& name) {
     return StatusCode::FAILURE;
   }
 
+  if(m_do3Dcorrection && (m_doSequentialResidual || m_doMuOnly || m_doNPVOnly || m_doNJetOnly ) ){
+    ATH_MSG_FATAL("3D correction incompatible with any other PU correction. Please turn off any PU residual options.");
+    return StatusCode::FAILURE;
+
+  }
+  
   m_jetStartScale = m_config->GetValue("PileupStartingScale","JetConstitScaleMomentum");
   ATH_MSG_INFO("JetPileupCorrection: Starting scale: " << m_jetStartScale);
   if ( m_jetStartScale == "DO_NOT_USE" ) {
@@ -102,7 +111,11 @@ StatusCode JetPileupCorrection::initializeTool(const std::string& name) {
   if ( m_useFull4vectorArea ) ATH_MSG_INFO("  Full 4-vector jet area correction is activated."); 
   //ATH_MSG_INFO(" \n");
 
-  if ( m_doResidual ) { 
+  if(m_do3Dcorrection){
+    m_residual3DCorr = new PUCorrection::PU3DCorrectionHelper();
+    m_residual3DCorr->loadParameters(m_config->GetValue("PU3DCorrection.constants", "pu3DResidualsConstants.root") );
+    
+  }else if ( m_doResidual ) { 
     std::string suffix = "_Residual";
     m_residualOffsetCorr = new ResidualOffsetCorrection(name+suffix,m_config,m_jetAlgo,m_calibAreaTag,m_isData,m_dev);
     m_residualOffsetCorr->msg().setLevel( this->msg().level() );
@@ -144,7 +157,17 @@ StatusCode JetPileupCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetE
   const double rho = jetEventInfo.rho();
   ATH_MSG_VERBOSE("    Rho = " << rho);
 
-  if(m_useFull4vectorArea) {
+  ATH_MSG_INFO(" AAAAAAAA "<< m_do3Dcorrection );
+  
+  if(m_do3Dcorrection){
+    int NPV = jetEventInfo.NPV();
+    int mu  = jetEventInfo.mu();
+    
+    double pt_calib= m_residual3DCorr->correctedPt(pT_det,  eta_det, jetareaP4.Pt(), rho, mu, NPV ) ;
+    xAOD::JetFourMom_t calibP4 = jetStartP4 * pt_calib/pT_det;
+    jet.setJetP4( calibP4 );
+    
+  } else if(m_useFull4vectorArea) {
     ATH_MSG_VERBOSE("  Applying area-subtraction calibration to jet " << jet.index() << " with pT = " << 0.001*jet.pt() << " GeV");
     //subtract rho * the jet area from the jet
     xAOD::JetFourMom_t calibP4 = jetStartP4 - rho*jetareaP4;
@@ -250,9 +273,8 @@ StatusCode JetPileupCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetE
 	}
 	calibP4 = jetStartP4*area_SF;
       }
-    } else
-      calibP4 = jetStartP4*area_SF;
-
+    } else calibP4 = jetStartP4*area_SF;
+    
     //Attribute to track if a jet has received the origin correction
     jet.setAttribute<int>("OriginCorrected",m_doOrigin);
     //Attribute to track if a jet has received the pileup subtraction (always true if this code was run)
diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
new file mode 100644
index 000000000000..981fb1bcf0ac
--- /dev/null
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -0,0 +1,237 @@
+// this file is -*- C++ -*-
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef JETCALIBTOOLS_PURESIDUAL3DCORRECTION_H
+#define JETCALIBTOOLS_PURESIDUAL3DCORRECTION_H
+////////////////////////////////////////////////////////////////
+/// \class PU3DCorrectionHelper
+///
+/// PU correction helper for "3D" method
+/// This self class allows :
+///  - to load the 3D calibration constants from a filename
+///  - to calculate the corrected pT (or momentum scale factor) given pt, eta, rho, mu,NPV
+/// 
+/// This class is self contained in a single header file : this is to ease its re-usage in
+/// the development phases of calib derivation frameworks.
+///
+
+#include "TH2D.h"
+#include "TFile.h"
+#include "TAxis.h"
+#include <stdio.h>
+#include <vector>
+
+namespace PUCorrection {
+  
+  struct PU3DCorrectionHelper {
+
+    ~PU3DCorrectionHelper(){
+      if(m_inputFile) delete m_inputFile; // this should also delete all objects owned by this file (?)
+      if(m_etaBins) delete m_etaBins;
+    }
+
+    /// Main function which returns the corrected pT
+    float correctedPt(float pt, float eta, float area, float rho, int mu, int NPV ) const {
+      float areaCorr = area*rho*m_rhoEnergyScale;
+
+      float  pt_ref = pt ;
+      float calibration3D = correction3D(pt_ref - areaCorr ,
+					 eta,
+					 mu,
+					 NPV);
+      pt_ref =  pt_ref - areaCorr - calibration3D;
+      float deltaPt = deltaPtCorrection( pt_ref, eta );
+
+      return pt -areaCorr - calibration3D + deltaPt;      
+    }
+
+    /// same as above but returns the ration pT_corrected/pT_uncorrected
+    float correctionFactor(float pt, float eta, float area, float rho, int mu, int NPV ) const {
+      float ptCorr = correctedPt(pt,eta,area,rho,mu,NPV);
+      return ptCorr/pt;
+    }
+
+    
+    
+    /// calculate the mu,NPV dependent part of the correction
+    float correction3D(float pt, float eta , int mu, int NPV) const {
+      int muNPVbin = m_ref3DHisto->FindBin(mu, NPV);
+      int etaBin = m_etaBins->FindFixBin(fabs(eta)) - 1;
+      float t0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
+      if(t0 <= -999.9) {
+	muNPVbin = m_closestNonEmpty[etaBin][muNPVbin];
+      }
+      //std::cout << " etaBin "<< etaBin << "  muNPVbin "<< muNPVbin  << std::endl;
+      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
+      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
+      
+      if(m_use3Dp2) {
+	float p2= m_3Dp2_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin) ;
+	return p0+p1*pt+p2*log(pt); 
+      }
+      return p0+p1*pt;            
+    }
+
+    
+    /// calculate the mu,NPV dependent part of the correction (this is only used for tests and validation)    
+    float correction3D_noextrap(float pt, float eta , int mu, int NPV) const {
+      int muNPVbin = m_ref3DHisto->FindBin(mu, NPV);
+      int etaBin = m_etaBins->FindFixBin(fabs(eta)) - 1;
+      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
+      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);       
+
+      //      std::cout << " etaBin "<< etaBin << "  muNPVbin "<< muNPVbin  << "   p0="<<p0<<std::endl;
+      
+      if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
+      
+      if(m_use3Dp2) {
+	float p2= m_3Dp2_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin) ;
+	if ( p2<=-999.9 ) return 0;
+	return p0+p1*log(pt-p2);	
+      }
+      return p0+p1*pt;            
+    }
+
+
+    // Just a test to see if we get smoother calib by doing an interpolation at point (mu,NPV), not used yet
+    float correction3D_interp(float pt, float eta , int mu, int NPV) const {
+      int etaBin = m_etaBins->FindFixBin(fabs(eta)) - 1;
+      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->Interpolate(mu, NPV);
+      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->Interpolate(mu,NPV);
+      
+      if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
+      
+      if(m_use3Dp2) {
+	float p2= m_3Dp2_vs_muNPV[ etaBin ]->Interpolate(mu,NPV) ;
+	if ( p2<=-999.9 ) return 0;
+	return p0+p1*log(pt-p2);	
+      }
+      return p0+p1*pt;            
+    }
+
+    
+    float deltaPtCorrection(float pt, float eta) const {
+      int etabin = m_Dptp0_vs_eta->FindBin(fabs(eta)) ;
+      //std::cout << " delta etabin "<< etabin<< std::endl;
+      float p0 = m_Dptp0_vs_eta->GetBinContent(etabin);
+      float p1 = m_Dptp1_vs_eta->GetBinContent(etabin);
+      return p0+pt*p1;
+    }
+    
+
+    
+    /// Loads the calib constants from histograms in TFile named fileName. 
+    void loadParameters(const std::string & fileName,
+			const std::string &param3D_name = "param3D",
+			const std::string &paramDelta_name = "paramDeltaPt",
+			const std::string &etaBins_name = "etaBins"
+			){
+      m_inputFile = new TFile(fileName.c_str() );
+      std::vector<float> * etaBins_v = (std::vector<float>*)m_inputFile->Get(etaBins_name.c_str());
+      std::vector<double> tmp(etaBins_v->begin(), etaBins_v->end() );
+      m_etaBins = new TAxis( tmp.size()-1, tmp.data() );
+      TList *param3D_l = (TList*) m_inputFile->Get(param3D_name.c_str());
+
+      TList *param3D_p0 = (TList*) param3D_l->At(0);
+      m_3Dp0_vs_muNPV.resize( param3D_p0->GetSize() );
+      TList *param3D_p1 = (TList*) param3D_l->At(1);
+      m_3Dp1_vs_muNPV.resize( param3D_p1->GetSize() );
+
+      TList *param3D_p2 = nullptr;
+      if(m_use3Dp2) {
+	param3D_p2 = (TList*) param3D_l->At(2);
+	m_3Dp2_vs_muNPV.resize( param3D_p2->GetSize() );
+      }
+
+      for(size_t i=0 ; i<(etaBins_v->size()-1); i++){
+	m_3Dp0_vs_muNPV[i] = (TH2D*)param3D_p0->At(i);
+	m_3Dp1_vs_muNPV[i] = (TH2D*)param3D_p1->At(i);
+	if(m_use3Dp2) m_3Dp2_vs_muNPV[i] = (TH2D*)param3D_p2->At(i);
+      }
+      m_ref3DHisto = m_3Dp0_vs_muNPV[0];
+      
+      TList* paramDelta_l = (TList*) m_inputFile->Get(paramDelta_name.c_str());
+      m_Dptp0_vs_eta = (TH1F*) paramDelta_l->At(0);
+      m_Dptp1_vs_eta = (TH1F*) paramDelta_l->At(1);      
+
+      setupClosestNonEmptyBins();
+    }
+
+
+
+    void setupClosestNonEmptyBins(){
+      // Pre calculate the positions of the closest non empty bins 
+
+      // we have a map (bin -> non-empty bin) for each eta slice :
+      m_closestNonEmpty.resize( m_3Dp0_vs_muNPV.size() );
+      for(size_t etabin=0;  etabin< m_closestNonEmpty.size() ;etabin++ ){
+
+	TH2D *refHisto =  m_3Dp0_vs_muNPV[etabin] ;
+	int nTot = refHisto->GetNcells();
+	TAxis * xax = refHisto->GetXaxis();
+	TAxis * yax = refHisto->GetYaxis();
+	float xscale = 1./(xax->GetXmax()-xax->GetXmin()); xscale *= xscale;
+	float yscale = 1./(yax->GetXmax()-yax->GetXmin()); yscale *= yscale;
+	int nbinX = xax->GetNbins();
+	int nbinY = yax->GetNbins();
+	// initialize the map with -1 :
+	m_closestNonEmpty[etabin].resize( nTot, -1 );
+
+	// loop over each bin
+	for(int xi=1;xi<nbinX+1;xi++) for(int yi=1;yi<nbinY+1;yi++) {
+	    int bi = refHisto->GetBin(xi,yi);
+	    if(refHisto->GetBinContent( bi ) >-999.) continue; // non empty bin, we don't need to find anything for it.
+
+	    int clBin = -1;
+	    float x0 = xax->GetBinCenter(xi);
+	    float y0 = yax->GetBinCenter(yi);
+	    float minDr2=1e10;
+	    // loop over other bins to find the closest non-empty :
+	    for(int xj=1;xj<nbinX+1;xj++) for(int yj=1;yj<nbinY+1;yj++) {
+		int bj = refHisto->GetBin(xj,yj);
+		if(refHisto->GetBinContent( bj ) <=-999.) continue; // this is an empty bin
+		float x = xax->GetBinCenter(xj);
+		float y = yax->GetBinCenter(yj);
+		float dr2 = (x0-x)*(x0-x)*xscale+(y0-y)*(y0-y)*yscale;
+		if(dr2<minDr2){ minDr2 = dr2; clBin = bj;}
+	      }
+	    m_closestNonEmpty[etabin][bi] = clBin;
+	  }
+
+      }
+    }
+
+
+    
+    //************************************
+    // data members
+    
+    // 3D corrections parameters :
+    TAxis *m_etaBins = nullptr;
+    std::vector<TH2D*> m_3Dp0_vs_muNPV;
+    std::vector<TH2D*> m_3Dp1_vs_muNPV;
+    std::vector<TH2D*> m_3Dp2_vs_muNPV;
+    TH2D* m_ref3DHisto = nullptr;
+    bool m_use3Dp2=true;
+
+    // DeltaPt corrections parameters :
+    TH1F* m_Dptp0_vs_eta=nullptr;
+    TH1F* m_Dptp1_vs_eta=nullptr;
+
+    // keep a pointer to the file from which we read constants from.
+    TFile* m_inputFile = nullptr;
+
+
+    //
+    float m_maxPt=170.0 ; // GeV !!
+    float m_rhoEnergyScale = 0.001; // 0.001 when rho is given in MeV. 
+
+    // ***************
+    // 
+    std::vector< std::vector<int> > m_closestNonEmpty;
+  };
+
+}
+
+#endif
-- 
GitLab


From 2c5b2197477652b3caab30c46bee45e5086276ce Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Tue, 7 Apr 2020 14:52:56 +0200
Subject: [PATCH 2/9] update calib for PFlow jets including PU 3D residuals

---
 Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py b/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py
index d2171455e9d1..88a07b444449 100644
--- a/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py
+++ b/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py
@@ -52,12 +52,12 @@ class JetRecCalibrationFinder:
     "triggerTrim"     : "JES_MC15recommendation_FatJet_June2015.config",
     # Note that these are not available in the most recent calibration area
     "reco"            : "JES_MC15cRecommendation_May2016_rel21.config",
-    "pflow"           : "JES_MC15cRecommendation_PFlow_Aug2016_rel21.config"
+    "pflow"           : "JES_MC16Recommendation_Consolidated_PFlow_Apr2019_Rel21.config"
   }
 
   # Default the calibration area tag to that used for T0 reconstruction for consistency
   # This is based on the initial R21 production in 2016
-  def find(self, alg, rad, inpin, seq, configkeyin, evsprefix, calibareatag="00-04-77"):
+  def find(self, alg, rad, inpin, seq, configkeyin, evsprefix, calibareatag="00-04-82"):
     from JetCalibTools.JetCalibToolsConf import JetCalibrationTool
     from JetRec.JetRecStandardToolManager import jtm
     inp = inpin
-- 
GitLab


From 71009eb4ea06cc021d0f8d8f6aff1894ebabd82c Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Tue, 7 Apr 2020 17:27:49 +0200
Subject: [PATCH 3/9] fix minor point from MR review

---
 .../Root/JetPileupCorrection.cxx              |  1 -
 .../Root/PUResidual3DCorrection.h             | 22 +++++++++----------
 2 files changed, 10 insertions(+), 13 deletions(-)

diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
index aa49d5defbf5..d49f6099c016 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
@@ -157,7 +157,6 @@ StatusCode JetPileupCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetE
   const double rho = jetEventInfo.rho();
   ATH_MSG_VERBOSE("    Rho = " << rho);
 
-  ATH_MSG_INFO(" AAAAAAAA "<< m_do3Dcorrection );
   
   if(m_do3Dcorrection){
     int NPV = jetEventInfo.NPV();
diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index 981fb1bcf0ac..a49351c5e84d 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -32,7 +32,7 @@ namespace PUCorrection {
     }
 
     /// Main function which returns the corrected pT
-    float correctedPt(float pt, float eta, float area, float rho, int mu, int NPV ) const {
+    float correctedPt(float pt, float eta, float area, float rho, float mu, int NPV ) const {
       float areaCorr = area*rho*m_rhoEnergyScale;
 
       float  pt_ref = pt ;
@@ -47,7 +47,7 @@ namespace PUCorrection {
     }
 
     /// same as above but returns the ration pT_corrected/pT_uncorrected
-    float correctionFactor(float pt, float eta, float area, float rho, int mu, int NPV ) const {
+    float correctionFactor(float pt, float eta, float area, float rho, float mu, int NPV ) const {
       float ptCorr = correctedPt(pt,eta,area,rho,mu,NPV);
       return ptCorr/pt;
     }
@@ -55,14 +55,14 @@ namespace PUCorrection {
     
     
     /// calculate the mu,NPV dependent part of the correction
-    float correction3D(float pt, float eta , int mu, int NPV) const {
+    float correction3D(float pt, float eta , float mu, int NPV) const {
       int muNPVbin = m_ref3DHisto->FindBin(mu, NPV);
-      int etaBin = m_etaBins->FindFixBin(fabs(eta)) - 1;
+      int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1;
       float t0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
       if(t0 <= -999.9) {
 	muNPVbin = m_closestNonEmpty[etaBin][muNPVbin];
       }
-      //std::cout << " etaBin "<< etaBin << "  muNPVbin "<< muNPVbin  << std::endl;
+
       float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
       float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
       
@@ -75,13 +75,12 @@ namespace PUCorrection {
 
     
     /// calculate the mu,NPV dependent part of the correction (this is only used for tests and validation)    
-    float correction3D_noextrap(float pt, float eta , int mu, int NPV) const {
+    float correction3D_noextrap(float pt, float eta , float mu, int NPV) const {
       int muNPVbin = m_ref3DHisto->FindBin(mu, NPV);
-      int etaBin = m_etaBins->FindFixBin(fabs(eta)) - 1;
+      int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1;
       float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
       float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);       
 
-      //      std::cout << " etaBin "<< etaBin << "  muNPVbin "<< muNPVbin  << "   p0="<<p0<<std::endl;
       
       if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
       
@@ -95,8 +94,8 @@ namespace PUCorrection {
 
 
     // Just a test to see if we get smoother calib by doing an interpolation at point (mu,NPV), not used yet
-    float correction3D_interp(float pt, float eta , int mu, int NPV) const {
-      int etaBin = m_etaBins->FindFixBin(fabs(eta)) - 1;
+    float correction3D_interp(float pt, float eta , float mu, int NPV) const {
+      int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1;
       float p0 = m_3Dp0_vs_muNPV[ etaBin ]->Interpolate(mu, NPV);
       float p1 = m_3Dp1_vs_muNPV[ etaBin ]->Interpolate(mu,NPV);
       
@@ -112,8 +111,7 @@ namespace PUCorrection {
 
     
     float deltaPtCorrection(float pt, float eta) const {
-      int etabin = m_Dptp0_vs_eta->FindBin(fabs(eta)) ;
-      //std::cout << " delta etabin "<< etabin<< std::endl;
+      int etabin = m_Dptp0_vs_eta->FindBin(std::abs(eta)) ;
       float p0 = m_Dptp0_vs_eta->GetBinContent(etabin);
       float p1 = m_Dptp1_vs_eta->GetBinContent(etabin);
       return p0+pt*p1;
-- 
GitLab


From 07d9b70fa6def3a95aa8a6ad85a959378371cfcd Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Wed, 8 Apr 2020 16:30:46 +0200
Subject: [PATCH 4/9] More protection and refinements for 3D correction

---
 .../Jet/JetCalibTools/Root/JetPileupCorrection.cxx    | 10 ++++++++--
 .../Jet/JetCalibTools/Root/PUResidual3DCorrection.h   | 11 ++++++++---
 .../Jet/JetRec/python/JetRecCalibrationFinder.py      |  4 ++--
 3 files changed, 18 insertions(+), 7 deletions(-)

diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
index d49f6099c016..f491c614c2b0 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
@@ -114,7 +114,12 @@ StatusCode JetPileupCorrection::initializeTool(const std::string& name) {
   if(m_do3Dcorrection){
     m_residual3DCorr = new PUCorrection::PU3DCorrectionHelper();
     m_residual3DCorr->loadParameters(m_config->GetValue("PU3DCorrection.constants", "pu3DResidualsConstants.root") );
-    
+    m_residual3DCorr->m_rhoEnergyScale = m_config->GetValue("PU3DCorrection.rhoEnergyScale", 0.001);
+    m_residual3DCorr->m_pTEnergyScale = m_config->GetValue("PU3DCorrection.pTEnergyScale", 0.001);
+    ATH_MSG_INFO("Pile-up 3D correction. Configured with :");
+    ATH_MSG_INFO("  calib constants file="<< m_config->GetValue("PU3DCorrection.constants", "pu3DResidualsConstants.root") );
+    ATH_MSG_INFO("  rho scale ="<<m_residual3DCorr->m_rhoEnergyScale );
+    ATH_MSG_INFO("  pT scale ="<<m_residual3DCorr->m_pTEnergyScale);    
   }else if ( m_doResidual ) { 
     std::string suffix = "_Residual";
     m_residualOffsetCorr = new ResidualOffsetCorrection(name+suffix,m_config,m_jetAlgo,m_calibAreaTag,m_isData,m_dev);
@@ -163,7 +168,8 @@ StatusCode JetPileupCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetE
     int mu  = jetEventInfo.mu();
     
     double pt_calib= m_residual3DCorr->correctedPt(pT_det,  eta_det, jetareaP4.Pt(), rho, mu, NPV ) ;
-    xAOD::JetFourMom_t calibP4 = jetStartP4 * pt_calib/pT_det;
+    double scaleF = pt_calib < 0 ? 0.01*m_GeV/pT_det : pt_calib/pT_det;
+    xAOD::JetFourMom_t calibP4 = jetStartP4 * scaleF;
     jet.setJetP4( calibP4 );
     
   } else if(m_useFull4vectorArea) {
diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index a49351c5e84d..8c8af60e53ae 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -35,7 +35,7 @@ namespace PUCorrection {
     float correctedPt(float pt, float eta, float area, float rho, float mu, int NPV ) const {
       float areaCorr = area*rho*m_rhoEnergyScale;
 
-      float  pt_ref = pt ;
+      float  pt_ref = pt*m_pTEnergyScale ;
       float calibration3D = correction3D(pt_ref - areaCorr ,
 					 eta,
 					 mu,
@@ -43,7 +43,7 @@ namespace PUCorrection {
       pt_ref =  pt_ref - areaCorr - calibration3D;
       float deltaPt = deltaPtCorrection( pt_ref, eta );
 
-      return pt -areaCorr - calibration3D + deltaPt;      
+      return (pt*m_pTEnergyScale -areaCorr - calibration3D + deltaPt)/m_pTEnergyScale;      
     }
 
     /// same as above but returns the ration pT_corrected/pT_uncorrected
@@ -54,8 +54,10 @@ namespace PUCorrection {
 
     
     
-    /// calculate the mu,NPV dependent part of the correction
+    /// calculate the mu,NPV dependent part of the correction.
+    /// IMPORTANT : the pt must be given in GeV 
     float correction3D(float pt, float eta , float mu, int NPV) const {
+      pt = pt < 1 ? 1 : pt;
       int muNPVbin = m_ref3DHisto->FindBin(mu, NPV);
       int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1;
       float t0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
@@ -75,6 +77,7 @@ namespace PUCorrection {
 
     
     /// calculate the mu,NPV dependent part of the correction (this is only used for tests and validation)    
+    /// IMPORTANT : the pt must be given in GeV 
     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;
@@ -110,6 +113,7 @@ namespace PUCorrection {
     }
 
     
+    /// IMPORTANT : the pt must be given in GeV 
     float deltaPtCorrection(float pt, float eta) const {
       int etabin = m_Dptp0_vs_eta->FindBin(std::abs(eta)) ;
       float p0 = m_Dptp0_vs_eta->GetBinContent(etabin);
@@ -224,6 +228,7 @@ namespace PUCorrection {
     //
     float m_maxPt=170.0 ; // GeV !!
     float m_rhoEnergyScale = 0.001; // 0.001 when rho is given in MeV. 
+    float m_pTEnergyScale = 0.001; // 0.001 when pT is given in MeV. 
 
     // ***************
     // 
diff --git a/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py b/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py
index 88a07b444449..d2171455e9d1 100644
--- a/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py
+++ b/Reconstruction/Jet/JetRec/python/JetRecCalibrationFinder.py
@@ -52,12 +52,12 @@ class JetRecCalibrationFinder:
     "triggerTrim"     : "JES_MC15recommendation_FatJet_June2015.config",
     # Note that these are not available in the most recent calibration area
     "reco"            : "JES_MC15cRecommendation_May2016_rel21.config",
-    "pflow"           : "JES_MC16Recommendation_Consolidated_PFlow_Apr2019_Rel21.config"
+    "pflow"           : "JES_MC15cRecommendation_PFlow_Aug2016_rel21.config"
   }
 
   # Default the calibration area tag to that used for T0 reconstruction for consistency
   # This is based on the initial R21 production in 2016
-  def find(self, alg, rad, inpin, seq, configkeyin, evsprefix, calibareatag="00-04-82"):
+  def find(self, alg, rad, inpin, seq, configkeyin, evsprefix, calibareatag="00-04-77"):
     from JetCalibTools.JetCalibToolsConf import JetCalibrationTool
     from JetRec.JetRecStandardToolManager import jtm
     inp = inpin
-- 
GitLab


From cba70b8d514057f472d76d22babf2b82cd2fe3c4 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Thu, 9 Apr 2020 17:07:14 +0200
Subject: [PATCH 5/9] move test functions and add comments

---
 .../Root/PUResidual3DCorrection.h             | 90 +++++++++++--------
 1 file changed, 54 insertions(+), 36 deletions(-)

diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index 8c8af60e53ae..485495bd07b2 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -76,42 +76,6 @@ namespace PUCorrection {
     }
 
     
-    /// calculate the mu,NPV dependent part of the correction (this is only used for tests and validation)    
-    /// IMPORTANT : the pt must be given in GeV 
-    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;
-      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
-      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);       
-
-      
-      if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
-      
-      if(m_use3Dp2) {
-	float p2= m_3Dp2_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin) ;
-	if ( p2<=-999.9 ) return 0;
-	return p0+p1*log(pt-p2);	
-      }
-      return p0+p1*pt;            
-    }
-
-
-    // Just a test to see if we get smoother calib by doing an interpolation at point (mu,NPV), not used yet
-    float correction3D_interp(float pt, float eta , float mu, int NPV) const {
-      int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1;
-      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->Interpolate(mu, NPV);
-      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->Interpolate(mu,NPV);
-      
-      if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
-      
-      if(m_use3Dp2) {
-	float p2= m_3Dp2_vs_muNPV[ etaBin ]->Interpolate(mu,NPV) ;
-	if ( p2<=-999.9 ) return 0;
-	return p0+p1*log(pt-p2);	
-      }
-      return p0+p1*pt;            
-    }
-
     
     /// IMPORTANT : the pt must be given in GeV 
     float deltaPtCorrection(float pt, float eta) const {
@@ -233,6 +197,60 @@ namespace PUCorrection {
     // ***************
     // 
     std::vector< std::vector<int> > m_closestNonEmpty;
+
+
+
+
+
+
+
+
+
+
+
+
+    // *******************************************************
+    // function belows are not used in the correction evaluation but have proven useful for tests during developments
+    // of the calibration methods. We keep them here just in case.
+    // 
+
+    /// calculate the mu,NPV dependent part of the correction (this is only used for tests and validation)    
+    /// IMPORTANT : the pt must be given in GeV 
+    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;
+      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);
+      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin);       
+
+      
+      if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
+      
+      if(m_use3Dp2) {
+	float p2= m_3Dp2_vs_muNPV[ etaBin ]->GetBinContent(muNPVbin) ;
+	if ( p2<=-999.9 ) return 0;
+	return p0+p1*log(pt-p2);	
+      }
+      return p0+p1*pt;            
+    }
+
+
+    // Just a test to see if we get smoother calib by doing an interpolation at point (mu,NPV), not used yet
+    float correction3D_interp(float pt, float eta , float mu, int NPV) const {
+      int etaBin = m_etaBins->FindFixBin(std::abs(eta)) - 1;
+      float p0 = m_3Dp0_vs_muNPV[ etaBin ]->Interpolate(mu, NPV);
+      float p1 = m_3Dp1_vs_muNPV[ etaBin ]->Interpolate(mu,NPV);
+      
+      if ( (p0<= -999.9) || (p1<=-999.9) ) return 0;
+      
+      if(m_use3Dp2) {
+	float p2= m_3Dp2_vs_muNPV[ etaBin ]->Interpolate(mu,NPV) ;
+	if ( p2<=-999.9 ) return 0;
+	return p0+p1*log(pt-p2);	
+      }
+      return p0+p1*pt;            
+    }
+    
+
   };
 
 }
-- 
GitLab


From 2b402b24b9e22053e44fbb4adae91b40b7a0a4a1 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Tue, 14 Apr 2020 12:06:51 +0200
Subject: [PATCH 6/9] minor fix following MR reviews

---
 Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx | 2 +-
 .../Jet/JetCalibTools/Root/PUResidual3DCorrection.h           | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
index f491c614c2b0..a1805e14ad87 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
@@ -165,7 +165,7 @@ StatusCode JetPileupCorrection::calibrateImpl(xAOD::Jet& jet, JetEventInfo& jetE
   
   if(m_do3Dcorrection){
     int NPV = jetEventInfo.NPV();
-    int mu  = jetEventInfo.mu();
+    float mu  = jetEventInfo.mu();
     
     double pt_calib= m_residual3DCorr->correctedPt(pT_det,  eta_det, jetareaP4.Pt(), rho, mu, NPV ) ;
     double scaleF = pt_calib < 0 ? 0.01*m_GeV/pT_det : pt_calib/pT_det;
diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index 485495bd07b2..11fddd49a11a 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -152,7 +152,7 @@ namespace PUCorrection {
 	    int clBin = -1;
 	    float x0 = xax->GetBinCenter(xi);
 	    float y0 = yax->GetBinCenter(yi);
-	    float minDr2=1e10;
+	    float minDr2=1e10; // just pick a bigger value than any distance in the (mu,NPV) plan
 	    // loop over other bins to find the closest non-empty :
 	    for(int xj=1;xj<nbinX+1;xj++) for(int yj=1;yj<nbinY+1;yj++) {
 		int bj = refHisto->GetBin(xj,yj);
@@ -160,7 +160,7 @@ namespace PUCorrection {
 		float x = xax->GetBinCenter(xj);
 		float y = yax->GetBinCenter(yj);
 		float dr2 = (x0-x)*(x0-x)*xscale+(y0-y)*(y0-y)*yscale;
-		if(dr2<minDr2){ minDr2 = dr2; clBin = bj;}
+		if(dr2<minDr2){ minDr2 = dr2; clBin = bj;} // found a closest-bin candidate
 	      }
 	    m_closestNonEmpty[etabin][bi] = clBin;
 	  }
-- 
GitLab


From db17f92879e535cb9811cf7668446585f4e9cc81 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Tue, 14 Apr 2020 23:36:33 +0200
Subject: [PATCH 7/9] replace pointers by unique_ptr in PU 3D corrections

---
 .../CalibrationMethods/JetPileupCorrection.h  |  2 +-
 .../Root/JetPileupCorrection.cxx              |  3 +-
 .../Root/PUResidual3DCorrection.h             | 43 ++++++++-----------
 3 files changed, 21 insertions(+), 27 deletions(-)

diff --git a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h
index a91f3bdade9d..b1cdff2fc91d 100644
--- a/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/JetCalibTools/CalibrationMethods/JetPileupCorrection.h
@@ -61,7 +61,7 @@ class JetPileupCorrection
   bool m_useFull4vectorArea;
   ResidualOffsetCorrection * m_residualOffsetCorr;
 
-  PUCorrection::PU3DCorrectionHelper *m_residual3DCorr=nullptr;
+  std::unique_ptr<PUCorrection::PU3DCorrectionHelper> m_residual3DCorr;
   
   bool m_doOnlyResidual;
   
diff --git a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
index a1805e14ad87..1b9750f01feb 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
+++ b/Reconstruction/Jet/JetCalibTools/Root/JetPileupCorrection.cxx
@@ -38,7 +38,6 @@ JetPileupCorrection::JetPileupCorrection(const std::string& name, TEnv * config,
 JetPileupCorrection::~JetPileupCorrection() {
 
   if(m_residualOffsetCorr) delete m_residualOffsetCorr;
-  if(m_residual3DCorr) delete m_residual3DCorr;
 }
 
 
@@ -112,7 +111,7 @@ StatusCode JetPileupCorrection::initializeTool(const std::string& name) {
   //ATH_MSG_INFO(" \n");
 
   if(m_do3Dcorrection){
-    m_residual3DCorr = new PUCorrection::PU3DCorrectionHelper();
+    m_residual3DCorr.reset( new PUCorrection::PU3DCorrectionHelper() ) ;
     m_residual3DCorr->loadParameters(m_config->GetValue("PU3DCorrection.constants", "pu3DResidualsConstants.root") );
     m_residual3DCorr->m_rhoEnergyScale = m_config->GetValue("PU3DCorrection.rhoEnergyScale", 0.001);
     m_residual3DCorr->m_pTEnergyScale = m_config->GetValue("PU3DCorrection.pTEnergyScale", 0.001);
diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index 11fddd49a11a..5fe08c421bbc 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -27,8 +27,7 @@ namespace PUCorrection {
   struct PU3DCorrectionHelper {
 
     ~PU3DCorrectionHelper(){
-      if(m_inputFile) delete m_inputFile; // this should also delete all objects owned by this file (?)
-      if(m_etaBins) delete m_etaBins;
+      
     }
 
     /// Main function which returns the corrected pT
@@ -93,11 +92,11 @@ namespace PUCorrection {
 			const std::string &paramDelta_name = "paramDeltaPt",
 			const std::string &etaBins_name = "etaBins"
 			){
-      m_inputFile = new TFile(fileName.c_str() );
-      std::vector<float> * etaBins_v = (std::vector<float>*)m_inputFile->Get(etaBins_name.c_str());
+      TFile tmpF(fileName.c_str() );
+      std::vector<float> * etaBins_v = (std::vector<float>*)tmpF.Get(etaBins_name.c_str());
       std::vector<double> tmp(etaBins_v->begin(), etaBins_v->end() );
-      m_etaBins = new TAxis( tmp.size()-1, tmp.data() );
-      TList *param3D_l = (TList*) m_inputFile->Get(param3D_name.c_str());
+      m_etaBins.reset( new TAxis( tmp.size()-1, tmp.data() ) );
+      TList *param3D_l = (TList*) tmpF.Get(param3D_name.c_str());
 
       TList *param3D_p0 = (TList*) param3D_l->At(0);
       m_3Dp0_vs_muNPV.resize( param3D_p0->GetSize() );
@@ -111,16 +110,15 @@ namespace PUCorrection {
       }
 
       for(size_t i=0 ; i<(etaBins_v->size()-1); i++){
-	m_3Dp0_vs_muNPV[i] = (TH2D*)param3D_p0->At(i);
-	m_3Dp1_vs_muNPV[i] = (TH2D*)param3D_p1->At(i);
-	if(m_use3Dp2) m_3Dp2_vs_muNPV[i] = (TH2D*)param3D_p2->At(i);
+	m_3Dp0_vs_muNPV[i].reset((TH2D*)param3D_p0->At(i));
+	m_3Dp1_vs_muNPV[i].reset((TH2D*)param3D_p1->At(i));
+	if(m_use3Dp2) m_3Dp2_vs_muNPV[i].reset( (TH2D*)param3D_p2->At(i) );
       }
-      m_ref3DHisto = m_3Dp0_vs_muNPV[0];
+      m_ref3DHisto = m_3Dp0_vs_muNPV[0].get();
       
-      TList* paramDelta_l = (TList*) m_inputFile->Get(paramDelta_name.c_str());
-      m_Dptp0_vs_eta = (TH1F*) paramDelta_l->At(0);
-      m_Dptp1_vs_eta = (TH1F*) paramDelta_l->At(1);      
-
+      TList* paramDelta_l = (TList*) tmpF.Get(paramDelta_name.c_str());
+      m_Dptp0_vs_eta.reset( (TH1F*) paramDelta_l->At(0) );
+      m_Dptp1_vs_eta.reset( (TH1F*) paramDelta_l->At(1) ) ;      
       setupClosestNonEmptyBins();
     }
 
@@ -133,7 +131,7 @@ namespace PUCorrection {
       m_closestNonEmpty.resize( m_3Dp0_vs_muNPV.size() );
       for(size_t etabin=0;  etabin< m_closestNonEmpty.size() ;etabin++ ){
 
-	TH2D *refHisto =  m_3Dp0_vs_muNPV[etabin] ;
+	TH2D *refHisto =  m_3Dp0_vs_muNPV[etabin].get() ;
 	int nTot = refHisto->GetNcells();
 	TAxis * xax = refHisto->GetXaxis();
 	TAxis * yax = refHisto->GetYaxis();
@@ -174,19 +172,16 @@ namespace PUCorrection {
     // data members
     
     // 3D corrections parameters :
-    TAxis *m_etaBins = nullptr;
-    std::vector<TH2D*> m_3Dp0_vs_muNPV;
-    std::vector<TH2D*> m_3Dp1_vs_muNPV;
-    std::vector<TH2D*> m_3Dp2_vs_muNPV;
+    std::unique_ptr<TAxis> m_etaBins ;
+    std::vector<std::unique_ptr<TH2D> > m_3Dp0_vs_muNPV;
+    std::vector<std::unique_ptr<TH2D> > m_3Dp1_vs_muNPV;
+    std::vector<std::unique_ptr<TH2D> > m_3Dp2_vs_muNPV;
     TH2D* m_ref3DHisto = nullptr;
     bool m_use3Dp2=true;
 
     // DeltaPt corrections parameters :
-    TH1F* m_Dptp0_vs_eta=nullptr;
-    TH1F* m_Dptp1_vs_eta=nullptr;
-
-    // keep a pointer to the file from which we read constants from.
-    TFile* m_inputFile = nullptr;
+    std::unique_ptr<TH1F> m_Dptp0_vs_eta=nullptr;
+    std::unique_ptr<TH1F> m_Dptp1_vs_eta=nullptr;
 
 
     //
-- 
GitLab


From 5aed5430cde809199879fcf0b1650c0193f590ca Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Wed, 15 Apr 2020 00:14:46 +0200
Subject: [PATCH 8/9] Use SetDirectory() on histos for 3D PU corrections

---
 .../Jet/JetCalibTools/Root/PUResidual3DCorrection.h      | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index 5fe08c421bbc..a6610498f171 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -111,14 +111,21 @@ namespace PUCorrection {
 
       for(size_t i=0 ; i<(etaBins_v->size()-1); i++){
 	m_3Dp0_vs_muNPV[i].reset((TH2D*)param3D_p0->At(i));
+	m_3Dp0_vs_muNPV[i]->SetDirectory(nullptr);
 	m_3Dp1_vs_muNPV[i].reset((TH2D*)param3D_p1->At(i));
-	if(m_use3Dp2) m_3Dp2_vs_muNPV[i].reset( (TH2D*)param3D_p2->At(i) );
+	m_3Dp1_vs_muNPV[i]->SetDirectory(nullptr);
+	if(m_use3Dp2) {
+	  m_3Dp2_vs_muNPV[i].reset( (TH2D*)param3D_p2->At(i) );
+	  m_3Dp2_vs_muNPV[i]->SetDirectory(nullptr);
+	}
       }
       m_ref3DHisto = m_3Dp0_vs_muNPV[0].get();
       
       TList* paramDelta_l = (TList*) tmpF.Get(paramDelta_name.c_str());
       m_Dptp0_vs_eta.reset( (TH1F*) paramDelta_l->At(0) );
+      m_Dptp0_vs_eta->SetDirectory(nullptr);
       m_Dptp1_vs_eta.reset( (TH1F*) paramDelta_l->At(1) ) ;      
+      m_Dptp1_vs_eta->SetDirectory(nullptr);
       setupClosestNonEmptyBins();
     }
 
-- 
GitLab


From 77f24d5415925fed571827cacd1c57702bc94db2 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Delsart <delsart@in2p3.fr>
Date: Wed, 15 Apr 2020 20:27:41 +0200
Subject: [PATCH 9/9] use TFile::Open in PU 3d correction

---
 .../Jet/JetCalibTools/Root/PUResidual3DCorrection.h       | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
index a6610498f171..52aa7dc88705 100644
--- a/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
+++ b/Reconstruction/Jet/JetCalibTools/Root/PUResidual3DCorrection.h
@@ -92,11 +92,11 @@ namespace PUCorrection {
 			const std::string &paramDelta_name = "paramDeltaPt",
 			const std::string &etaBins_name = "etaBins"
 			){
-      TFile tmpF(fileName.c_str() );
-      std::vector<float> * etaBins_v = (std::vector<float>*)tmpF.Get(etaBins_name.c_str());
+      std::unique_ptr<TFile> tmpF(TFile::Open( fileName.c_str() ));
+      std::vector<float> * etaBins_v = (std::vector<float>*)tmpF->Get(etaBins_name.c_str());
       std::vector<double> tmp(etaBins_v->begin(), etaBins_v->end() );
       m_etaBins.reset( new TAxis( tmp.size()-1, tmp.data() ) );
-      TList *param3D_l = (TList*) tmpF.Get(param3D_name.c_str());
+      TList *param3D_l = (TList*) tmpF->Get(param3D_name.c_str());
 
       TList *param3D_p0 = (TList*) param3D_l->At(0);
       m_3Dp0_vs_muNPV.resize( param3D_p0->GetSize() );
@@ -121,7 +121,7 @@ namespace PUCorrection {
       }
       m_ref3DHisto = m_3Dp0_vs_muNPV[0].get();
       
-      TList* paramDelta_l = (TList*) tmpF.Get(paramDelta_name.c_str());
+      TList* paramDelta_l = (TList*) tmpF->Get(paramDelta_name.c_str());
       m_Dptp0_vs_eta.reset( (TH1F*) paramDelta_l->At(0) );
       m_Dptp0_vs_eta->SetDirectory(nullptr);
       m_Dptp1_vs_eta.reset( (TH1F*) paramDelta_l->At(1) ) ;      
-- 
GitLab