ATLAS Offline Software
CaloSwCalibHitsCalibration.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /********************************************************************
6 
7 NAME: CaloSwCalibHitsCalibration.cxx
8 PACKAGE: offline/Calorimeter/CaloSwCalibHitsCalibration
9 
10 AUTHORS: L. Carminati
11 CREATED: January 4, 2007
12 
13 PURPOSE: Calibrate SW em clusters with coefficients extracted
14  from Calibration Hits.
15  Calibration constants derived by D. Banfi, L.C., L.
16  Mandelli and E. Nebot Del Busto
17  base class: CaloClusterCorrectionCommon
18  Uses ToolWithConstants to get correction constants
19 
20 Updated: January 17, 2008 (LC)
21  Latest parameterizations from E. Nebot Del Busto for the
22  endcap calibration
23 
24 ********************************************************************/
25 
26 #include "CaloSwCalibHitsCalibration.h"
27 #include "CLHEP/Units/SystemOfUnits.h"
28 
29 
30 using xAOD::CaloCluster;
31 using CLHEP::TeV;
32 
33 
56 void CaloSwCalibHitsCalibration::makeTheCorrection
57 (const Context& myctx,
58  CaloCluster* cluster,
59  const CaloDetDescrElement* /*elt*/,
60  float eta,
61  float adj_eta,
62  float /*phi*/,
63  float /*adj_phi*/,
64  CaloSampling::CaloSample /*samp*/) const
65 {
66  // ??? In principle, we should use adj_eta for the interpolation
67  // and range checks. However, the v2 corrections were derived
68  // using regular eta instead.
69  float the_aeta;
70  if (m_use_raw_eta (myctx))
71  the_aeta = std::abs (adj_eta);
72  else
73  the_aeta = std::abs (eta);
74 
75  const float etamax = m_etamax (myctx);
76  if (the_aeta >= etamax) return;
77 
78  const float eta_start_crack = m_eta_start_crack (myctx);
79  const float eta_end_crack = m_eta_end_crack (myctx);
80  int si;
81  if (the_aeta < eta_start_crack)
82  si = 0;
83  else if (the_aeta > eta_end_crack)
84  si = 1;
85  else {
86  // No corrections are applied for the crack region.
87  return;
88  }
89 
90  const CxxUtils::Array<3> correction = m_correction (myctx);
91 
92 // the_aeta = std::abs (cluster->eta() );
93 
94 
95  ATH_MSG_DEBUG( "************************************************************************************************" << endmsg);
96  ATH_MSG_DEBUG( " USING CALIBHITS CALIBRATION " << endmsg);
97  ATH_MSG_DEBUG( " Tool Name " << name() << endmsg);
98  ATH_MSG_DEBUG( "************************************************************************************************" << endmsg);
99 
100  unsigned int shape[] = {2};
101  CaloRec::WritableArrayData<1> interp_barriers (shape);
102  interp_barriers[0] = eta_start_crack;
103  interp_barriers[1] = eta_end_crack;
104 
105  int ibin = (static_cast<int> (the_aeta / etamax * 100)) ;
106 
107  int ibin_frontCorr = ibin;
108  if (m_fix_v6_pathologies (myctx)) {
109  if (the_aeta>2.35) ibin_frontCorr = (static_cast<int> (2.35 / etamax * 100)) ;
110  }
111 
112 // -------------------------------------------------------------
113 // Load calibration coefficients
114 // -------------------------------------------------------------
115 
116  CaloRec::Array<1> acc = correction[0][ibin];
117  CaloRec::Array<1> ooc = correction[1][ibin];
118  CaloRec::Array<1> lleak = correction[2][ibin];
119  CaloRec::Array<1> froffset = correction[3][ibin_frontCorr];
120  CaloRec::Array<1> frslope = correction[4][ibin_frontCorr];
121  CaloRec::Array<1> sec = correction[5][ibin_frontCorr];
122 
123  ATH_MSG_DEBUG( "Check etas -------------------------------------------------------------------"<< endmsg);
124  ATH_MSG_DEBUG( "Eta --> " << the_aeta << " Bin --> " << ibin <<" Cluster eta = " << cluster->eta() << endmsg);
125  ATH_MSG_DEBUG( "ETA = " << std::abs (adj_eta) << " ADJ_ETA = " << std::abs (eta)<< endmsg);
126 
127  ATH_MSG_DEBUG( "Check calibration coefficients -----------------------------------------------"<< endmsg);
128  ATH_MSG_DEBUG( "Accordion : " << acc[0] <<" " << acc[1]
129  << " " << acc[2] << " " << acc[3] << endmsg);
130  ATH_MSG_DEBUG( "OutOfCOne : " << ooc[0] <<" " << ooc[1]
131  << " " << ooc[2] << " " << ooc[3] << endmsg);
132  ATH_MSG_DEBUG( "Leakage : " << lleak[0] <<" " << lleak[1]
133  << " " << lleak[2] << " " << lleak[3] << endmsg);
134  ATH_MSG_DEBUG( "Front offset: " << froffset[0] <<" "
135  <<froffset[1] << " " << froffset[2] << " " << froffset[3] <<endmsg);
136  ATH_MSG_DEBUG( "Front Slope : " << frslope[0] <<" " << frslope[1]
137  << " " << frslope[2] << " " << frslope[3] << endmsg);
138  ATH_MSG_DEBUG( "Second order: " << sec[0] << " " << sec[1]
139  << " " << sec[2] << " " << sec[3] << endmsg);
140 
141  static const CaloSampling::CaloSample samps[2][4] = {
142  { CaloSampling::PreSamplerB,
143  CaloSampling::EMB1,
144  CaloSampling::EMB2,
145  CaloSampling::EMB3 },
146  { CaloSampling::PreSamplerE,
147  CaloSampling::EME1,
148  CaloSampling::EME2,
149  CaloSampling::EME3 }
150  };
151 
152  // Compute longitudinal barycenter: this is a very old and approximate
153  // parametrization that needs to be updated.
154  double shower_lbary = CaloClusterCorr::CaloSwCalibHitsShowerDepth::depth (the_aeta,
155  eta_start_crack,
156  eta_end_crack,
157  m_sampling_depth(myctx),
158  etamax,
159  cluster, msg() );
160  if (shower_lbary == 0) return;
161 
162  if (shower_lbary < 5. || shower_lbary > 25.) {
163  shower_lbary =15.;
164  ATH_MSG_DEBUG( " replace pathological depth by 15 X0" << endmsg);
165  }
166 
167  // Compute total energy in the accordion (eacc_base) and
168  // presampler (eps_base)
169 
170  double eacc_base = 0;
171  for (int sampling=1; sampling<4; sampling++) {
172  eacc_base += cluster->eSample(samps[si][sampling]);
173  ATH_MSG_DEBUG( "Barrel/endcap = " << si << " Sampling = " << sampling << " Energy -->> " << cluster->eSample(samps[si][sampling]) << endmsg);
174  }
175  double eps_base = cluster->eSample (samps[si][0]);
176 
177  ATH_MSG_DEBUG( "E accordion base --->>>> " << eacc_base << " Eps base " << eps_base << endmsg);
178 
179  // Add a protection against large longitudinal barycenter for clusters from
180  // hadrons which may cause over-calibration. A energy dependent upper limit
181  // on the long bary is introduced (from 20 at 0 to 23 at 1 TeV)
182 
183  double depth_max = 20. + (eacc_base+eps_base)*(3./TeV) ;
184  ATH_MSG_DEBUG( "Raw energy ---->> " << (eacc_base+eps_base) << endmsg) ;
185  ATH_MSG_DEBUG( "Bary max for this event ---->> " << depth_max
186  << endmsg) ;
187 
188  if ( shower_lbary > depth_max ) {
189  shower_lbary = 15.;
190  //shower_lbary = depth_max;
191  ATH_MSG_DEBUG( " replace pathological depth by 15 X0 "
192  << endmsg);
193  }
194 
195 
196 
197 // -------------------------------------------------------------
198 // Now calibrate the accordion
199 // -------------------------------------------------------------
200 
201  double acc_corr =
202  acc[1] + acc[2] * shower_lbary + acc[3] * shower_lbary * shower_lbary ;
203 
204  double e_out_perc = 0;
205  if (the_aeta < eta_start_crack) {
206  e_out_perc = ooc[1] + ooc[2] * shower_lbary + ooc[3] * shower_lbary * shower_lbary ;
207  }
208  else if (the_aeta > eta_end_crack) {
209  e_out_perc = ooc[1] + ooc[2] * shower_lbary + ooc[3] / shower_lbary ;
210  }
211 
212  double e_acc_reco=acc_corr*(eacc_base )*(1+(e_out_perc)*0.01);
213 
214 // -------------------------------------------------------------
215 // Now estimate the longitudinal leakage
216 // -------------------------------------------------------------
217 
218  double e_leak_perc = 0;
219 /*
220  if (the_aeta < eta_start_crack) {
221  e_leak_perc = lleak[1] + lleak[2] * shower_lbary + lleak[3] * exp(shower_lbary);
222  } else if (the_aeta > eta_end_crack) {
223  e_leak_perc = lleak[1] * shower_lbary + lleak[2] * exp(shower_lbary);
224  }
225 */
226  e_leak_perc = lleak[3] + lleak[1] * shower_lbary + lleak[2] * exp(shower_lbary);
227 
228  if (e_leak_perc < 0 ) e_leak_perc = 0.;
229  if (e_leak_perc > 100.) e_leak_perc = 100.;
230  double e_leak_reco = e_leak_perc * (e_acc_reco)*0.01;
231 
232 // -------------------------------------------------------------
233 // Now estimate the energy lost in front of the calorimeter
234 // -------------------------------------------------------------
235 
236  // raw_energy = e_acc_reco ;
237 
238  double raw_energy=e_acc_reco*1e-3;
239  double e_front_reco = eps_base;
240 
241  if (raw_energy <= 1) {
242  // Protect against log() going negative in pow();
243  // also protect against sqrt() of a negative number and div-by-zero.
244  }
245  else if (the_aeta < 1.8) {
246 
247  if (the_aeta < eta_start_crack) {
248  double WpsOff = froffset[1] + froffset[2] * raw_energy +
249  froffset[3] * raw_energy * raw_energy;
250  double WpsSlo = frslope[1] * pow(log(raw_energy),
251  static_cast<double>(frslope[2])) +
252  frslope[3] *sqrt( raw_energy );
253  e_front_reco=WpsOff + WpsSlo*(eps_base );
254 
255  ATH_MSG_DEBUG( " raw event " << raw_energy << endmsg);
256  ATH_MSG_DEBUG( " froffset coeff " << froffset[1] << " " << froffset[2] << " " << froffset[3] << endmsg);
257  ATH_MSG_DEBUG( " frslope coeff " << frslope[1] << " " << frslope[2] << " " << frslope[3] << endmsg);
258  ATH_MSG_DEBUG( " WpsOff,WpsSlo " << WpsOff << " " << WpsSlo << endmsg);
259  ATH_MSG_DEBUG( " eps_base, efront_reco " << eps_base << " " << e_front_reco << endmsg);
260  }
261  else{
262 
263  if (raw_energy<20.) raw_energy=20.;
264 
265  double WpsOff = froffset[1] + froffset[2] * raw_energy +
266  froffset[3] * sqrt(raw_energy);
267  double WpsSlo = frslope[1] + frslope[2] * log(raw_energy) +
268  frslope[3] * sqrt(raw_energy);
269  double WpsSlo2 = sec[1] + sec[2] * raw_energy -
270  sec[3] / (raw_energy * raw_energy) ;
271  e_front_reco=WpsOff + WpsSlo*(eps_base) + WpsSlo2*(eps_base)*(eps_base);
272  if (e_front_reco<0.) e_front_reco= eps_base;
273 
274  ATH_MSG_DEBUG( " raw energy " << raw_energy << endmsg);
275  ATH_MSG_DEBUG( "p1 " << froffset[1] << " "
276  << froffset[2] << " " << froffset[3] << " " << WpsOff << endmsg);
277  ATH_MSG_DEBUG( "p2 " << frslope[1] << " "
278  << frslope[2] << " " << frslope[3] << " " << WpsSlo << endmsg);
279  ATH_MSG_DEBUG( "p3 " << sec[1] << " "
280  << sec[2] << " " << sec[3] << " " << WpsSlo2 << endmsg);
281  ATH_MSG_DEBUG( " WpsOff, WpsSlo, WpsSlo2 " << WpsOff << " " << WpsSlo << " " << WpsSlo2 << endmsg);
282  ATH_MSG_DEBUG( " eps_base, efront_reco " << eps_base << " " << e_front_reco << endmsg);
283 
284  }
285 
286  }
287  else {
288 
289  double p1 = froffset[1] + froffset[2] * raw_energy +
290  froffset[3]* raw_energy * raw_energy ;
291  double p2 = frslope[1] + frslope[2] * raw_energy +
292  frslope[3] * raw_energy * raw_energy ;
293  double p3 = sec[1] + sec[2] * raw_energy +
294  sec[3] * raw_energy * raw_energy ;
295 
296  ATH_MSG_DEBUG( "p1 " << froffset[1] << " "
297  << froffset[2] << " " << froffset[3] << endmsg);
298  ATH_MSG_DEBUG( "p2 " << frslope[1] << " "
299  << frslope[2] << " " << frslope[3] << endmsg);
300  ATH_MSG_DEBUG( "p3 " << sec[1] << " "
301  << sec[2] << " " << sec[3] << endmsg);
302 
303  e_front_reco= (p1 + p2 * shower_lbary + p3 * shower_lbary * shower_lbary);
304  if (e_front_reco<0.) e_front_reco=eps_base;
305  }
306 
307 // -------------------------------------------------------------
308 // Now compute the total energy and finally update the cluster energies
309 // -------------------------------------------------------------
310 
311  double e_calo_reco =e_front_reco + e_leak_reco + e_acc_reco;
312 
313  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Final reco energy ---------------------- " << e_calo_reco << endmsg);
314  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Front ---------------------- " << e_front_reco << endmsg);
315  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Accordion ------------------ " << e_acc_reco << endmsg);
316  ATH_MSG_DEBUG( "CaloSwCalibrationHits::out of cone ---------------- " << acc_corr*(eacc_base )*(e_out_perc)*0.01 << endmsg);
317  ATH_MSG_DEBUG( "CaloSwCalibrationHits::Leakage -------------------- " << e_leak_reco << endmsg);
318 
319 
320 
321 /*
322  -------------------------------------------------------------
323  Conventions discussed in Dec 2007 larg week to fill energies in the samplings:
324  E = total energy including all corrections
325  E0 = total energy in front of the accordion :
326  presampler energy + energy lost in front + energy lost between PS and strips
327  E1 = energy in the strips with no out of cone corrections
328  E2 = energy in the middle with no out of cone corrections
329  E3 = energy in the back with no out of cone correction + longitudinal leakage
330 
331  Please note that E != E0+E1+E2+E3
332  -------------------------------------------------------------
333 */
334  if (m_updateSamplingEnergies (myctx))
335  {
336  // presampler
337 
338  cluster->setEnergy (samps[si][0], e_front_reco );
339 
340  // strips and middle
341 
342  for (int sampling=1; sampling<=2; sampling++){
343  cluster->setEnergy (samps[si][sampling], cluster->eSample(samps[si][sampling]) * acc_corr );
344  }
345 
346  // back
347 
348  cluster->setEnergy (samps[si][3], cluster->eSample (samps[si][3]) * acc_corr + e_leak_reco);
349  }
350 // total energy
351 
352  double e_temp = 0;
353  for (int nl = 0 ; nl< 4 ; nl++) e_temp += cluster->eSample (samps [si][nl]);
354 
355  ATH_MSG_DEBUG( "---------- Sum of the sampling energy --- >> " << e_temp << " EcaloReco = " << e_calo_reco << endmsg);
356 
357  cluster->setE (e_calo_reco);
358 
359  ATH_MSG_DEBUG( "CaloSwCalibHitsCalibration::Energy after correction --> " << cluster->e() << endmsg);
360 
361 }
362 
CaloCluster::eSample
double eSample(sampling_type sampling) const
Retrieve energy in a given sampling.
Definition: Calorimeter/CaloEvent/CaloEvent/CaloCluster.h:975
constants.EMB1
int EMB1
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:53
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
CaloDetDescrElement
This class groups all DetDescr information related to a CaloCell. Provides a generic interface for al...
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:66
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
python.SystemOfUnits.TeV
int TeV
Definition: SystemOfUnits.py:158
CaloSwCalibHitsCalibration::m_correction
Constant< CxxUtils::Array< 3 > > m_correction
Definition: CaloSwCalibHitsCalibration.h:73
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
xAOD::CaloCluster
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
Definition: Event/xAOD/xAODCaloEvent/xAODCaloEvent/CaloCluster.h:19
CaloSwCalibHitsCalibration::m_eta_end_crack
Constant< float > m_eta_end_crack
Definition: CaloSwCalibHitsCalibration.h:76
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
constants.EMB2
int EMB2
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:54
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
CaloSwCalibHitsCalibration::m_updateSamplingEnergies
Constant< bool > m_updateSamplingEnergies
Definition: CaloSwCalibHitsCalibration.h:80
CxxUtils::Array< 3 >
constants.EME1
int EME1
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:55
CxxUtils::WritableArrayData
Definition: Control/CxxUtils/CxxUtils/Array.h:778
CaloCluster
Principal data class for CaloCell clusters.
Definition: Calorimeter/CaloEvent/CaloEvent/CaloCluster.h:79
AthenaPoolTestRead.acc
acc
Definition: AthenaPoolTestRead.py:16
CaloSwCalibHitsCalibration::m_sampling_depth
Constant< CxxUtils::Array< 2 > > m_sampling_depth
Definition: CaloSwCalibHitsCalibration.h:74
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
CaloCell_ID_FCS::EME3
@ EME3
Definition: FastCaloSim_CaloCell_ID.h:26
CaloClusterCorr::CaloSwCalibHitsShowerDepth::depth
static double depth(const float aeta, const float start_crack, const float end_crack, const CaloRec::Array< 2 > &sampling_depth, const float etamax, const xAOD::CaloCluster *cluster, MsgStream &log)
Calculate the depth of the cluster.
Definition: CaloSwCalibHitsShowerDepth.cxx:36
CaloCluster::setE
virtual void setE(double e)
Set energy.
Definition: Calorimeter/CaloEvent/CaloEvent/CaloCluster.h:767
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
CaloSwCalibHitsCalibration::m_etamax
Constant< float > m_etamax
Definition: CaloSwCalibHitsCalibration.h:77
CaloSwCalibHitsCalibration::m_use_raw_eta
Constant< bool > m_use_raw_eta
Definition: CaloSwCalibHitsCalibration.h:78
CaloCluster::eta
virtual double eta() const
Retrieve eta independent of signal state.
Definition: Calorimeter/CaloEvent/CaloEvent/CaloCluster.h:755
CaloCell_ID_FCS::PreSamplerE
@ PreSamplerE
Definition: FastCaloSim_CaloCell_ID.h:23
CaloCell_ID_FCS::PreSamplerB
@ PreSamplerB
Definition: FastCaloSim_CaloCell_ID.h:19
CaloSwCalibHitsCalibration::makeTheCorrection
virtual void makeTheCorrection(const Context &myctx, xAOD::CaloCluster *cluster, const CaloDetDescrElement *elt, float eta, float adj_eta, float phi, float adj_phi, CaloSampling::CaloSample samp) const override
Virtual function for the correction-specific code.
Definition: CaloSwCalibHitsCalibration.cxx:57
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
CaloSwCalibHitsCalibration::m_fix_v6_pathologies
Constant< bool > m_fix_v6_pathologies
Definition: CaloSwCalibHitsCalibration.h:79
CaloCluster::e
virtual double e() const
Retrieve energy independent of signal state.
Definition: Calorimeter/CaloEvent/CaloEvent/CaloCluster.h:753
CaloUtils::ToolConstantsContext
Context object for retrieving ToolConstant values.
Definition: ToolWithConstants.h:61
CaloSwCalibHitsCalibration::m_eta_start_crack
Constant< float > m_eta_start_crack
Definition: CaloSwCalibHitsCalibration.h:75
CaloSwCalibHitsCalibration.h
CaloCell_ID_FCS::EMB3
@ EMB3
Definition: FastCaloSim_CaloCell_ID.h:22
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
constants.EME2
int EME2
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:56
Generated on Tue Jul 2 2024 21:11:54 for ATLAS Offline Software by   doxygen 1.8.18