7 Dumps the beam spot parameters to a text (ascii) file
9 __author__ =
'Peter Loscutoff'
11 __usage__ =
'plotBeamSpotVxVal.py [options] bs.nt.root split.nt.root'
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')
26 (options,args) = parser.parse_args()
28 parser.error(
'wrong number of command line arguments')
32 os.unsetenv(
'DISPLAY')
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
40 from array
import array
42 canvas = TCanvas(
"PlotBsVxVal",
"PlotBsVxVal",1000,500)
43 canvas.SetFillColor(0)
46 bsfile = TFile(args[0])
47 splitfile = TFile(args[1])
59 plotrange =
int(options.plotrange) + 1
61 plotSV = TH1D(
"p(s)SV",
"Pull vs # Tracks",plotrange*2,-.25,plotrange-.25)
63 for i
in range(0,plotrange*2):
64 nameSV =
"plotSV%s" % i
65 subplotSVs[nameSV]=TH1D(nameSV,nameSV,fitbins,-1*outlierCut,outlierCut)
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) )
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())
83 plotSV.SetBinContent(bin,fitsigma)
84 plotSV.SetBinError(bin,fserror)
87 gPad.SetLeftMargin(.15)
88 gPad.SetBottomMargin(.15)
90 plotSV.SetXTitle(
"Average # Tracks in Split Vertices")
91 plotSV.SetYTitle(
"Vertex Error Scale Factor")
94 plotSV.SetMarkerStyle(1)
97 legend = TLegend(.60,.75,.92,.92)
98 legend.AddEntry(plotSV,
"K Factor from Splitting")
99 legend.SetFillColor(0)
100 legend.SetBorderSize(0)
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)
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):
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)
143 plotX = TH1D(
"p(s)X",
"k vs Lumi Block",lumibins,lumiedge)
144 plotBS = TH1D(
"p(s)bs",
"k vs Lumi Block",lumibins,lumiedge)
146 for i
in range(0,21):
147 splitWeight.append(0.0)
148 fullWeight.append(0.0)
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
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
171 for i
in range(0,len(splitWeight)):
172 if splitWeight[i] != 0
and fullWeight[i] != 0:
173 splitWeight[i] = fullWeight[i]/splitWeight[i]
176 sptree = gDirectory.Get(
"splitVertex")
177 entries = sptree.GetEntriesFast()
178 for entry
in range(entries):
179 noprint = sptree.GetEntry(entry)
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"]
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 ):
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)
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:
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)
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:
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])
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])
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:
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
260 if gDirectory.Get(
'BeamSpotNt'):
261 bstree = gDirecotry.Get(
'BeamSpotNt')
262 elif gDirectory.Get(
'Beamspot/Beamspots'):
263 bstree = gDirectory.Get(
'Beamspot/Beamspots')
265 entries = bstree.GetEntriesFast()
266 for i
in range(entries):
267 noprint = bstree.GetEntry(i)
269 fserror = sqrt(bstree.kk)
270 if fitsigma > 0
and fserror > 0:
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
284 print (
"Splitting: ",vsavg,
" +- ",vserr)
285 print (
"Beam spot: ",bsavg,
" +- ",bserr)
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)
300 plotBS.SetMarkerStyle(21)
301 plotBS.SetLineWidth(3)
302 plotBS.Draw(
"same,e")
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)
313 for o
in options.output.split(
','):
315 canvas.SaveAs(basename+what+o)
320 os.environ[
'PYTHONINSPECT'] =
'1'