Merge branch 'ShowerDepthToolFixes-main-20230821' into 'main'

CP::ShowerDepthTool Updates, main branch (2023.08.21.)

See merge request !65183
Conflicts:
	PhysicsAnalysis/ElectronPhotonID/IsolationCorrections/Root/ShowerDepthTool.cxx

PhysicsAnalysis/ElectronPhotonID/IsolationCorrections/IsolationCorrections/ShowerDepthTool.h

 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef ELECTRONISOLATIONSELECTION_SHOWERDEPTHTOOL_H
 #define ELECTRONISOLATIONSELECTION_SHOWERDEPTHTOOL_H
 
-// ROOT include(s):
-#include <TString.h>
+// System include(s).
+#include <memory>
 #include <optional>
+#include <utility>
 
 // Forward declaration(s):
 class TH1;
@@ -16,52 +17,47 @@ namespace CP {
 
   class ShowerDepthTool{
   public :
-    ShowerDepthTool(); 
-    ~ShowerDepthTool(); 
+    ShowerDepthTool();
+    ~ShowerDepthTool();
     /// Function initialising the tool
-    bool initialize();      
-      
+    bool initialize();
+
     /** Shower depth (in mm) on EM1 vs. eta, considering misalignments **/
-    float getCorrectedShowerDepthEM1(const float& etas1,const float& phi,const bool& isData = true) const;
+    float getCorrectedShowerDepthEM1(float etas1, float phi, bool isData = true) const;
 
     /** Shower depth (in mm) on EM2 vs. eta, considering misalignments **/
-    float getCorrectedShowerDepthEM2(const float& etas2,const float& phi,const bool& isData = true) const;
+    float getCorrectedShowerDepthEM2(float etas2, float phi, bool isData = true) const;
 
     /** Return the shower depth in R,Z considering misalignments **/
-    std::pair<float, float> getCorrectedRZ(const float& eta,const float& phi,const bool& isData = true,const int& sampling = 1) const;
-      
+    std::pair<float, float> getCorrectedRZ(float eta, float phi, bool isData = true, int sampling = 1) const;
+
     /** Return the calorimeter displacement in R(Z) for barrel (endcap) **/
-    float getRZCorrection(const float& eta,const float& phi,const bool& isData = true) const;     
-      
+    float getRZCorrection(float eta, float phi, bool isData = true) const;
+
     /** Eta direction from zvertex to the shower in the given sampling **/
-    float getCorrectedEtaDirection(const float& zvertex,const float& eta,const float& phi,const bool& isData=true,const int& sampling = 1) const;
-      
+    float getCorrectedEtaDirection(float zvertex, float eta, float phi, bool isData=true, int sampling = 1) const;
+
     /** Eta direction from samplings 1 and 2 (pointing) **/
-    std::optional<float> getCaloPointingEta(const float& etas1,const float& etas2,const float& phi,const bool& isData=true) const;
+    std::optional<float> getCaloPointingEta(float etas1, float etas2, float phi, bool isData=true) const;
 
     /** Shower depth (in mm) vs. eta on EM1 **/
-    static float getShowerDepthEM1(const float& etas1) ;
-      
+    static float getShowerDepthEM1(float etas1);
+
     /** Shower depth (in mm) vs. eta on EM2 **/
-    static float getShowerDepthEM2(const float& etas2) ;
+    static float getShowerDepthEM2(float etas2);
 
     /** Shower depth in R,Z for the given sampling **/
-    static std::pair<float, float> getRZ(const float& eta,const int& sampling) ;
+    static std::pair<float, float> getRZ(float eta, int sampling);
 
-    static float getEtaDirection(const float& zvertex,const float& R,const float& z) ;
-      
-      
-  protected:
+    static float getEtaDirection(float zvertex, float R, float z);
 
-    /** Return TH1* from file given fileName, histoName **/
-    static TH1* getHistoFromFile(const TString& fileName,const TString& histoName);      
-      
   private:
-    TH1 *m_hData;
-    TH1 *m_hMC;    
-    std::string m_configFileName;
-    TString m_dataHistoName;
-    TString m_mcHistoName;
+
+    /** Return TH1* from file given fileName, histoName **/
+    static std::unique_ptr<TH1> getHistoFromFile(const char* fileName, const char* histoName);
+
+    std::unique_ptr<TH1> m_hData;
+    std::unique_ptr<TH1> m_hMC;
   };
 
 } // namespace CP

PhysicsAnalysis/ElectronPhotonID/IsolationCorrections/Root/ShowerDepthTool.cxx

@@ -2,69 +2,62 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
+// Local include(s).
 #include "IsolationCorrections/ShowerDepthTool.h"
 
+// Project include(s).
+#include "AsgMessaging/MessageCheck.h"
+#include "PathResolver/PathResolver.h"
+
+// ROOT include(s).
 #include <TFile.h>
 #include <TH1.h>
-#include <TSystem.h>
 
+// System include(s).
 #include <cmath>
 
-#include <cmath>
-#include <cstdlib>
 #include <string>
 
+namespace CP{
 
-#include <map>
-#include <memory>
+  // Calibration file / histogram name(s).
+  static const char* const CONFIG_FILE_NAME = "ElectronIsolationSelection/v1/CaloDeltaRZ.root";
+  static const char* const DATA_HISTO_NAME = "hData";
+  static const char* const MC_HISTO_NAME = "hMC";
 
-#include "PathResolver/PathResolver.h"
+  // Make the messaging functions available.
+  ANA_MSG_SOURCE(ShowerDepthToolMessaging, "CP::ShowerDepthTool");
+  using namespace ShowerDepthToolMessaging;
 
-using std::make_pair;
+  ShowerDepthTool::ShowerDepthTool() {}
 
-namespace CP{
+  ShowerDepthTool::~ShowerDepthTool() {}
 
-  ShowerDepthTool::ShowerDepthTool() : 
-    m_hData(nullptr),
-    m_hMC(nullptr),
-    m_configFileName("ElectronIsolationSelection/v1/CaloDeltaRZ.root"),
-    m_dataHistoName("hData"),
-    m_mcHistoName("hMC")
-  {          
-  }
+  bool ShowerDepthTool::initialize()
+  {
+    const std::string filename = PathResolverFindCalibFile( CONFIG_FILE_NAME );
 
-  ShowerDepthTool::~ShowerDepthTool() {
-    delete m_hData;
-    delete m_hMC;
+    m_hData = getHistoFromFile( filename.c_str(), DATA_HISTO_NAME );
+    m_hMC = getHistoFromFile( filename.c_str(), MC_HISTO_NAME );
+
+    return (m_hData && m_hMC);
   }
-      
-  bool ShowerDepthTool::initialize() 
-  {   
-    std::string filename = PathResolverFindCalibFile( m_configFileName );
-       
-    TString Tfilename( filename );    
-
-    m_hData = getHistoFromFile( Tfilename , m_dataHistoName );
-    m_hMC = getHistoFromFile( Tfilename , m_mcHistoName );   
-    
-    return !(m_hData == nullptr || m_hMC == nullptr);
-  }      
-   
+
   /** Shower depth (in mm) on EM1 vs. eta, considering misalignments **/
-  float ShowerDepthTool::getCorrectedShowerDepthEM1(const float& etas1,const float& phi,const bool& isData) const
+  float ShowerDepthTool::getCorrectedShowerDepthEM1(float etas1, float phi, bool isData) const
   {
-    return getShowerDepthEM1(etas1) - getRZCorrection(etas1, phi, isData); 
+    return getShowerDepthEM1(etas1) - getRZCorrection(etas1, phi, isData);
   }
 
   /** Shower depth (in mm) on EM2 vs. eta, considering misalignments **/
-  float ShowerDepthTool::getCorrectedShowerDepthEM2(const float& etas2,const float& phi,const bool& isData) const
+  float ShowerDepthTool::getCorrectedShowerDepthEM2(float etas2, float phi, bool isData) const
   {
-    return getShowerDepthEM2(etas2) - getRZCorrection(etas2, phi, isData); 
+    return getShowerDepthEM2(etas2) - getRZCorrection(etas2, phi, isData);
   }
 
   /** Return the shower depth on sampling 1 given etas1. From:
   https://svnweb.cern.ch/trac/atlasoff/browser/Calorimeter/CaloDetDescr/trunk/src/CaloDepthTool.cxx#L347 **/
-  float ShowerDepthTool::getShowerDepthEM1(const float& etas1) 
+  float ShowerDepthTool::getShowerDepthEM1(float etas1)
   {
     float radius, aetas1 = std::fabs(etas1);
     if (aetas1 < 0.8)
@@ -80,7 +73,7 @@ namespace CP{
 
   /** Return the shower depth on sampling 2 given etas2. From:
   https://svnweb.cern.ch/trac/atlasoff/browser/Calorimeter/CaloDetDescr/trunk/src/CaloDepthTool.cxx#L347 **/
-  float ShowerDepthTool::getShowerDepthEM2(const float& etas2) 
+  float ShowerDepthTool::getShowerDepthEM2(float etas2)
   {
     float radius, aetas2 = std::fabs(etas2);
     if (aetas2 < 1.425) // Barrel, my definition
@@ -95,21 +88,21 @@ namespace CP{
   }
 
 
-  float ShowerDepthTool::getCorrectedEtaDirection(const float& zvertex, 
-                                                  const float& eta,
-                                                  const float& phi, 
-                                                  const bool& isData,
-                                                  const int& sampling) const 
+  float ShowerDepthTool::getCorrectedEtaDirection(float zvertex,
+                                                  float eta,
+                                                  float phi,
+                                                  bool isData,
+                                                  int sampling) const
   {
     std::pair<float, float> RZ = getCorrectedRZ(eta, phi, isData, sampling);
     return getEtaDirection(zvertex, RZ.first, RZ.second);
   }
 
-  std::pair<float,float> ShowerDepthTool::getRZ(const float& eta,const int& sampling) 
+  std::pair<float,float> ShowerDepthTool::getRZ(float eta, int sampling)
   {
     if ((sampling != 1 && sampling != 2) || (std::fabs(eta)>10))
     {
-//       ATH_MSG_INFO( "Invalid sampling: " << sampling );
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Invalid sampling, eta: " << sampling << ", " << eta);
       return std::make_pair(0., 0.);
     }
     float depth = (sampling == 1 ? getShowerDepthEM1(eta) : getShowerDepthEM2(eta) );
@@ -119,7 +112,7 @@ namespace CP{
   }
 
 
-  std::optional<float> ShowerDepthTool::getCaloPointingEta(const float& etas1,const float& etas2,const float& phi,const bool& isData) const
+  std::optional<float> ShowerDepthTool::getCaloPointingEta(float etas1, float etas2, float phi, bool isData) const
   {
     std::pair<float, float> RZ1 = getCorrectedRZ(etas1, phi, isData, 1);
     std::pair<float, float> RZ2 = getCorrectedRZ(etas2, phi, isData, 2);
@@ -132,15 +125,14 @@ namespace CP{
   }
 
 
-
-  std::pair<float, float> ShowerDepthTool::getCorrectedRZ(const float& eta, 
-                                                          const float& phi,
-                                                          const bool& isData,
-                                                          const int& sampling) const 
+  std::pair<float, float> ShowerDepthTool::getCorrectedRZ(float eta,
+                                                          float phi,
+                                                          bool isData,
+                                                          int sampling) const
   {
     if ((sampling != 1 && sampling != 2) || (std::fabs(eta)>10))
     {
-//       ATH_MSG_INFO( "Invalid sampling: " << sampling );
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Invalid sampling, eta: " << sampling << ", " << eta);
       return std::make_pair(0., 0.);
     }
     float depth = (sampling == 1 ? getCorrectedShowerDepthEM1(eta, phi, isData) :
@@ -152,39 +144,72 @@ namespace CP{
 
 
   /** Return the calorimeter displacement in R(Z) for barrel (endcap) **/
-  float ShowerDepthTool::getRZCorrection(const float& eta,const float& phi,const bool& isData) const
+  float ShowerDepthTool::getRZCorrection(float eta, float phi, bool isData) const
   {
-    TH1* histo = (isData ? m_hData : m_hMC);
-    if (!histo)
+    // Get the correct histogram.
+    const TH1* histo = (isData ? m_hData.get() : m_hMC.get());
+    if (!histo) {
       return 0;
-    constexpr float epsilon=1e-5;
-    if (std::fabs(eta)==2.5)
-      return histo->Interpolate(eta*(1-epsilon), phi);
-    else
-      return histo->Interpolate(eta, phi);
+    }
+    // Make sure that we can perform the interpolation in both eta and phi.
+    // Note that std::numeric_limits<float>::epsilon() is just not large enough
+    // for the following. :-(
+    static constexpr float epsilon = 1e-6f;
+    const Int_t etaBin = histo->GetXaxis()->FindFixBin(eta);
+    if (etaBin < 1) {
+      const float etaOld = eta;
+      eta = histo->GetXaxis()->GetBinLowEdge(1) + epsilon;
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Using eta " << eta << " instead of "
+                          << etaOld);
+    }
+    else if (etaBin > histo->GetNbinsX()) {
+      const float etaOld = eta;
+      eta = histo->GetXaxis()->GetBinUpEdge(histo->GetNbinsX()) - epsilon;
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Using eta " << eta << " instead of "
+                          << etaOld);
+    }
+    const Int_t phiBin = histo->GetYaxis()->FindFixBin(phi);
+    if (phiBin < 1) {
+      const float phiOld = phi;
+      phi = histo->GetYaxis()->GetBinLowEdge(1) + epsilon;
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Using phi " << phi << " instead of "
+                          << phiOld);
+    }
+    else if (phiBin > histo->GetNbinsY()) {
+      const float phiOld = phi;
+      phi = histo->GetYaxis()->GetBinUpEdge(histo->GetNbinsY()) - epsilon;
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Using phi " << phi << " instead of "
+                          << phiOld);
+    }
+    // Get the correction as an interpolation.
+    return histo->Interpolate(eta, phi);
   }
-  
-  float ShowerDepthTool::getEtaDirection(const float& zvertex,const float& R,const float& z) 
+
+
+  float ShowerDepthTool::getEtaDirection(float zvertex, float R, float z)
   {
     return std::asinh( (z- zvertex)/R );
-  }  
+  }
 
 
-  TH1* ShowerDepthTool::getHistoFromFile(const TString& fileName,const TString& histoName)
+  std::unique_ptr<TH1> ShowerDepthTool::getHistoFromFile(const char* fileName, const char* histoName)
   {
-    std::unique_ptr<TFile> f(TFile::Open(fileName));
+    std::unique_ptr<TFile> f(TFile::Open(fileName, "READ"));
     if (!f){
-      return nullptr;
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING,
+                          "Could not open file: \"" << fileName << "\"");
+      return {};
     }
     TH1 *h = dynamic_cast<TH1*>( f->Get(histoName) );
     if (!h){
-      f->Close();
-      return nullptr;
+      ANA_MSG_LVL_SERIOUS(MSG::WARNING,
+                          "Could not get histogram: \"" << histoName
+                          << "\" from file: \"" << fileName << "\"");
+      return {};
     }
     //The file we be deleted so use SetDirectory
     h->SetDirectory(nullptr);
-    f->Close();
-    return h;
+    return std::unique_ptr<TH1>(h);
   }
 
 } // namespace CP

PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/Root/PhotonVertexHelpers.cxx

 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
 */
 
 // Local includes
@@ -15,6 +15,9 @@
 // Asg tools
 #include "IsolationCorrections/ShowerDepthTool.h"
 
+// ROOT include(s).
+#include <TString.h>
+
 namespace xAOD {
 namespace PVHelpers {
 
@@ -51,7 +54,7 @@ getZCommonAndError(const xAOD::EventInfo* eventInfo,
   // Beam position is the base for zCommon weighted average
   double beamPosZ = eventInfo->beamPosZ();
   double beamPosSigmaZ = eventInfo->beamPosSigmaZ();
-  
+
   if( beamPosSigmaZ == 0 )
     beamPosSigmaZ = 10;