ATLAS Offline Software
MuidCaloEnergyTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // MuidCaloEnergyTool
7 // AlgTool for estimating the energy loss of a muon and its error.
8 // The estimation can be based on the parameterization of the energy loss,
9 // on the measurement of the calorimeter or a combination of the two (the
10 // hybrid method).
11 //
12 // Please refer to:
13 // K. Nikolopoulos, D. Fassouliotis, C. Kourkoumelis, and A. Poppleton,
14 // "Event-by-Event Estimate of Muon Energy Loss in ATLAS",
15 // IEEE Trans. Nucl. Sci., vol. 54, no. 5, pt. 2, pp. 1792-1796, Oct. 2007.
16 //
18 
19 
20 #include "MuidCaloEnergyTool.h"
21 
22 #include <cmath>
23 
24 #include "AthenaKernel/Units.h"
25 #include "MuidEvent/CaloMeas.h"
32 #include "muonEvent/CaloEnergy.h"
33 
34 namespace Units = Athena::Units;
35 
36 
37 namespace Rec {
38 
39  MuidCaloEnergyTool::MuidCaloEnergyTool(const std::string& type, const std::string& name, const IInterface* parent) :
41  m_cosmics(false),
42  m_energyLossMeasurement(true),
43  m_forceIsolationFailure(false),
44  m_FSRtreatment(true),
45  m_MOPparametrization(true),
46  m_trackIsolation(false),
47  m_emEtCut(2.5 * Units::GeV),
48  m_emF1Cut(0.15),
49  m_emMinEnergy(2. * Units::GeV),
50  m_hecMinEnergy(10. * Units::GeV),
51  m_maxNTracksIso(2),
52  m_minFinalEnergy(0. * Units::GeV),
53  m_minMuonPt(15. * Units::GeV) {
54  declareInterface<IMuidCaloEnergy>(this);
55  declareProperty("Cosmics", m_cosmics);
56  declareProperty("EnergyLossMeasurement", m_energyLossMeasurement);
57  declareProperty("ForceIsolationFailure", m_forceIsolationFailure);
58  declareProperty("FSRtreatment", m_FSRtreatment);
59  declareProperty("MopParametrization", m_MOPparametrization);
60  declareProperty("TrackIsolation", m_trackIsolation);
61  declareProperty("EmEtCut", m_emEtCut);
62  declareProperty("EmF1Cut", m_emF1Cut);
63  declareProperty("EmMinEnergy", m_emMinEnergy);
64  declareProperty("HecMinEnergy", m_hecMinEnergy);
65  declareProperty("MaxNTracksIso", m_maxNTracksIso);
66  declareProperty("MinFinalEnergy", m_minFinalEnergy);
67  declareProperty("MinMuonPt", m_minMuonPt);
68  }
69 
71 
72  //<<<<<< PUBLIC MEMBER FUNCTION DEFINITIONS >>>>>>
73 
76  ATH_MSG_INFO("Initializing MuidCaloEnergyTool - measured calorimeter energy deposition ");
77  } else {
78  ATH_MSG_INFO("Initializing MuidCaloEnergyTool - parametrized calorimeter energy deposition ");
79  }
80 
81  // get the Tools
82  ATH_CHECK(m_caloParamTool.retrieve());
83  ATH_MSG_DEBUG("Retrieved tool " << m_caloParamTool);
84 
86  ATH_CHECK(m_caloMeasTool.retrieve());
87  ATH_MSG_DEBUG("Retrieved tool " << m_caloMeasTool);
88 
89  if (m_trackIsolation) {
90  ATH_CHECK(m_trackIsolationTool.retrieve());
91  ATH_MSG_DEBUG("Retrieved tool " << m_trackIsolationTool);
92  } else {
93  ATH_MSG_WARNING(" Using energy measurement without trackIsolation ");
94  m_trackIsolationTool.disable();
95  }
96  } else {
97  m_trackIsolationTool.disable();
98  m_caloMeasTool.disable();
99  }
100 
101  return StatusCode::SUCCESS;
102  }
103 
105  ATH_MSG_INFO("Finalizing MuidCaloEnergyTool."
106  << " Tracks used mean: " << m_countMean << ", tracks used mop: " << m_countMop
107  << ", tracks used measurement: " << m_countMeasurement);
108 
109  return StatusCode::SUCCESS;
110  }
111  std::unique_ptr<CaloEnergy> MuidCaloEnergyTool::energyLoss(const EventContext& ctx, double trackMomentum, double eta,
112  double phi) const {
113  // debug
114  ATH_MSG_VERBOSE("Muon with : p = " << trackMomentum / Units::GeV << " Phi = " << phi << " Eta = " << eta);
115  std::unique_ptr<CaloEnergy> caloEnergy;
117  // Energy Dep/Iso from calorimeters (projective assumption)
118  std::unique_ptr<CaloMeas> caloMeas{m_caloMeasTool->energyMeasurement(ctx, eta, phi, eta, phi)};
119  if (caloMeas) { caloEnergy = measurement(ctx, trackMomentum, eta, phi, *caloMeas); }
120  }
121 
122  if (!caloEnergy) {
123  if (m_MOPparametrization) {
124  ++m_countMop;
125  caloEnergy.reset(m_caloParamTool->mopParametrizedEnergy(trackMomentum, eta, phi));
126  ATH_MSG_VERBOSE("Selected the Mop Parametrization value! ==> ");
127  } else {
128  ++m_countMean;
129  caloEnergy.reset(m_caloParamTool->meanParametrizedEnergy(trackMomentum, eta, phi));
130  ATH_MSG_VERBOSE("Selected the Mean Parametrization value! ==> ");
131  }
132  }
133 
134  if (msgLvl(MSG::DEBUG)) {
135  std::string eLossType = " no Calo !!";
136  switch (caloEnergy->energyLossType()) {
137  case CaloEnergy::Parametrized: eLossType = "Parametrized"; break;
138  case CaloEnergy::NotIsolated: eLossType = "NotIsolated "; break;
139  case CaloEnergy::MOP: eLossType = "MOP "; break;
140  case CaloEnergy::Tail: eLossType = "Tail "; break;
141  case CaloEnergy::FSRcandidate: eLossType = "FSRcandidate"; break;
142  default: break;
143  };
144  ATH_MSG_DEBUG(std::setiosflags(std::ios::fixed)
145  << " energyLoss with"
146  << " momentum =" << std::setw(6) << std::setprecision(1) << trackMomentum / Units::GeV
147  << " phi =" << std::setw(6) << std::setprecision(3) << phi << " eta =" << std::setw(6) << std::setprecision(3)
148  << eta << ". CaloEnergy: deltaE = " << std::setw(8) << std::setprecision(3) << caloEnergy->deltaE() / Units::GeV
149  << " +" << std::setw(5) << std::setprecision(3) << caloEnergy->sigmaPlusDeltaE() / Units::GeV << " -"
150  << std::setw(5) << std::setprecision(3) << caloEnergy->sigmaMinusDeltaE() / Units::GeV << " (" << std::setw(5)
151  << std::setprecision(3) << caloEnergy->sigmaDeltaE() / Units::GeV << ") GeV, CaloEnergy::Type " << eLossType);
152  }
153 
154  return caloEnergy;
155  }
156  std::unique_ptr<Trk::TrackStateOnSurface> MuidCaloEnergyTool::trackStateOnSurface(const EventContext& ctx,
157  const Trk::TrackParameters& middleParameters,
158  const Trk::TrackParameters* innerParameters,
159  const Trk::TrackParameters* outerParameters) const {
160  ATH_MSG_VERBOSE("Muon with : p = " << middleParameters.momentum().mag() / Units::GeV << " Phi = "
161  << middleParameters.position().phi() << " Eta = " << middleParameters.position().eta());
162 
163  std::unique_ptr<CaloEnergy> caloEnergy;
165  // energy deposition according to the calo measurement
166  const double eta = middleParameters.position().eta();
167  const double phi = middleParameters.position().phi();
168  const double etaEM = innerParameters ? innerParameters->position().eta() : eta;
169  const double phiEM = innerParameters ? innerParameters->position().phi() : phi;
170  const double etaHad = outerParameters ? outerParameters->position().eta() : eta;
171  const double phiHad = outerParameters ? outerParameters->position().phi() : phi;
172  std::unique_ptr<CaloMeas> caloMeas{m_caloMeasTool->energyMeasurement(ctx, etaEM, phiEM, etaHad, phiHad)};
173  if (caloMeas) { caloEnergy = measurement(ctx, middleParameters.momentum().mag(), eta, phi, *caloMeas); }
174  }
175 
176  if (!caloEnergy) {
177  // parametrized energy deposition
178  ATH_MSG_DEBUG("Calculate parametrized energy loss");
179  caloEnergy =
180  energyLoss(ctx, middleParameters.momentum().mag(), middleParameters.position().eta(), middleParameters.position().phi());
182  if (caloEnergy->deltaE() > 8. * Units::GeV && middleParameters.momentum().mag() < 500. * Units::GeV)
183  ATH_MSG_WARNING(" high parametrized energy loss: " << caloEnergy->deltaE() / Units::GeV << " GeV"
184  << " for p " << middleParameters.momentum().mag() / Units::GeV
185  << " GeV"
186  << " eta " << middleParameters.position().eta());
187  }
188  // create MEOT owning CaloEnergy
189  std::bitset<Trk::MaterialEffectsBase::NumberOfMaterialEffectsTypes> typePattern(0);
194 
195  std::unique_ptr<Trk::TrackParameters> middle_clone = middleParameters.uniqueClone();
197  const CaloEnergy::EnergyLossType eLossType = caloEnergy->energyLossType();
198  const double dE = caloEnergy->deltaE() / Units::GeV;
199  const double sigPlus_dE = caloEnergy->sigmaPlusDeltaE() / Units::GeV;
200  const double sigMinus_dE = caloEnergy->sigmaMinusDeltaE() / Units::GeV;
201  const double sig_dE = caloEnergy->sigmaDeltaE() / Units::GeV;
202  //
203  auto materialEffects =
204  std::make_unique<Trk::MaterialEffectsOnTrack>(
205  0.,
206  std::move(caloEnergy),
207  middle_clone->associatedSurface(),
208  typePattern);
209 
210  // create TSOS
211  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> pattern(0);
213 
214  // debugging
215  if (msgLvl(MSG::DEBUG)) {
216  std::string eLossTypeStr = " no Calo !!";
217  switch (eLossType) {
218  case CaloEnergy::Parametrized: eLossTypeStr = "Parametrized"; break;
219  case CaloEnergy::NotIsolated: eLossTypeStr = "NotIsolated "; break;
220  case CaloEnergy::MOP: eLossTypeStr = "MOP "; break;
221  case CaloEnergy::Tail: eLossTypeStr = "Tail "; break;
222  case CaloEnergy::FSRcandidate: eLossTypeStr = "FSRcandidate"; break;
223  default: break;
224  };
225  ATH_MSG_DEBUG(std::setiosflags(std::ios::fixed)
226  << " trackStateOnSurface with"
227  << " momentum =" << std::setw(6) << std::setprecision(1) << middleParameters.momentum().mag() / Units::GeV
228  << " phi =" << std::setw(6) << std::setprecision(3) << middleParameters.position().phi()
229  << " eta =" << std::setw(6) << std::setprecision(3) << middleParameters.position().eta()
230  << ". CaloEnergy: deltaE = " << std::setw(8) << std::setprecision(3) << dE
231  << " +" << std::setw(5) << std::setprecision(3) << sigPlus_dE << " -"
232  << std::setw(5) << std::setprecision(3) << sigMinus_dE << " (" << std::setw(5)
233  << std::setprecision(3) << sig_dE << ") GeV, CaloEnergy::Type " << eLossTypeStr);
234  }
235  return std::make_unique<Trk::TrackStateOnSurface>(
236  nullptr, std::move(middle_clone), std::move(materialEffects), pattern);
237  }
238 
239  //<<<<<< PRIVATE MEMBER FUNCTION DEFINITIONS >>>>>>
240 
241  std::unique_ptr<CaloEnergy> MuidCaloEnergyTool::measurement(const EventContext& ctx, double trackMomentum, double eta, double phi,
242  const CaloMeas& caloMeas) const {
243  // Mean Energy Loss parametrization
244  std::unique_ptr<CaloEnergy> caloParamMean{m_caloParamTool->meanParametrizedEnergy(trackMomentum, eta, phi)};
245  // Mop Energy Loss parametrization
246  std::unique_ptr<CaloEnergy> caloParamMop{m_caloParamTool->mopParametrizedEnergy(trackMomentum, eta, phi)};
247  // Mip Energy Loss
248  std::unique_ptr<CaloEnergy> caloParamMip{m_caloParamTool->meanParametrizedEnergy(10. * Units::GeV, eta, phi)};
249  // Mop Energy Loss peak
250  std::unique_ptr<CaloEnergy> mopPeak{m_caloParamTool->mopPeakEnergy(trackMomentum, eta, phi)};
251 
252  // // flag upper hemisphere cosmic
253  // bool cosmic = false;
254  // if (caloParamMop->deltaE() < 0.) cosmic = true;
255 
256  // mop energy deposition
257  double MopLoss = std::abs(caloParamMop->deltaE());
258  // mop energy deposition uncertainty
259  double MopError = mopPeak->sigmaDeltaE();
260  // mop energy deposition corrected
261  double MopLossCorrected = paramCorrection(trackMomentum, eta, MopLoss, MopError);
262 
263  // percentage of inert material
264  const double InertMaterial = m_caloParamTool->x0mapInertMaterial(eta);
265  // percentage of em calorimeter material
266  const double EmMaterial = m_caloParamTool->x0mapEmMaterial(eta);
267  // percentage of hec calorimeter material
268  const double HECMaterial = m_caloParamTool->x0mapHecMaterial(eta);
269  // correction for the inert material
270  double MaterialCorrection = InertMaterial * MopLossCorrected;
271  // scale to get mop loss in em calo
272  double MopLossEm = MopLoss * m_caloParamTool->emMopFraction(eta);
273  // fraction of Tile used for the measurement
274  const double TileMeasurementMaterial = caloMeas.Tile_SamplingFraction();
275  // fraction of LArHEC used for the measurement
276  const double LArHECMeasurementMaterial = caloMeas.LArHEC_SamplingFraction();
277  // fraction of LArEM used for the measurement
278  const double LArEmMeasurementMaterial = caloMeas.LArEM_SamplingFraction();
279  // Measured energy deposition in Tile
280  const double TileEnergy = caloMeas.Tile_EnergyMeasured();
281  // Measured energy deposition in E/M
282  const double EmEnergy = caloMeas.LArEM_EnergyMeasured();
283 
284  // Correction for forward calorimetry
285  double ForwardHECCorrection = 0.;
286  if (std::abs(eta) > 2. && caloMeas.LArHEC_EnergyMeasured() > 100.)
287  ForwardHECCorrection = (1. - LArHECMeasurementMaterial) * HECMaterial * MopLossCorrected;
288  const double LArHECEnergy = caloMeas.LArHEC_EnergyMeasured() + ForwardHECCorrection; // Measured energy deposition in LArHEC
289 
290  double TotalMeasuredEnergy = TileEnergy + EmEnergy + LArHECEnergy;
291 
292  ATH_MSG_VERBOSE("Energy Deposition:Tile= " << TileEnergy / Units::GeV << ",LArHEC= " << LArHECEnergy / Units::GeV
293  << ",EM= " << EmEnergy / Units::GeV << " ForwardHECCorrection "
294  << ForwardHECCorrection / Units::GeV << " HECMaterial " << HECMaterial
295  << " MopLossCorrected " << MopLossCorrected / Units::GeV);
296 
297  bool bHEC = false; // performed HEC measurement?
298  bool bEM = false; // performed Em measurement?
299 
300  // If muon isolated, and no significant measurement is made then use the mop parameterization, else the mean
301  if (std::abs(eta) < 1.4) {
302  if (LArHECEnergy + TileEnergy > 0.1 * MopLoss * HECMaterial) bHEC = true;
303  } else if (std::abs(eta) > 1.8) {
304  if (LArHECEnergy + TileEnergy > 0.2 * MopLoss * HECMaterial) bHEC = true;
305  } else {
306  if (LArHECEnergy + TileEnergy > 0.25 * MopLoss * HECMaterial) bHEC = true;
307  }
308  if (EmEnergy > 0.5 * MopLoss * EmMaterial) bEM = true;
309 
310  if (!bHEC) {
311  // TotalMeasuredEnergy -= TileEnergy + LArHECEnergy;
312  // MaterialCorrection += HECMaterial * MopLossCorrected;
313  }
314  if (!bEM) {
315  // TotalMeasuredEnergy -= EmEnergy;
316  // MaterialCorrection += EmMaterial * MopLossCorrected;
317  }
318 
319  double MeasCorrected = TotalMeasuredEnergy + MaterialCorrection;
320  // Need to calculate the corresponding mip energy deposition
321  // Muons of 10 GeV are already in the relativistic rise region
322  // in order to obtain the mip deposition from the mean energy deposition of 10 GeV muons
323  // should divide by approximately 1.4 (Review of Particle Physics Figure 27.3 p.243)
324  const double IonizationLoss = (1. / 1.4) * std::abs(caloParamMip->deltaE());
325 
326  double eOverMipCorrectionEm = 0.;
327  double eOverMipCorrectionHEC = 0.;
328 
329  // Etrue = emip * Emeas
330  // -DE = Emeas - Etrue = Etrue ( 1./emip -1.)
331  if (bEM) {
332  const double emipEM = 0.78;
333  eOverMipCorrectionEm = -(1. / emipEM - 1.) * IonizationLoss * EmMaterial * LArEmMeasurementMaterial;
334  if (EmEnergy + eOverMipCorrectionEm < 0.) eOverMipCorrectionEm = 0.;
335  }
336  if (bHEC) {
337  const double emipTile = 0.86;
338  const double emipLAr = 0.94;
339  const double HECEnergy = TileEnergy + LArHECEnergy;
340  const double eOverMipCorrectionTile =
341  -(1. / emipTile - 1.) * TileEnergy / HECEnergy * IonizationLoss * HECMaterial * TileMeasurementMaterial;
342  const double eOverMipCorrectionLAr =
343  -(1. / emipLAr - 1.) * LArHECEnergy / HECEnergy * IonizationLoss * HECMaterial * LArHECMeasurementMaterial;
344  eOverMipCorrectionHEC = eOverMipCorrectionTile + eOverMipCorrectionLAr;
345  if (LArHECEnergy + TileEnergy + eOverMipCorrectionHEC < 0.) eOverMipCorrectionHEC = 0.;
346  }
347  const double eOverMipCorrection = eOverMipCorrectionEm + eOverMipCorrectionHEC;
348 
349  // additional offset from high-statistics Z->mumu MC (measured by Peter K
350  // 30/11/2011)
351  constexpr double fix1FromPeter[26] = {0.424104, 0.479637, 0.483419, 0.490242, 0.52806, 0.573582, 0.822098,
352  0.767301, 0.809919, 0.658745, 0.157187, 0.413214, 0.771074, 0.61815,
353  0.350113, 0.322785, 0.479294, 0.806183, 0.822161, 0.757731, -0.0857186,
354  -0.0992693, -0.0492252, 0.0650174, 0.261538, 0.360413};
355  // (update from Peter K 09/12/2011)
356  constexpr double fix2FromPeter[26] = {-0.647703, -0.303498, -0.268645, -0.261292, -0.260152, -0.269253, -0.266212,
357  -0.240837, -0.130172, -0.111638, -0.329423, -0.321011, -0.346050, -0.305592,
358  -0.313293, -0.317111, -0.428393, -0.524839, -0.599547, -0.464013, -0.159663,
359  -0.140879, -0.0975618, 0.0225352, 0.0701925, -0.24778};
360  int ieta = static_cast<int>(std::abs(eta) / 0.10);
361  if (ieta > 25) ieta = 25;
362  double FinalMeasuredEnergy = MeasCorrected + eOverMipCorrection + (fix1FromPeter[ieta] + fix2FromPeter[ieta]) * Units::GeV;
363 
364  ATH_MSG_VERBOSE("Sum of cells " << (TileEnergy + EmEnergy + LArHECEnergy) / Units::GeV << " Total energy deposition "
365  << TotalMeasuredEnergy / Units::GeV << " corrected energy deposition " << MeasCorrected / Units::GeV
366  << " e/mip corre " << FinalMeasuredEnergy / Units::GeV << std::endl
367  << "\nFinal Energy Measurement = "
368  << FinalMeasuredEnergy / Units::GeV
369  //<< "\nMean Energy Deposition = " << MeanLoss/Units::GeV
370  //<< " - " << MeanErrorLeft/Units::GeV << " + "<< MeanErrorRight/Units::GeV
371  << "\nMop Energy Deposition = " << MopLoss / Units::GeV << " +- "
372  << MopError / Units::GeV
373  //<< "\nOld parametrization energy= " << m_caloParamOld->deltaE()/Units::GeV
374  //<< " +- " << m_caloParamOld->sigmaDeltaE()/Units::GeV
375  //<< "\nTrack Momentum = " << trackMomentum/Units::GeV
376  //<< " Eta= " << eta << " Phi= " << phi
377  << std::endl
378  << "Final Meas = " << FinalMeasuredEnergy / Units::GeV << " Mop Dep = " << MopLoss / Units::GeV
379  << " +- " << MopError / Units::GeV);
380 
381  const double HECIso = caloMeas.Tile_Isolation() + caloMeas.LArHEC_Isolation();
382  const double EmIso = caloMeas.LArEM_Isolation();
383  const double sinTheta = 1. / std::cosh(eta);
384  const double pT = trackMomentum * sinTheta * Units::MeV;
385 
386  constexpr double OneOver85 = 1. / 85;
387  constexpr double FifteenGeV = 15. * Units::GeV;
388  const double EmCut = m_emMinEnergy + OneOver85 * (pT - FifteenGeV);
389  const double HECCut = m_hecMinEnergy;
390  const double pTCut = m_minMuonPt;
391  bool PassCut = true;
392  if (m_forceIsolationFailure || EmIso > EmCut || HECIso > HECCut || pT < pTCut || FinalMeasuredEnergy < m_minFinalEnergy)
393  PassCut = false;
394 
395  int nTracks = 0;
396  // double tracksEnergy = 0.;
397  if (m_trackIsolation) {
398  std::pair<int, double> inner = m_trackIsolationTool->trackIsolation(ctx, eta, phi);
399  nTracks = inner.first;
400  // tracksEnergy = inner.second - maxP;
401  if (pT < 100. * Units::GeV && nTracks > m_maxNTracksIso) PassCut = false;
402  }
403 
404  ATH_MSG_VERBOSE("pT= " << pT / Units::GeV << ",HECIso= " << HECIso / Units::GeV << ",EmIso= " << EmIso / Units::GeV
405  << ", nTracks= " << nTracks << ",PassCut= " << PassCut);
406 
408  std::unique_ptr<CaloEnergy> caloEnergy;
409 
410  // choose between lossTypes MOP, Tail, FSR and NotIsolated according
411  // to measured energy, isolation cut and Et in em
412  if (FinalMeasuredEnergy < MopLoss + 2. * MopError && FinalMeasuredEnergy > m_minFinalEnergy) {
413  ++m_countMop;
414  caloEnergy.swap(mopPeak);
415  } else if (PassCut) {
416  // tail offset from high-statistics Z->mumu MC (measured by Peter K 09/12/2011),
417  // but next we try to separate any FSR contribution from the Landau tail
418  double F1 = 0.;
419  if (caloMeas.LArEM_EnergyMeasured() > m_emEtCut) F1 = caloMeas.LArEM_FirstCompartmentEnergy() / caloMeas.LArEM_EnergyMeasured();
420  ATH_MSG_VERBOSE(" start Tail and FSR treatment: Et in e.m. " << EmEnergy * sinTheta / Units::GeV << " F1 ratio " << F1);
421  if (!m_FSRtreatment || EmEnergy * sinTheta < m_emEtCut || F1 < m_emF1Cut) {
423  double FinalEnergyErrorPlus = 0.50 * std::sqrt(FinalMeasuredEnergy / Units::GeV) * Units::GeV;
424  double FinalEnergyErrorMinus = FinalEnergyErrorPlus;
425 
426  // overall also have 50% resolution in EC rather than the 70% naively expected from LArHEC
427  if (LArHECEnergy > 1. * Units::GeV) {
428  FinalEnergyErrorMinus = FinalEnergyErrorPlus = 0.50 * std::sqrt(FinalMeasuredEnergy / Units::GeV) * Units::GeV;
429  }
430  double FinalEnergyError = 0.5 * (FinalEnergyErrorMinus + FinalEnergyErrorPlus);
431  lossType = CaloEnergy::Tail;
432  caloEnergy = std::make_unique<CaloEnergy>(FinalMeasuredEnergy, FinalEnergyError, FinalEnergyErrorMinus,
433  FinalEnergyErrorPlus, lossType);
434  ATH_MSG_VERBOSE(" CaloEnergy Tail : " << FinalMeasuredEnergy);
435  } else {
436  // significant e.m. deposit and high F1
437  double MopErrorEm = MopError * m_caloParamTool->emMopFraction(eta);
438  double FinalMeasuredEnergyNoEm = FinalMeasuredEnergy - EmEnergy + MopLossEm;
439  // lossType = CaloEnergy::FSRcandidate;
440 
441  ATH_MSG_VERBOSE(" CaloEnergy FSR : EmEnergy " << EmEnergy << " FinalMeasuredEnergy " << FinalMeasuredEnergy
442  << " MopLossEm " << MopLossEm << " MopErrorEm " << MopErrorEm
443  << " FinalMeasuredEnergyNoEm " << FinalMeasuredEnergyNoEm);
444  if (FinalMeasuredEnergyNoEm < MopLoss + 2. * MopError) {
445  // small hadronic energy deposit: like NotIsolated (MOP or Mean according to configuration)
446  lossType = CaloEnergy::NotIsolated;
447  if (m_MOPparametrization) {
448  ++m_countMop;
449  caloEnergy =
450  std::make_unique<CaloEnergy>(caloParamMop->deltaE(), caloParamMop->sigmaDeltaE(),
451  caloParamMop->sigmaMinusDeltaE(), caloParamMop->sigmaPlusDeltaE(), lossType);
452  } else {
453  ++m_countMean;
454  caloEnergy =
455  std::make_unique<CaloEnergy>(caloParamMean->deltaE(), caloParamMean->sigmaDeltaE(),
456  caloParamMean->sigmaMinusDeltaE(), caloParamMean->sigmaPlusDeltaE(), lossType);
457  }
458  ATH_MSG_VERBOSE(" CaloEnergy FSR : Small deposit, FinalMeasuredEnergyNoEm " << FinalMeasuredEnergyNoEm
459  << " using Eloss " << caloEnergy->deltaE());
460  } else {
461  // significant hadronic energy deposit
463  lossType = CaloEnergy::FSRcandidate;
464  double FinalEnergyErrorNoEm = 0.50 * std::sqrt(FinalMeasuredEnergyNoEm / Units::GeV) * Units::GeV;
465  FinalEnergyErrorNoEm = std::hypot(FinalEnergyErrorNoEm, MopErrorEm);
466  caloEnergy = std::make_unique<CaloEnergy>(FinalMeasuredEnergyNoEm, FinalEnergyErrorNoEm, FinalEnergyErrorNoEm,
467  FinalEnergyErrorNoEm, lossType);
468  ATH_MSG_VERBOSE(" CaloEnergy FSR : Large deposit, FinalMeasuredEnergyNoEm " << FinalMeasuredEnergyNoEm);
469  }
470  }
471  } else {
472  // lossType NotIsolated: MOP or Mean according to configuration
473  if (m_MOPparametrization) {
474  ++m_countMop;
475  caloEnergy = std::make_unique<CaloEnergy>(caloParamMop->deltaE(), caloParamMop->sigmaDeltaE(),
476  caloParamMop->sigmaMinusDeltaE(), caloParamMop->sigmaPlusDeltaE(), lossType);
477  } else {
478  ++m_countMean;
479  caloEnergy = std::make_unique<CaloEnergy>(caloParamMean->deltaE(), caloParamMean->sigmaDeltaE(),
480  caloParamMean->sigmaMinusDeltaE(), caloParamMean->sigmaPlusDeltaE(), lossType);
481  }
482  }
483 
484  caloEnergy->set_fsrCandidateEnergy(MopLossEm);
485  return caloEnergy;
486  }
487 
489  const double pT = trackMomentum / std::cosh(eta) / Units::GeV; // pt in GeV
490  double a = 0.;
491  double b = 0.;
492  if (std::abs(eta) < 1.) {
493  a = 0.02255;
494  b = 7.708e-5;
495  } else if (std::abs(eta) > 1. && std::abs(eta) < 2.) {
496  a = 0.04198;
497  b = 8.912e-5;
498  } else {
499  a = 0.02181;
500  b = 7.282e-5;
501  }
502  return std::hypot(a, (b * pT));
503  }
504 
505  double MuidCaloEnergyTool::paramCorrection(double trackMomentum, double eta, double MopLoss, double MopLossSigma) {
506  const double Nsigma = 2.;
507  double MSres = muSpecResolParam(trackMomentum, eta);
508  double MSsigma = trackMomentum * MSres;
509  double sigma = std::hypot(MSsigma, MopLossSigma);
510  double sum = 0.;
511  double weight = 0.;
512  double xlow = MopLoss - Nsigma * sigma;
513  if (xlow < 0.) xlow = 0.;
514  double xup = MopLoss + Nsigma * sigma;
515  int Na = 50;
516  const double inv_Na = 1. / static_cast<double>(Na);
517  for (int j = 0; j < Na; ++j) {
518  double x = xlow + j * (xup - xlow) * inv_Na;
519  double w = landau(x, MopLoss, MopLossSigma, true);
520  sum += x * w;
521  weight += w;
522  }
523  double MopStat = sum / weight;
524  return MopStat;
525  }
526 
527  double MuidCaloEnergyTool::landau(double x, double mpv, double sigma, bool norm) {
528  // The LANDAU function with mpv(most probable value) and sigma.
529  // This function has been adapted from the CERNLIB routine G110 denlan.
530  // If norm=kTRUE (default is kFALSE) the result is divided by sigma
531 
532  constexpr double p1[5] = {0.4259894875, -0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
533  constexpr double q1[5] = {1.0, -0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
534 
535  constexpr double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
536  constexpr double q2[5] = {1.0, 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
537 
538  constexpr double p3[5] = {0.1788544503, 0.09359161662, 0.006325387654, 0.00006611667319, -0.000002031049101};
539  constexpr double q3[5] = {1.0, 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
540 
541  constexpr double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
542  constexpr double q4[5] = {1.0, 106.8615961, 337.6496214, 2016.712389, 1597.063511};
543 
544  constexpr double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
545  constexpr double q5[5] = {1.0, 156.9424537, 3745.310488, 9834.698876, 66924.28357};
546 
547  constexpr double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
548  constexpr double q6[5] = {1.0, 651.4101098, 56974.73333, 165917.4725, -2815759.939};
549 
550  constexpr double a1[3] = {0.04166666667, -0.01996527778, 0.02709538966};
551 
552  constexpr double a2[2] = {-1.845568670, -4.284640743};
553 
554  if (sigma <= 0) return 0.;
555  double v = (x - mpv) / sigma;
556  double u, ue, us, den;
557  if (v < -5.5) {
558  u = exp(v + 1.0);
559  if (u < 1e-10) return 0.0;
560  ue = std::exp(-1 / u);
561  us = std::sqrt(u);
562  den = 0.3989422803 * (ue / us) * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
563  } else if (v < -1) {
564  u = std::exp(-v - 1);
565  den = exp(-u) * std::sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) /
566  (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
567  } else if (v < 1) {
568  den = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) /
569  (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
570  } else if (v < 5) {
571  den = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) /
572  (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
573  } else if (v < 12) {
574  u = 1 / v;
575  den = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) /
576  (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
577  } else if (v < 50) {
578  u = 1 / v;
579  den = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) /
580  (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
581  } else if (v < 300) {
582  u = 1 / v;
583  den = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) /
584  (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
585  } else {
586  u = 1 / (v - v * log(v) / (v + 1));
587  den = u * u * (1 + (a2[0] + a2[1] * u) * u);
588  }
589  if (!norm) return den;
590  return den / sigma;
591  }
592 
593 } // namespace Rec
str::trackMomentum
const std::string trackMomentum
Definition: BTagTrackIpAccessor.cxx:15
Rec::MuidCaloEnergyTool::m_forceIsolationFailure
bool m_forceIsolationFailure
Definition: MuidCaloEnergyTool.h:89
CaloEnergy::Tail
@ Tail
Definition: CaloEnergy.h:43
mergePhysValFiles.pattern
pattern
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:25
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::EnergyLoss::sigmaMinusDeltaE
double sigmaMinusDeltaE() const
returns the negative side
python.HION7.pTCut
int pTCut
Definition: HION7.py:72
CalculateHighPtTerm.pT
pT
Definition: ICHEP2016/CalculateHighPtTerm.py:57
Trk::TrackStateOnSurface::CaloDeposit
@ CaloDeposit
This TSOS contains a CaloEnergy object.
Definition: TrackStateOnSurface.h:135
EnergyLoss.h
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:18
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
Rec::MuidCaloEnergyTool::paramCorrection
static double paramCorrection(double trackMomentum, double eta, double MopLoss, double MopSigma)
Definition: MuidCaloEnergyTool.cxx:505
Rec::CaloMeas::Tile_Isolation
double Tile_Isolation(void) const
Definition: CaloMeas.h:64
CaloEnergy::EnergyLossType
EnergyLossType
Calo Energy Loss Type Parametrized : reconstruction configured to use the parametrization w/o looking...
Definition: CaloEnergy.h:43
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Rec::MuidCaloEnergyTool::trackStateOnSurface
std::unique_ptr< Trk::TrackStateOnSurface > trackStateOnSurface(const EventContext &ctx, const Trk::TrackParameters &middleParameters, const Trk::TrackParameters *innerParameters, const Trk::TrackParameters *outerParameters) const override
IMuidCaloEnergy interface: TrackStateOnSurface for parameters and energyLoss at the calorimeter mid-s...
Definition: MuidCaloEnergyTool.cxx:156
Rec::MuidCaloEnergyTool::m_minFinalEnergy
double m_minFinalEnergy
Definition: MuidCaloEnergyTool.h:100
OfflineHitType::InertMaterial
@ InertMaterial
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Rec::MuidCaloEnergyTool::m_countMean
std::atomic_int m_countMean
Definition: MuidCaloEnergyTool.h:104
Trk::ParametersBase::associatedSurface
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
Trk::ParametersBase::uniqueClone
std::unique_ptr< ParametersBase< DIM, T > > uniqueClone() const
clone method for polymorphic deep copy returning unique_ptr; it is not overriden, but uses the existi...
Definition: ParametersBase.h:97
Rec::CaloMeas::LArEM_EnergyMeasured
double LArEM_EnergyMeasured(void) const
Definition: CaloMeas.h:84
Rec::MuidCaloEnergyTool::m_caloMeasTool
ToolHandle< IMuidCaloEnergyMeas > m_caloMeasTool
Definition: MuidCaloEnergyTool.h:70
Rec::MuidCaloEnergyTool::m_emEtCut
double m_emEtCut
Definition: MuidCaloEnergyTool.h:95
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
Trk::EnergyLoss::sigmaDeltaE
double sigmaDeltaE() const
returns the symmatric error
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
Rec::MuidCaloEnergyTool::m_energyLossMeasurement
bool m_energyLossMeasurement
Definition: MuidCaloEnergyTool.h:88
Rec::MuidCaloEnergyTool::m_MOPparametrization
bool m_MOPparametrization
Definition: MuidCaloEnergyTool.h:91
Rec::MuidCaloEnergyTool::muSpecResolParam
static double muSpecResolParam(double trackMomentum, double eta)
Definition: MuidCaloEnergyTool.cxx:488
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
Rec::CaloMeas::LArHEC_EnergyMeasured
double LArHEC_EnergyMeasured(void) const
Definition: CaloMeas.h:72
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:77
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
Rec::MuidCaloEnergyTool::m_emF1Cut
double m_emF1Cut
Definition: MuidCaloEnergyTool.h:96
python.SystemOfUnits.MeV
float MeV
Definition: SystemOfUnits.py:172
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
Rec::MuidCaloEnergyTool::initialize
StatusCode initialize() override
Definition: MuidCaloEnergyTool.cxx:74
IMuidTrackIsolation.h
MaterialEffectsOnTrack.h
Rec::MuidCaloEnergyTool::finalize
StatusCode finalize() override
Definition: MuidCaloEnergyTool.cxx:104
Rec::MuidCaloEnergyTool::m_cosmics
bool m_cosmics
Definition: MuidCaloEnergyTool.h:87
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
Rec::CaloMeas::LArEM_FirstCompartmentEnergy
double LArEM_FirstCompartmentEnergy(void) const
Definition: CaloMeas.h:88
python.SystemOfUnits.us
float us
Definition: SystemOfUnits.py:149
CaloEnergy::energyLossType
CaloEnergy::EnergyLossType energyLossType(void) const
Accessor methods.
Definition: CaloEnergy.h:162
Rec::MuidCaloEnergyTool::m_maxNTracksIso
int m_maxNTracksIso
Definition: MuidCaloEnergyTool.h:99
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
Rec::CaloMeas
Definition: CaloMeas.h:19
Rec::MuidCaloEnergyTool::landau
static double landau(double x, double mpv, double sigma, bool norm)
Definition: MuidCaloEnergyTool.cxx:527
Rec
Name: MuonSpContainer.h Package : offline/Reconstruction/MuonIdentification/muonEvent.
Definition: FakeTrackBuilder.h:10
Rec::CaloMeas::LArEM_Isolation
double LArEM_Isolation(void) const
Definition: CaloMeas.h:92
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
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition: AthCommonDataStore.h:145
Rec::MuidCaloEnergyTool::m_caloParamTool
ToolHandle< IMuidCaloEnergyParam > m_caloParamTool
Definition: MuidCaloEnergyTool.h:75
Trk::EnergyLoss::deltaE
double deltaE() const
returns the
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Rec::MuidCaloEnergyTool::m_FSRtreatment
bool m_FSRtreatment
Definition: MuidCaloEnergyTool.h:90
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Rec::MuidCaloEnergyTool::measurement
std::unique_ptr< CaloEnergy > measurement(const EventContext &ctx, double trackMomentum, double eta, double phi, const CaloMeas &caloMeas) const
Definition: MuidCaloEnergyTool.cxx:241
Trk::ParametersBase
Definition: ParametersBase.h:55
CaloEnergy::FSRcandidate
@ FSRcandidate
Definition: CaloEnergy.h:43
MuidCaloEnergyTool.h
CaloMeas.h
Athena::Units
Definition: Units.h:45
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
Rec::CaloMeas::Tile_EnergyMeasured
double Tile_EnergyMeasured(void) const
Definition: CaloMeas.h:60
IMuidCaloEnergyParam.h
Rec::CaloMeas::LArHEC_Isolation
double LArHEC_Isolation(void) const
Definition: CaloMeas.h:76
Rec::MuidCaloEnergyTool::m_countMeasurement
std::atomic_int m_countMeasurement
Definition: MuidCaloEnergyTool.h:105
Units.h
Wrapper to avoid constant divisions when using units.
Rec::CaloMeas::Tile_SamplingFraction
double Tile_SamplingFraction(void) const
Definition: CaloMeas.h:68
Trk::MaterialEffectsBase::EnergyLossEffects
@ EnergyLossEffects
contains energy loss corrections
Definition: MaterialEffectsBase.h:48
python.PyAthena.v
v
Definition: PyAthena.py:154
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Rec::MuidCaloEnergyTool::m_trackIsolationTool
ToolHandle< IMuidTrackIsolation > m_trackIsolationTool
Definition: MuidCaloEnergyTool.h:80
CaloEnergy::set_fsrCandidateEnergy
void set_fsrCandidateEnergy(const float fs)
FSR Candidate Energy.
Definition: CaloEnergy.h:118
a
TList * a
Definition: liststreamerinfos.cxx:10
CaloEnergy::Parametrized
@ Parametrized
Definition: CaloEnergy.h:43
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
DEBUG
#define DEBUG
Definition: page_access.h:11
CaloEnergy::NotIsolated
@ NotIsolated
Definition: CaloEnergy.h:43
Rec::MuidCaloEnergyTool::MuidCaloEnergyTool
MuidCaloEnergyTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: MuidCaloEnergyTool.cxx:39
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
Rec::MuidCaloEnergyTool::m_countMop
std::atomic_int m_countMop
Definition: MuidCaloEnergyTool.h:106
CaloEnergy::MOP
@ MOP
Definition: CaloEnergy.h:43
AthAlgTool
Definition: AthAlgTool.h:26
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:198
Rec::MuidCaloEnergyTool::~MuidCaloEnergyTool
~MuidCaloEnergyTool(void)
TRTCalib_cfilter.p3
p3
Definition: TRTCalib_cfilter.py:132
IMuidCaloEnergyMeas.h
Rec::CaloMeas::LArEM_SamplingFraction
double LArEM_SamplingFraction(void) const
Definition: CaloMeas.h:96
Rec::MuidCaloEnergyTool::m_minMuonPt
double m_minMuonPt
Definition: MuidCaloEnergyTool.h:101
Trk::EnergyLoss::sigmaPlusDeltaE
double sigmaPlusDeltaE() const
returns the positive side
Rec::MuidCaloEnergyTool::m_trackIsolation
bool m_trackIsolation
Definition: MuidCaloEnergyTool.h:92
TrackStateOnSurface.h
Rec::MuidCaloEnergyTool::m_hecMinEnergy
double m_hecMinEnergy
Definition: MuidCaloEnergyTool.h:98
Rec::CaloMeas::LArHEC_SamplingFraction
double LArHEC_SamplingFraction(void) const
Definition: CaloMeas.h:80
Rec::MuidCaloEnergyTool::m_emMinEnergy
double m_emMinEnergy
Definition: MuidCaloEnergyTool.h:97
Rec::MuidCaloEnergyTool::energyLoss
std::unique_ptr< CaloEnergy > energyLoss(const EventContext &ctx, double trackMomentum, double eta, double phi) const override
IMuidCaloEnergy interface: to get the total energyLoss in the calorimeters.
Definition: MuidCaloEnergyTool.cxx:111
CaloEnergy.h