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