ATLAS Offline Software
plotBeamSpotVxVal.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
4 
5 
6 """
7 Dumps the beam spot parameters to a text (ascii) file
8 """
9 __author__ = 'Peter Loscutoff'
10 __version__ = '$Id $'
11 __usage__ = 'plotBeamSpotVxVal.py [options] bs.nt.root split.nt.root'
12 
13 import sys, os
14 import time
15 from math import sqrt
16 
17 # Argument parsing
18 from optparse import OptionParser
19 parser = OptionParser(usage=__usage__, version=__version__)
20 parser.add_option('-b', '--batch', dest='batch', action='store_true', default=False, help='run in batch mode')
21 parser.add_option('-o', '--output', dest='output', default='.gif', help='comma-separated list of output files or formats')
22 parser.add_option('-H', '--test', dest='test', default=False, help='a test variable (NYI)')
23 parser.add_option('', '--time', dest='time', default=False, help='use time instead of lumiblocks as an axis (NYI)')
24 parser.add_option('-t', '--tracks', dest='plotrange', default=100, help='maximum number of tracks/vertex to plot')
25 
26 (options,args) = parser.parse_args()
27 if len(args) < 2:
28  parser.error('wrong number of command line arguments')
29 
30 # Setup ROOT
31 if options.batch:
32  os.unsetenv('DISPLAY')
33 
34 import math
35 from math import sqrt,sin,cos,floor,fabs
36 from ROOT import TFile,TH1D,gDirectory,TTree,TCanvas,TLegend,gStyle,gROOT,TMatrixDSym,gPad
37 from InDetBeamSpotExample import ROOTUtils
38 from InDetBeamSpotExample.BeamSpotData import *
39 import time, calendar
40 from array import array
41 
42 canvas = TCanvas("PlotBsVxVal","PlotBsVxVal",1000,500)
43 canvas.SetFillColor(0)
44 canvas.Divide(2,1)
45 
46 bsfile = TFile(args[0])
47 splitfile = TFile(args[1])
48 
49 # Define our tools to get the time and fill number:
50 #from COOLUtils import RunTimeQuery
51 #c = RunTimeQuery()
52 #from CoolConvUtilities import AtlCoolTool
53 #tool = AtlCoolTool.AtlCoolTool("COOLOFL_DCS/COMP200")
54 
55 subplotSVs = {}
56 ROOTUtils.setStyle()
57 
58 fitbins = 3001
59 plotrange = int(options.plotrange) + 1
60 outlierCut = 20
61 plotSV = TH1D("p(s)SV","Pull vs # Tracks",plotrange*2,-.25,plotrange-.25)
62 
63 for i in range(0,plotrange*2):
64  nameSV = "plotSV%s" % i
65  subplotSVs[nameSV]=TH1D(nameSV,nameSV,fitbins,-1*outlierCut,outlierCut)
66 
67 splitfile.cd()
68 tree = gDirectory.Get("splitVertex")
69 entries = tree.GetEntriesFast()
70 for entry in range(entries):
71  noprint = tree.GetEntry(entry)
72  binSV = (tree.tracks_odd+tree.tracks_even)
73  nameSV = "plotSV%i" % binSV
74  if (binSV < (plotrange-1)*2 and tree.tracks_odd >= 4 and tree.tracks_even >= 4):
75  subplotSVs[nameSV].Fill( (tree.x_odd - tree.x_even)/sqrt(tree.c00_odd+tree.c00_even) )
76 
77 for i in range(0,plotrange*2):
78  nameSV = "plotSV%s" % i
79  if subplotSVs[nameSV].Integral() > 2:
80  fitsigma = subplotSVs[nameSV].GetRMS()
81  fserror = (subplotSVs[nameSV].GetRMS())/sqrt(subplotSVs[nameSV].Integral())
82  bin = i+1
83  plotSV.SetBinContent(bin,fitsigma)
84  plotSV.SetBinError(bin,fserror)
85 
86 canvas.cd(1)
87 gPad.SetLeftMargin(.15)
88 gPad.SetBottomMargin(.15)
89 
90 plotSV.SetXTitle("Average # Tracks in Split Vertices")
91 plotSV.SetYTitle("Vertex Error Scale Factor")
92 #plotSV.SetMarkerColor(4)
93 #plotSV.SetMarkerStyle(26)
94 plotSV.SetMarkerStyle(1)
95 plotSV.Draw("e")
96 
97 legend = TLegend(.60,.75,.92,.92)
98 legend.AddEntry(plotSV,"K Factor from Splitting")
99 legend.SetFillColor(0)
100 legend.SetBorderSize(0)
101 
102 legend.Draw()
103 
104 subplots = {}
105 
106 bstable = {}
107 
108 splitWeight = []
109 fullWeight = []
110 
111 lumiedge = array('d')
112 lumibins = 0
113 lumiend = 0
114 
115 bsfile.cd()
116 if gDirectory.Get("BeamSpotNt"):
117  bstree = gDirectory.Get("BeamSpotNt")
118 elif gDirectory.Get("Beamspot/Beamspots"):
119  bstree = gDirectory.Get("Beamspot/Beamspots")
120 entries = bstree.GetEntriesFast()
121 for i in range(entries):
122  noprint = bstree.GetEntry(i)
123  if bstree.k != 0:
124  lumibegin = bstree.lbStart
125  lumiedge.append(lumibegin)
126  lumiend = bstree.lbEnd
127  name = str(lumibegin) + "to" + str(lumiend)
128  lumibins = lumibins + 1
129  subplots[name]=TH1D(name,name,fitbins,-100,100)
130  for lb in range(lumibegin,lumiend+1):
131  bstable[lb] = {}
132  bstable[lb]["x0"] = bstree.x0
133  bstable[lb]["y0"] = bstree.y0
134  bstable[lb]["ax"] = bstree.ax
135  bstable[lb]["ay"] = bstree.ay
136  bstable[lb]["z"] = bstree.z
137  bstable[lb]["sx"] = bstree.sx
138  bstable[lb]["sy"] = bstree.sy
139  bstable[lb]["sz"] = bstree.sz
140  bstable[lb]["rhoxy"] = bstree.rhoxy
141 lumiedge.append(lumiend)
142 
143 plotX = TH1D("p(s)X","k vs Lumi Block",lumibins,lumiedge)
144 plotBS = TH1D("p(s)bs","k vs Lumi Block",lumibins,lumiedge)
145 
146 for i in range(0,21):
147  splitWeight.append(0.0)
148  fullWeight.append(0.0)
149 
150 splitfile.cd()
151 sptree = gDirectory.Get("splitVertex")
152 entries = sptree.GetEntriesFast()
153 for entry in range(entries):
154  noprint = sptree.GetEntry(entry)
155  splitBin = int(.5*(sptree.tracks_odd + sptree.tracks_even) )
156  if splitBin < 11: splitWeight[splitBin] = splitWeight[splitBin] + 1.0
157  else: splitWeight[11] = splitWeight[11] + 1.0
158 
159 bsfile.cd()
160 if gDirectory.Get('Vertices'):
161  vxtree = gDirectory.Get('Vertices')
162 elif gDirectory.Get('Beamspot/Vertices'):
163  vxtree = gDirectory.Get('Beamspot/Vertices')
164 entries = vxtree.GetEntriesFast()
165 for entry in range(entries):
166  noprint = vxtree.GetEntry(entry)
167  fullBin = vxtree.nTracks
168  if fullBin < 11: fullWeight[fullBin] = fullWeight[fullBin] + 1.0
169  else: fullWeight[11] = fullWeight[11] + 1.0
170 
171 for i in range(0,len(splitWeight)):
172  if splitWeight[i] != 0 and fullWeight[i] != 0:
173  splitWeight[i] = fullWeight[i]/splitWeight[i]
174 
175 splitfile.cd()
176 sptree = gDirectory.Get("splitVertex")
177 entries = sptree.GetEntriesFast()
178 for entry in range(entries):
179  noprint = sptree.GetEntry(entry)
180  lumi = sptree.lumi
181  if not lumi in bstable.keys(): continue
182  bs_x0 = bstable[lumi]["x0"]
183  bs_y0 = bstable[lumi]["y0"]
184  bs_a = bstable[lumi]["ax"]
185  bs_b = bstable[lumi]["ay"]
186  bs_z = bstable[lumi]["z"]
187  bs_sx = bstable[lumi]["sx"]
188  bs_sy = bstable[lumi]["sy"]
189  bs_sz = bstable[lumi]["sz"]
190  bs_rho = bstable[lumi]["rhoxy"]
191 
192  name = ""
193  for k in subplots.keys():
194  range = k.split("to")
195  if lumi >= int(range[0]) and lumi <= int(range[1]): name = k
196  if name == "": continue
197  splitBin = int(.5*(sptree.tracks_odd + sptree.tracks_even) )
198  if (sptree.tracks_odd >= 4 and sptree.tracks_even >= 4 ):
199 
200  cov = array('d')
201  cov.append(bs_sx*bs_sx + sptree.c00_odd)
202  cov.append(bs_sx*bs_sy*bs_rho + sptree.c01_odd)
203  cov.append(bs_sx*bs_sy*bs_rho + sptree.c01_odd)
204  cov.append(bs_sy*bs_sy + sptree.c11_odd)
205  covmat = TMatrixDSym(2)
206  covmat.Use(2,cov)
207  covmat.Invert()
208  if (covmat(0,0)*(sptree.x_odd - (bs_x0 + sptree.z_odd*bs_a))**2 +
209  covmat(1,1)*(sptree.y_odd - (bs_y0 + sptree.z_odd*bs_b))**2 +
210  2*covmat(0,1)*(sptree.x_odd - (bs_x0 + sptree.z_odd*bs_a))*(sptree.y_odd - (bs_y0 + sptree.z_odd*bs_b))) > outlierCut:
211  continue
212  cov = array('d')
213  cov.append(bs_sx*bs_sx + sptree.c00_even)
214  cov.append(bs_sx*bs_sy*bs_rho + sptree.c01_even)
215  cov.append(bs_sx*bs_sy*bs_rho + sptree.c01_even)
216  cov.append(bs_sy*bs_sy + sptree.c11_even)
217  covmat = TMatrixDSym(2)
218  covmat.Use(2,cov)
219  covmat.Invert()
220  if (covmat(0,0)*(sptree.x_even - (bs_x0 + sptree.z_even*bs_a))**2 +
221  covmat(1,1)*(sptree.y_even - (bs_y0 + sptree.z_even*bs_b))**2 +
222  2*covmat(0,1)*(sptree.x_even - (bs_x0 + sptree.z_even*bs_a))*(sptree.y_even - (bs_y0 + sptree.z_even*bs_b))) > outlierCut:
223  continue
224 
225  if splitBin < 11:
226  subplots[name].Fill( (sptree.x_odd - sptree.x_even)/sqrt(sptree.c00_odd+sptree.c00_even) , splitWeight[splitBin])
227  subplots[name].Fill( (sptree.y_odd - sptree.y_even)/sqrt(sptree.c11_odd+sptree.c11_even) , splitWeight[splitBin])
228  else:
229  subplots[name].Fill( (sptree.x_odd - sptree.x_even)/sqrt(sptree.c00_odd+sptree.c00_even) , splitWeight[11])
230  subplots[name].Fill( (sptree.y_odd - sptree.y_even)/sqrt(sptree.c11_odd+sptree.c11_even) , splitWeight[11])
231 
232 min = .8
233 max = 1.2
234 
235 sumw=0
236 sumwx=0
237 
238 for name in subplots.keys():
239  if subplots[name].Integral() > 100:
240  midpoint = int(name.split("to")[0])+.5
241  fitsigma = subplots[name].GetRMS()
242  fserror = (subplots[name].GetRMS())/sqrt(subplots[name].Integral())
243  if fitsigma > 0 and fserror > 0:
244  sumw += 1/fserror
245  sumwx += fitsigma/fserror
246  globalbin = plotBS.FindBin(midpoint)
247  bin = plotBS.GetBin(globalbin)
248  plotX.SetBinContent(bin,fitsigma)
249  plotX.SetBinError(bin,fserror)
250  if fitsigma*1.3 > max: max = fitsigma*1.3
251  if fitsigma*.8 < min and fitsigma != 0: min = fitsigma*.8
252 
253 vsavg = sumwx/sumw
254 vserr = 1.0/sumw
255 
256 sumw=0
257 sumwx=0
258 
259 bsfile.cd()
260 if gDirectory.Get('BeamSpotNt'):
261  bstree = gDirecotry.Get('BeamSpotNt')
262 elif gDirectory.Get('Beamspot/Beamspots'):
263  bstree = gDirectory.Get('Beamspot/Beamspots')
264 #bstree = gDirectory.Get("BeamSpotNt")
265 entries = bstree.GetEntriesFast()
266 for i in range(entries):
267  noprint = bstree.GetEntry(i)
268  fitsigma = bstree.k
269  fserror = sqrt(bstree.kk)
270  if fitsigma > 0 and fserror > 0:
271  sumw += 1/fserror
272  sumwx += fitsigma/fserror
273  midpoint = bstree.lbStart+.5
274  globalbin = plotBS.FindBin(midpoint)
275  bin = plotBS.GetBin(globalbin)
276  plotBS.SetBinContent(bin,fitsigma)
277  plotBS.SetBinError(bin,fserror)
278  if fitsigma*1.3 > max: max = fitsigma*1.3
279  if fitsigma*.8 < min and fitsigma != 0: min = fitsigma*.8
280 
281 bsavg = sumwx/sumw
282 bserr = 1.0/sumw
283 
284 print ("Splitting: ",vsavg," +- ",vserr)
285 print ("Beam spot: ",bsavg," +- ",bserr)
286 
287 canvas.cd(2)
288 
289 gPad.SetLeftMargin(.15)
290 gPad.SetBottomMargin(.15)
291 plotX.SetXTitle("Lumi Block")
292 plotX.SetYTitle("Vertex Error Scale Factor")
293 plotX.SetMaximum(max)
294 plotX.SetMinimum(min)
295 plotX.SetLineColor(33)
296 plotX.SetMarkerColor(33)
297 plotX.SetMarkerStyle(1)
298 plotX.SetFillColor(33)
299 plotX.Draw("e2")
300 plotBS.SetMarkerStyle(21)
301 plotBS.SetLineWidth(3)
302 plotBS.Draw("same,e")
303 
304 legend = TLegend(.50,.70,.80,.85)
305 legend.AddEntry(plotBS,"k from beamspot")
306 legend.AddEntry(plotX,"k from vx split")
307 legend.SetFillColor(0)
308 legend.SetBorderSize(0)
309 
310 legend.Draw()
311 
312 if options.output:
313  for o in options.output.split(','):
314  if o[0]=='.':
315  canvas.SaveAs(basename+what+o)
316  else:
317  canvas.SaveAs(o)
318 
319 import os
320 os.environ['PYTHONINSPECT'] = '1'
321 
322 
BeamSpotData
ROOTUtils.setStyle
def setStyle(style=None)
Definition: roofit/ROOTUtils.py:413
plotBeamSpotVxVal.covmat
covmat
Definition: plotBeamSpotVxVal.py:205
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
array
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
str
Definition: BTagTrackIpAccessor.cxx:11
Generated on Mon Sep 22 2025 21:16:34 for ATLAS Offline Software by   doxygen 1.8.18