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