Belle II Software light-2603-ina
test6_CPVResolutionVertexTagV.py
1#!/usr/bin/env python3
2
3
10
11"""
12<header>
13 <input>CPVToolsOutput.root</input>
14 <output>test6_CPVResolutionVertexTagV.root</output>
15 <contact>Paul Feichtinger; paul.feichtinger@ijs.si</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
24import ROOT
25import math
26import array
27from operator import itemgetter
28
29
30PATH = "."
31
32workingFiles = ["../CPVToolsOutput.root"]
33
34limDeltaT = 5
35limDeltaTErr = 3.0
36limZSig = 0.03
37limZTag = 0.03
38
39treename = "B0tree"
40
41# Output Validation file
42outputFile = ROOT.TFile("test6_CPVResolutionVertexTagV.root", "RECREATE")
43
44# Values to be watched
45outputNtuple = 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
51outputNtuple.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.")
54outputNtuple.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.")
58outputNtuple.SetAlias('Contact', "paul.feichtinger@ijs.si")
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)"
62VXDReqs = ["PXD0", "PXD2"]
63
64
65fitResults = []
66fitResultsForNtuple = []
67numberOfEntries = []
68
69for 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(f"c1_{VXDReq}", "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', 'paul.feichtinger@ijs.si'))
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', 'paul.feichtinger@ijs.si'))
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', 'paul.feichtinger@ijs.si'))
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', 'paul.feichtinger@ijs.si'))
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 argSet,
240 Import=tdat,
241 Cut=cut)
242
243 if VXDReq == 'PXD1' or VXDReq == 'PXD2':
244 fitDataDTErr = ROOT.RooDataSet(
245 "data",
246 "data",
247 ROOT.RooArgSet(B0_isSignal, B0_Jpsi_mu0_nPXDHits, B0_Jpsi_mu1_nPXDHits, deltaTErr),
248 Import=tdat,
249 Cut=f'{cut} && DeltaTErr >= {deltaTErr.getMin()} && DeltaTErr <= {deltaTErr.getMax()}')
250 elif VXDReq == 'SVD1' or VXDReq == 'SVD2':
251 fitDataDTErr = ROOT.RooDataSet(
252 "data",
253 "data",
254 ROOT.RooArgSet(B0_isSignal, B0_Jpsi_mu0_nSVDHits, B0_Jpsi_mu1_nSVDHits, deltaTErr),
255 Import=tdat,
256 Cut=f'{cut} && DeltaTErr >= {deltaTErr.getMin()} && DeltaTErr <= {deltaTErr.getMax()}')
257 else:
258 fitDataDTErr = ROOT.RooDataSet(
259 "data",
260 "data",
261 ROOT.RooArgSet(B0_isSignal, deltaTErr),
262 Import=tdat,
263 Cut=f'{cut} && DeltaTErr >= {deltaTErr.getMin()} && DeltaTErr <= {deltaTErr.getMax()}')
264
265 # fitData.append(data)
266
267 fitDataDT = ROOT.RooDataSet("fitDataDT", "fitDataDT", ROOT.RooArgSet(DT))
268 fitDataSigZ = ROOT.RooDataSet("fitDataSigZ", "fitDataSigZ", ROOT.RooArgSet(DSigZ))
269 fitDataTagZ = ROOT.RooDataSet("fitDataTagZ", "fitDataTagZ", ROOT.RooArgSet(DTagZ))
270
271 for i in range(data.numEntries()):
272
273 row = data.get(i)
274
275 tDT = row.getRealValue("DeltaT", 0, ROOT.kTRUE) - row.getRealValue("mcDeltaT", 0, ROOT.kTRUE)
276 if abs(tDT) < limDeltaT:
277 DT.setVal(tDT)
278 fitDataDT.add(ROOT.RooArgSet(DT))
279
280 tDSigZ = row.getRealValue("z", 0, ROOT.kTRUE) - row.getRealValue("mcDecayVertexZ", 0, ROOT.kTRUE)
281
282 if abs(tDSigZ) < limZSig:
283 DSigZ.setVal(tDSigZ)
284 fitDataSigZ.add(ROOT.RooArgSet(DSigZ))
285
286 tDTagZ = row.getRealValue("TagVz", 0, ROOT.kTRUE) - row.getRealValue("mcTagVz", 0, ROOT.kTRUE)
287
288 if abs(tDTagZ) < limZTag:
289 DTagZ.setVal(tDTagZ)
290 fitDataTagZ.add(ROOT.RooArgSet(DTagZ))
291
292 fitDataDTErr.Print()
293 fitDataDT.Print()
294 fitDataSigZ.Print()
295 fitDataTagZ.Print()
296 numberOfEntries.append(data.numEntries())
297
298 # Fit and plot of the DeltaT Error and DeltaTRECO - DeltaTMC
299
300 Mu1 = ROOT.RooRealVar("Mu1", "Mu1", 0., -limDeltaT, limDeltaT)
301 Mu2 = ROOT.RooRealVar("Mu2", "Mu2", 0., -limDeltaT, limDeltaT)
302 Mu3 = ROOT.RooRealVar("Mu3", "Mu3", 0., -limDeltaT, limDeltaT)
303 Sigma1 = ROOT.RooRealVar("Sigma1", "Sigma1", 0.4, 0., limDeltaT)
304 Sigma2 = ROOT.RooRealVar("Sigma2", "Sigma2", 2.4, 0., limDeltaT)
305 Sigma3 = ROOT.RooRealVar("Sigma3", "Sigma3", 0.8, 0., limDeltaT)
306 frac1 = ROOT.RooRealVar("frac1", "frac1", 0.6, 0.0, 1.)
307 frac2 = ROOT.RooRealVar("frac2", "frac2", 0.1, 0.0, 1.)
308
309 g1 = ROOT.RooGaussModel("g1", "g1", DT, Mu1, Sigma1)
310 g2 = ROOT.RooGaussModel("g2", "g2", DT, Mu2, Sigma2)
311 g3 = ROOT.RooGaussModel("g3", "g3", DT, Mu3, Sigma3)
312
313 argset1 = ROOT.RooArgSet(g1)
314 argset2 = ROOT.RooArgSet(g2)
315 argset3 = ROOT.RooArgSet(g3)
316
317 model = ROOT.RooAddModel("model", "model", ROOT.RooArgList(g1, g2, g3), ROOT.RooArgList(frac1, frac2))
318
319 DT.setRange("fitRange", -limDeltaT, limDeltaT)
320
321 fitRes = model.fitTo(fitDataDT, ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
322
323 fitRes.Print("v")
324
325 deltaTErr.Print("v")
326
327 resFrame = DT.frame()
328 fitDataDT.plotOn(resFrame)
329
330 meanCBS = ROOT.RooRealVar("meanCBS", "meanCBS", 0.4, 0.3, 1, "ps")
331 sigmaCBS = ROOT.RooRealVar("sigmaCBS", "sigmaCBS", 0.03, 0.01, 0.1, "ps")
332 alphaCBS = ROOT.RooRealVar("alphaCBS", "alphaCBS", -0.4, -1, 0, "")
333 nCBS = ROOT.RooRealVar("nCBS", "nCBS", 3, 0.5, 5, "")
334 dtErrCBS = ROOT.RooCBShape("dtErrGen", "dtErrGen", deltaTErr, meanCBS, sigmaCBS, alphaCBS, nCBS)
335
336 MuErr1 = ROOT.RooRealVar("MuErr1", "MuErr1", 0.3, 0.2, 0.6, "ps")
337 SigmaErr1 = ROOT.RooRealVar("SigmaErr1", "SigmaErr1", 0.03, 0.01, 0.07, "ps")
338 gErr1 = ROOT.RooGaussModel("gErr1", "gErr1", deltaTErr, MuErr1, SigmaErr1)
339 fracErr1 = ROOT.RooRealVar("fracErr1", "fracErr1", 0.45, 0.0, 0.7)
340
341 MuErr2 = ROOT.RooRealVar("MuErr2", "MuErr2", 0.24, 0.2, 0.4, "ps")
342 SigmaErr2 = ROOT.RooRealVar("SigmaErr2", "SigmaErr2", 0.03, 0.01, 0.08, "ps")
343 gErr2 = ROOT.RooGaussModel("gErr2", "gErr2", deltaTErr, MuErr2, SigmaErr2)
344 fracErr2 = ROOT.RooRealVar("fracErr2", "fracErr2", 0.2, 0.0, 0.5)
345
346 modelTErr = ROOT.RooAddModel(
347 "modelErr", "modelErr", ROOT.RooArgList(
348 dtErrCBS, gErr1, gErr2), ROOT.RooArgList(
349 fracErr1, fracErr2))
350
351 if VXDReq == 'PXD0' or VXDReq == 'PXD1' or VXDReq == 'PXD2':
352 CBSFitRes = modelTErr.fitTo(fitDataDTErr, ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
353 CBSFitRes.Print("v")
354
355 dtErrCBS.Print()
356
357 argset1TErr = ROOT.RooArgSet(dtErrCBS)
358 argset2TErr = ROOT.RooArgSet(gErr1)
359 argset3TErr = ROOT.RooArgSet(gErr2)
360
361 resFrameDtErr = deltaTErr.frame()
362 fitDataDTErr.plotOn(resFrameDtErr)
363 modelTErr.plotOn(resFrameDtErr)
364 modelTErr.plotOn(
365 resFrameDtErr,
366 ROOT.RooFit.Components(argset1TErr),
367 ROOT.RooFit.LineColor(ROOT.kRed + 2),
368 ROOT.RooFit.LineWidth(4))
369 modelTErr.plotOn(
370 resFrameDtErr,
371 ROOT.RooFit.Components(argset2TErr),
372 ROOT.RooFit.LineColor(ROOT.kGreen + 3),
373 ROOT.RooFit.LineWidth(4))
374 modelTErr.plotOn(
375 resFrameDtErr,
376 ROOT.RooFit.Components(argset3TErr),
377 ROOT.RooFit.LineColor(ROOT.kMagenta + 2),
378 ROOT.RooFit.LineWidth(4))
379
380 f1 = frac1.getVal()
381 f2 = frac2.getVal()
382 f3 = 1 - f1 - f2
383
384 argsetList = []
385 argsetList.append([argset1, f1])
386 argsetList.append([argset2, f2])
387 argsetList.append([argset3, f3])
388
389 argsetList = sorted(argsetList, key=itemgetter(1))
390
391 model.plotOn(resFrame, ROOT.RooFit.LineWidth(3))
392 model.plotOn(resFrame, ROOT.RooFit.Components(argsetList[2][0]), ROOT.RooFit.LineColor(ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
393 model.plotOn(
394 resFrame, ROOT.RooFit.Components(
395 argsetList[1][0]), ROOT.RooFit.LineColor(
396 ROOT.kGreen + 3), ROOT.RooFit.LineWidth(4))
397 model.plotOn(
398 resFrame, ROOT.RooFit.Components(
399 argsetList[0][0]), ROOT.RooFit.LineColor(
400 ROOT.kMagenta + 2), ROOT.RooFit.LineWidth(4))
401
402 resFrame.SetTitle("")
403 sXtitle = "#Deltat - Gen. #Delta t / ps"
404 resFrame.GetXaxis().SetTitle(sXtitle)
405 resFrame.GetXaxis().SetTitleSize(0.05)
406 resFrame.GetXaxis().SetLabelSize(0.045)
407 resFrame.GetYaxis().SetTitleSize(0.05)
408 resFrame.GetYaxis().SetTitleOffset(1.5)
409 resFrame.GetYaxis().SetLabelSize(0.045)
410
411 shift = Mu1.getVal() * f1 + Mu2.getVal() * f2 + Mu3.getVal() * f3
412 resolution = Sigma1.getVal() * f1 + Sigma2.getVal() * f2 + Sigma3.getVal() * f3
413
414 shiftErr = math.sqrt((Mu1.getError() * f1)**2 + (Mu2.getError() * f2)**2 + (Mu3.getError() * f3)**2)
415 resolutionErr = math.sqrt((Sigma1.getError() * f1)**2 + (Sigma2.getError() * f2)**2 + (Sigma3.getError() * f3)**2)
416
417 if shiftErr < 0.01:
418 shiftErr = 0.01
419
420 if resolutionErr < 0.01:
421 resolutionErr = 0.01
422
423 Number = f'{int((f1 + f2) * fitDataDT.numEntries()):d}'
424
425 c1.cd()
426 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
427 Pad.SetLeftMargin(0.15)
428 Pad.SetBottomMargin(0.15)
429 Pad.Draw()
430 Pad.cd()
431 resFrame.Draw()
432 legend = ROOT.TLegend(0.59, 0.6, 0.9, 0.9)
433 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(Number))
434 legend.AddEntry(0, f'#splitline{{#mu_{{#Delta t}} = {shift:4.2f}}}{{ #pm {shiftErr:4.2f} ps}}')
435 legend.AddEntry(0, f'#splitline{{#sigma_{{#Delta t}} = {resolution:4.2f}}}{{ #pm {resolutionErr:4.2f} ps}}')
436 legend.SetTextSize(0.054)
437 legend.SetFillColorAlpha(ROOT.kWhite, 0)
438 legend.Draw()
439 Pad.Update()
440 nPlot = PATH + "/test6_CPVResDeltaT" + VXDReq + ".pdf"
441 c1.SaveAs(nPlot)
442 c1.Clear()
443
444 iResult.append([f'mu = {shift:4.2f} +- {shiftErr:4.2f} ps',
445 f'sigma = {resolution:4.2f} +- {resolutionErr:4.2f} ps'])
446 fitResultsForNtuple.append(shift)
447 fitResultsForNtuple.append(resolution)
448
449 resFrameDtErr.SetTitle("")
450 sXtitleLandau = "#sigma_{#Deltat} / ps"
451 resFrameDtErr.GetXaxis().SetTitle(sXtitleLandau)
452 resFrameDtErr.GetXaxis().SetTitleSize(0.05)
453 resFrameDtErr.GetXaxis().SetLabelSize(0.045)
454 resFrameDtErr.GetYaxis().SetTitleSize(0.05)
455 resFrameDtErr.GetYaxis().SetTitleOffset(1.5)
456 resFrameDtErr.GetYaxis().SetLabelSize(0.045)
457
458 c1.cd()
459 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
460 Pad.SetLeftMargin(0.15)
461 Pad.SetBottomMargin(0.15)
462 Pad.Draw()
463 Pad.cd()
464 resFrameDtErr.Draw()
465 legend = ROOT.TLegend(0.59, 0.6, 0.9, 0.9)
466 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(Number))
467 legend.AddEntry(0, f'#splitline{{#mu_{{#Delta t}} = {meanCBS.getVal():4.2f}}}{{ #pm '
468 f'{meanCBS.getError():4.2f} ps}}') # '{:>6}'.format(Shift)
469 legend.AddEntry(0, f'#splitline{{#sigma_{{#Delta t}} = {sigmaCBS.getVal():4.2f}}}'
470 f'{{ #pm {sigmaCBS.getError():4.2f} ps}}') # '{:>4}'.format(Resol)
471 legend.SetTextSize(0.054)
472 legend.SetFillColorAlpha(ROOT.kWhite, 0)
473 # legend.Draw()
474 Pad.Update()
475 nPlot = PATH + "/test6_CPVResDeltaTError" + VXDReq + ".pdf"
476 c1.SaveAs(nPlot)
477 c1.Clear()
478
479 fitRes.Clear()
480 model.Clear()
481
482 # Fit of Delta Z for B0_sig
483
484 Mu1SigZ = ROOT.RooRealVar("Mu1SigZ", "Mu1SigZ", -9e-06, -limZSig, limZSig)
485 Mu2SigZ = ROOT.RooRealVar("Mu2SigZ", "Mu2SigZ", -1.8e-05, -limZSig, limZSig)
486 Mu3SigZ = ROOT.RooRealVar("Mu3SigZ", "Mu3SigZ", 6e-05, -limZSig, limZSig)
487 Sigma1SigZ = ROOT.RooRealVar("Sigma1SigZ", "Sigma1SigZ", 3e-03, 1e-6, limZSig)
488 Sigma2SigZ = ROOT.RooRealVar("Sigma2SigZ", "Sigma2SigZ", 1.4e-03, 1e-6, limZSig)
489 Sigma3SigZ = ROOT.RooRealVar("Sigma3SigZ", "Sigma3SigZ", 0.01, 1e-6, limZSig)
490 frac1SigZ = ROOT.RooRealVar("frac1SigZ", "frac1SigZ", 0.2, 0.0, 0.5)
491 frac2SigZ = ROOT.RooRealVar("frac2SigZ", "frac2SigZ", 0.75, 0.5, 1.)
492
493 g1SigZ = ROOT.RooGaussian("g1", "g1", DSigZ, Mu1SigZ, Sigma1SigZ)
494 g2SigZ = ROOT.RooGaussian("g2", "g2", DSigZ, Mu2SigZ, Sigma2SigZ)
495 g3SigZ = ROOT.RooGaussian("g3", "g3", DSigZ, Mu3SigZ, Sigma3SigZ)
496
497 argset1SigZ = ROOT.RooArgSet(g1SigZ)
498 argset2SigZ = ROOT.RooArgSet(g2SigZ)
499 argset3SigZ = ROOT.RooArgSet(g3SigZ)
500
501 modelSigZ = ROOT.RooAddPdf(
502 "modelSigZ", "modelSigZ", ROOT.RooArgList(g1SigZ, g2SigZ, g3SigZ),
503 ROOT.RooArgList(frac1SigZ, frac2SigZ))
504
505 DSigZ.setRange("fitRange", -limZSig, limZSig)
506
507 fitResSigZ = modelSigZ.fitTo(fitDataSigZ, ROOT.RooFit.NumCPU(1), ROOT.RooFit.Save())
508
509 fitResSigZ.Print("v")
510
511 resFrameSigZ = DSigZ.frame()
512
513 f1SigZ = frac1SigZ.getVal()
514 f2SigZ = frac2SigZ.getVal()
515 f3SigZ = 1 - f1SigZ - f2SigZ
516
517 argsetList = []
518 argsetList.append([argset1SigZ, f1SigZ])
519 argsetList.append([argset2SigZ, f2SigZ])
520 argsetList.append([argset3SigZ, f3SigZ])
521
522 argsetList = sorted(argsetList, key=itemgetter(1))
523
524 fitDataSigZ.plotOn(resFrameSigZ)
525 modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.LineWidth(4))
526 modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.Components(argsetList[2][0]), ROOT.RooFit.LineColor(ROOT.kRed + 2),
527 ROOT.RooFit.LineWidth(4))
528 modelSigZ.plotOn(resFrameSigZ, ROOT.RooFit.Components(argsetList[1][0]), ROOT.RooFit.LineColor(ROOT.kGreen + 3),
529 ROOT.RooFit.LineWidth(4))
530 modelSigZ.plotOn(
531 resFrameSigZ,
532 ROOT.RooFit.Components(argsetList[0][0]),
533 ROOT.RooFit.LineColor(
534 ROOT.kMagenta + 2),
535 ROOT.RooFit.LineWidth(4))
536 resFrameSigZ.SetTitle("")
537
538 sXtitleSigZ = "B0_Z - Gen. B0_Z / cm"
539
540 resFrameSigZ.GetXaxis().SetTitle(sXtitleSigZ)
541 resFrameSigZ.GetXaxis().SetLimits(-limZSig, limZSig)
542 resFrameSigZ.GetXaxis().SetTitleSize(0.05)
543 resFrameSigZ.GetXaxis().SetLabelSize(0.045)
544 resFrameSigZ.GetYaxis().SetTitleSize(0.05)
545 resFrameSigZ.GetYaxis().SetTitleOffset(1.5)
546 resFrameSigZ.GetYaxis().SetLabelSize(0.045)
547
548 shiftSigZ = Mu1SigZ.getVal() * f1SigZ + Mu2SigZ.getVal() * f2SigZ + Mu3SigZ.getVal() * f3SigZ
549 resolutionSigZ = Sigma1SigZ.getVal() * f1SigZ + Sigma2SigZ.getVal() * f2SigZ + Sigma3SigZ.getVal() * f3SigZ
550
551 shiftSigZ = shiftSigZ * 10000
552 resolutionSigZ = resolutionSigZ * 10000
553
554 shiftErrSigZ = math.sqrt((Mu1SigZ.getError() * f1SigZ)**2 + (Mu2SigZ.getError() * f2SigZ) ** 2 +
555 (Mu3SigZ.getError() * f3SigZ)**2) * 10000
556 resolutionErrSigZ = math.sqrt((Sigma1SigZ.getError() * f1SigZ)**2 + (Sigma2SigZ.getError() * f2SigZ)**2 +
557 (Sigma3SigZ.getError() * f3SigZ)**2) * 10000
558
559 c1.cd()
560 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
561 Pad.SetLeftMargin(0.15)
562 Pad.SetBottomMargin(0.15)
563 Pad.Draw()
564 Pad.cd()
565 resFrameSigZ.Draw()
566 legend = ROOT.TLegend(0.63, 0.65, 0.9, 0.9)
567 # NumbrSigZ = '{:d}'.format(int((f1+f2)*fitDataSigZ.numEntries()))
568 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(NumbrSigZ))
569
570 legend.AddEntry(
571 0,
572 f'#splitline{{#mu_{{#Delta z}} = {shiftSigZ:1.1f}}}'
573 f'{{ #pm {shiftErrSigZ:1.1f} #mum}}')
574 legend.AddEntry(
575 0,
576 f'#splitline{{#sigma_{{#Delta z}} = {resolutionSigZ:1.1f}}}'
577 f'{{ #pm {resolutionErrSigZ:1.1f} #mum}}')
578
579 legend.SetTextSize(0.05)
580 legend.SetFillColorAlpha(ROOT.kWhite, 0)
581 legend.Draw()
582 Pad.Update()
583 nPlot = PATH + "/test6_CPVResDeltaZsig" + VXDReq + ".pdf"
584 c1.SaveAs(nPlot)
585 c1.Clear()
586
587 fitResSigZ.Clear()
588 modelSigZ.Clear()
589
590 iResult.append([f'mu = {shiftSigZ:4.1f} +- {shiftErrSigZ:3.1f} mum',
591 f'sigma = {resolutionSigZ:4.1f} +- {resolutionErrSigZ:3.1f} mum'])
592 fitResultsForNtuple.append(shiftSigZ)
593 fitResultsForNtuple.append(resolutionSigZ)
594
595 # Fit of Delta z for B0_Tag
596
597 Mu1TagZ = ROOT.RooRealVar("Mu1TagZ", "Mu1TagZ", 0., -limZTag, limZTag)
598 Mu2TagZ = ROOT.RooRealVar("Mu2TagZ", "Mu2TagZ", 0., -limZTag, limZTag)
599 Mu3TagZ = ROOT.RooRealVar("Mu3TagZ", "Mu3TagZ", 0., -limZTag, limZTag)
600 Sigma1TagZ = ROOT.RooRealVar("Sigma1TagZ", "Sigma1TagZ", 2.51877e-02, 1e-6, limZTag)
601 Sigma2TagZ = ROOT.RooRealVar("Sigma2TagZ", "Sigma2TagZ", 1.54011e-02, 1e-6, limZTag)
602 Sigma3TagZ = ROOT.RooRealVar("Sigma3TagZ", "Sigma3TagZ", 1.61081e-02, 1e-6, limZTag)
603 frac1TagZ = ROOT.RooRealVar("frac1TagZ", "frac1TagZ", 1.20825e-01, 0.0, 1.)
604 frac2TagZ = ROOT.RooRealVar("frac2TagZ", "frac2TagZ", 1.10840e-01, 0.0, 1.)
605
606 g1TagZ = ROOT.RooGaussian("g1", "g1", DTagZ, Mu1TagZ, Sigma1TagZ)
607 g2TagZ = ROOT.RooGaussian("g2", "g2", DTagZ, Mu2TagZ, Sigma2TagZ)
608 g3TagZ = ROOT.RooGaussian("g3", "g3", DTagZ, Mu3TagZ, Sigma3TagZ)
609
610 argset1TagZ = ROOT.RooArgSet(g1TagZ)
611 argset2TagZ = ROOT.RooArgSet(g2TagZ)
612 argset3TagZ = ROOT.RooArgSet(g3TagZ)
613
614 modelTagZ = ROOT.RooAddPdf("modelTagZ", "modelTagZ", ROOT.RooArgList(
615 g1TagZ, g2TagZ, g3TagZ), ROOT.RooArgList(frac1TagZ, frac2TagZ))
616
617 DTagZ.setRange("fitRange", -limZTag, limZTag)
618
619 fitResTagZ = modelTagZ.fitTo(fitDataTagZ, ROOT.RooFit.NumCPU(1), ROOT.RooFit.Save())
620 fitResTagZ.Print("v")
621
622 resFrameTagZ = DTagZ.frame()
623
624 f1TagZ = frac1TagZ.getVal()
625 f2TagZ = frac2TagZ.getVal()
626 f3TagZ = 1 - f1TagZ - f2TagZ
627
628 argsetList = []
629 argsetList.append([argset1TagZ, f1TagZ])
630 argsetList.append([argset2TagZ, f2TagZ])
631 argsetList.append([argset3TagZ, f3TagZ])
632
633 argsetList = sorted(argsetList, key=itemgetter(1))
634
635 fitDataTagZ.plotOn(resFrameTagZ)
636 modelTagZ.plotOn(resFrameTagZ, ROOT.RooFit.LineWidth(4))
637 modelTagZ.plotOn(
638 resFrameTagZ, ROOT.RooFit.Components(
639 argsetList[2][0]), ROOT.RooFit.LineColor(
640 ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
641 modelTagZ.plotOn(
642 resFrameTagZ,
643 ROOT.RooFit.Components(argsetList[1][0]),
644 ROOT.RooFit.LineColor(ROOT.kGreen + 3),
645 ROOT.RooFit.LineWidth(4))
646 modelTagZ.plotOn(
647 resFrameTagZ,
648 ROOT.RooFit.Components(argsetList[0][0]),
649 ROOT.RooFit.LineColor(ROOT.kMagenta + 2),
650 ROOT.RooFit.LineWidth(4))
651 resFrameTagZ.SetTitle("")
652
653 sXtitleTagZ = "B0_TagVz - Gen. B0_TagVz / cm"
654
655 resFrameTagZ.GetXaxis().SetTitle(sXtitleTagZ)
656 resFrameTagZ.GetXaxis().SetTitleSize(0.05)
657 resFrameTagZ.GetXaxis().SetLabelSize(0.045)
658 resFrameTagZ.GetYaxis().SetTitleSize(0.05)
659 resFrameTagZ.GetYaxis().SetTitleOffset(1.5)
660 resFrameTagZ.GetYaxis().SetLabelSize(0.045)
661
662 shiftTagZ = Mu1TagZ.getVal() * f1TagZ + Mu2TagZ.getVal() * f2TagZ + Mu3TagZ.getVal() * f3TagZ
663 resolutionTagZ = Sigma1TagZ.getVal() * f1TagZ + Sigma2TagZ.getVal() * f2TagZ + Sigma3TagZ.getVal() * f3TagZ
664
665 shiftTagZ = shiftTagZ * 10000
666 resolutionTagZ = resolutionTagZ * 10000
667
668 shiftErrTagZ = math.sqrt((Mu1TagZ.getError() * f1TagZ)**2 + (Mu2TagZ.getError() * f2TagZ) ** 2 +
669 (Mu3TagZ.getError() * f3TagZ)**2) * 10000
670 resolutionErrTagZ = math.sqrt((Sigma1TagZ.getError() * f1TagZ)**2 + (Sigma2TagZ.getError() * f2TagZ)**2 +
671 (Sigma3TagZ.getError() * f3TagZ)**2) * 10000
672
673 c1.cd()
674 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
675 Pad.SetLeftMargin(0.15)
676 Pad.SetBottomMargin(0.15)
677 Pad.Draw()
678 Pad.cd()
679 resFrameTagZ.Draw()
680 legend = ROOT.TLegend(0.64, 0.65, 0.9, 0.9)
681 # NumbrTagZ = '{:d}'.format(int((f1+f2)*fitDataTagZ.numEntries()))
682 # legend.AddEntry(0, 'Entries' + '{:>11}'.format(NumbrTagZ))
683
684 legend.AddEntry(0, f'#splitline{{#mu_{{#Delta z}} = {shiftTagZ:1.1f}'
685 f'}}{{ #pm {shiftErrTagZ:1.1f} #mum}}')
686 legend.AddEntry(0, f'#splitline{{#sigma_{{#Delta z}} = {resolutionTagZ:1.1f}'
687 f'}}{{ #pm {resolutionErrTagZ:1.1f} #mum}}')
688 legend.SetTextSize(0.05)
689 legend.SetFillColorAlpha(ROOT.kWhite, 0)
690 legend.Draw()
691 Pad.Update()
692 nPlot = PATH + f"/test6_CPVResDeltaZtag{VXDReq}.pdf"
693 c1.SaveAs(nPlot)
694 c1.Clear()
695
696 iResult.append([f'mu = {shiftTagZ:4.1f} +- {shiftErrTagZ:3.1f} mum',
697 f'sigma = {resolutionTagZ:4.1f} +- {resolutionErrTagZ:3.1f} mum'])
698 fitResultsForNtuple.append(shiftTagZ)
699 fitResultsForNtuple.append(resolutionTagZ)
700
701 fitResults.append(iResult)
702
703fitResultsForNtuple.append(float((numberOfEntries[1] / numberOfEntries[0]) * 100))
704outputNtuple.Fill(array.array('f', fitResultsForNtuple))
705outputNtuple.Write()
706outputFile.Close()
707
708print('*********************** FIT RESULTS ***************************')
709print('* *')
710print('* WEIGHTED AVERAGES OF THE PARAMETERS *')
711print('* OF THE FITTED 3 GAUSSIAN *')
712print('* *')
713print('**************** WITHOUT PXD HIT REQUIREMENT ******************')
714print('* *')
715print('* DeltaT - Gen. DeltaT *')
716print('* *')
717print('* ' + fitResults[0][0][0] + ' ' + fitResults[0][0][1] + ' *')
718print('* *')
719print('* SigZ - Gen. SigZ *')
720print('* *')
721print('* ' + fitResults[0][1][0] + ' ' + fitResults[0][1][1] + ' *')
722print('* *')
723print('* TagZ - Gen. TagZ *')
724print('* *')
725print('* ' + fitResults[0][2][0] + ' ' + fitResults[0][2][1] + ' *')
726print('* *')
727print('********REQUIRING BOTH MUON TRACKS TO HAVE A PXD HIT***********')
728print('* *')
729print('* Efficiency *')
730print('* *')
731print(f'* N_{VXDReqs[1]}/N_{VXDReqs[0]} = {numberOfEntries[1]}/{numberOfEntries[0]} = '
732 f'{(numberOfEntries[1] / numberOfEntries[0]) * 100:3.2f}% *')
733print('* *')
734print('* DeltaT - Gen. DeltaT *')
735print('* *')
736print('* ' + fitResults[1][0][0] + ' ' + fitResults[1][0][1] + ' *')
737print('* *')
738print('* SigZ - Gen. SigZ *')
739print('* *')
740print('* ' + fitResults[1][1][0] + ' ' + fitResults[1][1][1] + ' *')
741print('* *')
742print('* TagZ - Gen. TagZ *')
743print('* *')
744print('* ' + fitResults[1][2][0] + ' ' + fitResults[1][2][1] + ' *')
745print('* *')
746print('***************************************************************')