Belle II Software development
deltaTVertexTagVResolution.py
1#!/usr/bin/env python
2
3
10
11
25
26
27from basf2 import B2FATAL
28import ROOT
29import sys
30import glob
31import math
32from operator import itemgetter
33
34if len(sys.argv) != 4:
35 sys.exit(
36 "Must provide 3 arguments: [input_sim_file] or ['input_sim_file*'] wildcards, [treeName] " +
37 "and [VXD_requirement (PXD or SVD)]")
38
39PATH = "."
40
41workingFile = sys.argv[1]
42workingFiles = glob.glob(str(workingFile))
43
44limDeltaT = 5
45limDeltaTErr = 3.0
46limZSig = 0.03
47limZTag = 0.03
48
49treeName = str(sys.argv[2])
50
51# No PXD hit required: PXD0. At least one PXD (SVD) hit for one of the muon tracks: PXD1 (SVD1).
52# Hit required for both muon tracks: PXD2 (SVD2)"
53VXDReqs = []
54if str(sys.argv[3]) == "PXD":
55 VXDReqs = ["PXD0", "PXD2"]
56elif str(sys.argv[3]) == "SVD":
57 VXDReqs = ["SVD0", "SVD2"]
58else:
59 B2FATAL('Not available VXD requirement " + str(sys.argv[3]) + ". Available are "PXD" and "SVD".')
60
61ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
62
63fitResults = []
64numberOfEntries = []
65
66for VXDReq in VXDReqs:
67 iResult = []
68
69 tdat = ROOT.TChain(treeName)
70
71 for iFile in workingFiles:
72 tdat.AddFile(iFile)
73
74 B0_DeltaT = ROOT.RooRealVar("DeltaT", "DeltaT", 0.)
75 deltaTErr = ROOT.RooRealVar("DeltaTErr", "DeltaTErr", 0, limDeltaTErr, "ps")
76 B0_TruthDeltaT = ROOT.RooRealVar("MCDeltaT", "MCDeltaT", 0.)
77
78 B0_qrMC = ROOT.RooRealVar("qrMC", "qrMC", 0., -100, 100)
79 B0_isSignal = ROOT.RooRealVar("isSignal", "isSignal", 0., 0., 512.)
80
81 B0_TagVz = ROOT.RooRealVar("TagVz", "TagVz", 0., -100, 100, "cm")
82 B0_TruthTagVz = ROOT.RooRealVar("mcTagVz", "mcTagVz", 0., -100, 100, "cm")
83
84 B0_Z = ROOT.RooRealVar("z", "z", 0., -100, 100, "cm")
85 B0_TruthZ = ROOT.RooRealVar("mcZ", "mcZ", 0., -100, 100, "cm")
86
87 B0_Jpsi_mu0_nPXDHits = ROOT.RooRealVar("Jpsi_mu_0_nPXDHits", "Jpsi_mu_0_nPXDHits", 0., -10., 100.)
88 B0_Jpsi_mu1_nPXDHits = ROOT.RooRealVar("Jpsi_mu_1_nPXDHits", "Jpsi_mu_1_nPXDHits", 0., -10., 100.)
89 B0_Jpsi_mu0_nSVDHits = ROOT.RooRealVar("Jpsi_mu_0_nSVDHits", "Jpsi_mu_0_nSVDHits", 0., -10., 100.)
90 B0_Jpsi_mu1_nSVDHits = ROOT.RooRealVar("Jpsi_mu_1_nSVDHits", "Jpsi_mu_1_nSVDHits", 0., -10., 100.)
91
92 # B0_Jpsi_mu0_firstSVDLayer = ROOT.RooRealVar("B0_Jpsi_mu0__firstSVDLayer", "B0_pi0_e1__firstSVDLayer", 1., -10., 100.)
93 # B0_Jpsi_mu1_firstSVDLayer = ROOT.RooRealVar("B0_Jpsi_mu1__firstSVDLayer", "B0_pi0_e0__firstSVDLayer", 1., -10., 100.)
94
95 DT = ROOT.RooRealVar("DT", "DT", 0., -limDeltaT, limDeltaT, "ps")
96 DSigZ = ROOT.RooRealVar("DSigZ", "DSigZ", 0., -limZSig, limZSig, "cm")
97 # DSigZ = ROOT.RooFormulaVar("DSigZ", "DSigZ", "@@0-@@1", ROOT.RooArgList(B0_Z, B0_TruthZ))
98 DTagZ = ROOT.RooRealVar("DSigZ", "DSigZ", 0., -limZTag, limZTag, "cm")
99 # DTagZ = ROOT.RooFormulaVar("DTagZ", "DTagZ", "@@0-@@1", ROOT.RooArgList(B0_TagVz, B0_TruthTagVz))
100
101 histo_DeltaT = ROOT.TH1F('B0_DeltaT_' + VXDReq, 'Residual of DeltaT',
102 100, -limDeltaT, limDeltaT)
103 histo_DeltaTErr = ROOT.TH1F('B0_DeltaTErr_' + VXDReq, 'Residual of DeltaZsig',
104 100, 0, 30)
105 histo_DeltaZSig = ROOT.TH1F('B0_DeltaZsig_' + VXDReq, 'Residual of DeltaZsig',
106 100, -limZSig, limZSig)
107 histo_DeltaZTag = ROOT.TH1F('B0_DeltaZtag_' + VXDReq, 'Residual of DeltaZsig',
108 100, -limZTag, limZTag)
109
110 cut = "isSignal == 1 " # + "&& abs(B0_DeltaTErr)< " + str(limDeltaTErr) + " "
111
112 if VXDReq == 'PXD1':
113 cut = cut + "&& (Jpsi_mu_0_nPXDHits> 0 || Jpsi_mu_1_nPXDHits> 0) "
114 if VXDReq == 'PXD2':
115 cut = cut + "&& Jpsi_mu_0_nPXDHits> 0 && Jpsi_mu_1_nPXDHits> 0 "
116
117 if VXDReq == 'SVD1':
118 cut = cut + "&& (Jpsi_mu_0_SVDHits> 0 || Jpsi_mu_1_nSVDHits> 0) "
119 if VXDReq == 'SVD2':
120 cut = cut + "&& Jpsi_mu_0_nSVDHits> 0 && Jpsi_mu_1_nSVDHits> 0 "
121
122 argSet = ROOT.RooArgSet(
123 B0_DeltaT,
124 B0_TruthDeltaT,
125 B0_TagVz,
126 B0_TruthTagVz,
127 B0_Z,
128 B0_TruthZ)
129
130 argSet.add(B0_Jpsi_mu0_nPXDHits)
131 argSet.add(B0_Jpsi_mu1_nPXDHits)
132 argSet.add(B0_Jpsi_mu0_nSVDHits)
133 argSet.add(B0_Jpsi_mu1_nSVDHits)
134
135 argSet.add(B0_isSignal)
136
137 data = ROOT.RooDataSet(
138 "data",
139 "data",
140 tdat,
141 argSet,
142 cut)
143
144 if VXDReq == 'PXD1' or VXDReq == 'PXD2':
145 fitDataDTErr = ROOT.RooDataSet("data", "data", tdat, ROOT.RooArgSet(
146 B0_isSignal, B0_Jpsi_mu0_nPXDHits, B0_Jpsi_mu1_nPXDHits, deltaTErr), cut)
147 elif VXDReq == 'SVD1' or VXDReq == 'SVD2':
148 fitDataDTErr = ROOT.RooDataSet("data", "data", tdat, ROOT.RooArgSet(
149 B0_isSignal, B0_Jpsi_mu0_nSVDHits, B0_Jpsi_mu1_nSVDHits, deltaTErr), cut)
150 else:
151 fitDataDTErr = ROOT.RooDataSet("data", "data", tdat, ROOT.RooArgSet(B0_isSignal, deltaTErr), cut)
152
153 # fitData.append(data)
154
155 fitDataDT = ROOT.RooDataSet("fitDataDT", "fitDataDT", ROOT.RooArgSet(DT))
156 fitDataSigZ = ROOT.RooDataSet("fitDataSigZ", "fitDataSigZ", ROOT.RooArgSet(DSigZ))
157 fitDataTagZ = ROOT.RooDataSet("fitDataTagZ", "fitDataTagZ", ROOT.RooArgSet(DTagZ))
158
159 for i in range(data.numEntries()):
160
161 row = data.get(i)
162
163 tDT = row.getRealValue("DeltaT", 0, ROOT.kTRUE) - row.getRealValue("MCDeltaT", 0, ROOT.kTRUE)
164 if abs(tDT) < limDeltaT:
165 DT.setVal(tDT)
166 fitDataDT.add(ROOT.RooArgSet(DT))
167
168 tDSigZ = row.getRealValue("z", 0, ROOT.kTRUE) - row.getRealValue("mcZ", 0, ROOT.kTRUE)
169
170 if abs(tDSigZ) < limZSig:
171 DSigZ.setVal(tDSigZ)
172 fitDataSigZ.add(ROOT.RooArgSet(DSigZ))
173
174 tDTagZ = row.getRealValue("TagVz", 0, ROOT.kTRUE) - row.getRealValue("mcTagVz", 0, ROOT.kTRUE)
175
176 if abs(tDTagZ) < limZTag:
177 DTagZ.setVal(tDTagZ)
178 fitDataTagZ.add(ROOT.RooArgSet(DTagZ))
179
180 fitDataDTErr.Print()
181 fitDataDT.Print()
182 fitDataSigZ.Print()
183 fitDataTagZ.Print()
184 numberOfEntries.append(data.numEntries())
185
186# Fit and plot of the DeltaT Error and DeltaTRECO - DeltaTMC
187
188 Mu1 = ROOT.RooRealVar("Mu1", "Mu1", 0., -limDeltaT, limDeltaT)
189 Mu2 = ROOT.RooRealVar("Mu2", "Mu2", 0., -limDeltaT, limDeltaT)
190 Mu3 = ROOT.RooRealVar("Mu3", "Mu3", 0., -limDeltaT, limDeltaT)
191 Sigma1 = ROOT.RooRealVar("Sigma1", "Sigma1", 1.88046e+00, 0., limDeltaT)
192 Sigma2 = ROOT.RooRealVar("Sigma2", "Sigma2", 3.40331e+00, 0., limDeltaT)
193 Sigma3 = ROOT.RooRealVar("Sigma3", "Sigma3", 8.23171e-01, 0., limDeltaT)
194 frac1 = ROOT.RooRealVar("frac1", "frac1", 5.48703e-01, 0.0, 1.)
195 frac2 = ROOT.RooRealVar("frac2", "frac2", 2.60604e-01, 0.0, 1.)
196
197 g1 = ROOT.RooGaussModel("g1", "g1", DT, Mu1, Sigma1)
198 g2 = ROOT.RooGaussModel("g2", "g2", DT, Mu2, Sigma2)
199 g3 = ROOT.RooGaussModel("g3", "g3", DT, Mu3, Sigma3)
200
201 argset1 = ROOT.RooArgSet(g1)
202 argset2 = ROOT.RooArgSet(g2)
203 argset3 = ROOT.RooArgSet(g3)
204
205 model = ROOT.RooAddModel("model", "model", ROOT.RooArgList(g1, g2, g3), ROOT.RooArgList(frac1, frac2))
206
207 DT.setRange("fitRange", -limDeltaT, limDeltaT)
208
209 fitRes = model.fitTo(
210 fitDataDT,
211 ROOT.RooFit.Minos(ROOT.kFALSE), ROOT.RooFit.Extended(ROOT.kFALSE),
212 ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
213
214 fitRes.Print("v")
215
216 deltaTErr.Print("v")
217
218 resFrame = DT.frame()
219 fitDataDT.plotOn(resFrame)
220
221 meanCBS = ROOT.RooRealVar("meanCBS", "meanCBS", 5.37602e-01, 0.5, 1, "ps")
222 sigmaCBS = ROOT.RooRealVar("sigmaCBS", "sigmaCBS", 8.16334e-02, 0, 0.1, "ps")
223 alphaCBS = ROOT.RooRealVar("alphaCBS", "alphaCBS", -4.85571e-01, -1, 0, "")
224 nCBS = ROOT.RooRealVar("nCBS", "nCBS", 1.86325e+00, 0.5, 3, "")
225 dtErrCBS = ROOT.RooCBShape("dtErrGen", "dtErrGen", deltaTErr, meanCBS, sigmaCBS, alphaCBS, nCBS)
226
227 MuErr1 = ROOT.RooRealVar("MuErr1", "MuErr1", 4.12399e-01, 0.3, 0.6, "ps")
228 SigmaErr1 = ROOT.RooRealVar("SigmaErr1", "SigmaErr1", 5.41152e-02, 0. - 1, 0.07, "ps")
229 gErr1 = ROOT.RooGaussModel("gErr1", "gErr1", deltaTErr, MuErr1, SigmaErr1)
230 fracErr1 = ROOT.RooRealVar("fracErr1", "fracErr1", 7.50810e-01, 0.0, 1.)
231
232 MuErr2 = ROOT.RooRealVar("MuErr2", "MuErr2", 3.26658e-01, 0.2, 0.4, "ps")
233 SigmaErr2 = ROOT.RooRealVar("SigmaErr2", "SigmaErr2", 3.66794e-02, 0.01, 0.08, "ps")
234 gErr2 = ROOT.RooGaussModel("gErr2", "gErr2", deltaTErr, MuErr2, SigmaErr2)
235 fracErr2 = ROOT.RooRealVar("fracErr2", "fracErr2", 1.82254e-01, 0.0, 1.)
236
237 modelTErr = ROOT.RooAddModel(
238 "modelErr", "modelErr", ROOT.RooArgList(
239 dtErrCBS, gErr1, gErr2), ROOT.RooArgList(
240 fracErr1, fracErr2))
241
242 if VXDReq == 'PXD0' or VXDReq == 'PXD1' or VXDReq == 'PXD2':
243 CBSFitRes = modelTErr.fitTo(
244 fitDataDTErr,
245 ROOT.RooFit.Minos(ROOT.kFALSE), ROOT.RooFit.Extended(ROOT.kFALSE),
246 ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
247
248 dtErrCBS.Print()
249
250 argset1TErr = ROOT.RooArgSet(dtErrCBS)
251 argset2TErr = ROOT.RooArgSet(gErr1)
252 argset3TErr = ROOT.RooArgSet(gErr2)
253
254 resFrameDtErr = deltaTErr.frame()
255 fitDataDTErr.plotOn(resFrameDtErr)
256 modelTErr.plotOn(resFrameDtErr)
257 modelTErr.plotOn(
258 resFrameDtErr,
259 ROOT.RooFit.Components(argset1TErr),
260 ROOT.RooFit.LineColor(ROOT.kRed + 2),
261 ROOT.RooFit.LineWidth(4))
262 modelTErr.plotOn(
263 resFrameDtErr,
264 ROOT.RooFit.Components(argset2TErr),
265 ROOT.RooFit.LineColor(ROOT.kGreen + 3),
266 ROOT.RooFit.LineWidth(4))
267 modelTErr.plotOn(
268 resFrameDtErr,
269 ROOT.RooFit.Components(argset3TErr),
270 ROOT.RooFit.LineColor(ROOT.kMagenta + 2),
271 ROOT.RooFit.LineWidth(4))
272
273 f1 = frac1.getVal()
274 f2 = frac2.getVal()
275 f3 = 1 - f1 - f2
276
277 argsetList = []
278 argsetList.append([argset1, f1])
279 argsetList.append([argset2, f2])
280 argsetList.append([argset3, f3])
281
282 argsetList = sorted(argsetList, key=itemgetter(1))
283
284 model.plotOn(resFrame, ROOT.RooFit.LineWidth(3))
285 model.plotOn(resFrame, ROOT.RooFit.Components(argsetList[2][0]), ROOT.RooFit.LineColor(ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
286 model.plotOn(
287 resFrame, ROOT.RooFit.Components(
288 argsetList[1][0]), ROOT.RooFit.LineColor(
289 ROOT.kGreen + 3), ROOT.RooFit.LineWidth(4))
290 model.plotOn(
291 resFrame, ROOT.RooFit.Components(
292 argsetList[0][0]), ROOT.RooFit.LineColor(
293 ROOT.kMagenta + 2), ROOT.RooFit.LineWidth(4))
294
295 resFrame.SetTitle("")
296 sXtitle = "#Deltat - Gen. #Delta t / ps"
297 resFrame.GetXaxis().SetTitle(sXtitle)
298 resFrame.GetXaxis().SetTitleSize(0.05)
299 resFrame.GetXaxis().SetLabelSize(0.045)
300 resFrame.GetYaxis().SetTitleSize(0.05)
301 resFrame.GetYaxis().SetTitleOffset(1.5)
302 resFrame.GetYaxis().SetLabelSize(0.045)
303
304 shift = Mu1.getVal() * f1 + Mu2.getVal() * f2 + Mu3.getVal() * f3
305 resolution = Sigma1.getVal() * f1 + Sigma2.getVal() * f2 + Sigma3.getVal() * f3
306
307 shiftErr = math.sqrt((Mu1.getError() * f1)**2 + (Mu2.getError() * f2)**2 + (Mu3.getError() * f3)**2)
308 resolutionErr = math.sqrt((Sigma1.getError() * f1)**2 + (Sigma2.getError() * f2)**2 + (Sigma3.getError() * f3)**2)
309
310 if shiftErr < 0.01:
311 shiftErr = 0.01
312
313 if resolutionErr < 0.01:
314 resolutionErr = 0.01
315
316 Number = f'{int((f1 + f2) * fitDataDT.numEntries()):d}'
317
318 c1 = ROOT.TCanvas("c1", "c1", 1400, 1100)
319 c1.cd()
320 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
321 Pad.SetLeftMargin(0.15)
322 Pad.SetBottomMargin(0.15)
323 Pad.Draw()
324 Pad.cd()
325 resFrame.Draw()
326 legend = ROOT.TLegend(0.59, 0.6, 0.9, 0.9)
327 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(Number))
328 legend.AddEntry(
329 0,
330 '#splitline{#mu_{#Delta t} =' +
331 f'{shift: 4.2f}' +
332 '}{ #pm ' +
333 f'{shiftErr:4.2f}' +
334 ' ps}')
335 legend.AddEntry(0, '#splitline{#sigma_{#Delta t} =' + f'{resolution: 4.2f}' + '}{ #pm ' +
336 f'{resolutionErr:4.2f}' + ' ps}')
337 legend.SetTextSize(0.054)
338 legend.SetFillColorAlpha(ROOT.kWhite, 0)
339 legend.Draw()
340 Pad.Update()
341 nPlot = PATH + "/test6_CPVResDeltaT" + VXDReq + ".pdf"
342 c1.SaveAs(nPlot)
343 c1.Clear()
344
345 iResult.append(['mu = ' + f'{shift: 4.3f}' + ' +- ' + f'{shiftErr:4.3f}' + ' ps',
346 'sigma =' + f'{resolution: 4.3f}' + ' +- ' + f'{resolutionErr:4.3f}' + ' ps'])
347
348 resFrameDtErr.SetTitle("")
349 sXtitleLandau = "#sigma_{#Deltat} / ps"
350 resFrameDtErr.GetXaxis().SetTitle(sXtitleLandau)
351 resFrameDtErr.GetXaxis().SetTitleSize(0.05)
352 resFrameDtErr.GetXaxis().SetLabelSize(0.045)
353 resFrameDtErr.GetYaxis().SetTitleSize(0.05)
354 resFrameDtErr.GetYaxis().SetTitleOffset(1.5)
355 resFrameDtErr.GetYaxis().SetLabelSize(0.045)
356
357 c1 = ROOT.TCanvas("c1", "c1", 1400, 1100)
358 c1.cd()
359 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
360 Pad.SetLeftMargin(0.15)
361 Pad.SetBottomMargin(0.15)
362 Pad.Draw()
363 Pad.cd()
364 resFrameDtErr.Draw()
365 legend = ROOT.TLegend(0.59, 0.6, 0.9, 0.9)
366 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(Number))
367 legend.AddEntry(0, '#splitline{#mu_{#Delta t} =' + f'{meanCBS.getVal(): 4.2f}' + '}{ #pm ' +
368 f'{meanCBS.getError():4.2f}' + ' ps}') # '{:>6}'.format(Shift)
369 legend.AddEntry(0, '#splitline{#sigma_{#Delta t} =' + f'{sigmaCBS.getVal(): 4.2f}' +
370 '}{ #pm ' + f'{sigmaCBS.getError():4.2f}' + ' ps}') # '{:>4}'.format(Resol)
371 legend.SetTextSize(0.054)
372 legend.SetFillColorAlpha(ROOT.kWhite, 0)
373 # legend.Draw()
374 Pad.Update()
375 nPlot = PATH + "/test6_CPVResDeltaTError" + VXDReq + ".pdf"
376 c1.SaveAs(nPlot)
377 c1.Clear()
378
379 fitRes.Clear()
380 model.Clear()
381
382 # Fit of Delta Z for B0_sig
383
384 Mu1SigZ = ROOT.RooRealVar("Mu1SigZ", "Mu1SigZ", -6.03806e-06, -limZSig, limZSig)
385 Mu2SigZ = ROOT.RooRealVar("Mu2SigZ", "Mu2SigZ", 1.45755e-05, -limZSig, limZSig)
386 Mu3SigZ = ROOT.RooRealVar("Mu3SigZ", "Mu3SigZ", -1.84464e-04, -limZSig, limZSig)
387 Sigma1SigZ = ROOT.RooRealVar("Sigma1SigZ", "Sigma1SigZ", 4.03530e-03, 0., limZSig)
388 Sigma2SigZ = ROOT.RooRealVar("Sigma2SigZ", "Sigma2SigZ", 1.73995e-03, 0., limZSig)
389 Sigma3SigZ = ROOT.RooRealVar("Sigma3SigZ", "Sigma3SigZ", 2.18176e-02, 0., limZSig)
390 frac1SigZ = ROOT.RooRealVar("frac1SigZ", "frac1SigZ", 2.08032e-01, 0.0, 1.)
391 frac2SigZ = ROOT.RooRealVar("frac2SigZ", "frac2SigZ", 7.80053e-01, 0.0, 1.)
392
393 g1SigZ = ROOT.RooGaussian("g1", "g1", DSigZ, Mu1SigZ, Sigma1SigZ)
394 g2SigZ = ROOT.RooGaussian("g2", "g2", DSigZ, Mu2SigZ, Sigma2SigZ)
395 g3SigZ = ROOT.RooGaussian("g3", "g3", DSigZ, Mu3SigZ, Sigma3SigZ)
396
397 argset1SigZ = ROOT.RooArgSet(g1SigZ)
398 argset2SigZ = ROOT.RooArgSet(g2SigZ)
399 argset3SigZ = ROOT.RooArgSet(g3SigZ)
400
401 modelSigZ = ROOT.RooAddPdf(
402 "modelSigZ", "modelSigZ", ROOT.RooArgList(g1SigZ, g2SigZ, g3SigZ),
403 ROOT.RooArgList(frac1SigZ, frac2SigZ))
404
405 DSigZ.setRange("fitRange", -limZSig, limZSig)
406
407 fitResSigZ = modelSigZ.fitTo(
408 fitDataSigZ, ROOT.RooFit.Minos(ROOT.kFALSE),
409 ROOT.RooFit.Extended(ROOT.kFALSE), ROOT.RooFit.NumCPU(1), ROOT.RooFit.Save())
410
411 resFrameSigZ = DSigZ.frame()
412
413 f1SigZ = frac1SigZ.getVal()
414 f2SigZ = frac2SigZ.getVal()
415 f3SigZ = 1 - f1SigZ - f2SigZ
416
417 argsetList = []
418 argsetList.append([argset1SigZ, f1SigZ])
419 argsetList.append([argset2SigZ, f2SigZ])
420 argsetList.append([argset3SigZ, f3SigZ])
421
422 argsetList = sorted(argsetList, key=itemgetter(1))
423
424 fitDataSigZ.plotOn(resFrameSigZ)
425 modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.LineWidth(4))
426 modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.Components(argsetList[2][0]), ROOT.RooFit.LineColor(ROOT.kRed + 2),
427 ROOT.RooFit.LineWidth(4))
428 modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.Components(argsetList[1][0]), ROOT.RooFit.LineColor(ROOT.kGreen + 3),
429 ROOT.RooFit.LineWidth(4))
430 modelSigZ.plotOn(
431 resFrameSigZ,
432 ROOT.RooFit.Components(argsetList[0][0]),
433 ROOT.RooFit.LineColor(
434 ROOT.kMagenta + 2),
435 ROOT.RooFit.LineWidth(4))
436 resFrameSigZ.SetTitle("")
437
438 sXtitleSigZ = "z - mcZ / cm"
439
440 resFrameSigZ.GetXaxis().SetTitle(sXtitleSigZ)
441 resFrameSigZ.GetXaxis().SetLimits(-limZSig, limZSig)
442 resFrameSigZ.GetXaxis().SetTitleSize(0.05)
443 resFrameSigZ.GetXaxis().SetLabelSize(0.045)
444 resFrameSigZ.GetYaxis().SetTitleSize(0.05)
445 resFrameSigZ.GetYaxis().SetTitleOffset(1.5)
446 resFrameSigZ.GetYaxis().SetLabelSize(0.045)
447
448 shiftSigZ = Mu1SigZ.getVal() * f1SigZ + Mu2SigZ.getVal() * f2SigZ + Mu3SigZ.getVal() * f3SigZ
449 resolutionSigZ = Sigma1SigZ.getVal() * f1SigZ + Sigma2SigZ.getVal() * f2SigZ + Sigma3SigZ.getVal() * f3SigZ
450
451 shiftSigZ = shiftSigZ * 10000
452 resolutionSigZ = resolutionSigZ * 10000
453
454 shiftErrSigZ = math.sqrt((Mu1SigZ.getError() * f1SigZ)**2 + (Mu2SigZ.getError() * f2SigZ) ** 2 +
455 (Mu3SigZ.getError() * f3SigZ)**2) * 10000
456 resolutionErrSigZ = math.sqrt((Sigma1SigZ.getError() * f1SigZ)**2 + (Sigma2SigZ.getError() * f2SigZ)**2 +
457 (Sigma3SigZ.getError() * f3SigZ)**2) * 10000
458
459 cSig = ROOT.TCanvas("c1", "c1", 1400, 1100)
460 cSig.cd()
461 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
462 Pad.SetLeftMargin(0.15)
463 Pad.SetBottomMargin(0.15)
464 Pad.Draw()
465 Pad.cd()
466 resFrameSigZ.Draw()
467 legend = ROOT.TLegend(0.63, 0.65, 0.9, 0.9)
468 # NumbrSigZ = '{:d}'.format(int((f1+f2)*fitDataSigZ.numEntries()))
469 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(NumbrSigZ))
470
471 legend.AddEntry(
472 0,
473 '#splitline{#mu_{#Delta z} =' +
474 f'{shiftSigZ: 1.1f}' +
475 '}{ #pm ' +
476 f'{shiftErrSigZ:1.1f}' +
477 ' #mum}')
478 legend.AddEntry(
479 0,
480 '#splitline{#sigma_{#Delta z} =' +
481 f'{resolutionSigZ: 1.1f}' +
482 '}{ #pm ' +
483 f'{resolutionErrSigZ:1.1f}' +
484 ' #mum}')
485
486 legend.SetTextSize(0.05)
487 legend.SetFillColorAlpha(ROOT.kWhite, 0)
488 legend.Draw()
489 Pad.Update()
490 nPlot = PATH + "/test6_CPVResDeltaZsig" + VXDReq + ".pdf"
491 cSig.SaveAs(nPlot)
492 cSig.Clear()
493
494 fitResSigZ.Clear()
495 modelSigZ.Clear()
496
497 iResult.append(['mu = ' + f'{shiftSigZ:^5.3f}' + ' +- ' + f'{shiftErrSigZ:^4.3f}' + ' mum',
498 'sigma = ' + f'{resolutionSigZ:^4.3f}' + ' +- ' + f'{resolutionErrSigZ:^4.3f}' + ' mum'])
499
500 # Fit of Delta z for B0_Tag
501
502 Mu1TagZ = ROOT.RooRealVar("Mu1TagZ", "Mu1TagZ", 0., -limZTag, limZTag)
503 Mu2TagZ = ROOT.RooRealVar("Mu2TagZ", "Mu2TagZ", 0., -limZTag, limZTag)
504 Mu3TagZ = ROOT.RooRealVar("Mu3TagZ", "Mu3TagZ", 0., -limZTag, limZTag)
505 Sigma1TagZ = ROOT.RooRealVar("Sigma1TagZ", "Sigma1TagZ", 2.51877e-02, 0., limZTag)
506 Sigma2TagZ = ROOT.RooRealVar("Sigma2TagZ", "Sigma2TagZ", 1.54011e-02, 0., limZTag)
507 Sigma3TagZ = ROOT.RooRealVar("Sigma3TagZ", "Sigma3TagZ", 1.61081e-02, 0., limZTag)
508 frac1TagZ = ROOT.RooRealVar("frac1TagZ", "frac1TagZ", 1.20825e-01, 0.0, 1.)
509 frac2TagZ = ROOT.RooRealVar("frac2TagZ", "frac2TagZ", 1.10840e-01, 0.0, 1.)
510
511 g1TagZ = ROOT.RooGaussian("g1", "g1", DTagZ, Mu1TagZ, Sigma1TagZ)
512 g2TagZ = ROOT.RooGaussian("g2", "g2", DTagZ, Mu2TagZ, Sigma2TagZ)
513 g3TagZ = ROOT.RooGaussian("g3", "g3", DTagZ, Mu3TagZ, Sigma3TagZ)
514
515 argset1TagZ = ROOT.RooArgSet(g1TagZ)
516 argset2TagZ = ROOT.RooArgSet(g2TagZ)
517 argset3TagZ = ROOT.RooArgSet(g3TagZ)
518
519 modelTagZ = ROOT.RooAddPdf("modelTagZ", "modelTagZ", ROOT.RooArgList(
520 g1TagZ, g2TagZ, g3TagZ), ROOT.RooArgList(frac1TagZ, frac2TagZ))
521
522 DTagZ.setRange("fitRange", -limZTag, limZTag)
523
524 fitResTagZ = modelTagZ.fitTo(
525 fitDataTagZ, ROOT.RooFit.Minos(
526 ROOT.kFALSE), ROOT.RooFit.Extended(
527 ROOT.kFALSE), ROOT.RooFit.NumCPU(1), ROOT.RooFit.Save())
528
529 fitResTagZ.Print("v")
530 resFrameTagZ = DTagZ.frame()
531
532 f1TagZ = frac1TagZ.getVal()
533 f2TagZ = frac2TagZ.getVal()
534 f3TagZ = 1 - f1TagZ - f2TagZ
535
536 argsetList = []
537 argsetList.append([argset1TagZ, f1TagZ])
538 argsetList.append([argset2TagZ, f2TagZ])
539 argsetList.append([argset3TagZ, f3TagZ])
540
541 argsetList = sorted(argsetList, key=itemgetter(1))
542
543 fitDataTagZ.plotOn(resFrameTagZ)
544 modelTagZ.plotOn(resFrameTagZ, ROOT.RooFit.LineWidth(4))
545 modelTagZ.plotOn(
546 resFrameTagZ, ROOT.RooFit.Components(
547 argsetList[2][0]), ROOT.RooFit.LineColor(
548 ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
549 modelTagZ.plotOn(
550 resFrameTagZ,
551 ROOT.RooFit.Components(argsetList[1][0]),
552 ROOT.RooFit.LineColor(ROOT.kGreen + 3),
553 ROOT.RooFit.LineWidth(4))
554 modelTagZ.plotOn(
555 resFrameTagZ,
556 ROOT.RooFit.Components(argsetList[0][0]),
557 ROOT.RooFit.LineColor(ROOT.kMagenta + 2),
558 ROOT.RooFit.LineWidth(4))
559 resFrameTagZ.SetTitle("")
560
561 sXtitleTagZ = "TagVz - mcTagVz / cm"
562
563 resFrameTagZ.GetXaxis().SetTitle(sXtitleTagZ)
564 resFrameTagZ.GetXaxis().SetTitleSize(0.05)
565 resFrameTagZ.GetXaxis().SetLabelSize(0.045)
566 resFrameTagZ.GetYaxis().SetTitleSize(0.05)
567 resFrameTagZ.GetYaxis().SetTitleOffset(1.5)
568 resFrameTagZ.GetYaxis().SetLabelSize(0.045)
569
570 shiftTagZ = Mu1TagZ.getVal() * f1TagZ + Mu2TagZ.getVal() * f2TagZ + Mu3TagZ.getVal() * f3TagZ
571 resolutionTagZ = Sigma1TagZ.getVal() * f1TagZ + Sigma2TagZ.getVal() * f2TagZ + Sigma3TagZ.getVal() * f3TagZ
572
573 shiftTagZ = shiftTagZ * 10000
574 resolutionTagZ = resolutionTagZ * 10000
575
576 shiftErrTagZ = math.sqrt((Mu1TagZ.getError() * f1TagZ)**2 + (Mu2TagZ.getError() * f2TagZ) ** 2 +
577 (Mu3TagZ.getError() * f3TagZ)**2) * 10000
578 resolutionErrTagZ = math.sqrt((Sigma1TagZ.getError() * f1TagZ)**2 + (Sigma2TagZ.getError() * f2TagZ)**2 +
579 (Sigma3TagZ.getError() * f3TagZ)**2) * 10000
580
581 cTag = ROOT.TCanvas("c1", "c1", 1400, 1100)
582 cTag.cd()
583 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
584 Pad.SetLeftMargin(0.15)
585 Pad.SetBottomMargin(0.15)
586 Pad.Draw()
587 Pad.cd()
588 resFrameTagZ.Draw()
589 legend = ROOT.TLegend(0.64, 0.65, 0.9, 0.9)
590 # NumbrTagZ = '{:d}'.format(int((f1+f2)*fitDataTagZ.numEntries()))
591 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(NumbrTagZ))
592
593 legend.AddEntry(0, '#splitline{#mu_{#Delta z} =' + f'{shiftTagZ: 1.1f}' +
594 '}{ #pm ' + f'{shiftErrTagZ: 1.1f}' + ' #mum}')
595 legend.AddEntry(0, '#splitline{#sigma_{#Delta z} =' + f'{resolutionTagZ: 1.1f}' +
596 '}{ #pm ' + f'{resolutionErrTagZ: 1.1f}' + ' #mum}')
597 legend.SetTextSize(0.05)
598 legend.SetFillColorAlpha(ROOT.kWhite, 0)
599 legend.Draw()
600 Pad.Update()
601 nPlot = PATH + "/test6_CPVResDeltaZtag" + VXDReq + ".pdf"
602 cTag.SaveAs(nPlot)
603 cTag.Clear()
604
605 iResult.append(['mu = ' + f'{shiftTagZ:^5.3f}' + ' +- ' + f'{shiftErrTagZ:^4.3f}' + ' mum',
606 'sigma = ' + f'{resolutionTagZ:^4.3f}' + ' +- ' + f'{resolutionErrTagZ:^4.3f}' + ' mum'])
607
608 fitResults.append(iResult)
609
610
611print('*********************** FIT RESULTS ***************************')
612print('* *')
613print('* WEIGHTED AVERAGES OF THE PARAMETERS *')
614print('* OF THE FITTED 3 GAUSSIAN *')
615print('* *')
616print('**************** WITHOUT VXD HIT REQUIREMENT ******************')
617print('* *')
618print('* DeltaT - Gen. DeltaT *')
619print('* *')
620print('* ' + fitResults[0][0][0] + ' ' + fitResults[0][0][1] + ' *')
621print('* *')
622print('* SigZ - Gen. SigZ *')
623print('* *')
624print('* ' + fitResults[0][1][0] + ' ' + fitResults[0][1][1] + ' *')
625print('* *')
626print('* TagZ - Gen. TagZ *')
627print('* *')
628print('* ' + fitResults[0][2][0] + ' ' + fitResults[0][2][1] + ' *')
629print('* *')
630print('********REQUIRING BOTH MUON TRACKS TO HAVE A ' + sys.argv[3] + ' HIT***********')
631print('* *')
632print('* Efficiency *')
633print('* *')
634print('* N_' + VXDReqs[1] + '/N_' + VXDReqs[0] + ' = ' + str(numberOfEntries[1]) + "/" + str(numberOfEntries[0]) + ' = ' +
635 f'{float(numberOfEntries[1] / numberOfEntries[0] * 100):^3.2f}' + '% *')
636print('* *')
637print('* DeltaT - Gen. DeltaT *')
638print('* *')
639print('* ' + fitResults[1][0][0] + ' ' + fitResults[1][0][1] + ' *')
640print('* *')
641print('* SigZ - Gen. SigZ *')
642print('* *')
643print('* ' + fitResults[1][1][0] + ' ' + fitResults[1][1][1] + ' *')
644print('* *')
645print('* TagZ - Gen. TagZ *')
646print('* *')
647print('* ' + fitResults[1][2][0] + ' ' + fitResults[1][2][1] + ' *')
648print('* *')
649print('***************************************************************')