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