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