Belle II Software development
2_analysis_pi0.py
1#!/usr/bin/env python3
2
3
10
11"""
12<header>
13 <input>../GenericB_GENSIMRECtoDST.dst.root</input>
14 <output>Pi0_Validation.root</output>
15 <contact>Mario Merola (mario.merola@na.infn.it)</contact>
16 <description>
17 Check the calibration of the ECL in the MC by determining the measured pi0 invariant mass.
18 </description>
19</header>
20"""
21
22import basf2
23import ROOT
24from modularAnalysis import cutAndCopyList, inputMdst
25from stdPi0s import stdPi0s
26from validation_tools.metadata import create_validation_histograms
27from validation_tools.metadata import validation_metadata_update
28from variables import variables as vm
29
30INPUT_FILENAME = "../GenericB_GENSIMRECtoDST.dst.root"
31OUTPUT_FILENAME = "Pi0_Validation.root"
32
33main = basf2.Path()
34inputMdst(INPUT_FILENAME, path=main)
35
36stdPi0s('all', path=main)
37cutAndCopyList('pi0:rec', 'pi0:all', 'daughter(0, E)>0.05 and daughter(1, E)>0.05', path=main)
38cutAndCopyList('pi0:mc', 'pi0:all', 'mcErrors<1', path=main)
39
40vm.addAlias('Mreco', 'M')
41
42
43create_validation_histograms(
44 main, OUTPUT_FILENAME, "pi0:rec",
45 variables_1d=[
46 (
47 "Mreco", 40, 0.08, 0.18,
48 "#pi^{0} reconstructed candidates, invariant mass",
49 "Eldar Ganiev <eldar.ganiev@desy.de>",
50 r"The $\pi^0$ invariant mass distribution with $E_{\gamma}>0.05\, \text{GeV}$",
51 r"Distribution should be peaking at the nominal $\pi^0$ mass.",
52 "M(#pi^{0}) [GeV/c^{2}]", "Candidates", "shifter"
53 ),
54 ],
55 description=r"$\pi^0$ reconstructed mass distribution",
56)
57
58
59vm.addAlias('Mmc', 'M')
60
61create_validation_histograms(
62 main, OUTPUT_FILENAME, "pi0:mc",
63 variables_1d=[
64 (
65 "Mmc", 40, 0.08, 0.18,
66 "#pi^{0} MC candidates, invariant mass",
67 "Eldar Ganiev <eldar.ganiev@desy.de>",
68 r"The $\pi^0$ invariant mass distribution for truth matched candidates",
69 r"Distribution should be peaking at the nominal $\pi^0$ mass.",
70 "M(#pi^{0}) [GeV/c^{2}]", "Candidates", "shifter"
71 ),
72 ],
73 description=r"$\pi^0$ MC mass distribution",
74)
75
76main.add_module('Progress')
77basf2.process(main)
78
79
80f = ROOT.TFile(OUTPUT_FILENAME)
81Mrecohist = f.Get('Mreco')
82Mmchist = f.Get('Mmc')
83
84
85mass = ROOT.RooRealVar("recomass", "m_{#gamma#gamma} [GeV/c^{2}]", 0.11, 0.15)
86
87h_pi0_reco = ROOT.RooDataHist("h_pi0_reco", "h_pi0_reco", ROOT.RooArgList(mass), Mrecohist)
88h_pi0_mc = ROOT.RooDataHist("h_pi0_mc", "h_pi0_mc", ROOT.RooArgList(mass), Mmchist)
89
90
91# pi0 signal PDF is a Gaussian (Crystal Ball also listed in case we want to switch)
92mean = ROOT.RooRealVar("mean", "mean", 0.125, 0.11, 0.15)
93sig1 = ROOT.RooRealVar("#sigma", "sig", 0.007, 0.002, 0.1)
94gau1 = ROOT.RooGaussian("gau1", "gau1", mass, mean, sig1)
95
96alphacb = ROOT.RooRealVar("alphacb", "alpha", 1.5, 0.1, 1.9)
97ncb = ROOT.RooRealVar("ncb", "n", 8) # ,2.,15)
98sigcb = ROOT.RooCBShape("sigcb", "sigcb", mass, mean, sig1, alphacb, ncb)
99
100# pi0 background PDF is a 2nd order Chebyshev
101b1 = ROOT.RooRealVar("b1", "b1", 0.1, -1, 1)
102a1 = ROOT.RooRealVar("a1", "a1", 0.1, -1, 1)
103bList = ROOT.RooArgList(a1, b1)
104bkg = ROOT.RooChebychev("bkg", "bkg", mass, bList)
105
106
107nsig = ROOT.RooRealVar("nsig", "nsig", 3000, 0, 1000000)
108nbkg = ROOT.RooRealVar("nbkg", "nbkg", 12000, 0, 1000000)
109
110
111totalPdf = ROOT.RooAddPdf("totalpdf", "", ROOT.RooArgList(gau1, bkg), ROOT.RooArgList(nsig, nbkg))
112
113
114output = ROOT.TFile("Pi0_Validation_ntuple.root", "recreate")
115
116# Store pi0 mass fit results to a tuple for comparison of mean and width among releases.
117outputNtuple = ROOT.TNtuple(
118 "pi0_mass",
119 "Pi0 mass fit results",
120 "mean:meanerror:width:widtherror:mean_MC:meanerror_MC:width_MC:widtherror_MC")
121
122
123ROOT.gROOT.SetBatch(True)
124canvas = ROOT.TCanvas("canvas", "pi0 mass fit", 1000, 600)
125canvas.Divide(2, 1)
126canvas.cd(1)
127
128# Fit to the reco mass
129totalPdf.fitTo(h_pi0_reco, ROOT.RooFit.Extended(True), ROOT.RooFit.Minimizer("Minuit2", "Migrad"))
130frame1 = mass.frame()
131h_pi0_reco.plotOn(frame1, ROOT.RooFit.Name("Hist"))
132frame1.SetMaximum(frame1.GetMaximum())
133totalPdf.plotOn(frame1, ROOT.RooFit.Name("curve"))
134totalPdf.plotOn(frame1, ROOT.RooFit.Components("gau1"), ROOT.RooFit.LineStyle(ROOT.kDashed), ROOT.RooFit.LineColor(ROOT.kRed))
135totalPdf.plotOn(frame1, ROOT.RooFit.Components("bkg"), ROOT.RooFit.LineStyle(3), ROOT.RooFit.LineColor(ROOT.kBlue))
136frame1.SetMaximum(Mrecohist.GetMaximum() * 1.5)
137frame1.GetXaxis().SetTitleOffset(1.4)
138frame1.GetYaxis().SetTitleOffset(1.5)
139meanval = mean.getVal()
140meanerror = mean.getError()
141width = sig1.getVal()
142widtherror = sig1.getError()
143frame1.Draw("")
144
145canvas.cd(2)
146
147# ---------------------------
148# Fit to the truth matched mass.
149# The same signal parametrisation as for the reco mass but only a 1st order Chebyshev polynomial are used.
150# Re-initialize the fit parameters.
151mean.setVal(0.125)
152sig1.setVal(0.007)
153nsig.setVal(2200)
154nbkg.setVal(700)
155a1.setVal(-0.5)
156
157bkg = ROOT.RooChebychev("bkg1", "bkg", mass, a1)
158totalPdf = ROOT.RooAddPdf("totalpdfMC", "", ROOT.RooArgList(gau1, bkg), ROOT.RooArgList(nsig, nbkg))
159
160# Fit to the truth matched mass
161totalPdf.fitTo(h_pi0_mc, ROOT.RooFit.Extended(True), ROOT.RooFit.Minimizer("Minuit2", "Migrad"))
162frame2 = mass.frame()
163h_pi0_mc.plotOn(frame2, ROOT.RooFit.Name("Hist"))
164frame2.SetMaximum(frame2.GetMaximum())
165totalPdf.plotOn(frame2, ROOT.RooFit.Name("curve"))
166totalPdf.plotOn(frame2, ROOT.RooFit.Components("gau1"), ROOT.RooFit.LineStyle(ROOT.kDashed), ROOT.RooFit.LineColor(ROOT.kRed))
167totalPdf.plotOn(frame2, ROOT.RooFit.Components("bkg1"), ROOT.RooFit.LineStyle(3), ROOT.RooFit.LineColor(ROOT.kBlue))
168frame2.SetMaximum(Mmchist.GetMaximum() * 1.5)
169frame2.GetXaxis().SetTitleOffset(1.4)
170frame2.GetYaxis().SetTitleOffset(1.5)
171meanval_mc = mean.getVal()
172meanerror_mc = mean.getError()
173width_mc = sig1.getVal()
174widtherror_mc = sig1.getError()
175frame2.Draw("")
176
177canvas.Write()
178
179outputNtuple.Fill(meanval, meanerror, width, widtherror, meanval_mc, meanerror_mc, width_mc, widtherror_mc)
180
181outputNtuple.Write()
182
183validation_metadata_update(
184 output,
185 "pi0_mass",
186 title="Pi0 mass fit results",
187 contact="eldar.ganiev@desy.de",
188 description="Fit to the invariant mass of the reconstructed and truth matched pi0s",
189 check="Consistent numerical fit results. Stable mean and width.",
190 metaoptions="shifter")
191
192
193output.Close()