ATLAS Offline Software
ShowerDepthUtil.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // Local include(s).
7 
8 // Project include(s).
11 
12 // ROOT include(s).
13 #include <TFile.h>
14 
15 // System include(s).
16 #include <cmath>
17 #include <string>
18 #include <stdexcept>
19 
20 namespace {
21 // Note that std::numeric_limits<float>::epsilon()
22 // is just not large enough for what we need. :-(
23 constexpr float epsilon = 1e-6;
24 // Calibration file / histogram name(s).
25 const char* const CONFIG_FILE_NAME =
26  "ElectronIsolationSelection/v1/CaloDeltaRZ.root";
27 const char* const DATA_HISTO_NAME = "hData";
28 const char* const MC_HISTO_NAME = "hMC";
29 
30 } // namespace
31 
32 namespace CP {
33 
34 // Make the messaging functions available.
35 ANA_MSG_SOURCE(ShowerDepthUtilMessaging, "CP::ShowerDepthUtil");
36 using namespace ShowerDepthUtilMessaging;
37 
39  const std::string filename = PathResolverFindCalibFile(CONFIG_FILE_NAME);
40 
41  m_hData = getHistoFromFile(filename.c_str(), DATA_HISTO_NAME);
42  m_hMC = getHistoFromFile(filename.c_str(), MC_HISTO_NAME);
43  if(!m_hData || !m_hMC){
44  throw std::runtime_error("Could not create instance of ShowerDepthUtil");
45  }
46 }
47 
49 float ShowerDepthUtil::getCorrectedShowerDepthEM1(float etas1, float phi,
50  bool isData) const {
51  return getShowerDepthEM1(etas1) - getRZCorrection(etas1, phi, isData);
52 }
53 
55 float ShowerDepthUtil::getCorrectedShowerDepthEM2(float etas2, float phi,
56  bool isData) const {
57  return getShowerDepthEM2(etas2) - getRZCorrection(etas2, phi, isData);
58 }
59 
66  float radius, aetas1 = std::abs(etas1);
67  if (aetas1 < 0.8) {
68  radius = (1558.859292 - 4.990838 * aetas1 - 21.144279 * aetas1 * aetas1);
69  } else if (aetas1 < 1.5) {
70  radius = (1522.775373 + 27.970192 * aetas1 - 21.104108 * aetas1 * aetas1);
71  } else {
72  // endcap so Z
73  radius = 3790.671754;
74  }
75  if (etas1 < 0. and aetas1 > 1.5) {
76  // negative endcap
77  return -radius;
78  }
79  return radius;
80 }
81 
88  float radius, aetas2 = std::abs(etas2);
89  if (aetas2 < 1.425) { // Barrel
90  radius = (1698.990944 - 49.431767 * aetas2 - 24.504976 * aetas2 * aetas2);
91  } else if (aetas2 < 1.5) { // EME2 in tool
92  radius = (8027.574119 - 2717.653528 * aetas2);
93  } else {
94  // endcap so Z
95  radius = (3473.473909 + 453.941515 * aetas2 - 119.101945 * aetas2 * aetas2);
96  }
97  if (etas2 < 0. and aetas2 > 1.5) {
98  // negative endcap
99  return -radius;
100  }
101  return radius;
102 }
103 
105  float phi, bool isData,
106  int sampling) const {
107  std::pair<float, float> RZ = getCorrectedRZ(eta, phi, isData, sampling);
108  return getEtaDirection(zvertex, RZ.first, RZ.second);
109 }
110 
111 std::pair<float, float> ShowerDepthUtil::getRZ(float eta, int sampling) {
112  if ((sampling != 1 && sampling != 2) || (std::abs(eta) > 10)) {
113  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
114  "Invalid sampling, eta: " << sampling << ", " << eta);
115  return std::make_pair(0., 0.);
116  }
117  float depth =
118  (sampling == 1 ? getShowerDepthEM1(eta) : getShowerDepthEM2(eta));
119  if (std::abs(eta) < 1.5)
120  return std::make_pair(depth, depth * std::sinh(eta));
121  return std::make_pair(depth / std::sinh(eta), depth);
122 }
123 
124 std::optional<float> ShowerDepthUtil::getCaloPointingEta(float etas1,
125  float etas2, float phi,
126  bool isData) const {
127  std::pair<float, float> RZ1 = getCorrectedRZ(etas1, phi, isData, 1);
128  std::pair<float, float> RZ2 = getCorrectedRZ(etas2, phi, isData, 2);
129  // Sanity check
130  if (std::abs(RZ2.first - RZ1.first) < epsilon) {
131  return std::nullopt;
132  }
133 
134  return {std::asinh((RZ2.second - RZ1.second) / (RZ2.first - RZ1.first))};
135 }
136 
137 std::pair<float, float> ShowerDepthUtil::getCorrectedRZ(float eta, float phi,
138  bool isData,
139  int sampling) const {
140  if ((sampling != 1 && sampling != 2) || (std::abs(eta) > 10)) {
141  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
142  "Invalid sampling, eta: " << sampling << ", " << eta);
143  return std::make_pair(0., 0.);
144  }
145  float depth = (sampling == 1 ? getCorrectedShowerDepthEM1(eta, phi, isData)
146  : getCorrectedShowerDepthEM2(eta, phi, isData));
147  if (std::abs(eta) < 1.5)
148  return std::make_pair(depth, depth * std::sinh(eta));
149  return std::make_pair(depth / std::sinh(eta), depth);
150 }
151 
153 float ShowerDepthUtil::getRZCorrection(float eta, float phi,
154  bool isData) const {
155  // Get the correct histogram.
156  const TH1* histo = (isData ? m_hData.get() : m_hMC.get());
157  if (!histo) {
158  return 0;
159  }
160  // Make sure that we can perform the interpolation in both eta and phi.
161  const Int_t etaBin = histo->GetXaxis()->FindFixBin(eta);
162  if (etaBin < 1) {
163  const float etaOld = eta;
164  eta = histo->GetXaxis()->GetBinLowEdge(1) + epsilon;
165  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
166  "Using eta " << eta << " instead of " << etaOld);
167  } else if (etaBin > histo->GetNbinsX()) {
168  const float etaOld = eta;
169  eta = histo->GetXaxis()->GetBinUpEdge(histo->GetNbinsX()) - epsilon;
170  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
171  "Using eta " << eta << " instead of " << etaOld);
172  }
173  const Int_t phiBin = histo->GetYaxis()->FindFixBin(phi);
174  if (phiBin < 1) {
175  const float phiOld = phi;
176  phi = histo->GetYaxis()->GetBinLowEdge(1) + epsilon;
177  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
178  "Using phi " << phi << " instead of " << phiOld);
179  } else if (phiBin > histo->GetNbinsY()) {
180  const float phiOld = phi;
181  phi = histo->GetYaxis()->GetBinUpEdge(histo->GetNbinsY()) - epsilon;
182  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
183  "Using phi " << phi << " instead of " << phiOld);
184  }
185  // Get the correction as an interpolation.
186  return histo->Interpolate(eta, phi);
187 }
188 
189 float ShowerDepthUtil::getEtaDirection(float zvertex, float R, float z) {
190  return std::asinh((z - zvertex) / R);
191 }
192 
193 std::unique_ptr<TH1> ShowerDepthUtil::getHistoFromFile(const char* fileName,
194  const char* histoName) {
195  std::unique_ptr<TFile> f(TFile::Open(fileName, "READ"));
196  if (!f) {
197  ANA_MSG_LVL_SERIOUS(MSG::WARNING,
198  "Could not open file: \"" << fileName << "\"");
199  return {};
200  }
201  TH1* h = dynamic_cast<TH1*>(f->Get(histoName));
202  if (!h) {
203  ANA_MSG_LVL_SERIOUS(MSG::WARNING, "Could not get histogram: \""
204  << histoName << "\" from file: \""
205  << fileName << "\"");
206  return {};
207  }
208  // The file will be deleted. We take ownership.
209  h->SetDirectory(nullptr);
210  return std::unique_ptr<TH1>(h);
211 }
212 
213 } // namespace CP
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
egammaParameters::depth
@ depth
pointing depth of the shower as calculated in egammaqgcld
Definition: egammaParamDefs.h:276
egammaParameters::zvertex
@ zvertex
pointing z at vertex reconstructed from the cluster
Definition: egammaParamDefs.h:270
ShowerDepthUtil.h
CP
Select isolated Photons, Electrons and Muons.
Definition: Control/xAODRootAccess/xAODRootAccess/TEvent.h:49
CP::ShowerDepthUtil::getCaloPointingEta
std::optional< float > getCaloPointingEta(float etas1, float etas2, float phi, bool isData=true) const
Eta direction from samplings 1 and 2 (pointing)
Definition: ShowerDepthUtil.cxx:124
xAOD::phi
setEt phi
Definition: TrigEMCluster_v1.cxx:29
xAOD::etaBin
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap etaBin
Definition: L2StandAloneMuon_v1.cxx:149
CP::ShowerDepthUtil::getShowerDepthEM2
static float getShowerDepthEM2(float etas2)
Shower depth (in mm) vs.
Definition: ShowerDepthUtil.cxx:87
z
#define z
h
CP::ShowerDepthUtil::getHistoFromFile
static std::unique_ptr< TH1 > getHistoFromFile(const char *fileName, const char *histoName)
Return TH1* from file given fileName, histoName.
Definition: ShowerDepthUtil.cxx:193
MessageCheck.h
macros for messaging and checking status codes
hist_file_dump.f
f
Definition: hist_file_dump.py:140
python.Dumpers.asinh
def asinh(x)
helper methods ---------------------------------------------------------—
Definition: Dumpers.py:88
CP::ShowerDepthUtil::getCorrectedEtaDirection
float getCorrectedEtaDirection(float zvertex, float eta, float phi, bool isData=true, int sampling=1) const
Eta direction from zvertex to the shower in the given sampling.
Definition: ShowerDepthUtil.cxx:104
CP::ShowerDepthUtil::getCorrectedShowerDepthEM1
float getCorrectedShowerDepthEM1(float etas1, float phi, bool isData=true) const
Shower depth (in mm) on EM1 vs.
Definition: ShowerDepthUtil.cxx:49
CP::ShowerDepthUtil::getRZ
static std::pair< float, float > getRZ(float eta, int sampling)
Shower depth in R,Z for the given sampling.
Definition: ShowerDepthUtil.cxx:111
CP::ShowerDepthUtil::getShowerDepthEM1
static float getShowerDepthEM1(float etas1)
Shower depth (in mm) vs.
Definition: ShowerDepthUtil.cxx:65
PathResolver.h
xAOD::phiBin
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setPhiMap phiBin
Definition: L2StandAloneMuon_v2.cxx:145
CP::ShowerDepthUtil::getEtaDirection
static float getEtaDirection(float zvertex, float R, float z)
Definition: ShowerDepthUtil.cxx:189
CP::ShowerDepthUtil::getCorrectedRZ
std::pair< float, float > getCorrectedRZ(float eta, float phi, bool isData=true, int sampling=1) const
Return the shower depth in R,Z considering misalignments.
Definition: ShowerDepthUtil.cxx:137
ANA_MSG_LVL_SERIOUS
#define ANA_MSG_LVL_SERIOUS(lvl, xmsg)
Macro used to print "serious" messages.
Definition: Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:270
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
CP::ShowerDepthUtil::getCorrectedShowerDepthEM2
float getCorrectedShowerDepthEM2(float etas2, float phi, bool isData=true) const
Shower depth (in mm) on EM2 vs.
Definition: ShowerDepthUtil.cxx:55
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:283
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:23
CP::ShowerDepthUtil::ShowerDepthUtil
ShowerDepthUtil()
Definition: ShowerDepthUtil.cxx:38
jobOptions.fileName
fileName
Definition: jobOptions.SuperChic_ALP2.py:39
plotBeamSpotCompare.histo
histo
Definition: plotBeamSpotCompare.py:414
CP::ANA_MSG_SOURCE
ANA_MSG_SOURCE(ShowerDepthUtilMessaging, "CP::ShowerDepthUtil")
CP::ShowerDepthUtil::getRZCorrection
float getRZCorrection(float eta, float phi, bool isData=true) const
Return the calorimeter displacement in R(Z) for barrel (endcap)
Definition: ShowerDepthUtil.cxx:153