ATLAS Offline Software
beamspotVtxAnalysis.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
4 
5 from optparse import OptionParser
6 import os, glob
7 import string
8 from math import sqrt
9 
10 __author__ = 'James Walder'
11 __version__ = '$Id$'
12 __usage__ = '%prog [options] '
13 
14 def makeListFromString(s,token=','):
15  names = s.split(token)
16  return names
17 
18 parser = OptionParser(usage=__usage__, version=__version__)
19 parser.add_option('-f','--files',dest='files', default='vtx.root', help='Root files for data vertex plots')
20 parser.add_option('-g','--bs',dest='beamspotFile', default='', help='Root files for data beamspot plots')
21 
22 parser.add_option('-d','--dir',dest='dir', default='Beamspot', help='Root directory for beamspots')
23 parser.add_option('-b','--batch', dest='batch', action='store_true', default=False, help='Run in batch mode')
24 
25 (options,args) = parser.parse_args()
26 
27 if options.batch:
28  os.unsetenv('DISPLAY')
29 
30 
31 from ROOT import TH1D, TCanvas
32 from ROOT import gSystem, gROOT, TFile, TH1D, TH2D, TCanvas, TTree, TChain
33 from ROOT import TGraphErrors
34 from ROOT import TMultiGraph
35 from ROOT import EColor
36 from ROOT import TLegend,TLine
37 from ROOT import TLatex,TGaxis, gFile
38 
39 legends = [] # holder
40 pdglines= [] # holder
41 labels = []
42 ghists = [] # additional hists
43 failed = [] # list of failed plots
44 graphs = []
45 canvases = []
46 
47 PointSize = 1.6
48 MaxHeightFactor = 1.5
49 
50 colours = [ EColor.kRed+1, EColor.kBlue+1, EColor.kGreen-8, EColor.kYellow+2]
51 markers = [20,25,22,27,23,28,30,21]
52 plotTextSize = 0.03
53 
54 gROOT.ProcessLine(".L AtlasStyle.C")
55 gROOT.ProcessLine(".L AtlasUtils.C")
56 gROOT.ProcessLine("SetAtlasStyle();")
57 
58 def findBeamspots(file,dir, pattern):
59  import re
60  from ROOT import gDirectory
61  file.cd(dir)
62  p = re.compile(pattern)
63  print p.pattern
64  n = gDirectory.GetListOfKeys().GetSize()
65  print "N: ", n
66  matched = []
67  for i in range(n):
68  name = gDirectory.GetListOfKeys().At(i).GetName()
69  if not p.match(name):
70  continue
71  print name
72  #print p.match(name)
73  h = gDirectory.Get(gDirectory.GetListOfKeys().At(i).GetName())
74  #print h, h.GetName()
75  matched.append(h)
76 
77  file.cd()
78  return matched
79 
80 def makePoint( tree,prefix ,xval='x0',xlow=-8, xhigh=8, xbins=100, doGausFit=True, cuts=""):
81  h = TH1D(prefix+"_"+xval, "h_"+tree.GetName()+"_"+xval, xbins, xlow,xhigh)
82  p = '(%(x)s) >> %(h)s' % { 'x':xval, 'h':h.GetName()}
83  h.SetXTitle(xval)
84  ylabel = "Entries / XXX"
85  n = h.GetBinWidth(1)
86  l = '%(n).2f'%{'n':n}
87  s = ylabel.replace("XXX",l)
88  h.SetYTitle(s)
89 
90  print p
91  tree.Draw( p , cuts);
92  h.Print()
93 
94  mean = h.GetMean()
95  meanErr = h.GetMeanError()
96  rms = h.GetRMS()
97  rmsErr = h.GetRMSError()
98  if doGausFit:
99  h.Fit("gaus","LL")
100  f = h.GetFunction("gaus")
101  mean = f.GetParameter(1)
102  meanErr = f.GetParError(1)
103  rms = f.GetParameter(2)
104  rmsErr = f.GetParError(2)
105  ghists.append(h)
106  print h.GetName(), mean, meanErr, rms, rmsErr
107  return (mean, meanErr, rms, rmsErr,h)
108 
109 
110 def AddText(xmin,ymin,text,size=0.015, font=None):
111  m = TLatex() # //l.SetTextAlign(12); l.SetTextSize(tsize);
112  m.SetNDC();
113  #m.SetTextFont(72);
114  m.SetTextColor(EColor.kBlack);
115  m.SetTextSize(size)
116  if font != None:
117  m.SetTextFont(font)
118  m.DrawLatex(xmin,ymin,text);
119  labels.append(m)
120 
121 
123  subs = { 'p':'.', 'm':'-' }
124  n = ""
125  for i in s:
126  c = i
127  if i in subs.keys():
128  c = subs[i]
129  n = n+c
130  return float(n)
131 
132 
133 def splitBeamspotTreeTitle(title, cutToken="CUT"):
134  t = title.split(cutToken)
135  namestring = t[0].strip('_')
136  cutstring = t[-1].strip('_')
137  cuts = cutstring.split('_')
138  output = {}
139  for i in range(0,len(cuts),2):
140  output[cuts[i]] = makeFloatFromString(cuts[i+1])
141 
142  return output,namestring
143 
144 def getTree(name, file, dir):
145  tr = file.Get(dir +"/" + name)
146  print type(tr)
147  return tr
148 
149 
150 # Start of main programs
151 
152 #Open the files
153 
154 files = makeListFromString(options.files)
155 if len(files) == 1:
156  fii = TFile(files[0])
157 else:
158  print "Multiple files not yet supported"
159  sys.exit()
160 
161 
162 
163 
164 def getObject(chain,name):
165  try:
166  h = chain.Get(name)
167  #h=chain.Get(name)
168  h.Print()
169  except:
170  print "No object in ", chain.GetName(), dir,name
171  h = None
172  return h
173 
174 
175 def AddStatText(h, minx=0.7,miny=0.7,title="",units="",pars=['Mean','RMS']):
176  n = h.GetEntries()
177  En = sqrt(n)
178 
179  m = h.GetMean()
180  Em = h.GetMeanError()
181 
182  r = h.GetRMS()
183  Er = h.GetRMSError()
184 
185  step =0.033
186  ypos = miny;
187  if len(title) > 0:
188  AddText(minx, ypos, title,plotTextSize)
189  ypos = ypos - step
190  if "Entries" in pars:
191  nsig = 'N = %(v).0f #pm %(ev).0f' % { 'v': n, 'ev': En}
192  AddText(minx, ypos, nsig,plotTextSize)
193  ypos = ypos - step
194  if "Mean" in pars:
195  mass = '#mu= %(v).04f #pm %(ev).04f %(u)s'% { 'v': m, 'ev': Em, 'u':units}
196  AddText(minx, ypos, mass,plotTextSize)
197  ypos = ypos - step
198  if "RMS" in pars:
199  sigma= 'RMS = %(v)2.4f #pm %(ev)2.4f %(u)s'% { 'v': r, 'ev': Er, 'u':units}
200  AddText(minx, ypos, sigma,plotTextSize)
201  ypos = ypos - step
202 
203 
204 def makeHists( tree, name):
205  points= []
206 
207  cuts='nEvents>5000&&fitID==1&&fitStatus==1'
208 
209  xpos = makePoint( tree,name, xval='xc',xbins=100,xlow=-0.16,xhigh=-0.12,doGausFit=False,cuts=cuts)[-1]
210  ypos = makePoint( tree,name, xval='yc',xbins=100,xlow=0.99,xhigh=1.02,doGausFit=False,cuts=cuts)[-1]
211  zpos = makePoint( tree,name, xval='z', xbins=50,xlow=-2,xhigh=3,doGausFit=False,cuts=cuts)[-1]
212 
213  sx = makePoint( tree,name,xval='sx',xbins=100,xlow=0.06,xhigh=0.1,doGausFit=False,cuts=cuts)[-1]
214  sy = makePoint( tree,name, xval='sy',xbins=100,xlow=0.06,xhigh=0.1,doGausFit=False,cuts=cuts)[-1]
215  sz = makePoint( tree,name, xval='sz',xbins=100,xlow=55,xhigh=65,doGausFit=False,cuts=cuts)[-1]
216  ax = makePoint( tree,name, xval='ax',xbins=50,xlow=400e-6,xhigh=700e-6,doGausFit=False,cuts=cuts)[-1]
217  ay = makePoint( tree,name, xval='ay',xbins=50,xlow=-20e-6,xhigh=100e-6,doGausFit=False,cuts=cuts)[-1]
218  k = makePoint( tree,name, xval='k', xbins=100,xlow=0.9,xhigh=1.3,doGausFit=False,cuts=cuts)[-1]
219  rhoxy = makePoint( tree,name, xval='rhoxy',xbins=100,xlow=-0.3,xhigh=0.3,doGausFit=False,cuts=cuts)[-1]
220 
221  return (xpos, ypos, zpos, sx, sy, sz, ax, ay, k, rhoxy)
222 
223 
224 
225 def makeCanvas(name, hists, xlabel, ylabel, doStats=True,doLegend=True, units="", logy=False,setMinZero=True,drawOptions="",histLabels=[]):
226  canv = TCanvas(name, "",800,600)
227 
228  min=10e6
229  max=-10e6
230  for i in hists:
231  m = i.GetMaximum()
232  max = ( m if m>max else max)
233  m = i.GetMinimum()
234  min = (m if m < min else min)
235 
236  if setMinZero and not logy:
237  min=0
238  if min<=0 and logy:
239  min=0.01
240  #if min > 1 and logy:
241  # min=0.1
242 
243  c=0
244  m=0
245  isFirst=True
246  for i in hists:
247  i.SetMarkerColor(colours[c])
248  i.SetLineColor(colours[c])
249  i.SetMarkerStyle(markers[m])
250  c = c+1
251  if c == len(colours):
252  m = m+1
253  c=0
254  if m == len(markers):
255  m=0
256 
257  if isFirst:
258  i.SetMinimum( min)
259  i.SetMaximum( max*1.5)
260  i.SetXTitle(xlabel)
261  if "XXX" not in ylabel :
262  print "XX"
263  s = ylabel
264  else:
265  print "YY"
266  n = i.GetBinWidth(1)
267  l = '%(n).4f'%{'n':n}
268  s = ylabel.replace("XXX",l)
269  i.SetYTitle(s)
270  i.Draw(drawOptions)
271  isFirst = False
272  else:
273  i.Draw(drawOptions+"sames")
274  if logy:
275  canv.SetLogy(1)
276  if doStats:
277  pos = 0.9
278  for i in range(len(hists)):
279  if len(histLabels) == 0:
280  AddStatText(hists[i], minx=0.7, miny=pos, title = hists[i].GetName())
281  else:
282  AddStatText(hists[i], minx=0.7, miny=pos, title = histLabels[i])
283  pos = pos -0.15
284  if doLegend:
285  leg = TLegend(0.2,0.7,0.4,0.9)
286  leg.SetFillColor(0);
287  leg.SetFillStyle(0);
288  leg.SetBorderSize(0);
289  for i in range(len(hists)):
290  if len(histLabels) == 0:
291  leg.AddEntry(hists[i], hists[i].GetName())
292  else:
293  leg.AddEntry(hists[i], histLabels[i])
294 
295  legends.append(leg)
296  leg.Draw()
297 
298  return canv
299 
300 
301 def getObjects(file, prefix,names):
302  o=[]
303  for i in names:
304  o.append(getObject(file,prefix+i))
305  return o
306 
307 def makeStandardPlots(file, prefix,names=[],nameLabels=[], doMC=False, logy=True):
308 
309  cs = []
310  cs.append(makeCanvas(prefix+"_nTracks",getObjects(file,"hnTracks_",names)
311  , xlabel="nTracks", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
312 
313  cs.append(makeCanvas(prefix+"_chi2",getObjects(file,"hchi2_",names)
314  , xlabel="chi2", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
315  cs.append(makeCanvas(prefix+"_prob",getObjects(file,"hprob_",names)
316  , xlabel="prob", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
317  cs.append(makeCanvas(prefix+"_sumpt",getObjects(file,"hsumpt_",names)
318  , xlabel="sumpt", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
319 
320  cs.append(makeCanvas(prefix+"_x",getObjects(file,"hvx_",names)
321  ,xlabel="x [mm]", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
322  cs.append(makeCanvas(prefix+"_y",getObjects(file,"hvy_",names)
323  ,xlabel="y [mm]", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
324  cs.append(makeCanvas(prefix+"_z",getObjects(file,"hvz_",names)
325  ,xlabel="z [mm]", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
326 
327  cs.append(makeCanvas(prefix+"_Ex",getObjects(file,"hvEx_",names)
328  ,xlabel="#sigma(x) [mm]", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
329  cs.append(makeCanvas(prefix+"_Ey",getObjects(file,"hvEy_",names)
330  ,xlabel="#sigma(y) [mm]", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
331  cs.append(makeCanvas(prefix+"_Ez",getObjects(file,"hvEz_",names)
332  ,xlabel="#sigma(z) [mm]", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
333 
334  if doMC:
335  cs.append(makeCanvas(prefix+"_Pullx",getObjects(file,"hvxPullMC_",names)
336  ,xlabel="Pull x", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
337  cs.append(makeCanvas(prefix+"_Pully",getObjects(file,"hvyPullMC_",names)
338  ,xlabel="Pull y", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
339  cs.append(makeCanvas(prefix+"_Pullz",getObjects(file,"hvzPullMC_",names)
340  ,xlabel="Pull z", ylabel = "Entries/ XXX", logy=logy,histLabels=nameLabels))
341 
342 
343  return cs
344 
345 def makeDeltaPlots(file, prefix,names=[],nameLabels=[]):
346 
347  cs = []
348  cs.append(makeCanvas(prefix+"_DnTracks",getObjects(file,"hDnTracks_",names)
349  , xlabel="#DeltanTracks", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
350  cs.append(makeCanvas(prefix+"_Dsumpt",getObjects(file,"hDsumpt_",names)
351  , xlabel="#Deltasumpt", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
352 
353  cs.append(makeCanvas(prefix+"_Dx",getObjects(file,"hvDx_",names)
354  ,xlabel="#Deltax [mm]", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
355  cs.append(makeCanvas(prefix+"_Dy",getObjects(file,"hvDy_",names)
356  ,xlabel="#Deltay [mm]", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
357  cs.append(makeCanvas(prefix+"_Dz",getObjects(file,"hvDz_",names)
358  ,xlabel="#Deltaz [mm]", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
359 
360  cs.append(makeCanvas(prefix+"_DEx",getObjects(file,"hvDEx_",names)
361  ,xlabel="#Delta#sigma(x) [mm]", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
362  cs.append(makeCanvas(prefix+"_DEy",getObjects(file,"hvDEy_",names)
363  ,xlabel="#Delta#sigma(y) [mm]", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
364  cs.append(makeCanvas(prefix+"_DEz",getObjects(file,"hvDEz_",names)
365  ,xlabel="#Delta#sigma(z) [mm]", ylabel = "Entries/ XXX", logy=True,histLabels=nameLabels))
366 
367 
368 
369  return cs
370 
371 
372 # For the analys
373 #Open the files
374 
375 files = makeListFromString(options.files)
376 if len(files) == 1:
377  fii = TFile(files[0])
378 else:
379  print "Multiple files not yet supported"
380  sys.exit()
381 
382 
383 #chainVtx = TChain("")
384 #for i in files:
385 # chainVtx.Add(i)
386 #print "Found: ", chainVtx.GetNtrees()
387 
388 fii = TFile(files[0])
389 
390 
391 cs = []
392 
393 #cs.extend( makeStandardPlots(fii,prefix="VtxAllCombRes", names=['AllDefault','AllRes'],nameLabels=['Combined','Resolved']) )
394 #cs.extend( makeStandardPlots(fii,prefix="VtxMatCombRes", names=['MatchedDefault','resMatched'],nameLabels=['Combined','Resolved']) )
395 
396 #cs.extend( makeStandardPlots(fii,prefix="VtxAllcombpixsct", names=['AllDefault','AllPix','AllSCT'],nameLabels=['Combined','Pix','SCT']) )
397 
398 cs.extend( makeStandardPlots(fii,prefix="VtxAllCombRespixsct", names=['AllDefault','AllRes','AllPix','AllSCT'],nameLabels=['Combined','Resolved','Pix','SCT']) )
399 cs.extend( makeStandardPlots(fii,prefix="VtxMatcombrespixsct", names=['MatchedDefault','resMatched','pixMatched','sctMatched'],nameLabels=['Combined','resMatched','Pix','SCT'],doMC=True) )
400 cs.extend( makeStandardPlots(fii,prefix="VtxMatcombrespixsctNoLog", names=['MatchedDefault','resMatched','pixMatched','sctMatched'],nameLabels=['Combined','resMatched','Pix','SCT'],doMC=True,logy=False) )
401 
402 cs.extend( makeDeltaPlots(fii, prefix="VtxDpixcomb", names=['pixMatchedComb'],nameLabels=['Comb-Pix']) )
403 cs.extend( makeDeltaPlots(fii, prefix="VtxDrescomb", names=['resMatchedComb'],nameLabels=['Comb-Res']) )
404 cs.extend( makeDeltaPlots(fii, prefix="VtxDsctcomb", names=['sctMatchedComb'],nameLabels=['Comb-SCT']) )
405 cs.extend( makeDeltaPlots(fii, prefix="VtxDpixsct", names=['sctMatchedPix'], nameLabels=['Pix-SCT']) )
406 
407 cs.extend( makeDeltaPlots(fii, prefix="VtxDrespixsctcomb", names=['resMatchedComb','pixMatchedComb','sctMatchedComb'],nameLabels=['Res-Comb','Pix-Comb','SCT-Comb']) )
408 
409 
410 # def makeCanvas(name, hists, xlabel, ylabel, doStats=True,doLegend=True, units="", logy=False,setMinZero=True,drawOptions="",histLabels=[]):
411 
412 
413 
414 if len(options.beamspotFile):
415  bsfii = TFile(options.beamspotFile)
416 
417  bsfii.cd()
418  bs = findBeamspots(bsfii,"Beamspot","Beamspots.*")
419 
420  type = ['Adaptive','InDetPriVxFinderFastFinder','InDetPriVxFinderFullFinder']
421  for t in type:
422  bsAdaptive = findBeamspots(bsfii,"Beamspot","Beamspots.*"+t+".*")
423 
424  hx=[]
425  hy=[]
426  hz=[]
427  hsx=[]
428  hsy=[]
429  hsz=[]
430  hax=[]
431  hay=[]
432  hk=[]
433  hrhoxy=[]
434 
435 
436 
437  for i in bsAdaptive:
438  if "TRT" in i.GetName():
439  continue
440  if "SCT" in i.GetName():
441  name = "SCT"
442  elif "Pix" in i.GetName():
443  name = "Pix"
444  elif "Combined" in i.GetName():
445  name = "Combined"
446  elif "Res" in i.GetName():
447  name = "Resolved"
448  else:
449  name = i.GetName()
450 
451 
452  hists = makeHists(i,name)
453  hx.append(hists[0])
454  hy.append(hists[1])
455  hz.append(hists[2])
456  hsx.append(hists[3])
457  hsy.append(hists[4])
458  hsz.append(hists[5])
459  hax.append(hists[6])
460  hay.append(hists[7])
461  hk.append(hists[8])
462  hrhoxy.append(hists[9])
463 
464  cs.append( makeCanvas(t+"_x",hx,"x [mm]","Entries / (XXX mm)"))
465  cs.append( makeCanvas(t+"_y",hy,"y [mm]","Entries / (XXX mm)"))
466  cs.append( makeCanvas(t+"_z",hz,"z [mm]","Entries / (XXX mm)"))
467  cs.append( makeCanvas(t+"_sx",hsx,"sx [mm]","Entries / (XXX mm)"))
468  cs.append( makeCanvas(t+"_sy",hsy,"sy [mm]","Entries / (XXX mm)"))
469  cs.append( makeCanvas(t+"_sz",hsz,"sz [mm]","Entries / (XXX mm)"))
470  cs.append( makeCanvas(t+"_ax",hax,"ax [mm]","Entries / (XXX mm)"))
471  cs.append( makeCanvas(t+"_ay",hay,"ay [mm]","Entries / (XXX mm)"))
472  cs.append( makeCanvas(t+"_k",hk,"k [mm]","Entries / (XXX mm)"))
473  cs.append( makeCanvas(t+"_rhoxy",hrhoxy,"rhoxy [mm]","Entries / (XXX mm)"))
474 
475 for i in cs:
476  i.SaveAs("plots/"+i.GetName()+".eps")
477  i.SaveAs("plots/"+i.GetName()+".png")
478 
479 
480 
481 s = raw_input('--> ')
482 exit()
483 
beamspotVtxAnalysis.AddText
def AddText(xmin, ymin, text, size=0.015, font=None)
Definition: beamspotVtxAnalysis.py:110
makeHists
Definition: makeHists.py:1
beamspotVtxAnalysis.makePoint
def makePoint(tree, prefix, xval='x0', xlow=-8, xhigh=8, xbins=100, doGausFit=True, cuts="")
Definition: beamspotVtxAnalysis.py:80
beamspotVtxAnalysis.makeHists
def makeHists(tree, name)
Definition: beamspotVtxAnalysis.py:204
beamspotVtxAnalysis.makeFloatFromString
def makeFloatFromString(s)
Definition: beamspotVtxAnalysis.py:122
beamspotVtxAnalysis.splitBeamspotTreeTitle
def splitBeamspotTreeTitle(title, cutToken="CUT")
Definition: beamspotVtxAnalysis.py:133
beamspotVtxAnalysis.getObject
def getObject(chain, name)
Definition: beamspotVtxAnalysis.py:164
beamspotVtxAnalysis.getObjects
def getObjects(file, prefix, names)
Definition: beamspotVtxAnalysis.py:301
beamspotVtxAnalysis.type
list type
Definition: beamspotVtxAnalysis.py:420
beamspotVtxAnalysis.AddStatText
def AddStatText(h, minx=0.7, miny=0.7, title="", units="", pars=['Mean', 'RMS'])
Definition: beamspotVtxAnalysis.py:175
beamspotVtxAnalysis.makeCanvas
def makeCanvas(name, hists, xlabel, ylabel, doStats=True, doLegend=True, units="", logy=False, setMinZero=True, drawOptions="", histLabels=[])
Definition: beamspotVtxAnalysis.py:225
beamspotVtxAnalysis.getTree
def getTree(name, file, dir)
Definition: beamspotVtxAnalysis.py:144
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
beamspotVtxAnalysis.findBeamspots
def findBeamspots(file, dir, pattern)
Definition: beamspotVtxAnalysis.py:58
calibdata.exit
exit
Definition: calibdata.py:236
beamspotVtxAnalysis.makeListFromString
def makeListFromString(s, token=',')
Definition: beamspotVtxAnalysis.py:14
beamspotVtxAnalysis.makeStandardPlots
def makeStandardPlots(file, prefix, names=[], nameLabels=[], doMC=False, logy=True)
Definition: beamspotVtxAnalysis.py:307
beamspotVtxAnalysis.makeDeltaPlots
def makeDeltaPlots(file, prefix, names=[], nameLabels=[])
Definition: beamspotVtxAnalysis.py:345
readCCLHist.float
float
Definition: readCCLHist.py:83