Belle II Software  release-08-01-10
deltaTVertexTagVResolution.py
1 #!/usr/bin/env python
2 
3 
10 
11 
25 
26 
27 from basf2 import B2FATAL
28 import ROOT
29 import sys
30 import glob
31 import math
32 from operator import itemgetter
33 
34 if 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 
39 PATH = "."
40 
41 workingFile = sys.argv[1]
42 workingFiles = glob.glob(str(workingFile))
43 
44 limDeltaT = 5
45 limDeltaTErr = 3.0
46 limZSig = 0.03
47 limZTag = 0.03
48 
49 treeName = str(sys.argv[2])
50 
51 # No PXD hit equired: PXD0. At least one PXD (SVD) hit for one of the muon tracks: PXD1 (SVD1).
52 # Hit required for both muon tracks: PXD2 (SVD2)"
53 VXDReqs = []
54 if str(sys.argv[3]) == "PXD":
55  VXDReqs = ["PXD0", "PXD2"]
56 elif str(sys.argv[3]) == "SVD":
57  VXDReqs = ["SVD0", "SVD2"]
58 else:
59  B2FATAL('Not available VXD requirement " + str(sys.argv[3]) + ". Available are "PXD" and "SVD".')
60 
61 ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
62 
63 fitResults = []
64 numberOfEntries = []
65 
66 for 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  Numbr = '{:d}'.format(int((f1 + f2) * fitDataDT.numEntries()))
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(Numbr))
328  legend.AddEntry(
329  0,
330  '#splitline{#mu_{#Delta t} =' +
331  '{: 4.2f}'.format(shift) +
332  '}{ #pm ' +
333  '{:4.2f}'.format(shiftErr) +
334  ' ps}')
335  legend.AddEntry(0, '#splitline{#sigma_{#Delta t} =' + '{: 4.2f}'.format(resolution) + '}{ #pm ' +
336  '{:4.2f}'.format(resolutionErr) + ' 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 = ' + '{: 4.3f}'.format(shift) + ' +- ' + '{:4.3f}'.format(shiftErr) + ' ps',
346  'sigma =' + '{: 4.3f}'.format(resolution) + ' +- ' + '{:4.3f}'.format(resolutionErr) + ' 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(Numbr))
367  legend.AddEntry(0, '#splitline{#mu_{#Delta t} =' + '{: 4.2f}'.format(meanCBS.getVal()) + '}{ #pm ' +
368  '{:4.2f}'.format(meanCBS.getError()) + ' ps}') # '{:>6}'.format(Shift)
369  legend.AddEntry(0, '#splitline{#sigma_{#Delta t} =' + '{: 4.2f}'.format(sigmaCBS.getVal()) +
370  '}{ #pm ' + '{:4.2f}'.format(sigmaCBS.getError()) + ' 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  '{: 1.1f}'.format(shiftSigZ) +
475  '}{ #pm ' +
476  '{:1.1f}'.format(shiftErrSigZ) +
477  ' #mum}')
478  legend.AddEntry(
479  0,
480  '#splitline{#sigma_{#Delta z} =' +
481  '{: 1.1f}'.format(resolutionSigZ) +
482  '}{ #pm ' +
483  '{:1.1f}'.format(resolutionErrSigZ) +
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 = ' + '{:^5.3f}'.format(shiftSigZ) + ' +- ' + '{:^4.3f}'.format(shiftErrSigZ) + ' mum',
498  'sigma = ' + '{:^4.3f}'.format(resolutionSigZ) + ' +- ' + '{:^4.3f}'.format(resolutionErrSigZ) + ' 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} =' + '{: 1.1f}'.format(shiftTagZ) +
594  '}{ #pm ' + '{: 1.1f}'.format(shiftErrTagZ) + ' #mum}')
595  legend.AddEntry(0, '#splitline{#sigma_{#Delta z} =' + '{: 1.1f}'.format(resolutionTagZ) +
596  '}{ #pm ' + '{: 1.1f}'.format(resolutionErrTagZ) + ' #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 = ' + '{:^5.3f}'.format(shiftTagZ) + ' +- ' + '{:^4.3f}'.format(shiftErrTagZ) + ' mum',
606  'sigma = ' + '{:^4.3f}'.format(resolutionTagZ) + ' +- ' + '{:^4.3f}'.format(resolutionErrTagZ) + ' mum'])
607 
608  fitResults.append(iResult)
609 
610 
611 print('*********************** FIT RESULTS ***************************')
612 print('* *')
613 print('* WEIGHTED AVERAGES OF THE PARAMETERS *')
614 print('* OF THE FITTED 3 GAUSSIAN *')
615 print('* *')
616 print('**************** WITHOUT VXD HIT REQUIREMENT ******************')
617 print('* *')
618 print('* DeltaT - Gen. DeltaT *')
619 print('* *')
620 print('* ' + fitResults[0][0][0] + ' ' + fitResults[0][0][1] + ' *')
621 print('* *')
622 print('* SigZ - Gen. SigZ *')
623 print('* *')
624 print('* ' + fitResults[0][1][0] + ' ' + fitResults[0][1][1] + ' *')
625 print('* *')
626 print('* TagZ - Gen. TagZ *')
627 print('* *')
628 print('* ' + fitResults[0][2][0] + ' ' + fitResults[0][2][1] + ' *')
629 print('* *')
630 print('********REQUIRING BOTH MUON TRACKS TO HAVE A ' + sys.argv[3] + ' HIT***********')
631 print('* *')
632 print('* Efficiency *')
633 print('* *')
634 print('* N_' + VXDReqs[1] + '/N_' + VXDReqs[0] + ' = ' + str(numberOfEntries[1]) + "/" + str(numberOfEntries[0]) + ' = ' +
635  '{:^3.2f}'.format(float((numberOfEntries[1] / numberOfEntries[0]) * 100)) + '% *')
636 print('* *')
637 print('* DeltaT - Gen. DeltaT *')
638 print('* *')
639 print('* ' + fitResults[1][0][0] + ' ' + fitResults[1][0][1] + ' *')
640 print('* *')
641 print('* SigZ - Gen. SigZ *')
642 print('* *')
643 print('* ' + fitResults[1][1][0] + ' ' + fitResults[1][1][1] + ' *')
644 print('* *')
645 print('* TagZ - Gen. TagZ *')
646 print('* *')
647 print('* ' + fitResults[1][2][0] + ' ' + fitResults[1][2][1] + ' *')
648 print('* *')
649 print('***************************************************************')