ATLAS Offline Software
xAODHist.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 # Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 
4 import sys
5 import os
6 import ROOT
7 import argparse
8 
9 
10 def safeRetrieve(evt, typ, key):
11  if evt.contains(typ, key):
12  return evt.retrieve(typ, key)
13  print(f'WARNING: Cannot find object {typ}/{key}')
14  return []
15 
16 
17 def xAODHist(evt, phys=False, analysis=False, histfile=None):
18  # Set output HIST filename
19  if histfile:
20  histfilename = histfile
21  else:
22  histfilename = "hist.root"
23 
24  # Book histograms
25  hfile = ROOT.TFile( histfilename, 'RECREATE', 'ROOT file with histograms' )
26 
27  if phys:
28  hphys_nclus = ROOT.TH1F( 'phys_nclus', 'phys_nclus', 20, 0, 20 )
29  hphys_nidtracks = ROOT.TH1F( 'phys_nidtracks', 'phys_nidtracks', 20, 0, 20 )
30  hphys_ntautracks = ROOT.TH1F( 'phys_ntautracks', 'phys_ntautracks', 20, 0, 20 )
31  hphys_ntaus = ROOT.TH1F( 'phys_ntaus', 'phys_taus', 20, 0, 20 )
32  hphys_nmuons = ROOT.TH1F( 'phys_nmuons', 'phys_nmuons', 20, 0, 20 )
33  hphys_nelecs = ROOT.TH1F( 'phys_nelecs', 'phys_nelecs', 20, 0, 20 )
34  hphys_nphotons = ROOT.TH1F( 'phys_nphotons', 'phys_nphotons', 20, 0, 20 )
35  hphys_njets = ROOT.TH1F( 'phys_njets', 'phys_njets', 20, 0, 20 )
36  hphys_nfakeelectrons = ROOT.TH1F( 'phys_nfakeelectrons', 'phys_nfakeelectrons', 20, 0, 20 )
37  hphys_nfakephotons = ROOT.TH1F( 'phys_nfakephotons', 'phys_nfakephotons', 20, 0, 20 )
38 
39  if analysis:
40  hntaujet = ROOT.TH1F( 'ana_ntaujet', 'ana_ntaujet', 20, 0, 20 )
41  htaujetpt = ROOT.TH1F( 'ana_taujet_pt', 'ana_taujet_pt', 200, 0, 200 )
42  htaujeteta = ROOT.TH1F( 'ana_taujet_eta', 'ana_taujet_eta', 100, -5, 5 )
43  htaujetphi = ROOT.TH1F( 'ana_taujet_phi', 'ana_taujet_phi', 70, -3.5, 3.5 )
44 
45  hnelectron = ROOT.TH1F( 'ana_nelectron', 'ana_nelectron', 20, 0, 20 )
46  helectronpt = ROOT.TH1F( 'ana_electron_pt', 'ana_electron_pt', 200, 0, 200 )
47  helectroneta = ROOT.TH1F( 'ana_electron_eta', 'ana_electron_eta', 100, -5, 5 )
48  helectronphi = ROOT.TH1F( 'ana_electron_phi', 'ana_electron_phi', 70, -3.5, 3.5 )
49 
50  hnmuon = ROOT.TH1F( 'ana_nmuon', 'ana_nmuon', 20, 0, 20 )
51  hmuonpt = ROOT.TH1F( 'ana_muon_pt', 'ana_muon_pt', 200, 0, 200 )
52  hmuoneta = ROOT.TH1F( 'ana_muon_eta', 'ana_muon_eta', 100, -5, 5 )
53  hmuonphi = ROOT.TH1F( 'ana_muon_phi', 'ana_muon_phi', 70, -3.5, 3.5 )
54 
55  hnphoton = ROOT.TH1F( 'ana_nphoton', 'ana_nphoton', 20, 0, 20 )
56  hphotonpt = ROOT.TH1F( 'ana_photon_pt', 'ana_photon_pt', 200, 0, 200 )
57  hphotoneta = ROOT.TH1F( 'ana_photon_eta', 'ana_photon_eta', 100, -5, 5 )
58  hphotonphi = ROOT.TH1F( 'ana_photon_phi', 'ana_photon_phi', 70, -3.5, 3.5 )
59 
60  hnjet = ROOT.TH1F( 'ana_njet', 'ana_njet', 20, 0, 20 )
61  hjetpt = ROOT.TH1F( 'ana_jet_pt', 'ana_jet_pt', 200, 0, 200 )
62  hjeteta = ROOT.TH1F( 'ana_jet_eta', 'ana_jet_eta', 100, -5, 5 )
63  hjetphi = ROOT.TH1F( 'ana_jet_phi', 'ana_jet_phi', 70, -3.5, 3.5 )
64 
65  # Loop over events
66  for i in range(0, evt.getEntries()):
67  evt.getEntry(i)
68 
69  if phys:
70  clusters = safeRetrieve(
71  evt, "xAOD::CaloClusterContainer", "CaloCalTopoClusters")
72  nclus = len(clusters)
73  hphys_nclus.Fill( nclus )
74 
75  idTracks = safeRetrieve(evt,
76  "xAOD::TrackParticleContainer", "InDetTrackParticles")
77  nIdTracks = len(idTracks)
78  hphys_nidtracks.Fill( nIdTracks )
79 
80  tautracks = safeRetrieve(evt, "xAOD::TauTrackContainer", "TauTracks")
81  nTauTracks = len(tautracks)
82  hphys_ntautracks.Fill( nTauTracks )
83 
84  taus = safeRetrieve(evt, "xAOD::TauJetContainer", "TauJets")
85  nTaus = len(taus)
86  hphys_ntaus.Fill( nTaus )
87 
88  muons = safeRetrieve(evt, "xAOD::MuonContainer", "Muons")
89  nMuons = len(muons)
90  hphys_nmuons.Fill( nMuons )
91 
92  electrons = safeRetrieve(evt, "xAOD::ElectronContainer", "Electrons")
93  nElec = len(electrons)
94  hphys_nelecs.Fill( nElec )
95 
96  photons = safeRetrieve(evt, "xAOD::PhotonContainer", "Photons")
97  nPhot = len(photons)
98  hphys_nphotons.Fill( nPhot )
99 
100  jets = safeRetrieve(evt,"xAOD::JetContainer", "AntiKt4EMPFlowJets")
101  nJet = len(jets)
102  hphys_njets.Fill( nJet )
103 
104  nTrueElectrons = 0
105  nTruePhotons = 0
106  acc = ROOT.SG.ConstAccessor(
107  'ElementLink< xAOD::TruthParticleContainer>')('truthParticleLink')
108 
109  if nElec > 0 and acc.isAvailable(electrons.at(0)):
110  for i in range(nElec):
111  truthLink = acc(electrons.at(i))
112  if(truthLink.isValid()):
113  pdgId = truthLink.pdgId()
114  if abs(pdgId) == 11:
115  nTrueElectrons += 1
116 
117  if nPhot > 0 and acc.isAvailable(photons.at(0)):
118  for i in range(nPhot):
119  truthLink = acc(photons.at(i))
120  if(truthLink.isValid()):
121  pdgId = truthLink.pdgId()
122  if (pdgId == 22):
123  nTruePhotons += 1
124 
125  nFakeElectrons = nElec - nTrueElectrons
126  nFakePhotons = nPhot - nTruePhotons
127  hphys_nfakeelectrons.Fill( nFakeElectrons )
128  hphys_nfakephotons.Fill( nFakePhotons )
129 
130  # PHYSLITE input with Analysis containers
131  if analysis:
132  ana_taujets = safeRetrieve(evt,"xAOD::TauJetContainer", "AnalysisTauJets")
133  ana_ntaujet = len(ana_taujets)
134  hntaujet.Fill( ana_ntaujet )
135  for j in ana_taujets:
136  htaujetpt.Fill( j.pt()/1000. )
137  htaujeteta.Fill( j.eta() )
138  htaujetphi.Fill( j.phi() )
139 
140  ana_muons = safeRetrieve(evt,"xAOD::MuonContainer", "AnalysisMuons")
141  ana_nmuon = len(ana_muons)
142  hnmuon.Fill( ana_nmuon )
143  for j in ana_muons:
144  hmuonpt.Fill( j.pt()/1000. )
145  hmuoneta.Fill( j.eta() )
146  hmuonphi.Fill( j.phi() )
147 
148  ana_electrons = safeRetrieve(evt,"xAOD::ElectronContainer", "AnalysisElectrons")
149  ana_nelectron = len(ana_electrons)
150  hnelectron.Fill( ana_nelectron )
151  for j in ana_electrons:
152  helectronpt.Fill( j.pt()/1000. )
153  helectroneta.Fill( j.eta() )
154  helectronphi.Fill( j.phi() )
155 
156  ana_photons = safeRetrieve(evt,"xAOD::PhotonContainer", "AnalysisPhotons")
157  ana_nPhot = len(ana_photons)
158  hnphoton.Fill( ana_nPhot )
159  for j in ana_photons:
160  hphotonpt.Fill( j.pt()/1000. )
161  hphotoneta.Fill( j.eta() )
162  hphotonphi.Fill( j.phi() )
163 
164  ana_jets = safeRetrieve(evt,"xAOD::JetContainer", "AnalysisJets")
165  ana_nJet = len(ana_jets)
166  hnjet.Fill( ana_nJet )
167  for j in ana_jets:
168  hjetpt.Fill( j.pt()/1000. )
169  hjeteta.Fill( j.eta() )
170  hjetphi.Fill( j.phi() )
171 
172  # Store HIST in output file in different directories
173  if phys:
174  hfile.cd()
175  hfile.mkdir("PHYS")
176  hfile.cd("PHYS")
177  hphys_nclus.Write()
178  hphys_nidtracks.Write()
179  hphys_ntautracks.Write()
180  hphys_ntaus.Write()
181  hphys_nmuons.Write()
182  hphys_nelecs.Write()
183  hphys_nphotons.Write()
184  hphys_njets.Write()
185  hphys_nfakeelectrons.Write()
186  hphys_nfakephotons.Write()
187 
188  if analysis:
189  hfile.cd()
190  hfile.mkdir("AnalysisJets")
191  hfile.cd("AnalysisJets")
192  hnjet.Write()
193  hjetpt.Write()
194  hjeteta.Write()
195  hjetphi.Write()
196 
197  hfile.cd("..")
198  hfile.mkdir("AnalysisPhotons")
199  hfile.cd("AnalysisPhotons")
200  hnphoton.Write()
201  hphotonpt.Write()
202  hphotoneta.Write()
203  hphotonphi.Write()
204 
205  hfile.cd("..")
206  hfile.mkdir("AnalysisTauJets")
207  hfile.cd("AnalysisTauJets")
208  hntaujet.Write()
209  htaujetpt.Write()
210  htaujeteta.Write()
211  htaujetphi.Write()
212 
213  hfile.cd("..")
214  hfile.mkdir("AnalysisElectrons")
215  hfile.cd("AnalysisElectrons")
216  hnelectron.Write()
217  helectronpt.Write()
218  helectroneta.Write()
219  helectronphi.Write()
220 
221  hfile.cd("..")
222  hfile.mkdir("AnalysisMuons")
223  hfile.cd("AnalysisMuons")
224  hnmuon.Write()
225  hmuonpt.Write()
226  hmuoneta.Write()
227  hmuonphi.Write()
228 
229  hfile.Close()
230 
231  return
232 
233 def main():
234  parser = argparse.ArgumentParser(
235  description="Extracts a few basic quantities from the xAOD file and dumps them into a hist ROOT file")
236  parser.add_argument("xAODFile", nargs='?', type=str,
237  help="xAOD filename", action="store")
238  parser.add_argument("--phys", help="Create histogram file from DAOD_PHYS or AOD variables",
239  action="store_true", default=False)
240  parser.add_argument("--analysis", help="Create histogram file from DAOD_PHYSLITE n/pt/eta/phi variables",
241  action="store_true", default=False)
242  parser.add_argument("--outputHISTFile", help="histogram output filename",
243  action="store", default=None)
244  parser.add_argument("--inputisESD", help="Set if input is ESD",
245  action="store_true", default=False)
246 
247  args = parser.parse_args()
248 
249  if len(sys.argv) < 2:
250  parser.print_help()
251  sys.exit(1)
252 
253  # Check input file existance
254  if not os.access(args.xAODFile, os.R_OK):
255  print("ERROR, can't access file {}".format(args.xAODFile))
256  sys.exit(1)
257 
258  # Create single inputfile
259  filelist = args.xAODFile
260 
261  # Setup TEvent object and add inputs
262  evt = ROOT.POOL.TEvent(
263  ROOT.POOL.TEvent.kPOOLAccess if args.inputisESD else ROOT.POOL.TEvent.kClassAccess)
264  stat = evt.readFrom(filelist)
265  if not stat:
266  print("ERROR, failed to open file {} with POOL.TEvent".format(
267  args.xAODFile))
268  sys.exit(1)
269  pass
270 
271  xAODHist(evt, args.phys, args.analysis, args.outputHISTFile)
272 
273  return 0
274 
275 if __name__ == "__main__":
276 
277  main()
vtune_athena.format
format
Definition: vtune_athena.py:14
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
AthenaPoolTestRead.acc
acc
Definition: AthenaPoolTestRead.py:16
python.xAODHist.main
def main()
Definition: xAODHist.py:233
print
void print(char *figname, TCanvas *c1)
Definition: TRTCalib_StrawStatusPlots.cxx:26
python.xAODHist.xAODHist
def xAODHist(evt, phys=False, analysis=False, histfile=None)
Definition: xAODHist.py:17
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
python.xAODHist.safeRetrieve
def safeRetrieve(evt, typ, key)
Definition: xAODHist.py:10