ATLAS Offline Software
Loading...
Searching...
No Matches
Functions | Variables
MCAST_Tester.cxx File Reference
#include <cstdlib>
#include <map>
#include <memory>
#include <string>
#include <unordered_map>
#include <TChain.h>
#include <TError.h>
#include <TFile.h>
#include <TTree.h>
#include "xAODRootAccess/Init.h"
#include "xAODRootAccess/TEvent.h"
#include "xAODRootAccess/TStore.h"
#include "AsgMessaging/Check.h"
#include "AsgMessaging/MessageCheck.h"
#include "AsgMessaging/StatusCode.h"
#include "AsgTools/StandaloneToolHandle.h"
#include "MuonAnalysisInterfaces/IMuonCalibrationAndSmearingTool.h"
#include "MuonAnalysisInterfaces/IMuonSelectionTool.h"
#include "PATInterfaces/SystematicRegistry.h"
#include "PATInterfaces/SystematicVariation.h"
#include "xAODCore/ShallowAuxContainer.h"
#include "xAODCore/ShallowCopy.h"
#include "xAODCore/tools/IOStats.h"
#include "xAODCore/tools/ReadStats.h"
#include "xAODEventInfo/EventInfo.h"
#include "xAODMuon/MuonContainer.h"
#include "AthContainers/ConstAccessor.h"
#include <MuonMomentumCorrections/CalibContainer.h>
#include <MuonMomentumCorrections/CalibInitializer.h>
#include "CxxUtils/checker_macros.h"

Go to the source code of this file.

Functions

int main (int argc, char *argv[])

Variables

 ATLAS_NO_CHECK_FILE_THREAD_SAFETY

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 50 of file MCAST_Tester.cxx.

50 {
51
52 const char* APP_NAME = argv[0];
53
54 // setup for ANA_CHECK()
55 using namespace msgMMC;
56 ANA_CHECK_SET_TYPE(int);
57
58 bool useCorrectedCopy = true;
59
60 // The application's name:
61
62 // Check if we received a file name:
63 if (argc < 2) {
64 Error(APP_NAME, "==============================================");
65 Error(APP_NAME, "No file name received!");
66 Error(APP_NAME, "Usage: $> %s [xAOD file name]", APP_NAME);
67 Error(APP_NAME, " $> %s [xAOD file name]", APP_NAME);
68 Error(APP_NAME, " $> %s [xAOD file name] -n X | X = number of events you want to run on", APP_NAME);
69 Error(APP_NAME, " $> %s [xAOD file name] -event X | X = specific number of the event to run on - for debugging", APP_NAME);
70 Error(APP_NAME, "==============================================");
71 return 1;
72 }
73
75 //::: parse the options
77 std::string options;
78 for (int i = 0; i < argc; i++) { options += (argv[i]); }
79
80 int Ievent = -1;
81
82 if (options.find("-event") != std::string::npos) {
83 for (int ipos = 0; ipos < argc; ipos++) {
84 if (std::string(argv[ipos]).compare("-event") == 0) {
85 Ievent = atoi(argv[ipos + 1]);
86 Info(APP_NAME, "Argument (-event) : Running only on event # %i", Ievent);
87 break;
88 }
89 }
90 }
91
92 int nEvents = -1;
93 if (options.find("-n") != std::string::npos) {
94 for (int ipos = 0; ipos < argc; ipos++) {
95 if (std::string(argv[ipos]).compare("-n") == 0) {
96 nEvents = atoi(argv[ipos + 1]);
97 Info(APP_NAME, "Argument (-n) : Running on NEvents = %i", nEvents);
98 break;
99 }
100 }
101 }
102
103 TString limitSys = "";
104 if (options.find("-s") != std::string::npos) {
105 for (int ipos = 0; ipos < argc; ipos++) {
106 if (std::string(argv[ipos]).compare("-s") == 0) {
107 limitSys = TString(argv[ipos + 1]);
108 break;
109 }
110 }
111 }
112
113 bool doSys = false;
114 if (options.find("-doSys") != std::string::npos) {
115 doSys = true;
116 }
118 //::: initialize the application and get the event
120 ANA_CHECK(xAOD::Init(APP_NAME));
121
122 //::: Open the input file:
123 std::string fileName = argv[1];
124 Info(APP_NAME, "Opening file: %s", fileName.c_str());
125 std::unique_ptr<TFile> ifile(TFile::Open(fileName.c_str(), "READ"));
126 if (!ifile) Error(APP_NAME, "Cannot find file %s", fileName.c_str());
127
128 std::unique_ptr<TChain> chain = std::make_unique<TChain>("CollectionTree", "CollectionTree");
129 chain->Add(fileName.c_str());
130
131 //::: Create a TEvent object:
132 xAOD::TEvent event((TTree*)chain.get(), xAOD::TEvent::kAthenaAccess);
133 Info(APP_NAME, "Number of events in the file: %i", static_cast<int>(event.getEntries()));
134
135 //::: Create a transient object store. Needed for the tools.
136 xAOD::TStore store;
137
138 //::: Decide how many events to run over:
139 Long64_t entries = event.getEntries();
140 if (fileName.find("data") != std::string::npos) entries = 2000;
141 // if (fileName.find("mc20") != std::string::npos) entries = 2000;
142
143 if(doSys) entries = 200;
144
146 //::: MuonCalibrationAndSmearingTool
147 // setup the tool handle as per the
148 // recommendation by ASG - https://twiki.cern.ch/twiki/bin/view/AtlasProtected/AthAnalysisBase#How_to_use_AnaToolHandle
150
151 // bool is run3 geo
152 bool isRun3Geo = false;
153 if(fileName.find("mc21") != std::string::npos) isRun3Geo = true;
154 if(fileName.find("data22") != std::string::npos) isRun3Geo = true;
155
156 asg::StandaloneToolHandle<CP::IMuonCalibrationAndSmearingTool> corrTool;
157 corrTool.setProperty("calibMode", 1).ignore();
158 corrTool.setProperty("IsRun3Geo", isRun3Geo).ignore();
159 // corrTool.setProperty("doExtraSmearing", true).ignore();
160 // corrTool.setProperty("do2StationsHighPt", false).ignore();
161
162 corrTool.setTypeAndName("CP::MuonCalibTool/NewMCT");
163
164
165
166 bool isDebug = false;
167 if (nEvents >= 0 || Ievent >= 0) isDebug = true;
168
169 if (isDebug) corrTool.setProperty("OutputLevel", MSG::VERBOSE).ignore();
170 //::: retrieve the tool
171 StatusCode sc = corrTool.retrieve();
172 if (sc.isFailure()) {
173 Error(APP_NAME, "Cannot retrieve MuonCorrectionTool");
174 return 1;
175 }
176
177 //::: retrieve the tool
178 sc = corrTool.retrieve();
179 if (sc.isFailure()) {
180 Error(APP_NAME, "Cannot retrieve MuonCorrectionTool");
181 return 1;
182 }
183
184
186 //::: MuonSelectionTool
187 // setup the tool handle as per the
188 // recommendation by ASG - https://twiki.cern.ch/twiki/bin/view/AtlasProtected/AthAnalysisBase#How_to_use_AnaToolHandle
190 //::: create the tool handle
191 asg::StandaloneToolHandle<CP::IMuonSelectionTool> selTool;
192 selTool.setTypeAndName("CP::MuonSelectionTool/MuonSelectionTool");
193
194 //::: set the properties
195 selTool.setProperty("IsRun3Geo", isRun3Geo).ignore();
196 selTool.setProperty("MaxEta", 2.5).ignore();
197 selTool.setProperty("MuQuality", (int)xAOD::Muon::Loose).ignore(); // corresponds to 0=Tight, 1=Medium, 2=Loose, 3=VeryLoose, 4=HighPt, 5=LowPtEfficiency
198
199 //::: retrieve the tool
200 if (selTool.retrieve().isFailure() || sc.isFailure()) {
201 Error(APP_NAME, "Cannot retrieve MuonSelectionTool");
202 return 1;
203 }
204
206 //::: Systematics initialization
208 std::vector<CP::SystematicSet> sysList;
209 if(limitSys.Length() == 0)sysList.push_back(CP::SystematicSet());
210
211 if(doSys)
212 {
213 const CP::SystematicRegistry& registry = CP::SystematicRegistry::getInstance();
214 const CP::SystematicSet& recommendedSystematics = registry.recommendedSystematics();
215 for (CP::SystematicSet::const_iterator sysItr = recommendedSystematics.begin(); sysItr != recommendedSystematics.end(); ++sysItr) {
216 if(!TString(sysItr->name()).Contains(limitSys) && limitSys.Length() > 0) continue;
217 sysList.push_back(CP::SystematicSet());
218 sysList.back().insert(*sysItr);
219 }
220 }
221
222 std::vector<CP::SystematicSet>::const_iterator sysListItr;
223
224 msgMMC::ANA_MSG_INFO("main() - Systematics are ");
225 for (sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr) msgMMC::ANA_MSG_INFO(sysListItr->name());
226
227 // branches to be stored
228 Float_t InitPtCB(0.), InitPtID(0.), InitPtMS(0.);
229 Float_t CorrPtCB(0.), CorrPtID(0.), CorrPtMS(0.);
230 Float_t Eta(0.), Phi(0.), Charge(0.);
231 Float_t ExpResoCB(0.), ExpResoID(0.), ExpResoMS(0.);
232 long long unsigned int eventNum(0);
233
234 // output file
235 TFile* outputFile = TFile::Open("output.root", "recreate");
236
237 // it contains a single TTree per systematic variation
238 std::unordered_map<CP::SystematicSet, TTree*> sysTreeMap;
239 for (sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr) {
240 // create new tree for the systematic in question
241 std::string treeName = "test_tree_" + (sysListItr->name().size() == 0 ? "NOMINAL" : sysListItr->name());
242 TTree* sysTree = new TTree(treeName.c_str(), "test tree for MCAST");
243
244 // add branches
245 sysTree->Branch("InitPtCB", &InitPtCB, "InitPtCB/F");
246 sysTree->Branch("InitPtID", &InitPtID, "InitPtID/F");
247 sysTree->Branch("InitPtMS", &InitPtMS, "InitPtMS/F");
248 sysTree->Branch("CorrPtCB", &CorrPtCB, "CorrPtCB/F");
249 sysTree->Branch("CorrPtID", &CorrPtID, "CorrPtID/F");
250 sysTree->Branch("CorrPtMS", &CorrPtMS, "CorrPtMS/F");
251 sysTree->Branch("Eta", &Eta, "Eta/F");
252 sysTree->Branch("Phi", &Phi, "Phi/F");
253 sysTree->Branch("Charge", &Charge, "Charge/F");
254 sysTree->Branch("ExpResoCB", &ExpResoCB, "ExpResoCB/F");
255 sysTree->Branch("ExpResoID", &ExpResoID, "ExpResoID/F");
256 sysTree->Branch("ExpResoMS", &ExpResoMS, "ExpResoMS/F");
257 sysTree->Branch("eventNum", &eventNum);
258
259 sysTreeMap[*sysListItr] = sysTree;
260 }
261
263 //::: Loop over the events
265 for (Long64_t entry = 0; entry < entries; ++entry) {
266 if (nEvents != -1 && entry > nEvents) break;
267 // Tell the object which entry to look at:
268 event.getEntry(entry);
269
270 // Print some event information
271 const xAOD::EventInfo* evtInfo = 0;
272 ANA_CHECK(event.retrieve(evtInfo, "EventInfo"));
273 if (Ievent != -1 && static_cast<int>(evtInfo->eventNumber()) != Ievent) { continue; }
274
275 eventNum = evtInfo->eventNumber();
276
277 //::: Get the Muons from the event:
278 const xAOD::MuonContainer* muons = 0;
279 ANA_CHECK(event.retrieve(muons, "Muons"));
280
281 //::: Loop over systematics
282 for (sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr) {
283 // create a shallow copy of the muons container
284 std::pair<xAOD::MuonContainer*, xAOD::ShallowAuxContainer*> muons_shallowCopy = xAOD::shallowCopyContainer(*muons);
285
286 xAOD::MuonContainer* muonsCorr = muons_shallowCopy.first;
287
288 if (isDebug) {
289 Info(APP_NAME, "-----------------------------------------------------------");
290 Info(APP_NAME, "Looking at %s systematic", (sysListItr->name()).c_str());
291 }
292 //::: Check if systematic is applicable to the correction tool
293 if (corrTool->applySystematicVariation(*sysListItr) != StatusCode::SUCCESS) {
294 Error(APP_NAME, "Cannot configure muon calibration tool for systematic");
295 }
296
297 //::: Loop over muon container
298 for (auto muon : *muonsCorr) {
299 //::: Select "good" muons:
300 if (!selTool->accept(*muon)) {
301 if (isDebug) Info(APP_NAME, "This muon doesn't pass the ID hits quality cuts");
302 continue;
303 }
304
305 //::: Should be using correctedCopy here, testing behaviour of applyCorrection though
306 InitPtCB = muon->pt();
307 InitPtID = -999;
308 if (muon->inDetTrackParticleLink().isValid()) {
309 const ElementLink<xAOD::TrackParticleContainer>& id_track = muon->inDetTrackParticleLink();
310 InitPtID = (!id_track) ? 0 : (*id_track)->pt();
311 }
312 InitPtMS = -999;
313 if (muon->extrapolatedMuonSpectrometerTrackParticleLink().isValid()) {
314 const ElementLink<xAOD::TrackParticleContainer>& ms_track = muon->extrapolatedMuonSpectrometerTrackParticleLink();
315 InitPtMS = (!ms_track) ? 0 : (*ms_track)->pt();
316 }
317
318 Eta = muon->eta();
319 Phi = muon->phi();
320 Charge = muon->charge();
321
322 //::: Print some info about the selected muon:
323 if (isDebug) Info(APP_NAME, "Selected muon: eta = %g, phi = %g, pt = %g", muon->eta(), muon->phi(), muon->pt() / 1e3);
324
325 float ptCB = 0;
326 if (muon->primaryTrackParticleLink().isValid()) {
327 const ElementLink<xAOD::TrackParticleContainer>& cb_track = muon->primaryTrackParticleLink();
328 ptCB = (!cb_track) ? 0 : (*cb_track)->pt();
329 } else {
330 if (isDebug)
331 Info(APP_NAME, "Missing primary track particle link for --> CB %g, author: %d, type: %d", ptCB, muon->author(),
332 muon->muonType());
333 }
334 float ptID = 0;
335 if (muon->inDetTrackParticleLink().isValid()) {
336 const ElementLink<xAOD::TrackParticleContainer>& id_track = muon->inDetTrackParticleLink();
337 ptID = (!id_track) ? 0 : (*id_track)->pt();
338 }
339 float ptME = 0;
340 if (muon->extrapolatedMuonSpectrometerTrackParticleLink().isValid()) {
341 const ElementLink<xAOD::TrackParticleContainer>& ms_track = muon->extrapolatedMuonSpectrometerTrackParticleLink();
342 ptME = (!ms_track) ? 0 : (*ms_track)->pt();
343 }
344
345 if (isDebug)
346 Info(APP_NAME, "--> CB %g, ID %g, ME %g, author: %d, type: %d", ptCB / 1e3, ptID / 1e3, ptME / 1e3, muon->author(),
347 muon->muonType());
348
349 // either use the correctedCopy call or correct the muon object itself
350 static const SG::ConstAccessor<float> InnerDetectorPtAcc ("InnerDetectorPt");
351 static const SG::ConstAccessor<float> MuonSpectrometerPtAcc ("MuonSpectrometerPt");
352 if (useCorrectedCopy) {
353 // ::: Create a calibrated muon:
354 xAOD::Muon* mu = 0;
355 if (!corrTool->correctedCopy(*muon, mu)) {
356 Error(APP_NAME, "Cannot really apply calibration nor smearing");
357 continue;
358 }
359 CorrPtCB = mu->pt();
360 CorrPtID = InnerDetectorPtAcc (*mu);
361 CorrPtMS = MuonSpectrometerPtAcc (*mu);
362
363
364 sysTreeMap[*sysListItr]->Fill();
365 //::: Delete the calibrated muon:
366 delete mu;
367 }
368 else {
369 if (!corrTool->applyCorrection(*muon)) {
370 Error(APP_NAME, "Cannot really apply calibration nor smearing");
371 continue;
372 }
373 CorrPtCB = muon->pt();
374 CorrPtID = InnerDetectorPtAcc (*muon);
375 CorrPtMS = MuonSpectrometerPtAcc (*muon);
376 ExpResoCB = corrTool->expectedResolution("CB", *muon, true);
377 ExpResoID = corrTool->expectedResolution("ID", *muon, true);
378 ExpResoMS = corrTool->expectedResolution("MS", *muon, true);
379 if (isDebug)
380 Info(APP_NAME, "Calibrated muon: eta = %g, phi = %g, pt(CB) = %g, pt(ID) = %g, pt(MS) = %g", muon->eta(),
381 muon->phi(), muon->pt() / 1e3, CorrPtID / 1e3, CorrPtMS / 1e3);
382 if (isDebug)
383 Info(APP_NAME, " expReso : ExpResoCB = %g , ExpResoID = %g , ExpResoMS = %g", ExpResoCB, ExpResoID, ExpResoMS);
384 sysTreeMap[*sysListItr]->Fill();
385 }
386 }
387 if (isDebug) Info(APP_NAME, "-----------------------------------------------------------");
388
389 delete muons_shallowCopy.first;
390 delete muons_shallowCopy.second;
391 }
392 //::: Close with a message:
393 if (entry % 5 == 0)
394 Info(APP_NAME,
395 "===>>> done processing event #%i, run #%i %i events processed so far <<<===", static_cast<int>(evtInfo->eventNumber()),
396 static_cast<int>(evtInfo->runNumber()), static_cast<int>(entry + 1));
397 }
398
399 for (sysListItr = sysList.begin(); sysListItr != sysList.end(); ++sysListItr) { sysTreeMap[*sysListItr]->Write(); }
400
401 //::: Close output file
402 outputFile->Close();
403
404 xAOD::IOStats::instance().stats().printSmartSlimmingBranchList();
405
406 //::: Return gracefully:
407 return 0;
408}
APP_NAME
#define APP_NAME
Definition BoostedXbbTag.cxx:23
ANA_CHECK
#define ANA_CHECK(EXP)
check whether the given expression was successful
Definition Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:324
ANA_CHECK_SET_TYPE
#define ANA_CHECK_SET_TYPE(TYPE)
set the type for ANA_CHECK to report failures
Definition Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:314
sc
static Double_t sc
Definition LArPhysWaveHECTool.cxx:37
Phi
@ Phi
Definition RPCdef.h:8
Eta
@ Eta
Definition RPCdef.h:8
CP::SystematicRegistry
This module implements the central registry for handling systematic uncertainties with CP tools.
Definition SystematicRegistry.h:25
CP::SystematicRegistry::recommendedSystematics
const SystematicSet & recommendedSystematics() const
returns: the recommended set of systematics
Definition SystematicRegistry.cxx:38
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::end
const_iterator end() const
description: const iterator to the end of the set
Definition SystematicSet.h:59
CP::SystematicSet::const_iterator
std::set< SystematicVariation >::const_iterator const_iterator
Definition SystematicSet.h:52
CP::SystematicSet::begin
const_iterator begin() const
description: const iterator to the beginning of the set
Definition SystematicSet.h:55
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition ConstAccessor.h:55
asg::StandaloneToolHandle
an "initializing" ToolHandle for stand-alone applications
Definition StandaloneToolHandle.h:44
asg::StandaloneToolHandle::setProperty
StatusCode setProperty(const std::string &name, T2 &&value)
Definition StandaloneToolHandle.h:105
asg::StandaloneToolHandle::retrieve
StatusCode retrieve()
initialize the tool, will succeed if the tool was already initialized
Definition StandaloneToolHandle.h:147
asg::StandaloneToolHandle::setTypeAndName
void setTypeAndName(const std::string &typeAndName)
Definition StandaloneToolHandle.h:101
xAOD::EventInfo_v1::runNumber
uint32_t runNumber() const
The current event's run number.
xAOD::EventInfo_v1::eventNumber
uint64_t eventNumber() const
The current event's event number.
xAOD::IOStats::stats
ReadStats & stats()
Access the object belonging to the current thread.
Definition IOStats.cxx:17
xAOD::IOStats::instance
static IOStats & instance()
Singleton object accessor.
Definition IOStats.cxx:11
xAOD::ReadStats::printSmartSlimmingBranchList
void printSmartSlimmingBranchList(bool autoIncludeLinks=false) const
Print the accessed variables, formatted for smart slimming.
xAOD::TEvent
Tool for accessing xAOD files outside of Athena.
Definition Control/xAODRootAccess/xAODRootAccess/TEvent.h:59
xAOD::TEvent::kAthenaAccess
@ kAthenaAccess
Access containers/objects like Athena does.
Definition Control/xAODRootAccess/xAODRootAccess/TEvent.h:74
xAOD::TStore
A relatively simple transient store for objects created in analysis.
Definition TStore.h:45
nEvents
const int nEvents
Definition fbtTestBasics.cxx:78
entries
double entries
Definition listroot.cxx:49
Amg::compare
std::pair< int, int > compare(const AmgSymMatrix(N) &m1, const AmgSymMatrix(N) &m2, double precision=1e-9, bool relative=false)
compare two matrices, returns the indices of the first element that fails the condition,...
Definition EventPrimitivesHelpers.h:109
CaloNoise_fillDB.mu
mu
Definition CaloNoise_fillDB.py:51
CxxUtils::atoi
int atoi(std::string_view str)
Helper functions to unpack numbers decoded in string into integers and doubles The strings are requir...
Definition Control/CxxUtils/Root/StringUtils.cxx:85
L1Topo::Error
Error
The different types of error that can be flagged in the L1TopoRDO.
Definition Error.h:16
SGTest::store
TestStore store
Definition TestStore.cxx:23
ZDCMsg::Info
@ Info
Definition ZDCMsg.h:20
test_pyathena.pt
pt
Definition test_pyathena.py:11
xAOD::Init
StatusCode Init(const char *appname)
Function initialising ROOT/PyROOT for using the ATLAS EDM.
Definition Init.cxx:31
xAOD::EventInfo
EventInfo_v1 EventInfo
Definition of the latest event info version.
Definition IEventInfoCnvTool.h:18
xAOD::shallowCopyContainer
std::pair< std::unique_ptr< T >, std::unique_ptr< ShallowAuxContainer > > shallowCopyContainer(const T &cont, const EventContext &ctx)
Function making a shallow copy of a constant container.
Definition ShallowCopy.h:110
xAOD::Muon
Muon_v1 Muon
Reference the current persistent version:
Definition Event/xAOD/xAODMuon/xAODMuon/Muon.h:13
xAOD::MuonContainer
MuonContainer_v1 MuonContainer
Definition of the current "Muon container version".
Definition Event/xAOD/xAODMuon/xAODMuon/MuonContainer.h:14

Variable Documentation

◆ ATLAS_NO_CHECK_FILE_THREAD_SAFETY

ATLAS_NO_CHECK_FILE_THREAD_SAFETY

Definition at line 47 of file MCAST_Tester.cxx.

Generated on for ATLAS Offline Software by doxygen 1.14.0