ATLAS Offline Software
Loading...
Searching...
No Matches
MuonCalibIntScaleSmearTool.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5// Framework include(s):
6#include "AsgDataHandles/ReadHandle.h"
7#include "PathResolver/PathResolver.h"
8
9// Local include(s):
10#include <cmath>
11#include "TRandom3.h"
12
13#include "MuonMomentumCorrections/MuonCalibTool.h"
14#include "MuonMomentumCorrections/MuonCalibIntScaleSmearTool.h"
15#include "MuonMomentumCorrections/EnumDef.h"
16#include <MuonMomentumCorrections/CalibInitializer.h>
17
18namespace ADD_CMODE(CP)
19{
20
21 MuonCalibIntScaleSmearTool::MuonCalibIntScaleSmearTool(const std::string &name) : asg::AsgTool(name)
22 {}
23
24 StatusCode MuonCalibIntScaleSmearTool::initialize()
25 {
26 // Greet the user:
27 ATH_MSG_INFO("Initializing MuonCalibIntScaleSmearTool");
28
29 // Get the m_eventinfo container
30 ATH_CHECK(m_eventInfo.initialize());
31
32
33 // Load the constants
34 for(const auto & year: MCP::dataYearList)
35 {
36 for(const auto & param: m_paramList)
37 {
38 m_IDparams[year][param] = MCP::CalibInitializer::createScaleResCorrMap(year, MCP::TrackType::ID, m_release, param);
39 m_MEparams[year][param] = MCP::CalibInitializer::createScaleResCorrMap(year, MCP::TrackType::ME, m_release, param);
40 m_CBparams[year][param] = MCP::CalibInitializer::createScaleResCorrMap(year, MCP::TrackType::CB, m_release, param);
41 }
42
43 m_IDExpectedResparams[year] = MCP::CalibInitializer::createExpectedResMap(year, MCP::TrackType::ID, m_release);
44 m_MEExpectedResparams[year] = MCP::CalibInitializer::createExpectedResMap(year, MCP::TrackType::ME, m_release);
45 }
46
47
48 m_currentParameters = nullptr;
49 // Init the systematics
50 m_Parameters.initialize(affectingSystematics(), [this](const CP::SystematicSet &systConfig, ParameterSetScaleSmear &param)
51 { return calcSystematicVariation(systConfig, param); });
52 if (!applySystematicVariation(CP::SystematicSet()))
53 {
54 ATH_MSG_ERROR("Unable to run with no systematic");
55 return StatusCode::FAILURE;
56 }
57 CP::SystematicRegistry &registry = CP::SystematicRegistry::getInstance();
58 if (registry.registerSystematics(*this) != StatusCode::SUCCESS)
59 {
60 ATH_MSG_ERROR("Unkown systematic list");
61 return StatusCode::FAILURE;
62 }
63 // Return gracefully:
64 return StatusCode::SUCCESS;
65 }
66
67 CP::CorrectionCode MuonCalibIntScaleSmearTool::applyCorrection(MCP::MuonObj& mu) const
68 {
69 // Do nothing, if it is data
70 if(mu.CB.isData) return CP::CorrectionCode::Ok;
71
72 // Get the correction constants for each year
73 auto IDcorrConstants = getConstants(mu.ID, m_IDparams.at(mu.ID.year), m_currentParameters->m_IDparams);
74 auto MEcorrConstants = getConstants(mu.ME, m_MEparams.at(mu.ME.year), m_currentParameters->m_MEparams);
75 auto CBcorrConstants = getConstants(mu.CB, m_CBparams.at(mu.CB.year), m_currentParameters->m_CBparams);
76
77 double corrIDpT = getCorrectedPt(mu, mu.ID, IDcorrConstants);
78 double corrMEpT = getCorrectedPt(mu, mu.ME, MEcorrConstants);
79 double corrCBpT = getCorrectedPt(mu, mu.CB, CBcorrConstants);
80
81 //Smearing can be positive and negative since it depends on Gaussians centered in 0 with sigma 1.
82 //However, the resolution correction factor depends on 1/(1+smearing)
83 //and therefore if the smearing is less than -1 the pT turns out to have a negative value.
84 //Physically it means that the resolution has caused the pT to be reconstructed with a flip of the charge.
85 //Therefore in addition to correcting the pT, a flip of the charge is also done through the calib_charge variable.
86
87 // First compute the calibrated CB pT under the recombination scheme with uncalibrated ID & ME pT
88 double corrCBpTWithIDME = getCorrectedCBPtWithIDMSComb(mu, IDcorrConstants, MEcorrConstants);
89
90 // Write the pT into the object (if negative pT, multiply it by -1)
91 //< 0.1 because sometimes the pT is -0.
92 mu.ID.calib_pt = corrIDpT * ((corrIDpT < -0.1) ? -1 : 1);
93 mu.ME.calib_pt = corrMEpT * ((corrMEpT < -0.1) ? -1 : 1);
94 mu.CB.calib_pt = corrCBpT * ((corrCBpT < -0.1) ? -1 : 1);
95 //charge calibration (flip of the charge if needed)
96 mu.ID.calib_charge = mu.ID.uncalib_charge * ((corrIDpT < -0.1) ? -1 : 1);
97 mu.ME.calib_charge = mu.ME.uncalib_charge * ((corrMEpT < -0.1) ? -1 : 1);
98 mu.CB.calib_charge = mu.CB.uncalib_charge * ((corrCBpT < -0.1) ? -1 : 1);
99
100 if(m_calibMode == MuonCalibTool::correctData_IDMS || m_calibMode == MuonCalibTool::notCorrectData_IDMS) {
101 mu.CB.calib_pt = corrCBpTWithIDME * ((corrCBpTWithIDME < -0.1) ? -1 : 1);
102 mu.CB.calib_charge = mu.CB.uncalib_charge * ((corrCBpTWithIDME < -0.1) ? -1 : 1);
103 }
104 else if(m_calibMode == MuonCalibTool::correctData_IDonly) {
105 mu.CB.calib_pt = mu.ID.calib_pt;
106 mu.CB.calib_charge = mu.ID.calib_charge;
107 }
108 else if(m_calibMode == MuonCalibTool::correctData_MSonly) {
109 mu.CB.calib_pt = mu.ME.calib_pt;
110 mu.CB.calib_charge = mu.ME.calib_charge;
111 }
112
113 // Return gracefully:
114 return CP::CorrectionCode::Ok;
115 }
116
117 std::map<MCP::ScaleSmearParam, double> MuonCalibIntScaleSmearTool::getConstants(const MCP::TrackCalibObj& trk, const std::map<MCP::ScaleSmearParam, ScaleSmearCorrConstMap>& constants, const std::map<MCP::ScaleSmearParam, double>& direction) const
118 {
119
120 std::map<MCP::ScaleSmearParam, double> calibConstants;
121
122 // Extra the constants from container into a simple map
123 for(const auto& param: m_paramList)
124 {
125 const auto & contantList = constants.at(param);
126
127 double val = contantList.at(MCP::ScaleResCorrection::Nominal)->getCalibConstant(trk);
128 if(direction.at(param) == MCP::SystVariation::Up)
129 {
130 val += contantList.at(MCP::ScaleResCorrection::SystErr__1up)->getCalibConstant(trk);
131 }
132 else if (direction.at(param) == MCP::SystVariation::Down)
133 {
134 val -= contantList.at(MCP::ScaleResCorrection::SystErr__1down)->getCalibConstant(trk);
135 }
136
137 // Remove unphysical smearing
138 if(param == MCP::ScaleSmearParam::r0 || param == MCP::ScaleSmearParam::r1 || param == MCP::ScaleSmearParam::r2)
139 {
140 val = std::max(0.0, val);
141 }
142
143 calibConstants[param] = val;
144 }
145
146 return calibConstants;
147 }
148
149
150 double MuonCalibIntScaleSmearTool::getCorrectedPt(const MCP::MuonObj& mu, const MCP::TrackCalibObj& trk, const std::map<MCP::ScaleSmearParam, double>& calibConstant) const
151 {
152 if(trk.calib_pt == 0) return 0;
153
154 // For debugging:: Todo add checking if in verbose mode
155 ATH_MSG_VERBOSE("MuonType: "<<MCP::toString(trk.type));
156 for(const auto& var: calibConstant) ATH_MSG_VERBOSE("var: "<<MCP::toString(var.first)<<" = "<<var.second);
157
158 double pT = trk.calib_pt;
159
160 // Calculate the denominator for the scale correction
161 double smearCorr = getSmearCorr(mu, trk, calibConstant.at(MCP::ScaleSmearParam::r0), calibConstant.at(MCP::ScaleSmearParam::r1), calibConstant.at(MCP::ScaleSmearParam::r2));
162 ATH_MSG_VERBOSE("smearCorr: "<<smearCorr);
163
164 // apply the full correction
165 double corrpT = (pT + pT * calibConstant.at(MCP::ScaleSmearParam::s1) + calibConstant.at(MCP::ScaleSmearParam::s0))/(1 + smearCorr);
166
167 return corrpT;
168 }
169
170 double MuonCalibIntScaleSmearTool::getSmearCorr(const MCP::MuonObj& mu, const MCP::TrackCalibObj& trk, double r0, double r1, double r2) const
171 {
172 if(trk.calib_pt == 0) return 1;
173
174 double pT = trk.calib_pt;
175 if(trk.type == MCP::TrackType::ID)
176 {
177 // ID applies a tan(theta) correction for r2 for high eta muons
178 double additional_weight = 1.;
179 if (std::abs(trk.eta) > 2) additional_weight = sinh(trk.eta);
180
181 return r1 * mu.rnd_g3 + r2 * mu.rnd_g4 * pT * additional_weight;
182 }
183 return r0 * mu.rnd_g0 / pT + r1 * mu.rnd_g1 + r2 * mu.rnd_g2 * pT;
184 }
185
186
187 double MuonCalibIntScaleSmearTool::getCorrectedCBPtWithIDMSComb(const MCP::MuonObj& mu, const std::map<MCP::ScaleSmearParam, double>& calibIDConstant, const std::map<MCP::ScaleSmearParam, double>& calibMEConstant) const
188 {
189 // calculate the relative of ID and ME for the corrections
190 double weightID = 0.5;
191 double weightME = 0.5;
192
193 double deltaCBME = mu.CB.calib_pt - mu.ME.calib_pt;
194 double deltaCBID = mu.CB.calib_pt - mu.ID.calib_pt;
195
196 if (mu.ME.calib_pt == 0)
197 {
198 weightID = 1.0;
199 weightME = 0.0;
200 }
201 else if (mu.ID.calib_pt == 0)
202 {
203 weightID = 0.0;
204 weightME = 1.0;
205 }
206 else if (mu.CB.calib_pt != 0)
207 {
208 if (std::abs(deltaCBME) > 0 || std::abs(deltaCBID) > 0)
209 {
210 double R = 1, Rplus = 1;
211 if (std::abs(deltaCBME) == std::abs(deltaCBID))
212 {
213 // do nothing
214 }
215 else if (std::abs(deltaCBME) != 0 &&
216 std::abs(deltaCBME) > std::abs(deltaCBID)) {
217 R = (-deltaCBID) / deltaCBME; /* R~wMS/wID */
218 Rplus = 1 + R;
219 if (Rplus != 0 && R > 0) {
220 weightID = 1 / Rplus;
221 weightME = R / Rplus;
222 }
223 }
224 else if (std::abs(deltaCBID) != 0 &&
225 std::abs(deltaCBME) < std::abs(deltaCBID)) {
226 R = (-deltaCBME) / (deltaCBID); /* R~wID/wMS */
227 Rplus = 1 + R;
228 if (Rplus != 0 && R > 0) {
229 weightID = R / Rplus;
230 weightME = 1 / Rplus;
231 }
232 }
233 }
234 }
235 // If no correction was found
236 if(weightID == 0.5 && weightME == 0.5)
237 {
238 if(mu.expectedPercentResME<std::numeric_limits<double>::epsilon() ||
239 mu.CB.calib_pt<std::numeric_limits<double>::epsilon() || mu.expectedPercentResID<std::numeric_limits<double>::epsilon() ){
240 ATH_MSG_VERBOSE("Potential FPE caught! Continue with averaging ID+MS");
241 ATH_MSG_VERBOSE("mu.expectedPercentResME = "<<mu.expectedPercentResME);
242 ATH_MSG_VERBOSE("mu.expectedPercentResID = "<<mu.expectedPercentResID);
243 ATH_MSG_VERBOSE("mu.CB.calib_pt = "<<mu.CB.calib_pt);
244 } else {
245 double wME = mu.ME.calib_pt / mu.CB.calib_pt / std::pow(mu.expectedPercentResME, 2);
246 double wID = mu.ID.calib_pt / mu.CB.calib_pt / std::pow(mu.expectedPercentResID, 2);
247 weightID = wID / (wME + wID);
248 weightME = wME / (wME + wID);
249 }
250 }
251
252
253 // Calculate the denominator for the scale correction
254 double smearIDCorr = getSmearCorr(mu, mu.ID, calibIDConstant.at(MCP::ScaleSmearParam::r0), calibIDConstant.at(MCP::ScaleSmearParam::r1), calibIDConstant.at(MCP::ScaleSmearParam::r2));
255 double smearMECorr = getSmearCorr(mu, mu.ME, calibMEConstant.at(MCP::ScaleSmearParam::r0), calibMEConstant.at(MCP::ScaleSmearParam::r1), calibMEConstant.at(MCP::ScaleSmearParam::r2));
256 double smearCorr = weightID * smearIDCorr + weightME * smearMECorr;
257 double scaleCB = 0;
258
259 // apply the full correction
260 if(weightID == 0) scaleCB = (calibMEConstant.at(MCP::ScaleSmearParam::s1) + calibMEConstant.at(MCP::ScaleSmearParam::s0) / mu.ME.calib_pt) * weightME;
261 else if(weightME == 0) scaleCB = calibIDConstant.at(MCP::ScaleSmearParam::s1) * weightID ;
262 else scaleCB = (calibIDConstant.at(MCP::ScaleSmearParam::s1) * weightID + (calibMEConstant.at(MCP::ScaleSmearParam::s1) + calibMEConstant.at(MCP::ScaleSmearParam::s0) / mu.ME.calib_pt) * weightME);
263
264 ATH_MSG_VERBOSE("mu.ID.calib_pt: "<<mu.ID.calib_pt);
265 ATH_MSG_VERBOSE("mu.ME.calib_pt: "<<mu.ME.calib_pt);
266 ATH_MSG_VERBOSE("mu.CB.calib_pt: "<<mu.CB.calib_pt);
267 ATH_MSG_VERBOSE("mu.expectedPercentResME: "<<mu.expectedPercentResME);
268 ATH_MSG_VERBOSE("mu.expectedPercentResID: "<<mu.expectedPercentResID);
269
270 ATH_MSG_VERBOSE("weightID: "<<weightID);
271 ATH_MSG_VERBOSE("weightME: "<<weightME);
272 ATH_MSG_VERBOSE("smearIDCorr: "<<smearIDCorr);
273 ATH_MSG_VERBOSE("smearMECorr: "<<smearMECorr);
274 ATH_MSG_VERBOSE("smearCorr: "<<smearCorr);
275 ATH_MSG_VERBOSE("scaleCB: "<<scaleCB);
276
277 double pT = mu.CB.calib_pt;
278 double corrpT = (pT + pT * scaleCB)/(1 + smearCorr);
279
280 return corrpT;
281 }
282
283
284 bool MuonCalibIntScaleSmearTool::isAffectedBySystematic(const CP::SystematicVariation &systematic) const
285 {
286 CP::SystematicSet sys = affectingSystematics();
287 return sys.find(systematic) != sys.end();
288 }
289
290 CP::SystematicSet MuonCalibIntScaleSmearTool::affectingSystematics() const
291 {
292 CP::SystematicSet result;
294 // Resolution systematics
296 if (m_calibMode == MuonCalibTool::correctData_CB || m_calibMode == MuonCalibTool::notCorrectData_CB || m_sysScheme == "AllSys") {
297 // CB systematics
298 result.insert(CP::SystematicVariation("MUON_CB", 1));
299 result.insert(CP::SystematicVariation("MUON_CB", -1));
300 }
301 if (m_calibMode == MuonCalibTool::correctData_IDMS || m_calibMode == MuonCalibTool::notCorrectData_IDMS || m_calibMode == MuonCalibTool::correctData_IDonly || m_sysScheme == "AllSys") {
302 // ID systematics
303 result.insert(CP::SystematicVariation("MUON_ID", 1));
304 result.insert(CP::SystematicVariation("MUON_ID", -1));
305 }
306 if (m_calibMode == MuonCalibTool::correctData_IDMS || m_calibMode == MuonCalibTool::notCorrectData_IDMS || m_calibMode == MuonCalibTool::correctData_MSonly || m_sysScheme == "AllSys") {
307 // MS systematics
308 result.insert(CP::SystematicVariation("MUON_MS", 1));
309 result.insert(CP::SystematicVariation("MUON_MS", -1));
310 }
311
315 if (m_sysScheme == "Corr_Scale") {
316 result.insert(CP::SystematicVariation("MUON_SCALE", 1));
317 result.insert(CP::SystematicVariation("MUON_SCALE", -1));
318 }
319 else if (m_sysScheme == "Decorr_Scale" || m_sysScheme == "AllSys") {
320 // Either doing direct calib of CB or asking for all the sys
321 if (m_calibMode == MuonCalibTool::correctData_CB || m_calibMode == MuonCalibTool::notCorrectData_CB || m_sysScheme == "AllSys") {
322 result.insert(CP::SystematicVariation("MUON_SCALE_CB", 1));
323 result.insert(CP::SystematicVariation("MUON_SCALE_CB", -1));
324
325 result.insert(CP::SystematicVariation("MUON_SCALE_CB_ELOSS", 1));
326 result.insert(CP::SystematicVariation("MUON_SCALE_CB_ELOSS", -1));
327 }
328
329 // Either not doing direct calib of CB or asking for IDonly calib or asking for all the sys
330 if (m_calibMode == MuonCalibTool::correctData_IDMS || m_calibMode == MuonCalibTool::notCorrectData_IDMS || m_calibMode == MuonCalibTool::correctData_IDonly || m_sysScheme == "AllSys") {
331 result.insert(CP::SystematicVariation("MUON_SCALE_ID", 1));
332 result.insert(CP::SystematicVariation("MUON_SCALE_ID", -1));
333 }
334
335 // Either not doing direct calib of CB or asking for MEonly calib or asking for all the sys
336 if (m_calibMode == MuonCalibTool::correctData_IDMS || m_calibMode == MuonCalibTool::notCorrectData_IDMS || m_calibMode == MuonCalibTool::correctData_MSonly || m_sysScheme == "AllSys") {
337 result.insert(CP::SystematicVariation("MUON_SCALE_MS", 1));
338 result.insert(CP::SystematicVariation("MUON_SCALE_MS", -1));
339
340 result.insert(CP::SystematicVariation("MUON_SCALE_MS_ELOSS", 1));
341 result.insert(CP::SystematicVariation("MUON_SCALE_MS_ELOSS", -1));
342 }
343 }
344
345
346
347
348 return result;
349 }
350
351 CP::SystematicSet MuonCalibIntScaleSmearTool::recommendedSystematics() const { return affectingSystematics(); }
352
353 StatusCode MuonCalibIntScaleSmearTool::calcSystematicVariation(const CP::SystematicSet &systConfig, ParameterSetScaleSmear &param) const
354 {
355 param.m_IDparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Default;
356 param.m_IDparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Default;
357 param.m_IDparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Default;
358 param.m_IDparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Default;
359 param.m_IDparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Default;
360
361 param.m_MEparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Default;
362 param.m_MEparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Default;
363 param.m_MEparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Default;
364 param.m_MEparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Default;
365 param.m_MEparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Default;
366
367 param.m_CBparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Default;
368 param.m_CBparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Default;
369 param.m_CBparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Default;
370 param.m_CBparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Default;
371 param.m_CBparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Default;
372
373
374 // ID systematics
375 CP::SystematicVariation syst = systConfig.getSystematicByBaseName("MUON_ID");
376
377 if (syst == CP::SystematicVariation("MUON_ID", 1)) {
378 param.m_IDparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Up;
379 param.m_IDparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Up;
380 param.m_IDparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Up;
381 } else if (syst == CP::SystematicVariation("MUON_ID", -1)) {
382 param.m_IDparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Down;
383 param.m_IDparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Down;
384 param.m_IDparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Down;
385 } else if (!syst.empty()) return StatusCode::FAILURE;
386
387 // MS systematics
388 syst = systConfig.getSystematicByBaseName("MUON_MS");
389
390 if (syst == CP::SystematicVariation("MUON_MS", 1)) {
391 param.m_MEparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Up;
392 param.m_MEparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Up;
393 param.m_MEparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Up;
394 } else if (syst == CP::SystematicVariation("MUON_MS", -1)) {
395 param.m_MEparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Down;
396 param.m_MEparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Down;
397 param.m_MEparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Down;
398 } else if (!syst.empty()) return StatusCode::FAILURE;
399
400 // CB systematics
401 syst = systConfig.getSystematicByBaseName("MUON_CB");
402
403 if (syst == CP::SystematicVariation("MUON_CB", 1)) {
404 param.m_CBparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Up;
405 param.m_CBparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Up;
406 param.m_CBparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Up;
407 } else if (syst == CP::SystematicVariation("MUON_CB", -1)) {
408 param.m_CBparams[MCP::ScaleSmearParam::r0] = MCP::SystVariation::Down;
409 param.m_CBparams[MCP::ScaleSmearParam::r1] = MCP::SystVariation::Down;
410 param.m_CBparams[MCP::ScaleSmearParam::r2] = MCP::SystVariation::Down;
411 } else if (!syst.empty()) return StatusCode::FAILURE;
412
413 // Scale systematics
414 syst = systConfig.getSystematicByBaseName("MUON_SCALE");
415
416 if (syst == CP::SystematicVariation("MUON_SCALE", 1)) {
417 param.m_IDparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Down;
418 param.m_IDparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Down;
419
420 param.m_MEparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Down;
421 param.m_MEparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Down;
422
423 param.m_CBparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Down;
424 param.m_CBparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Down;
425
426 } else if (syst == CP::SystematicVariation("MUON_SCALE", -1)) {
427 param.m_IDparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Up;
428 param.m_IDparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Up;
429
430 param.m_MEparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Up;
431 param.m_MEparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Up;
432
433 param.m_CBparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Up;
434 param.m_CBparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Up;
435 } else if (!syst.empty()) return StatusCode::FAILURE;
436
437 // Split scale ID/MS/EGloss
438 syst = systConfig.getSystematicByBaseName("MUON_SCALE_ID");
439
440 if (syst == CP::SystematicVariation("MUON_SCALE_ID", 1)) {
441 param.m_IDparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Down;
442 param.m_IDparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Down;
443 } else if (syst == CP::SystematicVariation("MUON_SCALE_ID", -1)) {
444 param.m_IDparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Up;
445 param.m_IDparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Up;
446 } else if (!syst.empty()) return StatusCode::FAILURE;
447
448 syst = systConfig.getSystematicByBaseName("MUON_SCALE_MS");
449
450 if (syst == CP::SystematicVariation("MUON_SCALE_MS", 1)) {
451 param.m_MEparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Down;
452 } else if (syst == CP::SystematicVariation("MUON_SCALE_MS", -1)) {
453 param.m_MEparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Up;
454 } else if (!syst.empty()) return StatusCode::FAILURE;
455
456 syst = systConfig.getSystematicByBaseName("MUON_SCALE_MS_ELOSS");
457
458 if (syst == CP::SystematicVariation("MUON_SCALE_MS_ELOSS", 1)) {
459 param.m_MEparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Down;
460 } else if (syst == CP::SystematicVariation("MUON_SCALE_MS_ELOSS", -1)) {
461 param.m_MEparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Up;
462 } else if (!syst.empty()) return StatusCode::FAILURE;
463
464 syst = systConfig.getSystematicByBaseName("MUON_SCALE_CB");
465
466 if (syst == CP::SystematicVariation("MUON_SCALE_CB", 1)) {
467 param.m_CBparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Down;
468 } else if (syst == CP::SystematicVariation("MUON_SCALE_CB", -1)) {
469 param.m_CBparams[MCP::ScaleSmearParam::s1] = MCP::SystVariation::Up;
470 } else if (!syst.empty()) return StatusCode::FAILURE;
471
472 syst = systConfig.getSystematicByBaseName("MUON_SCALE_CB_ELOSS");
473
474 if (syst == CP::SystematicVariation("MUON_SCALE_CB_ELOSS", 1)) {
475 param.m_CBparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Down;
476 } else if (syst == CP::SystematicVariation("MUON_SCALE_CB_ELOSS", -1)) {
477 param.m_CBparams[MCP::ScaleSmearParam::s0] = MCP::SystVariation::Up;
478 } else if (!syst.empty()) return StatusCode::FAILURE;
479
480 return StatusCode::SUCCESS;
481 }
482
483 StatusCode MuonCalibIntScaleSmearTool::applySystematicVariation(const CP::SystematicSet &systConfig)
484 {
485 return m_Parameters.get(systConfig, m_currentParameters);
486 }
487
488
489 double MuonCalibIntScaleSmearTool::getExpectedResolution(const int &DetType, double pT, double eta, double phi, MCP::DataYear year, bool addMCCorrectionSmearing) const
490 {
491
492 const auto & IDcorrConstants = m_IDparams.at(year);
493 const auto & MEcorrConstants = m_MEparams.at(year);
494 const auto & IDExpectedResConstants = m_IDExpectedResparams.at(year);
495 const auto & MEExpectedResConstants = m_MEExpectedResparams.at(year);
496
497 MCP::TrackType type = MCP::TrackType::ID;
498 if(DetType == MCP::DetectorType::MS) type = MCP::TrackType::ME;
499
500 auto trk = MCP::TrackCalibObj(type, pT * GeVtoMeV, eta, phi, year, false);
501
502
503 double expRes = 0.;
504 if (DetType == MCP::DetectorType::MS)
505 {
506 if (pT == 0) return 1e12;
507
508 double p0 = 0;
509 double p1 = 0;
510 double p2 = 0;
511
512 if(!addMCCorrectionSmearing)
513 {
514 p0 = MEExpectedResConstants.at(MCP::ExpectedResParam::r0)->getCalibConstant(trk);
515 p1 = MEExpectedResConstants.at(MCP::ExpectedResParam::r1)->getCalibConstant(trk);
516 p2 = MEExpectedResConstants.at(MCP::ExpectedResParam::r2)->getCalibConstant(trk);
517 }
518 else
519 {
520 double expectedP0 = MEExpectedResConstants.at(MCP::ExpectedResParam::r0)->getCalibConstant(trk);
521 double expectedP1 = MEExpectedResConstants.at(MCP::ExpectedResParam::r1)->getCalibConstant(trk);
522 double expectedP2 = MEExpectedResConstants.at(MCP::ExpectedResParam::r2)->getCalibConstant(trk);
523
524 double r0 = MEcorrConstants.at(MCP::ScaleSmearParam::r0).at(MCP::ScaleResCorrection::Nominal)->getCalibConstant(trk);
525 double r1 = MEcorrConstants.at(MCP::ScaleSmearParam::r1).at(MCP::ScaleResCorrection::Nominal)->getCalibConstant(trk);
526 double r2 = MEcorrConstants.at(MCP::ScaleSmearParam::r2).at(MCP::ScaleResCorrection::Nominal)->getCalibConstant(trk);
527
528 p0 = std::sqrt(std::pow(expectedP0, 2) + std::pow(r0, 2));
529 p1 = std::sqrt(std::pow(expectedP1, 2) + std::pow(r1, 2));
530 p2 = std::sqrt(std::pow(expectedP2, 2) + std::pow(r2, 2));
531 }
532
533 expRes = std::sqrt(std::pow(p0 / pT, 2) + std::pow(p1, 2) + std::pow(p2 * pT, 2));
534
535 }
536 else if (DetType == MCP::DetectorType::ID)
537 {
538 if (pT == 0) return 1e12;
539
540 double p1 = 0;
541 double p2 = 0;
542
543 if(!addMCCorrectionSmearing)
544 {
545 p1 = IDExpectedResConstants.at(MCP::ExpectedResParam::r1)->getCalibConstant(trk);
546 p2 = IDExpectedResConstants.at(MCP::ExpectedResParam::r2)->getCalibConstant(trk);
547 double p2Tan = IDExpectedResConstants.at(MCP::ExpectedResParam::r2tan2)->getCalibConstant(trk);
548 if(p2Tan) p2 = p2Tan;
549 }
550 else
551 {
552 double expectedP1 = IDExpectedResConstants.at(MCP::ExpectedResParam::r1)->getCalibConstant(trk);
553 double expectedP2 = IDExpectedResConstants.at(MCP::ExpectedResParam::r2)->getCalibConstant(trk);
554 double p2Tan = IDExpectedResConstants.at(MCP::ExpectedResParam::r2tan2)->getCalibConstant(trk);
555 if(p2Tan) expectedP2 = p2Tan;
556
557 double r1 = IDcorrConstants.at(MCP::ScaleSmearParam::r1).at(MCP::ScaleResCorrection::Nominal)->getCalibConstant(trk);
558 double r2 = IDcorrConstants.at(MCP::ScaleSmearParam::r2).at(MCP::ScaleResCorrection::Nominal)->getCalibConstant(trk);
559
560 p1 = std::sqrt(std::pow(expectedP1, 2) + std::pow(r1, 2));
561 p2 = std::sqrt(std::pow(expectedP2, 2) + std::pow(r2, 2));
562
563 if(p2Tan) p2 = p2 * std::sinh(eta) * std::sinh(eta);
564 }
565
566 expRes = std::sqrt(std::pow(p1, 2) + std::pow(p2 * pT, 2));
567 return expRes;
568 }
569 else
570 {
571 ATH_MSG_ERROR("wrong DetType in input " << DetType);
572 return 0.;
573 }
574
575 return expRes;
576 }
577
578
579} // namespace ADD_CMODE(CP)
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x,...)
Definition AthMsgStreamMacros.h:47
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x,...)
Definition AthMsgStreamMacros.h:42
ATH_MSG_INFO
#define ATH_MSG_INFO(x,...)
Definition AthMsgStreamMacros.h:45
ReadHandle.h
Handle class for reading from StoreGate.
CP::CorrectionCode
Return value from object correction CP tools.
Definition CorrectionCode.h:31
CP::CorrectionCode::Ok
@ Ok
The correction was done successfully.
Definition CorrectionCode.h:38
CP::SystematicRegistry
This module implements the central registry for handling systematic uncertainties with CP tools.
Definition SystematicRegistry.h:25
CP::SystematicRegistry::getInstance
static SystematicRegistry & getInstance()
Get the singleton instance of the registry for the curren thread.
Definition SystematicRegistry.cxx:25
CP::SystematicSet
Class to wrap a set of SystematicVariations.
Definition SystematicSet.h:31
CP::SystematicSet::getSystematicByBaseName
SystematicVariation getSystematicByBaseName(const std::string &basename) const
description: get the first systematic matching basename
Definition SystematicSet.cxx:170
CP::SystematicVariation::empty
bool empty() const
returns: whether this is an empty systematic, i.e.
Definition SystematicVariation.cxx:291
asg::AsgTool
Base class for the dual-use tool implementation classes.
Definition AsgTool.h:47
MCP::CalibInitializer::createExpectedResMap
std::map< ExpectedResParam, std::shared_ptr< CalibContainer > > createExpectedResMap(DataYear dataYear, TrackType type, const std::string &recommendationPath)
Definition CalibInitializer.cxx:74
MCP::CalibInitializer::createScaleResCorrMap
std::map< ScaleResCorrection, std::shared_ptr< CalibContainer > > createScaleResCorrMap(DataYear dataYear, TrackType type, const std::string &recommendationPath, ScaleSmearParam param)
Definition CalibInitializer.cxx:41
MCP::dataYearList
static constexpr std::array< MCP::DataYear, 8 > dataYearList
Definition EnumDef.h:34
MCP::toString
std::string toString(TrackType trkType)
Definition EnumDef.h:39
MCP::DataYear
DataYear
Definition EnumDef.h:28
MCP::TrackType
TrackType
Definition EnumDef.h:13
type
MCP::TrackCalibObj
Basic object to cache all relevant information from the track.
Definition MuonObj.h:74
MCP::TrackCalibObj::eta
const double eta
Value of the track-eta.
Definition MuonObj.h:157
MCP::TrackCalibObj::calib_pt
double calib_pt
Smeared track pt.
Definition MuonObj.h:155
MCP::TrackCalibObj::type
const TrackType type
Flag telling the code whether this is CB/ME/ID.
Definition MuonObj.h:149
Generated on for ATLAS Offline Software by doxygen 1.17.0