Belle II Software  release-08-01-10
test6_CPVResolutionVertexTagV.py
1 #!/usr/bin/env python
2 
3 
10 
11 """
12 <header>
13  <input>CPVToolsOutput.root</input>
14  <output>test6_CPVResolutionVertexTagV.root</output>
15  <contact>Yo Sato; yosato@post.kek.jp</contact>
16  <description>This validation script performs a fit of the values of DeltaT, DeltaTErr, Dz for B0_sig and Deltaz for B0_tag.
17  The DeltaT and DeltaZ distributions are fitted with 3 Gaussian functions.
18  DeltaTErr is fitted with a CBShape function and two Gaussians.
19  PXD requirement PXD0 means no hit is required.
20  PXD requirement PXD2 means both muon tracks from B->J/Psi Ks are required.</description>
21 </header>
22 """
23 
24 import ROOT
25 import math
26 import array
27 from operator import itemgetter
28 
29 
30 PATH = "."
31 
32 workingFiles = ["../CPVToolsOutput.root"]
33 
34 limDeltaT = 5
35 limDeltaTErr = 3.0
36 limZSig = 0.03
37 limZTag = 0.03
38 
39 treename = "B0tree"
40 
41 # Output Validation file
42 outputFile = ROOT.TFile("test6_CPVResolutionVertexTagV.root", "RECREATE")
43 
44 # Values to be watched
45 outputNtuple = ROOT.TNtuple(
46  "Vertex_TagV_Resols",
47  "Weighted averages of the resolution params of Vertex and TagV Tools",
48  "mu_DT_PXD0:sigma_DT_PXD0:mu_DZsig_PXD0:sigma_DZsig_PXD0:mu_DZtag_PXD0:sigma_DZtag_PXD0:" +
49  "mu_DT_PXD2:sigma_DT_PXD2:mu_DZsig_PXD2:sigma_DZsig_PXD2:mu_DZtag_PXD2:sigma_DZtag_PXD2:PXD2_PXD0_Eff")
50 
51 outputNtuple.SetAlias('Description', "These are the weighted averages of the mean and the standard deviation " +
52  "of the residuals for DeltaT, DeltaZsig and DeltaZtag. The fit is performed with 3 Gaussian functions. " +
53  "The units are ps for DeltaT and microns for DeltaZ.")
54 outputNtuple.SetAlias(
55  'Check',
56  "These parameters should not change drastically. Since the nightly reconstruction validation runs " +
57  "on the same input file (which changes only from release to release), the values between builds should be the same.")
58 outputNtuple.SetAlias('Contact', "yosato@post.kek.jp")
59 
60 # No PXD hit required: PXD0. At least one PXD (SVD) hit for one of the muon tracks: PXD1 (SVD1).
61 # Hit required for both muon tracks: PXD2 (SVD2)"
62 VXDReqs = ["PXD0", "PXD2"]
63 
64 
65 fitResults = []
66 fitResultsForNtuple = []
67 numberOfEntries = []
68 
69 for VXDReq in VXDReqs:
70  iResult = []
71 
72  tdat = ROOT.TChain(treename)
73 
74  for iFile in workingFiles:
75  tdat.AddFile(iFile)
76 
77  evt_no = ROOT.RooRealVar("evt_no", "evt_no", 0., -10., 10000000.)
78 
79  B0_DeltaT = ROOT.RooRealVar("DeltaT", "DeltaT", 0.)
80  deltaTErr = ROOT.RooRealVar("DeltaTErr", "DeltaTErr", 0, limDeltaTErr, "ps")
81  B0_TruthDeltaT = ROOT.RooRealVar("mcDeltaT", "mcDeltaT", 0.)
82 
83  B0_qrMC = ROOT.RooRealVar("qrMC", "qrMC", 0., -100, 100)
84  B0_isSignal = ROOT.RooRealVar("isSignal", "isSignal", 0., 0., 512.)
85 
86  B0_TagVz = ROOT.RooRealVar("TagVz", "TagVz", 0., -100, 100, "cm")
87  B0_TruthTagVz = ROOT.RooRealVar("mcTagVz", "mcTagVz", 0., -100, 100, "cm")
88 
89  B0_Z = ROOT.RooRealVar("z", "z", 0., -100, 100, "cm")
90  B0_TruthZ = ROOT.RooRealVar("mcDecayVertexZ", "mcDecayVertexZ", 0., -100, 100, "cm")
91 
92  B0_Jpsi_mu0_nPXDHits = ROOT.RooRealVar("Jpsi_mu_0_nPXDHits", "Jpsi_mu_0_nPXDHits", 0., -10., 100.)
93  B0_Jpsi_mu1_nPXDHits = ROOT.RooRealVar("Jpsi_mu_1_nPXDHits", "Jpsi_mu_1_nPXDHits", 0., -10., 100.)
94  B0_Jpsi_mu0_nSVDHits = ROOT.RooRealVar("Jpsi_mu_0_nSVDHits", "Jpsi_mu_0_nSVDHits", 0., -10., 100.)
95  B0_Jpsi_mu1_nSVDHits = ROOT.RooRealVar("Jpsi_mu_1_nSVDHits", "Jpsi_mu_1_nSVDHits", 0., -10., 100.)
96 
97  # B0_Jpsi_mu0_firstSVDLayer = ROOT.RooRealVar("B0_Jpsi_mu0__firstSVDLayer", "B0_pi0_e1__firstSVDLayer", 1., -10., 100.)
98  # B0_Jpsi_mu1_firstSVDLayer = ROOT.RooRealVar("B0_Jpsi_mu1__firstSVDLayer", "B0_pi0_e0__firstSVDLayer", 1., -10., 100.)
99 
100  DT = ROOT.RooRealVar("DT", "DT", 0., -limDeltaT, limDeltaT, "ps")
101  DSigZ = ROOT.RooRealVar("DSigZ", "DSigZ", 0., -limZSig, limZSig, "cm")
102  # DSigZ = ROOT.RooFormulaVar("DSigZ", "DSigZ", "@@0-@@1", ROOT.RooArgList(B0_Z, B0_TruthZ))
103  DTagZ = ROOT.RooRealVar("DSigZ", "DSigZ", 0., -limZTag, limZTag, "cm")
104  # DTagZ = ROOT.RooFormulaVar("DTagZ", "DTagZ", "@@0-@@1", ROOT.RooArgList(B0_TagVz, B0_TruthTagVz))
105 
106  histo_DeltaT = ROOT.TH1F('B0_DeltaT_' + VXDReq, 'Residual of DeltaT',
107  100, -limDeltaT, limDeltaT)
108  histo_DeltaTErr = ROOT.TH1F('B0_DeltaTErr_' + VXDReq, 'Residual of DeltaZsig',
109  100, 0, 30)
110  histo_DeltaZSig = ROOT.TH1F('B0_DeltaZsig_' + VXDReq, 'Residual of DeltaZsig',
111  100, -limZSig, limZSig)
112  histo_DeltaZTag = ROOT.TH1F('B0_DeltaZtag_' + VXDReq, 'Residual of DeltaZsig',
113  100, -limZTag, limZTag)
114 
115  cut = "isSignal == 1 " # + "&& abs(B0_DeltaTErr)< " + str(limDeltaTErr) + " "
116 
117  if VXDReq == 'PXD1':
118  cut = cut + "&& (Jpsi_mu_0_nPXDHits> 0 || Jpsi_mu_1_nPXDHits> 0) "
119  if VXDReq == 'PXD2':
120  cut = cut + "&& Jpsi_mu_0_nPXDHits> 0 && Jpsi_mu_1_nPXDHits> 0 "
121 
122  if VXDReq == 'SVD1':
123  cut = cut + "&& (Jpsi_mu_0_SVDHits> 0 || Jpsi_mu_1_nSVDHits> 0) "
124  if VXDReq == 'SVD2':
125  cut = cut + "&& Jpsi_mu_0_nSVDHits> 0 && Jpsi_mu_1_nSVDHits> 0 "
126 
127  ROOT.gROOT.SetBatch(True)
128  c1 = ROOT.TCanvas("c1", "c1", 1400, 1100)
129 
130  tdat.Draw("DeltaT - mcDeltaT >> B0_DeltaT_" + VXDReq, cut)
131  tdat.Draw("DeltaTErr >> B0_DeltaTErr_" + VXDReq, cut)
132  tdat.Draw("z - mcDecayVertexZ >> B0_DeltaZsig_" + VXDReq, cut)
133  tdat.Draw("TagVz - mcTagVz >> B0_DeltaZtag_" + VXDReq, cut)
134 
135  # Validation Plot 1
136  histo_DeltaT.GetXaxis().SetLabelSize(0.04)
137  histo_DeltaT.GetYaxis().SetLabelSize(0.04)
138  histo_DeltaT.GetYaxis().SetTitleOffset(0.7)
139  histo_DeltaT.GetXaxis().SetTitleOffset(0.7)
140  histo_DeltaT.GetXaxis().SetTitleSize(0.06)
141  histo_DeltaT.GetYaxis().SetTitleSize(0.07)
142  histo_DeltaT.SetTitle('DeltaT Residual for PXD requirement ' + VXDReq + '; #Deltat - Gen. #Deltat / ps ; Events')
143 
144  histo_DeltaT.GetListOfFunctions().Add(
145  ROOT.TNamed('Description', 'DeltaT Residual for PXD requirement ' + VXDReq +
146  '. PXD0 means no PXD hit required. PXD2 means both muon tracks are required to have a PXD hit.'))
147  histo_DeltaT.GetListOfFunctions().Add(ROOT.TNamed('Check', 'Std. Dev. and Mean should not change drastically.'))
148  histo_DeltaT.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
149  histo_DeltaT.Write()
150 
151  # Validation Plot 2
152  histo_DeltaTErr.GetXaxis().SetLabelSize(0.04)
153  histo_DeltaTErr.GetYaxis().SetLabelSize(0.04)
154  histo_DeltaTErr.GetYaxis().SetTitleOffset(0.7)
155  histo_DeltaTErr.GetXaxis().SetTitleOffset(0.7)
156  histo_DeltaTErr.GetXaxis().SetTitleSize(0.06)
157  histo_DeltaTErr.GetYaxis().SetTitleSize(0.07)
158  histo_DeltaTErr.SetTitle('DeltaT error for PXD requirement ' + VXDReq + ' ; #sigma_{#Deltat} / ps ; Events')
159 
160  histo_DeltaTErr.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'logy'))
161  histo_DeltaTErr.GetListOfFunctions().Add(
162  ROOT.TNamed('Description', 'DeltaT error for PXD requirement ' + VXDReq +
163  '. PXD0 means no PXD hit required. PXD2 means both muon tracks are required to have a PXD hit.'))
164  histo_DeltaTErr.GetListOfFunctions().Add(
165  ROOT.TNamed(
166  'Check',
167  'Std. Dev. and Mean should not change drastically. Peaks after 2.6 ps should not increase.'))
168  histo_DeltaTErr.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
169  histo_DeltaTErr.Write()
170 
171  # Validation Plot 3
172  histo_DeltaZSig.GetXaxis().SetLabelSize(0.04)
173  histo_DeltaZSig.GetYaxis().SetLabelSize(0.04)
174  histo_DeltaZSig.GetYaxis().SetTitleOffset(0.7)
175  histo_DeltaZSig.GetXaxis().SetTitleOffset(0.7)
176  histo_DeltaZSig.GetXaxis().SetTitleSize(0.06)
177  histo_DeltaZSig.GetYaxis().SetTitleSize(0.07)
178  histo_DeltaZSig.SetTitle('DeltaZ Residual on signal side for requirement ' + VXDReq + '; B0_Z - Gen. B0_Z / cm ; Events')
179 
180  histo_DeltaZSig.GetListOfFunctions().Add(
181  ROOT.TNamed('Description', 'DeltaZ Residual on signal side for PXD requirement ' + VXDReq +
182  '. PXD0 means no PXD hit required. PXD2 means both muon tracks are required to have a PXD hit.'))
183  histo_DeltaZSig.GetListOfFunctions().Add(ROOT.TNamed('Check', 'Std. Dev. and Mean should not change drastically.'))
184  histo_DeltaZSig.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
185  histo_DeltaZSig.Write()
186 
187  # Validation Plot 4
188  histo_DeltaZTag.GetXaxis().SetLabelSize(0.04)
189  histo_DeltaZTag.GetYaxis().SetLabelSize(0.04)
190  histo_DeltaZTag.GetYaxis().SetTitleOffset(0.7)
191  histo_DeltaZTag.GetXaxis().SetTitleOffset(0.7)
192  histo_DeltaZTag.GetXaxis().SetTitleSize(0.06)
193  histo_DeltaZTag.GetYaxis().SetTitleSize(0.07)
194  histo_DeltaZTag.SetTitle('DeltaZ Residual on tag side for requirement ' + VXDReq + '; B0_TagVz - Gen. B0_TagVz / cm; Events')
195 
196  histo_DeltaZTag.GetListOfFunctions().Add(
197  ROOT.TNamed('Description', 'DeltaZ Residual on tag side for PXD requirement ' + VXDReq +
198  '. PXD0 means no PXD hit required. PXD2 means both muon tracks are required to have a PXD hit.'))
199  histo_DeltaZTag.GetListOfFunctions().Add(ROOT.TNamed('Check', 'Std. Dev. and Mean should not change drastically.'))
200  histo_DeltaZTag.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
201  histo_DeltaZTag.Write()
202 
203  c1.Clear()
204 
205  argSet = ROOT.RooArgSet(
206  # B0_mcTagPDG,
207  B0_DeltaT,
208  B0_TruthDeltaT,
209  # B0_mcErrors,
210  # evt_no,
211  B0_TagVz,
212  B0_TruthTagVz,
213  B0_Z,
214  B0_TruthZ)
215 
216  argSet.add(B0_Jpsi_mu0_nPXDHits)
217  argSet.add(B0_Jpsi_mu1_nPXDHits)
218  argSet.add(B0_Jpsi_mu0_nSVDHits)
219  argSet.add(B0_Jpsi_mu1_nSVDHits)
220 
221  argSet.add(B0_isSignal)
222 
223  # argSet.add(B0_X)
224  # argSet.add(B0_TruthX)
225  # argSet.add(B0_Y)
226  # argSet.add(B0_TruthY)
227 
228  # argSet.add(B0_TagVx)
229  # argSet.add(B0_TruthTagVx)
230  # argSet.add(B0_TagVy)
231  # argSet.add(B0_TruthTagVy)
232 
233  # argSet.add(deltaTErr)
234  # argSet.add(B0_qrMC)
235 
236  data = ROOT.RooDataSet(
237  "data",
238  "data",
239  tdat,
240  argSet,
241  cut)
242 
243  if VXDReq == 'PXD1' or VXDReq == 'PXD2':
244  fitDataDTErr = ROOT.RooDataSet("data", "data", tdat, ROOT.RooArgSet(
245  B0_isSignal, B0_Jpsi_mu0_nPXDHits, B0_Jpsi_mu1_nPXDHits, deltaTErr),
246  f'{cut} && DeltaTErr >= {deltaTErr.getMin()} && DeltaTErr <= {deltaTErr.getMax()}')
247  elif VXDReq == 'SVD1' or VXDReq == 'SVD2':
248  fitDataDTErr = ROOT.RooDataSet("data", "data", tdat, ROOT.RooArgSet(
249  B0_isSignal, B0_Jpsi_mu0_nSVDHits, B0_Jpsi_mu1_nSVDHits, deltaTErr),
250  f'{cut} && DeltaTErr >= {deltaTErr.getMin()} && DeltaTErr <= {deltaTErr.getMax()}')
251  else:
252  fitDataDTErr = ROOT.RooDataSet("data", "data", tdat, ROOT.RooArgSet(B0_isSignal, deltaTErr),
253  f'{cut} && DeltaTErr >= {deltaTErr.getMin()} && DeltaTErr <= {deltaTErr.getMax()}')
254 
255  # fitData.append(data)
256 
257  fitDataDT = ROOT.RooDataSet("fitDataDT", "fitDataDT", ROOT.RooArgSet(DT))
258  fitDataSigZ = ROOT.RooDataSet("fitDataSigZ", "fitDataSigZ", ROOT.RooArgSet(DSigZ))
259  fitDataTagZ = ROOT.RooDataSet("fitDataTagZ", "fitDataTagZ", ROOT.RooArgSet(DTagZ))
260 
261  for i in range(data.numEntries()):
262 
263  row = data.get(i)
264 
265  tDT = row.getRealValue("DeltaT", 0, ROOT.kTRUE) - row.getRealValue("mcDeltaT", 0, ROOT.kTRUE)
266  if abs(tDT) < limDeltaT:
267  DT.setVal(tDT)
268  fitDataDT.add(ROOT.RooArgSet(DT))
269 
270  tDSigZ = row.getRealValue("z", 0, ROOT.kTRUE) - row.getRealValue("mcDecayVertexZ", 0, ROOT.kTRUE)
271 
272  if abs(tDSigZ) < limZSig:
273  DSigZ.setVal(tDSigZ)
274  fitDataSigZ.add(ROOT.RooArgSet(DSigZ))
275 
276  tDTagZ = row.getRealValue("TagVz", 0, ROOT.kTRUE) - row.getRealValue("mcTagVz", 0, ROOT.kTRUE)
277 
278  if abs(tDTagZ) < limZTag:
279  DTagZ.setVal(tDTagZ)
280  fitDataTagZ.add(ROOT.RooArgSet(DTagZ))
281 
282  fitDataDTErr.Print()
283  fitDataDT.Print()
284  fitDataSigZ.Print()
285  fitDataTagZ.Print()
286  numberOfEntries.append(data.numEntries())
287 
288  # Fit and plot of the DeltaT Error and DeltaTRECO - DeltaTMC
289 
290  Mu1 = ROOT.RooRealVar("Mu1", "Mu1", 0., -limDeltaT, limDeltaT)
291  Mu2 = ROOT.RooRealVar("Mu2", "Mu2", 0., -limDeltaT, limDeltaT)
292  Mu3 = ROOT.RooRealVar("Mu3", "Mu3", 0., -limDeltaT, limDeltaT)
293  Sigma1 = ROOT.RooRealVar("Sigma1", "Sigma1", 0.4, 0., limDeltaT)
294  Sigma2 = ROOT.RooRealVar("Sigma2", "Sigma2", 2.4, 0., limDeltaT)
295  Sigma3 = ROOT.RooRealVar("Sigma3", "Sigma3", 0.8, 0., limDeltaT)
296  frac1 = ROOT.RooRealVar("frac1", "frac1", 0.6, 0.0, 1.)
297  frac2 = ROOT.RooRealVar("frac2", "frac2", 0.1, 0.0, 1.)
298 
299  g1 = ROOT.RooGaussModel("g1", "g1", DT, Mu1, Sigma1)
300  g2 = ROOT.RooGaussModel("g2", "g2", DT, Mu2, Sigma2)
301  g3 = ROOT.RooGaussModel("g3", "g3", DT, Mu3, Sigma3)
302 
303  argset1 = ROOT.RooArgSet(g1)
304  argset2 = ROOT.RooArgSet(g2)
305  argset3 = ROOT.RooArgSet(g3)
306 
307  model = ROOT.RooAddModel("model", "model", ROOT.RooArgList(g1, g2, g3), ROOT.RooArgList(frac1, frac2))
308 
309  DT.setRange("fitRange", -limDeltaT, limDeltaT)
310 
311  fitRes = model.fitTo(fitDataDT, ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
312 
313  fitRes.Print("v")
314 
315  deltaTErr.Print("v")
316 
317  resFrame = DT.frame()
318  fitDataDT.plotOn(resFrame)
319 
320  meanCBS = ROOT.RooRealVar("meanCBS", "meanCBS", 0.4, 0.3, 1, "ps")
321  sigmaCBS = ROOT.RooRealVar("sigmaCBS", "sigmaCBS", 0.03, 0.01, 0.1, "ps")
322  alphaCBS = ROOT.RooRealVar("alphaCBS", "alphaCBS", -0.4, -1, 0, "")
323  nCBS = ROOT.RooRealVar("nCBS", "nCBS", 3, 0.5, 5, "")
324  dtErrCBS = ROOT.RooCBShape("dtErrGen", "dtErrGen", deltaTErr, meanCBS, sigmaCBS, alphaCBS, nCBS)
325 
326  MuErr1 = ROOT.RooRealVar("MuErr1", "MuErr1", 0.3, 0.2, 0.6, "ps")
327  SigmaErr1 = ROOT.RooRealVar("SigmaErr1", "SigmaErr1", 0.03, 0.01, 0.07, "ps")
328  gErr1 = ROOT.RooGaussModel("gErr1", "gErr1", deltaTErr, MuErr1, SigmaErr1)
329  fracErr1 = ROOT.RooRealVar("fracErr1", "fracErr1", 0.45, 0.0, 0.7)
330 
331  MuErr2 = ROOT.RooRealVar("MuErr2", "MuErr2", 0.24, 0.2, 0.4, "ps")
332  SigmaErr2 = ROOT.RooRealVar("SigmaErr2", "SigmaErr2", 0.03, 0.01, 0.08, "ps")
333  gErr2 = ROOT.RooGaussModel("gErr2", "gErr2", deltaTErr, MuErr2, SigmaErr2)
334  fracErr2 = ROOT.RooRealVar("fracErr2", "fracErr2", 0.2, 0.0, 0.5)
335 
336  modelTErr = ROOT.RooAddModel(
337  "modelErr", "modelErr", ROOT.RooArgList(
338  dtErrCBS, gErr1, gErr2), ROOT.RooArgList(
339  fracErr1, fracErr2))
340 
341  if VXDReq == 'PXD0' or VXDReq == 'PXD1' or VXDReq == 'PXD2':
342  CBSFitRes = modelTErr.fitTo(fitDataDTErr, ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
343  CBSFitRes.Print("v")
344 
345  dtErrCBS.Print()
346 
347  argset1TErr = ROOT.RooArgSet(dtErrCBS)
348  argset2TErr = ROOT.RooArgSet(gErr1)
349  argset3TErr = ROOT.RooArgSet(gErr2)
350 
351  resFrameDtErr = deltaTErr.frame()
352  fitDataDTErr.plotOn(resFrameDtErr)
353  modelTErr.plotOn(resFrameDtErr)
354  modelTErr.plotOn(
355  resFrameDtErr,
356  ROOT.RooFit.Components(argset1TErr),
357  ROOT.RooFit.LineColor(ROOT.kRed + 2),
358  ROOT.RooFit.LineWidth(4))
359  modelTErr.plotOn(
360  resFrameDtErr,
361  ROOT.RooFit.Components(argset2TErr),
362  ROOT.RooFit.LineColor(ROOT.kGreen + 3),
363  ROOT.RooFit.LineWidth(4))
364  modelTErr.plotOn(
365  resFrameDtErr,
366  ROOT.RooFit.Components(argset3TErr),
367  ROOT.RooFit.LineColor(ROOT.kMagenta + 2),
368  ROOT.RooFit.LineWidth(4))
369 
370  f1 = frac1.getVal()
371  f2 = frac2.getVal()
372  f3 = 1 - f1 - f2
373 
374  argsetList = []
375  argsetList.append([argset1, f1])
376  argsetList.append([argset2, f2])
377  argsetList.append([argset3, f3])
378 
379  argsetList = sorted(argsetList, key=itemgetter(1))
380 
381  model.plotOn(resFrame, ROOT.RooFit.LineWidth(3))
382  model.plotOn(resFrame, ROOT.RooFit.Components(argsetList[2][0]), ROOT.RooFit.LineColor(ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
383  model.plotOn(
384  resFrame, ROOT.RooFit.Components(
385  argsetList[1][0]), ROOT.RooFit.LineColor(
386  ROOT.kGreen + 3), ROOT.RooFit.LineWidth(4))
387  model.plotOn(
388  resFrame, ROOT.RooFit.Components(
389  argsetList[0][0]), ROOT.RooFit.LineColor(
390  ROOT.kMagenta + 2), ROOT.RooFit.LineWidth(4))
391 
392  resFrame.SetTitle("")
393  sXtitle = "#Deltat - Gen. #Delta t / ps"
394  resFrame.GetXaxis().SetTitle(sXtitle)
395  resFrame.GetXaxis().SetTitleSize(0.05)
396  resFrame.GetXaxis().SetLabelSize(0.045)
397  resFrame.GetYaxis().SetTitleSize(0.05)
398  resFrame.GetYaxis().SetTitleOffset(1.5)
399  resFrame.GetYaxis().SetLabelSize(0.045)
400 
401  shift = Mu1.getVal() * f1 + Mu2.getVal() * f2 + Mu3.getVal() * f3
402  resolution = Sigma1.getVal() * f1 + Sigma2.getVal() * f2 + Sigma3.getVal() * f3
403 
404  shiftErr = math.sqrt((Mu1.getError() * f1)**2 + (Mu2.getError() * f2)**2 + (Mu3.getError() * f3)**2)
405  resolutionErr = math.sqrt((Sigma1.getError() * f1)**2 + (Sigma2.getError() * f2)**2 + (Sigma3.getError() * f3)**2)
406 
407  if shiftErr < 0.01:
408  shiftErr = 0.01
409 
410  if resolutionErr < 0.01:
411  resolutionErr = 0.01
412 
413  Numbr = '{:d}'.format(int((f1 + f2) * fitDataDT.numEntries()))
414 
415  c1.cd()
416  Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
417  Pad.SetLeftMargin(0.15)
418  Pad.SetBottomMargin(0.15)
419  Pad.Draw()
420  Pad.cd()
421  resFrame.Draw()
422  legend = ROOT.TLegend(0.59, 0.6, 0.9, 0.9)
423  # legend.AddEntry(0, 'Entries' + '{:>11}'.format(Numbr))
424  legend.AddEntry(0, f'#splitline{{#mu_{{#Delta t}} = {shift:4.2f}}}{{ #pm {shiftErr:4.2f} ps}}')
425  legend.AddEntry(0, f'#splitline{{#sigma_{{#Delta t}} = {resolution:4.2f}}}{{ #pm {resolutionErr:4.2f} ps}}')
426  legend.SetTextSize(0.054)
427  legend.SetFillColorAlpha(ROOT.kWhite, 0)
428  legend.Draw()
429  Pad.Update()
430  nPlot = PATH + "/test6_CPVResDeltaT" + VXDReq + ".pdf"
431  c1.SaveAs(nPlot)
432  c1.Clear()
433 
434  iResult.append([f'mu = {shift:4.2f} +- {shiftErr:4.2f} ps',
435  f'sigma = {resolution:4.2f} +- {resolutionErr:4.2f} ps'])
436  fitResultsForNtuple.append(shift)
437  fitResultsForNtuple.append(resolution)
438 
439  resFrameDtErr.SetTitle("")
440  sXtitleLandau = "#sigma_{#Deltat} / ps"
441  resFrameDtErr.GetXaxis().SetTitle(sXtitleLandau)
442  resFrameDtErr.GetXaxis().SetTitleSize(0.05)
443  resFrameDtErr.GetXaxis().SetLabelSize(0.045)
444  resFrameDtErr.GetYaxis().SetTitleSize(0.05)
445  resFrameDtErr.GetYaxis().SetTitleOffset(1.5)
446  resFrameDtErr.GetYaxis().SetLabelSize(0.045)
447 
448  c1.cd()
449  Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
450  Pad.SetLeftMargin(0.15)
451  Pad.SetBottomMargin(0.15)
452  Pad.Draw()
453  Pad.cd()
454  resFrameDtErr.Draw()
455  legend = ROOT.TLegend(0.59, 0.6, 0.9, 0.9)
456  # legend.AddEntry(0, 'Entries' + '{:>11}'.format(Numbr))
457  legend.AddEntry(0, f'#splitline{{#mu_{{#Delta t}} = {meanCBS.getVal():4.2f}}}{{ #pm '
458  f'{meanCBS.getError():4.2f} ps}}') # '{:>6}'.format(Shift)
459  legend.AddEntry(0, f'#splitline{{#sigma_{{#Delta t}} = {sigmaCBS.getVal():4.2f}}}'
460  f'{{ #pm {sigmaCBS.getError():4.2f} ps}}') # '{:>4}'.format(Resol)
461  legend.SetTextSize(0.054)
462  legend.SetFillColorAlpha(ROOT.kWhite, 0)
463  # legend.Draw()
464  Pad.Update()
465  nPlot = PATH + "/test6_CPVResDeltaTError" + VXDReq + ".pdf"
466  c1.SaveAs(nPlot)
467  c1.Clear()
468 
469  fitRes.Clear()
470  model.Clear()
471 
472  # Fit of Delta Z for B0_sig
473 
474  Mu1SigZ = ROOT.RooRealVar("Mu1SigZ", "Mu1SigZ", -9e-06, -limZSig, limZSig)
475  Mu2SigZ = ROOT.RooRealVar("Mu2SigZ", "Mu2SigZ", -1.8e-05, -limZSig, limZSig)
476  Mu3SigZ = ROOT.RooRealVar("Mu3SigZ", "Mu3SigZ", 6e-05, -limZSig, limZSig)
477  Sigma1SigZ = ROOT.RooRealVar("Sigma1SigZ", "Sigma1SigZ", 3e-03, 1e-6, limZSig)
478  Sigma2SigZ = ROOT.RooRealVar("Sigma2SigZ", "Sigma2SigZ", 1.4e-03, 1e-6, limZSig)
479  Sigma3SigZ = ROOT.RooRealVar("Sigma3SigZ", "Sigma3SigZ", 0.01, 1e-6, limZSig)
480  frac1SigZ = ROOT.RooRealVar("frac1SigZ", "frac1SigZ", 0.2, 0.0, 0.5)
481  frac2SigZ = ROOT.RooRealVar("frac2SigZ", "frac2SigZ", 0.75, 0.5, 1.)
482 
483  g1SigZ = ROOT.RooGaussian("g1", "g1", DSigZ, Mu1SigZ, Sigma1SigZ)
484  g2SigZ = ROOT.RooGaussian("g2", "g2", DSigZ, Mu2SigZ, Sigma2SigZ)
485  g3SigZ = ROOT.RooGaussian("g3", "g3", DSigZ, Mu3SigZ, Sigma3SigZ)
486 
487  argset1SigZ = ROOT.RooArgSet(g1SigZ)
488  argset2SigZ = ROOT.RooArgSet(g2SigZ)
489  argset3SigZ = ROOT.RooArgSet(g3SigZ)
490 
491  modelSigZ = ROOT.RooAddPdf(
492  "modelSigZ", "modelSigZ", ROOT.RooArgList(g1SigZ, g2SigZ, g3SigZ),
493  ROOT.RooArgList(frac1SigZ, frac2SigZ))
494 
495  DSigZ.setRange("fitRange", -limZSig, limZSig)
496 
497  fitResSigZ = modelSigZ.fitTo(fitDataSigZ, ROOT.RooFit.NumCPU(1), ROOT.RooFit.Save())
498 
499  fitResSigZ.Print("v")
500 
501  resFrameSigZ = DSigZ.frame()
502 
503  f1SigZ = frac1SigZ.getVal()
504  f2SigZ = frac2SigZ.getVal()
505  f3SigZ = 1 - f1SigZ - f2SigZ
506 
507  argsetList = []
508  argsetList.append([argset1SigZ, f1SigZ])
509  argsetList.append([argset2SigZ, f2SigZ])
510  argsetList.append([argset3SigZ, f3SigZ])
511 
512  argsetList = sorted(argsetList, key=itemgetter(1))
513 
514  fitDataSigZ.plotOn(resFrameSigZ)
515  modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.LineWidth(4))
516  modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.Components(argsetList[2][0]), ROOT.RooFit.LineColor(ROOT.kRed + 2),
517  ROOT.RooFit.LineWidth(4))
518  modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.Components(argsetList[1][0]), ROOT.RooFit.LineColor(ROOT.kGreen + 3),
519  ROOT.RooFit.LineWidth(4))
520  modelSigZ.plotOn(
521  resFrameSigZ,
522  ROOT.RooFit.Components(argsetList[0][0]),
523  ROOT.RooFit.LineColor(
524  ROOT.kMagenta + 2),
525  ROOT.RooFit.LineWidth(4))
526  resFrameSigZ.SetTitle("")
527 
528  sXtitleSigZ = "B0_Z - Gen. B0_Z / cm"
529 
530  resFrameSigZ.GetXaxis().SetTitle(sXtitleSigZ)
531  resFrameSigZ.GetXaxis().SetLimits(-limZSig, limZSig)
532  resFrameSigZ.GetXaxis().SetTitleSize(0.05)
533  resFrameSigZ.GetXaxis().SetLabelSize(0.045)
534  resFrameSigZ.GetYaxis().SetTitleSize(0.05)
535  resFrameSigZ.GetYaxis().SetTitleOffset(1.5)
536  resFrameSigZ.GetYaxis().SetLabelSize(0.045)
537 
538  shiftSigZ = Mu1SigZ.getVal() * f1SigZ + Mu2SigZ.getVal() * f2SigZ + Mu3SigZ.getVal() * f3SigZ
539  resolutionSigZ = Sigma1SigZ.getVal() * f1SigZ + Sigma2SigZ.getVal() * f2SigZ + Sigma3SigZ.getVal() * f3SigZ
540 
541  shiftSigZ = shiftSigZ * 10000
542  resolutionSigZ = resolutionSigZ * 10000
543 
544  shiftErrSigZ = math.sqrt((Mu1SigZ.getError() * f1SigZ)**2 + (Mu2SigZ.getError() * f2SigZ) ** 2 +
545  (Mu3SigZ.getError() * f3SigZ)**2) * 10000
546  resolutionErrSigZ = math.sqrt((Sigma1SigZ.getError() * f1SigZ)**2 + (Sigma2SigZ.getError() * f2SigZ)**2 +
547  (Sigma3SigZ.getError() * f3SigZ)**2) * 10000
548 
549  c1.cd()
550  Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
551  Pad.SetLeftMargin(0.15)
552  Pad.SetBottomMargin(0.15)
553  Pad.Draw()
554  Pad.cd()
555  resFrameSigZ.Draw()
556  legend = ROOT.TLegend(0.63, 0.65, 0.9, 0.9)
557  # NumbrSigZ = '{:d}'.format(int((f1+f2)*fitDataSigZ.numEntries()))
558  # legend.AddEntry(0, 'Entries' + '{:>11}'.format(NumbrSigZ))
559 
560  legend.AddEntry(
561  0,
562  f'#splitline{{#mu_{{#Delta z}} = {shiftSigZ:1.1f}}}'
563  f'{{ #pm {shiftErrSigZ:1.1f} #mum}}')
564  legend.AddEntry(
565  0,
566  f'#splitline{{#sigma_{{#Delta z}} = {resolutionSigZ:1.1f}}}'
567  f'{{ #pm {resolutionErrSigZ:1.1f} #mum}}')
568 
569  legend.SetTextSize(0.05)
570  legend.SetFillColorAlpha(ROOT.kWhite, 0)
571  legend.Draw()
572  Pad.Update()
573  nPlot = PATH + "/test6_CPVResDeltaZsig" + VXDReq + ".pdf"
574  c1.SaveAs(nPlot)
575  c1.Clear()
576 
577  fitResSigZ.Clear()
578  modelSigZ.Clear()
579 
580  iResult.append([f'mu = {shiftSigZ:4.1f} +- {shiftErrSigZ:3.1f} mum',
581  f'sigma = {resolutionSigZ:4.1f} +- {resolutionErrSigZ:3.1f} mum'])
582  fitResultsForNtuple.append(shiftSigZ)
583  fitResultsForNtuple.append(resolutionSigZ)
584 
585  # Fit of Delta z for B0_Tag
586 
587  Mu1TagZ = ROOT.RooRealVar("Mu1TagZ", "Mu1TagZ", 0., -limZTag, limZTag)
588  Mu2TagZ = ROOT.RooRealVar("Mu2TagZ", "Mu2TagZ", 0., -limZTag, limZTag)
589  Mu3TagZ = ROOT.RooRealVar("Mu3TagZ", "Mu3TagZ", 0., -limZTag, limZTag)
590  Sigma1TagZ = ROOT.RooRealVar("Sigma1TagZ", "Sigma1TagZ", 2.51877e-02, 1e-6, limZTag)
591  Sigma2TagZ = ROOT.RooRealVar("Sigma2TagZ", "Sigma2TagZ", 1.54011e-02, 1e-6, limZTag)
592  Sigma3TagZ = ROOT.RooRealVar("Sigma3TagZ", "Sigma3TagZ", 1.61081e-02, 1e-6, limZTag)
593  frac1TagZ = ROOT.RooRealVar("frac1TagZ", "frac1TagZ", 1.20825e-01, 0.0, 1.)
594  frac2TagZ = ROOT.RooRealVar("frac2TagZ", "frac2TagZ", 1.10840e-01, 0.0, 1.)
595 
596  g1TagZ = ROOT.RooGaussian("g1", "g1", DTagZ, Mu1TagZ, Sigma1TagZ)
597  g2TagZ = ROOT.RooGaussian("g2", "g2", DTagZ, Mu2TagZ, Sigma2TagZ)
598  g3TagZ = ROOT.RooGaussian("g3", "g3", DTagZ, Mu3TagZ, Sigma3TagZ)
599 
600  argset1TagZ = ROOT.RooArgSet(g1TagZ)
601  argset2TagZ = ROOT.RooArgSet(g2TagZ)
602  argset3TagZ = ROOT.RooArgSet(g3TagZ)
603 
604  modelTagZ = ROOT.RooAddPdf("modelTagZ", "modelTagZ", ROOT.RooArgList(
605  g1TagZ, g2TagZ, g3TagZ), ROOT.RooArgList(frac1TagZ, frac2TagZ))
606 
607  DTagZ.setRange("fitRange", -limZTag, limZTag)
608 
609  fitResTagZ = modelTagZ.fitTo(fitDataTagZ, ROOT.RooFit.NumCPU(1), ROOT.RooFit.Save())
610  fitResTagZ.Print("v")
611 
612  resFrameTagZ = DTagZ.frame()
613 
614  f1TagZ = frac1TagZ.getVal()
615  f2TagZ = frac2TagZ.getVal()
616  f3TagZ = 1 - f1TagZ - f2TagZ
617 
618  argsetList = []
619  argsetList.append([argset1TagZ, f1TagZ])
620  argsetList.append([argset2TagZ, f2TagZ])
621  argsetList.append([argset3TagZ, f3TagZ])
622 
623  argsetList = sorted(argsetList, key=itemgetter(1))
624 
625  fitDataTagZ.plotOn(resFrameTagZ)
626  modelTagZ.plotOn(resFrameTagZ, ROOT.RooFit.LineWidth(4))
627  modelTagZ.plotOn(
628  resFrameTagZ, ROOT.RooFit.Components(
629  argsetList[2][0]), ROOT.RooFit.LineColor(
630  ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
631  modelTagZ.plotOn(
632  resFrameTagZ,
633  ROOT.RooFit.Components(argsetList[1][0]),
634  ROOT.RooFit.LineColor(ROOT.kGreen + 3),
635  ROOT.RooFit.LineWidth(4))
636  modelTagZ.plotOn(
637  resFrameTagZ,
638  ROOT.RooFit.Components(argsetList[0][0]),
639  ROOT.RooFit.LineColor(ROOT.kMagenta + 2),
640  ROOT.RooFit.LineWidth(4))
641  resFrameTagZ.SetTitle("")
642 
643  sXtitleTagZ = "B0_TagVz - Gen. B0_TagVz / cm"
644 
645  resFrameTagZ.GetXaxis().SetTitle(sXtitleTagZ)
646  resFrameTagZ.GetXaxis().SetTitleSize(0.05)
647  resFrameTagZ.GetXaxis().SetLabelSize(0.045)
648  resFrameTagZ.GetYaxis().SetTitleSize(0.05)
649  resFrameTagZ.GetYaxis().SetTitleOffset(1.5)
650  resFrameTagZ.GetYaxis().SetLabelSize(0.045)
651 
652  shiftTagZ = Mu1TagZ.getVal() * f1TagZ + Mu2TagZ.getVal() * f2TagZ + Mu3TagZ.getVal() * f3TagZ
653  resolutionTagZ = Sigma1TagZ.getVal() * f1TagZ + Sigma2TagZ.getVal() * f2TagZ + Sigma3TagZ.getVal() * f3TagZ
654 
655  shiftTagZ = shiftTagZ * 10000
656  resolutionTagZ = resolutionTagZ * 10000
657 
658  shiftErrTagZ = math.sqrt((Mu1TagZ.getError() * f1TagZ)**2 + (Mu2TagZ.getError() * f2TagZ) ** 2 +
659  (Mu3TagZ.getError() * f3TagZ)**2) * 10000
660  resolutionErrTagZ = math.sqrt((Sigma1TagZ.getError() * f1TagZ)**2 + (Sigma2TagZ.getError() * f2TagZ)**2 +
661  (Sigma3TagZ.getError() * f3TagZ)**2) * 10000
662 
663  c1.cd()
664  Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
665  Pad.SetLeftMargin(0.15)
666  Pad.SetBottomMargin(0.15)
667  Pad.Draw()
668  Pad.cd()
669  resFrameTagZ.Draw()
670  legend = ROOT.TLegend(0.64, 0.65, 0.9, 0.9)
671  # NumbrTagZ = '{:d}'.format(int((f1+f2)*fitDataTagZ.numEntries()))
672  # legend.AddEntry(0, 'Entries' + '{:>11}'.format(NumbrTagZ))
673 
674  legend.AddEntry(0, f'#splitline{{#mu_{{#Delta z}} = {shiftTagZ:1.1f}'
675  f'}}{{ #pm {shiftErrTagZ:1.1f} #mum}}')
676  legend.AddEntry(0, f'#splitline{{#sigma_{{#Delta z}} = {resolutionTagZ:1.1f}'
677  f'}}{{ #pm {resolutionErrTagZ:1.1f} #mum}}')
678  legend.SetTextSize(0.05)
679  legend.SetFillColorAlpha(ROOT.kWhite, 0)
680  legend.Draw()
681  Pad.Update()
682  nPlot = PATH + f"/test6_CPVResDeltaZtag{VXDReq}.pdf"
683  c1.SaveAs(nPlot)
684  c1.Clear()
685 
686  iResult.append([f'mu = {shiftTagZ:4.1f} +- {shiftErrTagZ:3.1f} mum',
687  f'sigma = {resolutionTagZ:4.1f} +- {resolutionErrTagZ:3.1f} mum'])
688  fitResultsForNtuple.append(shiftTagZ)
689  fitResultsForNtuple.append(resolutionTagZ)
690 
691  fitResults.append(iResult)
692 
693 fitResultsForNtuple.append(float((numberOfEntries[1] / numberOfEntries[0]) * 100))
694 outputNtuple.Fill(array.array('f', fitResultsForNtuple))
695 outputNtuple.Write()
696 outputFile.Close()
697 
698 print('*********************** FIT RESULTS ***************************')
699 print('* *')
700 print('* WEIGHTED AVERAGES OF THE PARAMETERS *')
701 print('* OF THE FITTED 3 GAUSSIAN *')
702 print('* *')
703 print('**************** WITHOUT PXD HIT REQUIREMENT ******************')
704 print('* *')
705 print('* DeltaT - Gen. DeltaT *')
706 print('* *')
707 print('* ' + fitResults[0][0][0] + ' ' + fitResults[0][0][1] + ' *')
708 print('* *')
709 print('* SigZ - Gen. SigZ *')
710 print('* *')
711 print('* ' + fitResults[0][1][0] + ' ' + fitResults[0][1][1] + ' *')
712 print('* *')
713 print('* TagZ - Gen. TagZ *')
714 print('* *')
715 print('* ' + fitResults[0][2][0] + ' ' + fitResults[0][2][1] + ' *')
716 print('* *')
717 print('********REQUIRING BOTH MUON TRACKS TO HAVE A PXD HIT***********')
718 print('* *')
719 print('* Efficiency *')
720 print('* *')
721 print(f'* N_{VXDReqs[1]}/N_{VXDReqs[0]} = {numberOfEntries[1]}/{numberOfEntries[0]} = '
722  f'{(numberOfEntries[1] / numberOfEntries[0]) * 100:3.2f}% *')
723 print('* *')
724 print('* DeltaT - Gen. DeltaT *')
725 print('* *')
726 print('* ' + fitResults[1][0][0] + ' ' + fitResults[1][0][1] + ' *')
727 print('* *')
728 print('* SigZ - Gen. SigZ *')
729 print('* *')
730 print('* ' + fitResults[1][1][0] + ' ' + fitResults[1][1][1] + ' *')
731 print('* *')
732 print('* TagZ - Gen. TagZ *')
733 print('* *')
734 print('* ' + fitResults[1][2][0] + ' ' + fitResults[1][2][1] + ' *')
735 print('* *')
736 print('***************************************************************')