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