ATLAS Offline Software
EgammaHitsCalibration.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // ********************************************************************
6 //
7 // NAME: EgammaHitsCalibration.cxx
8 // PACKAGE: Trigger/TrigTools/EgammaHitsCalibration.cxx
9 //
10 // AUTHOR: Werner Spolidoro Freund (wsfreund@cern.ch)
11 //
12 // REFERENCES: Tool to provide Hits Calibration to HLT (BASED ON OFFLINE)
13 //
14 // Updated: January 14, 2011
15 // Modifications on EgammaHitsCalibration updated for offline rev 30069
16 //
17 // ********************************************************************
18 
19 #include "EgammaHitsCalibration.h"
20 //#include "TrigCaloEvent/TrigEMCluster.h"
21 #include "xAODTrigCalo/TrigEMCluster.h"
22 #include "GaudiKernel/MsgStream.h"
23 // Need Cosh
24 #include <math.h>
25 
26 
27 StatusCode EgammaHitsCalibration::initialize(){
28 
29  CHECK (base_class::initialize());
30 
31  ATH_MSG_DEBUG( "Initialize Tool : " << name() );
32 
33 
34 
35  return StatusCode::SUCCESS;
36 }
37 
38 StatusCode EgammaHitsCalibration::finalize(){
39  ATH_MSG_DEBUG( "Finalize Tool : " << name() );
40  return StatusCode::SUCCESS;
41 }
42 
43 void EgammaHitsCalibration::makeCorrection(xAOD::TrigEMCluster* clus,
44  const void * /*param*/) const{
45 
46  const float etamax = m_etamax();
47  const float eta_start_crack = m_eta_start_crack();
48  const float eta_end_crack = m_eta_end_crack();
49  const CxxUtils::Array<3> correction = m_correction();
50 
51  float the_aeta = (clus->eta());
52 
53  the_aeta = (the_aeta>0)?(the_aeta):(- the_aeta);
54 
55  if (the_aeta >= etamax) return;
56 
57 #ifndef NDEBUG
58  ATH_MSG_DEBUG( "************************************************************************************************" );
59  ATH_MSG_DEBUG( " USING CALIBHITS CALIBRATION " );
60  ATH_MSG_DEBUG( "************************************************************************************************" );
61 #endif
62  int si;
63  if (the_aeta < eta_start_crack)
64  si = 0; //barrel
65  else if (the_aeta > eta_end_crack)
66  si = 1; //endcap
67  else {
68  // No corrections are applied for the crack region.
69  return;
70  }
71 
72 
73  int ibin = (static_cast<int> (the_aeta / etamax * 100));
74  int ibin_frontCorr = ibin;
75  if (m_fix_v6_pathologies()) {
76  if ( the_aeta>2.35)
77  ibin_frontCorr = (static_cast<int> (235. / etamax ));
78  }
79 
80 #ifndef NDEBUG
81  ATH_MSG_DEBUG( "Check etas -------------------------------------------------------------------" );
82  ATH_MSG_DEBUG( "Eta --> " << the_aeta << " Bin --> " << ibin <<" Cluster eta = " << clus->eta() );
83 
84  ATH_MSG_DEBUG( "Check calibration coefficients -----------------------------------------------" );
85  ATH_MSG_DEBUG( "Accordion : " << correction[0][ibin][0] <<" " << correction[0][ibin][1] << " " << correction[0][ibin][2] << " " << correction[0][ibin][3] );
86  ATH_MSG_DEBUG( "OutOfCOne : " << correction[1][ibin][0] <<" " << correction[1][ibin][1] << " " << correction[1][ibin][2] << " " << correction[1][ibin][3] );
87  ATH_MSG_DEBUG( "Leakage : " << correction[2][ibin][0] <<" " << correction[2][ibin][1] << " " << correction[2][ibin][2] << " " << correction[2][ibin][3] );
88  ATH_MSG_DEBUG( "Front offset: " << correction[3][ibin_frontCorr][0] <<" " <<correction[3][ibin_frontCorr][1] << " " << correction[3][ibin_frontCorr][2] << " " << correction[3][ibin_frontCorr][3] );
89  ATH_MSG_DEBUG( "Front Slope : " << correction[4][ibin_frontCorr][0] <<" " << correction[4][ibin_frontCorr][1] << " " << correction[4][ibin_frontCorr][2] << " " << correction[4][ibin_frontCorr][3] );
90  ATH_MSG_DEBUG( "Second order: " << correction[5][ibin_frontCorr][0] << " " << correction[5][ibin_frontCorr][1] << " " << correction[5][ibin_frontCorr][2] << " " << correction[5][ibin_frontCorr][3] );
91 #endif
92 
93  EgammaHitsShowerDepth showerDepth;
94  float shower_lbary = showerDepth.depth(the_aeta,
95  eta_start_crack,
96  eta_end_crack,
97  m_sampling_depth(),
98  etamax,
99  clus);
100  if (!shower_lbary)
101  return;
102 
103  if (shower_lbary < 5. || shower_lbary > 25.){
104  shower_lbary = 15.;
105 #ifndef NDEBUG
106  ATH_MSG_DEBUG( "replace pathological depth by 15 X0" );
107 #endif
108  }
109 
110  // Compute total energy in the Accordion (eacc_base) and
111  // presampler (eps_base)
112 
113  float eacc_base = 0;
114  for (int sampling=1; sampling<4; sampling++) {
115  eacc_base += clus->energy(EgammaHitsShowerDepth::m_samps[si][sampling]);
116 #ifndef NDEBUG
117  ATH_MSG_DEBUG( "Barrel/endcap = " << si << " Sampling = " << sampling << " Energy -->> " << clus->energy(EgammaHitsShowerDepth::m_samps[si][sampling]) );
118 #endif
119  }
120  float eps_base = clus->energy (EgammaHitsShowerDepth::m_samps[si][0]);
121 
122 #ifndef NDEBUG
123  ATH_MSG_DEBUG( "E accordion base --->>>> " << eacc_base << " Eps base " << eps_base );
124 #endif
125 
126 
127  // Add a protection against large longitudinal barycenter for clusters from
128  // hadrons which may cause over-calibration. A energy dependent upper limit
129  // on the long bary is introduced (from 20 at 0 to 23 at 1 TeV)
130 
131  float depth_max = 20. + 3. * (eacc_base+eps_base) * 1e-6; // divided by TeV
132 #ifndef NDEBUG
133  ATH_MSG_DEBUG( "Raw energy ---->> " << (eacc_base+eps_base) );
134  ATH_MSG_DEBUG( "Bary max for this event ---->> " << depth_max );
135 #endif
136 
137  if ( shower_lbary > depth_max ) {
138  shower_lbary = 15.;
139  //shower_lbary = depth_max;
140 #ifndef NDEBUG
141  ATH_MSG_DEBUG( " replace pathological depth by 15 X0 " );
142 #endif
143  }
144  // -------------------------------------------------------------
145  // Now calibrate the Accordion
146  // -------------------------------------------------------------
147 
148  float acc_corr = correction[0][ibin][1] + correction[0][ibin][2] * shower_lbary + correction[0][ibin][3] * shower_lbary * shower_lbary ;
149 
150  float e_out_perc = 0;
151  if (the_aeta < eta_start_crack) {
152  e_out_perc = correction[1][ibin][1] + correction[1][ibin][2] * shower_lbary + correction[1][ibin][3] * shower_lbary * shower_lbary ;
153  }
154  else {
155  e_out_perc = correction[1][ibin][1] + correction[1][ibin][2] * shower_lbary + correction[1][ibin][3] / shower_lbary ;
156  }
157 
158  float e_acc_reco=acc_corr*(eacc_base )*(1+(e_out_perc)/100);
159 
160  // -------------------------------------------------------------
161  // Now estimate the longitudinal leakage
162  // -------------------------------------------------------------
163 
164  float e_leak_perc = 0;
165  e_leak_perc = correction[2][ibin][3] + correction[2][ibin][1] * shower_lbary + correction[2][ibin][2] * exp(shower_lbary);
166 
167  if (e_leak_perc < 0 ) e_leak_perc = 0.;
168  if (e_leak_perc > 100.) e_leak_perc = 100.;
169  float e_leak_reco = e_leak_perc * (e_acc_reco)/100;
170 
171  // -------------------------------------------------------------
172  // Now estimate the energy lost in front of the calorimeter
173  // -------------------------------------------------------------
174 
175  float raw_energy=e_acc_reco*.001;
176  float e_front_reco = eps_base;
177 
178  if (raw_energy <= 1) {
179  // Protect against log() going negative in pow();
180  // also protect against sqrt() of a negative number and div-by-zero.
181  }
182  else if (the_aeta < 1.8) {
183 
184  if (the_aeta < eta_start_crack) {
185  float WpsOff = correction[3][ibin_frontCorr][1] + correction[3][ibin_frontCorr][2] * raw_energy + correction[3][ibin_frontCorr][3] * raw_energy * raw_energy;
186  float WpsSlo = correction[4][ibin_frontCorr][1] * pow(log(raw_energy),correction[4][ibin_frontCorr][2]) + correction[4][ibin_frontCorr][3] *sqrt( raw_energy );
187  e_front_reco=WpsOff + WpsSlo*(eps_base);
188 
189 #ifndef NDEBUG
190  ATH_MSG_DEBUG( " raw event " << raw_energy );
191  ATH_MSG_DEBUG( " froffset coeff " << correction[3][ibin_frontCorr][1] << " " << correction[3][ibin_frontCorr][2] << " " << correction[3][ibin_frontCorr][3] );
192  ATH_MSG_DEBUG( " frslope coeff " << correction[4][ibin_frontCorr][1] << " " << correction[4][ibin_frontCorr][2] << " " << correction[4][ibin_frontCorr][3] );
193  ATH_MSG_DEBUG( " WpsOff,WpsSlo " << WpsOff << " " << WpsSlo );
194  ATH_MSG_DEBUG( " eps_base, efront_reco " << eps_base << " " << e_front_reco );
195 #endif
196  }
197  else{
198 
199 
200  if (raw_energy<20.) raw_energy=20.;
201 
202  float WpsOff = correction[3][ibin_frontCorr][1] + correction[3][ibin_frontCorr][2] * raw_energy + correction[3][ibin_frontCorr][3] * sqrt(raw_energy);
203  float WpsSlo = correction[4][ibin_frontCorr][1] + correction[4][ibin_frontCorr][2] * log(raw_energy) + correction[4][ibin_frontCorr][3] * sqrt(raw_energy);
204  float WpsSlo2 = correction[5][ibin_frontCorr][1] + correction[5][ibin_frontCorr][2] * raw_energy - correction[5][ibin_frontCorr][3] / (raw_energy * raw_energy) ;
205  e_front_reco=WpsOff + WpsSlo*(eps_base) + WpsSlo2*(eps_base)*(eps_base);
206  if (e_front_reco<0.) e_front_reco= eps_base;
207 #ifndef NDEBUG
208  ATH_MSG_DEBUG( " raw energy " << raw_energy );
209  ATH_MSG_DEBUG( "p1 " << correction[3][ibin_frontCorr][1] << " " << correction[3][ibin_frontCorr][2] << " " << correction[3][ibin_frontCorr][3] << " " << WpsOff );
210  ATH_MSG_DEBUG( "p2 " << correction[4][ibin_frontCorr][1] << " " << correction[4][ibin_frontCorr][2] << " " << correction[4][ibin_frontCorr][3] << " " << WpsSlo );
211  ATH_MSG_DEBUG( "p3 " << correction[5][ibin_frontCorr][1] << " " << correction[5][ibin_frontCorr][2] << " " << correction[5][ibin_frontCorr][3] << " " << WpsSlo2 );
212  ATH_MSG_DEBUG( " WpsOff, WpsSlo, WpsSlo2 " << WpsOff << " " << WpsSlo << " " << WpsSlo2 );
213  ATH_MSG_DEBUG( " eps_base, efront_reco " << eps_base << " " << e_front_reco );
214 #endif
215 
216  }
217 
218  }
219  else {
220  float p1 = correction[3][ibin_frontCorr][1] + correction[3][ibin_frontCorr][2] * raw_energy + correction[3][ibin_frontCorr][3] * raw_energy * raw_energy ;
221  float p2 = correction[4][ibin_frontCorr][1] + correction[4][ibin_frontCorr][2] * raw_energy + correction[4][ibin_frontCorr][3] * raw_energy * raw_energy ;
222  float p3 = correction[5][ibin_frontCorr][1] + correction[5][ibin_frontCorr][2] * raw_energy + correction[5][ibin_frontCorr][3] * raw_energy * raw_energy ;
223 
224 #ifndef NDEBUG
225  ATH_MSG_DEBUG( "p1 " << correction[3][ibin_frontCorr][1] << " " << correction[3][ibin_frontCorr][2] << " " << correction[3][ibin_frontCorr][3] );
226  ATH_MSG_DEBUG( "p2 " << correction[4][ibin_frontCorr][1] << " " << correction[4][ibin_frontCorr][2] << " " << correction[4][ibin_frontCorr][3] );
227  ATH_MSG_DEBUG( "p3 " << correction[5][ibin_frontCorr][1] << " " << correction[5][ibin_frontCorr][2] << " " << correction[5][ibin_frontCorr][3] );
228 #endif
229 
230  e_front_reco= (p1 + p2 * shower_lbary + p3 * shower_lbary * shower_lbary);
231  if (e_front_reco<0.) e_front_reco=eps_base;
232  }
233 
234  // -------------------------------------------------------------
235  // Now compute the total energy and finally update the cluster energies
236  // -------------------------------------------------------------
237 
238  float e_calo_reco =e_front_reco + e_leak_reco + e_acc_reco;
239 
240 #ifndef NDEBUG
241  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Final reco energy ---------------------- " << e_calo_reco );
242  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Front ---------------------- " << e_front_reco );
243  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Accordion ------------------ " << e_acc_reco );
244  ATH_MSG_DEBUG( "CaloSwCalibrationHits::out of cone ---------------- " << acc_corr*(eacc_base )*(e_out_perc)/100 );
245  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Leakage -------------------- " << e_leak_reco );
246 #endif
247 
248 
249 
250  /*
251  -------------------------------------------------------------
252  Conventions discussed in Dec 2007 larg week to fill energies in the samplings:
253  E = total energy including all corrections
254  E0 = total energy in front of the accordion :
255  presampler energy + energy lost in front + energy lost between PS and strips
256  E1 = energy in the strips with no out of cone corrections
257  E2 = energy in the middle with no out of cone corrections
258  E3 = energy in the back with no out of cone correction + longitudinal leakage
259 
260  Please note that E != E0+E1+E2+E3
261  -------------------------------------------------------------
262  */
263 
264  // presampler
265  clus->setEnergy (EgammaHitsShowerDepth::m_samps[si][0], e_front_reco );
266 
267  // strips and middle
268  for (int sampling=1; sampling<=2; sampling++){
269  clus->setEnergy (EgammaHitsShowerDepth::m_samps[si][sampling], clus->energy(EgammaHitsShowerDepth::m_samps[si][sampling]) * acc_corr );
270  }
271 
272  // back
273  clus->setEnergy (EgammaHitsShowerDepth::m_samps[si][3], clus->energy (EgammaHitsShowerDepth::m_samps[si][3]) * acc_corr + e_leak_reco);
274 
275 #ifndef NDEBUG
276  // total energy
277  float e_temp = 0;
278  for (int nl = 0 ; nl< 4 ; nl++) e_temp += clus->energy (EgammaHitsShowerDepth::m_samps[si][nl]);
279  ATH_MSG_DEBUG( "---------- Sum of the sampling energy --- >> " << e_temp << " EcaloReco = " << e_calo_reco );
280  ATH_MSG_DEBUG( "CaloSwCalibHitsCalibration Energy after correction --> " << clus->energy() );
281 #endif
282 
283  clus->setEnergy(e_calo_reco);
284 
285  clus->setEt(clus->energy()/cosh(the_aeta));
286 
287 }
288 
xAOD::TrigEMCluster_v1::eta
float eta() const
get Eta (calibrated)
xAOD::TrigEMCluster_v1::setEnergy
void setEnergy(float energy)
set Energy (calibrated)
initialize
void initialize()
Definition: run_EoverP.cxx:894
TrigEMCluster.h
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
EgammaHitsCalibration::initialize
virtual StatusCode initialize() override
Initialization of the tool.
Definition: EgammaHitsCalibration.cxx:27
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
EgammaHitsCalibration::m_etamax
Constant< float > m_etamax
Definition: EgammaHitsCalibration.h:51
ParticleJetTools::p3
Amg::Vector3D p3(const xAOD::TruthVertex *p)
Definition: ParticleJetLabelCommon.cxx:55
tools.zlumi_mc_cf.correction
def correction(mu, runmode, campaign, run=None)
Definition: zlumi_mc_cf.py:4
EgammaHitsShowerDepth
Definition: EgammaHitsShowerDepth.h:15
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
CxxUtils::Array< 3 >
EgammaHitsCalibration::makeCorrection
virtual void makeCorrection(xAOD::TrigEMCluster *, const void *v=nullptr) const override
method to perform the correction.
Definition: EgammaHitsCalibration.cxx:43
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
EgammaHitsCalibration::m_correction
Constant< CxxUtils::Array< 3 > > m_correction
Definition: EgammaHitsCalibration.h:47
EgammaHitsCalibration::m_sampling_depth
Constant< CxxUtils::Array< 2 > > m_sampling_depth
Definition: EgammaHitsCalibration.h:48
EgammaHitsCalibration::finalize
virtual StatusCode finalize() override
Finalization of the tool.
Definition: EgammaHitsCalibration.cxx:38
EgammaHitsCalibration::m_eta_start_crack
Constant< float > m_eta_start_crack
Definition: EgammaHitsCalibration.h:49
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
EgammaHitsShowerDepth::depth
float depth(const float &aeta, const float start_crack, const float end_crack, const CxxUtils::Array< 2 > &sampling_depth, const float etamax, const xAOD::TrigEMCluster *cluster) const
Calculate the depth of the cluster.
Definition: EgammaHitsShowerDepth.cxx:44
EgammaHitsCalibration::m_eta_end_crack
Constant< float > m_eta_end_crack
Definition: EgammaHitsCalibration.h:50
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
EgammaHitsCalibration::m_fix_v6_pathologies
Constant< bool > m_fix_v6_pathologies
Definition: EgammaHitsCalibration.h:54
xAOD::TrigEMCluster_v1::setEt
void setEt(float)
set Et (calibrated)
EgammaHitsCalibration.h
xAOD::TrigEMCluster_v1
Description of a trigger EM cluster.
Definition: TrigEMCluster_v1.h:28
xAOD::TrigEMCluster_v1::energy
float energy() const
get Energy (calibrated)
EgammaHitsShowerDepth::m_samps
static const CaloSampling::CaloSample m_samps[2][4]
Definition: EgammaHitsShowerDepth.h:24
Generated on Fri Jul 5 2024 21:14:48 for ATLAS Offline Software by   doxygen 1.8.18