ATLAS Offline Software
ICHEP2016/ParsePunchthroughInput.py
Go to the documentation of this file.
1 # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
2 
3 from ROOT import *
4 from array import array
5 import glob
6 import re
7 import math
8 import ProviderHistoHelpers
9 
10 def GetKeyNames(self,dir=""):
11  self.cd(dir)
12  return [key.GetName() for key in gDirectory.GetListOfKeys()]
13 TFile.GetKeyNames = GetKeyNames
14 
15 jetDefDict = {
16  'AntiKt4EMTopo' : 'AntiKt4Topo_EMJES',
17 # 'AntiKt6TopoEM' : 'AntiKt6Topo_EMJES',
18  'AntiKt4LCTopo' : 'AntiKt4Topo_LCJES',
19 # 'AntiKt6LCTopo' : 'AntiKt6Topo_LCJES'
20  }
21 
22 SystematicNameList = ['PunchThrough_MC15','PunchThrough_AFII']
23 
24 def ReadPunchthroughHistograms(dirName):
25  if not dirName.endswith("/"):
26  dirName = dirName + "/"
27 
28  # Run over the one file
29  rootFileList = sorted(glob.glob(dirName+"*.root"))
30  fileDict = {}
31  if len(rootFileList) != 1:
32  print "Found a number of root files not equal to 1 in dir:",dirName
33  print "Assume separate files for EM, LC"
34  for aJetDef in jetDefDict.keys() :
35  for file in rootFileList :
36  if aJetDef in file :
37  if aJetDef not in fileDict.keys() :
38  fileDict[aJetDef] = file
39  else :
40  print "Oh no! More than one file corresponding to this jet definition:",aJetDef
41  return None
42  else :
43  print "Found only one root file: assume contents can be used for all jet collections."
44  for aJetDef in jetDefDict.keys() :
45  fileDict[aJetDef] = rootFileList[0]
46 
47  histos = {}
48  for aJetDef,aJetDefProv in jetDefDict.iteritems():
49  rootFile = TFile(fileDict[aJetDef],"READ")
50  histos[aJetDefProv] = {}
51 
52  for histName in rootFile.GetKeyNames():
53  if aJetDef not in histName: continue
54  for aKey in SystematicNameList:
55  if "_MC" in aKey :
56  histName = aKey+"_"+aJetDef
57  histo = rootFile.Get(histName)
58  if histo is None:
59  print "Failed to get histogram:",histName
60  return None
61  histo.SetName(aKey+"_"+aJetDefProv+"badaxes")
62  histo.SetDirectory(0)
63 
64  # In Lydia's histogram:
65  # x axis is pT.
66  # y axis is eta
67  # z axis is NSegments
68  # Last two axes need to get swapped to be compatible with 8TeV!!!!
69 
70  # In Christopher's histogram,
71  # these are back to being as requested.
72 
73  pTBinArray = []
74  nPTBins = histo.GetXaxis().GetNbins()
75  for bin in range(1,nPTBins+2) :
76  binEdge = histo.GetXaxis().GetBinLowEdge(bin)
77  pTBinArray.append(binEdge)
78  pTBinArray = array('d',pTBinArray)
79 
80  etaBinArray = []
81  if "EM" in aJetDef :
82  nEtaBins = histo.GetYaxis().GetNbins()
83  for bin in range(1,nEtaBins+2) :
84  binEdge = histo.GetYaxis().GetBinLowEdge(bin)
85  etaBinArray.append(binEdge)
86  else :
87  nEtaBins = histo.GetZaxis().GetNbins()
88  for bin in range(1,nEtaBins+2) :
89  binEdge = histo.GetZaxis().GetBinLowEdge(bin)
90  etaBinArray.append(binEdge)
91  etaBinArray = array('d',etaBinArray)
92 
93  nSegmentsBinArray = []
94  if "EM" in aJetDef :
95  nNSegBins = histo.GetZaxis().GetNbins()
96  for bin in range(1,nNSegBins+2) :
97  binEdge = histo.GetZaxis().GetBinLowEdge(bin)
98  nSegmentsBinArray.append(binEdge)
99  else :
100  nNSegBins = histo.GetYaxis().GetNbins()
101  for bin in range(1,nNSegBins+2) :
102  binEdge = histo.GetYaxis().GetBinLowEdge(bin)
103  nSegmentsBinArray.append(binEdge)
104 
105  nSegmentsBinArray = array('d',nSegmentsBinArray)
106 
107  outhisto = TH3D(aKey+"_"+aJetDefProv,aKey+"_"+aJetDefProv,nPTBins,pTBinArray,nNSegBins,nSegmentsBinArray,nEtaBins,etaBinArray)
108  outhisto.SetDirectory(0)
109 
110  # Fix histo margins to make sure it only applies the calibration where it should.
111  # Loop over all bins in x, y, z.
112  for xbin in range(histo.GetNbinsX()+2) :
113  for ybin in range(histo.GetNbinsY()+2) :
114  for zbin in range(histo.GetNbinsZ()+2) :
115  # Permitted regions are only:
116  # - pT > 50 GeV
117  # - eta < 2.7
118  # - NSegments > 20
119 
120  if "EM" in aJetDef :
121  if histo.GetXaxis().GetBinCenter(xbin) < 50 or \
122  histo.GetYaxis().GetBinCenter(ybin) > 2.7 or \
123  histo.GetZaxis().GetBinCenter(zbin) < 20 :
124  histo.SetBinContent(xbin,ybin,zbin,0.0)
125  outhisto.SetBinContent(xbin,zbin,ybin,0.0)
126  else :
127  outhisto.SetBinContent(xbin,zbin,ybin,histo.GetBinContent(xbin,ybin,zbin))
128 
129  elif "LC" in aJetDef :
130  if histo.GetXaxis().GetBinCenter(xbin) < 50 or \
131  histo.GetZaxis().GetBinCenter(zbin) > 2.7 or \
132  histo.GetYaxis().GetBinCenter(ybin) < 20 :
133  histo.SetBinContent(xbin,ybin,zbin,0.0)
134  outhisto.SetBinContent(xbin,ybin,zbin,0.0)
135  else :
136  outhisto.SetBinContent(xbin,ybin,zbin,histo.GetBinContent(xbin,ybin,zbin))
137 
138 
139  histos[aJetDefProv][aKey] = outhisto
140 
141  # Add a blank histogram for AFII corresponding to this
142  else :
143  outhisto = TH3D(aKey+"_"+aJetDefProv,aKey+"_"+aJetDefProv,2,15.0,4500.0,2,0.0,1000.0,2,-4.5,4.5)
144  outhisto.SetDirectory(0)
145  histos[aJetDefProv][aKey] = outhisto
146 
147  # Done reading, close the file
148  rootFile.Close()
149 
150  return histos
151 
ParsePunchthroughInput.ReadPunchthroughHistograms
def ReadPunchthroughHistograms(dirName)
Definition: Final2012/ParsePunchthroughInput.py:31
TH3D
Definition: rootspy.cxx:505
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
ParsePunchthroughInput.GetKeyNames
def GetKeyNames(self, dir="")
Definition: Final2012/ParsePunchthroughInput.py:10
DerivationFramework::TriggerMatchingUtils::sorted
std::vector< typename T::value_type > sorted(T begin, T end)
Helper function to create a sorted vector from an unsorted one.
array
Generated on Sat Oct 5 2024 21:23:31 for ATLAS Offline Software by   doxygen 1.8.18