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