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>
10 Check the calibration of the ECL in the MC by determining the measured pi0 invariant mass.
15 INPUT_FILENAME =
"../GenericB_GENSIMRECtoDST.dst.root"
16 OUTPUT_FILENAME =
"Pi0_Validation.root"
21 from modularAnalysis
import inputMdst, reconstructDecay, matchMCTruth, cutAndCopyList
22 from stdPhotons
import stdPhotons
23 from stdPi0s
import stdPi0s
26 from variables
import variables
as vm
29 inputMdst(
'default', INPUT_FILENAME, path=main)
32 matchMCTruth(
'pi0:all', path=main)
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)
37 vm.addAlias(
'Mreco',
'M')
40 create_validation_histograms(
41 main, OUTPUT_FILENAME,
"pi0:rec",
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"
52 description=
r"$\pi^0$ reconstructed mass distribution",
56 vm.addAlias(
'Mmc',
'M')
58 create_validation_histograms(
59 main, OUTPUT_FILENAME,
"pi0:mc",
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"
70 description=
r"$\pi^0$ MC mass distribution",
74 print(basf2.statistics)
77 f = ROOT.TFile(OUTPUT_FILENAME)
78 Mrecohist = f.Get(
'Mreco')
79 Mmchist = f.Get(
'Mmc')
82 mass = ROOT.RooRealVar(
"recomass",
"m_{#gamma#gamma} [GeV/c^{2}]", 0.11, 0.15)
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)
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)
93 alphacb = ROOT.RooRealVar(
"alphacb",
"alpha", 1.5, 0.1, 1.9)
94 ncb = ROOT.RooRealVar(
"ncb",
"n", 8)
95 sigcb = ROOT.RooCBShape(
"sigcb",
"sigcb", mass, mean, sig1, alphacb, ncb)
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)
104 nsig = ROOT.RooRealVar(
"nsig",
"nsig", 1000, 0, 1000000)
105 nbkg = ROOT.RooRealVar(
"nbkg",
"nbkg", 1000, 0, 1000000)
108 totalPdf = ROOT.RooAddPdf(
"totalpdf",
"", ROOT.RooArgList(sigcb, bkg), ROOT.RooArgList(nsig, nbkg))
111 output = ROOT.TFile(
"Pi0_Validation_ntuple.root",
"recreate")
114 outputNtuple = ROOT.TNtuple(
116 "Pi0 mass fit results",
117 "mean:meanerror:width:widtherror")
120 canvas = ROOT.TCanvas(
"canvas",
"pi0 mass fit", 1000, 600)
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))
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()
144 outputNtuple.Fill(meanval, meanerror, width, widtherror)
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)
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))
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()
186 outputNtuple.Fill(meanval_mc, meanerror_mc, width_mc, widtherror_mc)
190 validation_metadata_update(
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")